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_]:=Sin[x]
fmin=0;
fmax=Pi/2;

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

n=5;
xs=N[Table[(fmax-fmin)*i/n+fmin,{i,0,n}]];
x0=fmin;
x1=xs[[2]];
x2=xs[[3]];
xn=fmax;
p1=Plot[f[x],{x,fmin,fmax}];
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:

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,fmin,fmax}]
Plot[f[x]-HermiteSpline[x,xs,coefs],{x,fmin,fmax}]

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

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