hp41programs

Double Integrals for the HP-41
 

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.
 
 

01  LBL "DIN"  02  1.6 03  STO 07 04  .15 05  SQRT 06  STO 08  07  RCL 02 08  RCL 01 09  STO 10 10  - 11  RCL 03         12  STO 15 13  / 14  STO 09 15  2 16  / 17  ST+ 10 18  CLX 19  STO 06 20  LBL 01

21  RCL 10 22  RCL 09 23  RCL 08 24  * 25  - 26  XEQ 02  27  ST+ 06 28  RCL 10  29  XEQ 02 30  RCL 07         31  * 32  ST+ 06 33  RCL 10 34  RCL 09 35  ST+ 10 36  RCL 08 37  * 38  + 39  XEQ 02 40  ST+ 06

41  DSE 15 42  GTO 01 43  GTO 04 44  LBL 02 45  STO 13 46  XEQ IND 04 47  STO 12  48  RCL 13 49  XEQ IND 05 50  RCL 12 51  - 52  RCL 03 53  STO 16 54  / 55  STO 11 56  2 57  / 58  ST+ 12 59  CLX 60  STO 14

61  LBL 03 62  RCL 12 63  RCL 11 64  RCL 08 65  * 66  - 67  RCL 13 68  XEQ IND 00 69  ST+ 14 70  RCL 12 71  RCL 13 72  XEQ IND 00 73  RCL 07 74  * 75  ST+ 14 76  RCL 12 77  RCL 11 78  ST+ 12 79  RCL 08 80  *

81  + 82  RCL 13 83  XEQ IND 00 84  ST+ 14 85  DSE 16 86  GTO 03 87  RCL 14 88  RCL 11 89  * 90  RTN 91  LBL 04 92  RCL 06 93  RCL 09 94  * 95  3.6 96  X^2 97  / 98  STO 06 99  END

 
   ( 140 bytes / SIZE 017 )
 
 

STACK	INPUTS	OUTPUTS
X	/	I

 
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
 
 
 

01 LBL "DIN"  02 1.6  03 STO 07  04 .15  05 SQRT  06 STO 08   07 RCL 02  08 RCL 01  09 STO 10  10 -  11 RCL 03          12 STO 16  13 /  14 STO 09  15 2  16 STO 14  17 /  18 ST+ 10  19 CLX

20 STO 05  21 GTO 01  22 LBL 02  23 STO 13  24 XEQ IND 04  25 STO 12   26 -  27 RCL 03  28 STO 17  29 /  30 STO 11  31 RCL 14  32 /  33 ST+ 12  34 CLX  35 STO 15  36 LBL 03  37 RCL 12  38 RCL 11

39 RCL 08  40 *  41 -  42 RCL 13  43 XEQ IND 00  44 ST+ 15  45 RCL 12  46 RCL 13  47 XEQ IND 00  48 RCL 07  49 *  50 ST+ 15  51 RCL 12  52 RCL 11  53 ST+ 12  54 RCL 08  55 *  56 +  57 RCL 13

58 XEQ IND 00  59 ST+ 15  60 DSE 17  61 GTO 03  62 RCL 11  63 RCL 15  64 *  65 RTN  66 LBL 01  67 RCL 10  68 RCL 09  69 RCL 08  70 *  71 -  72 XEQ 02  73 ST+ 05  74 RCL 10  75 XEQ 02  76 RCL 07

77 *  78 ST+ 05  79 RCL 10  80 RCL 09  81 ST+ 10  82 RCL 08  83 *  84 +  85 XEQ 02  86 ST+ 05  87 DSE 16  88 GTO 01  89 RCL 05          90 RCL 09  91 *  92 3.6  93 X^2  94 /  95 STO 05  96 END

 
    ( 137 bytes / SIZE 018 )
 
 

STACK	INPUT	OUTPUT
X	/	Integral

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

01 LBL "S"  02 STO Z  03 ST* Z  04 X<>Y  05 X^2  06 ST+ Z  07 X<> L  08 *  09 LN1+X  10 X<>Y  11 SQRT  12 /  13 RTN  14 LBL "U"  15 X^2  16 LASTX  17 RTN  18 END

 
    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
 
 

01 LBL "DIN4"  02 3  03 STO Y  04 4.8  05 SQRT  06 ST+ Z  07 -  08 28  09 ST/ Z  10 /  11 SQRT  12 STO 07          13 X<>Y  14 SQRT  15 STO 08  16 .5  17 STO Y  18 43.2  19 SQRT  20 1/X  21 ST+ Z  22 -  23 STO 19  24 /

25 STO 15  26 RCL 02  27 RCL 01  28 STO 10  29 -  30 RCL 03  31 STO 16  32 /  33 STO 09  34 2  35 STO 14  36 /  37 ST+ 10  38 CLX  39 STO 05  40 GTO 01  41 LBL 02  42 STO 13  43 XEQ IND 04  44 STO 12   45 -  46 RCL 03  47 STO 17  48 /

49 STO 11  50 RCL 14  51 /  52 ST+ 12  53 CLX  54 STO 18  55 LBL 03  56 RCL 12  57 RCL 11  58 RCL 08  59 *  60 -  61 RCL 13  62 XEQ IND 00  63 ST+ 18  64 RCL 12          65 RCL 11  66 RCL 07  67 *  68 -  69 RCL 13  70 XEQ IND 00  71 RCL 15  72 *

73 ST+ 18  74 RCL 12  75 RCL 11  76 RCL 07  77 *  78 +  79 RCL 13  80 XEQ IND 00  81 RCL 15  82 *  83 ST+ 18  84 RCL 12  85 RCL 11  86 ST+ 12  87 RCL 08  88 *  89 +  90 RCL 13  91 XEQ IND 00  92 ST+ 18  93 DSE 17  94 GTO 03  95 RCL 11  96 RCL 18

97 *  98 RTN  99 LBL 01 100 RCL 10 101 RCL 09 102 RCL 08 103 * 104 - 105 XEQ 02 106 ST+ 05 107 RCL 10 108 RCL 09        109 RCL 07 110 * 111 - 112 XEQ 02 113 RCL 15  114 * 115 ST+ 05 116 RCL 10 117 RCL 09 118 RCL 07 119 * 120 +

121 XEQ 02 122 RCL 15 123 * 124 ST+ 05 125 RCL 10 126 RCL 09 127 ST+ 10 128 RCL 08 129 * 130 + 131 XEQ 02 132 ST+ 05 133 DSE 16 134 GTO 01 135 RCL 05        136 RCL 09  137 * 138 4 139 / 140 RCL 19 141 X^2 142 * 143 STO 05 144 END

 
     ( 201 bytes / SIZE 020 )
 
 

STACK	INPUT	OUTPUT
X	/	Integral

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

01 LBL "S"  02 STO Z  03 ST* Z  04 X<>Y  05 X^2  06 ST+ Z  07 X<> L  08 *  09 LN1+X  10 X<>Y  11 SQRT  12 /  13 RTN  14 LBL "U"  15 X^2  16 LASTX  17 RTN  18 END

 
    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++" ...............