Jump to content

Talk:Thiele's interpolation formula

Page contents not supported in other languages.
From Wikipedia, the free encyclopedia

The article on Thiele's interpolation formula duplicates the form given in Abramowitz and Stegun's "Handbook of Mathematical Functions" - as does Eric Weisstein's. However, it is not particularly easy to implement, from the description given, even by writing recursively. In the (antique) version of "Maple" to which I have access, there is a good non-recursive working implementation of Thiele's algorithm. (Dognose (as Anthony Burgess wrote it in one novel), I've used it to good effect often enough.) If anyone is interested, I can pass on the code (copyright abuse permitting).

Hair Commodore 20:58, 2 November 2006 (UTC)[reply]

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" 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 reduces to the ratio of two polynomials, which may have factors in common - thus yielding a 0/0 form.)

Comments? Hair Commodore 18:00, 1 February 2007 (UTC)[reply]

Algol 68???!? Really?

Arghman (talk) 17:43, 27 June 2015 (UTC)[reply]

Algol -> Simula

Algol -> Ada

Algol -> Pascal -> C#

Not too bad :) 2001:4643:EBFE:0:A11B:257F:77A1:A7F (talk) 10:05, 23 March 2016 (UTC)[reply]