RatFunInter

# Rational Function Interpolation for the HP-41

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)   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 = 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    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 (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      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 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
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