# hp41programs

Multiplinteg Multiple Integrals for the HP-41

Overview

1°)  Gauss-Legendre 3-point Formula

a)  Program#1
b)  Program#2

2°)  Gauss-Legendre 4-point Formula
3°)  Gauss-Legendre 5-point Formula
4°)  Gauss-Legendre 2-point Formula
5°)  Constant Limits of Integration
6°)  Monte Carlo Integration

a)  Example1
b)  Example2

1°)  Gauss-Legendre 3-point Formula

a)  Program#1

-The following program calculates:    §ab §u1(x1)v1(x1)  §u2(x1,x2)v2(x1,x2) .........  §u(n-1)(x1,...,x(n-1))v(n-1)(x1,...,x(n-1))   f(x1,x2,....;xn)  dx1dx2....dxn   ( n < 7 )

i-e integrals , double integrals , ...... , up to sextuple integrals, provided the subroutines do not contain any XEQ.

-This limitation is only due to the fact that the return stack can accomodate up to 6 pending return addresses.

-The 3 point Gauss-Legendre formula is applied to n sub-intervals in the direction of all axis:  §-11 f(x).dx = ( 5.f(-(3/5)1/2) + 8.f(0) + 5.f((3/5)1/2) )/9
-Far from any singularity, this formula is of order 6.
-Linear changes of arguments transform the intervals into the standard interval  [ 1 ; 1 ]

Data Registers:       R00 = m =  number of  §    ;   R01 = x1 ;  R02 = x2 ; .......... ;  R06 = x6
R07 = 7 ; R08 = 4 ; R09 = 2 ; R10 = 1.6 ; R11 = 3.6 ; R12 = 0.151/2 ; R13 = 0.61/2   ;    R14 thru R17 = control numbers

R18 = n = number of sub-intervals

R19 = f name                                  R25 = u1 name         R32 = u2 name         .....................    R53 = u5 name
R20 = a                                          R26 = v1 name         R33 = v2 name        .....................     R54 = v5 name
R21 = b                                          R27 = u1(x1)            R34 = u2(x1;x2)       .....................     R55 = u5(x1;x2;...;x5)
R22 = (b-a)/n                                 R28 = v1(x1)            R35 = v2(x1;x2)       .....................      R56 = v5(x1;x2;...;x5)
R23 = index                                   R29 = (v1-u1)/n        R36 = (v2-u2)/n        .....................     R57 = (v5-u5)/n
R24 = partial sum,                          R30 = index             R37 = index             .....................      R58 = index
and, finally, the integral                          R31 = partial sum     R38 = partial sum    ......................      R59 = partial sum

Flags: /
Subroutines:  The functions  f ;  u1 ; v1 ; u2 ; v2 ; ....... are to be computed assuming  x1 is in R01 ; x2 is in R02 ; ........

-Lines 02 to 47 are useful to store the various inputs in the proper registers.
-If you initialize registers  R00 , R18 , R19 , R20 , R21   ;   R25 , R26  ;  R32 , R33  ;  ....before executing "MIT" , these lines may be deleted.

-The append character is denoted  ~

 01  LBL "MIT"   02  " A=?"   03  PROMPT   04  STO 20   05  " B=?"   06  PROMPT   07  STO 21   08  " FNAME?"   09  AON   10  STOP   11  ASTO 19   12  1   13  STO 00    14  25   15  FIX 0   16  CF 29   17  LBL 00    18  CF 23   19  " U"   20  ARCL 00   21  "~NAME?"   22  STOP   23  FC?C 23   24  GTO 05   25  ASTO IND X   26  1   27  +   28  " V"   29  ARCL 00   30  "~NAME?"   31  STOP   32  ASTO IND X 33  6   34  +   35  ISG 00   36  CLX   37  GTO 00   38  LBL 05   39  AOFF   40  FIX 9   41  SF 29   42  " N=?"   43  PROMPT    44  CLA   45  LBL 10    46  STO 18    47  RCL 00   48   E3   49  /   50  STO 14            51  15   52  STO 15   53  16   54  STO 16   55  17   56  STO 17   57  RCL 21   58  RCL 20   59  -   60  STO 22   61  7   62  STO 07   63  4   64  STO 08 65  SQRT   66  STO 09   67  1.6   68  STO 10   69  +   70  STO 11   71  .15   72  SQRT   73  STO 12   74  ST+ X   75  STO 13   76  LBL 01    77  ISG 14   78  GTO 02   79  DSE 14   80  CLX   81  GTO IND 19   82  LBL 02   83  RCL 07   84  ST+ 15   85  ST+ 16   86  ST+ 17   87  RCL 15   88  RCL 08   89  -   90  RCL IND X   91  SIGN   92  X#0?   93  GTO 03   94  XEQ IND L   95  RCL 15   96  RCL 09 97  -   98  X<>Y   99  STO IND Y 100  CHS 101  STO IND 15 102  DSE Y 103  RCL IND Y 104  XEQ IND X 105  RCL 15 106  DSE X 107  X<>Y 108  STO IND Y 109  ST+ IND 15 110  LBL 03  111  RCL 18 112  ST/ IND 15 113  STO IND 16 114  CLX 115  STO IND 17 116  LBL 04 117  RCL 15 118  RCL IND 16 119  RCL 09 120  ST- Z 121  1/X 122  - 123  RCL IND 15 124  * 125  RCL IND Y 126  + 127  STO IND 14 128  XEQ 01 129  RCL 10 130  * 131  ST+ IND 17 132  RCL IND 15 133  RCL 12 134  * 135  ST- IND 14 136  XEQ 01 137  ST+ IND 17 138  RCL IND 15 139  RCL 13  140  * 141  ST+ IND 14 142  XEQ 01 143  ST+ IND 17 144  DSE IND 16 145  GTO 04 146  RCL IND 15 147  RCL 11  148  / 149  ST* IND 17 150  RCL IND 17 151  RCL 07 152  ST- 15 153  ST- 16 154  ST- 17 155  SIGN 156  ST-14 157  X<>Y 158  RTN 159  END

( 283 bytes / SIZE 18+7m )

 STACK INPUT OUTPUT X / I

I = R24

Example:    Evaluate   I  =  §13  §xx^2  §x+yx.y  §zx+z   ln(x2+y/z+t) dx.dy.dz.dt

-First, we load the 7 required subroutines ( I've used one-character global labels to speed up execution, namely "T" , "U" , "V" , "W" , "X" , "Y" , "Z" ):

LBL "T"             f(x,y,z,t) = ln(x2+y/z+t)
RCL 01
X^2
RCL 02
RCL 03
/
+
RCL 04
+
LN
RTN
LBL "U"           u1(x) = x
RCL 01
RTN
LBL "V"           v1(x) = x2
RCL 01
X^2
RTN
LBL "W"          u2(x,y) = x+y
RCL 01
RCL 02
+
RTN
LBL "X"          v2(x,y) = x.y
RCL 01
RCL 02
*
RTN
LBL "Y"          u3(x,y,z) = z
RCL 03
RTN
LBL "Z"          v3(x,y,z) = x+z
RCL 01
RCL 03
+
RTN

-We execute SIZE 046 ( or greater ) and

XEQ "MIT"   >>>>    " A=?"
1   R/S       >>>>     " B=?"
3   R/S       >>>>     " FNAME?"
T   R/S      >>>>      " U1NAME?"
U   R/S     >>>>      " V1NAME?"
V   R/S     >>>>      " U2NAME?"
W  R/S     >>>>      " V2NAME?"
X   R/S     >>>>      " U3NAME?"
Y   R/S     >>>>      " V3NAME?"
Z   R/S     >>>>      " U4NAME?"     press  R/S   without any entry when all function names have been keyed in
R/S     >>>>      " N=?"                ( number of sub-intervals )  let's try     n = 1

1   R/S     >>>>      the result is obtained about 3mn18s later:     I = 160.452315     in X-register and in register R24

-To recalculate I with another n-value, key in this value and press XEQ 10

for example,

with  n = 2     2  XEQ 10  yields   160.631496
and   n = 4     4  XEQ 10  yields   160.634273

Notes:

-If n is multiplied by 2, execution time is approximately multiplied by 16 for quadruple integrals , by 64 for sextuple integrals.
-We can use Lagrange interpolation formula to improve our results by extrapolation to the limit.
-In this example, we find  I = 160.634317

b)  Program#2

-The following variant is much longer but significantly faster ( more exactly less slow... )

-Flags F01 to F05 are used to define single, double, .... , sextuple integrals
and instead of 2 subroutines to calculate for instance u(x) & v(x),
a unique subroutine must return v(x) in Y-register and u(x) in X-register.

-There is no PROMPT because it's much easier to remember what data registers must be initialize.

Data Registers:           •  R00 = f name               ( Registers R00 & R11 thru R12+n are to be initialized before executing "MI3" )

R01 = x1 ,  R02 = x2 , .......... ,  R06 = x6            R07 to R09: temp     R10 = Integral

•  R11 = a                                    •  R14 = u1v1 name               R19 to R17+4n:  temp
•  R12 = b                                    •  R15 = u2v2 name
•  R13 = m = Nb of subintervals   •  R16 = u3v3 name
•  R17 = u4v4 name
•  R18 = u5v5 name

Flags:  F01 to F05          SF 01 for single integrals
CF 01  SF 02  for double integrals
CF 01  CF 02  SF 03  for triple integrals
CF 01  CF 02  CF 03  SF 04  for quadruple integrals
CF 01  CF 02  CF 03  CF 04  SF 05  for quintuple integrals
CF 01  CF 02  CF 03  CF 04  CF 05 for sextuple integrals

Subroutines:  The functions  f ,  u1v1 , u2v2 , ....... are to be computed assuming  x1 is in R01 , x2 is in R02 , ........

f  returns f(x1,x2,....;xn) in X-register
u1v1 returns  v1(x1)  in Y-register  &  u1(x1)  in X-register
u2v2 returns  v2(x1,x2)  in Y-register  &  u2(x1,x2)  in X-register
..................................................................................................

un-1vn-1 returns  vn-1(x1,x2,....,xn)  in Y-register  &  un-1(x1,x2,....,xn)  in X-register

-Line 21 is a three-byte GTO 01

 01 LBL "MI3"  02 RCL 12  03 RCL 11  04 STO 19  05 -  06 RCL 13  07 STO 20  08 /  09 STO 21  10 2  11 STO 07  12 /  13 ST+ 19  14 CLX  15 STO 10  16 .15  17 SQRT  18 STO 09  19 1.6  20 STO 08  21 GTO 01  22 LBL 02  23 STO 01  24 FS? 01  25 GTO IND 00   26 XEQ IND 14  27 STO 22  28 -  29 RCL 13  30 STO 23  31 /  32 STO 24  33 RCL 07  34 /  35 ST+ 22  36 CLX  37 STO 25  38 LBL 03  39 RCL 22  40 RCL 24 41 RCL 09  42 *  43 -  44 XEQ 04  45 ST+ 25  46 RCL 22  47 XEQ 04  48 RCL 08  49 *  50 ST+ 25  51 RCL 22  52 RCL 24  53 ST+ 22  54 RCL 09  55 *  56 +  57 XEQ 04  58 ST+ 25  59 DSE 23  60 GTO 03  61 RCL 24  62 RCL 25  63 *  64 RTN  65 LBL 04  66 STO 02  67 FS? 02  68 GTO IND 00   69 XEQ IND 15  70 STO 26  71 -  72 RCL 13  73 STO 27  74 /  75 STO 28  76 RCL 07  77 /  78 ST+ 26  79 CLX  80 STO 29 81 LBL 05  82 RCL 26  83 RCL 28  84 RCL 09  85 *  86 -  87 XEQ 06  88 ST+ 29  89 RCL 26  90 XEQ 06  91 RCL 08  92 *  93 ST+ 29  94 RCL 26  95 RCL 28  96 ST+ 26  97 RCL 09  98 *  99 + 100 XEQ 06 101 ST+ 29 102 DSE 27 103 GTO 05 104 RCL 28 105 RCL 29 106 * 107 RTN 108 LBL 06 109 STO 03 110 FS? 03 111 GTO IND 00 112 XEQ IND 16 113 STO 30 114 - 115 RCL 13 116 STO 31 117 / 118 STO 32 119 RCL 07 120 / 121 ST+ 30 122 CLX 123 STO 33 124 LBL 07 125 RCL 30 126 RCL 32 127 RCL 09 128 * 129 - 130 XEQ 08 131 ST+ 33 132 RCL 30 133 XEQ 08 134 RCL 08 135 * 136 ST+ 33 137 RCL 30 138 RCL 32 139 ST+ 30 140 RCL 09 141 * 142 + 143 XEQ 08 144 ST+ 33 145 DSE 31 146 GTO 07 147 RCL 32 148 RCL 33 149 * 150 RTN 151 LBL 08 152 STO 04 153 FS? 04 154 GTO IND 00 155 XEQ IND 17 156 STO 34 157 - 158 RCL 13 159 STO 35 160 / 161 STO 36 162 RCL 07 163 / 164 ST+ 34 165 CLX 166 STO 37 167 LBL 09 168 RCL 34 169 RCL 36 170 RCL 09 171 * 172 - 173 XEQ 10 174 ST+ 37 175 RCL 34 176 XEQ 10 177 RCL 08 178 * 179 ST+ 37 180 RCL 34 181 RCL 36 182 ST+ 34 183 RCL 09 184 * 185 + 186 XEQ 10 187 ST+ 37 188 DSE 35 189 GTO 09 190 RCL 36          191 RCL 37 192 * 193 RTN 194 LBL 10 195 STO 05 196 FS? 05 197 GTO IND 00 198 XEQ IND 18 199 STO 38 200 - 201 RCL 13 202 STO 39 203 / 204 STO 40 205 RCL 07 206 / 207 ST+ 38 208 CLX 209 STO 41 210 LBL 11 211 RCL 38 212 RCL 40 213 RCL 09 214 * 215 - 216 STO 06 217 XEQ IND 00 218 ST+ 41 219 RCL 38 220 STO 06 221 XEQ IND 00 222 RCL 08 223 * 224 ST+ 41 225 RCL 38 226 RCL 40 227 ST+ 38 228 RCL 09 229 * 230 + 231 STO 06 232 XEQ IND 00 233 ST+ 41 234 DSE 39 235 GTO 11 236 RCL 40 237 RCL 41 238 * 239 RTN 240 LBL 01 241 RCL 19 242 RCL 21 243 RCL 09 244 * 245 - 246 XEQ 02 247 ST+ 10 248 RCL 19 249 XEQ 02 250 RCL 08 251 * 252 ST+ 10 253 RCL 19 254 RCL 21 255 ST+ 19 256 RCL 09 257 * 258 + 259 XEQ 02 260 ST+ 10 261 DSE 20 262 GTO 01 263 RCL 10 264 RCL 21 265 * 266 1.005 267 LBL 00 268 FS? IND X    269 GTO 12 270 ISG X 271 GTO 00 272 LBL 12 273 INT 274 3.6 275 X<>Y 276 Y^X 277 / 278 STO 10 279 END

( 457 bytes / SIZE 18+4n )

 STACK INPUT OUTPUT X / I

I = R10

Example1:     Evaluate   I  =  §13  §xx^2  §x+yx.y  §zx+z   ln(x2+y/z+t) dx.dy.dz.dt

-Now, we only have to store 4 subroutines:

 01 LBL "S"  02 RCL 01  03 X^2  04 RCL 02  05 RCL 03  06 /  07 +  08 RCL 04  09 +  10 LN  11 RTN  12 LBL "Y"  13 RCL 01  14 X^2  15 LASTX  16 RTN  17 LBL "Z"  18 RCL 02  19 RCL 02  20 RCL 01  21 ST* Z  22 +  23 RTN  24 LBL "T"  25 RCL 01  26 RCL 03  27 ST+ Y  28 RTN

"S"  ASTO 00

"Y"  ASTO 14
"Z"  ASTO 15
"T"  ASTO 16

CF 01  CF 02  CF 03  SF 04  ( quadruple integral )

1  STO 11
3  STO 12    and if we choose  m = 1  subinterval:

1  STO 13   XEQ "MI3"  >>>>   I = 160.4523151 = R10                                     ---Execution time = 2m02s---

-To recalculate I with another m-value, store this value in R13 and  R/S

with  m = 2     2  STO 13  R/S   >>>>   160.631496
and   m = 4     4  STO 13  R/S   >>>>   160.634273

Note:

-The execution time was 3m18s with the 1st version.
-So, we have saved more than 1 minute.

Example2:     Evaluate   I  =  §11.3  §xx^2  §yx.y  §zyz   §uz.u  §vu.v   Ln(1+x+y+z+u+v+w) dx.dy.dz.du.dv.dw

 01 LBL "S"  02 RCL 01  03 RCL 02  04 +  05 RCL 03  06 +  07 RCL 04  08 + 09 RCL 05  10 +  11 RCL 06  12 +  13 LN1+X  14 RTN  15 LBL "Y"  16 RCL 01 17 X^2  18 LASTX  19 RTN  20 LBL "Z"  21 RCL 01  22 RCL 02  23 ST* Y  24 RTN 25 LBL "U"  26 RCL 02  27 RCL 03  28 ST* Y  29 RTN  30 LBL "V"  31 RCL 03  32 RCL 04 33 ST* Y  34 RTN  35 LBL "W"  36 RCL 04  37 RCL 05  38 ST* Y  39 RTN  40 END

"S"  ASTO 00

"Y"  ASTO 14
"Z"  ASTO 15
"U"  ASTO 16
"V"  ASTO 17
"W" ASTO 18

1    STO 11
1.3  STO 12

CF 01  CF 02  CF 03  CF 04  CF 05        ( sextuple integral )

-If you choose  m = 1  subinterval

1  STO 13   XEQ "MI3"   >>>>    I = 0.119433886                                  ---Execution time = 18m33s---

-Likewise   2  STO 13  R/S  >>>>    I = 0.126416371
4  STO 13  R/S  >>>>    I = 0.126694346

-With m = 2 & 4, I've used V41 or free42binary

Notes:

-With Lagrange formula to use extrapolation to the limit, we obtain  I ~ 0.126699797

-If you insert the subroutines inside the program itself, and if you replace the XEQ IND & GTO IND by local XEQ,
the execution time becomes 16m22s  for m = 1  ( provided the GTOs & XEQs have been compiled )

2°)  Gauss-Legendre 4-point Formula

-Using 4-point Gauss-Legendre formula will produce of course a better precision.

§-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 & R11 thru R12+n are to be initialized before executing "MI4" )

R01 = x1 ,  R02 = x2 , .......... ,  R06 = x6            R07 to R09: temp     R10 = Integral

•  R11 = a                                    •  R14 = u1v1 name               R19 to R19+4n:  temp
•  R12 = b                                    •  R15 = u2v2 name
•  R13 = m = Nb of subintervals   •  R16 = u3v3 name
•  R17 = u4v4 name
•  R18 = u5v5 name

Flags:  F01 to F05          SF 01 for single integrals
CF 01  SF 02  for double integrals
CF 01  CF 02  SF 03  for triple integrals
CF 01  CF 02  CF 03  SF 04  for quadruple integrals
CF 01  CF 02  CF 03  CF 04  SF 05  for quintuple integrals
CF 01  CF 02  CF 03  CF 04  CF 05 for sextuple integrals

Subroutines:  The functions  f ,  u1v1 , u2v2 , ....... are to be computed assuming  x1 is in R01 , x2 is in R02 , ........

f  returns f(x1,x2,....;xn) in X-register
u1v1 returns  v1(x1)  in Y-register  &  u1(x1)  in X-register
u2v2 returns  v2(x1,x2)  in Y-register  &  u2(x1,x2)  in X-register
..................................................................................................

un-1vn-1 returns  vn-1(x1,x2,....,xn)  in Y-register  &  un-1(x1,x2,....,xn)  in X-register

-Line 41 is a three-byte GTO 01

 01 LBL "MI4"  02 RCL 12  03 RCL 11  04 STO 20  05 -  06 RCL 13  07 STO 21  08 /  09 STO 22  10 2  11 STO 10  12 /  13 ST+ 20  14 CLX  15 STO 19  16 3  17 STO Y  18 4.8  19 SQRT  20 ST+ Z  21 -  22 28  23 ST/ Z  24 /  25 SQRT  26 STO 07  27 X<>Y  28 SQRT  29 STO 08  30 .5  31 ENTER  32 STO 23  33 43.2  34 SQRT  35 1/X  36 ST+ Z  37 -  38 ST* 23  39 /  40 STO 09  41 GTO 01  42 LBL 02  43 STO 01  44 FS? 01  45 GTO IND 00   46 XEQ IND 14  47 STO 24  48 -  49 RCL 13  50 STO 25  51 /  52 STO 26  53 RCL 10  54 / 55 ST+ 24  56 CLX  57 STO 27  58 LBL 03  59 RCL 24  60 RCL 26  61 RCL 08  62 *  63 -  64 XEQ 04  65 ST+ 27  66 RCL 24  67 RCL 26  68 RCL 07  69 *  70 -  71 XEQ 04  72 RCL 09  73 *  74 ST+ 27  75 RCL 24  76 RCL 26  77 RCL 07  78 *  79 +  80 XEQ 04  81 RCL 09  82 *  83 ST+ 27  84 RCL 24  85 RCL 26  86 ST+ 24  87 RCL 08  88 *  89 +  90 XEQ 04  91 ST+ 27  92 DSE 25  93 GTO 03  94 RCL 26  95 RCL 27  96 *  97 RTN  98 LBL 04  99 STO 02 100 FS? 02 101 GTO IND 00 102 XEQ IND 15 103 STO 28 104 - 105 RCL 13 106 STO 29 107 / 108 STO 30 109 RCL 10 110 / 111 ST+ 28 112 CLX 113 STO 31 114 LBL 05 115 RCL 28 116 RCL 30 117 RCL 08 118 * 119 - 120 XEQ 06 121 ST+ 31 122 RCL 28 123 RCL 30 124 RCL 07 125 * 126 - 127 XEQ 06 128 RCL 09 129 * 130 ST+ 31 131 RCL 28 132 RCL 30 133 RCL 07 134 * 135 + 136 XEQ 06 137 RCL 09 138 * 139 ST+ 31 140 RCL 28 141 RCL 30 142 ST+ 28 143 RCL 08 144 * 145 + 146 XEQ 06 147 ST+ 31 148 DSE 29 149 GTO 05 150 RCL 30 151 RCL 31 152 * 153 RTN 154 LBL 06 155 STO 03 156 FS? 03 157 GTO IND 00 158 XEQ IND 16 159 STO 32 160 - 161 RCL 13 162 STO 33 163 / 164 STO 34 165 RCL 10 166 / 167 ST+ 32 168 CLX 169 STO 35 170 LBL 07 171 RCL 32 172 RCL 34 173 RCL 08 174 * 175 - 176 XEQ 08 177 ST+ 35 178 RCL 32 179 RCL 34 180 RCL 07 181 * 182 - 183 XEQ 08 184 RCL 09 185 * 186 ST+ 35 187 RCL 32 188 RCL 34 189 RCL 07 190 * 191 + 192 XEQ 08 193 RCL 09 194 * 195 ST+ 35 196 RCL 32 197 RCL 34 198 ST+ 32 199 RCL 08 200 * 201 + 202 XEQ 08 203 ST+ 35 204 DSE 33 205 GTO 07 206 RCL 34 207 RCL 35 208 * 209 RTN 210 LBL 08 211 STO 04 212 FS? 04 213 GTO IND 00 214 XEQ IND 17 215 STO 36 216 - 217 RCL 13 218 STO 37 219 / 220 STO 38 221 RCL 10 222 / 223 ST+ 36 224 CLX 225 STO 39 226 LBL 09 227 RCL 36 228 RCL 38 229 RCL 08 230 * 231 - 232 XEQ 10 233 ST+ 39 234 RCL 36 235 RCL 38 236 RCL 07 237 * 238 - 239 XEQ 10 240 RCL 09 241 * 242 ST+ 39 243 RCL 36 244 RCL 38 245 RCL 07 246 * 247 + 248 XEQ 10 249 RCL 09 250 * 251 ST+ 39 252 RCL 36 253 RCL 38 254 ST+ 36 255 RCL 08 256 * 257 + 258 XEQ 10 259 ST+ 39 260 DSE 37 261 GTO 09 262 RCL 38 263 RCL 39 264 * 265 RTN 266 LBL 10 267 STO 05 268 FS? 05 269 GTO IND 00 270 XEQ IND 18 271 STO 40 272 - 273 RCL 13 274 STO 41 275 / 276 STO 42 277 RCL 10 278 / 279 ST+ 40 280 CLX 281 STO 43 282 LBL 11 283 RCL 40 284 RCL 42 285 RCL 08 286 * 287 - 288 STO 06 289 XEQ IND 00 290 ST+ 43 291 RCL 40 292 RCL 42 293 RCL 07 294 * 295 - 296 STO 06 297 XEQ IND 00 298 RCL 09 299 * 300 ST+ 43 301 RCL 40 302 RCL 42 303 RCL 07 304 * 305 + 306 STO 06 307 XEQ IND 00 308 RCL 09 309 * 310 ST+ 43 311 RCL 40 312 RCL 42 313 ST+ 40 314 RCL 08 315 * 316 + 317 STO 06 318 XEQ IND 00 319 ST+ 43 320 DSE 41 321 GTO 11 322 RCL 42 323 RCL 43 324 * 325 RTN 326 LBL 01 327 RCL 20 328 RCL 22 329 RCL 08 330 * 331 - 332 XEQ 02 333 ST+ 19 334 RCL 20 335 RCL 22 336 RCL 07 337 * 338 - 339 XEQ 02 340 RCL 09 341 * 342 ST+ 19 343 RCL 20 344 RCL 22 345 RCL 07 346 * 347 + 348 XEQ 02 349 RCL 09 350 * 351 ST+ 19 352 RCL 20 353 RCL 22 354 ST+ 20 355 RCL 08 356 * 357 + 358 XEQ 02 359 ST+ 19 360 DSE 21 361 GTO 01 362 RCL 19 363 RCL 22 364 * 365 RCL 23 366 1.005 367 LBL 00 368 FS? IND X    369 GTO 12 370 ISG X 371 GTO 00 372 LBL 12 373 INT 374 Y^X 375 * 376 STO 10 377 END

( 599 bytes / SIZE 20+4n )

 STACK INPUT OUTPUT X / I

I = R10

Example:     Let's take again the 2nd example above:   I  =  §11.3  §xx^2  §yx.y  §zyz   §uz.u  §vu.v   Ln(1+x+y+z+u+v+w) dx.dy.dz.du.dv.dw

 01 LBL "S"  02 RCL 01  03 RCL 02  04 +  05 RCL 03  06 +  07 RCL 04  08 + 09 RCL 05  10 +  11 RCL 06  12 +  13 LN1+X  14 RTN  15 LBL "Y"  16 RCL 01 17 X^2  18 LASTX  19 RTN  20 LBL "Z"  21 RCL 01  22 RCL 02  23 ST* Y  24 RTN 25 LBL "U"  26 RCL 02  27 RCL 03  28 ST* Y  29 RTN  30 LBL "V"  31 RCL 03  32 RCL 04 33 ST* Y  34 RTN  35 LBL "W"  36 RCL 04  37 RCL 05  38 ST* Y  39 RTN  40 END

"S"  ASTO 00

"Y"  ASTO 14
"Z"  ASTO 15
"U"  ASTO 16
"V"  ASTO 17
"W" ASTO 18

1    STO 11
1.3  STO 12

CF 01  CF 02  CF 03  CF 04  CF 05        ( sextuple integral )

-If you choose  m = 1  subinterval

1  STO 13   XEQ "MI4"   >>>>    I = 0.126248363                                  ---Execution time = 10s---    with V41 in turbo mode

-Likewise   2  STO 13  R/S  >>>>    I = 0.126696185
4  STO 13  R/S  >>>>    I = 0.126700205

-With m = 2 & 4, I've used free42binary

Notes:

-With Lagrange formula to use extrapolation to the limit, we obtain  I ~ 0.126700221

( 4-point Gauss-Legendre formula suggests that the error is roughly inversely proportional to m8:     I = Im + k/m8 )

-With a real HP41, even with m = 1, execution time is probably about  1h40m !

3°)  Gauss-Legendre 5-point Formula

-An even better accuracy may be obtained by Gauss-Legendre 5-point formula

Data Registers:           •  R00 = f name               ( Registers R00 & R11 thru R12+n are to be initialized before executing "MI5" )

R01 = x1 ,  R02 = x2 , .......... ,  R06 = x6            R07 to R09: temp     R10 = Integral

•  R11 = a                                    •  R14 = u1v1 name               R19 to R43:  temp
•  R12 = b                                    •  R15 = u2v2 name
•  R13 = m = Nb of subintervals   •  R16 = u3v3 name
•  R17 = u4v4 name
•  R18 = u5v5 name

Flags:  F01 to F05          SF 01 for single integrals
CF 01  SF 02  for double integrals
CF 01  CF 02  SF 03  for triple integrals
CF 01  CF 02  CF 03  SF 04  for quadruple integrals
CF 01  CF 02  CF 03  CF 04  SF 05  for quintuple integrals
CF 01  CF 02  CF 03  CF 04  CF 05 for sextuple integrals

Subroutines:  The functions  f ,  u1v1 , u2v2 , ....... are to be computed assuming  x1 is in R01 , x2 is in R02 , ........

f  returns f(x1,x2,....;xn) in X-register
u1v1 returns  v1(x1)  in Y-register  &  u1(x1)  in X-register
u2v2 returns  v2(x1,x2)  in Y-register  &  u2(x1,x2)  in X-register
..................................................................................................

un-1vn-1 returns  vn-1(x1,x2,....,xn)  in Y-register  &  un-1(x1,x2,....,xn)  in X-register

-Line 24 is a three-byte GTO 01

 01 LBL "IM5"  02 RCL 12  03 RCL 11  04 STO 19  05 -  06 RCL 13  07 STO 20  08 /  09 STO 21  10 2  11 STO 43  12 /  13 ST+ 19  14 .453089923  15 STO 07  16 .2692346551  17 STO 08  18 2.020153476  19 STO 09  20 2.401115807  21 STO 10  22 CLX  23 STO 42  24 GTO 01  25 LBL 02  26 STO 01  27 FS? 01  28 GTO IND 00  29 XEQ IND 14  30 STO 22          31 -  32 RCL 13  33 STO 23  34 /  35 STO 24  36 RCL 43  37 /  38 ST+ 22  39 CLX  40 STO 25  41 LBL 03  42 RCL 22  43 RCL 24  44 RCL 07  45 *  46 -  47 XEQ 04  48 ST+ 25  49 RCL 22  50 RCL 24  51 RCL 08  52 *  53 -  54 XEQ 04  55 RCL 09  56 * 57 ST+ 25  58 RCL 22  59 XEQ 04  60 RCL 10  61 *  62 ST+ 25  63 RCL 22  64 RCL 24  65 RCL 08  66 *  67 +  68 XEQ 04  69 RCL 09  70 *  71 ST+ 25  72 RCL 22  73 RCL 24  74 ST+ 22  75 RCL 07  76 *  77 +  78 XEQ 04  79 ST+ 25  80 DSE 23  81 GTO 03  82 RCL 24  83 RCL 25  84 *  85 RTN  86 LBL 04  87 STO 02  88 FS? 02  89 GTO IND 00  90 XEQ IND 15  91 STO 26  92 -  93 RCL 13  94 STO 27  95 /  96 STO 28  97 RCL 43  98 /  99 ST+ 26 100 CLX 101 STO 29 102 LBL 05 103 RCL 26 104 RCL 28 105 RCL 07 106 * 107 - 108 XEQ 06 109 ST+ 29 110 RCL 26 111 RCL 28 112 RCL 08 113 * 114 - 115 XEQ 06 116 RCL 09 117 * 118 ST+ 29 119 RCL 26 120 XEQ 06 121 RCL 10 122 * 123 ST+ 29 124 RCL 26 125 RCL 28 126 RCL 08 127 * 128 + 129 XEQ 06 130 RCL 09 131 * 132 ST+ 29 133 RCL 26 134 RCL 28 135 ST+ 26 136 RCL 07 137 * 138 + 139 XEQ 06 140 ST+ 29 141 DSE 27 142 GTO 05 143 RCL 28 144 RCL 29 145 * 146 RTN 147 LBL 06 148 STO 03 149 FS? 03 150 GTO IND 00 151 XEQ IND 16 152 STO 30 153 - 154 RCL 13 155 STO 31 156 / 157 STO 32 158 RCL 43 159 / 160 ST+ 30 161 CLX 162 STO 33 163 LBL 07 164 RCL 30 165 RCL 32 166 RCL 07 167 * 168 - 169 XEQ 08 170 ST+ 33 171 RCL 30 172 RCL 32 173 RCL 08 174 * 175 - 176 XEQ 08 177 RCL 09 178 * 179 ST+ 33 180 RCL 30 181 XEQ 08 182 RCL 10 183 * 184 ST+ 33 185 RCL 30 186 RCL 32 187 RCL 08 188 * 189 + 190 XEQ 08 191 RCL 09 192 * 193 ST+ 33 194 RCL 30 195 RCL 32 196 ST+ 30 197 RCL 07 198 * 199 + 200 XEQ 08 201 ST+ 33 202 DSE 31 203 GTO 07 204 RCL 32 205 RCL 33 206 * 207 RTN 208 LBL 08 209 STO 04 210 FS? 04 211 GTO IND 00 212 XEQ IND 17 213 STO 34 214 - 215 RCL 13 216 STO 35 217 / 218 STO 36 219 RCL 43 220 / 221 ST+ 34 222 CLX 223 STO 37 224 LBL 09 225 RCL 34 226 RCL 36 227 RCL 07 228 * 229 - 230 XEQ 10 231 ST+ 37 232 RCL 34 233 RCL 36 234 RCL 08 235 * 236 - 237 XEQ 10 238 RCL 09 239 * 240 ST+ 37 241 RCL 34 242 XEQ 10 243 RCL 10 244 * 245 ST+ 37 246 RCL 34 247 RCL 36 248 RCL 08 249 * 250 + 251 XEQ 10 252 RCL 09 253 * 254 ST+ 37 255 RCL 34 256 RCL 36 257 ST+ 34 258 RCL 07 259 * 260 + 261 XEQ 10 262 ST+ 37 263 DSE 35 264 GTO 09 265 RCL 36 266 RCL 37 267 * 268 RTN 269 LBL 10 270 STO 05 271 FS? 05 272 GTO IND 00 273 XEQ IND 18 274 STO 38 275 - 276 RCL 13 277 STO 39 278 / 279 STO 40 280 RCL 43 281 / 282 ST+ 38 283 CLX 284 STO 41 285 LBL 11 286 RCL 38 287 RCL 40 288 RCL 07 289 * 290 - 291 STO 06 292 XEQ IND 00 293 ST+ 41 294 RCL 38 295 RCL 40 296 RCL 08 297 * 298 - 299 STO 06 300 XEQ IND 00 301 RCL 09 302 * 303 ST+ 41 304 RCL 38 305 STO 06 306 XEQ IND 00 307 RCL 10 308 * 309 ST+ 41 310 RCL 38 311 RCL 40 312 RCL 08 313 * 314 + 315 STO 06 316 XEQ IND 00 317 RCL 09 318 * 319 ST+ 41 320 RCL 38 321 RCL 40 322 ST+ 38 323 RCL 07 324 * 325 + 326 STO 06 327 XEQ IND 00 328 ST+ 41 329 DSE 39 330 GTO 11 331 RCL 40 332 RCL 41 333 * 334 RTN 335 LBL 01 336 RCL 19 337 RCL 21 338 RCL 07 339 * 340 - 341 XEQ 02 342 ST+ 42 343 RCL 19 344 RCL 21 345 RCL 08 346 * 347 - 348 XEQ 02 349 RCL 09 350 * 351 ST+ 42 352 RCL 19 353 XEQ 02 354 RCL 10 355 * 356 ST+ 42 357 RCL 19 358 RCL 21 359 RCL 08 360 * 361 + 362 XEQ 02 363 RCL 09 364 * 365 ST+ 42 366 RCL 19 367 RCL 21 368 ST+ 19 369 RCL 07 370 * 371 + 372 XEQ 02 373 ST+ 42 374 DSE 20 375 GTO 01 376 RCL 21 377 RCL 42         378 * 379 .1184634425 380 1.005 381 LBL 00 382 FS? IND X 383 GTO 12 384 ISG X 385 GTO 00 386 LBL 12 387 INT 388 Y^X 389 * 390 STO 10 391 END

( 679 bytes / SIZE 044 )

 STACK INPUT OUTPUT X / I

I = R10

Example:     Let's take again the 2nd example above:   I  =  §11.3  §xx^2  §yx.y  §zyz   §uz.u  §vu.v   Ln(1+x+y+z+u+v+w) dx.dy.dz.du.dv.dw

 01 LBL "S"  02 RCL 01  03 RCL 02  04 +  05 RCL 03  06 +  07 RCL 04  08 + 09 RCL 05  10 +  11 RCL 06  12 +  13 LN1+X  14 RTN  15 LBL "Y"  16 RCL 01 17 X^2  18 LASTX  19 RTN  20 LBL "Z"  21 RCL 01  22 RCL 02  23 ST* Y  24 RTN 25 LBL "U"  26 RCL 02  27 RCL 03  28 ST* Y  29 RTN  30 LBL "V"  31 RCL 03  32 RCL 04 33 ST* Y  34 RTN  35 LBL "W"  36 RCL 04  37 RCL 05  38 ST* Y  39 RTN  40 END

"S"  ASTO 00

"Y"  ASTO 14
"Z"  ASTO 15
"U"  ASTO 16
"V"  ASTO 17
"W" ASTO 18

1    STO 11
1.3  STO 12

CF 01  CF 02  CF 03  CF 04  CF 05        ( sextuple integral )

-If you choose  m = 1  subinterval

1  STO 13   XEQ "MI5"   >>>>    I = 0.126686001                ( With V41 )

-Likewise   2  STO 13  R/S  >>>>    I = 0.126700195

Notes:

-The exact values of the constants in lines 14-16-18-20-379  are respectively:

(1/6) SQRT ( 5 + sqrt(40/7) )
(1/6) SQRT ( 5 - sqrt(40/7) )
( 322 + sqrt (11830) ) / ( 322 - sqrt (11830) )
512 / ( 322 - sqrt (11830) )
( 322 - sqrt (11830) ) / 1800

4°)  Gauss-Legendre 2-point Formula

-On the other hand, we'll have a faster & shorter routine if we use Gauss-Legendre 2-point formula, but of course with a much lower precision:

Data Registers:           •  R00 = f name               ( Registers R00 & R11 thru R12+n are to be initialized before executing "MI2" )

R01 = x1 ,  R02 = x2 , .......... ,  R06 = x6            R07 to R09: temp     R10 = Integral

•  R11 = a                                    •  R14 = u1v1 name               R19 to R40:  temp
•  R12 = b                                    •  R15 = u2v2 name
•  R13 = m = Nb of subintervals   •  R16 = u3v3 name
•  R17 = u4v4 name
•  R18 = u5v5 name

Flags:  F01 to F05          SF 01 for single integrals
CF 01  SF 02  for double integrals
CF 01  CF 02  SF 03  for triple integrals
CF 01  CF 02  CF 03  SF 04  for quadruple integrals
CF 01  CF 02  CF 03  CF 04  SF 05  for quintuple integrals
CF 01  CF 02  CF 03  CF 04  CF 05 for sextuple integrals

Subroutines:  The functions  f ,  u1v1 , u2v2 , ....... are to be computed assuming  x1 is in R01 , x2 is in R02 , ........

f  returns f(x1,x2,....;xn) in X-register
u1v1 returns  v1(x1)  in Y-register  &  u1(x1)  in X-register
u2v2 returns  v2(x1,x2)  in Y-register  &  u2(x1,x2)  in X-register
..................................................................................................

un-1vn-1 returns  vn-1(x1,x2,....,xn)  in Y-register  &  un-1(x1,x2,....,xn)  in X-register

-Line 19 is a three-byte GTO 01

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

( 393 bytes / SIZE 041 )

 STACK INPUT OUTPUT X / I

I = R10

Example:     Let's take again the 2nd example above:   I  =  §11.3  §xx^2  §yx.y  §zyz   §uz.u  §vu.v   Ln(1+x+y+z+u+v+w) dx.dy.dz.du.dv.dw

 01 LBL "S"  02 RCL 01  03 RCL 02  04 +  05 RCL 03  06 +  07 RCL 04  08 + 09 RCL 05  10 +  11 RCL 06  12 +  13 LN1+X  14 RTN  15 LBL "Y"  16 RCL 01 17 X^2  18 LASTX  19 RTN  20 LBL "Z"  21 RCL 01  22 RCL 02  23 ST* Y  24 RTN 25 LBL "U"  26 RCL 02  27 RCL 03  28 ST* Y  29 RTN  30 LBL "V"  31 RCL 03  32 RCL 04 33 ST* Y  34 RTN  35 LBL "W"  36 RCL 04  37 RCL 05  38 ST* Y  39 RTN  40 END

"S"  ASTO 00

"Y"  ASTO 14
"Z"  ASTO 15
"U"  ASTO 16
"V"  ASTO 17
"W" ASTO 18

1    STO 11
1.3  STO 12

CF 01  CF 02  CF 03  CF 04  CF 05        ( sextuple integral )

-If you choose  m = 1  subinterval

1  STO 13   XEQ "MI2"   >>>>    I = 0.074398572

-Likewise   2  STO 13  R/S  >>>>    I = 0.118113746
4  STO 13  R/S  >>>>    I = 0.125963544
8  STO 13  R/S  >>>>    I = 0.126650084

Note:

-Even with n = 8 , the result is not very accurate...

5°)  Constant Limits of Integration

-It's of course simpler to evaluate  I =   §ab §cd .........  §kl   f(x1,x2,....,xn)  dx1dx2....dxn

-"MI2" can be used, for instance, with an m-point Gaussian formula and it works up to n = 10

Data Registers:           •  R00 = n                       ( Register R00 & Rbb thru Ree are to be initialized before executing "MI2" )

R01 = x1      R02 = x2     ..................       Rnn = xn
R11 = S1      R12 = S2     ..................   R10+nn = Sn          ( Sj = partial sums )
R21 to R30 = control numbers ,  R31 = 2.m

( Rbb-1 = 1 )  •  Rbb = x1 ,  •  Rbb+1 = w1 ,  .......................... ,  •  Ree-1 = xm ,  •  Ree = wm

x1 w1 x2 w2 .......... xm wm  = abscissas and weights of an integration formula of your choice.
Flags: /
Subroutine:  A program that computes   f(x1,x2,....;xn)   with  x1,x2,....;xn  in  R01 , R02 , ........... , Rnn

-Lines 140-148-157-165-174 are three-byte  GTO 06  GTO 07  GTO 08  GTO 09  GTO 10

 01  LBL "MI2"   02  STO 10   03  X<>Y   04  ENTER^   05  INT   06  DSE X   07  ENTER^   08  DSE X   09  X<>Y   10   E3   11  /   12  +   13  30.02   14  RDN   15  LBL 00    16  STO IND T    17  DSE T   18  GTO 00   19  SIGN   20  ST+ L   21  STO IND L    22  RCL 10    23  20.02   24  +   25  R^   26  LBL 14   27  STO IND Y   28  DSE Y   29  GTO 14   30  INT 31  LASTX   32  FRC   33   E3   34  *   35  X<>Y   36  -   37  1   38  +   39  STO 31    40  11.02   41  CLRGX   42  GTO IND 10   43  LBL 10    44  RCL IND 30   45  STO 10   46  CLX    47  STO 19   48  LBL 09   49  RCL IND 29   50  STO 09   51  CLX   52  STO 18   53  LBL 08   54  RCL IND 28   55  STO 08   56  CLX   57  STO 17   58  LBL 07   59  RCL IND 27   60  STO 07 61  CLX   62  STO 16   63  LBL 06   64  RCL IND 26   65  STO 06    66  CLX   67  STO 15   68  LBL 05   69  RCL IND 25   70  STO 05   71  CLX   72  STO 14   73  LBL 04    74  RCL IND 24   75  STO 04   76  CLX   77  STO 13   78  LBL 03   79  RCL IND 23   80  STO 03   81  CLX   82  STO 12   83  LBL 02   84  RCL IND 22   85  STO 02    86  CLX   87  STO 11   88  LBL 01   89  RCL IND 21   90  STO 01 91  XEQ IND 00   92  ISG 21   93  RCL IND 21   94  *   95  ST+ 11   96  ISG 21   97  GTO 01   98  RCL 11   99  ISG 22 100  RCL IND 22 101  * 102  ST+ 12 103  RCL 31  104  ST- 21 105  ISG 22 106  GTO 02 107  ST- 22 108  RCL 12 109  ISG 23 110  RCL IND 23 111  * 112  ST+ 13 113  ISG 23 114  GTO 03 115  RCL 13  116  ISG 24 117  RCL IND 24 118  * 119  ST+ 14 120  RCL 31 121  ST- 23 122  ISG 24 123  GTO 04 124  ST- 24 125  RCL 14 126  ISG 25 127  RCL IND 25 128  * 129  ST+ 15 130  ISG 25 131  GTO 05 132  RCL 15  133  ISG 26 134  RCL IND 26 135  * 136  ST+ 16 137  RCL 31 138  ST- 25 139  ISG 26 140  GTO 06 141  ST- 26 142  RCL 16 143  ISG 27 144  RCL IND 27 145  * 146  ST+ 17 147  ISG 27 148  GTO 07 149  RCL 17 150  ISG 28 151  RCL IND 28 152  * 153  ST+ 18 154  RCL 31 155  ST- 27 156  ISG 28 157  GTO 08 158  ST- 28 159  RCL 18  160  ISG 29 161  RCL IND 29 162  * 163  ST+ 19 164  ISG 29 165  GTO 09 166  RCL 19 167  ISG 30 168  RCL IND 30 169  * 170  ST+ 20 171  RCL 31 172  ST- 29 173  ISG 30 174  GTO 10 175  RCL 21 176  RCL 20 177  END

( 302 bytes / SIZE ??? )

 STACK INPUTS OUTPUTS Y bbb.eee bbb.eee X n I

where  bbb.eee  is the control number of the list of abscissas & weights  { x1 w1 x2 w2 .......... xm wm }    with   bbb > 32
and         n = number of variables  ,  n < 11

Example:     Evaluate     I  =  §01  §01  §01  §01   1 / ( 1 + x + y + z + t )  dx dy dz dt

-If you choose for instance the Gauss-Legendre 6-point formula, store the 12 following coefficients into R33 to R44:

R33 = 0.2386191861          R39 = -0.2386191861
R34 = 0.4679139346          R40 =  0.4679139346
R35 = 0.6612093865          R41 = -0.6612093865
R36 = 0.3607615730          R42 =  0.3607615730
R37 = 0.9324695142          R43 = -0.9324695142
R38 = 0.1713244924          R44 =  0.1713244924

-However, we have to make a change of variable so that all the limits of integration become  -1 , +1

x' = 2 x - 1 , y' = 2 y - 1 , z' = 2 z - 1 , t' = 2 t - 1  whence

I = ( 1/16 )  §-11  §-11  §-11  §-11   2 / ( 6 + x' + y' + z' + t' )  dx' dy' dz' dt'
I = ( 1/8 )  §-11  §-11  §-11  §-11   1 / ( 6 + x' + y' + z' + t' )  dx' dy' dz' dt'

 01  LBL "T"  02  RCL 01  03  RCL 02  04  +  05  RCL 03  06  +  07  RCL 04  08  +  09  RCL 45  10  +  11  1/X  12  RTN

6   STO 45  ( Line 09  RCL 45  is much faster than a digit entry line )

"T"  ASTO 00

33.044  ENTER^
4       XEQ "MI2"  >>>>   J ~  2.777151459           ---Execution time = 20m18s---

-Therefore,     I = J / 8 ~  0.347143932

Notes:

-A good emulator in turbo mode is of course very useful !
-You can check that  §01  §01  §01  §01  §01  §01   1 / ( 1 + x + y + z + u + v + w )  dx dy dz du dv dw ~  0.258610350
-But even with an emulator, n = 10 would probably impose an m-point formula with m < 6

-For example  I = §01 ............. §01  1 / ( x1 + ............. + x10 )  dx1 ............... dx10

-With the 2-point Gauss-Legendre formula  "IM2" gives           I ~ 87.45154095 / 512 = 0.170803791   in less than 4 seconds  ( with V41 )
-With the 3-point Gauss-Legendre formula  "IM2" returns        I ~ 87.45682557 / 512 = 0.170814112   in 2mn10s

6°)  Monte Carlo Integration

a)  Example1

-Here, we estimate an integral by the formula

I = § § §Volume  f(x,y,z) dx dy dz ~  Volume < f >            where    < f > = (1/N) SUM  f(xi,yj,zk)     with  N pseudo-random points   (xi,yj,zk)

and the standard deviation   s = Volume  [ ( < f2 > - < f >2 ) / N ]1/2  is also computed

>>>  s is only a rough estimation of the "probable" error.

-As an example, "MC1" calculates the integral  I  =  §01  §01  §01   1 / ( 1 + x + y + z )  dx dy dz
-Here, the volume V = 1

Data Registers:           •  R00 = random seed                ( Register R00 is to be initialized before executing "MC1" )

R01 = SUM f(xi,yj,zk)                 R03 = 9821               R05 = N
R02 = SUM [ f(xi,yj,zk) ]2           R04 = 0.211327        R06 = N , N-1 , ......... , 0

Flags: /
Subroutines: /

 01  LBL "MC1"  02  STO 05       03  STO 06  04  9821  05  STO 03   06  .211327  07  STO 04  08  CLX  09  STO 01  10  STO 02  11  LBL 01  12  RCL 00 13  RCL 03  14  *  15  RCL 04   16  +  17  FRC  18  ENTER^  19  STO 00       20  RCL 03  21  *  22  RCL 04  23  +  24  FRC 25  STO 00  26  ST+ Y  27  RCL 03   28  *  29  RCL 04  30  +  31  FRC  32  STO 00       33  +  34  1  35  +  36  1/X 37  ST+ 01  38  X^2  39  ST+ 02  40  DSE 06  41  GTO 01  42  RCL 02       43  RCL 01  44  RCL 05  45  ST/ Z  46  /  47  STO Z  48  X^2 49  -  50  RCL 05       51  /  52  SQRT  53  X<>Y  54  RTN  55  ST+ 05  56  STO 06   57  GTO 01  58  END

( 84 bytes / SIZE 007 )

 STACK INPUTS OUTPUTS Y / s X N I

where  N  is the number of sample points ,  I  =  §01  §01  §01   1 / ( 1 + x + y + z )  dx dy dz
and    s  =  standard deviation  =   [ ( < f2 > - < f >2 ) / N ]1/2

Example:    With the random seed   r = 0.1   STO 00

•  1000  XEQ "MC1"  >>>>   I = 0.41977                              ---Execution time = 16m44s--
X<>Y   s = 0.00302

•  9000        R/S         >>>>    I = 0.418162                            ( i-e with 10000 "random" points )
X<>Y    s = 0.000938

•  90000      R/S         >>>>    I = 0.417901                            ( i-e with 100000 "random" points )
X<>Y    s = 0.000295

-The exact result is  I = 0.417972076.....

Notes:

-The previous computations are kept when you press R/S
-With 9000  R/S  and  90000  R/S , I've used V41 in turbo mode...

b)  Example2

-In the example above, all the pseudo-random points belong to the domain of integration
-If it is not the case, we can choose a superset of the domain of integration.

-For instance, to calculate the triple integral over the sphere  S:  x2 + y2 + z2 <= 1

I = § § §S  dx dy dz / ( 1 +  x2 + y2 + z2 )

we choose at random N points (x,y,z)  in the cube  -1 < x < 1 , -1 < y < 1 , -1 < z < 1
but we only keep the M points inside the sphere

Data Registers:           •  R00 = random seed                  ( Register R00 is to be initialized before executing "MC2" )

R01 = SUM f(xi,yj,zk)                 R03 = 9821               R05 = N                                    R07 = M
R02 = SUM [ f(xi,yj,zk) ]2           R04 = 0.211327        R06 = N , N-1 , ......... , 0          R08 = 1

Flags: /
Subroutines: /

 01  LBL "MC2"  02  STO 05       03  STO 06  04  9821  05  STO 03   06  .211327  07  STO 04  08  CLX  09  STO 01   10  STO 02  11  STO 07  12  SIGN  13  STO 08  14  LBL 01  15  RCL 00  16  RCL 03  17  * 18  RCL 04  19  +  20  FRC  21  STO 00   22  ST+ X  23  RCL 08       24  -  25  X^2  26  RCL 00   27  RCL 03  28  *  29  RCL 04  30  +  31  FRC  32  STO 00  33  ST+ X  34  RCL 08 35  -  36  X^2  37  +  38  RCL 00   39  RCL 03       40  *  41  RCL 04   42  +  43  FRC  44  STO 00  45  ST+ X  46  RCL 08  47  -  48  X^2  49  +  50  RCL 08  51  XY  74  PI  75  4  76  *  77  3  78  /  79  ST* Z  80  *  81  RTN  82  ST+ 05  83  STO 06       84  GTO 01  85  END

( 116 bytes / SIZE 009 )

 STACK INPUTS OUTPUTS Y / s X N I

where  N  is the number of sample points ,  I  =  § § § x^2+y^2+z^2<1  dx dy dz / ( 1 +  x2 + y2 + z2 )
and    s  =  standard deviation  =  V [ ( < f2 > - < f >2 ) / M ]1/2

Example:    With the same random seed  r = 0.1   STO 00

•  1000  XEQ "MC2"  >>>>   I = 2.683650                              ---Execution time = 26mn--
X<>Y   s = 0.021178

and R07 = M = 547

•  9000        R/S         >>>>    I = 2.691756                            ( i-e with N = 10000 "random" points )
X<>Y    s = 0.006716

and R07 = M = 5307

•  90000       R/S        >>>>    I = 2.695752                            ( i-e with N = 100000 "random" points )
X<>Y    s = 0.002134

and R07 = M = 52464

-The exact result may be be found by a change of variables ( spherical coordinates )
-It yields   I = 4.PI ( 1 - PI/4 ) ~  2.696766213.....

Notes:

-The previous computations are kept when you press R/S
-With 9000  R/S  and  90000  R/S , I've used V41 in turbo mode...

-The ratio M/N is an approximation of the volume of the sphere divided by the volume of the cube
i-e  (4.PI/3) / 8 = PI/6 ~ 0.5236

-"MC1" & '"MC2" employ the same random number generator:  rn+1 = FRC ( 9821 rn + 0.211327 )

-Better - and more complicated - methods are described in "Numerical Recipes"

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