Dblinteg

# Triple 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

-"TIN" evaluates:  I = §ab §u(x)v(x)§w(x;y)t(x;y)   f (x;y;z) dx dy dz                       ( § = Integral )

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

•  R01 = a                                   R08 = I               R09 to R23: 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).
•  R06 = the name of the subroutine which calculates w(x;y)
•  R07 = the name of the subroutine which calculates t(x;y)

Flags:  /
Subroutines:  5 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).
-Another one takes x and y from the X and Y registers respectively and calculates w(x;y)
-The last one takes x and y from the X and Y registers respectively and calculates t(x;y)

N.B:

-When the program stops, the result is in X-register and in R08.
-Line 43 is a synthetic three-byte GTO ( but it can be replaced by a two-byte GTO  because this line is executed only once ).

 01  LBL "TIN"    02  1.6   03  STO 09   04  .15   05  SQRT   06  STO 10   07  RCL 02   08  RCL 01   09  STO 12   10  -   11  RCL 03   12  STO 21   13  /   14  STO 11           15  2   16  /   17  ST+ 12   18  CLX   19  STO 08    20  LBL 01   21  RCL 12   22  RCL 11   23  RCL 10   24  *   25  -   26  XEQ 02   27  ST+ 08   28  RCL 12   29  XEQ 02   30  RCL 09   31  * 32  ST+ 08   33  RCL 12   34  RCL 11   35  ST+ 12   36  RCL 10   37  *   38  +   39  XEQ 02   40  ST+ 08   41  DSE 21   42  GTO 01   43  GTO 06   44  LBL 02   45  STO 15   46  XEQ IND 04   47  STO 14   48  RCL 15   49  XEQ IND 05   50  RCL 14   51  -   52  RCL 03   53  STO 22   54  /   55  STO 13   56  2   57  /   58  ST+ 14   59  CLX   60  STO 19   61  LBL 03   62  RCL 14 63  RCL 13   64  RCL 10   65  *   66  -   67  XEQ 04   68  ST+ 19   69  RCL 14   70  XEQ 04   71  RCL 09   72  *   73  ST+ 19   74  RCL 14   75  RCL 13   76  ST+ 14   77  RCL 10   78  *   79  +   80  XEQ 04   81  ST+ 19   82  DSE 22   83  GTO 03   84  RCL 19   85  RCL 13   86  *   87  RTN   88  LBL 04   89  STO 17   90  RCL 15   91  XEQ IND 06   92  STO 18   93  RCL 17 94  RCL 15   95  XEQ IND 07   96  RCL 18   97  -   98  RCL 03   99  STO 23 100  / 101  STO 16 102  2 103  / 104  ST+ 18 105  CLX 106  STO 20 107  LBL 05 108  RCL 18 109  RCL 16 110  RCL 10 111  * 112  - 113  RCL 17 114  RCL 15 115  XEQ IND 00 116  ST+ 20 117  RCL 18 118  RCL 17 119  RCL 15 120  XEQ IND 00 121  RCL 09 122  * 123  ST+ 20 124  RCL 18 125  RCL 16 126  ST+ 18 127  RCL 10 128  * 129  + 130  RCL 17 131  RCL 15 132  XEQ IND 00 133  ST+ 20 134  DSE 23 135  GTO 05 136  RCL 20 137  RCL 16 138  * 139  RTN 140  LBL 06 141  RCL 08 142  RCL 11 143  * 144  RCL 09 145  2 146  + 147  3 148  Y^X 149  / 150  STO 08 151  END

( 226 bytes / SIZE 024 )

 STACK INPUT OUTPUT X / Integral

Example:      Evaluate   I = §12 §xx^2 §x+yxy  xyz / sqrt ( x2 + y2 + z2 ) dx dy dz

The 5 required subroutines are, for instance:
( with global labels of 6 characters maximum )

 01  LBL "FF"  02  ENTER^  03  X^2  04  R^  05  ST* Z  06  X^2  07  +  08  R^  09  ST* Z  10  X^2  11  +  12  SQRT  13  /  14  RTN  15  LBL "U"  16  RTN  17  LBL "V"  18  X^2  19  RTN  20  LBL "W"  21  +  22  RTN  23  LBL "T"  24  *  25  RTN

( All the calculations fit into the stack, but if the functions are very complicated,
one may of course use data registers R24 ...etc... )

Then we initialize registers R00 thru R07:

"FF"   ASTO 00
1        STO 01
2        STO 02
"U"    ASTO 04
"V"    ASTO 05
"W"    ASTO 06
"T"    ASTO 07

n = 1      1  STO 03    XEQ "TIN"  >>>>    I1 =  0.765014888
n = 2      2  STO 03         R/S         >>>>    I2 =  0.770640690    ( in 4mn 29s )
n = 4      4  STO 03         R/S         >>>>    I4 =  0.770731245
n = 8      8  STO 03         R/S         >>>>    I8 =  0.770732669    ( in 4h02mn )                   ( the exact value is 0.7707326901 to ten places )

-Unfortunately, with triple integrals, when n is multiplied by 2, execution time is roughly multiplied by 8 !

-An HP-48GX ( in 10 FIX format ) gives 0.7707326900   in 26 hours with the built-in  ò function !

In 7 FIX format, the result   0.7707326903   is obtained in 3h49mn.
In 6 FIX format, the result   0.7707326923   is obtained in 1h38mn.

-Thus,  "TIN" is not so good as the ò function of the HP-48, but the HP-41 still bears comparison ...

-Execution time can be reduce by keying the 5 subroutines inside the "TIN" program itself ( with local labels ).
The program will run about 17% faster.

-Extrapolation to the limit may also be used to improve accuracy:

There are some theoretical reasons to think that the error is roughly inversely proportional to n6:     I = In + k/n6
Applying this formula with n = 4 and n = 8 leads to a system that yields:    I = 0.7707326916     the error is only 1.5 10 -9 .

b) Program#2

-It's in fact preferable to use just 2 subroutines instead of 4 to calculate u(x) & v(x)  and  w(x,y) & t(x,y)
-We can also use a flag to evaluate also double integrals with the same program:

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

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

Flags: F02                        CF 02 = Triple Integral
SF 02 = Double Integral

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
A subroutine that takes x in X-register & y in Y-register and returns t(x,y) in Y-register & w(x,y) in X-register.

-Line 21 is a three-byte GTO 01

 01 LBL "DTI"  02 .15  03 SQRT  04 STO 09  05 1.6  06 STO 11  07 RCL 02  08 RCL 01  09 STO 21  10 -  11 RCL 03  12 STO 17  13 /  14 STO 10  15 2  16 STO 06  17 /  18 ST+ 21  19 CLX  20 STO 22  21 GTO 01  22 LBL 02  23 STO 14  24 XEQ IND 04  25 STO 13  26 -  27 RCL 03  28 STO 18  29 / 30 STO 12  31 RCL 06  32 /  33 ST+ 13  34 CLX  35 STO 16  36 LBL 03  37 RCL 13  38 RCL 12  39 RCL 09  40 *  41 -  42 XEQ 04  43 ST+ 16  44 RCL 13  45 XEQ 04  46 RCL 11          47 *  48 ST+ 16  49 RCL 13  50 RCL 12   51 ST+ 13  52 RCL 09  53 *  54 +  55 XEQ 04  56 ST+ 16  57 DSE 18  58 GTO 03 59 RCL 12  60 RCL 16  61 *  62 RTN  63 LBL 04  64 STO 15  65 RCL 14  66 FS? 02  67 GTO IND 00  68 XEQ IND 05  69 STO 07  70 -  71 RCL 03  72 STO 20  73 /  74 STO 08  75 RCL 06  76 /  77 ST+ 07  78 CLX  79 STO 19  80 LBL 05  81 RCL 07  82 RCL 08  83 RCL 09  84 *  85 -  86 RCL 15  87 RCL 14 88 XEQ IND 00  89 ST+ 19  90 RCL 07  91 RCL 15  92 RCL 14  93 XEQ IND 00  94 RCL 11  95 *  96 ST+ 19  97 RCL 07  98 RCL 08  99 ST+ 07 100 RCL 09 101 * 102 + 103 RCL 15 104 RCL 14 105 XEQ IND 00 106 ST+ 19 107 DSE 20 108 GTO 05 109 RCL 08 110 RCL 19 111 * 112 RTN 113 LBL 01 114 RCL 21 115 RCL 10 116 RCL 09 117 * 118 - 119 XEQ 02 120 ST+ 22 121 RCL 21 122 XEQ 02 123 RCL 11 124 * 125 ST+ 22 126 RCL 21 127 RCL 10 128 ST+ 21 129 RCL 09 130 * 131 + 132 XEQ 02 133 ST+ 22 134 DSE 17 135 GTO 01 136 RCL 10         137 RCL 22 138 * 139 3.6 140 FC? 02 141 ST/ Y 142 X^2 143 / 144 STO 06 145 END

( 216 bytes / SIZE 023 )

 STACK INPUT OUTPUT X / Integral

Example1:      Evaluate   I = §12 §xx^2 §x+yx^2.y  [ Ln(1+x.y.z)] / sqrt ( x2 + y2 + z2 ) dx dy dz

 01 LBL "S"  02 STO T  03 ST* T  04 X<>Y  05 X^2  06 ST+ T  07 X<> L  08 *  09 X<>Y  10 X^2  11 ST+ Z  12 X<> L  13 *  14 LN1+X  15 X<>Y  16 SQRT  17 /  18 RTN  19 LBL "U"  20 X^2  21 LASTX  22 RTN  23 LBL "W"  24 STO Z  25 ST* Z  26 X<>Y  27 ST* Z  28 +  29 RTN  30 END

CF 02   ( triple integral )

S  ASTO 00

U  ASTO 04
V  ASTO 05

1   STO 01
2   STO 02

-With n = 1    1  STO 03   XEQ "DTI"   >>>>   I = 1.226398672                     ---Execution time = 48s---
-With n = 2    2  STO 03        R/S          >>>>   I = 1.226795831                     ---Execution time = 5m26s---
-With n = 4    4  STO 03        R/S          >>>>   I = 1.226799708

Example2:      Evaluate   I = §12 §xx^2   [ Ln(1+x.y)] / sqrt ( x2 + y2 ) 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

SF 02  ( double integral )

S  ASTO 00

U  ASTO 04

1   STO 01
2   STO 02

-With n = 1    1  STO 03   XEQ "DTI"   >>>>   I = 0.456387227                            ---Execution time = 15s---
-With n = 2    2  STO 03        R/S          >>>>   I = 0.456373589                            ---Execution time = 49s---
-With n = 4    4  STO 03        R/S          >>>>   I = 0.456373361                            ---Execution time = 184s---

2°)  Gauss-Legendre 4-point Formula

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

§-11 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 R05 are to be initialized before executing "DTI4" )

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

Flags: F02                        CF 02 = Triple Integral
SF 02 = Double Integral

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
A subroutine that takes x in X-register & y in Y-register and returns t(x,y) in Y-register & w(x,y) in X-register.

-Line 40 is a three-byte GTO 01

 01 LBL "DTI4"  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 24  24 /  25 STO 09  26 RCL 02  27 RCL 01  28 STO 20  29 -  30 RCL 03 31 STO 17  32 /  33 STO 22  34 2  35 STO 10  36 /  37 ST+ 20  38 CLX  39 STO 16  40 GTO 01  41 LBL 02  42 STO 06  43 XEQ IND 04  44 STO 13  45 -  46 RCL 03  47 STO 19  48 /  49 STO 14  50 RCL 10  51 /  52 ST+ 13  53 CLX  54 STO 18  55 LBL 03  56 RCL 13  57 RCL 14  58 RCL 08  59 *  60 - 61 XEQ 04  62 ST+ 18  63 RCL 13  64 RCL 14  65 RCL 07  66 *  67 -  68 XEQ 04  69 RCL 09  70 *  71 ST+ 18  72 RCL 13          73 RCL 14  74 RCL 07  75 *  76 +  77 XEQ 04  78 RCL 09  79 *  80 ST+ 18  81 RCL 13  82 RCL 14  83 ST+ 13  84 RCL 08  85 *  86 +  87 XEQ 04  88 ST+ 18  89 DSE 19  90 GTO 03 91 RCL 14  92 RCL 18  93 *  94 RTN  95 LBL 04  96 STO 15  97 RCL 06  98 FS? 02  99 GTO IND 00 100 XEQ IND 05 101 STO 11 102 - 103 RCL 03 104 STO 21 105 / 106 STO 12 107 RCL 10 108 / 109 ST+ 11 110 CLX 111 STO 23 112 LBL 05 113 RCL 11 114 RCL 12 115 RCL 08 116 * 117 - 118 RCL 15 119 RCL 06 120 XEQ IND 00 121 ST+ 23 122 RCL 11 123 RCL 12 124 RCL 07 125 * 126 - 127 RCL 15 128 RCL 06 129 XEQ IND 00 130 RCL 09 131 * 132 ST+ 23 133 RCL 11 134 RCL 12 135 RCL 07 136 * 137 + 138 RCL 15 139 RCL 06 140 XEQ IND 00 141 RCL 09 142 * 143 ST+ 23 144 RCL 11 145 RCL 12 146 ST+ 11 147 RCL 08 148 * 149 + 150 RCL 15 151 RCL 06 152 XEQ IND 00 153 ST+ 23 154 DSE 21 155 GTO 05 156 RCL 12 157 RCL 23 158 * 159 RTN 160 LBL 01 161 RCL 20 162 RCL 22 163 RCL 08 164 * 165 - 166 XEQ 02 167 ST+ 16 168 RCL 20 169 RCL 22 170 RCL 07 171 * 172 - 173 XEQ 02 174 RCL 09 175 * 176 ST+ 16 177 RCL 20 178 RCL 22 179 RCL 07 180 * 181 + 182 XEQ 02 183 RCL 09 184 * 185 ST+ 16 186 RCL 20         187 RCL 22 188 ST+ 20 189 RCL 08 190 * 191 + 192 XEQ 02 193 ST+ 16 194 DSE 17 195 GTO 01 196 RCL 16 197 RCL 22 198 * 199 RCL 24 200 2 201 / 202 FC? 02 203 ST* Y 204 X^2 205 * 206 STO 06 207 END

( 302 bytes / SIZE 025 )

 STACK INPUT OUTPUT X / Integral

Example1:      Evaluate   I = §12 §xx^2 §x+yx^2.y  [ Ln(1+x.y.z)] / sqrt ( x2 + y2 + z2 ) dx dy dz

 01 LBL "S"  02 STO T  03 ST* T  04 X<>Y  05 X^2  06 ST+ T  07 X<> L  08 *  09 X<>Y  10 X^2  11 ST+ Z  12 X<> L  13 *  14 LN1+X  15 X<>Y  16 SQRT  17 /  18 RTN  19 LBL "U"  20 X^2  21 LASTX  22 RTN  23 LBL "W"  24 STO Z  25 ST* Z  26 X<>Y  27 ST* Z  28 +  29 RTN  30 END

CF 02   ( triple integral )

S  ASTO 00

U  ASTO 04
V  ASTO 05

1   STO 01
2   STO 02

-With n = 1    1  STO 03   XEQ "DTI"   >>>>   I = 1.226803750                    ---Execution time = 1m51s---
-With n = 2    2  STO 03        R/S          >>>>   I = 1.226799889
-With n = 4    4  STO 03        R/S          >>>>   I = 1.226799707

Example2:      Evaluate   I = §12 §xx^2   [ Ln(1+x.y)] / sqrt ( x2 + y2 ) 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

SF 02  ( double integral )

S  ASTO 00

U  ASTO 04

1   STO 01
2   STO 02

-With n = 1    1  STO 03   XEQ "DTI"   >>>>   I = 0.456373416                          ---Execution time = 26s---
-With n = 2    2  STO 03        R/S          >>>>   I = 0.456373357                          ---Execution time = 90s---
-With n = 4    4  STO 03        R/S          >>>>   I = 0.456373358

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