# 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, add  CLX  STO 03  after line 06 and of course, choose 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:

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