Math 310, Fall 1999
Bueler

A Neville's algorithm example in Matlab


Neville's algorithm as given in the book is easy to implement in Matlab, if one realizes that there is a quick way to deal with the "0 vs. 1" indexing issue.  That is, arrays in other languages are frequently indexed from  i=0  to  i=n.  This is what the book assumes.  But since the Matlab language has the concept of a matrix so firmly imbedded, all indices start at   i=1.

So, to write Algorithm 3.1 in Matlab, one trick is to add one to indices, but keep the "for" loops the same.  That is, the only line that needs to change is the main line: 

        Q(i+1,j+1) = ((xx-x(i-j+1))*Q(i+1,j) - (xx-x(i+1))*Q(i,j))/...
                     (x(i+1)-x(i-j+1));
Again, note that the only change is that the indices in this line have all been increased by one.  The loops are still: 
  for i = 1:n
     for j = 1:i
There is a better way to write Algorithm 3.1.  Namely, to compute each new column of the Neville's algorithm table, so that the program doesn't have to store away the whole table as a matrix, but rather just a vector which is a column.  This is what the code below does.  There is a disadvantage--the program can not be modified to add nodes as needed, as suggested on page 119.  It should be saved as "nev.m": 
function y = nev(xx,n,x,Q)
% Neville's algorithm as a function (save as "nev.m")
% 
% inputs:
%    n = order of interpolation (n+1 = # of points)
%    x(1),...,x(n+1)    x coords
%    Q(1),...,Q(n+1)    y coords
%    xx=evaluation point for interpolating polynomial p
%
% output:  p(xx)

for i = n:-1:1
   for j = 1:i
      Q(j) = (xx-x(j))*Q(j+1) - (xx-x(j+n+1-i))*Q(j);
      Q(j) = Q(j)/(x(j+n+1-i)-x(j));
   end
end

y = Q(1);


Now, to run this on problem I. in assignment #2 would look like: 

» x=1940:10:1990

x =

        1940        1950        1960        1970        1980        1990

» Q=[132165 151326 179323 203302 226542 249633]

Q =

      132165      151326      179323      203302      226542      249633

» nev(1940,5,x,Q) % just a check

ans =

      132165

» nev(1930,5,x,Q)

ans =

      169649

» nev(2000,5,x,Q)

ans =

      251654

»
But to actually plot out some polynomials, as in problem II. on assignment #2, you probably want to actually write a new program, for instance:
  
% save in something like "probII.m"
%    --the name doesn't matter since it's not a function
% plot out the n=2, 4, 8, 16 polynomials which interpolate the 
%    function f(x) = 1/(1+25x^2) on the interval [-1,1] using 
%    equally spaced points

xx=[-1:.01:1];  % plot a couple hundred points in each case

% first plot the original function
yy=1./(1+25*xx.^2); 
plot(xx,yy);
hold on;               % this makes the plots that follow pile up


for n=[2 4 8 16]
   x=-1:2/n:1;         % (check this works for n=4, for instance)
   Q=1./(1+25*x.^2);   % get the interpolation points
   for i=1:201
      yy(i)=nev(xx(i),n,x,Q);  %interpolate
   end;
   plot(xx,yy);
end;

%  inessential embellishments:
axis([-1 1 -2 2]);     % experimentation showed this was nice
plot(x,Q,'o');         % show the 16 interpolation points--last values of x and Q
title('Polynomial interpolation of y=1/(1+25*x^2)--the Runge example'); 

hold off;              % so future plots start over


Run it--it produces a nice picture.

Another implementation of Neville's algorithm can be found in Numerical Recipes, where there is an attempt to estimate error within the algorithm: Section 3.1 of Numerical Recipes.