The Museum of HP Calculators

Rational Function Interpolation for the HP-41

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