(* Simple polynomial interpolation *) (* To run on command line: math -run "<{Hue[1], Thickness[0.01]}, PlotLabel->"Original function", DisplayFunction -> Identity ] (* Set data points -- nodes. *) n=4; h = (b-a)/n; (* define data points as values of F[x]. Try n=2,...7 *) data = Table[N[{x, F[x]}], {x, a, b, h}] fit = InterpolatingPolynomial[data, x] plotdata = ListPlot[data, PlotLabel -> "Data points", Prolog-> AbsolutePointSize [5] ] plotfit = Plot[fit, {x, a, b}, PlotStyle->{Hue[0.5], Thickness[0.005]}, PlotLabel-> "interpolation", DisplayFunction -> Identity ] Show[plotfunction, plotdata, plotfit, PlotLabel -> FontForm["red - original f(x) , blue - interpolation ", {"Courier", 15}], Prolog-> AbsolutePointSize[5], DisplayFunction -> $DisplayFunction ] Plot[F[x], {x, a, b}, PlotStyle->{Hue[1], Thickness[0.01]}, PlotLabel->"Original function to be interpolated"] Print[] Print[] Print[ "Interpolating polynomial Pn = " fit] Print[] Print[" Or, in simplified form, Pn = " Expand[Simplify[fit]]] Print[] (* So, what is the minimal degree of Pn neccessary for a decent interpolation of this function? If it has 2 extrema, then derivative=0 at two points => (Pn)' = P2, and so the minimum power is P3. Try P4 and look at the breakdown of the coefficients. How large is the coefficient infront of x^4? What does that mean? The general principle is "Like likes like". So, for the above function, do you want to try odd or even degree polynomials? *) Example II comment end *) (* As is with any method, one can run into trouble *) (* Example III comment begins *) F[x_] := 1/(1 + 25*x*x); (* Runge's function *) a = -1.0; b= 1.0; (* range *) plotfunction = Plot[F[x], {x, a, b}, PlotStyle->{Hue[1], Thickness[0.01]}, PlotLabel->"Original function", DisplayFunction -> Identity ] n=20; h = (b-a)/n; (* define data points as values of F[x] *) (* data = Table[N[{x, F[x]}], {x, a, b, h}] *) data = Table[N[{0.5*((a+b) + (b-a)*Cos[(i/n)*Pi]), F[0.5*((a+b) + (b-a)*Cos[(i/n)*Pi])] } ], {i, 0, n}] fit = InterpolatingPolynomial[data, x] Print[" Interpolating function, Pn = " Expand[Simplify[fit]]] plotdata = ListPlot[data, PlotLabel -> "Data points", Prolog-> AbsolutePointSize[5] ] plotfit = Plot[fit, {x, a, b}, PlotStyle->{Hue[0.5], Thickness[0.005]}, PlotLabel-> "interpolation", DisplayFunction -> Identity ] Show[plotfunction, plotdata, plotfit, PlotLabel -> FontForm["red - original f(x), blue - interpolation ", {"Courier", 15}], Prolog-> AbsolutePointSize[5], DisplayFunction -> $DisplayFunction ] *) (* Possible fixes below *) (* Just cut and paste into the appropriate place in the above example *) (* ------------------------------------------------------- *) Example III. End comments *) (* Fix 1. Increase n *) (* Generally BAD idea *) (* Fix 2. Extend range of data. Approximates well if deep inside max range *) (* Fix 3. Use non-uniform or "unequally spaced" data points ("nodes", "sampling") *) (* e.g. n Chebyshev points xi = cos[Pi*(2i +1)/(2n +2)], i <= n*) (* data = Table[N[{0.5*((a+b) + (b-a)*Cos[(i/n)*Pi]), F[0.5*((a+b) + (b-a)*Cos[(i/n)*Pi])] } ], {i, 0, n}] *) (* Fix 4. Use more appropriate MODEL (basis functions). Gaussians in OUR case *) (* fit = Fit[data, Table[Exp[-i*1.0*x^2], {i, 1, n}], x] *) Clear[data] Clear[fit] (* Thanks to Niteesh Bharara Fall 2003 for fixing the multiple plot problem *)