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

This notebook demonstrates polynomial least squares approximation with B-splines.
Figure 8.12 was generated with this notebook.

Compare the results here with those in leastsq_pol.nb

Input: Knot sequence called 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.  The v_i live within the u-knot vector.

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

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

(* Set a zero tolerance -- anything smaller than this is considered to be zero *)
 eps=5.0^-6;


(* Definition of a B-spline *)
(* It is composed of a left and right span *)

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]]))];


(* Recursive definition of a B-spline *)
(* First procedure is the terminating condition *)
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;




(* Define knots for a B-spline *)
(* The last knot is made a little larger so that we won't get a zero divide when evaluating at the very end of the domain of definition *)
(* In figure 8.12 we add knots 2,6 and then 1 to look at an improve approximation *)
u={0,0,0,0,4,8,12,16,20,20,20,20.00001};

(* Degree = 3 *)
n=3;

(* Where we will evaluate the function f to create data values *) 
v={0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20};

ldata=Length[v]; 

(* The number of parameter intervals and degree determines the number of control points *)
k=Length[u]-3;

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


(* A B-spline is define as a product of b-spline polynomials and coefficients d_i *)
bsfunct[n_,u_,d_,l_,x_]:=(∑_(i=0)^l)d[[i+1]] bspline[n,u,i,x];

(* Create the generalized Vandermonde matrix for a degree 3 B-spline *)
(* The B-spline functions are evaluated at each knot *)
Vander[n_,l_,u_,ldata_,v_]:=Table[bspline[3,u,i,v[[j]]],{j,1,ldata},{i,0,l}]; 



(* Test the Vandermonde matrix *)
(* Vander[n,l,u,ldata,v]  *)

(* Define the function from which to take values *)
f[x_]:=N[(Sqrt[x])+((Sin[x])/10)];

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

(* Create data by evaluating the function at the evaluation points *)
data[ldata_]:=Table[f[v[[j]]],{j,1,ldata}];

(* Create the b-spline coefficients *)
bvals[n_,l_,u_,ldata_,v_]:=Inverse[Mat[n,l,u,ldata,v]].(Transpose[Vander[n,l,u,ldata,v]].data[ldata]); 

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

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

(* Just to check ... print out the b-spline  at a given parameter = 4.4 *)
 (* bsfunct[n,u,d,l,4.4] *)

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


(* Create a postscript plot *)
(Export["bspleast4.eps",PLT];)