CSE 477 / G. Farin
Interpolation

This explores two iterative approaches. A set of knots is prescribed and then a set of data points is defined with the same number as elements. The interpolation problem is solved once we find a Bezier polygon bez such that the corresponding curve passes through the data points.

n=6;
iter=5n;
knot=Table[t,{t,0.0,1.0,1/n}]; (* uniform knots *)
data=Table[{t^1.5,3Sin[6t]},{t,0.0,10.0,10/n}];

curve[t_]:=Sum[BernsteinBasis[n,i,t]*bez[[i+1]],{i,0,n}];
curvej[j_,t_]:=curve[t]-BernsteinBasis[n,j,t]*bez[[j]]; (*curvej: leave out j-th term *)

center=Sum[data[[i]],{i,1,n+1}]/(n+1);
dists=Table[Norm[data[[i]]-center],{i,1n+1}];
radius=Max[dists];

(* PIA *)
bez=data;
pia1=TimeUsed[];Do[(errors=Table[curve[knot[[j]]]-data[[j]],{j,1,n+1}];
Do[bez[[j]]=bez[[j]]-errors[[j]],{j,1,n+1}]),
{k,1,iter}];
piatime=TimeUsed[]-pia1;

errnorms=Table[Norm[curve[knot[[j]]]-data[[j]]],{j,1,n+1}];
maxerr=Max[errnorms];


dplot=Graphics[{Green,PointSize[0.015],Point[data]}];
points=Graphics[{PointSize[0.015],Point[bez]}];
bplot=ListLinePlot[bez];

cplot=ParametricPlot[curve[t],{t,0,1},PlotStyle->{Thick,Red}];
Show[{bplot,cplot,dplot,points}]
Print["max error PIA:  ",maxerr, "  time: ",piatime]

(*********************************)
(* Gauss-Seidel: *)
bez=data; (*initial guess*)
time1=TimeUsed[];
Do[Do[bez[[j]]= (data[[j]]-curvej[j,knot[[j]]])/BernsteinBasis[n,j,knot[[j]]],{j,2,n}],
{k,1,iter}]
gstime=TimeUsed[]-time1;
errnorms=Table[Norm[curve[knot[[j]]]-data[[j]]],{j,1,n+1}];
maxerr=Max[errnorms];


dplot=Graphics[{Green,PointSize[0.015],Point[data]}];
points=Graphics[{PointSize[0.015],Point[bez]}];
bplot=ListLinePlot[bez];

cplot=ParametricPlot[curve[t],{t,0,1},PlotStyle->{Thick,Red}];
Show[{bplot,dplot,points,cplot}]
Print["max error Gauss-Seidel:  ",maxerr, " time: ",gstime]

(*********************************)
(* Vandermonde: *)

time3=TimeUsed[];
vandermonde=Table[BernsteinBasis[n,i-1,knot[[j]]],{j,1,n+1},{i,1,n+1}];
bez=LinearSolve[vandermonde,data];
vantime=TimeUsed[]-time3;

errnorms=Table[Norm[curve[knot[[j]]]-data[[j]]],{j,1,n+1}];
maxerr=Max[errnorms];


dplot=Graphics[{Green,PointSize[0.015],Point[data]}];
points=Graphics[{PointSize[0.015],Point[bez]}];
bplot=ListLinePlot[bez];

cplot=ParametricPlot[curve[t],{t,0,1},PlotStyle->{Thick,Red}];
Show[{bplot,dplot,points,cplot}]
Print["max error Vandermonde:  ",maxerr, " time: ",vantime]