Matlab programs for Interpolation and Extrapolation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function y = aitkend2(x)
% y = aitkend2(x) applies Aitken's Delta-Squared convergence-
% acceleration process to a column x(:) of consecutive elements
% of a presumably convergent sequence to obtain a column y(:)
% whose elements converge faster to x's limit L whenever the
% ratios of consecutive terms of x-L tend to a limit K ~= 1 .
% To diminish roundoff's contribution, exploit the identity
% y-L = aitkend2(x-L) using any estimate L of x's limit.
% Length(y) = length(x)-2 . The process may be repeated but may
% amplify errors in x by up to ((|K|+1)/|K-1|)^2 each time.
% (C) 5 Feb. 2002 W. Kahan
x = x(:) ; n = length(x)-1 ;
d = diff(x) ;
y = x(2:n) - d(1:n-1).*d(2:n)./diff(d) ;
function [df, ed] = divdiff(n,f,x,ef)
% df = divdiff(n,f,x) = array of divided differences of f(:) over x(:)
% from 0-th to the n-th divided differences. If f(:) = F(x(:)) ,
% df(i,j+1) = (Diff^j) F({x(i), x(i+1), ..., x(i+j)}) for j = 0:n
% = (Deriv^j) F( somewhere amidst x(i), ..., x(i+j) )/j! ;
% and df figures in a polynomial
% df(i,1) + (z-x(i))*(df(i,2) + (z-x(i+1))*(df(i,3) + (z-x(i+2))*(df(i,...
% that interpolates F(z) at z = x(i), x(i+1), x(i+2), x(i+3), ... .
% When F(z) is a k-th degree polynomial, df(i,j+1) = 0 for j > k ;
% but increasing j amplifies f(:)'s error ever more into df(:,j+1) .
% The elements of x(:) must be distinct lest divisions by zero occur.
% Except for roundoff, the value of df(i,j+1) should be independent of
% the order of the rows in [x(i:i+j), f(i:i+j)] ; when x is real the
% accuracy of df is best if x(:) is in sorted (monotonic) order.
% Undefineable entries of df are left as NaNs.
%
% [df, ed] = divdiff(n,f,x,ef) provides an optional array ed of error-
% bounds for the error df inherits from errors as big as ef in f ;
% but if optional input ef is omitted it is assumed to be ones(...)
% and then ed is an array of factors by which errors in f may be
% amplified. Bounds are valid ONLY for real x in a monotonic order.
% (C) 10 Feb. 1998, W. Kahan
x = x(:) ; f = f(:) ; nx = length(x) ; nf = length(f) ;
if nf~=nx,
lengthx = nx , lengthf = nf ,
error('divdiff(n,f,x) requires length(f(:)) = length(x(:)) .')
end
NaNs = NaN ;
df = [ f, NaNs(ones(nx,1),ones(1,n)) ] ;
for j = 1:min(n,nx-1) ,
m = nx-j ; k = m+1 ;
df([1:m],j+1) = diff(df(1:k, j)) ./ ( x(1+j:nx) - x(1:m) ) ;
end
if nargout > 1, % ... begin computation of ed :
if ~all(isreal(x)),
error('[df,ed] = divdiff(n,f,x,...) requires real x .')
end
if nargin < 4,
ef = ones(nf, 1) ;
else
ef = abs(ef(:)) ; lef = length(ef) ;
if (lef == 1), ef = ef(ones(nf,1)) ;
elseif (lef ~= nf),
lengthf = nf, lengthef = lef,
error('divdiff(n,f,x,ef) requires length(ef) = length(f) .')
end, end
m = cumprod(-ones(nf,1)) ; % ... = [-1 1 -1 1 -1 1 ...]' .
t = sort(x) ;
if ~(all(t==x)|all(flipud(t)==x)),
error('[df,ed] = divdiff(n,f,x,...) requires monotonic x .')
end
ed = abs( divdiff(n, ef.*m, x) ) ; % ... called recursively!
end
function tf = thldiff(n,f,x)
% tf = thldiff(n,f,x) = array of Thiele's reciprocal differences of
% f(:) over x(:) from 0-th to n-th reciprocal differences; if
% f(:) = F(x(:)) then, for j = 0:n ,
% tf(i,j+1) = ReciprocalDifference^j F( {x(i), x(i+1), ..., x(i+j)} )
% and tf figures in a continued fraction
% z - x(i)
% tf(i,1) + -------- z - x(i+1)
% tf(i,2) + ---------- z - x(i+2)
% tf(i,3)-tf(i,1) + ---------- z - x(i+3)
% tf(i,4)-tf(i,2) + ----------
% tf(i,5)-...
% that interpolates F(z) at z = x(i), x(i+1), x(i+2), x(i+3), ... .
% But increasing j amplifies f(:)'s error ever more into tf(:,j+1) .
% If F(z) is a rational function then a column of tf should become
% infinite, followed by NaNs (unless some occur earlier by accident);
% terminate F's continued fraction just before infinity. All elements
% of x(:) must be distinct lest divisions 0/0 occur. But for roundoff
% and accidental NaNs, tf(i,j+1) would be independent of the order of
% the rows of [x(i:i+j), f(i:i+j)] ; misbehavior is least likely when
% x(:) is in sorted (monotonic) order. Undefineable entries of tf are
% left as NaNs. (C) 3 Feb. 2002, W. Kahan
x = x(:) ; f = f(:) ; nx = length(x) ; nf = length(f) ;
if (nf ~= nx),
lengthx = nx , lengthf = nf ,
error('thldiff(n,f,x) requires vectors f and x of equal lengths.')
end
N = round(n) ;
if (N ~= n)|(n < 1)|(n >= nx),
N = n , lengthx = nx ,
error('thldiff(N,f,x) requires length(x) > integer N > 0 .')
end
NaNs = NaN ;
tf = [ f, NaNs(ones(nx,1),ones(1,n)) ] ;
tf([1:nx-1], 2) = diff(x)./diff(f) ;
for j = 2:n,
m = nx-j ; k = m+1 ;
tf([1:m],j+1) = (x(1+j:nx)-x(1:m))./diff(tf(1:k,j)) + tf(2:k,j-1) ;
end
function [y,ey] = polyintp(z, n, f, x, ef)
% y = polyintp(z, n, fa, xa) is an array of values at x = z of
% polynomials in x of degrees 1, 2, 3, ..., n that interpolate
% given values fa = f(xa) of a function f(x) at an array xa(:)
% of arguments x . Specifically, y(i,j) is the value at x = z
% of the polynomial of degree j (or less) in x that takes the
% values fa(k) at x = xa(k) for k = i:i+j . This works only
% for j = 1:min(n,length(fa)-i) , beyond which y(i,j) is NaN.
% Entries xa(:) must be distinct and, for best accuracy, in a
% monotonic (sorted) order. Using Neville's algorithm, polyintp
% would satisfy the identity y-Y == polyintp(z,n,fa-Y,xa) for any
% scalar Y but for roundoff, which can be attenuated by choosing
% any Y near enough the desired value(s) of y to diminish them.
% Values y(i,j) for different i agree to the extent that f(x)
% is approximated well near x = z by a polynomial of degree j ;
% but as j increases so do the effects of error in array fa(:) .
%
% [y,ey] = polyintp(z,n,fa,xa,ef) yields an optional array ey of
% error-bounds for the error y inherits from errors as big as ef
% in fa ; but if optional array ef is omitted it is assumed to be
% ones(...) and then ey is an array of factors by which errors in
% fa may be amplified. Bounds are valid ONLY if z is real and
% x is real and in a monotonic order. (C) 10 Feb. 1998, W. Kahan
x = x(:) ; f = f(:) ; lx = length(x) ; lf = length(f) ;
if (lf ~= lx), lengthx = lx, lengthf = lf,
error(' Polyintp(z,n,x,f) requires length(x) = length(f) .')
end
NaNs = NaN ; y = [f, NaNs(ones(lf,1), ones(1,n))] ;
for j = 1:n ,
i = 1:lx-j ;
y(i,j+1) = y(i,j) + (y(i+1,j)-y(i,j)).*(z - x(i))./(x(i+j)-x(i)) ;
end
y = y(:, 2:n+1) ;
if nargout > 1, % ... begin computation of ey :
if ~(all(isreal(x)) & isreal(z)),
error('[y,ey] = polyintp(z,n,f,x,...) requires real x and z.')
end
if nargin < 5,
ef = ones(lf, 1) ;
else
ef = abs(ef(:)) ; lef = length(ef) ;
if (lef == 1), ef = ef(ones(lf,1)) ;
elseif (lef ~= lf),
lengthf = lf, lengthef = lef,
error('polyintp(z,n,f,x,ef) needs length(ef) = length(f) .')
end, end
m = cumprod(-ones(lf,1)) ; % ... = [-1 1 -1 1 -1 1 ...]' .
t = sort(x) ;
if ~(all(t==x)|all(flipud(t)==x)),
error('[y,ey] = polyintp(z,n,f,x,...) requires monotonic x .')
end
m = m.*(2*(x>z) - 1) ;
ey = abs( polyintp(z,n,ef.*m,x) ) ; % ... called recursively!
end
function y = ratlintp(z, n, f, x)
% y = ratlintp(z, n, fa, xa) is an array of values at x = z of
% rational functions of x that interpolate given values fa = f(xa)
% of some function f(x) at ever more points of an array xa(:) of
% arguments x . Specifically, y(i,j) = Rj(z) for some rational
% function Rj(x) that takes values Rj(xa(k)) = fa(k) over the
% range k = i:i+j while j = 1:min(n,length(fa)-i) , beyond which
% y(i,j) is NaN. Rj(x) is intended to be a ratio Rj = Pj/Qj of
% polynomials with degree(Pj(x)) = floor(j/2) = j - degree(Qj(x)) ,
% but their degrees may turn out smaller. Entries xa(:) must be
% distinct and, for best results, sorted. Values y(i,j) for
% different i agree to the extent that f(x) is approximated well
% near x = z by a rational function Rj(x) of the kind intended,
% but as j increases so do the effects of errors in array fa(:) .
% And sometimes accidents introduce NaN spuriously into array y .
% Ratlintp(inf, nm xa, fa) is computed right but is rarely useful;
% and sometimes y(i,j) is useful for only even j or only odd j .
% Sometimes approximations by ratlintp(2/z, n, 2 ./xa, fa) or by
% 2 ./ratlintp(z, n, xa, 2 ./fa) etc. may come closer to f(z) .
% The formula used derives from ch. 2.2.3 of "Intro. to Numerical
% Analysis" by J. Stoer & R. Bulirsch. W. Kahan (C) 1 Feb. 2002
x = x(:) ; f = f(:) ; lx = length(x) ; lf = length(f) ;
if (lf ~= lx), lengthx = lx, lengthf = lf,
error(' Ratlintp(z,n,x,f) requires length(x) = length(f) .')
end
NaNs = NaN ; y = [zeros(lf,1), f, NaNs(ones(lf,1), ones(1,n))] ;
if (abs(z) == inf),
for j = 1:min(n,lx-1) ,
i = 1:lx-j ;
y(i, j+2) = y(i+1, j) ;
end
else
for j = 1:min(n,lx-1) ,
i = 1:lx-j ; i1 = 1:lx+1-j ;
q = (z - x(i+j)).*( y(i+1,j+1) - y(i+1,j) ) ;
q = q./( (z - x(i)) .* ( y(i,j+1) - y(i+1,j) ) - q ) ;
y(i,j+2) = y(i+1,j+1) + q.*diff(y(i1,j+1)) ;
end
end
y = y(:, 3:n+2) ;