Splines

Here we are going to construct Splines by solving the appropriate equations.

Lets get a list of points to work with.

list={{2,1},{4,3},{5,6},{7,3}};
plot1=ListPlot[list,PlotStyle->{Red,PointSize[0.05]},PlotRange->All]

n=Length[list]-1;
Do[x[k]=list[[k+1,1]];y[k]=list[[k+1,2]],{k,0,n}]

Now we construct the cubic polynomials to find.

S[i_,xx_]=a[i]+b[i] (xx-x[i])+c[i] (((xx-x[i]))^2)+d[i] (((xx-x[i]))^3)

Now, lets set up the equations.  The first two sets are for the values of S.  The third and fourth equate the first and second derivatives where the polynomials meet.  The last set of equations give the two addition restraints required (in our case, the free splines).

equ1=Table[S[j,x[j]]==y[j],{j,0,n-1}]

equ2=Table[S[j,x[j+1]]==y[j+1],{j,0,n-1}]

equ3=Table[Derivative[0,1][S][j,x[j+1]]==Derivative[0,1][S][j+1,x[j+1]],{j,0,n-2}]

equ4=Table[Derivative[0,2][S][j,x[j+1]]==Derivative[0,2][S][j+1,x[j+1]],{j,0,n-2}]

equ5={Derivative[0,2][S][0,x[0]]==0,Derivative[0,2][S][n-1,x[n]]==0}

equs=Simplify[Flatten[{equ1,equ2,equ3,equ4,equ5}]]

Now lets solve the equations.

solution=Solve[equs]

Do[S1[j,x_]=S[j,x]/.solution[[1]],{j,0,n-1}]
Table[N[S1[i,x]],{i,0,n-1}]

Lets plot the solution and the points it goes through.

curve[j_]:=Plot[S1[j,xx],{xx,x[j],x[j+1]}];
Show[plot1,Table[curve[j],{j,0,n-1}]]

Let us do the same with a longer list of points, and do it all in one go.  These points are generated from sin(x), with a bit of noise.

list=Table[{x,Sin[x]+Random[Real,{-0.1,0.1}]},{x,-3.0,3.0,0.5}];
plot1=ListPlot[list,PlotStyle->{Red,PointSize[0.02]},PlotRange->All];
plot2=Plot[Sin[x],{x,-3,3},PlotStyle->Green];
n=Length[list]-1;
Do[x[k]=list[[k+1,1]];y[k]=list[[k+1,2]],{k,0,n}]
S[i_,xx_]=a[i]+b[i] (xx-x[i])+c[i] (((xx-x[i]))^2)+d[i] (((xx-x[i]))^3);
equ1=Table[S[j,x[j]]==y[j],{j,0,n-1}];
equ2=Table[S[j,x[j+1]]==y[j+1],{j,0,n-1}];
equ3=Table[Derivative[0,1][S][j,x[j+1]]==Derivative[0,1][S][j+1,x[j+1]],{j,0,n-2}];
equ4=Table[Derivative[0,2][S][j,x[j+1]]==Derivative[0,2][S][j+1,x[j+1]],{j,0,n-2}];
equ5={Derivative[0,2][S][0,x[0]]==0,Derivative[0,2][S][n-1,x[n]]==0};
equs=Simplify[Flatten[{equ1,equ2,equ3,equ4,equ5}]];
solution=Solve[equs];
Do[S1[j,x_]=S[j,x]/.solution[[1]],{j,0,n-1}]
curve[j_]:=Plot[S1[j,xx],{xx,x[j],x[j+1]}];
Show[plot1,plot2,Table[curve[j],{j,0,n-1}]]

By the way, the algorithm in the text book is much faster than the computations done here.  This is because they solve the equations more intelligently than Mathematica does.  They use techniques from Chapter 6.