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:iThere 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.