type newtdd

%Program 3.1 Newton Divided Difference Interpolation Method
%Computes coefficients of interpolating polynomial
%Input: x and y are vectors containing the x and y coordinates
%       of the n data points
%Output: coefficients c of interpolating polynomial in nested form
%Use with nest.m to evaluate interpolating polynomial
function c=newtdd(x,y,n)
for j=1:n
  v(j,1)=y(j);        % Fill in y column of Newton triangle
end
for i=2:n             % For column i,
  for j=1:n+1-i       % fill in column from top to bottom
    v(j,i)=(v(j+1,i-1)-v(j,i-1))/(x(j+i-1)-x(j));
  end
end
for i=1:n
  c(i)=v(1,i);        % Read along top of triangle
end                   % for output coefficients


type nest.m

%Program 0.1 Nested multiplication
%Evaluates polynomial from nested form using Horner's method
%Input: degree d of polynomial,
%       array of d+1 coefficients (constant term first),
%       x-coordinate x at which to evaluate, and
%       array of d base points b, if needed
%Output: value y of polynomial at x
function y=nest(d,c,x,b)
if nargin<4, b=zeros(d,1); end 
y=c(d+1);
for i=d:-1:1
  y = y.*(x-b(i))+c(i);
end


n=4;
x0=(pi/2)*(0:(n-1))/(n-1);
y0=sin(x0);
c=newtdd(x0,y0,n);
x=-1:.01:3;
y=nest(n-1,c,x,x0);
c

c =

         0    0.9549   -0.2443   -0.1139

plot(x,y,'r',x0,y0,'go',x,sin(x),'b:')
grid
figure(2);semilogy(x,abs(y-sin(x)))
n=8;
x0=(pi/2)*(0:(n-1))/(n-1);
y0=sin(x0);
c=newtdd(x0,y0,n);
y=nest(n-1,c,x,x0);
figure(1)
plot(x,y,'r',x0,y0,'go',x,sin(x),'b:')
grid
figure(2);semilogy(x,abs(y-sin(x)))
grid
sin(1)

ans =

    0.8415

format long
sin(1)

ans =

   0.84147098480790

y=nest(n-1,c,1,x0);
y

y =

   0.84147099092917

c

c =

  Columns 1 through 3 

                  0   0.99162858425603  -0.11079437204152

  Columns 4 through 6 

  -0.15632638990518   0.01792733674552   0.00698236874874

  Columns 7 through 8 

  -0.00085512051116  -0.00013825261116

diary off