For Shira's Matlab version, click here.

function inter(X,Y) // find size of data, minimum and maximum x n=size(X,2); least=X(1); most=X(1); for i=2:n if (X(i)<least) least=X(i) end if (X(i)>most) most=X(i) end end // make table of divided differences table=zeros(n,n); table(:,1)=Y'; for i=2:n for j=i:n table(j,i)=(table(j,i-1)-table(j-1,i-1))/(X(j)-X(j-i+1)); end end // make graph of polynomial using Horner's method and plot x=[0:100]*(most-least)/100+least; y=table(n,n); for i=1:(n-1) y=y.*(x-X(n-i))+table(n-i,n-i); end xbasc(0) plot2d(X,Y,style=-4) plot2d(x,y,style=4) endfunction

Here are some examples of runs:

X=[10 9 8 4 2 3] // the x-coordinates need not be ordered Y=[0.1 0.3 0.4 0.9 0.8 0.5] inter(X,Y)

X=[0 1 2 3 4 5 6] Y=[1 1 1 1 1 1 1] // straight line output expected inter(X,Y)

X=[0 1 2 3 4 5 6] Y=[1 1.1 1 1.1 1 1.1 1] inter(X,Y)

Polynomial interpolation can give highly oscillatory results if the x-coordinates are equally spaced apart. The classic example of this is known as "Runge's example". We try to do polynomial interpolation of equally spaced points on the curve y=1/(1+x^2). Below are results using 11 and 21 points:

X=[-5:5] Y=(X.*X+1).^(-1) inter(X,Y)

X=[-10:10] Y=(X.*X+1).^(-1) inter(X,Y)

To cure this we need to squash the points up a bit towards the end of the ends of the interpolation interval. Here is a better way to do it, using 11 and 21 respectively:

X=5*sin(%pi*[-5:5]/10) Y=(X.*X+1).^(-1) inter(X,Y)

X=10*sin(%pi*[-10:10]/20) Y=(X.*X+1).^(-1) inter(X,Y)