Interpolating Polynomials

Suppose we are given n+1 points:

  ((Subscript[x, 0]),(Subscript[y, 0])),  ((Subscript[x, 1]),(Subscript[y, 1])),  ((Subscript[x, 2]),(Subscript[y, 2])), . . . ,  ((Subscript[x, n]),(Subscript[y, n])).

The idea is to find the unique polynomial of degree n:

  p(x)  =  (Subscript[a, 0]) + (Subscript[a, 1]) x + (Subscript[a, 2]) (x^2) + ... + (Subscript[a, n]) (x^n)

such that p((Subscript[x, i])) = (Subscript[y, i]) for all i.

Let us consider such a list.

list={{1,3},{2,5},{4,2},{6,7}};
plot1=ListPlot[list,PlotStyle->{Red,PointSize[0.05]},PlotRange->All]

Method 1 - Brute Force

The first method is to solve for the coefficients.

Clear[p,a]
n=Length[list]-1;
Do[x[k]=list[[k+1,1]];y[k]=list[[k+1,2]],{k,0,n}]
p[x_]=(∑_(k=0)^n)a[k] (x^k)

equations=Table[p[x[k]]==y[k],{k,0,n}]

solutions=Solve[equations]

p1[x_]=p[x]/.solutions[[1]]

plot2=Plot[p1[xx],{xx,x[0]-1,x[n]+1}];
Show[plot1,plot2]

Method 2 - Lagrange Polynomials

Here we use Lagrange Polynomials.

n=Length[list]-1;
Do[x[k]=list[[k+1,1]];y[k]=list[[k+1,2]],{k,0,n}]
Do[Lag[k,xx_]=((((∏_(j=0)^(k-1))(xx-x[j])) ((∏_(j=k+1)^n)(xx-x[j])))/(((∏_(j=0)^(k-1))(x[k]-x[j])) ((∏_(j=k+1)^n)(x[k]-x[j])))),{k,0,n}]
Plot[Evaluate[Table[Lag[k,xx],{k,0,n}]],{xx,x[0]-1,x[n]+1},
	          PlotStyle->Table[RGBColor[(k/n),0,1-(k/n)],{k,0,n}]]

p2[x_]=(∑_(k=0)^n)y[k] Lag[k,x]

Let's check that the two methods agree:

Simplify[p1[x]-p2[x]]

Method 3 - Use Mathematica

p3[x_]=InterpolatingPolynomial[list,x]

Simplify[p3[x]-p1[x]]

How good are interpolating polynomials?

Just how good are interpolating polynomials at "joining the dots"?  Let us illustrate this with two examples.  In the first one, we will try to match

    f(x) = sin(x).

list=Table[{x,Sin[x]},{x,-3.0,3.0,0.5}];
plot1=ListPlot[list,PlotStyle->{PointSize[0.02]},DisplayFunction->Identity];
p[x_]=Simplify[InterpolatingPolynomial[list,x]];
plot2=Plot[{Sin[x],p[x]},{x,-6,6},PlotStyle->{Red,Blue}];
Show[plot1,plot2]

Not bad, eh?  The function seemed to fit even outside the range -3 to 3 !!!

Now lets suppose that the data is a bit "noisy", that is, let's add a small random componant to the data.

list=Table[{x,Sin[x]+Random[Real,{-0.1,0.1}]},{x,-3.0,3.0,0.5}];
plot1=ListPlot[list,PlotStyle->{PointSize[0.02]},DisplayFunction->Identity];
p[x_]=Simplify[InterpolatingPolynomial[list,x]];
plot2=Plot[{Sin[x],p[x]},{x,-3.03,3.03},PlotStyle->{Red,Blue}];
Show[plot1,plot2]

Not very good!!!

Now lets see how splines do.  To be honest, I don't really
know how mathematica does these splines.  But they make 
a valuable comparison.

p2[x_]=Interpolation[list][x];
plot2=Plot[{Sin[x],p[x],p2[x]},{x,-3,3},PlotStyle->{Red,Blue,Green}];
Show[plot1,plot2]

In the following example, the `noise' is created by floating point errors.

list=Table[{x,Sin[x]},{x,-3.0,3.0,0.07}];
p[x_]=Expand[InterpolatingPolynomial[list,x]];
Plot[{Sin[x],p[x]},{x,-3.03,3.03},PlotStyle->{Red,Blue},PlotRange->{-4,4}]