The Museum of HP Calculators
The Museum of HP Calculators
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.
-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