hp41programs

Lagrange Lagrange's Interpolation Formula for the HP-41
 

Overview
 

 1°) One Dimension Problem

   a) Univariate Interpolation
   b) Derivative of the Lagrange Polynomial
   c) Integration of the Lagrange Polynomial

 2°) Bivariate Interpolation

   a) non-synthetic program
   b) synthetic program
 

1°) One Dimension Problem
 

    a) 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 eee+1 )
 
 

    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  yields   I2 = 6.912809880
  n = 4  yields   I4 = 6.521610110
  n = 8  yields   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:

- 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 )
-  n is usually doubled at every step, but it's not necessary here.
-  The same idea can be applied to any other problem of this kind.
 

    b) Derivative of the Lagrange Polynomial
 

-"DLAGR" takes x in X-register and returns  L'(x)  in X-register.
 

Data Registers:

-The same 2n numbers   x1 , y1 , x2 , y2 , ........ , xn , yn  are to be stored 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:  /
 
 

 01  LBL "DLAGR"
 02  CLA
 03  STO M           
 04  X<>Y
 05  STO Q
 06  STO a
 07  LBL 01
 08  RCL a
 09  STO O
 10  CLX
 11  LBL 02
 12  RCL O
 13  RCL Q
 14  X=Y?
 15  GTO 02
 16  X<> Z
 17  RCL a
 18  STO N           
 19  SIGN
 20  LBL 03
 21  RCL N
 22  RCL O
 23  X=Y?
 24  GTO 03
 25  CLX
 26  RCL Q
 27  X=Y?
 28  GTO 03
 29  RDN
 30  CLX
 31  RCL M           
 32  RCL IND N
 33  -
 34  *
 35  ENTER^
 36  ENTER^
 37  LBL 03
 38  R^
 39  R^
 40  ISG N
 41  ISG N
 42  GTO 03
 43  +
 44  STO Z
 45  LBL 02
 46  X<> Z
 47  ISG O
 48  ISG O
 49  GTO 02          
 50  RCL a
 51  X<>Y
 52  LBL 04
 53  RCL IND Q
 54  RCL IND Z
 55  -
 56  X=0?
 57  SIGN
 58  /
 59  ISG Y
 60  ISG Y
 61  GTO 04
 62  ISG Q
 63  RCL IND Q
 64  *
 65  ST+ P
 66  ISG Q
 67  GTO 01          
 68  RCL a
 69  RCL M
 70  SIGN
 71  X<> P
 72  CLA
 73  END

 
     ( 123 bytes / SIZE eee+1 )
 
 

    STACK      INPUTS    OUTPUTS
         Y       bbb.eee      bbb.eee
         X            x        L' (x)
         L            /          x

 
Example:

  Evaluate  y(PI)  for the collocation polynomial that takes the values prescribed below:
 
 

    xi     0     1     2     4     7
    yi     3     2     4     6     5

 
-If, 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^
    PI    XEQ "DLAGR"  >>>>  L' (PI) = 0.857353871                                  ---Execution time = 67s---

Notes:

-Synthetic register a is used, so this program cannot be called as more than a first level subroutine.
-Registers M N O P Q a  may of course be replaced by standard registers like R00 thru R05
-In this case, choose  bbb > 005.
 

    c) Integration of the Lagrange Polynomial
 

-This program is also listed in "Numerical Integration for the HP-41"

 n data points are given and stored in registers  R01 thru R2n
 
 

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

 
-The program below calculates its coefficients c0 , c1 , ........... , cn-1  when L(x) is written:

    L(x) = c0 + c1 ( x-x1 ) + c2 ( x-x1 )2 + .................. + cn-1 ( x-x1 )n-1

-First, x1 is subtracted from all other xi 's
-We have obviously  c0 = y1  and  the other yi 's  are replaced by  ( yi - y1 ) / ( xi - x1 )
-Then, we find the value of the new Lagrange polynomial at zero, thus producing c1 which is stored in register R04
-The process is repeated until  cn-1 is computed and stored in register R2n
-Finally, lines 42 to 57 evaluate  §x1xn  L(x).dx
 

Data Registers:       •  R00 = n = number of points                                ( Registers R00 thru R2n  are to be initialized before executing "LGI" )

                                  •  R01 = x1     •  R03 = x2      .......  •  R2n-1 = xn
                                  •  R02 = y1     •  R04 = y2      .......   •  R2n   = yn

  >>>>    when the program stops,  R02 = c0 = y1 ,  R04 = c1 , ............... ,  R2n = cn-1

Flags:  /
Subroutine:   "LAGR"  ( cf paragrph 1-a) )
 
 

01  LBL "LGI"
02  RCL 00
03  ST+ X
04  STO Y
05  RCL 01         
06  LBL 00
07  DSE Z
08  ST- IND Z 
09  DSE Z
10  GTO 00
11  CLX
12   E3
13  /
14  ISG X
15  STO 01
16  RCL 02
17  GTO 02         
18  LBL 01
19  2
20  ST+ 01
21  RCL 01
22  0
23  XEQ "LAGR"
24  LBL 02
25  ISG Y
26  STO IND Y
27  ISG Y
28  FS? 30
29  GTO 04         
30  LBL 03
31  RCL IND Y 
32  ISG Z
33  X<>Y
34  ST- IND Z
35  X<>Y
36  ST/ IND Z
37  RDN
38  ISG Y
39  GTO 03         
40  GTO 01
41  LBL 04
42  RCL 00
43  ST+ X
44  0
45  LBL 05
46  RCL IND Y 
47  RCL 00
48  /
49  +
50  RCL IND 01 
51  *
52  2
53  ST- Z
54  RDN
55  DSE 00
56  GTO 05         
57  END

 
   ( 98 bytes / SIZE 2n+1 )
 
 

      STACK        INPUT      OUTPUT
           X              /   §x1xn  L(x).dx

 
Example:     Calculate  §18  L(x).dx   from the following values
 
 

      x      1     2.4      4     5.2      7      8
      y      1      4      6      5      4      2

 
-Store these 12 numbers into registers

     R01  R03  R05  R07  R09  R11   ( x )
     R02  R04  R06  R08  R10  R12    ( y )    respectively

-There are 6 points:     6   STO 00      XEQ "LGI"     >>>>    §18  L(x).dx  =   29.61789480  ( in 43 seconds )

-We have also the coefficients of the Lagrange polynomial ( as a sum of powers of ( x-x1 ) )  in registers  R02  R04  .............  R12
-Here, x1 = 1 and we find:

  L(x) =  1 - 0.362103178 ( x-1 ) + 3.623795356 ( x-1 )2 - 1.661873944 ( x-1 )3 + 0.272598127 ( x-1 )4 - 0.015381483 ( x-1 )5
 

Remarks:

-For large n-values, the coefficients of the Lagrange polynomial are often great in magnitude
  and of alternate signs, thus leading to inaccurate results.
-Moreover, the coefficients themselves are not obtained very accurately if n is too great.

-For instance, with the following data:
 
 

  x   1   3   5   7   9  11  13  15  17  19  21  23  25  27  29  31  33  35  37  39
  y   2   4   6   8  10  12  14  16  18  20  22  24  26  28  30  32  34  36  38  40

 
 "LGI"  yields  797.9971774  wheras the exact result is 798

Execution time t:
 
 

       n        t
       4      15s
       6      43s
      10   2mn51s
      20  21mn13s

 

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          
79  SIGN 
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... )

            >>>   if, for instance, you choose bb = 11 ,  store these  19 numbers  into

     x \ y |   2    3    4    6                                                          |   R14   R18   R22   R26
--------------------------                                              ---------------------------------
       1   |   4    3    3    5                                                  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 table 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
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                =                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
 

References:

[1]   Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]   Press, Vetterling, Teukolsky, Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4