(* Simple polynomial interpolation *)

(* First, try something really simple that will surely work to demonstrate the 
basic concept *)

(* EXAMPLE I *)

(* 3 points on a parabola, {x0, f(x0)}, {x1, f(x1)}, {x2, f(x2)} *) 
(* f(x) = x^2 in this example *)
(* Don't forget to uncomment the appropriate part *)

(*  example I comment begin. Uncomment begin and end, comment out other examples


Poly[x_] = InterpolatingPolynomial[{ {-1, 1}, {0, 0}, {1, 1} }, x] 


Print[]
Print[]

Print[ "Interpolating polynomial Pn =   " Poly[x]] 
Print[]
Print[" Or, in simplified form, Pn = " Simplify[Poly[x]]]

Print[]

Plot[Poly[x], {x, -2, 2}] 

Example I comment end *)

(* EXAMPLE II What minimal power n of Pn we need to decently interpolate a 
smooth function? By "decently" I mean that the main features are reproduced, 
such as all maxs and mins. *) 

(* Example II comment begin 

F[x_] := Sin[x]; 

a = Pi/2; b = (7/2)*Pi; 

plotfunction = Plot[F[x], {x, a, b}, PlotStyle->{Hue[1], Thickness[0.01]}, PlotLabel->"Original function", DisplayFunction -> Identity ]

n=2; 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 *)



F[x_] := 1/(1 + x*x); (* Runge's function *) 


a = -5.0; b= 5.0; (* range *) 

plotfunction = Plot[F[x], {x, a, b}, PlotStyle->{Hue[1], Thickness[0.01]}, PlotLabel->"Original function", DisplayFunction -> Identity ]

n=2; h = (b-a)/n;  (* define data points as values of F[x] *)



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 ]  



(* Possible fixes below *)
(* Just cut and paste into the appropriate place in the above example *)

(* ------------------------------------------------------- *)









(* 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*0.1*x^2], {i, 1, n}], x] 

*)


 

Clear[data]

Clear[fit]

(* Thanks to Niteesh Bharara for fixing the multiple plot problem    *)