Mathematical Principles for Scientific Computing and Visualization  
http://www.farinhansford.com/books/scv
by Gerald Farin and Dianne Hansford

This notebook demonstrates interpolation with B-splines.
Figure 8.10 was generated with this notebook.

Compare the results here with those in Taylor.nb

Input: Knot sequence u. Data values (v_i,f(v_i))_i i=1,n ... Here we have selected 20 values from f(v) = sqrt(v) + sin(v)/10. 

Output: the least squares approximating B-spline to the input values.
Here we are using a degree 3 B-spline, and we are using five polynomial segments.
(See comments below -- you can add more segments by adding knots and see the improvement in the approximation.)


See Chapter 8: Data Fitting for details.

Updated July 2008

(SetDirectory[ToFileName[Extract["FileName"/.NotebookInformation[EvaluationNotebook[]],{1},FrontEnd`FileName]]];) 

(* Define a zero tolerance *)
eps=5.`*^-6;

(* B-spline evaluation *)
spanleft[n_,knots_,i_,x_]:=If[Abs[knots[[i+n+1]]-knots[[i+1]]]<eps,0,((x-knots[[i+1]])/(knots[[i+n+1]]-knots[[i+1]]))];

spanright[n_,knots_,i_,x_]:=If[Abs[knots[[i+n+2]]-knots[[i+2]]]<eps,0,((knots[[i+n+2]]-x)/(knots[[i+n+2]]-knots[[i+2]]))]

bspline[n_,knots_,i_,x_]:=If[x≥knots[[i+1]]&&x<knots[[i+2]],1,0];

bspline[n_,knots_,i_,x_]:=spanleft[n,knots,i,x] bspline[n-1,knots,i,x]+spanright[n,knots,i,x] bspline[n-1,knots,i+1,x]/;n>0;


(* Degree 3 B-spline *)
n=3;

(* Knot sequence *)
(* Play with adding a knot at 13.5 -- it will improve the results *)
u={0,0,0,0,4,8,12,16,20,20,20,29};

(* Places to create data *)
v={0,2,4,8,12,13.5, 16,18,20};

ldata=Length[v];

k=Length[u]-3;

(* Number of control points / coefficients *)
l=k-n+1;

(* B-spline curve *)
bsfunct[n_,u_,d_,l_,x_]:=(∑_(i=0)^l)d[[i+1]] bspline[n,u,i,x];

(* Set up design matrix *)
(* Evaluation of B-splines at each data's knot value *)
(Vander[n_,l_,u_,ldata_,v_]:=Table[bspline[3,u,i,v[[j]]],{j,1,ldata},{i,0,l}]) 

(* Function from which we will sample *)
(f[x_]:=N[(Sqrt[x])+((Sin[x])/10)];) 

(* Form matrix part of normal equations *)
(Mat[n_,l_,u_,ldata_,v_]:=Transpose[Vander[n,l,u,ldata,v]].Vander[n,l,u,ldata,v];) 


(* Evauate function at associated knot values *)
(data[ldata_]:=Table[f[v[[j]]],{j,1,ldata}];) 


(* Solve for coefficients of B-spline curve  *)(bvals[n_,l_,u_,ldata_,v_]:=Inverse[Mat[n,l,u,ldata,v]].(Transpose[Vander[n,l,u,ldata,v]].data[ldata]);) 

(* d stores the B-spline coefficients *)
(d=N[bvals[n,l,u,ldata,v]];) 

(Print["B-spline coeffs d: ",d];) 

(* Organize the knot vector locations and associated function values in a table *)
vals=Table[{v[[i]],f[v[[i]]]},{i,1,ldata}];

Print["ldata = ",ldata,"  vals = ",vals];

(* Plot the points in vals *)

ListPlot[vals,PlotStyle->PointSize[0.03`],BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11}]


PLT=Plot[{f[x],bsfunct[n,u,d,l,x]},{x,v[[1]],v[[ldata]]},PlotStyle->{{Thickness[0.01`],GrayLevel[0.5`]},{Thickness[0.01`]}},Axes->False,BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11}]

showit=Show[{points,PLT}]

(* Output to a postscript file *)
(* Export["splint1a.eps",showit]; *)