spatial/r3: add Gradient
MODIFY
spatial/r3/mat_test.go
spatial/r3/mat_test.go
@@ -235,6 +235,7 @@ func BenchmarkQuat(b *testing.B) {
var scalarFields = []struct {
// This is the scalar field function.
field func(p Vec) float64
+ gradient func(p Vec) Vec
hessian func(p Vec) *Mat
}{
{
@@ -243,6 +244,9 @@ var scalarFields = []struct {
z2 := p.Z * p.Z
return 4*p.X*p.X*p.X + 5*p.Y*p.Y + 3*z2*z2
},
+ gradient: func(p Vec) Vec {
+ return Vec{X: 12 * p.X * p.X, Y: 10 * p.Y, Z: 12 * p.Z * p.Z * p.Z}
+ },
hessian: func(p Vec) *Mat {
return NewMat([]float64{
24 * p.X, 0, 0,
@@ -257,6 +261,13 @@ var scalarFields = []struct {
y2 := p.Y * p.Y
return math.Cos(p.X) * math.Sin(p.Z) * y2 * y2
},
+ gradient: func(p Vec) Vec {
+ y3 := p.Y * p.Y * p.Y
+ y4 := p.Y * y3
+ sx, cx := math.Sincos(p.X)
+ sz, cz := math.Sincos(p.Z)
+ return Vec{X: -y4 * sx * sz, Y: 4 * y3 * cx * sz, Z: y4 * cx * cz}
+ },
hessian: func(p Vec) *Mat {
y3 := p.Y * p.Y * p.Y
y4 := y3 * p.Y

MODIFY
spatial/r3/vector.go
spatial/r3/vector.go
@@ -97,6 +97,19 @@ func Divergence(p, step Vec, field func(Vec) Vec) float64 {
return 0.5 * (divx + divy + divz)
}
+// Gradient returns the gradient of the scalar field at the point p,
+// approximated using finite differences with the given step sizes.
+func Gradient(p, step Vec, field func(Vec) float64) Vec {
+ dx := Vec{X: step.X}
+ dy := Vec{Y: step.Y}
+ dz := Vec{Z: step.Z}
+ return Vec{
+ X: field(Add(p, dx)) - field(Sub(p, dx)),
+ Y: field(Add(p, dy)) - field(Sub(p, dy)),
+ Z: field(Add(p, dz)) - field(Sub(p, dz)),
+ }
+}
+
// 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
@@ -304,6 +304,25 @@ func TestDivergence(t *testing.T) {
}
}
+func TestGradient(t *testing.T) {
+ const (
+ tol = 1e-12
+ h = 1e-4
+ )
+ step := Vec{X: h, Y: h, Z: h}
+ rnd := rand.New(rand.NewSource(1))
+ for _, test := range scalarFields {
+ for i := 0; i < 30; i++ {
+ p := randomVec(rnd)
+ got := Gradient(p, step, test.field)
+ want := test.gradient(p)
+ if vecApproxEqual(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 Gradient