Module for Cubic Splines

Numerical Methods using Mathematica  
Complementary software supplements for   
Numerical Methods textbooks  
John H. Mathews and Kurtis Fink, 2004

Cubic Spline Interpolant

Definition (Cubic Spline).  Suppose that  {((Subscript[x, k]), (Subscript[y, k]))}_(k = 0)^n  are  n+1  points, where  a = (Subscript[x, 0]) < (Subscript[x, 0]) <  ...  < (Subscript[x, n]) = b.  The function  S (x)  is called a cubic spline if there exists  n  cubic polynomials  (Subscript[S, k]) (x)  with coefficients  (Subscript[s, k, 0]),(Subscript[s, k, 1]),(Subscript[s, k, 2]), and  (Subscript[s, k, 3])  that satisfy the properties:

I.  	S (x)  =  (Subscript[S, k]) (x) = (Subscript[s, k, 0]) + (Subscript[s, k, 1]) (x - 
                (Subscript[x, k])) + (Subscript[s, k, 2]) (((x - (Subscript[x, k])))^2) + (Subscript[s, k, 
                3]) (((x - (Subscript[x, k])))^3)  
	for  x∈[(Subscript[x, k]), (Subscript[x, k + 1])]  and  k=0,1,...,n - 1.  

II.  	S ((Subscript[x, k]))  =  (Subscript[y, k])  for  k = 0, 1,  ... , n.  
	The spline passes through each data point.

III.  	(Subscript[S, k]) ((Subscript[x, k + 1])) = (Subscript[S, k + 
                1]) ((Subscript[x, k + 1]))  for  k = 0, 1,  ... , n - 2. 
	The spline forms a continuous function over [a,b].

IV.  	\(\(S\_k\%'\) \((x\_\(k + 1\))\) = \(S\_\(k + 
                1\)\%'\) \((x\_\(k + 1\))\)\)  for  k = 0, 1,  ... , n - 2. 
	The spline forms a smooth function.

IV.  	\(\(S\_k\%''\) \((x\_\(k + 1\))\) = \(S\_\(k + 
                1\)\%''\) \((x\_\(k + 1\))\)\)  for  k = 0, 1,  ... , n - 2. 
	The second derivative is continuous.

Lemma (Natural Spline).  There exists a unique cubic spline with the free boundary conditions  S'' (a) = 0  and  S'' (b) = 0.

Proof  Cubic Splines  Cubic Splines

Remark.  The natural spline is the curve obtained by forcing a flexible elastic rod through the points but letting the slope at the ends be free to equilibrate to the position that minimizes the oscillatory behavior of the curve.  It is useful for fitting a curve to experimental data that is significant to several significant digits.

Computer Programs  Cubic Splines  Cubic Splines

Program (Natural Cubic Spline).  To construct and evaluate the cubic spline interpolant S (x) for the  n + 1 data points  ((Subscript[x, 0]), (Subscript[y, 0])), ((Subscript[x, 1]), (Subscript[y, 1])),  ... , ((Subscript[x, n]), (Subscript[y, n])),  using the free boundary conditions  S'' (a) = 0  and  S'' (b) = 0.

Mathematica Subroutine (Natural Cubic Spline).

NaturalSpline[XY0_] := Module[{XY = XY0}, 
Differences :=  
                  Module[{k},  
n  =  
                      Length[XY] - 
                        1;  
X  =  (Subscript[(Transpose[XY]), [[1]]]);  
Y  =  (Subscript[(Transpose[XY]), [[2]]]);  
h = d  =  
                        Table[0, {n}];  
m  =  
                      Table[0, {n + 
                            1}];  
a = b = c = v  =  Table[0, {n - 1}];  
s  =  
                      Table[0, {n}, {4}];  
(Subscript[h, [[1]]])  =  
                      (Subscript[X, [[2]]])  -  
                        (Subscript[X, [[1]]]);  
(Subscript[d, [[1]]])  =  (((Subscript[Y, [[2]]])  -  (Subscript[Y, [[1]]]))/(Subscript[h, [[1]]]));  
For[ k = 2,  
                      k ≤ 
                        n,  k++,  
(Subscript[h, [[k]]])  =  
                        (Subscript[X, [[k + 1]]])  -  
                          (Subscript[X, [[k]]]);  
(Subscript[d, [[k]]])  =  (((Subscript[Y, [[k + 1]]])  -  (Subscript[Y, [[k]]]))/(Subscript[h, [[k]]]));  
(Subscript[a, [[k - 1]]])  =  
                        (Subscript[h, [[k]]]);  
(Subscript[b, [[k - 1]]])  =  
                        2 ((Subscript[h, [[k - 1]]])  +  
                              (Subscript[h, [[k]]]));  
(Subscript[c, [[k - 1]]])  =  
                        (Subscript[h, [[k]]]);  
(Subscript[v, [[k - 1]]])  =  
                        6 ((Subscript[d, [[k]]])  -  
                              (Subscript[d, [[k - 1]]]))]; ]; 
Print[h];
Print[d];
TriDiagonal := 
                  Module[{k, 
                      t},  
(Subscript[m, [[1]]])  =  
                      0;  
(Subscript[m, [[n + 1]]])  =  0; 
For[ k = 2,  
                      k ≤ 
                        n - 1,  k++,  
t  =  
                        ((Subscript[a, [[k - 1]]])/(Subscript[b, [[k - 1]]]));  
(Subscript[b, [[k]]])   =   (Subscript[b, [[k]]])  -  
                          t  (Subscript[c, [[k - 1]]]);  
(Subscript[v, [[k]]])   =   (Subscript[v, [[k]]])  -  
                          t  (Subscript[v, [[k - 1]]]); ];  
(Subscript[m, [[n]]])   =   ((Subscript[v, [[n - 1]]])/(Subscript[b, [[n - 1]]]));  
For[ k = n - 2,  
                      1 ≤ 
                        k,  k--,  
(Subscript[m, [[k + 1]]])  =  (((Subscript[v, [[k]]])  -  (Subscript[c, [[k]]]) (Subscript[m, [[k + 2]]]))/(Subscript[b, [[k]]])); ]; ]; 
ComputeCoeff := 
                  Module[{k},  
For[ k = 1,  
                        k ≤ 
                          n,  k++,  
(Subscript[s, [[k, 1]]])  =  
                          (Subscript[Y, [[k]]]);  
(Subscript[s, [[k, 2]]])  =  
                          (Subscript[d, [[k]]])  -  
                            (1/6) (Subscript[h, [[k]]]) (2 (Subscript[m, [[k]]]) + 
                                  (Subscript[m, [[k + 1]]]));  
(Subscript[s, [[k, 3]]])  =  
                          ((Subscript[m, [[k]]])/2);  
(Subscript[s, [[k, 4]]])  =  (((Subscript[m, [[k + 1]]])  -  (Subscript[m, [[k]]]))/(6 (Subscript[h, [[k]]]))); ]; ]; 
CS[t_] := 
                  Module[{j},  
For[ j = 1,  
                      j ≤ 
                        n,  j++,  
If[ 
                          (Subscript[X, [[j]]]) ≤ t  &&  
                            t < 
                              (Subscript[X, [[j + 1]]]), k = j ]; ];  
If[
                      t < (Subscript[X, [[1]]]), k = 1 ];  
If[
                      (Subscript[X, [[n + 1]]]) ≤ t, 
                      k = n ];  
w  =  
                      t  -  
                        (Subscript[X, [[k]]]);  
Return[ (((Subscript[s, [[k, 4]]]) w + 
                                    (Subscript[s, [[k, 3]]])) w + 
                              (Subscript[s, [[k, 2]]])) w + 
                        (Subscript[s, [[k, 1]]]) ]; ]; 
 (*The main part of the program starts here*) ;  
Differences; 
TriDiagonal; 
ComputeCoeff; ];

Example 1.  Construct the natural cubic spline for the points (0, 1), (1, 0), (2, 0), (3, 1), (4, 2), (5, 
        2), (6, 1)  that has the endpoint constraints   S'' (0)=0  and  S''(6) = 0.

Solution 1.

XY = {{0. , 1. }, {1. , 0. }, {2. , 0. }, {3. , 1. }, {4. , 
            2. }, {5. , 2. }, {6. , 1. }};  
NaturalSpline[XY];  
Needs["Graphics`Colors`"];   
dots = 
        ListPlot[XY, PlotStyle -> {Red, PointSize[0.02]}, 
          DisplayFunction -> Identity];  
gr = 
        Plot[CS[x], {x, 0, 6}, PlotStyle -> {Magenta}, 
          DisplayFunction -> Identity];  
Show[gr, dots, AxesLabel -> {"x", "y    "}, 
        DisplayFunction -> $DisplayFunction];  
Print["Spline y = ", (Subscript[S, i]), "[x]"];  
Print["data points = ", XY];

If more information is needed, we can look at the coefficients that form the cubic spline.

\(\(\(\(Print["The spline coefficients are"];\)\ \ \ 

    MatrixForm[Chop[s]]\)\(\ \ \ \)\)\)

2nd Solution:

s1[x_]  =  
        a1 (x^3)  +  b1 (x^2)  +  c1 x  +  d1;  
s2[x_]  =  
        a2 (x^3)  +  b2 (x^2)  +  c2 x  +  d2;   
s3[x_]  =  
        a3 (x^3)  +  b3 (x^2)  +  c3 x  +  d3;   
s4[x_]  =  
        a4 (x^3)  +  b4 (x^2)  +  c4 x  +  
          d4;   
s5[x_]  =  
        a5 (x^3)  +  b5 (x^2)  +  c5 x  +  d5;   
s6[x_]  =  
        a6 (x^3)  +  b6 (x^2)  +  c6 x  +  
          d6;   
Print["S1[x_] = ", s1[x] ];  
Print["S2[x_] = ", s2[x] ];  
Print["S3[x_] = ", s3[x] ];  
Print["S4[x_] = ", s4[x] ];  
Print["S5[x_] = ", s5[x] ];  
Print["S6[x_] = ", s6[x] ];

eq1 = s1[0] == 1; eq2 = s1[1] == 0; eq3 = s2[1] == 0; eq4 = s2[2] == 0; 
    eq5 = s3[2] == 0; eq6 = s3[3] == 1; eq7 = s4[3] == 1; 
    eq8 = s4[4] == 2; eq9 = s5[4] == 2; eq10 = s5[5] == 2; eq11 = s6[5] == 2; 
    eq12 = s6[6] == 1;
eqf13 = s1'[1] == s2'[1]; eqf14 = s2'[2] == s3'[2]; 
    eqf15 = s3'[3] == s4'[3]; eqf16 = s4'[4] == s5'[4]; 
    eqf17 = s5'[5] == s6'[5];
eqs18 = s1''[1] == s2''[1]; eqs19 = s2''[2] == s3''[2]; 
    eqs20 = s3''[3] == s4''[3]; eqs21 = s4''[4] == s5''[4]; 
    eqs22 = s5''[5] == s6''[5];

eqleftf = s1'[0] == 0; 
    eqrightf = s6'[6] == 0;

eqns  =  { eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8 , eq9 , 
          eq10 , eq11 , eq12, eqf13, eqf14, eqf15, eqf16, eqf17, eqs18, 
          eqs19, eqs20, eqs21, eqs22,  eqleftf, 
          eqrightf};

vars = {a1, b1, c1, d1, a2, b2, c2, d2, a3, b3, c3, d3, a4, b4, c4, 
          d4, a5, b5, c5, d5, a6, b6, c6, d6};

solset  =  Solve[eqns, vars];   
Print[TableForm[eqns]];   
Print[solset];

f1[x_]  =  ReplaceAll[ s1[x], solset[[1]] ];  
f2[x_]  =  ReplaceAll[ s2[x], solset[[1]]];  
f3[x_]  =  
        ReplaceAll[ s3[x], solset[[1]]]; 
f4[x_]  =  
        ReplaceAll[ s4[x], solset[[1]]]; 
f5[x_]  =  
        ReplaceAll[ s5[x], solset[[1]]]; 
f6[x_]  =  
        ReplaceAll[ s6[x], solset[[1]]]; 
\(\(Print[\*""\!\(f\_1\)[x] = "", f1[x]\ ];\)\ \ \)
\(\(Print[\*""\!\(f\_2\)[x] = "", 
        f2[x]\ ];\)\ \ \)
\(\(Print[\*""\!\(f\_3\)[x] = "", 
        f3[x]\ ];\)\ \)
\(\(Print[\*""\!\(f\_4\)[x] = "", 
        f4[x]\ ];\)\)
\(\ \(\(Print[\*""\!\(f\_5\)[x] = "", 
      f5[x]\ ]\)\(;\)\(\ \)\(Print[\*""\!\(f\_6\)[x] = "", 
      f6[x]\ ]\)\(;\)\(\ \)\)\)

xy = { {0, f1[0]}, {1, f1[1]}, {2, f2[2]}, {3, f3[3]}, {4, 
            f4[4]}, {5, f5[5]}, {6, f6[6]} };  
dots  =  
        ListPlot[xy, PlotStyle -> {Red, PointSize[0.02]} , 
          DisplayFunction -> Identity]; 
graph1  =  
        Plot[f1[x], {x, 0, 1}, PlotStyle -> Green, 
          DisplayFunction -> Identity]; 
graph2  =  
        Plot[f2[x], {x, 1, 2}, PlotStyle -> Blue, 
          DisplayFunction -> Identity];
graph3  =  
        Plot[f3[x], {x, 2, 3}, PlotStyle -> Orange, 
          DisplayFunction -> Identity];
graph4  =  
        Plot[f4[x], {x, 3, 4}, PlotStyle -> Green, 
          DisplayFunction -> Identity];
graph5 =  
        Plot[f5[x], {x, 4, 5}, PlotStyle -> Blue, 
          DisplayFunction -> Identity];
graph6  =  
        Plot[f6[x], {x, 5, 6}, PlotStyle -> Orange, 
          DisplayFunction -> Identity];
graph12  =  
        Show[graph1, graph2, graph3, graph4, graph5, graph6, dots, 
          AspectRatio -> (2.2/3), 
          DisplayFunction -> $DisplayFunction];
\(\(Print[\*""\!\(f\_1\)[x] = "", f1[x]\ ];\)\ \ \)
\(\(Print[\*""\!\(f\_2\)[x] = "", 
        f2[x]\ ];\)\ \ \)
\(\(Print[\*""\!\(f\_3\)[x] = "", 
        f3[x]\ ];\)\ \)
\(\(Print[\*""\!\(f\_4\)[x] = "", 
        f4[x]\ ];\)\)
\(\ \(\(Print[\*""\!\(f\_5\)[x] = "", 
      f5[x]\ ]\)\(;\)\(\ \)\(Print[\*""\!\(f\_6\)[x] = "", 
      f6[x]\ ]\)\(;\)\(\ \)\)\)

Example 2.  Construct the natural cubic spline for the points (0, 1), (1, 0), (2, 0), (3, 1), (4, 2), (5, 
        2), (6, 1)  that has the endpoint constraints   S'' (0)=0  and  S''(6) = 0.  
Use Mathematica's procedure SplineFit.

Solution 2.

XY = {{0. , 1. }, {1. , 0. }, {2. , 0. }, {3. , 1. }, {4. , 
            2. }, {5. , 2. }, {6. , 1. }};  
Needs["NumericalMath`SplineFit`"];  
spline = SplineFit[XY, Cubic];  
ParametricPlot[spline[t], {t, 0, 6}, PlotRange -> All, 
        Compiled -> False, PlotStyle -> Blue];

Remark.  There are five popular types of splines: natural spline, clamped spline, extrapolated spline, parabolically terminated spline, endpoint curvature adjusted spline.
When Mathematica constructs a cubic spline it uses the "natural cubic spline."

Clamped Spline.

Lemma (Clamped Spline).  There exists a unique cubic spline with the first derivative boundary conditions  S' (a) = (Subscript[d, 0])  and  S' (b) = (Subscript[d, n]).

Proof  Cubic Splines  Cubic Splines

A property of clamped cubic splines.

	A practical feature of splines is the minimum of the oscillatory behavior they possess.  Consequently, among all functions f(x) which are twice continuously differentiable on [a,b] and interpolate a given set data points {((Subscript[x, k]), f ((Subscript[x, k])))}_(k = 0)^n , the cubic spline has "less wiggle."  The next result explains this phenomenon. 

Theorem (Minimum property of clamped cubic splines).  Assume that  f ∈ (C^2)[a, b]  and  S (x)  is the unique clamped cubic spline interpolant for  f (x) which passes through {((Subscript[x, k]), f ((Subscript[x, k])))}_(k = 0)^n and satisfies the clamped end conditions  S' (a) = f' (a)  and  S' (b) = f' (b).   Then

		\(\(\(∫\_a\%b\ \(\((S'' \((x)\))\)\^2\) dx\)\(\ \ \)\(≤\)\(\ \ \)\)\)(∫_a^b) (((f'' (x)))^2) dx.

Proof  Cubic Splines  Cubic Splines

Computer Programs  Cubic Splines  Cubic Splines

Program (Clamped Cubic Spline).  To construct and evaluate the cubic spline interpolant  S(x)  for the  n+1  data points  ((Subscript[x, 0]), (Subscript[y, 0])), ((Subscript[x, 1]), (Subscript[y, 1])),  ... , ((Subscript[x, n]), (Subscript[y, n])),  using the first derivative boundary conditions  S' (a) = (Subscript[d, 0])  and  S' (b) = (Subscript[d, n]).

ClampedSpline[XY0_, D1_, DN1_] := Module[{XY = XY0, dx1 = D1, dxn1 = DN1}, 
Differences :=  
                  Module[{k},  
n  =  
                      Length[XY] - 
                        1;  
X  =  (Subscript[(Transpose[XY]), [[1]]]);  
Y  =  (Subscript[(Transpose[XY]), [[2]]]);  
h = d  =  
                        Table[0, {n}];  
m  =  
                      Table[0, {n + 
                            1}];  
a = b = c = v  =  Table[0, {n - 1}];  
s  =  
                      Table[0, {n}, {4}];  
(Subscript[h, [[1]]])  =  
                      (Subscript[X, [[2]]])  -  
                        (Subscript[X, [[1]]]);  
(Subscript[d, [[1]]])  =  (((Subscript[Y, [[2]]])  -  (Subscript[Y, [[1]]]))/(Subscript[h, [[1]]]));  
For[ k = 2,  
                      k ≤ 
                        n,  k++,  
(Subscript[h, [[k]]])  =  
                        (Subscript[X, [[k + 1]]])  -  
                          (Subscript[X, [[k]]]);  
(Subscript[d, [[k]]])  =  (((Subscript[Y, [[k + 1]]])  -  (Subscript[Y, [[k]]]))/(Subscript[h, [[k]]]));  
(Subscript[a, [[k - 1]]])  =  
                        (Subscript[h, [[k]]]);  
(Subscript[b, [[k - 1]]])  =  
                        2 ((Subscript[h, [[k - 1]]])  +  
                              (Subscript[h, [[k]]]));  
(Subscript[c, [[k - 1]]])  =  
                        (Subscript[h, [[k]]]);  
(Subscript[v, [[k - 1]]])  =  
                        6 ((Subscript[d, [[k]]])  -  
                              (Subscript[d, [[k - 1]]]))]; ];  
TriDiagonal := 
                  Module[{k, 
                      t},  
(Subscript[b, [[1]]])  =  
                      (Subscript[b, [[1]]])  -  
                        ((Subscript[h, [[1]]])/2);  
(Subscript[v, [[1]]])  =  (Subscript[v, [[1]]])  -  
                        3 ((Subscript[d, [[1]]])  -  
                              dx1);  
(Subscript[b, [[n - 1]]])  =  
                      (Subscript[b, [[n - 1]]]) - 
                        ((Subscript[h, [[n]]])/2);  
(Subscript[v, [[n - 1]]])  =  
                      (Subscript[v, [[n - 1]]]) - 
                        3 (dxn1 - 
                              (Subscript[d, [[n]]]));  
For[ k = 2,  
                      k ≤ 
                        n - 1,  k++,  
t  =  
                        ((Subscript[a, [[k - 1]]])/(Subscript[b, [[k - 1]]]));  
(Subscript[b, [[k]]])   =   (Subscript[b, [[k]]])  -  
                          t  (Subscript[c, [[k - 1]]]);  
(Subscript[v, [[k]]])   =   (Subscript[v, [[k]]])  -  
                          t  (Subscript[v, [[k - 1]]]); ];  
(Subscript[m, [[n]]])   =   ((Subscript[v, [[n - 1]]])/(Subscript[b, [[n - 1]]]));  
For[ k = n - 2,  
                      1 ≤ 
                        k,  k--,  
(Subscript[m, [[k + 1]]])  =  (((Subscript[v, [[k]]])  -  (Subscript[c, [[k]]]) (Subscript[m, [[k + 2]]]))/(Subscript[b, [[k]]])); ]; ];  
ComputeCoeff := 
                  Module[{k},  
For[ k = 1,  
                        k ≤ 
                          n,  k++,  
(Subscript[m, [[1]]])  =  ((3 ((Subscript[d, [[1]]])  -  dx1))/(Subscript[h, [[1]]]))  -  
                            ((Subscript[m, [[2]]])/2);  
(Subscript[m, [[n + 1]]])  =  ((3 (dxn1  -  (Subscript[d, [[n]]])))/(Subscript[h, [[n]]]))  -  
                            ((Subscript[m, [[n]]])/2);  
(Subscript[s, [[k, 1]]])  =  
                          (Subscript[Y, [[k]]]);  
(Subscript[s, [[k, 2]]])  =  
                          (Subscript[d, [[k]]])  -  
                            (1/6) (Subscript[h, [[k]]]) (2 (Subscript[m, [[k]]]) + 
                                  (Subscript[m, [[k + 1]]]));  
(Subscript[s, [[k, 3]]])  =  
                          ((Subscript[m, [[k]]])/2);  
(Subscript[s, [[k, 4]]])  =  (((Subscript[m, [[k + 1]]])  -  (Subscript[m, [[k]]]))/(6 (Subscript[h, [[k]]]))); ]; ];  
CS[t_] := 
                  Module[{j},  
For[ j = 1,  
                      j ≤ 
                        n,  j++,  
If[ 
                          (Subscript[X, [[j]]]) ≤ t  &&  
                            t < (Subscript[X, [[j + 1]]]), k = j ]; ];  
If[
                      t < (Subscript[X, [[1]]]), k = 1 ];  
If[
                      (Subscript[X, [[n + 1]]]) ≤ t, 
                      k = n ];  
w  =  
                      t  -  
                        (Subscript[X, [[k]]]);  
Return[ (((Subscript[s, [[k, 4]]]) w + 
                                    (Subscript[s, [[k, 3]]])) w + 
                              (Subscript[s, [[k, 2]]])) w + 
                        (Subscript[s, [[k, 1]]]) ]; ]; 
 (*The main part of the program starts here*) ;  
Differences; 
TriDiagonal; 
ComputeCoeff; ];

Example 3.  Construct the clamped cubic spline for the points (0, 1), (1, 0), (2, 0), (3, 1), (4, 2), (5, 
        2), (6, 1)  that has the endpoint constraints   S' (0)=-0.2  and  S'(6) = -2.5.

Solution 3.

XY = {{0. , 1. }, {1. , 0. }, {2. , 0. }, {3. , 1. }, {4. , 
            2. }, {5. , 2. }, {6. , 1. }};  
ClampedSpline[
        XY, -0.2, -2.5];  
Needs["Graphics`Colors`"];  
dots = 
        ListPlot[XY, PlotStyle -> {Red, PointSize[0.02]}, 
          DisplayFunction -> Identity];  
gr = 
        Plot[CS[x], {x, 0, 6}, PlotStyle -> {Magenta}, 
          DisplayFunction -> Identity];  
Show[gr, dots, DisplayFunction -> $DisplayFunction, 
        AxesLabel -> {"x", "y    "}];  
Print["Spline y = ", (Subscript[S, i]), "[x]"];  
Print["data points = ", XY];

If more information is needed, we can look at the coefficients that form the cubic spline.

Print["The spline coefficients are"];  
MatrixForm[Chop[s]]

Research Experience for Undergraduates

Cubic Splines  Cubic Splines  Internet hyperlinks to web sites and a bibliography of articles.