MS-C1650 - Numerical Analysis
Week 18: Lecture 5: Piecewise Interpolation

Harri Hakula, 2021

Hermite Interpolation (Textbook version)

Let us implement the method in its pure format.

HermiteInterval[{{x0_?NumericQ,f0_?NumericQ,df0_?NumericQ},{x1_?NumericQ,f1_?NumericQ,df1_?NumericQ}}]:=Block[{h,α,C,x},
h=x1-x0;
C=f0;
α=(3/h^2)(df0+df1)+(6/h^3)(f0-f1);
{
Simplify[-(df0/h)Integrate[t-x1,{t,x0,x}]+(df1/h)Integrate[t-x0,{t,x0,x}]+α Integrate[(t-x0)(t-x1),{t,x0,x}]+C],
x0≤x≤x1
}
]

HermiteInterpolation[data:{{_?NumericQ,_?NumericQ,_?NumericQ}..}]:=Block[{int},
Piecewise[HermiteInterval/@Partition[data,2,1]]
]

Examples

Case 1

fun=Function[ξ,3(ξ^4)-2(ξ^2)+Sin[ξ]+1];
Dfun=Function[ξ,Evaluate[D[fun[ξ],ξ]]];
DDfun=Function[ξ,Evaluate[D[fun[ξ],{ξ,2}]]];

xp=Subdivide[-1,1,3]

data=Transpose[{xp,fun[xp],Dfun[xp]}];

hpi=HermiteInterpolation[data]

Plot[{fun[x],hpi},{x,-1,1},Epilog->{Red,PointSize[0.02],Point[Transpose[{xp,fun[xp]}]]}]

Plot[{Dfun[x],Evaluate[D[hpi,x]]},{x,-1,1},Epilog->{Red,PointSize[0.02],Point[Transpose[{xp,Dfun[xp]}]]}]

Plot[{Dfun[x],Evaluate[D[hpi,x]]},{x,-1,1},Epilog->{Red,PointSize[0.02],Point[Transpose[{xp,Dfun[xp]}]]},PlotRange->All]

Plot[{DDfun[x],Evaluate[D[hpi,{x,2}]]},{x,-1,1},Epilog->{Red,PointSize[0.02],Point[Transpose[{xp,DDfun[xp]}]]},PlotRange->All]

Examples

Case 2

fun=Function[x,Cos[4x]Exp[2x]];
Dfun=Function[x,Evaluate[D[fun[x],x]]];
DDfun=Function[x,Evaluate[D[fun[x],{x,2}]]];

xp=Subdivide[-1,1,2]

data=Transpose[{xp,fun[xp],Dfun[xp]}];

hpi=HermiteInterpolation[data]

Plot[{fun[x],hpi},{x,-1,1}]

xp=Subdivide[-1,1,3]

data=Transpose[{xp,fun[xp],Dfun[xp]}];

hpi=HermiteInterpolation[data]

Plot[{fun[x],hpi},{x,-1,1},Epilog->{Red,PointSize[0.02],Point[Transpose[{xp,fun[xp]}]]}]

Plot[{Dfun[x],Evaluate[D[hpi,x]]},{x,-1,1},Epilog->{Red,PointSize[0.02],Point[Transpose[{xp,Dfun[xp]}]]}]

Plot[{Dfun[x],Evaluate[D[hpi,x]]},{x,-1,1},Epilog->{Red,PointSize[0.02],Point[Transpose[{xp,Dfun[xp]}]]},PlotRange->All]

Plot[{DDfun[x],Evaluate[D[hpi,{x,2}]]},{x,-1,1},Epilog->{Red,PointSize[0.02],Point[Transpose[{xp,DDfun[xp]}]]},PlotRange->All]

Hermite Interpolation (Non-standard)

Let us interpolate the function f(x) = 3(x^4)-2(x^2)+Sin[x]+1 over the interval x∈[-1,1].

Plot[3(x^4)-2(x^2)+Sin[x]+1,{x,-1,1}]

Choose points and find the values and the derivatives.

fun=Function[x,3(x^4)-2(x^2)+Sin[x]+1];
Dfun=Function[x,Evaluate[D[fun[x],x]]];
DDfun=Function[x,Evaluate[D[fun[x],{x,2}]]];

x=Subdivide[-1,1,6]

y=fun/@x

dy=Dfun/@x

data=MapThread[List[List[#1],#2,#3]&,{x,y,dy}]

if=Interpolation[data,Method->"Hermite",InterpolationOrder->1];

Plot[{fun[x],if[x]},{x,-1,1}]

if/@x
N[fun/@x]
Function[ξ,Evaluate[D[if[ξ],ξ]]]/@x
N[Dfun/@x]

N[Function[ξ,Evaluate[D[if[ξ],{ξ,2}]]]/@x]
N[DDfun/@x]

Splines (Natural spline, textbook version)

Let us implement the natural spline in its pure format.

SplineInterval[{
{x0_?NumericQ,f0_?NumericQ,z0_?NumericQ},{x1_?NumericQ,f1_?NumericQ,z1_?NumericQ}
}]:=Block[{x,h,c,d},
h=x1-x0;
c=(1/h)(f1-f0+(h^2/6)(z0-z1));
d=f0-(h^2/6)z0;
{
Expand[(1/6)(1/h)z0(x1-x)^3+(1/6)(1/h)z1(x-x0)^3+c(x-x0)+d],
x0≤x≤x1
}
]

SplineInterpolation[data:{{_?NumericQ,_?NumericQ}..}]:=Block[{z,h,f,A,b,dim},
dim=Length[data];
f=data[[All,2]];
z=ConstantArray[0,dim];
h=data[[2,1]]-data[[1,1]];
A=If[dim-2==1,
SparseArray[{Band[{1,1}]->2 h/3},{dim-2,dim-2}],
SparseArray[{Band[{1,1}]->2 h/3,Band[{1,2}]->h/6,Band[{2,1}]->h/6},{dim-2,dim-2}]
];
b=(1/h)(f[[1;;-3]]-2f[[2;;-2]]+f[[3;;-1]]);
z[[2;;-2]]=LinearSolve[A,N@b];
Piecewise[SplineInterval/@Partition[MapThread[Append,{data,z}],2,1]]
]

Slide

Case 1

fun=Function[x,3(x^4)-2(x^2)+Sin[x]+1];
Dfun=Function[x,Evaluate[D[fun[x],x]]];
DDfun=Function[x,Evaluate[D[fun[x],{x,2}]]];

xp=Subdivide[-1,1,3]

data=Transpose[{xp,fun[xp]}];

spi=SplineInterpolation[{{-1,1},{0,2},{1,-1}}]

spi=SplineInterpolation[data]

spi/.{x->-1/3}

Plot[{fun[x],spi},{x,-1,1}]

xp=Subdivide[-1,1,6]

data=Transpose[{xp,fun[xp]}];

spi=SplineInterpolation[data]

Plot[{fun[x],spi},{x,-1,1},Epilog->{Red,PointSize[0.02],Point[Transpose[{xp,fun[xp]}]]}]

Plot[{Dfun[x],Evaluate[D[spi,x]]},{x,-1,1},Epilog->{Red,PointSize[0.02],Point[Transpose[{xp,Dfun[xp]}]]}]

Plot[{DDfun[x],Evaluate[D[spi,{x,2}]]},{x,-1,1},Epilog->{Red,PointSize[0.02],Point[Transpose[{xp,DDfun[xp]}]]},PlotRange->All]

Slide

Case 2

fun=Function[x,Cos[4x]Exp[2x]];
Dfun=Function[x,Evaluate[D[fun[x],x]]];
DDfun=Function[x,Evaluate[D[fun[x],{x,2}]]];

xp=Subdivide[-1,1,3]

data=Transpose[{xp,fun[xp]}];

spi=SplineInterpolation[data]

Plot[{fun[x],spi},{x,-1,1}]

xp=Subdivide[-1,1,3]

data=Transpose[{xp,fun[xp]}];

spi=SplineInterpolation[data]

Plot[{fun[x],spi},{x,-1,1},Epilog->{Red,PointSize[0.02],Point[Transpose[{xp,fun[xp]}]]}]

Plot[{Dfun[x],Evaluate[D[spi,x]]},{x,-1,1},Epilog->{Red,PointSize[0.02],Point[Transpose[{xp,Dfun[xp]}]]}]

Plot[{DDfun[x],Evaluate[D[spi,{x,2}]]},{x,-1,1},Epilog->{Red,PointSize[0.02],Point[Transpose[{xp,DDfun[xp]}]]},PlotRange->All]

Spline  (Non-standard)

is=Interpolation[N@Transpose[{x,y}],Method->"Spline",InterpolationOrder->3];

Plot[{fun[x],if[x],is[x]},{x,-1,1}]

is/@x
N[fun/@x]
Function[ξ,Evaluate[D[is[ξ],ξ]]]/@x
N[Dfun/@x]

Function[ξ,Evaluate[D[is[ξ],{ξ,2}]]]/@x
N[DDfun/@x]