WIP: better spline

geom/cmd/splinetester.go

@@ -26,12 +26,17 @@ func main() {
 	if err := linear.Prepare(); err != nil {
 		log.Fatalf("linear.Prepare() = %v, want nil", err)
 	}
-	cubic := &geom.CubicSpline{Points: points}
+	cubic := &geom.CubicSpline{Points: points,
+		FixedPreslope:  true,
+		FixedPostslope: false,
+		Preslope:       -5,
+		//Postslope:      -4,
+	}
 	if err := cubic.Prepare(); err != nil {
 		log.Fatalf("cubic.Prepare() = %v, want nil", err)
 	}
 	// Produce interpolated points in CSV-like form.
-	for x := -8.0; x < 8.0; x += 0.125 {
+	for x := -8.0; x < 8.0; x += 0.0625 {
 		fmt.Printf("%f,%f,%f\n", x, linear.Interpolate(x), cubic.Interpolate(x))
 	}
 }

geom/spline.go

@@ -63,16 +63,25 @@ func (s *LinearSpline) Interpolate(x float64) float64 {
 	return s.Points[i-1].Y + (x-s.Points[i-1].X)*s.deriv[i-1]
 }
 
-// CubicSpline implements a natural cubic spline. A cubic spline interpolates
+// CubicSpline implements a cubic spline. A cubic spline interpolates
-// the given Points while ensuring first and second derivatives are continuous.
+// the given Points with cubic polynomials, ensuring first and second
+// derivatives along the whole spline are continuous.
 type CubicSpline struct {
 	Points []Float2
 
+	// If false, CubicSpline defines a natural cubic spline (the slopes at the
+	// endpoints are "free" and the moments at the ends are zero.)
+	// If true, Preslope (or Postslope, or both) is used to set the slope.
+	FixedPreslope, FixedPostslope bool
+
+	// Slope of line before and after spline, for extrapolation.
+	// If a natural cubic spline is being used, these are set by Prepare.
+	// If instead FixedPreslope or FixePostslope is true, these are read by
+	// Prepare to determine the moments.
+	Preslope, Postslope float64
+
 	// moments (second derivative at 1/6 scale) and intervals
 	m, h []float64
-
-	// slope of line before and after spline, for extrapolation
-	preslope, postslope float64
 }
 
 // Prepare sorts the points and computes internal information.
@@ -97,49 +106,117 @@ func (s *CubicSpline) Prepare() error {
 		}
 		s.h[i] = s.Points[i+1].X - s.Points[i].X
 	}
-	// Compute moments. m[0] and m[N-1] are chosen to be 0 (natural cubic spline).
+	// Compute moments. "moments" is a term from drafting that basically means
-	// Also these "moments" aren't the true values of the second derivatives
+	// the second derivative of the function. Points are known as "knots".
+	//
+	// Let's start with the natural cubic spline case, where m[0] and m[N-1] are
+	// chosen to be 0. I'll skip over the derivation of the equations below, but
+	// it follows from putting a cubic in each interval and having each one meet
+	// the knots and match second derivatives with its neighbors.
+	//
+	// Note: these "moments" aren't the true values of the second derivatives
 	// at the knots - they are calculated at 1/6th scale to avoid a multiply
 	// and divide by 6.
+	//
 	// Given:
+	//
 	//    ɣ(i) = 2.0 * (h[i-1] + h[i])
-	//    b(i) = ((Points[i+1].Y-Points[i].Y)/h[i] - (Points[i].Y-Points[i-1].Y)/h[i-1])
+	//    b(i) = ((Y[i+1]-Y[i])/h[i] - (Y[i]-Y[i-1])/h[i-1])
+	//
 	// we solve for m[i] in the equations:
+	//
 	//    h[i-1]*m[i-1] + ɣ(i)*m[i] + h[i]*m[i+1] = b(i)
+	//
 	// for i = 1...N-2.
 	//
 	// Written as a diagonally dominant tridiagonal matrix equation:
 	//
-	// [ɣ(1)  h[1]  0     0    ...    0      ] [  m[1]  ]   [  b(1)  ]
+	// [ ɣ(1)  h[1]     0     0     ...       0 ] [   m[1] ]   [   b(1) ]
-	// [h[1]  ɣ(2)  h[2]  0    ...    0      ] [  m[2]  ]   [  b(2)  ]
+	// [ h[1]  ɣ(2)  h[2]     0     ...       0 ] [   m[2] ]   [   b(2) ]
-	// [0     h[2]  ɣ(3)  h[3] ...    0      ] [  m[3]  ] = [  b(3)  ]
+	// [    0  h[2]  ɣ(3)  h[3]     ...       0 ] [   m[3] ] = [   b(3) ]
-	// [0     0     h[3]  ɣ(4) ...    ...    ] [  ...   ]   [  ...   ]
+	// [    0     0  h[3]  ɣ(4)     ...     ... ] [    ... ]   [    ... ]
-	// [...................... ...    h[N-3] ] [  ...   ]   [  ...   ]
+	// [  ...   ...   ...   ...     ...  h[N-3] ] [  .  .. ]   [    ... ]
-	// [0     0     ...   0    h[N-3] ɣ(N-2) ] [ m[N-2] ]   [ b(N-2) ]
+	// [    0     0   ...     0  h[N-3]  ɣ(N-2) ] [ m[N-2] ]   [ b(N-2) ]
 	//
 	// This is solvable in O(N) using simplified Gaussian elimination
 	// ("Thomas algorithm").
+	//
+	// For the fixed end-slopes case, we also need to derive m[0] and m[N-1]
+	// from the given end-slopes. Given:
+	//
+	//     b(0) = (Y[1] - Y[0]) / h[0] - Preslope
+	//   b(N-1) = Postslope - (Y[N-1] - Y[N-2]) / h[N-2]
+	//
+	// We solve two additional equations for the new unknowns m[0] and m[N-1]:
+	//
+	//       2*m[0]*h[0] + m[1]*h[0]       = b(0)
+	//     m[N-2]*h[N-2] + 2*m[N-1]*h[N-2] = b(N-1).
+	//
+	// Fortunately this is still a tridiagonal:
+	//
+	// [ 2h[0]  h[0]     0     0      0     ...        0 ] [   m[0] ]   [   b(0) ]
+	// [  h[0]  ɣ(1)  h[1]     0      0     ...        0 ] [   m[1] ]   [   b(1) ]
+	// [     0  h[1]  ɣ(2)  h[2]      0     ...        0 ] [   m[2] ]   [   b(2) ]
+	// [     0     0  h[2]  ɣ(3)   h[3]     ...      ... ] [   m[3] ] = [   b(3) ]
+	// [     0     0  h[3]  ɣ(4)    ...  h[N-3]        0 ] [    ... ]   [    ... ]
+	// [   ...   ...   ...   ... h[N-3]  ɣ(N-2)   h[N-2] ] [ m[N-2] ]   [ b(N-2) ]
+	// [     0     0     0   ...     0   h[N-2]  2h[N-2] ] [ m[N-1] ]   [ b(N-1) ].
+	//
+	// Fixing one end but leaving the other free leads to a mix of the two.
 
 	// Setup:
-	diag, upper, B := make([]float64, N-1), make([]float64, N-1), make([]float64, N-1)
+	diag, upper, B := make([]float64, N), make([]float64, N), make([]float64, N)
+	if s.FixedPreslope {
+		diag[0] = 2.0 * s.h[0]
+		upper[0] = s.h[0]
+		B[0] = (s.Points[1].Y-s.Points[0].Y)/s.h[0] - s.Preslope
+	}
 	for i := 1; i < N-1; i++ {
 		diag[i] = 2.0 * (s.h[i-1] + s.h[i])
 		upper[i] = s.h[i]
 		B[i] = (s.Points[i+1].Y-s.Points[i].Y)/s.h[i] - (s.Points[i].Y-s.Points[i-1].Y)/s.h[i-1]
 	}
+	if s.FixedPostslope {
+		diag[N-1] = 2.0 * s.h[N-2]
+		upper[N-1] = s.h[N-2]
+		B[N-1] = s.Postslope - (s.Points[N-1].Y-s.Points[N-2].Y)/s.h[N-2]
+	}
 	// Forward elimination:
+	if s.FixedPreslope {
+		// Use row 0 to eliminate lower h[0] from row 1.
+		// lower[1] = h[0]; diag[0] = 2h[0]; therefore lower[1]/diag[0] = 0.5.
+		diag[1] -= 0.5 * upper[0]
+		B[1] -= 0.5 * B[0]
+	}
 	for i := 2; i < N-1; i++ {
+		// Use row i-1 to eliminate lower h[i-1] from row i
 		t := s.h[i-1] / diag[i-1] // lower[i] / diag[i-1]
 		diag[i] -= t * upper[i-1]
 		B[i] -= t * B[i-1]
 	}
+	if s.FixedPostslope {
+		// Use row N-2 to eliminate lower h[N-2] from row N-1.
+		t := s.h[N-2] / diag[N-2]
+		diag[N-1] -= t * upper[N-2]
+		B[N-1] -= t * B[N-2]
+	}
 	// Back substitution:
+	if s.FixedPostslope {
+		s.m[N-1] = B[N-1] / diag[N-1]
+	}
 	for i := N - 2; i > 0; i-- {
 		s.m[i] = (B[i] - s.h[i]*s.m[i+1]) / diag[i]
 	}
-	// Pre- and post-slope:
+	if s.FixedPreslope {
-	s.preslope = -s.m[1]*s.h[0] + (s.Points[1].Y-s.Points[0].Y)/s.h[0]
+		s.m[0] = (B[0] - s.h[0]*s.m[1]) / diag[1]
-	s.postslope = s.m[N-2]*s.h[N-2] + (s.Points[N-1].Y-s.Points[N-2].Y)/s.h[N-2]
+	}
+	// Derive pre- and post-slope, if not fixed:
+	if !s.FixedPreslope {
+		s.Preslope = -s.m[1]*s.h[0] + (s.Points[1].Y-s.Points[0].Y)/s.h[0]
+	}
+	if !s.FixedPostslope {
+		s.Postslope = s.m[N-2]*s.h[N-2] + (s.Points[N-1].Y-s.Points[N-2].Y)/s.h[N-2]
+	}
 	return nil
 }
 
@@ -153,11 +230,11 @@ func (s *CubicSpline) Interpolate(x float64) float64 {
 	}
 	if x < s.Points[0].X {
 		// Comes before the start of the spline, extrapolate
-		return s.Points[0].Y + (x-s.Points[0].X)*s.preslope
+		return s.Points[0].Y + (x-s.Points[0].X)*s.Preslope
 	}
 	if x > s.Points[N-1].X {
 		// Comes after the end of the spline, extrapolate
-		return s.Points[N-1].Y + (x-s.Points[N-1].X)*s.postslope
+		return s.Points[N-1].Y + (x-s.Points[N-1].X)*s.Postslope
 	}
 	// Somewhere in the middle
 	i := sort.Search(N, func(n int) bool {

geom/spline_test.go

@@ -111,7 +111,7 @@ func TestCubicSplineOnePoint(t *testing.T) {
 	}
 }
 
-func TestCubicSpline(t *testing.T) {
+func TestNaturalCubicSpline(t *testing.T) {
 	s := &CubicSpline{
 		Points: []Float2{{-7, -2}, {-5, 1}, {-3, 0}, {-2, -3}, {0, 2}, {1, -5}, {3, -2}, {4, 4}},
 	}
@@ -156,3 +156,53 @@ func TestCubicSpline(t *testing.T) {
 		}
 	}
 }
+
+func TestFixedEndSlopesCubicSpline(t *testing.T) {
+	s := &CubicSpline{
+		Points:         []Float2{{-7, -2}, {-5, 1}, {-3, 0}, {-2, -3}, {0, 2}, {1, -5}, {3, -2}, {4, 4}},
+		FixedPreslope:  true,
+		FixedPostslope: true,
+		Preslope:       -1,
+		Postslope:      1,
+	}
+	if err := s.Prepare(); err != nil {
+		t.Errorf("s.Prepare() = %v, want nil", err)
+	}
+	tests := []struct {
+		x, want float64
+	}{
+		{x: -8, want: -3.648342225609756},
+		{x: -7.5, want: -2.824171112804878},
+		{x: -7, want: -2},
+		{x: -6.5, want: -1.180464581745427},
+		{x: -6, want: -0.3887433307926829},
+		{x: -5.5, want: 0.3473495855564025},
+		{x: -5, want: 1},
+		{x: -4.5, want: 1.5067079125381098},
+		{x: -4, want: 1.6662299923780488},
+		{x: -3.5, want: 1.2426370760289636},
+		{x: -3, want: 0},
+		{x: -2.5, want: -1.9368449885670733},
+		{x: -2, want: -3},
+		{x: -1.5, want: -1.855450886051829},
+		{x: -1, want: 0.45221989329268286},
+		{x: -0.5, want: 2.2837807259908534},
+		{x: 0, want: 2},
+		{x: 0.5, want: -1.229539824695122},
+		{x: 1, want: -5},
+		{x: 1.5, want: -6.734946646341463},
+		{x: 2, want: -6.406821646341463},
+		{x: 2.5, want: -4.6252858231707314},
+		{x: 3, want: -2},
+		{x: 3.5, want: 0.941477705792683},
+		{x: 4, want: 4},
+		{x: 4.5, want: 7.078029725609756},
+		{x: 5, want: 10.156059451219512},
+		{x: 5.5, want: 13.234089176829269},
+	}
+	for _, test := range tests {
+		if got := s.Interpolate(test.x); math.Abs(got-test.want) > 0.0000001 {
+			t.Errorf("s.Interpolate(%v) = %v, want %v", test.x, got, test.want)
+		}
+	}
+}