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

-Load for instance:
 
 

 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'

-Load the short routine below:
 
 

 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  X<Y?		  52  GTO 02
 53  ST+ 07
 54  +
 55  1/X
 56  ST+ 01
 57  X^2
 58  ST+ 02
 59  LBL 02
 60  DSE 06
 61  GTO 01
 62  RCL 02     
 63  RCL 01
 64  RCL 07
 65  ST/ Z
 66  /
 67  STO Z
 68  X^2		  69  -
 70  RCL 07
 71  /
 72  SQRT
 73  X<>Y
 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++" ...............