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) ;