Math 4401 Numerical Methods

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

Hermite Interpolation: Including Derivative Information

Algorithm: Data Set, Derivative supplied (Work from Original Definition)

Here are four data points that are the form ((Subscript[x, i]),(Subscript[y, i]),dy/(Subscript[dx, i]))

Clear[data,points]
data={{1.3,0.6200860,-0.5220232},{1.65,0.4554022,-0.5698959},{1.9,0.2818186,-0.5811571},{2.1,0.2818186,-0.5811571}};

points=Table[{Part[data,i,1],Part[data,i,2]},{i,1,Length[data]}]

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

Here is one way of creating the Hermite interpolating polynomial. The bit in red is a hack to get the derivative calculated and evaluated correctly.

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

deriv[n_,j_,x_]:=D[Lag[n,j,x],x]

HP[n_,x_]:=Simplify[(∑_(j=1)^n)yp[j]*(1-2(x-xp[j])deriv[n,j,xx]/.xx->xp[j])((Lag[n,j,x])^2)+(∑_(j=1)^n)dp[j]*(x-xp[j])*((Lag[n,j,x])^2)]

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

Here is the polynomial (n data points means you should use HP[n-1,x]):

HP[4,x]

Note that it is a 7th degree polynomial! There was twice as much information to encode in the polynomial (y value and y'), so we needed twice as much flexibility, which leads to a polynomial which is twice as long as the Lagrange interpolating polynomial on three points. Here is the Lagrange interpolating polynomial to compare:

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[Sum[yp[k]*Lag[n,k,x],{k,1,n}]]

LP[4,x]

Here is a plot of the Hermite (red) and Lagrange (green) interpolating polynomial:

data2={};
Table[data2=Append[data2,{xp[i],yp[i]}],{i,1,4}];

plotdata=ListPlot[data2,PlotStyle->{PointSize[0.03],Blue}];

plot1=Plot[{HP[4,x],LP[4,x]},{x,0,6},PlotRange->{All,{-1,1}},FrameLabel->{"x","y"},
Frame->True,
Axes->False,
PlotStyle->{{Thick,Red},{Thick,Darker[Green]}}];

Show[plot1,plotdata,ImageSize->{500,Automatic},BaseStyle->20]

Show[plot1,plotdata,
ImageSize->{500,Automatic},
PlotRange->{{1,2.3},{0.15,1}},BaseStyle->20]

The divided difference notation can help us understand this a bit better.

Algorithm: Data Set, Derivative supplied (Work from Divided Difference Representation)

Clear[data,points]
data={{1.3,0.6200860,-0.5220232},{1.65,0.4554022,-0.5698959},{1.9,0.2818186,-0.5811571},{2.1,0.2818186,-0.5811571}};

points=Table[{Part[data,i,1],Part[data,i,2]},{i,1,Length[data]}]

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

First, it is necessary to create the doubled data set:

Clear[datadb,z,fz,fpz]

len=2*Length[data]-1;
z={};
fz={};
fpz={};
Do[
z=Append[z,Part[data,k+1,1]];
fz=Append[fz,Part[data,k+1,2]];
fpz=Append[fpz,Part[data,k+1,3]];
z=Append[z,Part[data,k+1,1]];
fz=Append[fz,Part[data,k+1,2]];
fpz=Append[fpz,Part[data,k+1,3]],
{k,0,Length[data]-1}]
z
fz
fpz

Here is the divided difference scheme:

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

len=2*Length[data]-1;

(*This next line initializes the first divided differences*)
Do[divdif[{Part[z,n]}]=Part[fz,n],{n,1,len+1}];

(*This second line initializes half of the second divided differences, and computes the other half*)
Do[
divdif[{Part[z,n],Part[z,n+1]}]=
If[OddQ[n],Part[fpz,n],((divdif[{Part[z,n+1]}]-divdif[{Part[z,n]}])/(Part[z,n+1]-Part[z,n]))],
{n,1,len}];

(*The next line implements the divided differences recursively*)
divdif[{p_,q___,r_}]:=divdif[{p,q,r}]=((divdif[{q,r}]-divdif[{p,q}])/(r-p));

(*The next line creates the list and prod we need to create the Hermite approximating polynomial*)
list[0]={Part[z,1]};
Do[list[k]=Append[list[k-1],Part[z,k+1]],{k,1,len}];
prod[0]=(x-Part[z,1]);
Do[prod[k]=prod[k-1]*(x-Part[z,k+1]),{k,1,len}];

This is Equation at bottom of Page 133 in the text (old text!):

HP[x_]= divdif[list[0]]+(∑_(k=1)^len)divdif[list[k]]*prod[k-1];
Expand[%,x]

hpplot=Plot[HP[x],{x,0,15},
PlotRange->{{1,2.3},{0.15,1}},
FrameLabel->{"x","y"},
Frame->True,
Axes->False,
PlotStyle->{Thick,Red},
ImageSize->{500,Automatic}];

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

Which is exactly what we found before.

Algorithm: Data Set, no Derivative Provided: Be Careful!

If no derivative information is provided, you would be best to stick with Lagrange interpolation, as this simple example clearly shows.

Clear[data,points]
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}};

points=Table[{Part[data,i,1],Part[data,i,2]},{i,1,Length[data]}]

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

First, it is necessary to create the doubled data set. Here was have make up some derivative information from the given data set, which I've highlighted in red. This might seem like a good idea, but it doesn’t work very well, as we shall see.

Clear[datadb,z,fz,fpz]

len=2*Length[data]-1;
z={};
fz={};
fpz={};
Do[
z=Append[z,Part[data,k+1,1]];
fz=Append[fz,Part[data,k+1,2]];
fpz=If[k==0,Append[fpz,0.5],Append[fpz,((Part[data,k+1,2]-Part[data,k,2])/(Part[data,k+1,1]-Part[data,k,1]))]];
z=Append[z,Part[data,k+1,1]];
fz=Append[fz,Part[data,k+1,2]];
fpz=If[k==0,Append[fpz,0.5],Append[fpz,((Part[data,k+1,2]-Part[data,k,2])/(Part[data,k+1,1]-Part[data,k,1]))]],
{k,0,Length[data]-1}]
z
fz
fpz

This implements the divided difference technique to get the Hermite interpolating polynomial:

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

len=2*Length[data]-1;

Do[divdif[{Part[z,n]}]=Part[fz,n],{n,1,len+1}];

Do[
divdif[{Part[z,n],Part[z,n+1]}]=
If[OddQ[n],Part[fpz,n],((divdif[{Part[z,n+1]}]-divdif[{Part[z,n]}])/(Part[z,n+1]-Part[z,n]))],
{n,1,len}];

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

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

Here is our approximation:

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

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

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

Obviously, we have to be more careful about how we incorporate derivative information than we were here! If you have no derivative information, stick with the Lagrange Interpolating Polynomial.

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

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

LP[Length[data]-1,x];
Expand[%]

lpplot=Plot[LP[Length[data]-1,x],{x,0,15},
PlotRange->{{0,6},{1,3}},
FrameLabel->{"x","y"},
PlotStyle->{Green,Thick},
ImageSize->{500,Automatic}];

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

The Lagrangian worked much better in this case.

Algorithm: General Divided Difference Representation Formula for n=2 points

Here we do NOT define the divided difference table entries, just so we can see the structure at the end.

Clear[z,f,divdif,len,list,prod,LP]
n=2;
len=2*n-1;

varlist={};
Do[
varlist=Append[varlist,z[k]]; 
,{k,0,len}]

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

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

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.