The Museum of HP Calculators

# Lagrange's Interpolation Formula for the HP-41

This program is Copyright © 2003-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:

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:   M1 ( x1 ; y1 )  ,   M2 ( x2 ; y2 )  ,  ..........  ,  Mn ( xn ; yn )
There is a unique polynomial L of degree < n  such that:     L(xi) = yi       i = 1 ; 2 ; ..... ; n
-In other words, L is the collocation polynomial for these arguments.
( xi values can be unequally-spaced.)

-The following program calculates  L(x) = y1 L1(x) + y2 L2(x) + ..... + yn Ln(x)

where    Li(x) = ( x - x1 ) ...... ( x - xi-1 ) ( x - xi+1 ) ...... ( x - xn ) / [ ( xi - x1 ) ...... ( xi - xi-1 ) ( xi - xi+1 ) ( xi - xn ) ]
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   x1 , y1 , x2 , y2 , ........ , xn , yn  into contiguous registers,  starting at any register  Rbb:

 x1 x2 ........... xn y1 y2 ........... yn

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  xi and yi 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:

 xi 0 1 2 4 7 yi 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  =  §02  ex dx
The interval  [0;2]  is divided into n subintervals and

n = 1  yields  I1 = 8.389056101
n = 2  -----   I2 = 6.912809880
n = 4  -----   I4 = 6.521610110
n = 8  -----   I8 = 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/n2  , therefore, we can use Lagrange interpolation formula with

x1 = 1   x2 = 1/4   x3 = 1/16   x4 = 1/64
y1 = I1   y2 = I2      y3 = I4       y4 =  I8      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 e2 - 1 = 6.389056099 to ten places
and the error of our Lagrange approximation is only  3. 10-7  although I8 was only exact to 1 decimal!

N.B:   1- If you are using Simpson's rule instead of the trapezoidal rule, store  xn = 1/n4  : 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(xi,yj)   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 = y1                •  Rbb+2n+1 = y2       ..............................    •  Ree-n = ym
---------------------------------------------------------------------------------------------------------------------
•  Rbb         = x1 |   •  Rbb+n+1 = f(x1,y1)    •  Rbb+2n+2 = f(x1,y2)    ..............................     •  Ree-n+1 = f(x1,ym)
•  Rbb+1     = x2 |   •  Rbb+n+2 = f(x2,y1)    •  Rbb+2n+3 = f(x2,y2)    ..............................     •  Ree-n+2 = f(x2,ym)
....................    |     .............................................................................................................................................
|
•  Rbb+n-1 = xn  |   •  Rbb+2n  = f(xn,y1)     •  Rbb+3n+1 = f(xn,y2)   ...............................     •  Ree  = f(xn,ym)

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/103            ( Registers R00 thru Rnm+m+n are to be initialized before executing "BVI"  )

x  \  y |       •  Rnn+1 = y1             •  R2n+2 = y2       ..............................    •  Rnm+m = ym
---------------------------------------------------------------------------------------------------------------------
•  R01  = x1 |   •  Rnn+2 = f(x1,y1)    •  R2n+3 = f(x1,y2)    ..............................     •  Rnm+m+1 = f(x1,ym)
•  R02  = x2 |   •  Rnn+3 = f(x2,y1)    •  R2n+4 = f(x2,y2)    ..............................     •  Rnm+m+2 = f(x2,ym)
...........    |      .............................................................................................................................................
|
•  Rnn = xn  |   •  R2n+1  = f(xn,y1)   •  R3n+2 = f(xn,y2)   ..............................      •  Rnm+m+n  = f(xn,ym)

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