(10U
Thiele's Reciprocal Differences and Derivatives W. Kahan Feb. 2002
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
For any real function f(x) define ž[-1]f := 0 , ž[0]f(x) := f(x) ,
and for k = 1, 2, 3, ... in turn define the Reciprocal Difference
ž[k]f(x[0],x[1], ...,x[k]) :=
x[0] - x[k]
= ------------------------------------------------- +
ž[k-1]f(x[0],...,x[k-1]) - ž[k-1]f(x(1),...,x[k])
+ ž[k-2]f(x[1],...,x[k-1])
The value of ž[k]f(x[0],...,x[k]) turns out to be independent of the
order of distinct arguments except that some orders may produce NaN
from ì - ì . Usually the best order for x[0], x[1], x[2], ... is
monotonic. If now we abbreviate ž[k]f(x[0],...,[k]) = žf(0k) we
find that, with rare exceptions, the Continued Fraction
z - x[0]
žf(00) + -------- z - x[1]
žf(01) + -------- z - x[2]
žf(02)-žf(00) + -------- z - x[3]
žf(03)-žf(01) + --------
žf(04)-...
interpolates f(z) at z = x[0], x[1], x[2], ... . If f(z) is a
rational function, the continued fraction must terminate upon some
infinite žf(0k) although an accidental ì or NaN may interfere
earlier. For example take f(x) := 1/(1 + xý) ;
1/(1 + xý) = 1 + x/( -2 + (x-1)/( -« + (x+1)/( -2 + (x-2)/(-«) )))
so ž[3]f(-1, 0, 1, 2) = (-2) + (-2) = -4 , but actual computation
yields ž[3]f( 0, 1, -1, 2) = ì - ì , a NaN instead of -4 . Thus,
even if a rational function f(x) has a continued fraction like the
one above that interpolates it at z = x[0], x[1], x[2], ... ( and
there is no guarantee of that ), Thiele's reciprocal differences may
fail to deliver it. Sad. This risk is not mentioned by Maple V in
its help thiele ( readlib(thiele) first ).
Here is a Matlab program to compute reciprocal differences:
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:
% tf(i,j+1) = (RDiff^j) F({x(i), x(i+1), ..., x(i+j)}) for j = 0:n .
% If every f(i) = F(x(i)) then 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), ... .
% If F(z) is a rational function then some column of tf becomes all
% constant, possibly infinite, unless NaNs intervene. The elements
% of x(:) must be distinct lest divisions 0/0 occur. The value of a
% reciprocal difference would be independent of the order of the elements
% in the set {x(i), x(i+1), ..., x(i+j)} but for roundoff and NaNs.
% If real, tf is least likely to malfunction when array x is in sorted
% (monotonic) order. Undefineable entries of tf are left as NaNs.
% W. Kahan, 31 Jan. 2002
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:min(n,nx-1) ,
z = (x(1+j:nx)-x(1:nx-j))./diff(tf(1:nf+1-j,j)) + tf(2:nf+1-j,j-1) ;
tf([1:nx-j]',j+1) = z ;
end
.......................................................................
For a competent but brief account of rational interpolation and some of
its hazards, see ch. 2.2 of "Introduction to Numerical Analysis"
by J. Stoer & R. Bulirsch, translated by R. Bartels, W. Gautschi &
C. Witzgall (1980, Springer, New York) and citations therein.
.......................................................................
Analogous to reciprocal differences are reciprocal derivatives defined
thus: Given any sufficiently differentiable function f(x) , define
f[-1](x) := 0 , f[0](x) := f(x) , and for k = 1, 2, 3, ... in turn
f[k](x) := k/( d f[k-1](x) /dx ) + f[k-2](x) .
For example, when f[0](x) = 1/(1 + xý) we find
f[1](x) = -«(1 + xý)ý/x ,
f[2](x) = 1/(1 - 3xý) ,
f[3](x) = 4x(xý - 1)
f[4](x) = 0 .
In general, a rational function's reciprocal derivatives terminate at
some constant f[k](x) , but this is not easy to prove. If again we
abbreviate f0 := f[0](0) , f1 := f[1](0) , f2 := f[2](0) , ... then
with rare exceptions (when some fj is infinite prematurely) ...
x
f(x) = f0 + ---- x
f1 + --------- x
f2 - f0 + --------- x
f3 - f1 + --------- x
f4 - f2 + --------- x
f5 - f3 + ---------
f6 - f4 + ...
provides a Pad‚ continued fraction expansion of f(x) which matches
its Taylor series for all sufficiently tiny x but usually converges
over a bigger domain. Other ways exist to compute terminating Pad‚
continued fraction expansions that match the first several terms of a
Taylor series, but this way is convenient for functions f(x) that
satisfy certain simple differential equations or recurrences. For
example,
x
+û(1+x) = 1 + --- x for all x > -1
2 + --- x
2 + --- x
2 + ---
2 + ...
is convenient. So is
x
exp(x) = 1 + --- x
1 + --- x
-2 + --- x
-3 + --- x
2 + --- x
5 + --- x
-2 + --- x
-7 + --- x
2 + --- x
9 + --- x
-2 + ----
-11 + ... .
But all too often the complexity of reciprocal derivatives overwhelms
computerized algebra systems like DERIVE, for which code is appended:
" Reciprocal derivatives: f[-1]=0, f[0]=f, f[k]=k/f[k-1]'+f[k-2] ."
" RDIF(f,x,n) = [ f[0], f[1], ..., f[n] ] "
RDIF(f,x,n):=ELEMENT(ITERATES([ELEMENT(v_,1)+1,~
ELEMENT(v_,1)/DIF(ELEMENT(v_,2),x)+ELEMENT(v_,3),~
ELEMENT(v_,2)], v_,[1,f,0],n)`,2)
This code cannot cope with important instances like x/(1 + xý) and
x
tan(x) = --- xý
1 - --- xý
3 - --- xý
7 - --- xý
9 - --- xý
11 - ---
13 - ...
that lack Pad‚ continued fraction expansions with just x in every
numerator.