% Perform Newton Interpolation on function func
% evaluated at points x(1:n)
% function [c,err,NewtonTable,NewtonTableErr,efactors] = NewtonInterpEbnd(func,x,plt)
%
% Inputs:
% func = name of function to be interpolated
% x = array of n points at which to perform interpolation
% x MUST BE SORTED IN INCREASING ORDER
% plt = 1 to produce graphics, 0 not to
%
% Outputs:
% c = n+1 coefficients of interpolation polynomial
% p(x) = sum from i=0 to n c(i+1) * (x-x(1))*...*(x-x(i))
% err = maximum difference |p(x)-func(x)| at
% 1000 equally spaced points in interpolation interval
% NewtonTable = divided difference table (see page 355 of text)
% NewtonTableErr = bounds on accumumlated roundoff in NewtonTable
% efactors = factors in error bound:
% efactors(1) = estimate of max f^(n)(x)/n!
% efactors(2) = estimate roundoff error in efactor(1)
% efactors(3) = 1/n!
% if plot = 1,
% Figure 1: plot of func(x) (blue), p(x) (red) at 1000 equally space points
% Figure 2: plot of actual error |f(x)-p(x)| (blue dotted)
% Newton error estimate (black dashed)
% = (efactor(1)+efactor(2)) * prod {i=0 to n} (x-x{i})
% Newton error bound factor (magenta)
% = 1/n! prod {i=0 to n} (x-x(i))
% roundoff bound from evaluating p(x) (red)
% Error bound based on estimated roundoff error
% NewtonTableErr in computing NewtonTable (magenta)
%
function [c,err,NewtonTable,NewtonTableErr,efactors] = NewtonInterpEbnd(func,x,plt)
xmin = min(x);
xmax = max(x);
xx = xmin:(xmax-xmin)/999:xmax;
n = length(x);
x = reshape(x,n,1);
%
nold = n;
n = n + 10;
% Compute Divided Difference Table
% NewtonTable = zeros(n,n+1);
NewtonTable = zeros(n,nold+2);
% NewtonTableErr = zeros(n,n+1);
% NewtonTable(1:n,1) = reshape(x,n,1);
NewtonTable(1:nold,1) = x;
NewtonTable(nold+1:n,1) = rand(10,1)*(xmax-xmin) + xmin;
NewtonTable(1:n,2) = eval([func,'(NewtonTable(1:n,1))']);
NewtonTableErr(1:n,2) = eps*abs(NewtonTable(1:n,2));
% for j=1:n-1,
for j=1:nold,
NewtonTable(1:n-j,j+2) = ...
(NewtonTable(2:n-j+1,j+1) - NewtonTable(1:n-j,j+1)) ./ ...
(NewtonTable(j+1:n,1) - NewtonTable(1:n-j,1));
NewtonTableErr(1:n-j,j+2) = ...
( eps*abs(NewtonTable(2:n-j+1,j+1) - NewtonTable(1:n-j,j+1)) + ...
NewtonTableErr(2:n-j+1,j+1) + NewtonTableErr(1:n-j,j+1) ) ./ ...
(NewtonTable(j+1:n,1) - NewtonTable(1:n-j,1));
end
maxfn = max(abs(NewtonTable(:,nold+2)));
maxfnerr = max(abs(NewtonTableErr(:,nold+2)));
efactors = [maxfn,maxfnerr,1/gamma(nold+1)];
n = nold;
c = NewtonTable(1,2:n+1);
cbnd2f = NewtonTableErr(1,2:n+1);
cbnd2b = diag(NewtonTableErr(n:-1:1,:),1);
%
% Evaluate Newton Interpolation polynomial at xx, compute error,
% bound error from roundoff in forming and evaluating the
% Newton polynomial
p = c(n)*ones(size(xx));
pbnd1 = eps*abs(p);
pbnd2f = cbnd2f(n);
pbnd2b = cbnd2b(n);
for i = n:-1:2,
p = p.*(xx-x(i-1)) + c(i-1);
pbnd1 = pbnd1.*abs(xx-x(i-1)) + eps*abs(c(i-1));
pbnd2f = pbnd2f.*abs(xx-x(i-1)) + cbnd2f(i-1);
pbnd2b = pbnd2b.*abs(xx-x(n+1-(i-1))) + cbnd2b(i-1);
end
fxx = eval([func,'(xx)']);
fx = eval([func,'(x)']);
err = max(abs(fxx-p));
%
% Evaluate error due to roundoff in computing divided differences
%
kk=1;
% loop through all knots
for i = 1:n
% look at xx just less than know if j=1
% look at xx just greater than know if j=n
for j = 1:2
if ((i ~= 1 | j ~= 1) & (i ~= n | j ~= n)),
if j==1,
mn = (x(i-1)+x(i))/2; mx = x(i);
else
mn = x(i); mx = (x(i)+x(i+1))/2;
end
% compute range [kold,k] of xx() in range [mn,mx]
kold = kk+1; if (kold==2), kold =1; end
while(xx(kk)<mx), kk=kk+1; end
%% disp('i,j,kold,kk='), [i,j,kold,kk]
% sort indices of knots in order of increasing distance
% selected xx()
[dist,loc]=sort(abs(x-(mn+mx)/2));
%% disp('loc='),loc
% collecte divided differences in this order
cbndloc(1)=loc(1);
minl = loc(1);maxl=loc(1);
for k=2:n,
if loc(k)<minl,
cbndloc(k)=cbndloc(k-1)-1; minl=loc(k);
else
cbndloc(k)=cbndloc(k-1); maxl=loc(k);
end
end
%% disp('cbndloc='),cbndloc
for k=1:n, cbnddde(k)=NewtonTableErr(cbndloc(k),k+1); end
% evaluate error
pbnd2(kold:kk) = cbnddde(n);
for i = n:-1:2,
% p = p.*(xx-x(i-1)) + c(i-1);
% pbnd1 = pbnd1.*abs(xx-x(i-1)) + eps*abs(c(i-1));
% pbnd2f = pbnd2f.*abs(xx-x(i-1)) + cbnd2f(i-1);
% pbnd2b = pbnd2b.*abs(xx-x(n+1-(i-1))) + cbnd2b(i-1);
pbnd2(kold:kk) = pbnd2(kold:kk).*abs(xx(kold:kk)-x(loc(i-1))) ...
+ cbnddde(i-1);
end
%% disp('pause'),pause
end
end
end
%% size(xx),size(pbnd2),size(pbnd2b),size(pbnd2f)
%
% Compute factor of standard error bound for Newton Interpolation:
% 1/n! * prod (i=1 to n) (x - x(i))
%
errfactor = ones(size(xx))/gamma(n+1);
for i=1:n, errfactor = errfactor .* (xx-x(i)); end
errfactor = abs(errfactor);
% Plot
if (plt == 1)
figure(1)
hold off, clf
plot(xx,fxx,'b.',xx,p,'r',x,fx,'bx')
grid
title('f(x) = blue, (xi,f(xi)) = blue x, p(x) = red')
xlabel(['Maximum error = ',num2str(err)])
figure(2)
hold off, clf
semilogy(xx,abs(fxx-p),'b.',xx,pbnd1,'r',xx,errfactor,'m', ...
xx,(maxfn+maxfnerr)*errfactor*gamma(n+1),'k--',xx,pbnd2,'green'),
grid
title('|f(x) - p(x)| = blue (dotted), Newton error estimate = black dashed, Newton error factor = magenta')
xlabel('roundoff in p(x) = red, NewtonTable error = green')
v=axis;
v(1) = xmin; v(2) = xmax; v(3) = 1e-20*v(4); axis(v);
end