The Museum of HP Calculators


Rational Function Interpolation for the HP-41

This program is Copyright © 2006 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview


 

-The following program uses Thiele's formula:

-Given n data points  ( x1 , y1 ) , ( x2 , y2 ) , .................. , ( xn , yn ) ,   the rational function f(x) is evaluated by the continued fraction:

                       x - x1     x - x2      x - x3                                       x - xn-1
  f(x) =   y1 +  ------    -------    --------   ...............................  ---------
                       r12 +      r123 +     r1234 +                                      r12......n
 

  where  r1k = ( xk - x1 )/( yk - y1 )  ,    r12k = ( xk - x2 )/( r1k - r12 )  ,   r123k = ( xk - x3 )/( r12k - r123 )  , .............. ,  r12....n = ( xn - xn-1 )/( r12...n - r12...n-1 )
 

-Thus,  f(x) = p(x)/q(x)     with   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 = x1      Rn+1 = y1
                                        R02 = x2      Rn+2 = y2       R2n+1 = r12
                                         ---------        ----------         R2n+2 = r13     R3n = r123
                                         ---------        ----------         ------------       ---------

                                        Rnn = xn      R2n   = yn        R3n-1 = r1n     R4n-3 = r12n  ---------  R(n2+3n)/2 = r12.....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                ( synthetic )
  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                ( synthetic )
  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 (n2+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      into registers     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 r12...-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  y1 equals another y-value, there will be a DATA ERROR , so always choose the first point such that y1 is different from all other yk
-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  and  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 ( xi , yi ) 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 ( xi , 1/yi ) 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:        F. Scheid -  "Numerical Analysis" - Mc Graw Hill   ISBN 07-055197-9
                              Abramowitz and Stegun  -  "Handbook of Mathematical Functions"  -  Dover Publications  -  ISBN  0-486-61272-4
 

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall