Neville and Newton Interpolation

Save this as neville.m:

function p = neville(x,f,xx)
%
%  function p = neville(x,f,xx)
%
%  implements the interpolation algorithm of Neville
%
%  on entry:
%    x       abscisses, given as a column vector;
%    f       ordinates, given as a column vector;
%    xx      point where to evaluate the interpolating
%            polynomial through (x(i),f(i)).
%
%  on return:
%    p       last row of Neville's table where p(1) is
%            the value of the interpolator at xx.
%
%  note:
%    MATLAB does not support negative or zero indices,
%    so this algorithm differs from the one on paper.
%
%  example (pages 227-228 in Gerald-Wheatley):
%    x = [32.0 22.2 41.6 10.1 50.5]';
%    f = [0.52992 0.37784 0.66393 0.17537 0.63608]';
%    xx = 27.5;
%    p = neville(x,f,xx)
%
n = size(x,1);
p = f;
for i=1:n
  dx(i) = xx - x(i);
end;
for i=2:n
  for j=2:i
    p(i-j+1) = (dx(i)*p(i-j+1) - dx(i-j+1)*p(i-j+2))/(x(i-j+1)-x(i));
  end;
end;

Here is the output of the example from the text book:

>> x = [32.0 22.2 41.6 10.1 50.5]'

x =

  32.00000000000000
  22.20000000000000
  41.60000000000000
  10.10000000000000
  50.50000000000000

>> f = [0.52992 0.37784 0.66393 0.17537 0.63608]'

f =

   0.52992000000000
   0.37784000000000
   0.66393000000000
   0.17537000000000
   0.63608000000000

>> xx = 27.5

xx =

  27.50000000000000

>> p = neville(x,f,xx)

p =

   0.45753649919172
   0.47901170612965
   0.55843197307291
   0.37379460396040
   0.63608000000000

We can also test it on a random polynomial:

>> c = rand(1,6);      % random polynomial of degree five
>> x = rand(5,1);      % random abscisses
>> f = polyval(c,x)';  % evaluate polynomial
>> xx = rand(1,1);     % a random test point
>> y = polyval(c,xx)

y =

   0.93419464936959

>> p = neville(x,f,xx)'

p =

   0.93430718939273
   0.93882120509096
   0.96275050937465
   0.58504916015771
   0.19922880542495

>> p(1) - y

ans =

     1.125400231457574e-04

Since the interpolation problem has a unique solution with distinct abscisses, the interpolation with five points should lead to exactly the same five degree polynomial.

Save this as the file divdif.m:

function d = divdif(x,f)
%
%  function d = divdif(x,f)
%
%  returns the vector of divided differences for the
%  Newton form of the interpolating polynomial
%
%  on entry:
%    x       abscisses, given as a column vector;
%    f       ordinates, given as a column vector.
%
%  on return:
%    d       divided differences
%
%  note:
%    MATLAB does not support negative or zero indices,
%    so this algorithm differs from the one on paper.
%
%  example (page 234 in Gerald-Wheatley):
%    x = [1.1 2 3.5 5 7.1]';
%    f = [0.6981 1.4715 2.1287 2.0521 1.4480]';
%    d = divdif(x,f)
%
n = size(x,1);
d = f;
for i=2:n
  for j=1:i-1
    d(i) = (d(j) - d(i))/(x(j) - x(i));
  end;
end;

Following is the output on the example from the text book:

x = [1.1 2 3.5 5 7.1]'

x =

   1.10000000000000
   2.00000000000000
   3.50000000000000
   5.00000000000000
   7.10000000000000

>> f = [0.6981 1.4715 2.1287 2.0521 1.4480]'

f =

   0.69810000000000
   1.47150000000000
   2.12870000000000
   2.05210000000000
   1.44800000000000

>> d = divdif(x,f)

d =

   0.69810000000000
   0.85933333333333
  -0.17550000000000
   0.00318803418803
   0.00264985196358