Rational Function Interpolation for the HP-41

Overview
 

-The following program uses Thiele's formula:

-Given n data points  ( x₁ , y₁ ) , ( x₂ , y₂ ) , .................. , ( x_n , y_n ) ,   the rational function f(x) is evaluated by the continued fraction:

                       x - x₁     x - x₂      x - x₃                                       x - x_n-1
  f(x) =   y₁ +  ------    -------    --------   ...............................  ---------
                       r₁₂ +      r₁₂₃ +     r₁₂₃₄ +                                      r_12......n
 

  where  r_1k = ( x_k - x₁ )/( y_k - y₁ )  ,    r_12k = ( x_k - x₂ )/( r_1k - r₁₂ )  ,   r_123k = ( x_k - x₃ )/( r_12k - r₁₂₃ )  , .............. ,  r_12....n = ( x_n - x_n-1 )/( r_12...n - r_12...n-1 )
 

-Thus,  f(x) = p(x)/q(x)   where  deg(p) = deg(q)         if  n is odd
                                        and   deg(p) = deg(q) + 1   if n is even      ( the degrees may of course be smaller if some leading coefficients equal zero )

-However, if n is even, you may prefer:   deg(p) = deg(q) - 1
-In this case, simply set flag F02 and the above formula will be applied to 1/y
-But this will work only if all y-values are different from zero.

-"THIELE" uses synthetic registers M N O P Q  which can be replaced by any unused data registers.
 

Program Listing
 

Data Registers:           •  R00 = n                                                     ( Registers R00 thru R2n  are to be initialized before executing "THIELE"  )

                                      •  R01 = x₁    •  Rn+1 = y₁
                                      •  R02 = x₂    •  Rn+2 = y₂       R2n+1 = r₁₂
                                         ---------        ----------         R2n+2 = r₁₃     R3n = r₁₂₃
                                         ---------        ----------         ------------       ---------

                                      •  Rnn = x_n    •  R2n   = y_n        R3n-1 = r_1n     R4n-3 = r_12n  ---------  R(n²+3n)/2 = r_12.....n

Flags:  F01-F02-F10   Set flag F01 after executing the routine once: it will avoid to re-calculate the reciprocal differences.
                                     Set flag F02 if you prefer that the degree of the numerator is smaller than the degree of the denominator ( if n is even )
Subroutines: /
 
 

01  LBL "THIELE"   02  STO M   03  RCL 00   04  RCL 00   05  3   06  +   07  *   08  2   09  /   10  FS? 01   11  GTO 03   12  CF 10   13  FS? 02   14  SF 10   15  RCL 00   16   E3   17  /   18  STO O   19  .999   20  +

21  STO N   22  SIGN   23  RCL 00    24  +   25  STO P    26  STO Q   27  RCL 00             28  +   29  ISG Q   30  LBL 01   31  RCL N   32  INT   33  ISG X   34  CLX   35  RCL O   36  FRC   37  +   38  STO O   39  X<>Y   40  LBL 02

41  RCL IND N   42  RCL IND O   43  -   44  RCL IND P   45  FS? 10   46  1/X   47  RCL IND Q     48  FS? 10   49  1/X   50  -   51  /   52  STO IND Y   53  CLX   54  SIGN   55  ST+ Q   56  +   57  ISG O   58  GTO 02   59  CF 10   60  RCL Q

61  STO P    62  SIGN   63  ST+ Q   64  X<>Y   65  ISG N   66  GTO 01   67  LASTX   68  LBL 03   69  STO O   70  RCL 00   71  STO N   72  SIGN   73  ST- N   74  0   75  LBL 04   76  RCL IND O     77  +   78  RCL M   79  RCL IND N   80  -

81  X<>Y   82  /   83  X<>Y   84  1   85  +   86  ST- O   87  X<>Y   88  DSE N   89  GTO 04   90  RCL IND O     91  FS? 02   92  1/X   93  +   94  FS? 02   95  1/X   96  RCL M   97  SIGN   98  X<>Y   99  CLA 100  END

 
   ( 159 bytes SIZE (n²+3n+2)/2 )
 
 

STACK	INPUTS	OUTPUTS
X	x	f(x)
L	/	x

 
Example1:     Given  f(0) = 226 ,  f(2) = 58 ,  f(5) = 18 ,  f(10) = 6 ,  f(20) = 1     Calculate  f(1) ,  f(3)  ,  f(7)

-There are 5 points so   5  STO 00  and we store:

    0   226             R01   R06
    2    58              R02   R07
    5    18      in     R03   R08    respectively
   10    6               R04   R09
   20    1               R05   R10

-Then,  CF 01    1  XEQ "THIELE"   >>>>   f(1) =  109.4020   ( in 12 seconds )

-After executing "THIELE" once, the r_12...-values are in the required registers ( R11 R15 R18 R20 in this example ),
 and we can set flag F01 to speed up execution:

  SF 01   3   R/S  >>>>   f(3) =  35.8827   ( in 3 seconds )
              7   R/S  >>>>   f(7) =  10.8845

-Here, the number of points is odd ( n = 5 ), so setting or clearing flag F02 yields the same results ( provided no y-value equals zero )
 

Example2:     Given  f(0) = 226 ,  f(2) = 58 ,  f(10) = 6 ,  f(20) = 1      Evaluate again f(1) , f(3) , f(7)

-There are 4 points whence   4  STO 00  and

    0    STO 01    2   STO 02   10   STO 03   20  STO 04
  226  STO 05   58  STO 06    6    STO 07    1   STO 08

 1°)  With flag F02 clear CF 02, we find:

   CF 01     1  XEQ "THIELE"  >>>>   f(1) = 97.5178 ,  SF 01   3  R/S  >>>>   f(3) = 38.9928 ,  7  R/S  >>>>  f(7)  = 12.2710

 2°)  With flag F02 set SF 02, we find:

   CF 01     1  R/S  >>>>    f(1) = 100.1803 ,   SF 01   3  R/S  >>>>   f(3) = 37.8940 ,  7  R/S  >>>>  f(7)  = 11.5319
 

Notes:

-If  y₁ equals another y-value, there will be a DATA ERROR , so always choose the first point such that y₁ is different from all other y_k
-There is always a small risk that some denominators equal 0. If it ever happens, change the order of the data points and start again...

-This program only deals with the cases:

  deg(numerator) = deg(denominator)          if n is odd
  deg(numerator) = deg(denominator) +/-1  if n is even,  but there are many other possibilities.

-For instance, if you want to fit a set of data points ( x_i , y_i ) to a rational function of the type  1/p(x)  where p is a polynomial,
 you can use Lagrange's interpolation formula to the set ( x_i , 1/y_i ) and take the reciprocal of the results.

-In the first example above ( 5 points ) it gives  f(1) = 108.4804  ,   f(3) = 35.8694  ,  f(7) = 10.9406
 

References:

[1]  F. Scheid -  "Numerical Analysis" - Mc Graw Hill   ISBN 07-055197-9
[2]  Abramowitz and Stegun  -  "Handbook of Mathematical Functions"  -  Dover Publications  -  ISBN  0-486-61272-4