The Museum of HP Calculators


Multiple Integration for the HP-41

This program is Copyright © 2004 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

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

Program Listing
 

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.
 

  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?"     ( append 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  -----    160.634273

Note:   -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
 
 
 

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall