(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.