To do spline interpolation, we need a series of knots and a pair of slopes.  This is a demo program 

showing an interpolation for a 5 knot (four segment) spline curve. 

>	with( LinearAlgebra ): with( plots ):

 

With the knots matrix and the end point slopes, our spline becomes fully defined.   

>	knots := < <-3,2> | <-1,0.5> | <0,0> | <1,0.5> | <3,2> >; endSlopes := Vector( [-1,1] );

 

All we really care about is the coefficients for the cubic polynomials which make up our splines.   

With five knots, we have four cubic polynomials, each with four polynomials, giving us 16  

unknowns.  From the definition of a cubic spline, we build a matrix representing our relationships 

between the cubic polynomials and solve for the unknown coefficients.  We have:
 

Matrix A: Cubic polynomial equations
Vector xv: Coefficients for cubic polynomials (what we want to solve for)
Vector bv: Solutions (conditions) for each line in matrix A 

>	numKnots := ColumnDimension( knots ): unknowns := ( numKnots - 1 ) * 4: A := Matrix( unknowns, unknowns ): xv := Vector( unknowns ): bv := Vector( unknowns ):

 

These are some general helper structures which will clean up our later computations. 

>	cubicPoly := Vector( [1,x,x^2,x^3] ): cubicOneDiff := Vector( [0,1,2*x,3*x^2] ): cubicTwoDiff := Vector( [0,0,2,6*x] ): currentRow := 1:

 

For each segment, assert f_i( x_i ) = y_i 

>	for j from 1 to ( numKnots - 1 ) do    currentKnot := Column( knots, j ):    A[ currentRow, (j * 4 - 3)..(j * 4) ] := Transpose( eval( cubicPoly, x=currentKnot[ 1 ] ) ):    bv[ currentRow ] := currentKnot[ 2 ]:    currentRow := currentRow + 1: end do:

 

For each segment assert f_i( x_i+1) = y_i+1 

>	for j from 1 to ( numKnots - 1 ) do    currentKnot := Column( knots, j+1 ):    A[ currentRow, (j * 4 - 3)..(j * 4) ] := Transpose( eval( cubicPoly, x=currentKnot[ 1 ] ) ):    bv[ currentRow ] := currentKnot[ 2 ]:    currentRow := currentRow + 1: end do:

 

Where each segment meets, assert f_i ' ( x_i ) - f_i+1 ' ( x_i ) = 0 

>

for j from 1 to ( numKnots - 2 ) do    currentKnot := Column( knots, j+1 ):    A[ currentRow, (j * 4 - 3)..(j * 4) ] := Transpose( eval( cubicOneDiff, x=currentKnot[ 1 ] ) ):    A[ currentRow, (j * 4 + 1)..(j * 4 + 4) ] := (-1) * Transpose( eval( cubicOneDiff, x=currentKnot[ 1 ] ) ):    bv[ currentRow ] := 0:    currentRow := currentRow + 1: end do:

 

Where each segment meets, assert f_i '' (x_i) - f_i+1 '' (x_i) = 0 

>

for j from 1 to ( numKnots - 2 ) do    currentKnot := Column( knots, j+1 ):    A[ currentRow, (j * 4 - 3)..(j * 4) ] := Transpose( eval( cubicTwoDiff, x=currentKnot[ 1 ] ) ):    A[ currentRow, (j * 4 + 1)..(j * 4 + 4) ] := (-1) * Transpose( eval( cubicTwoDiff, x=currentKnot[ 1 ] ) ):    bv[ currentRow ] := 0:    currentRow := currentRow + 1: end do:

 

Finally, assert our end point slope conditions 

>

currentKnot := Column( knots, 1 ): A[ currentRow, 1..4 ] := Transpose( eval( cubicOneDiff, x=currentKnot[ 1 ] ) ): bv[ currentRow ] := endSlopes[ 1 ]: currentRow := currentRow + 1: currentKnot := Column( knots, numKnots ): A[ currentRow, (unknowns-3)..(unknowns) ] := Transpose( eval( cubicOneDiff, x=currentKnot[ 1 ] ) ): bv[ currentRow ] := endSlopes[ 2 ]:

 

Solve the coefficient matrix to find the coefficients for our components cubic polynomials. 

>	xv := A^(-1) . bv:

 

Represent our coefficients as cubic polynomials and build a piecwise function (the spline curve). 

>	f := 0: for j from 1 to ( numKnots - 1 ) do    f := piecewise( x < knots[ 1, j ], f,                    x >= knots[ 1, j ] and x < knots[ 1, j+1 ], xv[ (j * 4 - 3)..(j*4) ] . cubicPoly ): end do: f := simplify( f );

 

Plot our spline curve. 

>	( xMin, xMax ) := rtable_scanblock( Row( knots, 1 ), [1..numKnots],'Minimum','Maximum'): ( yMin, yMax ) := rtable_scanblock( Row( knots, 2 ), [1..numKnots],'Minimum','Maximum'): plot( f, x=(xMin-2)..(xMax+2), y=(yMin-2)..(yMax+2));

 

>