Package for Hermite Polynomial Interpolation
- To: mathgroup at yoda.physics.unc.edu
- Subject: Package for Hermite Polynomial Interpolation
- From: Gossett <gossett at bethel.edu>
- Date: Mon, 3 Aug 1992 10:47:25 -0500
This package does Hermite Polynomial Interpolation. It interpolates to a set
of function values and a set of first derivative values.
BeginPackage["Hermite`"]
HermiteP::usage = "HermiteP[Data_List,x_Symbol] or HermiteP[X_List,F_List,Fp_List,x_Symbol] produces a degree 2n Hermite polynomial in x, where Lenght[X] == n. It also plots the polynomial and the data. Data is a list of triples {x,f(x),f'(x)}. The alternate form takes three lists; the first contains the x values, the second contains the values f(x), and the third contains f'(x). Examples: HermiteP[{{1,2,3},{2,3,4},{3,4,5},{4,5,6}},x] or HermiteP[{1,2,3,4},{2,3,4,5},{3,4,5,6},x]."
(* This package may be freely used for non-commercial purposes. It has been
designed for educational purposes, but may be useful in other contexts.
Copyright June 1992
Dr. Eric Gossett
Bethel College
3900 Bethel Drive
St. Paul, MN 55112
gossett at bethel.edu
Version 1.0
This function is a procedural implementation of a divided difference algorithm. Reference: "Numerical Analysis" 4th edition, by Burden and Faires, PWS-Kent, 1989. It is quite slow for large data sets.
*)
Begin["`Private`"]
Unprotect[HermiteP]
HermiteP[Data_List,x_Symbol] :=
Module[{X,F,FP},
If[(Length[Dimensions[Data]]==2) && (Dimensions[Data][[2]]==3),
(* Then *)
X = Transpose[Data][[1]];
F = Transpose[Data][[2]];
FP = Transpose[Data][[3]];
HermiteP[X,F,FP,x],
(* Else *)
Print["Improper data: type ?HermiteP for proper format."]
]
]
HermiteP[X_List,F_List,Fp_List,x_Symbol] :=
Module[{i,j,n,z={},Q={},xprod,H},
n = Length[X];
If[(n != Length[F])||(n != Length[Fp]),
(* Then*)
Print["list sizes don't match"],
(* Else initialize 1st three columns for divided differences *)
For[i=1,i<=n,i++,
AppendTo[z,X[[i]]];
AppendTo[z,X[[i]]];
AppendTo[Q,{F[[i]]}];
AppendTo[Q,{F[[i]]}];
AppendTo[Q[[2i]],Fp[[i]]];
If[i != 1,
AppendTo[Q[[2i-1]],(Q[[2i-1,1]]-Q[[2i-2,1]])/(z[[2i-1]]-z[[2i-2]])],
];
]; (* end For i *)
(* Form the rest of the difference table *)
For[i=3,i<=2n,i++,
For[j=3,j<=i,j++,
AppendTo[Q[[i]],(Q[[i,j-1]]-Q[[i-1,j-1]])/(z[[i]]-z[[i-j+1]])];
]; (* end For j *)
]; (* end For i *)
(* form the polynomial *)
For[i=2;H=Q[[1,1]];xprod=1,i<=2n,i++,
xprod = xprod*(x - z[[i-1]]);
H = H + Q[[i,i]]*xprod;
]; (* end For i *)
(* plot the polynomial vs the data *)
Show[ListPlot[Transpose[{X,F}], PlotStyle->{PointSize[.015]},
DisplayFunction->Identity],
Plot[H,{x,z[[1]]-1,z[[2n]]+1}, DisplayFunction->Identity],
DisplayFunction->$DisplayFunction];
(* print the polynomial in divided difference form *)
Print[H];
(* retrun the function in standard form *)
Return[Expand[H]];
]; (* end If Lengths *)
]
End[]
Protect[HermiteP]
EndPackage[]