Tea Time Numerical Analysis
Section 3.2: Lagrange Polynomials

Andy Long, Spring, 2024

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):

p1=ListPlot[{{x0,y0},{x1,y1}}];
p2=Plot[y0 lagrange[x,x0,x1]+y1 lagrange[x,x1,x0],{x,x0,x1}];
Show[p1,p2]

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 *)
{x0,x1,x2}={50,60,75}
{y0,y1,y2}={180,210,230}
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

Let’s fit three points: (50,180), (75, 190), (60,210)

p1=ListPlot[{{x0,y0},{x1,y1},{x2,y2}}];
p2=Plot[
y0 lagrange2[x,x0,x1,x2]+
y1 lagrange2[x,x1,x0,x2]+
y2 lagrange2[x,x2,x0,x1],
{x,x0,x2}];
Show[p1,p2]

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]
{x0,x1,x2,x3}={50,75,60,70}
{y0,y1,y2,y3}={180,230,210,215}
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: (50,180), (75, 230), (60,210), (70,215)

p1=ListPlot[{{x0,y0},{x1,y1},{x2,y2},{x3,y3}}];
p2=Plot[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],{x,x0,x1}];
Show[p1,p2]

Lagrange cubic polynomials can be useful for design:

Manipulate[
Module[{},
f[x_]:=y0 lagrange3[x,x0,x1,x2,x3]+
y1 lagrange3[x,x1,x0,x2,x3]+
y2 lagrange3[x,x2,x0,x1,x3]+
y3 lagrange3[x,x3,x0,x1,x2];
Show[
Plot[f[x],{x,0,1},PlotRange->{0,200},PlotLabel->Together[f[x]]],
ListPlot[{{x0,y0},{x1,y1},{x2,y2},{x3,y3}},
PlotStyle->{PointSize->Large}]
]
],
{{y0,20},0,200},
{{y1,130},0,200},
{{y2,50},0,200},
{{y3,80},0,200},
{{x0,0},0,x1},
{{x1,1/3},x0,x2},
{{x2,2/3},x1,x3},
{{x3,1},x2,1}
]

But back to computational efficiency: we’re doing a lot more computation than we need to....

fitData[data_]:=
Module[{n=Length[data],trans=Transpose[data],xs,ys,window},
{xs,ys}=trans;
If[n==4,
Module[
{x0=xs[[1]],x1=xs[[2]],x2=xs[[3]],x3=xs[[4]],
y0=ys[[1]],y1=ys[[2]],y2=ys[[3]],y3=ys[[4]]},
f[x_]:=
y0 lagrange3[x,x0,x1,x2,x3]+
y1 lagrange3[x,x1,x0,x2,x3]+
y2 lagrange3[x,x2,x0,x1,x3]+
y3 lagrange3[x,x3,x0,x1,x2];
],
If[n==3,
Module[
{x0=xs[[1]],x1=xs[[2]],x2=xs[[3]],
y0=ys[[1]],y1=ys[[2]],y2=ys[[3]]},
f[x_]:=y0 lagrange2[x,x0,x1,x2]+
y1 lagrange2[x,x1,x0,x2]+
y2 lagrange2[x,x2,x0,x1];
],
If[n==2,
Module[
{x0=xs[[1]],x1=xs[[2]],
y0=ys[[1]],y1=ys[[2]]},
f[x_]:=y0 lagrange[x,x0,x1]+
y1 lagrange[x,x1,x0];
]
]
]
];
window=Max[ys]-Min[ys];
Show[
Plot[f[x],{x,Min[xs],Max[xs]},PlotRange->{Min[ys]+window*.5,Max[ys]-window*.5},PlotLabel->Together[f[x]]],
ListPlot[Transpose[{xs,ys}],PlotStyle->{PointSize->Large}]
]
]
(* #2 *)
data={{1,0.7},{0,0.8},{2,0.75},{3,0.5}};
fitData[data]

#3 and 4, p. 121

By the way, the  underlying function doesn’t need to be differentiable in order for us to make the estimate! Not sure why he added that, but maybe he is setting himself up for a problem later....

fitData[{{1,1},{2,1},{3,2}}]
f[1.3]
fitData[{{0,10},{30,58},{1029,-32}}]
f[1.3]
fitData[{{-10,10},{20,58},{1019,-32}}]
f[1.3]
fitData[{{5,15},{200,2},{10,15}}]
f[1.3]
fitData[{{-5,15},{-2,2},{3,15}}]
f[1.3]

#6, p. 121

data={{1,0.7},{2,0.75}};
fitData[data]
f[1.3]
data={{0,0.8},{1,0.7},{2,0.75}};
fitData[data]
f[1.3]
data={{0,0.8},{1,0.7},{2,0.75},{3,0.5}};
fitData[data]
f[1.3]

#7, p. 121

data={{2,0.8},{3.5,0.7}};
fitData[data]
f[3]
data={{2,0.8},{3.5,0.7},{4,0.75}};
fitData[data]
f[3]
data={{2,0.8},{3.5,0.7},{4,0.75},{5,0.5}};
fitData[data]
f[3]

#19, p. 122

data={{1,1},{2,3},{3,5},{4,7}};
fitData[data]

#22, p. 122

Manipulate[
Module[{y0=E^x0,y1=E^x1},
f[x_]:= y0 lagrange[x,x0,x1]+ y1 lagrange[x,x1,x0];
Show[
Plot[{E^x,f[x],E^x-f[x]},{x,0,1},PlotRange->{-1,E},PlotLabel->Together[f[x]]],
ListPlot[{{x0,y0},{x1,y1}},PlotStyle->{PointSize->Large}]
]
],
{{x0,0},0,x1},
{{x1,1},x0,1}
]

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

nevillesMethod[1.5,{{0,1},{1,E},{2,E^2}},False]
fakeExp[x_]:=nevillesMethod[x,{{0,1},{1,E},{2,E^2}},False]
Plot[{E^x,fakeExp[x]},{x,0,2}]

Farcsine:

Taylor polynomials fit arcsine with a polynomial, but one whose coefficients are determined by a single point and its higher derivatives (up to the nth derivative). So they do very well locally, creating tangent lines, parabolas,  cubics, etc. 

Lagrange polynomials take a different approach, fitting the function by matching at a number of function values (and ignoring higher derivatives). This certainly makes sense when one only has tabular data from an unknown (and perhaps unknowable) function; where any “derivatives” would simply be estimates, anyway.

One can easily imagine hybrid methods of these two ideas, provided the function is known. For example, cubic splines may use their four degrees of freedom to compute their coefficients by interpolating two points, but using two pieces of information at each point: function values and the derivatives at those points.

data=Table[{x,ArcSin[x]},{x,-1,1,.2}]; (* Let's fake ArcSin using 11 equally spaced points: *)
farcsine[x_]:=nevillesMethod[x,data,False]

Where does farcsine struggle? You may well guess that it’s at the endpoints....

{nevillesMethod[.99,data,True],ArcSin[.99]}
Plot[{ArcSin[x],farcsine[x]},{x,-1,1}]