The Museum of HP Calculators

Lagrange's Interpolation Formula 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:

1°) Univariate Interpolation
2°) Bivariate Interpolation        ( new )
   a) non-synthetic program        ( new )
   b) synthetic program                ( new )

1°) Univariate Interpolation

-Let be a set of n data points: M₁ ( x₁ ; y₁ ) , M₂ ( x₂ ; y₂ ) , .......... , M_n ( x_n ; y_n )
There is a unique polynomial L of degree < n such that: L(x_i) = y_i i = 1 ; 2 ; ..... ; n
-In other words, L is the collocation polynomial for these arguments.
( x_i values can be unequally-spaced.)

-The following program calculates L(x) = y₁ L₁(x) + y₂ L₂(x) + ..... + y_n L_n(x)

where L_i(x) = ( x - x₁ ) ...... ( x - x_i-1 ) ( x - x_i+1 ) ...... ( x - x_n ) / [ ( x_i - x₁ ) ...... ( x_i - x_i-1 ) ( x_i - x_i+1 ) ( x_i - x_n ) ]
are the Lagrange's multiplier functions.

-This formula is also applied to an extrapolation to the limit in example2.

Data Registers:

-You have to store the 2n numbers   x₁ , y₁ , x₂ , y₂ , ........ , x_n , y_n into contiguous registers, starting at any register Rbb:

       x₁         x₂    ...........        x_n

       y₁         y₂    ...........        y_n

                                  into

      Rbb    Rbb+2 .............    Ree-1

    Rbb+1    Rbb+3 .............      Ree

( ee = bb + 2n -1 )

Flags: /
Subroutines: /

-Thanks to synthetic programming, "LAGR" only uses the registers containing the x_i and y_i numbers,
but you can replace status registers M N O P Q by the standard R00 R01 R02 R03 R04 ( in this case, of course, bb > 04 ).

01 LBL "LAGR"
02 CLA
03 STO M
04 X<>Y
05 STO O
06 STO Q
07 LBL 01
08 RCL Q
09 STO N
10 SIGN
11 LBL 02
12 RCL M
13 RCL IND O
14 RCL IND N
15 ST- Z
16 -
17 X#0?
18 GTO 03
19 SIGN
20 STO Y
21 LBL 03
22 /
23 *
24 ISG N
25 ISG N
26 GTO 02
27 ISG O
28 RCL IND O
29 *
30 ST+ P
31 ISG O
32 GTO 01
33 RCL Q
34 RCL M
35 SIGN
36 X<> P
37 CLA
38 END

( 69 bytes / SIZE? )

    STACK      INPUTS    OUTPUTS

         Y       bbb.eee      bbb.eee

         X            x        L(x)

         L            /          x

Example1:

Evaluate y(3) and y(5) for the collocation polynomial that takes the values prescribed below:

x_i 0 1 2 4 7

y_i 3 2 4 6 5

For instance, you store these 10 numbers ( 0 3 1 2 2 4 4 6 7 5 ) into R01 thru R10, then

1.010 ENTER^
3 XEQ "LAGR" produces: y(3) = 5.847619048 ( in 14 seconds )

RDN 5 R/S ---------- y(5) = 4.523809520

Example2:

-Lagrange's polynomial can also be used for an extrapolation to the limit:
Suppose we are using the trapezoidal rule to evaluate the integral: I = §₀² e^x dx
The interval [0;2] is divided into n subintervals and

n = 1 yields I₁ = 8.389056101
n = 2 -----   I₂ = 6.912809880
n = 4 -----   I₄ = 6.521610110
n = 8 -----   I₈ = 6.422297820

We would like to know the result with n = infinity! i-e 1/n = 0.
The error of the trapezoidal rule is nearly proportional to 1/n² , therefore, we can use Lagrange interpolation formula with

x₁ = 1 x₂ = 1/4 x₃ = 1/16 x₄ = 1/64
y₁ = I₁ y₂ = I₂ y₃ = I₄ y₄ = I₈ and evaluate L (0)

If these 8 numbers are stored in R11 thru R18 for instance,

11.018 ENTER^
0 XEQ "LAGR" produces 6.389056387

Exact result is of course e² - 1 = 6.389056099 to ten places
and the error of our Lagrange approximation is only 3. 10^-7 although I₈ was only exact to 1 decimal!

N.B:   1- If you are using Simpson's rule instead of the trapezoidal rule, store x_n = 1/n⁴ : it's a 4th order method ( at least far from any singularity )
      2- n is usually doubled at every step, but it's not necessary here.
    3- The same idea can be applied to any other problem of this kind.

2°) Bivariate Interpolation

a) Non-synthetic program

-Now we have a grid of nxm values: f(x_i,y_j) with i = 1,2,......,n & j = 1,2,.....,m and we want to estimate f(x,y)

-"LAGR2" uses the Lagrange Polynomial in the directions of x-axis and y-axis.
-You have to store the data in registers Rbb to Ree as shown below ( bb > 10 ):

Data Registers: • R00 = bbb.eeenn ( Registers R00 and Rbb thru Ree are to be initialized before executing "LAGR2" )

R01 = x R03 = f(x,y)
R02 = y R04 thru R10: temp We must choose bb > 10

                                                 x \ y |        • Rbb+n = y₁                • Rbb+2n+1 = y₂       ..............................    • Ree-n = y_m
                                  ---------------------------------------------------------------------------------------------------------------------
                                • Rbb         = x₁ |   • Rbb+n+1 = f(x₁,y₁)    • Rbb+2n+2 = f(x₁,y₂)    ..............................     • Ree-n+1 = f(x₁,y_m)
                                • Rbb+1     = x₂ |   • Rbb+n+2 = f(x₂,y₁)    • Rbb+2n+3 = f(x₂,y₂)    ..............................     • Ree-n+2 = f(x₂,y_m)
                                   ....................    |     .............................................................................................................................................
                                                           |
                                • Rbb+n-1 = x_n |   • Rbb+2n = f(x_n,y₁)     • Rbb+3n+1 = f(x_n,y₂)   ...............................     • Ree = f(x_n,y_m)

Flags: /
Subroutines: /

01 LBL "LAGR2"
02 STO 01
03 X<>Y
04 STO 02
05 CLX
06 STO 03
07 RCL 00
08 ISG X
09 STO 07
10 INT
11 DSE X
12 .1
13 %
14 RCL 00
15 INT
16 +
17 STO 05
18 STO 06
19   E-5
20 ST+ 07
21 LBL 00
22 RCL 07
23 STO 08
24 STO 09
25 FRC
26 RCL 06
27 INT
28 +
29 ISG X
30 STO 10
31 CLX
32 STO 04
33 LBL 01
34 RCL IND 10
35 LBL 02
36 RCL 02
37 RCL IND 08
38 RCL IND 09
39 ST- Z
40 -
41 X#0?
42 GTO 02
43 SIGN
44 STO Y
45 LBL 02
46 /
47 *
48 ISG 09
49 GTO 02
50 ST+ 04
51 RCL 07
52 STO 09
53 ISG 10
54 CLX
55 ISG 08
56 GTO 01
57 RCL 05
58 STO 10
59 RCL 04
60 LBL 03
61 RCL 01
62 RCL IND 06
63 RCL IND 10
64 ST- Z
65 -
66 X#0?
67 GTO 03
68 SIGN
69 STO Y
70 LBL 03
71 /
72 *
73 ISG 10
74 GTO 03
75 ST+ 03
76 ISG 06
77 GTO 00
78 RCL 01                            Lines 78-79-80 are only useful to save y in Y-register and x in L-register
79 SIGN                               Otherwise, these lines may be deleted.
80 RCL 02
81 RCL 03
82 END

( 121 bytes / SIZE? )

    STACK      INPUTS    OUTPUTS

         Y            y          y

         X            x       f(x,y)

         L            /          x

Example: A function f(x,y) is only known by the following values ( meaning for example: f(1,4) = 3 f(4,6) = 9 ...etc... )

     x \ y |   2    3    4    6                                                                                                      |   R14   R18   R22   R26
--------------------------              if, for instance, you choose bb = 11                      ---------------------------------
       1   |   4    3    3    5                 store these 19 numbers into                                  R11 | R15   R19   R23   R27            respectively
       2   |   3    1    2    6                                                                                              R12 | R16   R20   R24   R28
       4   |   1    0    4    9                                                                                              R13 | R17   R21   R25   R29

-The control number of this tableau is 11.02903 so 11.02903 STO 00
( note that nn = 03 is the number of x-values which may be different from m = the number of y-values ( here m = 4 ) )

-Then, to compute, say f(3,5)

5 ENTER^
3 XEQ "LAGR2" >>>> f(3,5) ~ 5.8333 ( in 29 seconds )

-Likewise, to obtain f(1.6,2.7)

2.7 ENTER^
1.6 R/S >>>> f(1.6,2.7) ~ 1.9038

Note: Execution time is approximately proportional to nxm

b) Synthetic program

-Synthetic programming allows to store the data into registers R01 thru Rn.m+n+m
-However, synthetic register a is used, so this program cannot be called as more than a first level subroutine.

Data Registers: • R00 = n.mmm = n + m/10³ ( Registers R00 thru Rnm+m+n are to be initialized before executing "BVI" )

                                          x \ y |       • Rnn+1 = y₁             • R2n+2 = y₂       ..............................    • Rnm+m = y_m
                                  ---------------------------------------------------------------------------------------------------------------------
                                • R01 = x₁ |   • Rnn+2 = f(x₁,y₁)    • R2n+3 = f(x₁,y₂)    ..............................     • Rnm+m+1 = f(x₁,y_m)
                                • R02 = x₂ |   • Rnn+3 = f(x₂,y₁)    • R2n+4 = f(x₂,y₂)    ..............................     • Rnm+m+2 = f(x₂,y_m)
                                    ...........    |      .............................................................................................................................................
                                                   |
                                • Rnn = x_n |   • R2n+1 = f(x_n,y₁)   • R3n+2 = f(x_n,y₂)   ..............................      • Rnm+m+n = f(x_n,y_m)

Flags: /
Subroutines: /

01 LBL "BVI"
02 CLA
03 STO M
04 X<>Y
05 STO N
06 RCL 00
07 INT
08 STO P ( synthetic )
09 LBL 01
10 RCL 00
11 INT
12 LASTX
13 FRC
14 ENTER^
15 SIGN
16 10^X
17 CHS
18 10^X
19 SQRT
20 +
21 ST* Y
22 +
23 ISG X
24 STO Q
25 CLX
26 LBL 02
27 RCL Q
28 FRC
29 ISG X
30 STO a
31 SIGN
32 LBL 03
33 RCL N
34 RCL IND a
35 -
36 RCL IND Q
37 ST- L
38 X<> L
39 CHS
40 X#0?
41 GTO 03
42 SIGN
43 STO Y
44 LBL 03
45 /
46 *
47 ISG a
48 GTO 03
49 RCL P
50 RCL Q
51 +
52 RDN
53 RCL IND T
54 *
55 +
56 ISG Q
57 GTO 02
58 RCL 00
59 INT
60 STO Q
61 X<>Y
62 LBL 04
63 RCL M
64 RCL IND P
65 RCL IND Q
66 ST- Z
67 -
68 X#0?
69 GTO 04
70 SIGN
71 STO Y
72 LBL 04
73 /
74 *
75 DSE Q
76 GTO 04
77 ST+ O
78 DSE P
79 GTO 01
80 RCL M
81 SIGN
82 RCL N
83 RCL O
84 CLA
85 END

( 131 bytes / SIZE (n+1)(m+1) )

    STACK      INPUTS    OUTPUTS

         Y            y          y

         X            x       f(x,y)

         L            /          x

Example: With the same data as above:

     x \ y |   2    3    4    6                                                                                  |   R04   R08   R12   R16
--------------------------                                                                        --------------------------------
       1   |   4    3    3    5                 are to be stored into                           R01 | R05   R09   R13   R17            respectively
       2   |   3    1    2    6                                                                           R02 | R06   R10   R14   R18
       4   |   1    0    4    9                                                                           R03 | R07   R11   R15   R19

-Here, n = 3 and m = 4 so 3.004 STO 00

-Then, to compute f(3,5)

5 ENTER^
3 XEQ "BVI" >>>> f(3,5) ~ 5.8333 ( in 31 seconds )

-And to obtain f(1.6,2.7)

2.7 ENTER^
1.6 R/S >>>> f(1.6,2.7) ~ 1.9038

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