spatial/r3: add Jacobian
MODIFY
spatial/r3/mat.go
spatial/r3/mat.go
@@ -267,3 +267,28 @@ func (m *Mat) Hessian(p, step Vec, field func(Vec) float64) {
m.Set(2, 1, fyz)
m.Set(2, 2, fzz)
}
+
+// Jacobian sets the receiver to the Jacobian matrix of the vector field at the point p,
+// approximated using finite differences with the given step sizes.
+// Jacobian expects the field's first order partial
+// derivatives are all continuous for correct results.
+func (m *Mat) Jacobian(p, step Vec, field func(Vec) Vec) {
+ dx := Vec{X: step.X}
+ dy := Vec{Y: step.Y}
+ dz := Vec{Z: step.Z}
+
+ dfdx := Scale(0.5/step.X, Sub(field(Add(p, dx)), field(Sub(p, dx))))
+ dfdy := Scale(0.5/step.Y, Sub(field(Add(p, dy)), field(Sub(p, dy))))
+ dfdz := Scale(0.5/step.Z, Sub(field(Add(p, dz)), field(Sub(p, dz))))
+ m.Set(0, 0, dfdx.X)
+ m.Set(0, 1, dfdy.X)
+ m.Set(0, 2, dfdz.X)
+
+ m.Set(1, 0, dfdx.Y)
+ m.Set(1, 1, dfdy.Y)
+ m.Set(1, 2, dfdz.Y)
+
+ m.Set(2, 0, dfdx.Z)
+ m.Set(2, 1, dfdy.Z)
+ m.Set(2, 2, dfdz.Z)
+}

MODIFY
spatial/r3/mat_test.go
spatial/r3/mat_test.go
@@ -302,3 +302,24 @@ func TestMatHessian(t *testing.T) {
}
}
}
+
+func TestMatJacobian(t *testing.T) {
+ const (
+ tol = 1e-5
+ h = 8e-4
+ )
+ step := Vec{X: h, Y: h, Z: h}
+ rnd := rand.New(rand.NewSource(1))
+ for _, test := range vectorFields {
+ for i := 0; i < 1; i++ {
+ p := randomVec(rnd)
+ got := NewMat(nil)
+ got.Jacobian(p, step, test.field)
+ want := test.jacobian(p)
+ if !mat.EqualApprox(got, want, tol) {
+ t.Errorf("matrices not equal within tol\ngot: %v\nwant: %v",
+ mat.Formatted(got), mat.Formatted(want))
+ }
+ }
+ }
+}

MODIFY
spatial/r3/vector_test.go
spatial/r3/vector_test.go
@@ -259,28 +259,47 @@ func TestRotate(t *testing.T) {
var vectorFields = []struct {
field func(Vec) Vec
divergence func(Vec) float64
+ jacobian func(Vec) *Mat
}{
{
field: func(v Vec) Vec {
+ // (x*y*z, y*z, z*x)
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
},
+ jacobian: func(v Vec) *Mat {
+ return NewMat([]float64{
+ v.Y * v.Z, v.X * v.Z, v.X * v.Y,
+ 0, v.Z, v.Y,
+ v.Z, 0, v.X,
+ })
+ },
},
{
field: func(v Vec) Vec {
+ // (x*y*z*cos(y), y*z+sin(x), z*x*sin(y))
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,
+ 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
+ return v.X*sy + v.Y*v.Z*cy + v.Z
+ },
+ jacobian: func(v Vec) *Mat {
+ cx := math.Cos(v.X)
+ sy, cy := math.Sincos(v.Y)
+ return NewMat([]float64{
+ v.Y * v.Z * cy, v.X*v.Z*cy - v.X*v.Y*v.Z*sy, v.X * v.Y * cy,
+ cx, v.Z, v.Y,
+ v.Z * sy, v.X * v.Z * cy, v.X * sy,
+ })
},
},
}

spatial/r3: add Jacobian