Thiele's interpolation formula

In mathematics, Thiele's interpolation formula is an interpolation formula for a function "f"(•) of a single variable, due to the Danish mathematician Thorvald Nicolai Thiele. It is expressed as a continued fraction, where ρ represents the reciprocal difference:

:$f(x)\; =\; f(x\_1)\; +\; cfrac\{x-x\_1\}\{\; ho(x\_1,x\_2)\; +\; cfrac\{x-x\_2\}\{\; ho\_2(x\_1,x\_2,x\_3)\; -\; f(x\_1)\; +\; cfrac\{x-x\_3\}\{\; ho\_3(x\_1,x\_2,x\_3,x\_4)\; -\; ho(x\_1,x\_2)\; +\; cdots\}$

Here is a version of Thiele interpolation in Algol 68:

¢ The MODE of lx and ly here should really be a UNION of "something REAL", "something COMPLex", and "something SYMBOLIC" ... ¢ PROC thiele:=( [] REAL lx,ly, REAL x) REAL: BEGIN [] REAL xx=lx [@1] ,yy=ly [@1] ; INT n=UPB xx; IF UPB yy=n THEN ¢ Assuming that the values of xx are distinct ... ¢ [0:n-1,1:n] REAL p; p [0,] :=yy [] ; FOR i TO n-1 DO p [1,i] :=(xx [i] -xx [1+i] )/(p [0,i] -p [0,1+i] ) OD; FOR i FROM 2 TO n-1 DO FOR j TO n-i DO p [i,j] :=(xx [j] -xx [j+i] )/(p [i-1,j] -p [i-1,j+1] )+p [i-2,j+1] OD OD; REAL a:=0; FOR i FROM n-1 BY -1 TO 2 DO a:=(x-xx [i] )/(p [i,1] -p [i-2,1] +a) OD; y [1] +(x-xx [1] )/(p [1,1] +a) ELSE error ¢ Unequal length arrays supplied ¢ FI END;

Note that, although it works in most cases, it is sensitive to input values, especially those due to equally spaced abscissæ. (Essentially, in such a case, it may reduce to the ratio of two polynomials, which may have factors in common - thus yielding a 0/0 form.)

References

*