spatial/r3: add Divergence
MODIFY
spatial/r3/vector.go
spatial/r3/vector.go
@@ -83,6 +83,18 @@ func Cos(p, q Vec) float64 {
return Dot(p, q) / (Norm(p) * Norm(q))
}
+// Divergence returns the divergence of the vector field at the point p,
+// approximated using finite differences with the given step sizes.
+func Divergence(p, step Vec, field func(Vec) Vec) float64 {
+ sx := Vec{X: step.X}
+ divx := (field(Add(p, sx)).X - field(Sub(p, sx)).X) / step.X
+ sy := Vec{Y: step.Y}
+ divy := (field(Add(p, sy)).Y - field(Sub(p, sy)).Y) / step.Y
+ sz := Vec{Z: step.Z}
+ divz := (field(Add(p, sz)).Z - field(Sub(p, sz)).Z) / step.Z
+ return 0.5 * (divx + divy + divz)
+}
+
// minElem return a vector with the minimum components of two vectors.
func minElem(a, b Vec) Vec {
return Vec{

MODIFY
spatial/r3/vector_test.go
spatial/r3/vector_test.go
@@ -8,6 +8,8 @@ import (
"math"
"testing"
+ "golang.org/x/exp/rand"
+
"gonum.org/v1/gonum/floats/scalar"
"gonum.org/v1/gonum/mat"
)
@@ -254,6 +256,54 @@ func TestRotate(t *testing.T) {
}
}
+var vectorFields = []struct {
+ field func(Vec) Vec
+ divergence func(Vec) float64
+}{
+ {
+ field: func(v Vec) Vec {
+ return Vec{X: v.X * v.Y * v.Z, Y: v.Y * v.Z, Z: v.Z * v.X}
+ },
+ divergence: func(v Vec) float64 {
+ return v.X + v.Y*v.Z + v.Z
+ },
+ },
+ {
+ field: func(v Vec) Vec {
+ sx := math.Sin(v.X)
+ sy, cy := math.Sincos(v.Y)
+ return Vec{
+ X: v.X * v.Y * v.Z * cy,
+ Y: v.Y*v.Z + sx,
+ Z: v.Z * v.X / sy,
+ }
+ },
+ divergence: func(v Vec) float64 {
+ sy, cy := math.Sincos(v.Y)
+ return v.X/sy + v.Y*v.Z*cy + v.Z
+ },
+ },
+}
+
+func TestDivergence(t *testing.T) {
+ const (
+ tol = 1e-10
+ h = 1e-2
+ )
+ step := Vec{X: h, Y: h, Z: h}
+ rnd := rand.New(rand.NewSource(1))
+ for _, test := range vectorFields {
+ for i := 0; i < 30; i++ {
+ p := randomVec(rnd)
+ got := Divergence(p, step, test.field)
+ want := test.divergence(p)
+ if math.Abs(got-want) > tol {
+ t.Errorf("result out of tolerance. got %v, want %v", got, want)
+ }
+ }
+ }
+}
+
func vecDense(v Vec) *mat.VecDense {
return mat.NewVecDense(3, []float64{v.X, v.Y, v.Z})
}

spatial/r3: add Divergence