Mathematica Module for Hermite and Other Spline Interpolation
based on one from http://math.fullerton.edu/mathews (thanks!)
Andy Long
Numerical Analysis
Fall, 2009

Numerical Methods using Mathematica  
Complementary software supplements for   
Numerical Methods textbooks  
John H. Mathews and Kurtis Fink, 2004

Background for the Hermite Polynomial

	The cubic Hermite polynomial  p(x)  has the interpolative properties  p((Subscript[x, 0]))=f((Subscript[x, 0])),   p((Subscript[x, 1]))=f((Subscript[x, 1])),   p'((Subscript[x, 0]))=f'((Subscript[x, 0])),  and  p'((Subscript[x, 1]))=f'((Subscript[x, 1])); that is, p(x) fits both the function values and their derivatives at the endpoints of the interval  [(Subscript[x, 0]),(Subscript[x, 1])].  Hermite polynomials were studied by the French Mathematician Charles Hermite (1822-1901), and are referred to as a "clamped cubic," where "clamped" refers to the slope at the endpoints being fixed.  This situation is illustrated in the figure below.

Theorem (Cubic Hermite Polynomial).  If  f[x]  is continuous on the interval  [(Subscript[x, 0]),(Subscript[x, 1])], there exists a unique cubic polynomial  p[x] = a (x^3) + b (x^2) + c x + d  such that 
		   p[(Subscript[x, 0])] = f[(Subscript[x, 0])], 
		   p[(Subscript[x, 1])] = f[(Subscript[x, 1])], 
		p'[(Subscript[x, 0])] = f'[(Subscript[x, 0])],  
		p'[(Subscript[x, 2])] = f'[(Subscript[x, 2])].

Remark.  The cubic Hermite polynomial is a generalization of both the Taylor polynomial and Lagrange polynomial, and it is referred to as an "osculating polynomial."  Hermite polynomials can be generalized to higher degrees by requiring that they use of more nodes  {(Subscript[x, 0]),(Subscript[x, 1]),...,(Subscript[x, n])} and the extension to agreement at higher derivatives  Derivative[k][p][(Subscript[x, i])] =Derivative[k][f][(Subscript[x, i])]  for  i = 1,2,...,n  and    k = 1,2,...,(Subscript[m, i]).

Example 1.  Find the cubic Hermite polynomial or "clamped cubic" that satisfies  
		   p(1)=2
p'(1)=1
   p(3)=1
p'(3)=2

Solution 1. Follow the strategy we used in class, using divided differences:

x0=1
x1=3
dx=x1-x0
f0=2
f1=1
f00=1
f11=2
f01=(f1-f0)/dx
f001=(f01-f00)/dx
f011=(f11-f01)/dx
f0011=(f011-f001)/dx

p[x_]=f0+f00*(x-x0)+f001*(x-x0)^2+f0011*(x-x0)^2*(x-x1)
Plot[{p[x]},{x,x0, x1},PlotStyle->{Red}]

Solution 2. In contrast to the way we found the coefficients in class, we can proceed directly from the typical form of a polynomial as a sum of monomials, then use some linear algebra to find the cubic polynomial (the same one that we found above, only it will be expressed a little differently):

Enter the formula for a general cubic equation.

p[x_] = a (x^3) + b (x^2) + c x + d ; 
Print["p[x] = ",p[x] ];

Symbolic differentiation (integration too) is permitted with Mathematica.

Print["p'[x] = ",p'[x] ];

Set up four equations using the prescribed endpoint conditions. Then find the solution set to this linear system and store it in the variable solset.

eq1 =    p[1] == 2;  
eq2 = p'[1] == 1;  
eq3 =    p[3] == 1;  
eq4 = p'[3] == 2;  
eqns = {eq1,eq2,eq3,eq4};  
vars = {a,b,c,d};
solset = Solve[eqns,vars];  
Print[TableForm[eqns] ];  
Print[solset];

Use the solution given above for the coefficients and form the cubic function.  We must dig out one set of braces using  (Subscript[solset, [[1]]])  before we can use the ReplaceAll command.

(f[x_]=p[x]/.solset[[1]];) (xy={{1,f[1]},{3,f[3]}};) (dots=ListPlot[xy,PlotStyle->{Red,PointSize[0.02`]},DisplayFunction->Identity];) (gr=Plot[f[x],{x,1,3},PlotRange->{{0,3},{0,2.2`}},PlotStyle->{Magenta},AspectRatio->((2.2`)/3),DisplayFunction->Identity];) Show[gr,dots,DisplayFunction->$DisplayFunction] (Print["y = f[x] = ",f[x]];) (Print[" f[1] = ",f[1],",      f[3] = ",f[3]];) (Print["f'[1] = ",f'[1],",     f'[3] = ",f'[3]];)

The Clamped Cubic Spline

	A clamped cubic spline is obtained by forming a piecewise cubic function which passes through the given set of knots ((Subscript[x, 0]),(Subscript[y, 0])),((Subscript[x, 1]),(Subscript[y, 1])),...,((Subscript[x, n]),(Subscript[y, n]))  with the condition the function values, their derivatives and second derivatives of adjacent cubics agree at the interior nodes.  The endpoint conditions are S'((Subscript[x, 0]))=(Subscript[d, 0])  and  S'((Subscript[x, n]))=(Subscript[d, n]), where (Subscript[d, 0]) and (Subscript[d, n])  are given. 
	Once again it pays to use our linear algebra to find the coefficients, as we do in the following simple example, in which we find the coefficients for a pair of such cubics.

Example 2.  Find the "clamped cubic spline" that satisfies  
		      (Subscript[S, 1])(1) = 2 
   (Subscript[S, 1])'(1) = 1
      (Subscript[S, 1])(2) = 1 
      (Subscript[S, 1])(2) = (Subscript[S, 2])(2)
   (Subscript[S, 1])'(2) = (Subscript[S, 2])'(2)
(Subscript[S, 1])''(2) = (Subscript[S, 2])''(2)
      (Subscript[S, 2])(3) = 1
   (Subscript[S, 2])'(3) = 2

Solution.

Set up the formulas for the two cubic polynomials and form the equations to solve.  Subscripted variables are used if you prefer you can use ordinary variables  (Subscript[a, 1]),(Subscript[b, 1]),(Subscript[c, 1]),(Subscript[d, 1]),(Subscript[a, 2]),(Subscript[b, 2]),(Subscript[c, 2]),(Subscript[d, 2]).

S1[x_] = (Subscript[a, 1]) (x^3) + (Subscript[b, 1]) (x^2) + (Subscript[c, 1]) x + (Subscript[d, 1]);  
S2[x_] = (Subscript[a, 2]) (x^3) + (Subscript[b, 2]) (x^2) + (Subscript[c, 2]) x + (Subscript[d, 2]);  
Print["S1[x_] = ",S1[x] ];  
Print["S2[x_] = ",S2[x] ];

Set up eight equations using the prescribed endpoint conditions. Then find the solution set to this linear system and store it in the variable solset.

eq1 =      S1[1] == 2;   
eq2 =   S1'[1] == 1;   
eq3 =      S1[2] ==  1;   
eq4 =      S1[2] == S2[2];   
eq5 =   S1'[2] == S2'[2];   
eq6 =S1''[2] == S2''[2];   
eq7 =      S2[3] == 1;   
eq8 =   S2'[3] == 2;   
eqns = { eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8 };   
vars = { (Subscript[a, 1]),(Subscript[b, 1]),(Subscript[c, 1]),(Subscript[d, 1]),(Subscript[a, 2]),(Subscript[b, 2]),(Subscript[c, 2]),(Subscript[d, 2]) };   
solset = Solve[eqns,vars];
Print[TableForm[eqns]];   
Print[solset];

Use the solution given above for the coefficients and form the cubic functions.  Remember that we must dig out one set of braces using  (Subscript[solset, [[1]]])  before we can use the ReplaceAll command.

g1[x_] = ReplaceAll[ S1[x],(Subscript[solset, [[1]]]) ];  
g2[x_] = ReplaceAll[ S2[x],(Subscript[solset, [[1]]]) ];  
Print["(Subscript[g, 1])[x] = ",g1[x] ];  
Print["(Subscript[g, 2])[x] = ",g2[x] ];

Now graph the portion of each cubic in the interval over which it is to be used.  Then combine the two piecewise cubic graphs to form the spline.

(xy={{1,g1[1]},{2,g1[2]},{3,g2[3]}};) (dots=ListPlot[xy,PlotStyle->{Red,PointSize[0.02`]},DisplayFunction->Identity];) (graph1=Plot[g1[x],{x,1,2},PlotRange->{{0,3},{0,2.2`}},PlotStyle->Green,DisplayFunction->Identity];) (graph2=Plot[g2[x],{x,2,3},PlotRange->{{0,3},{0,2.2`}},PlotStyle->Blue,DisplayFunction->Identity];) (graph12=Show[graph1,graph2,dots,AspectRatio->((2.2`)/3),DisplayFunction->$DisplayFunction]) (Print["\!\(\*SubscriptBox[\(g\), \(1\)]\)[x] = ",g1[x]];) (Print["\!\(\*SubscriptBox[\(g\), \(2\)]\)[x] = ",g2[x]];)

Remark.  Note that the individual cubics were graphed over different intervals and then the graphs were combined. If they were plotted together over a common interval it would look different, as you can see below. Look at how nicely the curves match at the interior boundary, however:

(graph0=Plot[{g1[x],g2[x]},{x,0,4},PlotRange->{{0,4},{0,2.2`}},PlotStyle->{Green,Blue},DisplayFunction->Identity];) (graphs=Show[graph0,dots,AspectRatio->((2.2`)/4),DisplayFunction->$DisplayFunction]) (Print["\!\(\*SubscriptBox[\(g\), \(1\)]\)[x] = ",g1[x]];) (Print["\!\(\*SubscriptBox[\(g\), \(2\)]\)[x] = ",g2[x]];)

Remark.  It would be nice to have one piecewise cubic function S[z] that is used. The following formulas for  S[x]  uses the condition syntax  /;  for making a piecewise function.  Under the help menu we can find the following information about Condition.
patt  /;  test is a pattern which matches only if the evaluation of test yields True

Clear[S];  
S[x_/;1≤x≤2] = g1[x];   
S[x_/;2≤x≤3] = g2[x];

Caution.  We cannot print the formula for a piecewise function with the Print command.  It is only possible to interrogate the system and determine what rule it is storing for  S.  Hence, it was a good idea to always use two formulas g1 and g2 to define  S, since they can be printed.

?S

When a plot is made, it only uses real numbers specified in the conditions.

(graph=Plot[S[x],{x,1,3},PlotRange->{{0,4},{0,2.2`}},PlotStyle->Magenta,AspectRatio->((2.2`)/4),DisplayFunction->Identity];) Show[graph,dots,DisplayFunction->$DisplayFunction] (Print["S[x_/;1≤x≤2] = ",g1[x]];) (Print["S[x_/;2≤x≤3] = ",g2[x]];)

The clamped cubic spline forces the slope at the endpoints to be prescribed values.

The Natural Cubic Spline

	A natural cubic spline is obtained by forming a piecewise cubic function which passes through the given set of knots ((Subscript[x, 0]),(Subscript[y, 0])),((Subscript[x, 1]),(Subscript[y, 1])),...,((Subscript[x, n]),(Subscript[y, n])) with the condition the function values, their derivatives and second derivatives of adjacent cubics agree at the interior nodes.  The endpoint conditions are  S''((Subscript[x, 0]))=0  and  S''((Subscript[x, n]))=0. The natural cubic spline is said to be "a relaxed curve."

Example 3.  Find the "natural cubic spline" that satisfies  
		      (Subscript[S, 1])(1) = 2
 (Subscript[S, 1])''(1) = 0
      (Subscript[S, 1])(2) = 1
      (Subscript[S, 1])(2) = (Subscript[S, 2])(2)
   (Subscript[S, 1])'(2) = (Subscript[S, 2])'(2)
(Subscript[S, 1])''(2) = (Subscript[S, 2])''(2)
     (Subscript[S, 2])(3) = 1
(Subscript[S, 2])''(3) = 0

Solution.

Set up the formulas for the two cubic polynomials and form the equations to solve. Subscripted variables are used if you prefer you can use ordinary variables  (Subscript[a, 1]),(Subscript[b, 1]),(Subscript[c, 1]),(Subscript[d, 1]),(Subscript[a, 2]),(Subscript[b, 2]),(Subscript[c, 2]),(Subscript[d, 2]).

S1[x_] = (Subscript[a, 1]) (x^3) + (Subscript[b, 1]) (x^2) + (Subscript[c, 1]) x + (Subscript[d, 1]);  
S2[x_] = (Subscript[a, 2]) (x^3) + (Subscript[b, 2]) (x^2) + (Subscript[c, 2]) x + (Subscript[d, 2]);  
Print["S1[x_] = ",S1[x] ];  Print["S2[x_] = ",S2[x] ];

Set up eight equations using the prescribed endpoint conditions. Then find the solution set to this linear system and store it in the variable solset.

eq1 =     S1[1] == 2;   
eq2 =S1''[1] == 0;   
eq3 =     S1[2] ==  1;   
eq4 =     S1[2] == S2[2];   
eq5 =   S1'[2] == S2'[2];   
eq6 =S1''[2] == S2''[2];   
eq7 =     S2[3] == 1;   
eq8 =S2''[3] == 0;   
eqns = { eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8 };   
vars = { (Subscript[a, 1]),(Subscript[b, 1]),(Subscript[c, 1]),(Subscript[d, 1]),(Subscript[a, 2]),(Subscript[b, 2]),(Subscript[c, 2]),(Subscript[d, 2]) };   
solset = Solve[eqns,vars];   
Print[TableForm[eqns]];   
Print[solset];

Use the solution given above for the coefficients and form the cubic functions.  Remember that we must dig out one set of braces using  (Subscript[solset, [[1]]])  before we can use the ReplaceAll command.

g1[x_] = ReplaceAll[ S1[x],(Subscript[solset, [[1]]]) ];   
g2[x_] = ReplaceAll[ S2[x],(Subscript[solset, [[1]]]) ];   
Print["(Subscript[g, 1])[x] = ",g1[x] ];   
Print["(Subscript[g, 2])[x] = ",g2[x] ];

Now graph the portion of each cubic in the interval over which it is to be used. Then combine the two piecewise cubic graphs to form the spline.

(xy={{1,g1[1]},{2,g1[2]},{3,g2[3]}};) (dots=ListPlot[xy,PlotStyle->{Red,PointSize[0.02`]},DisplayFunction->Identity];) (graph1=Plot[g1[x],{x,1,2},PlotRange->{{0,3},{0,2.2`}},PlotStyle->Green,DisplayFunction->Identity];) (graph2=Plot[g2[x],{x,2,3},PlotRange->{{0,3},{0,2.2`}},PlotStyle->Blue,DisplayFunction->Identity];) (graph12=Show[graph1,graph2,dots,AspectRatio->((2.2`)/3),DisplayFunction->$DisplayFunction]) (Print["\!\(\*SubscriptBox[\(g\), \(1\)]\)[x] = ",g1[x]];) (Print["\!\(\*SubscriptBox[\(g\), \(2\)]\)[x] = ",g2[x]];)

Remark.  It would be nice to have one piecewise cubic function S[z] that is used. The following formulas for  S[x]  uses the condition syntax  /;  for making a piecewise function.  Under the help menu we can find the following information about Condition.
patt  /;  test is a pattern which matches only if the evaluation of test yields True

Clear[S];  
S[x_/;1≤x≤2] = g1[x];   
S[x_/;2≤x≤3] = g2[x];

Caution.  We cannot print the formula for a piecewise function with the Print command.  It is only possible to interrogate the system and determine what rule it is storing for  S.  Hence, it was a good idea to always use two formulas g1 and g2 to define  S, since they can be printed.

?S

When a plot is made, it only uses real numbers specified in the conditions.

(graph=Plot[S[x],{x,1,3},PlotRange->{{0,3},{0,2.2`}},PlotStyle->Magenta,AspectRatio->((2.2`)/3),DisplayFunction->Identity];) (Show[graph,dots,DisplayFunction->$DisplayFunction]) (Print["S[x_/;1≤x≤2] = ",g1[x]];) (Print["S[x_/;2≤x≤3] = ",g2[x]];)

Notice that the natural cubic spline is different from the clamped cubic spline, it is "a relaxed curve." (and happy too!)

Your Job, Should You Choose to Accept It (and it is homework, after all!):

	Now, having enjoyed your romp through a variety of splines, your job is to use Mathematica to do two things:
	1. Create a Hermite cubic spline of two cubics that fit the function f(x)=sin(x) and its derivative at the three x-values 0, Pi/4, Pi/2. You can do this via linear algebra, or "the hard way". In any event, you can use some of the tips learned via the procedures used for the clamped and natural cubic splines. Graph the function f and its interpolator on this interval, as well as their difference (and hence the error).
	2. Graph the derivative of the spline above. What do you expect it to look like?
	3. Construct the natural cubic spline that does the same (that is, repeating steps 1 and 2).
Good luck!