Homework for Interpolation
Andy Long

Problem 1: a Hermite problem 
p. 179, #8
p. 190, #4, 7
p. 221, #9, 17

Hermite Spline Interpolation

Andy Long
Northern Kentucky University
longa@nku.edu
3/26/15

Let’s interpolate Sine at enough points to allow us to compute Sine for any real value (we can use symmetry and periodicity to get the rest). That means we’ll work on the interval [0,Pi/2]. 

You can use any other function, on other limits, by simply replacing the values below:

f[x_]:=Log[x+1]
xmin=-.5;
xmax=.5;

Now we’ll work with n points over our limits (evenly placed from xmin to xmax):

n=3;
xs=N[Table[(xmax-xmin)*i/n+xmin,{i,0,n}]];
x0=xmin;
x1=xs[[2]];
x2=xs[[3]];
xn=xmax;
p1=Plot[f[x],{x,xmin,xmax}];
p2=ListPlot[Table[{xs[[i]],f[xs[[i]]]},{i,0,n+1}],PlotMarkers->{Automatic,Medium}];
Show[p1,p2]

This process uses divided differences: we can do it without this, but it’s worth augmenting divided differences so that it can handle the case when x-values are repeated (and divided differences become derivatives).

dividedDifference[list_]:=
Module[
{n=Length[list]},
If[n==1,
b=list[[1]][[2]],
If[n==2 && list[[1]][[1]]==list[[2]][[1]],
b=list[[1]][[3]],
b=(dividedDifference[list[[2;;n]]]-dividedDifference[list[[1;;n-1]]])/
(list[[n]][[1]]-list[[1]][[1]])];
];
b
]
f'[x0]
f'[x1]
dividedDifference[{{x0,f[x0],f'[x0]},{x0,f[x0],f'[x0]}}]
dividedDifference[{{x1,f[x1],f'[x1]},{x1,f[x1],f'[x1]}}]
dividedDifference[{{x0,f[x0],f'[x0]},{x1,f[x1],f'[x1]}}]

HermiteCoefs finds the four Hermite cubic coefficients to fit two points and the derivative of f at those points:

HermiteCoefs[f_,x1_,x2_]:=
Module[
{d21=x2-x1,y1=f[x1],y2=f[x2],y1p=f'[x1],y2p=f'[x2]},
dd=dividedDifference[{{x1,y1},{x2,y2}}];
{
y1,
y1p,
(dd-y1p)/d21,
(y1p+y2p-2*dd)/d21^2
}
]
HermiteCoefs[f,x0,x1]

HermiteCoefsList creates a list of the Hermite cubic coefficients for a set of abscissa:

HermiteCoefsList[f_,xs_]:=
Module[{n=Length[xs],i=1,result={}},
While[i<n,
result=Append[result,HermiteCoefs[f,xs[[i]],xs[[i+1]]]];
i=i+1;
];
result
]
coefs=HermiteCoefsList[f,xs]

HermitePoly is the Hermite cubic polynomial fitting two points, which evaluates at any value of x:

HermitePoly[x_,x1_,x2_,coefs_]:=
Module[{d21=x2-x1},
dx1=x-x1;
dx2=x-x2;
coefs[[1]]+dx1*(coefs[[2]]+dx1*(coefs[[3]]+dx2*coefs[[4]]))
]
Plot[{f[x],HermitePoly[x,x0,x1,HermiteCoefs[f,x0,x1]]},{x,x0,x1}]
Plot[f[x]-HermitePoly[x,x0,x1,HermiteCoefs[f,x0,x1]],{x,x0,x1}]

HermiteSpline is the Hermite cubic spline, which evaluates at any value of x. 

Here are the answers to my questions for the homework:

HermiteSpline[x_,xs_,coefs_]:=
Module[{i=2,n=Length[xs]-1},
While[x>xs[[i]] && i≤n,
i=i+1
];
HermitePoly[x,xs[[i-1]],xs[[i]],coefs[[i-1]]]
]
Plot[{f[x],HermiteSpline[x,xs,coefs]},{x,xmin,xmax}]
Plot[f[x]-HermiteSpline[x,xs,coefs],{x,xmin,xmax}]
HermiteSpline[.3,xs,coefs]
HermiteSpline[.4,xs,coefs]
(* See http://ceadserv1.nku.edu/longa//classes/mat360/days/day19.html for error *)
fp[x_]=Abs[D[f[x],{x,4}]]
(* The error in the polynomial part of the error is maximized at the midpoint of each interval -- that makes each of the differences (x-x1) and (x-x2) equal to half the interval width, or 1/6 *)
N[fp[-1/2+0/3]]*(1/6)^4
N[fp[-1/2+1/3]]*(1/6)^4
N[fp[-1/2+2/3]]*(1/6)^4

As we can see, extrapolation with this spline is dangerous:

Plot[{f[x],HermiteSpline[x,xs,coefs]},{x,2*xmin-(xmax-xmin),2xmax-xmin}]

p. 221, #9: we fit the three points with an interpolating polynomial

x1=10;
x2=11;
x3=15;
y1=70;
y2=69;
y3=68;
(* Now just for fun here's the interpolating polynomial: *)
a=dividedDifference[{{x2,y2},{x3,y3},{x1,y1}}]
b=dividedDifference[{{x2,y2},{x3,y3}}]-(x2-x3)*a
c=y3
qH[x_]:=y3+(x-x3)*(b+(x-x3)*a)
p1=Plot[qH[x],{x,x1,x3}];
p2=ListPlot[{{x1,y1},{x2,y2},{x3,y3}},PlotMarkers->{Automatic,Medium}];
Show[p1,p2]
(* here's the estimated value: *)
N[qH[14]]

Illustrating p. 190, #4, let’s show another inteprolating polynomial:

p1=Plot[{qH[x],qH[x]+.1*(x-x1)*(x-x2)*(x-x3)},{x,x1,x3}];
p2=ListPlot[{{x1,y1},{x2,y2},{x3,y3}},PlotMarkers->{Automatic,Medium}];
Show[p1,p2]

p. 221, #17: we fit the three points with an interpolating polynomial

x1=2;
x2=2.2;
x3=2.4;
y1=.30103;
y2=.34242;
y3=.38021;
(* here's the second derivative approximation: *)
2* dividedDifference[{{2,.30103,},{2.2,.34242},{2.4,.38021}}]
(* Now just for fun here's the interpolating polynomial: *)
a=dividedDifference[{{x2,y2},{x3,y3},{x1,y1}}]
b=dividedDifference[{{x2,y2},{x3,y3}}]-(x2-x3)*a
c=y3
qH[x_]:=y3+(x-x3)*(b+(x-x3)*a)
qH''[2]
p1=Plot[qH[x],{x,x1,x3}];
p2=ListPlot[{{x1,y1},{x2,y2},{x3,y3}},PlotMarkers->{Automatic,Medium}];
Show[p1,p2]