hp41programs

Overview
 

 1°) Gauss-Legendre 3-point Formula

  a) Program#1
  b) Program#2

 2°) Gauss-Legendre 4-point Formula
 
 

1°)  Gauss-Legendre 3-point Formula
 

     a) Program#1
 

We want to evaluate:   §_a^b §_u(x)^v(x)  f(x;y) dx dy                                                                  § =Integral

The idea behind the following program is to use the 3 point Gauss-Legendre formula
in both x-axis and y-axis. There are at least 3 reasons for this choice:
 1-The formula is relatively simple.
 2-It has a quite good accuracy ( better than Simpson's rule )
 3-It can be used even if the integrand has a singularity at the endpoints of the interval,
although in this case, convergence is of course much slower.

However, when the endpoints are infinite, it's necessary to make a change of argument ( like x = tan u )
to reduce the interval of integration to a bounded one.

Data Registers:        •  R00 = f name                           ( Registers R00 thru R05 are to be initialized before executing "DIN" )

         •  R01 = a                                   R06 = I               R09 to R16: temp
         •  R02 = b
         •  R03 = n = the number of subintervals into which the intervals of integration are to be divided.
         •  R04 = the name of the subroutine which calculates u(x).
         •  R05 = the name of the subroutine which calculates v(x).

Flags:  /
Subroutines:  3 subroutines must be keyed into program memory:

      -The first one must take x, y and z from the X-register,the Y-register and the Z-register respectively and calculate f(x;y;z).
      -The second takes x from the X-register and calculates u(x).
      -The third takes x from the X-register and calculates v(x).

N.B:  When the program stops, the result is in X-register and in R06.
 
 

 
   ( 140 bytes / SIZE 017 )
 
 

 
Example:   Evaluate  I = §₁² §_x^{x^2} sqrt(1+x⁴ y⁴) dx dy.
 

We must load the 3 needed subroutines into program memory ( any global labels, maximum of 6 characters ),
for instance:

 01  LBL "FF"
 02  *
 03  X^2
 04  X^2
 05  SIGN
 06  LAST X
 07  +
 08  SQRT
 09  RTN
 10  LBL "U"
 11  RTN
 12  LBL "V"
 13  X^2
 14  RTN

and then,

  "FF" ASTO 00
    1     STO 01
    2     STO 02
  "U"  ASTO 04
  "V"  ASTO 05

    n = 1     1  STO 03  XEQ "DIN"  >>>   15.45937082    in the X-register and in R06
    n = 2     2  STO 03      R/S           >>>   15.46673275
    n = 4     4  STO 03      R/S           >>>  15.46686031
    n = 8     8  STO 03      R/S           >>>   15.46686245  ... all the digits are correct!

With n = 8, execution time is about 7mn16s.
(When n is multiplied by 2, execution time is roughly multiplied by 4)
 

     b) Program#2
 

-In fact, it's preferable to use 1 subroutine instead of 2 to calculate u(x) & v(x)
 
 

Data Registers:           •  R00 = f name                            ( Registers R00 thru R04 are to be initialized before executing "DIN" )

                                      •  R01 = a                                    •  R04 =  uv name       R05 = I      R06 to R17: temp
                                      •  R02 = b
                                      •  R03 = n = Nb of subintervals

Flags: /
Subroutines:    A subroutine that takes x, y and z from the X-register,the Y-register and the Z-register respectively and returns f(x;y;z) in X-register.
                          A subroutine that takes x in X-register and returns v(x) in Y-register & u(x) in X-register
 
 
 

 
    ( 137 bytes / SIZE 018 )
 
 

 
Example:      Evaluate   I = §₁² §_x^{x^2}   [ Ln(1+x.y)] / sqrt ( x² + y² ) dx dy
 
 

 
    S  ASTO 00

    U  ASTO 04

    1   STO 01
    2   STO 02

-With n = 1    1  STO 03   XEQ "DIN"   >>>>   I = 0.456387227                            ---Execution time = 14s---
-With n = 2    2  STO 03        R/S          >>>>   I = 0.456373589
-With n = 4    4  STO 03        R/S          >>>>   I = 0.456373361
 

2°)  Gauss-Legendre 4-point Formula
 

-A better accuracy may of course be obtained by Gauss-Legendre 4-point formula:
 

     §_-1¹ f(x).dx ~ A.f(-a) + B.f(-b) + B.f(b) + A.f(a)

  with     a = sqrt[(3-sqrt(4.8))/7]      A = 0.5+1/sqrt(43.2)
             b = sqrt[(3+sqrt(4.8))/7]     B = 0.5-1/sqrt(43.2)
 

Data Registers:           •  R00 = f name                            ( Registers R00 thru R04 are to be initialized before executing "DIN4" )

                                      •  R01 = a                                    •  R04 =  uv name       R05 = I      R07 to R19: temp
                                      •  R02 = b
                                      •  R03 = n = Nb of subintervals

Flags: /
Subroutines:    A subroutine that takes x, y and z from the X-register,the Y-register and the Z-register respectively and returns f(x;y;z) in X-register.
                          A subroutine that takes x in X-register and returns v(x) in Y-register & u(x) in X-register
 
 

 
     ( 201 bytes / SIZE 020 )
 
 

 
Example:      Evaluate   I = §₁² §_x^{x^2}   [ Ln(1+x.y)] / sqrt ( x² + y² ) dx dy
 
 

 
    S  ASTO 00

    U  ASTO 04

    1   STO 01
    2   STO 02

-With n = 1    1  STO 03   XEQ "DIN4"   >>>>   I = 0.456373416                          ---Execution time = 25s---
-With n = 2    2  STO 03        R/S             >>>>   I = 0.456373357                          ---Execution time = 84s---
-With n = 4    4  STO 03        R/S             >>>>   I = 0.456373358                         ---Execution time = 313s---
 
 

References:

[1]   Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]   F. Scheid -  "Numerical Analysis" - Mc Graw Hill   ISBN 07-055197-9
[3]   Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes" in Fortran or in C++" ...............