Natural Cubic Spline

A natural cubic spline S(x) is a piecewise cubic polynomial function. Let a=(Subscript[x, 0]) < (Subscript[x, 1]) < ... <  (Subscript[x, n-1]) < (Subscript[x, n]) = b be a partition of the interval [a,b]. Then S(x) satisfies the following conditions:  
     (1) S(x) is defined on the interval [a, b];  
     (2) S(x) is a cubic polynomial on each interval [(Subscript[x, i]), (Subscript[x, i+1])],  i = 0, 1, ... n-1; 
     (3) S(x), S'(x), S''(x) are continuous on [a, b]; 
     (4) S''(a) = 0, S''(b) = 0. 
If a function satisifies the first three conditions, then it is a cubic spline. The fourth conditon is a characteristic of natural cubic spline.  

Given a set of data ((Subscript[x, 0]), (Subscript[y, 0])), ((Subscript[x, 1]), (Subscript[y, 1])) ..., ((Subscript[x, n]), (Subscript[y, n])), we can construct a natural cubic interpolation spline that satisifes a further interpolation condition: 
                                  S((Subscript[x, i])) = (Subscript[y, i]),         i = 0, 1, ..., n. 
The following is an algorithm for constructing a natural cubic interpolation spline.

Algorithm

Let XY ={((Subscript[x, 0]), (Subscript[y, 0])), ((Subscript[x, 1]), (Subscript[y, 1])) ..., ((Subscript[x, n]), (Subscript[y, n]))}  be the data, (Subscript[h, k]) =(Subscript[x, k+1])- (Subscript[x, k]). We will solve (Subscript[z, 1]), (Subscript[z, 2]),...., (Subscript[z, n-1]) from the tridiagonal system 
   (({{({{2((Subscript[h, 1])+(Subscript[h, 2])), (Subscript[h, 1])}, {(Subscript[h, 1]), 2((Subscript[h, 1])+(Subscript[h, 2]))}})({{ ,  }, {(Subscript[h, 2]),  }})
          ...                     ...          ....,   }, {  , ({{ 2((Subscript[h, n-2])+(Subscript[h, n-1])), (Subscript[h, n-2])}, {(Subscript[h, n-2]), 2((Subscript[h, n-1])+(Subscript[h, n]))}})}})) (({{ ({{(Subscript[z, 1])}, {(Subscript[z, 2])}}) }, {({{...}, {(Subscript[z, n-1])}})}})) = (({{ ({{(Subscript[v, 1])}, {(Subscript[v, 2])}}) }, {({{...}, {(Subscript[v, n-1])}})}}))        where    (Subscript[v, i]) = 6 ((((Subscript[y, i+1])-(Subscript[y, i]))/(Subscript[h, i]))-(((Subscript[y, i])-(Subscript[y, i-1]))/(Subscript[h, i-1]))). We set (Subscript[z, 0])= 0 and (Subscript[z, n])=0. 
The spline function S(x) is then given by the formula 
               S(x) = ((Subscript[z, i+1])/(6(Subscript[h, i])))(((x -(Subscript[t, i])))^3)+((Subscript[z, i])/(6(Subscript[h, i])))((((Subscript[t, i+1])-x))^3)+(((Subscript[y, i+1])/(Subscript[h, i]))-(((Subscript[h, i])(Subscript[z, i+1]))/6))(x -(Subscript[t, i]))+(((Subscript[y, i])/(Subscript[h, i]))-(((Subscript[h, i])(Subscript[z, i]))/6))((Subscript[t, i+1])-x)     for  (Subscript[t, i])≤x≤(Subscript[t, i+1]).

Algorithm

(*First we define the data and the tridiagonal system*)
Difference[XY0_] :=   
Block[{k},        
  XY = XY0; n = Length[XY]-1; X = Transpose[XY][[1]]; Y = Transpose[XY][[2]];
  h = Table[0,{n}]; d = Table[0,{n}]; z = Table[0,{n+1}]; a = Table[0,{n-1}]; 
  b = Table[0,{n-1}]; c = Table[0,{n-1}]; v = Table[0,{n-1}]; s = Table[0,{n},{4}];
  h[[1]] = X[[2]] - X[[1]]; d[[1]] = (Y[[2]]-Y[[1]])/h[[1]];
  For[k=2,k<=n,k++,
    h[[k]] = X[[k+1]] - X[[k]]; d[[k]] = (Y[[k+1]]-Y[[k]])/h[[k]];
    a[[k-1]] = h[[k]]; b[[k-1]] = 2 (h[[k-1]] + h[[k]]);
    c[[k-1]] = h[[k]]; v[[k-1]] = 6(d[[k]]-d[[k-1]])];];

(*next define the boundary conditons (Subscript[z, 0])=0 and (Subscript[z, n+1])=0*)
Coeffs:= Block[{}, z[[1]]  = 0; z[[n+1]] = 0];

(*This solves the tri-diagonal system of equations*)
TriDiagonal[a0_,b0_,c0_,v0_] :=   
Module[{k,t},
  a = a0; b=b0; c=c0; v=v0;
  For[k=2,k<=n-1,k++,  t = a[[k-1]]/b[[k-1]]; b[[k]] = b[[k]] - t c[[k-1]]; 
  v[[k]] = v[[k]] - t v[[k-1]];];  z[[n]] = v[[n-1]]/b[[n-1]];
  For[k=n-2,1<=k,k--,   z[[k+1]] = (v[[k]] -  c[[k]] z[[k+2]])/b[[k]];];];
     
(*We can then make up the spline function*)
ComputeCoeff[d_,h_,z_,y_] := (*coefficients of the cubic polynomials*)    
Block[{k},
  For[k=1,k<=n,k++,  s[[k,1]] = Y[[k]];  s[[k,2]] = d[[k]]-h[[k]](2 z[[k]]+  z[[k+1]])/6;
    s[[k,3]] = z[[k]]/2; s[[k,4]] = (z[[k+1]] - z[[k]])/(6 h[[k]]) ]; ];
P[m_,x_] :=    (*cubic polynomials that makes up the spline function*)
((s[[m,4]](x - X[[m]])+s[[m,3]])*(x-X[[m]])+s[[m,2]])(x - X[[m]])+ s[[m,1]] ;
Clear[Sm, SP]  (*Sm[x] is the m-th piece of the spline function, zero outside [(Subscript[x, m]),(Subscript[x, m+1])]*)
S1[x_] := If[x< X[[2]],P[1,x],0] ;
Sm[m_,x_] :=  If[X[[m]]<=x && x<X[[m+1]],P[m,x], 0]; 
Sn[n_,x_]:= If[X[[n]]<=x,P[n,x],0] ; 
SP[x_] := S1[x] +Sum[Sm[i,x], {i,2,n-1}] + Sn[n,x] (*SP[x] is the spline function*)
 
MakeSpline := (*This excutes all the above steps*)
Block[{},   Difference[XY]; TriDiagonal[a,b,c,v]; ComputeCoeff[d,h,z,Y] ];

To use the algorithm, all it needs is to have a given set of data XY in the form of  XY ={((Subscript[x, 0]), (Subscript[y, 0])), ((Subscript[x, 1]), (Subscript[y, 1])) ..., ((Subscript[x, n]), (Subscript[y, n]))}, then run   
                  MakeSpline
This commands executes the algorithm that solves the  tri-diagonal system and computes each piece of the cubic polynomial (denoted by P[m,x]) that makes up the spline function (denoted by SP[x]).

Example

Example 1.  Construct a natural cubic interpolation spline for the given data:

XY = {{0.,1.},{1.,0.},{2.,0.}, {3.,1.},{4.,2.}, {5.,2.},{6.,1.}}

The command MakeSpline will generate the natural cubic interpolaltion spline for this data set. The result is shown in the following plot:

MakeSpline;
dots=ListPlot[XY,PlotStyle->{PointSize[0.02`]},DisplayFunction->Identity];gr=Plot[SP[x],{x,0,6},DisplayFunction->Identity];txt ="Natural Cubic spline";Show[gr,dots,DisplayFunction->$DisplayFunction,PlotLabel->txt,AxesLabel->{"x","y"}]

The cubic polynomials that makes up the cubic spline function are printed out below:

Print[P[1,x]]; Print[P[2,x]]; Print[P[3,x]]; Print[P[4,x]]; Print[P[5,x]]; Print[P[6,x]];

We can plot all these polynomials together, which shows how they make up the spline function:

piece=Plot[{P[1,x], P[2,x],P[3,x],P[4,x],P[5,x],P[6,x]},{x, 0,6}, 
       PlotRange->{{0,6},{-0.2,2.2}},DisplayFunction->Identity];
Show[piece,dots,DisplayFunction->$DisplayFunction]

Example 2.  Construct a natural cubic interpolation spline for the given data:

XY = {{0.,0.8},{0.5,0.8},{1.,1.2},{8.,3.},{11.,6.},{18.,6.}, {20.5,3.5},{24.,3.}, {25.,0}}

MakeSpline;dots=ListPlot[XY,PlotStyle->{PointSize[0.02`]},DisplayFunction->Identity];gr=Plot[SP[x],{x,0,26},PlotRange->{{0,26},{0,8.1`}},AspectRatio->0.25`,DisplayFunction->Identity];txt ="Natural spline function";Show[gr,dots,DisplayFunction->$DisplayFunction,PlotLabel->txt,AxesLabel->{"x","y"}]

Example 3.  Construct a natural cubic interpolation spline that interpolate the function f (x) = 1/(1+25 (x^2)) over [-1,1] using equally spaced points: :

f[x_] = 1 /(1+ 25x^2);
XY= Table[{0.2 i, f[0.2 i]},{i, -5, 5}]

MakeSpline;dots=ListPlot[XY,PlotStyle->{PointSize[0.02`]},DisplayFunction->Identity];F=Plot[f[x],{x,-1,1},DisplayFunction->Identity];gr=Plot[SP[x],{x,-1,1},DisplayFunction->Identity]; Show[gr,dots,DisplayFunction->$DisplayFunction,PlotLabel->txt,AxesLabel->{"x","y"}]