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