Mat360 Homework

Tea Time Numerical Analysis
Section 3.2: Lagrange Polynomials

Andy Long, Spring, 2024

pp. 121-122, #1, 2 (find two polynomials!), 8, 9, 10, 13, 22

Preliminaries: 
Lagrange Polynomials make fitting points easy!

Lagrange’s linear functions are linear functions which are zero at one value of x, and one at the other:

lagrange[x_,x1_,x2_]:=(x-x2)/(x1-x2) (* 1 at x1, 0 at x2 *)
{x0,x1}={50,75}
{y0,y1}={180,230}
Show[
Plot[{lagrange[x,x0,x1],lagrange[x,x1,x0]},{x,x0,x1}],
ListPlot[
{{x0,0},{x1,0},{x0,1},{x1,1}},
PlotStyle->{PointSize->Large}]
]

Lagrange linear polynomials are just sums of two Lagrange linear functions that fit two points

Let’s use these Lagrange linear functions  fit two points: (50,180), (75, 230):

L1[x_,x0_,x1_,y0_,y1_]:=y0 lagrange[x,x0,x1]+y1 lagrange[x,x1,x0]
p1=ListPlot[{{x0,y0},{x1,y1}},PlotStyle->{PointSize->Large}];
p2=Plot[L1[x,x0,x1,y0,y1],{x,x0,x1}];
Show[p1,p2]

Exercise #1:

Lagrange quadratic functions are just products of Lagrange linear functions -- and they fit three points: two with function values zeros, and a function value of 1:

lagrange2[x_,x1_,x2_,x3_]:=lagrange[x,x1,x2]*lagrange[x,x1,x3] (* lagrange2 is 1 at x1, zero at x2 and x3 *)
(* (1,2),(1.5,−0.83),and (2.11,−1) For problem #1: *)
{x0,x1,x2}={1,1.5,2.11}
{y0,y1,y2}={2,-0.83,-1}
Show[
Plot[{lagrange2[x,x0,x1,x2],lagrange2[x,x1,x0,x2],lagrange2[x,x2,x0,x1]},{x,x0,x2}],
ListPlot[
{{x0,0},{x1,0},{x2,0},{x0,1},{x1,1},{x2,1}},
PlotStyle->{PointSize->Large}]
]

Lagrange quadratic polynomials are just sums of three Lagrange quadratic functions that each are zero at two points, and 1 at the other

Write down the Lagrange interpolating polynomial passing through (1, 2), (1.5, −0.83), and (2.11, −1).

L2[x_,x0_,x1_,x2_,y0_,y1_,y2_]:=y0 lagrange2[x,x0,x1,x2]+
y1 lagrange2[x,x1,x0,x2]+
y2 lagrange2[x,x2,x0,x1]
L2[x,x0,x1,x2]
p1=ListPlot[{{x0,y0},{x1,y1},{x2,y2}},PlotStyle->{PointSize->Large}];
p2=Plot[L2[x,x0,x1,x2,y0,y1,y2],{x,x0,x2}];
Show[p2,p1]

Exercise #2:

Lagrange cubic functions are also just products of Lagrange linear functions -- and they fit four points: three zeros and a 1 at the fourth point.

lagrange3[x_,x1_,x2_,x3_,x4_]:=lagrange[x,x1,x2]*lagrange[x,x1,x3]*lagrange[x,x1,x4]
(* (0,0),(1,2),(4,−3),and (10,−1) *)
{x0,x1,x2,x3}={0,1,4,10}
{y0,y1,y2,y3}={0,2,-3,-1}
Show[
Plot[{lagrange3[x,x0,x1,x2,x3],lagrange3[x,x1,x0,x2,x3],lagrange3[x,x2,x0,x1,x3],lagrange3[x,x3,x0,x1,x2]},{x,x0,x1}],
ListPlot[
{{x0,0},{x1,0},{x2,0},{x3,0},{x0,1},{x1,1},{x2,1},{x3,1}},
PlotStyle->{PointSize->Large}]
]

Lagrange cubic polynomials are just sums of four Lagrange cubic functions that each are zero at three points, and 1 at the other.

Let’s fit four points of Exercise 2: (0, 0), (1, 2), (4, −3), and (10, −1)

L3[x_,x0_,x1_,x2_,x3_,y0_,y1_,y2_,y3_]:=y0 lagrange3[x,x0,x1,x2,x3]+y1 lagrange3[x,x1,x2,x0,x3]+y2 lagrange3[x,x2,x0,x1,x3]+y3 lagrange3[x,x3,x2,x0,x1]
L3[x,x0,x1,x2,x3,y0,y1,y2,y3]
p1=ListPlot[{{x0,y0},{x1,y1},{x2,y2},{x3,y3}},PlotStyle->{PointSize->Large}];
p2=Plot[L3[x,x0,x1,x2,x3,y0,y1,y2,y3],{x,x0,x3}];
Show[p2,p1]

Exercise 8

a=0.43;
{x0,x1,x2,x3}={0,0.25,0.5,0.75};
{y0,y1,y2,y3}={1,1.64872,2.71828,4.48169};
fa1=L1[a,x1,x2,y1,y2]
fa2=L2[a,x1,x2,x3,y1,y2,y3]
fa3=L3[a,x0,x1,x2,x3,y0,y1,y2,y3]
p1=ListPlot[{{x0,y0},{x1,y1},{x2,y2},{x3,y3}},PlotStyle->{PointSize->Large}];
p1a=ListPlot[{{a,fa1}},PlotStyle->{Blue,PointSize->Large}];
p1b=ListPlot[{{a,fa2}},PlotStyle->{Red,PointSize->Large}];
p1c=ListPlot[{{a,fa3}},PlotStyle->{Yellow,PointSize->Large}];
p2=Plot[L3[x,x0,x1,x2,x3,y0,y1,y2,y3],{x,x0,x3}];
Show[p2,p1,p1a,p1b]

Let’s try to be smart about this and write a function to do each:

Ex8[a_,xpoints_,ypoints_,title_]:=
Module[{
x0=xpoints[[1]],x1=xpoints[[2]],x2=xpoints[[3]],x3=xpoints[[4]],
y0=ypoints[[1]],y1=ypoints[[2]],y2=ypoints[[3]],y3=ypoints[[4]],
fa1,fa2,fa3,p1,p1a,p1b,p2,p3,p4
},
fa1=L1[a,x0,x1,y0,y1];
fa2=L2[a,x0,x1,x2,y0,y1,y2];
fa3=L3[a,x0,x1,x2,x3,y0,y1,y2,y3];
Print[{fa1,fa2,fa3}];
p1=ListPlot[{{x0,y0},{x1,y1},{x2,y2},{x3,y3}},PlotStyle->{PointSize->Large}];
p1a=ListPlot[{{a,fa1}},PlotStyle->{Blue,PointSize->Large}];
p1b=ListPlot[{{a,fa2}},PlotStyle->{Red,PointSize->Large}];
p2=Plot[L1[x,x0,x1,y0,y1],{x,Min[xpoints],Max[xpoints]},PlotStyle->Blue];
p3=Plot[L2[x,x0,x1,x2,y0,y1,y2],{x,Min[xpoints],Max[xpoints]},PlotStyle->Red];
p4=Plot[L3[x,x0,x1,x2,x3,y0,y1,y2,y3],{x,Min[xpoints],Max[xpoints]},PlotStyle->Black,PlotLabel->title];
Show[p4,p3,p2,p1,p1a,p1b]
]

a=0.43;
xpoints={0.25,0.5,0.75,0};
ypoints={1.64872,2.71828,4.48169,1};
Ex8[a,xpoints,ypoints,"8(a)"]

a=0.18;
xpoints={0.1,0.2,0.3,0.4};
ypoints={-0.29004986,-0.56079734,-0.81401972,-1.0526302};
Ex8[a,xpoints,ypoints,"8(b)"]

a=2.26;
xpoints={2,2.5,1.5,1};
ypoints={−1.329,1.776,−2.569,1.654};
Ex8[a,xpoints,ypoints,"8(c)"]

a=11.26;
xpoints={11,12,10,13};
ypoints={−1.2352,−0.8765,−0.7865,0.0021};
Ex8[a,xpoints,ypoints,"8(d)"]

Exercise 9

Ex9[a_,true_,xpoints_,ypoints_,title_]:=
Module[{
x0=xpoints[[1]],x1=xpoints[[2]],x2=xpoints[[3]],
y0=ypoints[[1]],y1=ypoints[[2]],y2=ypoints[[3]],
fa1,fa2,fa3,p1,p1a,p1b,p2,p3,p4
},
fa1=L1[a,x0,x1,y0,y1];
Print[Simplify[L1[x,x0,x1,y0,y1]]];
fa2=L2[a,x0,x1,x2,y0,y1,y2];
Print[Simplify[L2[x,x0,x1,x2,y0,y1,y2]]];
Print[{fa1,fa2}];
Print[Abs[true-{fa1,fa2}]];
p1=ListPlot[{{x0,y0},{x1,y1},{x2,y2}},PlotStyle->{PointSize->Large}];
p1a=ListPlot[{{a,fa1}},PlotStyle->{Blue,PointSize->Large}];
p1b=ListPlot[{{a,fa2}},PlotStyle->{Red,PointSize->Large}];
p2=Plot[L1[x,x0,x1,y0,y1],{x,Min[xpoints],Max[xpoints]},PlotStyle->Blue];
p3=Plot[L2[x,x0,x1,x2,y0,y1,y2],{x,Min[xpoints],Max[xpoints]},PlotStyle->Red,PlotLabel->title];
Show[p3,p2,p1,p1a,p1b,PlotRange->All]
]

a=1.4;
xpoints={1.25,1.6,1};
f[x_]:=Sin[Pi x]
ypoints=f[xpoints];
Ex9[a,f[a],xpoints,ypoints,"9(a)"]

f[x_]:=CubeRoot[x-1]
ypoints=f[xpoints];
Ex9[a,f[a],xpoints,ypoints,"9(b)"]

f[x_]:=E^(2x-4)
ypoints=f[xpoints];
Ex9[a,f[a],xpoints,ypoints,"9(c)"]

f[x_]:=Log[10x]
ypoints=f[xpoints];
Ex9[a,f[a],xpoints,ypoints,"9(d)"]

Exercise 10

prodlist[a_,list_]:=Times@@(a-list)
prodlist[3,{3.1,6,4}]

Ex10[a_,f_,b_,c_,xpoints_,ypoints_,title_]:=
Module[{
x0=xpoints[[1]],x1=xpoints[[2]],x2=xpoints[[3]],
y0=ypoints[[1]],y1=ypoints[[2]],y2=ypoints[[3]],
fa1,fa2,fa3,p1,p1a,p1b,p2,p3,p4,
fpp=f''[a],fppp=f'''[a],product=prodlist[a,xpoints]
},
fa1=L1[a,x0,x1,y0,y1];
fa2=L2[a,x0,x1,x2,y0,y1,y2];
Print[{fa1,fa2}];
Print[Abs[f[a]-{fa1,fa2}]];
Print[Abs[{f''[b],f'''[c]/3*(a-x2)}/2!*(a-x0)*(a-x1)]];
p1=ListPlot[{{x0,y0},{x1,y1},{x2,y2}},PlotStyle->{PointSize->Large}];
p1a=ListPlot[{{a,fa1}},PlotStyle->{Blue,PointSize->Large}];
p1b=ListPlot[{{a,fa2}},PlotStyle->{Red,PointSize->Large}];
p2=Plot[L1[x,x0,x1,y0,y1],{x,Min[xpoints],Max[xpoints]},PlotStyle->Blue];
p3=Plot[L2[x,x0,x1,x2,y0,y1,y2],{x,Min[xpoints],Max[xpoints]},PlotStyle->Red];
p4=Plot[f[x],{x,Min[xpoints],Max[xpoints]},PlotStyle->Black,PlotLabel->title];
Show[p4,p3,p2,p1,p1a,p1b,PlotRange->All]
]

a=1.4;
xpoints={1.25,1.6,1};
f[x_]:=Sin[Pi x]
ypoints=f[xpoints];
Plot[{f''[x],f'''[x]},{x,Min[xpoints],Max[xpoints]},PlotStyle->{Blue,Black},PlotLabel->"Derivatives of f(x)"]
f[a]
b=3/2;
c=1;
Ex10[a,f,b,c,xpoints,ypoints,"10(a)"]

a=1.4;
xpoints={1.25,1.6,1};
f[x_]:=CubeRoot[x-1]
ypoints=f[xpoints];
Plot[{f''[x],f'''[x]},{x,Min[xpoints],Max[xpoints]},PlotStyle->{Blue,Black},PlotLabel->"Derivatives of f(x)"]
f[a]
b=1.25;
c=1;
Ex10[a,f,b,c,xpoints,ypoints,"10(b)"]

a=1.4;
xpoints={1.25,1.6,1};
f[x_]:=E^(2x-4)
ypoints=f[xpoints];
Plot[{f''[x],f'''[x]},{x,Min[xpoints],Max[xpoints]},PlotStyle->{Blue,Black},PlotLabel->"Derivatives of f(x)"]
f[a]
b=1.6;
c=1.6;
Ex10[a,f,b,c,xpoints,ypoints,"10(c)"]

a=1.4;
xpoints={1.25,1.6,1};
f[x_]:=Log[10x]
ypoints=f[xpoints];
Plot[{f''[x],f'''[x]},{x,Min[xpoints],Max[xpoints]},PlotStyle->{Blue,Black},PlotLabel->"Derivatives of f(x)"]
f[a]
b=1.25;
c=1;
Ex10[a,f,b,c,xpoints,ypoints,"10(d)"]

Exercise 13

Recursive definition of Interpolating polynomials: Neville’s method

Neville’s method is sort of the ideal limit of one of my favorite principles: that whenever you have two estimates for a quantity, you have a third: the average. In this case, the average of successive averages actually gives you the exact value at the end! We just slowly work up to the answer.

Each successive p[[i,j]] is a weighted average estimate of f(x), using the linear polynomial through two points.

nevillesMethod[x_,data_,verbose_]:=
Module[{n=Length[data],xs,ys,i,j,p},
p=ConstantArray[0,{n,n}];
{xs,ys}=Transpose[data];  (* The data is a list of coordinates *)
p[[All,1]]=ys; (* We initial the first column of p with the y-values *)
For[j=2,j<=n+1,j++,
For[i=1,i<=n+1-j,i++, 
p[[i,j]]= (* Success columns are filled with interpolated value, 
weighted linear interpolators of averaged y-values from the previous column. *)
((x-xs[[i]])p[[i+1,j-1]]-(x-xs[[i+j-1]])p[[i,j-1]])
/
(xs[[i+j-1]]-xs[[i]]);
];
];
If[verbose,Print[MatrixForm[p]]];
p[[1,n]]
]

Table 3.2, p. 119: Example of Neville’s method

a=1.5;
fapprox=nevillesMethod[a,{{0,1},{1,E},{2,E^2}},True]
fakeExp[x_]:=nevillesMethod[x,{{0,1},{1,E},{2,E^2}},False]
Show[Plot[{E^x,fakeExp[x]},{x,0,2},PlotLabel->"Fake e^x with a quadratic"],ListPlot[{{a,fapprox}},PlotStyle->{PointSize->Large}]
]

a=1.5;
f[x_]:=2^(x/2)
xpoints={0,1,2,3};
fapprox=nevillesMethod[a,Transpose[{xpoints,f[xpoints]}],True]
Plot[{f''''[x]},{x,Min[xpoints],Max[xpoints]},PlotStyle->{Blue,Black}]
f''''[3]/4!*prodlist[a,xpoints]
Abs[f[a]-fapprox]

Exercise 22

f[x_]:=E^(-x)

The absolute value of any derivative is going to be maximized at x=0, and equals 1:

Plot[f''[x],{x,0,1},PlotLabel->"All the derivative bounds look like this!"]

Picking any two points on the interval [0,1], we examine the function (x-x0)(x-x1). If both were chosen to be the same endpoint (e.g. x0=x1=0), then the function is x^2, and is maximized at 1, with value 1. That is the maximum possible.

Manipulate[
Show[
Plot[(x-x0)(x-x1),{x,0,1},PlotLabel->"Maximizing the quadratic part"],
ListPlot[{{x0,0},{x1,0}},PlotStyle->{PointSize->Large}]
],
{x0,0,x1},
{x1,x0,1}
]

Therefore the error bound is (1)/2!*1=1/2

1/2.0

Worst case scenarios involve using a tangent at either endpoint, and the absolute errors are within our bound (comfortably):

Plot[{f[x],1-x,E^(-1)-E^(-1)(x-1)},{x,0,1},PlotLabel->"with worst case scenarios...."]