hp41programs

Multiplinteg Multiple Integrals for the HP-41
 

Overview
 

 1°)  Gauss-Legendre 3-point Formula
 2°)  Constant Limits of Integration
 3°)  Monte Carlo Integration

   a)  Example1
   b)  Example2
 
 

1°)  Gauss-Legendre 3-point Formula
 

-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 )
 

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
 

2°)  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 "IM2"
  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
 

3°)  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++" ...............