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