Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Sept 2016
Lecture topic: Lagrange Interpolation

Lagrange Interpolating Polynomials

Lagrange Polynomials

Clear[n,k,x,xp,i,Lag]

Lag[n_,k_,x_]=((∏_(i=1)^(k-1))((x-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i])))

Lag[4,2,x]

LP[n_,x_]=(∑_(k=1)^n)(f[xp[k]]*Lag[n,k,x])

Example of Lagrange Polynomials: Continuous Function

Notice that I simply use a delayed assignment so the Lagrange interpolating polynomial is not defined until I set up the the values I want the function evaluated at:

Clear[n,k,x,xp,i,Lag]
Lag[n_,k_,x_]:=((∏_(i=1)^(k-1))((x-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i])))
LP[n_,x_]:=Simplify[(∑_(k=1)^n)(f[xp[k]]*Lag[n,k,x])]

len=6;
f[x_] = Cos[x]Exp[-x/2];
xp[1]=0.0; 
xp[2]=2.0;
xp[3]=1.5;
xp[4]=1.3;
xp[5]=0.4;
xp[6]=1.300001;

Here are the Lagrange interpolating polynomials for a few different values of n:

Table[{n,LP[n,x]},{n,1,len}]//TableForm

Here are plots for a few different values of n:

points={};
Do[
points=Append[points,Graphics[Disk[{xp[n],f[xp[n]]},0.05]]];

plot1=Plot[{LP[n,x],f[x]},{x,-2,5},PlotRange->{{-1,3},{-1,1.5}},Axes->True,PlotStyle->{{Red,Thick},{Blue,Thick}},ImageSize->{500,Automatic}];

plot2[n]=Show[plot1,points,BaseStyle->20];
,{n,1,len,1}]

Manipulate[plot2[n],{n,1,len,1}]

Example of Lagrange Polynomials: Data Set

Here we are given a data set rather than a function:

Clear[n,k,x,xp,i,Lag]

Lag[n_,k_,x_]:=((∏_(i=1)^(k-1))((x-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i])))
LP[n_,x_]:=Simplify[(∑_(k=1)^n)(yp[k]*Lag[n,k,x])]

len=6;
xp[1]=8.3;yp[1]=17.56492;
xp[2]=8.6;yp[2]=18.50515;
xp[3]=8.7;yp[3]=18.12091;
xp[4]=8.1; yp[4]=16.94410;
xp[5]=9.3;yp[5]=13.56492;
xp[6]=9.0;yp[6]=15.55893;

Here are the Lagrange interpolating polynomials for a few different values of n:

Table[{n,LP[n,x]},{n,1,len}]//TableForm

Here is a plot that shows how the interpolating polynomial changes as we add more data points (since the interpolating polynomial is forced to pass through all data points, it can change dramatically)

points={};
Do[
points=Append[points,Graphics[{Blue,Disk[{xp[n],yp[n]},{0.03,0.03/(2/6)}]}]];

plot1=Plot[LP[n,x],{x,8,10},PlotRange->{{8,10},{13,19}},Axes->True,PlotStyle->{{Red,Thick}},ImageSize->{500,Automatic}];

plot2[n]=Show[plot1,points,AspectRatio->1];
,{n,1,len,1}]

Manipulate[plot2[n],{n,1,len,1}]

Example of Lagrange Polynomials: Oscillations (Wiggles)

Here is a data set that shows some significant oscillation inside the region we wish to interpolate.

Clear[n,k,x,xp,i,Lag]

Lag[n_,k_,x_]:=((∏_(i=1)^(k-1))((x-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i])))
LP[n_,x_]:=Simplify[(∑_(k=1)^n)(yp[k]*Lag[n,k,x])]

xp[1]=0.9;
yp[1]=1.3;
xp[2]=1.3;
yp[2]=1.5;
xp[3]=1.9;
yp[3]=1.85;
xp[4]=2.1; 
yp[4]=2.1;
xp[5]=2.6; 
yp[5]=2.6;
xp[6]=3.0; 
yp[6]=2.7;
xp[7]=3.9; 
yp[7]=2.4;
xp[8]=4.4; 
yp[8]=2.15;
xp[9]=4.7; 
yp[9]=2.05;
xp[10]=5.0; 
yp[10]=2.1;
xp[11]=6.0; 
yp[11]=2.25;
xp[12]=7.0; 
yp[12]=2.3;
xp[13]=8.0; 
yp[13]=2.25;
xp[14]=9.2; 
yp[14]=1.95;
xp[15]=10.5; 
yp[15]=1.4;
xp[16]=11.3; 
yp[16]=0.9;
xp[17]=11.6; 
yp[17]=0.7;

data=Table[{xp[n],yp[n]},{n,1,17}];

dataplot=ListPlot[data,PlotStyle->{Blue,PointSize[0.02]},Frame->True,ImageSize->{500,Automatic}]

lpplot=Plot[LP[17,x],{x,0,15},PlotRange->{{0,15},{-2,6}},FrameLabel->{"x","y"},PlotStyle->{Red,Thick},ImageSize->{500,Automatic}];

Show[lpplot,dataplot,ImageSize->{500,Automatic}]

We shall see how to minimize the wiggles using Chebyshev polynomials later.

Example of Lagrange Polynomials: Continuous Function Module on equally spaced points in [-1,1]

This module creates the n-1 Lagrange interpolating polynomial  (Subscript[P, n-1])(x) on equally spaced points in the interval [-1,1] for a given continuous function. Notice that the result is returned as a function of t, which we can convert to x when we call the module.

LP[fcn_,n_,t_]:=Module[{a,b,h,xp,i,k,res},
a=-1.0;
b=1.0;
h=(b-a)/(n-1);
Do[xp[i]=a+h*(i-1),{i,1,n}];
res=Sum[(fcn/.x->xp[k])*((∏_(i=1)^(k-1))((t-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((t-xp[i])/(xp[k]-xp[i]))),{k,1,n}];
Expand[res]
]

LP[Cos[2x+3]^2/(1+15x^2),4,x]

Here is some graphical output, notice the problem with osciallations near the ends of the interval.

f[x_]=Cos[2x+3]^2/(1+15x^2)
len=13;
Do[
plot[n]=Plot[{LP[f[x],n,t],f[t]},{t,-1,1},PlotRange->{{-1,1},{0,1}},PlotStyle->{{Thick,Red},{Thick,Blue}}]
,{n,2,len}]
Table[{plot[n],plot[n+1],plot[n+2]},{n,2,len,3}];
Show[GraphicsGrid[%],ImageSize->{650,Automatic}]

Alternate Ways to Create Lagrange Interpolating Polynomials

Neville's Method: Creating Lagrange Interpolating Functions Recursively

Let's use Neville's method to construct a series of Lagrange Approximating Polynomials. Neville’s method is typically used on n+1 data points starting at (Subscript[x, 0])to (Subscript[x, n]), whereas what we used before was for n data points started at (Subscript[x, 1])to (Subscript[x, n]).

The recursion relations are given by:

(Subscript[q, i,0])=(Subscript[y, i]), 0≤i≤n
(Subscript[q, i,j])=(((x-(Subscript[x, i-j]))(Subscript[q, i,j-1])-(x-(Subscript[x, i]))(Subscript[q, i-1,j-1]))/((Subscript[x, i])-(Subscript[x, i-j]))),1≤i≤n,1≤j≤i

Note that the Wikipedia entry for Neville’s method uses an equivalent but different notation.

Clear[q,i,j,npts,x,xp]

f[x_] = Cos[x]Exp[-x/2];
xp[0]=0.0; 
xp[1]=0.5;
xp[2]=1.5;
xp[3]=2.0;
xp[4]=2.5;

npts= 5;

Do[q[i,0]=f[xp[i]],{i,0,npts-1}]

For[i=1,i≤npts-1,i=i+1,

For[j=1,j≤i,j=j+1,
q[i,j]=((x-xp[i-j])q[i,j-1]-(x-xp[i])q[i-1,j-1])/(xp[i]-xp[i-j])
]
]

Table[Simplify[q[i,i]],{i,0,4}]//TableForm

These are the same functions we get if we use our previous way of computing the Lagrange interpolating polynomial.

Clear[n,k,x,xp,i,Lag]

Lag[n_,k_,x_]:=((∏_(i=1)^(k-1))((x-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i])))
LP[n_,x_]:=Simplify[(∑_(k=1)^n)(yp[k]*Lag[n,k,x])]

len=5;
xp[1]=0.0;yp[1]=f[xp[1]]; 
xp[2]=0.5;yp[2]=f[xp[2]]; 
xp[3]=1.5;yp[3]=f[xp[3]]; 
xp[4]=2.0;yp[4]=f[xp[4]]; 
xp[5]=2.5;yp[5]=f[xp[5]]; 

Table[LP[i,x],{i,1,5}]//TableForm

Divided Differences

A Sequence (x,y)

Here is a set of data points, and we want to construct the Divided Difference Table for them.

Clear[p,datap]

datap={{1.0,0.7651977},{1.3,0.6200860},{1.6,0.4554022},{1.9,0.2818186},{2.2,0.1103623}}

Clear[p]

This the Mathematica syntax to get it done:

Clear[f,divdif]

Do[
divdif[{Part[datap,n,1]}]=Part[datap,n,2]
,{n,1,Length[datap]}];

divdif[{p_,q___,r_}]:=((divdif[{q,r}]-divdif[{p,q}])/((r-p)));

NOTE: Triple underscore means we can have zero or more elements in that position (q here). This is imperative for the divided difference evaluation!

First divided differences:

TableForm[
Table[divdif[{Part[datap,n-1,1],Part[datap,n,1]}],{n,2,Length[datap]}]
]

Second divided differences:

TableForm[
Table[divdif[{Part[datap,n-2,1],Part[datap,n-1,1],Part[datap,n,1]}],{n,3,Length[datap]}]
]

Third divided differences:

TableForm[
Table[divdif[{Part[datap,n-3,1],Part[datap,n-2,1],Part[datap,n-1,1],Part[datap,n,1]}],{n,4,Length[datap]}]
]

Fourth divided differences:

TableForm[
Table[divdif[{Part[datap,n-4,1],Part[datap,n-3,1],Part[datap,n-2,1],Part[datap,n-1,1],Part[datap,n,1]}],{n,5,Length[datap]}]
]

A Sequence with polynomial growth

Clear[p,datap]

datap= Table[{n,8n^4+2n^2+n+7},{n,0,10}]

Clear[p]

Clear[f,divdif]

Do[divdif[{Part[datap,n,1]}]=Part[datap,n,2],{n,1,Length[datap]}];

divdif[{p_,q___,r_}]:=((divdif[{q,r}]-divdif[{p,q}])/((r-p)));

First divided differences:

TableForm[Table[divdif[{Part[datap,n-1,1],Part[datap,n,1]}],{n,2,Length[datap]}]]

Second divided differences:

TableForm[Table[divdif[{Part[datap,n-2,1],Part[datap,n-1,1],Part[datap,n,1]}],{n,3,Length[datap]}]]

Third divided differences:

TableForm[Table[divdif[{Part[datap,n-3,1],Part[datap,n-2,1],Part[datap,n-1,1],Part[datap,n,1]}],{n,4,Length[datap]}]]

Fourth divided differences:

TableForm[Table[divdif[{Part[datap,n-4,1],Part[datap,n-3,1],Part[datap,n-2,1],Part[datap,n-1,1],Part[datap,n,1]}],{n,5,Length[datap]}]]

Fifth divided differences:

TableForm[Table[divdif[{Part[datap,n-5,1],Part[datap,n-4,1],Part[datap,n-3,1],Part[datap,n-2,1],Part[datap,n-1,1],Part[datap,n,1]}],{n,6,Length[datap]}]]

A sequence that grows asymptotically as α (x^n) will have the nth divided differences as α, and the (n+1)th divided differences will be zero. This can be useful for estimating the polynomial growth of sequences where the underlying function is unknown.

Newton's Interpolatory Divided Difference Formula

Data Set

Here is a data set we would like to interpolate on.

data={{0.9,1.3},{1.3,1.5},{1.9,1.85},{2.1,2.1},{2.6,2.6},{3.0,2.7},{3.9,2.4},{4.4,2.15}};

dataplot=ListPlot[data,PlotStyle->{Blue,PointSize[0.02]},Frame->True]

Algorithm

This sets up what we need to define Newton's Interpolatory Divided Difference Formula. There are some changes to the indexing due to labelling issues (0 to 16 or 1 to 17)

Clear[f,divdif,len,list,prod,LP]

len=Length[data]-1;

Do[divdif[{Part[data,n,1]}]=Part[data,n,2],{n,1,len+1}];

divdif[{p_,q___,r_}]:=divdif[{p,q,r}]=((divdif[{q,r}]-divdif[{p,q}])/(r-p));

list[0]={Part[data,1,1]};
Do[list[k]=Append[list[k-1],Part[data,k+1,1]],{k,1,len}];
prod[0]=(x-Part[data,1,1]);
Do[prod[k]=prod[k-1]*(x-Part[data,k+1,1]),{k,1,len}];

This is Equation 3.2 in the text:

LP[x_]= divdif[list[0]]+(∑_(k=1)^len)divdif[list[k]]*prod[k-1];

LP[x]

lpplot=Plot[LP[x],{x,0,15},PlotRange->{{0,6},{1,3}},FrameLabel->{"x","y"},Frame->True,Axes->False,PlotStyle->{Red,Thick}];

Show[lpplot,dataplot,ImageSize->{500,Automatic}]

Here is the old way of determining the Lagrange interpolating polynomial:

Clear[n,k,x,xp,i,Lag]
Lag[n_,k_,x_]:=((∏_(i=1)^(k-1))((x-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i])))
oldLP[n_,x_]:=Simplify[(∑_(k=1)^n)(yp[k]*Lag[n,k,x])]

Do[
xp[i]=data[[i,1]];
yp[i]=data[[i,2]];
,{i,1,Length[data]}]

oldLP[Length[data],x]

oldlpplot=Plot[oldLP[Length[data],x],{x,0,15},PlotRange->{{0,6},{1,3}},FrameLabel->{"x","y"},Frame->True,Axes->False,PlotStyle->{Red,Thick}];

Show[oldlpplot,dataplot,ImageSize->{500,Automatic}]

Simplify[oldLP[Length[data],x]-LP[x]]

Mathematica Built-in Command

In Mathematica, this procedure is automated using the command InterpolatingPolynomial

mathfit[x_]=InterpolatingPolynomial[data,x];
Plot[mathfit[x],{x,0,15},PlotRange->{{0,6},{1,3}},FrameLabel->{"x","y"},PlotStyle->{Red,Thick},ImageSize->{500,Automatic}]

References:

Heikki Ruskeepaa, Mathematica Navigator, Graphics and Methods of Applied Mathematics (1999) San Diego,Academic Press.
Richard Burden & J. Douglas Faires, Numerical Analysis, 7th ed. (2001) Pacific Grove, Brooks/Cole.