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 (and the problems with it).
Figure 8.6 was generated with this notebook.

Input: 
Knot values u_i and function values there: (u_i,f(u_i))_i i=1,L ... Here we have selected L=20 values from f(u) = sqrt(u) + sin(u)/10. 

Once we have the data, we need to choose:
Degree n for the approximation. (Here we are using degre 11.)


Output: the least squares approximating polynomial to the input values using the monomial basis and then using the Bernstein basis. 


In the first cell below, we see that the degree 11 monomial polynomial approximates fairly well near the start of the approximation, but it deviates greatly from the input funtion f at the end of the approximation. This is typical behavior for the monomial basis.

In the second cell, the degree 11 polynomial from the Bernstein basis looks like a great fit over the entire knot interval. As stated in the text, the condition numbers of V^T*V -- normal equation matrix -- is 10^33 for the monomial basis and 10^6 for the Bernstein basis.


See Chapter 8: Data Fitting for details.

Updated July 2008

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

(* ==== This section has the routines we will be using ====== *)

(* Create a monomial design matrix  -- see Equation (8.4) in the text *)
(* Rows of the matrix are formed by taking powers of the parameters u[i] *) 
(* These are the monomial polynomials evaluated at the input knot values *)

Vander[n_,L_,u_]:=Table[((u[[i]])^(j-1)),{i,1,L},{j,1,n+1}];


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

(* Define a degree n polynomial in monomial form with coefficients A[j] *)
 poly[n_,A_,x_]:=(∑_(j=1)^(n+1))A[[j]] (x^(j-1)); 

(* Define a matrix  created from V^T V -- see equation (8.5) *)
(* This is the matrix part of the normal equations *)
Mat[n_,L_,u_]:=Transpose[Vander[n,L,u]].Vander[n,L,u];

(* ==== End of the routines section *)

(* Create data by evaluating the function at the given, L,  knot values u_i *)
data[L_]:=Table[f[u[[i]]],{i,1,L}];

(* Find the coefficients of the polynomial by solving the normal equations *)
(* This is solving equation (8.5) *)
Coeffs[n_,L_,u_]:=Inverse[Mat[n,L,u]].(Transpose[Vander[n,L,u]].data[L]);

(* ==== Below  we construct our data and solve the problem ==== *)

(* Create 20 data values *)
L=20;
(* Degree 11 polynomial *)
n=11; 

(* Create uniform knot values u_i *)
u=Table[i,{i,1,L}]

(* Create the degree 11 approximation *)
(* The result is a set of coefficients defining the degree 11 least squares polynomial *)
(* We are solving the normal equations *)
 A=N[Coeffs[n,L,u]];

(* Plot the function and the degree 11 polynomial approximation *)
(* The plot is over x=1 to x=20 *)
PLT=Plot[{f[x],poly[n,A,x]},{x,u[[1]],u[[L]]},PlotStyle->{{Thickness[0.01]},{Thickness[0.01]}},BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11}, PlotRange->{{0,21},Automatic}]

 Export["sqrt_approx11.eps",PLT]

Using the Bernstein basis functions rather than the monomial basis functions.
Much better results!

(* Recursive definition of the Bernstein basis functions *)
Bern[n_,0,t_]:=(((1-t))^n);
Bern[n_,n_,t_]:=(t^n);
Bern[n_,i_,t_]:=Binomial[n,i] (((1-t))^(n-i)) (t^i);

(* Evaluate a set of Bernstein degree n polynomials evaluated at a given t-value *)

Bernvec[n_,t_]:=Table[Bern[n,i,t],{i,0,n}];

(* A Bezier polynomial is formed by multiplying the coefficients by the basis functions *)
 Bez[b_,n_,t_]:=Bernvec[n,t].b;

(* 20 knot values *)
 L=20;

(* Degree 11 *)
n=11;

(* Create uniform knot values *)
u=Table[i,{i,1,L}];

(* Form the design matrix using the Bernstein basis *)
(* See just after eqn (8.7) -- we form a row by evaluating each Bernstein function at a given knot vale *)
(* For Bernstein polynomials, the knots are assumed to be in [0,1] *)
Vander[n_,L_,u_]:=Table[Bern[n,j-1,((u[[i]]-u[[1]])/(u[[L]]-u[[1]]))],{i,1,L},{j,1,n+1}];


(* Evaluate the function -- again with the knots scaled to [0,1] *)
f[x_,u_]:=(Sqrt[((x-u[[1]])/(u[[L]]-u[[1]]))])+Sin[((x-u[[1]])/(u[[L]]-u[[1]]))]; 

(* Form the matrix part of the normal equations *)
 Mat[n_,L_,u_]:=Transpose[Vander[n,L,u]].Vander[n,L,u]; 

(* Create data at the knot values by evaluating our given function *)
data[L_]:=Table[f[u[[i]],u],{i,1,L}];

(* Create the coefficients of the approximating polynomial by solving the normal equations *)
Coeffs[n_,L_,u_]:=N[Inverse[Mat[n,L,u]].(Transpose[Vander[n,L,u]].data[L])];

(* Load the coefficients *)
A=N[Coeffs[n,L,u]];

PLT=Plot[{f[x,u],Bez[A,n,((x-u[[1]])/(u[[L]]-u[[1]]))]},{x,u[[1]],u[[L]]},PlotStyle->{{Thickness[0.02]},{Thickness[0.01]}},BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11}] 


(* Plot a postscript file *)
Export["sqrt_approx_bez_11.eps",PLT] ;