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@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[]