ECE 3340 - Numerical Fitting, Interpolation, Spline

Part B - Interpolation and Spline

1. Introduction

1.1 Example of polynomial fit

code

Manipulate[
If[!TrueQ[norderold==norder],
fitfunc=Fit[data1,Table[u^m,{m,0,norder}],u] ;
norderold=norder
,] ;
cursor=Graphics[{Dashed,Hue[0.85,1,1],Thickness[0.005],Line[{{x,-3.5},{x,5}}]}];
interpPoint=Graphics[{Hue[0.1,1.,1.],PointSize[0.04],Point[{x,fitfunc /. u->x}]}];
Show[ListLinePlot[data1,ImageSize->400,PlotStyle->{{RGBColor[1.,0.,0.],Thickness[0.01]},{RGBColor[0.,0.,1.],Thickness[0.01]}},Background->Hue[0.5,0.2,0.95],Frame->True]
,If[TrueQ[POLTFIT],{Plot[fitfunc,{u,0.,2.},PlotStyle-> {Hue[0.36,1,0.6],Thickness[0.007]}],interpPoint},Graphics[Point[{0.,0.}]]],cursor
,Graphics[{Blue,PointSize[0.025],Point /@ data1}]]
,{{x,0.5},0.,2},Row[{Control[{{norder,3},Range[Length[data1]-1],ControlType->SetterBar}]
,Spacer[20],Button[Dynamic[If[TrueQ[POLTFIT],,]]
,{ If[TrueQ[POLTFIT],POLTFIT=False,POLTFIT=True]}]}]
,Initialization:>{POLTFIT=False;
data1={{0.,-3.12},{0.2,0.44},{0.4,3.12},{0.6,3.2},{0.8,2.48},{1.,1.35},{1.2,0.49},{1.4,1.88},{1.6,2.61},{1.8,4.92}};}]

What value should we use for, say x=0.5? Or any x?

1.2 Basic ideas

When we try to estimate a value for a function at a point that is WITHIN the range of a the given points, we call that interpolation. When we try to estimate a value for a function at a point that is OUTSIDE the range, we call that extrapolation.

Question: which one, inter- or extra-polation is "safer"?

Here, we will learn 2 polynomial interpolation techniques: Newton & Lagrange.

1.3 Issues with interpolation

Key considerations when interpolating:
1. How accurate are the input points (data or calculation)?
2. Should "fitting" (interpolation or extrapolation) create excessive information beyond the input?

Analytical exact fit can be done with polynomial
y=(Subscript[a, 0])+(Subscript[a, 1]) x+(Subscript[a, 2]) (x^2)+...+(Subscript[a, n-1]) (x^(n-1))+(Subscript[a, n]) (x^n)

Do we really believe that the orange point represents a credible interpolated value?

2. Newton's and Lagrange’s polynomials

2.1 Newton’s polynomials (divided-difference)

2.1.1 Linear

f[x]=f[(Subscript[x, 0])] + (x-(Subscript[x, 0])) ((f[(Subscript[x, 1])]- f[(Subscript[x, 0])]     )/(((Subscript[x, 1])-(Subscript[x, 0]))))
f[(Subscript[x, 0])=1] =2.2 ; f[(Subscript[x, 1])=3.5] =5.1

x0=1 ; x1=3.5 ; (Subscript[f, 0]) = 2.2 ; (Subscript[f, 1]) =5.1 ;
f[x_]:=(Subscript[f, 0]) +(x-x0) (((Subscript[f, 1])-(Subscript[f, 0]))/(x1-x0))

aplot=Plot[f[x],{x,1.,3.5},   ImageSize->{300,300}
     , PlotStyle->{RGBColor[1.,0.,0.],Thickness[0.02]}
            
            ,Background->RGBColor[0.7,0.9,0.9]
            ,Frame->True
                
     ] ;
bplot=ListPlot[{{x0,(Subscript[f, 0])},{x1,(Subscript[f, 1])}}, PlotStyle->AbsolutePointSize[10]] ;
Show[aplot,bplot]

2.1.1 Quadratic fit

x0=1 ; x1=3.5 ; x2= 5.; (Subscript[f, 0]) = 2.2 ; (Subscript[f, 1]) =5.1 ; (Subscript[f, 2])=4.1 ;
bplot2=ListPlot[{{x0,(Subscript[f, 0])},{x1,(Subscript[f, 1])},{x2,(Subscript[f, 2])}}, PlotStyle->AbsolutePointSize[10]]

f(x)=(Subscript[b, 0]) + (Subscript[b, 1])  (x-(Subscript[x, 0]))  + (Subscript[b, 2])  (x-(Subscript[x, 0]))  (x-(Subscript[x, 1])) 
We have to determine the 3 unknown b's
f(x= (Subscript[x, 0])) = (Subscript[b, 0]) ;
f(x= (Subscript[x, 1])) = (Subscript[b, 0]) +(Subscript[b, 1])  ((Subscript[x, 1])-(Subscript[x, 0])) ;

ClearAll[f, b] ;
Solve[  (Subscript[b, 0]) + (Subscript[b, 1])*((Subscript[x, 1])-(Subscript[x, 0]))==f[(Subscript[x, 1])],{(Subscript[b, 1])}]

For the third term
f(x= (Subscript[x, 2])) = (Subscript[b, 0]) +(Subscript[b, 1])  ((Subscript[x, 2])-(Subscript[x, 0])) +(Subscript[b, 2])  ((Subscript[x, 2])-(Subscript[x, 0])) ((Subscript[x, 2])-(Subscript[x, 1])) ;
f( (Subscript[x, 2])) = f[(Subscript[x, 0])] +((f((Subscript[x, 1]))-f((Subscript[x, 0])))/((Subscript[x, 1])-(Subscript[x, 0]))) ((Subscript[x, 2])-(Subscript[x, 0])) +(Subscript[b, 2])  ((Subscript[x, 2])-(Subscript[x, 0])) ((Subscript[x, 2])-(Subscript[x, 1])) ;

Solve[ f[(Subscript[x, 0])] + ((f((Subscript[x, 1]))-f((Subscript[x, 0])))/((Subscript[x, 1])-(Subscript[x, 0])))((Subscript[x, 2])-(Subscript[x, 0])) + (Subscript[b, 2]) ((Subscript[x, 2])-(Subscript[x, 0]))((Subscript[x, 2])-(Subscript[x, 1]))==f[(Subscript[x, 2])]
	              ,{(Subscript[b, 2])}]

{{(Subscript[b, 2])->-((f((Subscript[x, 0]))-f((Subscript[x, 2]))+(((f((Subscript[x, 1]))-f((Subscript[x, 0]))) (((Subscript[x, 2])-(Subscript[x, 1]))+((Subscript[x, 1])-(Subscript[x, 0]))))/((Subscript[x, 1])-(Subscript[x, 0]))))/(((Subscript[x, 2])-(Subscript[x, 0])) ((Subscript[x, 2])-(Subscript[x, 1]))))}}

{{(Subscript[b, 2])->-(f((Subscript[x, 0]))-f((Subscript[x, 2]))+(((f((Subscript[x, 1]))-f((Subscript[x, 0]))) ((Subscript[x, 2])-(Subscript[x, 1])))/((Subscript[x, 1])-(Subscript[x, 0])))+(f((Subscript[x, 1]))-f((Subscript[x, 0]))))/(((Subscript[x, 2])-(Subscript[x, 0])) ((Subscript[x, 2])-(Subscript[x, 1])))}}

{{(Subscript[b, 2])->-(f((Subscript[x, 0]))-f((Subscript[x, 2]))+(((f((Subscript[x, 1]))-f((Subscript[x, 0]))) ((Subscript[x, 2])-(Subscript[x, 1])))/((Subscript[x, 1])-(Subscript[x, 0])))+(f((Subscript[x, 1]))-f((Subscript[x, 0]))))/(((Subscript[x, 2])-(Subscript[x, 0])) ((Subscript[x, 2])-(Subscript[x, 1])))}}

(Subscript[b, 2])->(((((f((Subscript[x, 2]))-f((Subscript[x, 1]))) )/((Subscript[x, 2])-(Subscript[x, 1])))-(((f((Subscript[x, 1]))-f((Subscript[x, 0]))) )/((Subscript[x, 1])-(Subscript[x, 0]))))/(((Subscript[x, 2])-(Subscript[x, 0])) ))

x0=1 ; x1=3.5 ; x2= 5.; (Subscript[f, 0]) = 2.2 ; (Subscript[f, 1]) =5.1 ; (Subscript[f, 2])=4.1 ;
f[x_]:=(Subscript[f, 0]) +(((Subscript[f, 1])-(Subscript[f, 0]))/(x1-x0))(x-x0) +((((((Subscript[f, 2])-(Subscript[f, 1])) )/(x2-x1))-((((Subscript[f, 1])-(Subscript[f, 0])) )/(x1-x0)))/((x2-x0) )) (x-x0) (x-x1)

aplot2=Plot[f[x],{x,1.,5},   ImageSize->{300,300}
     , PlotStyle->{RGBColor[1.,0.,0.],Thickness[0.01]}  ,Background->RGBColor[0.7,0.9,0.9]
            ,Frame->True] ;
Show[aplot2,bplot2]

2.1.3 Higher-order fit

For higher-order, the general format is
f(x)=(Subscript[b, 0]) + (Subscript[b, 1])  (x-(Subscript[x, 0]))  + (Subscript[b, 2])  (x-(Subscript[x, 0]))  (x-(Subscript[x, 1])) + ...
                                         (Subscript[b, n])  (x-(Subscript[x, 0]))  (x-(Subscript[x, 1])) ...  (x-(Subscript[x, n-1]))

Now we can generalize the formalism from the above:
(Subscript[b, 1]) = ((f((Subscript[x, 1]))-f((Subscript[x, 0])))/((Subscript[x, 1])-(Subscript[x, 0]))) ;

(Subscript[b, 2])= (((((f((Subscript[x, 2]))-f((Subscript[x, 1]))) )/((Subscript[x, 2])-(Subscript[x, 1])))-(((f((Subscript[x, 1]))-f((Subscript[x, 0]))) )/((Subscript[x, 1])-(Subscript[x, 0]))))/(((Subscript[x, 2])-(Subscript[x, 0])) ))= ((fd[(Subscript[x, 2]);(Subscript[x, 1])] - fd[(Subscript[x, 1]);(Subscript[x, 0])])/(((Subscript[x, 2])-(Subscript[x, 0]))))=fd[(Subscript[x, 2]);(Subscript[x, 1]); (Subscript[x, 0])] ;

(Subscript[b, 3]) = ((fd[(Subscript[x, 3]);(Subscript[x, 2]); (Subscript[x, 1])] - fd[(Subscript[x, 2]);(Subscript[x, 1]); (Subscript[x, 0])])/(((Subscript[x, 3])-(Subscript[x, 0]))))

2.2 Newton’s method: exercise in class

Use Newton's interpolation polynomials to approximate Sin[x] to the third-order for x=56 degrees...

Question 1: How many point we need?

Answer: For third order, we need 4 points: f[(Subscript[x, 0])], f[(Subscript[x, 1])], f[(Subscript[x, 2])], f[(Subscript[x, 3])]

Question 2: Write out the full expression for Newton's interpolation

...Answer:

f(x)=(Subscript[b, 0]) + (Subscript[b, 1])  (x-(Subscript[x, 0]))  + (Subscript[b, 2])  (x-(Subscript[x, 0]))  (x-(Subscript[x, 1])) +  (Subscript[b, 3])  (x-(Subscript[x, 0]))  (x-(Subscript[x, 1])) (x-(Subscript[x, 2]))

Question 3: Evaluate all the b coefficients given...

(Subscript[x, 0])=0; Sin[(Subscript[x, 0])]=0;
(Subscript[x, 1])=45 degrees; Sin[(Subscript[x, 1])]=((Sqrt[2])/2);
(Subscript[x, 2])=60 degrees; Sin[(Subscript[x, 2])]=((Sqrt[3])/2);
(Subscript[x, 3])=90 degrees; Sin[(Subscript[x, 3])]=1;

...Answer:

Write out formula for b0

b0=f[x0] /. f->Sin

b0   /. x0->0.

b0=%

Write out formula for b1

b1=((f[x1]-f[x0])/(x1-x0))

b1 /. f->Sin /. x1->(π/4)  /. x0->0

b1=%

Write out formula for b2

b2=((fd[x2,x1]-fd[x1,x0])/(x2-x0))

Write out the expression for fd

fd[u2_,u1_]:=((f[u2]-f[u1])/(u2-u1))

b2

Substitute known x1, x2, etc. for fd

fd[x2,x1]     /.  x1->(π/4) 
			                 /.  x2->(π/3)

((12 (f((π/3))-f((π/4))))/π) /. f->Sin

fd[x1,x0]     /.  x1->(π/4) 
			                 /.  x0->0

((4 (f((π/4))-f(0)))/π) /. f->Sin

Now that we have fd values, substitute them in for b2:
Write the b2 expression again

b2

((((12 (-(1/(Sqrt[2]))+((Sqrt[3])/2)))/π)-((2 (Sqrt[2]))/π))/(x2-x0))  /.  x2->(π/3) /. x0->0

Simplify[%]

b2=-((6 (4 (Sqrt[2])-3 (Sqrt[3])))/(π^2))

What else do we need?

...Answer:  we need b3
Write out the expression for b3:

b3=((fd3[x3,x2,x1]-fd3[x2,x1,x0])/(x3-x0))

fd3[u3_,u2_,u1_]:=((fd2[u3,u2]-fd2[u2,u1])/(u3-u1))

fd2[u2_,u1_]:=((f[u2]-f[u1])/(u2-u1))

b3

((f(x3)-f(x2))/(x3-x2)) /. f->Sin /.  x3->(π/2) /.  x2->(π/3)

((f(x2)-f(x1))/(x2-x1)) /. f->Sin /.  x2->(π/3) /.  x1->(π/4)

((((6 (1-((Sqrt[3])/2)))/π)-((12 (-(1/(Sqrt[2]))+((Sqrt[3])/2)))/π))/(x3-x1)) /.   x1->(π/4)  /.  x3->(π/2)

Simplify[%]

b3=((((12 (2+2 (Sqrt[2])-3 (Sqrt[3])))/(π^2))-b2)/(x3-x0))  /.   x0->0  /.  x3->(π/2)

Simplify[%]

b3=%

Let's write out all the coefficients

b0
b1
b2
b3

Question 4: Write the full Newton's interpolation expression with the coefficients found above

Answer:

b0+b1*(x-x0) + b2*(x-x0) *(x-x1)+b3*(x-x0) *(x-x1)*(x-x2)

b0+b1*(x-x0) + b2*(x-x0) *(x-x1)+b3*(x-x0) *(x-x1)*(x-x2)  /. x0->0 /.  x1->(π/4)  /.  x2->(π/3) /.  x3->(π/2)

Question 5: Do numerical evaluation for x=56 degrees.

f[x_]:=((12 (4+8 (Sqrt[2])-9 (Sqrt[3])) (x-(π/3)) (x-(π/4)) x)/(π^3))-((6 (4 (Sqrt[2])-3 (Sqrt[3])) (x-(π/4)) x)/(π^2))+((2 (Sqrt[2]) x)/π)

f[56.*Pi/180]

Compare with the more accurate answer

Sin[56*Pi/180.]

Relative error:

ScientificForm[ Abs[(0.82879-0.82904)/0.82904] ]

2.3 Newton’s method: Error estimate

From the above example, we see that the higher order the term that we include is, the more "accurate" the interpolation is supposed to be. Thus, an estimate of the error of the nth order is
(Subscript[R, n]) ~ (Subscript[f, n+1])(x) - (Subscript[f, n])(x) 
So,
(Subscript[R, n]) ~ f[(Subscript[x, n+1 ]);(Subscript[x, n ]); ... ; (Subscript[x, o])  ] (x-(Subscript[x, 0]))  (x-(Subscript[x, 1])) ...  (x-(Subscript[x, n]))   
But this is true only if the function is strongly convergence to the interpolation series! We can see that if we put "funny" numbers in, we'll get "funny" results. Thus, polynomial interpolation is a GIGO system (garbage-in, garbage-out). We must use it wisely.
Other interpolation such as Spline, Fourrier, and Wavelet are developed to alleviate these problems

2.4 Lagrange’s polynomials

Objective is the same as Newton, but just different algorithm: Looks for polynomial expressions that fit every point right away:

Therefore:
(Subscript[L, 0])[x]=1 at x=(Subscript[x, 0]); and  = 0 at other (Subscript[x, i])
(Subscript[L, 1])[x]=1 at x=(Subscript[x, 1]); and  = 0 at other (Subscript[x, i])
(Subscript[L, 2])[x]=1 at x=(Subscript[x, __?]); and  = 0 at other (Subscript[x, i])
So what are L[x]?
(Subscript[L, 0])[x]=(((x-(Subscript[x, 1]))(x-(Subscript[x, 2]))(x-(Subscript[x, 3]))...(x-(Subscript[x, n])))/(((Subscript[x, 0])-(Subscript[x, 1]))((Subscript[x, 0])-(Subscript[x, 2]))((Subscript[x, 0])-(Subscript[x, 3]))...((Subscript[x, 0])-(Subscript[x, n]))))
(Subscript[L, 1])[x]=(((x-(Subscript[x, 0]))(x-(Subscript[x, 2]))(x-(Subscript[x, 3]))...(x-(Subscript[x, n])))/(((Subscript[x, 1])-(Subscript[x, 0]))((Subscript[x, 1])-(Subscript[x, 2]))((Subscript[x, 1])-(Subscript[x, 3]))...((Subscript[x, 1])-(Subscript[x, n]))))
(Subscript[L, 2])[x]=(((x-(Subscript[x, __?]))(x-(Subscript[x, __?]))(x-(Subscript[x, __?]))...(x-(Subscript[x, __?])))/(((Subscript[x, __?])-(Subscript[x, 0]))((Subscript[x, __?])-(Subscript[x, 1]))((Subscript[x, __?])-(Subscript[x, __?]))...((Subscript[x, __?])-(Subscript[x, __?]))))

2.5 Lagrange’s method: Example (exercise in class)

Use Lagrange's interpolation polynomials to approximate Sin[x] to the third-order for x=56 degrees...

Question 1: How many points do we need?

Answer: For third order, we need 4 points: f[(Subscript[x, 0])], f[(Subscript[x, 1])], f[(Subscript[x, 2])], f[(Subscript[x, 3])]

Question 2: Write out all the L-polynomials given...

(Subscript[x, 0])=0; Sin[(Subscript[x, 0])]=0;
(Subscript[x, 1])=45 degrees; Sin[(Subscript[x, 1])]=((Sqrt[2])/2);
(Subscript[x, 2])=60 degrees; Sin[(Subscript[x, 2])]=((Sqrt[3])/2);
(Subscript[x, 3])=90 degrees; Sin[(Subscript[x, 3])]=1;

...Answer:

Write out formula for L0[x]

L0[x_]:=(((x-x1)  (x-x2)  (x-x3))/((x0-x1)  (x0-x2)  (x0-x3)))

Substitute known x1, x2, etc..

L0[x_]:=(((x-x1)  (x-x2)  (x-x3))/((x0-x1)  (x0-x2)  (x0-x3))) /. x0->0
                                                                                  /.  x1->(π/4) 
			                                                                       /.  x2->(π/3)
                                                                                   /.  x3->(π/2)

What do you get?

L0[x]

Do the same thing for L1[x]

L1[x_]:=(((x-x0)  (x-x2)  (x-x3))/((x1-x0)  (x1-x2)  (x1-x3)))/. x0->0
                                                                                  /.  x1->(π/4) 
			                                                                       /.  x2->(π/3)
                                                                                   /.  x3->(π/2)

What do you get?

L1[x]

What do you get for L2?

L2[x_]:=(((x-x0)  (x-x1)  (x-x3))/((x2-x0)  (x2-x1)  (x2-x3)))/. x0->0
                                                                                  /.  x1->(π/4) 
			                                                                       /.  x2->(π/3)
                                                                                   /.  x3->(π/2)

What do you get for L3?

L3[x_]:=(((x-x0)  (x-x1)  (x-x2))/((x3-x0)  (x3-x1)  (x3-x2)))/. x0->0
                                                                                  /.  x1->(π/4) 
			                                                                       /.  x2->(π/3)
                                                                                   /.  x3->(π/2)

L3[x]

Question 3: Write the full interpolation expression with the polynomials obtained from question 2, leave sin[x0], sin[x1],... unevaluated

Answer:

f[x0] L0[x] + f[x1] L1[x] + f[x2] L2[x] + f[x3] L3[x]

% /. f->Sin

Question 4: Do numerical evaluation for x=56 degrees.

Compare with the more accurate answer

Sin[56*Pi/180.]

Relative error:

ScientificForm[ Abs[(0.82879-0.82904)/0.82904] ]

Plot out

L0[x_]:=-((24 (x-(π/2)) (x-(π/3)) (x-(π/4)))/(π^3));
L1[x_]:=((192 x (x-(π/2)) (x-(π/3)))/(π^3));
L2[x_]:=-((216 x (x-(π/2)) (x-(π/4)))/(π^3));
L3[x_]:=((48 x (x-(π/3)) (x-(π/4)))/(π^3));

lplot=Plot[{ Sin[0.]*L0[x], Sin[Pi/4.]*L1[x],Sin[Pi/3.]* L2[x],Sin[Pi/2.]* L3[x], Sin[0.]*L0[x]+ Sin[Pi/4.]*L1[x]+Sin[Pi/3.]*L2[x]+Sin[Pi/2.]*L3[x] },{x,0.,Pi/2.}
, ImageSize->500, PlotStyle->{{Hue[0.,1,1.],Thickness[0.01]},
                                                 {Hue[0.3,1,0.5],Thickness[0.01]}
	                               	,   {Hue[0.6,1,1.],Thickness[0.01]}
		                                 ,   {Hue[0.8,1.,1],Thickness[0.01]}
			                              ,   {Hue[0.,0.,1],Thickness[0.01]}}
            ,Background->RGBColor[0.7,0.9,0.9]
            ,Frame->True,
	         GridLines->Automatic ];

pointplot=ListPlot[{{0.,Sin[0.]} ,{Pi/4., Sin[Pi/4.]} ,{Pi/3., Sin[Pi/3.]},{Pi/2., Sin[Pi/2.]}}, PlotStyle->{Red,AbsolutePointSize[15]},AspectRatio->1] ;

Show[lplot,pointplot]

3. Spline method

3.1 Overview discussion

(see ppt file)

3.2 BSpline

3.2.1 m=1

knot1={{0,0,1},{0,1/2,1},{0,1,1}};

Grid[{Table[Plot[Evaluate[BSplineBasis[{1,knot1[[i]]},0,x]],{x,0,1}, Frame-> True, PlotStyle->{ Thickness[0.01]
,Hue[((i-1)/4),1.,0.8]}
,Filling-> {1-> {0.,Hue[((i-1)/4),0.5,1,0.5]}}],{i,Length[knot1]}]}]

3.2.2 m=2

knot2={{0,0,0,1},{0,0,1/2,1},{0,1/2,1/2,1},{0,1/2,1,1}};

Grid[{Table[Plot[Evaluate[BSplineBasis[{2,knot2[[i]]},0,x]],{x,0,1}, Frame-> True, PlotStyle->{ Thickness[0.01]
,Hue[((i-1)/4),1.,0.8]}
,Filling-> {1-> {0.,Hue[((i-1)/4),0.5,1,0.5]}}],{i,Length[knot2]}]}]

3.2.3 m=3

knot3={{0,0,0,0,1/3},{0,0,0,1/3,2/3},
{0,0,1/3,2/3,1},{0,1/2,2/3,1,1}};

Grid[{Table[Plot[Evaluate[BSplineBasis[{3,knot3[[i]]},0,x]],{x,0,1}, Frame-> True, PlotStyle->{ Thickness[0.01]
,Hue[((i-1)/4),1.,0.8]}
,Filling-> {1-> {0.,Hue[((i-1)/4),0.5,1,0.5]}}],{i,Length[knot3]}]}]

Plot[BSplineBasis[{3,{0,1,2,3,4}/4},0,x],{x,0,1}, Frame-> True, PlotStyle->{ Thickness[0.01]
,Hue[((1-1)/4),1.,0.8]}
,Filling-> {1-> {0.,Hue[((1-1)/4),0.5,1,0.5]}}]

3.3 Comparison of various methods

For discussion:
- for many physical quantities, not enough to fit the function but also derivatives.
- spline is intrisically designed to optimize (minimize) derivative-dependent phenomena - example Lorentz lineshape.

3.4 BSpline applications

See example in HW.

See additional spline examples in Mathematica Demonstration webpage.