Dblinteg

Triple Integrals for the HP-41

Overview

1°) Gauss-Legendre 3-point Formula

a) Program#1
b) Program#2

2°) Gauss-Legendre 4-point Formula

1°) Gauss-Legendre 3-point Formula

a) Program#1

-"TIN" evaluates: I = §_a^b§_u(x)^v(x)§_w(x;y)^t(x;y) f (x;y;z) dx dy dz ( § = Integral )

Data Registers: • R00 = f name ( Registers R00 thru R07 are to be initialized before executing "TIN" )

         • R01 = a                                   R08 = I               R09 to R23: temp
         • R02 = b
         • R03 = n = the number of subintervals into which the intervals of integration are to be divided.
         • R04 = the name of the subroutine which calculates u(x).
         • R05 = the name of the subroutine which calculates v(x).
         • R06 = the name of the subroutine which calculates w(x;y)
         • R07 = the name of the subroutine which calculates t(x;y)

Flags: /
Subroutines: 5 subroutines must be keyed into program memory:

      -The first one must take x, y and z from the X-register,the Y-register and the Z-register respectively and calculate f(x;y;z).
      -The second takes x from the X-register and calculates u(x).
      -The third takes x from the X-register and calculates v(x).
      -Another one takes x and y from the X and Y registers respectively and calculates w(x;y)
      -The last one takes x and y from the X and Y registers respectively and calculates t(x;y)

N.B:

-When the program stops, the result is in X-register and in R08.
-Line 43 is a synthetic three-byte GTO ( but it can be replaced by a two-byte GTO because this line is executed only once ).

01 LBL "TIN"
02 1.6
03 STO 09
04 .15
05 SQRT
06 STO 10
07 RCL 02
08 RCL 01
09 STO 12
10 -
11 RCL 03
12 STO 21
13 /
14 STO 11
15 2
16 /
17 ST+ 12
18 CLX
19 STO 08
20 LBL 01
21 RCL 12
22 RCL 11
23 RCL 10
24 *
25 -
26 XEQ 02
27 ST+ 08
28 RCL 12
29 XEQ 02
30 RCL 09
31 *

32 ST+ 08
33 RCL 12
34 RCL 11
35 ST+ 12
36 RCL 10
37 *
38 +
39 XEQ 02
40 ST+ 08
41 DSE 21
42 GTO 01
43 GTO 06
44 LBL 02
45 STO 15
46 XEQ IND 04
47 STO 14
48 RCL 15
49 XEQ IND 05
50 RCL 14
51 -
52 RCL 03
53 STO 22
54 /
55 STO 13
56 2
57 /
58 ST+ 14
59 CLX
60 STO 19
61 LBL 03
62 RCL 14

63 RCL 13
64 RCL 10
65 *
66 -
67 XEQ 04
68 ST+ 19
69 RCL 14
70 XEQ 04
71 RCL 09
72 *
73 ST+ 19
74 RCL 14
75 RCL 13
76 ST+ 14
77 RCL 10
78 *
79 +
80 XEQ 04
81 ST+ 19
82 DSE 22
83 GTO 03
84 RCL 19
85 RCL 13
86 *
87 RTN
88 LBL 04
89 STO 17
90 RCL 15
91 XEQ IND 06
92 STO 18
93 RCL 17

94 RCL 15
95 XEQ IND 07
96 RCL 18
97 -
98 RCL 03
99 STO 23
100 /
101 STO 16
102 2
103 /
104 ST+ 18
105 CLX
106 STO 20
107 LBL 05
108 RCL 18
109 RCL 16
110 RCL 10
111 *
112 -
113 RCL 17
114 RCL 15
115 XEQ IND 00
116 ST+ 20
117 RCL 18
118 RCL 17
119 RCL 15
120 XEQ IND 00
121 RCL 09
122 *
123 ST+ 20
124 RCL 18

125 RCL 16
126 ST+ 18
127 RCL 10
128 *
129 +
130 RCL 17
131 RCL 15
132 XEQ IND 00
133 ST+ 20
134 DSE 23
135 GTO 05
136 RCL 20
137 RCL 16
138 *
139 RTN
140 LBL 06
141 RCL 08
142 RCL 11
143 *
144 RCL 09
145 2
146 +
147 3
148 Y^X
149 /
150 STO 08
151 END

( 226 bytes / SIZE 024 )

STACK	INPUT	OUTPUT
X	/	Integral

Example: Evaluate I = §₁²§_x^{x^2}§_x+y^xy xyz / sqrt ( x² + y² + z² ) dx dy dz

The 5 required subroutines are, for instance:
( with global labels of 6 characters maximum )

01 LBL "FF"
02 ENTER^
03 X^2
04 R^
05 ST* Z
06 X^2
07 +
08 R^
09 ST* Z
10 X^2
11 +
12 SQRT
13 /
14 RTN
15 LBL "U"
16 RTN
17 LBL "V"
18 X^2
19 RTN
20 LBL "W"
21 +
22 RTN
23 LBL "T"
24 *
25 RTN

( All the calculations fit into the stack, but if the functions are very complicated,
one may of course use data registers R24 ...etc... )

Then we initialize registers R00 thru R07:

    "FF"   ASTO 00
      1        STO 01
      2        STO 02
     "U"    ASTO 04
     "V"    ASTO 05
    "W"    ASTO 06
     "T"    ASTO 07

         n = 1      1 STO 03    XEQ "TIN" >>>>    I₁ = 0.765014888
         n = 2      2 STO 03         R/S         >>>>    I₂ = 0.770640690    ( in 4mn 29s )
         n = 4      4 STO 03         R/S         >>>>    I₄ = 0.770731245
         n = 8      8 STO 03         R/S         >>>>    I₈ = 0.770732669    ( in 4h02mn )                   ( the exact value is 0.7707326901 to ten places )

-Unfortunately, with triple integrals, when n is multiplied by 2, execution time is roughly multiplied by 8 !

-An HP-48GX ( in 10 FIX format ) gives 0.7707326900 in 26 hours with the built-in ò function !

In 7 FIX format, the result 0.7707326903 is obtained in 3h49mn.
In 6 FIX format, the result 0.7707326923 is obtained in 1h38mn.

-Thus, "TIN" is not so good as the ò function of the HP-48, but the HP-41 still bears comparison ...

-Execution time can be reduce by keying the 5 subroutines inside the "TIN" program itself ( with local labels ).
The program will run about 17% faster.

-Extrapolation to the limit may also be used to improve accuracy:

There are some theoretical reasons to think that the error is roughly inversely proportional to n⁶: I = I_n + k/n⁶
Applying this formula with n = 4 and n = 8 leads to a system that yields: I = 0.7707326916 the error is only 1.5 10^-9 .

b) Program#2

-It's in fact preferable to use just 2 subroutines instead of 4 to calculate u(x) & v(x) and w(x,y) & t(x,y)
-We can also use a flag to evaluate also double integrals with the same program:

Data Registers: • R00 = f name ( Registers R00 thru R05 are to be initialized before executing "DTI" )

                                      • R01 = a                                    • R04 = uv name       R06 = I      R07 to R23: temp
                                      • R02 = b                                    • R05 = tw name
                                      • R03 = n = Nb of subintervals

Flags: F02 CF 02 = Triple Integral
SF 02 = Double Integral

Subroutines:    A subroutine that takes x, y and z from the X-register,the Y-register and the Z-register respectively and returns f(x;y;z) in X-register.
                          A subroutine that takes x in X-register and returns v(x) in Y-register & u(x) in X-register
                          A subroutine that takes x in X-register & y in Y-register and returns t(x,y) in Y-register & w(x,y) in X-register.

-Line 21 is a three-byte GTO 01

01 LBL "DTI"
02 .15
03 SQRT
04 STO 09
05 1.6
06 STO 11
07 RCL 02
08 RCL 01
09 STO 21
10 -
11 RCL 03
12 STO 17
13 /
14 STO 10
15 2
16 STO 06
17 /
18 ST+ 21
19 CLX
20 STO 22
21 GTO 01
22 LBL 02
23 STO 14
24 XEQ IND 04
25 STO 13
26 -
27 RCL 03
28 STO 18
29 /

30 STO 12
31 RCL 06
32 /
33 ST+ 13
34 CLX
35 STO 16
36 LBL 03
37 RCL 13
38 RCL 12
39 RCL 09
40 *
41 -
42 XEQ 04
43 ST+ 16
44 RCL 13
45 XEQ 04
46 RCL 11
47 *
48 ST+ 16
49 RCL 13
50 RCL 12
51 ST+ 13
52 RCL 09
53 *
54 +
55 XEQ 04
56 ST+ 16
57 DSE 18
58 GTO 03

59 RCL 12
60 RCL 16
61 *
62 RTN
63 LBL 04
64 STO 15
65 RCL 14
66 FS? 02
67 GTO IND 00
68 XEQ IND 05
69 STO 07
70 -
71 RCL 03
72 STO 20
73 /
74 STO 08
75 RCL 06
76 /
77 ST+ 07
78 CLX
79 STO 19
80 LBL 05
81 RCL 07
82 RCL 08
83 RCL 09
84 *
85 -
86 RCL 15
87 RCL 14

88 XEQ IND 00
89 ST+ 19
90 RCL 07
91 RCL 15
92 RCL 14
93 XEQ IND 00
94 RCL 11
95 *
96 ST+ 19
97 RCL 07
98 RCL 08
99 ST+ 07
100 RCL 09
101 *
102 +
103 RCL 15
104 RCL 14
105 XEQ IND 00
106 ST+ 19
107 DSE 20
108 GTO 05
109 RCL 08
110 RCL 19
111 *
112 RTN
113 LBL 01
114 RCL 21
115 RCL 10
116 RCL 09

117 *
118 -
119 XEQ 02
120 ST+ 22
121 RCL 21
122 XEQ 02
123 RCL 11
124 *
125 ST+ 22
126 RCL 21
127 RCL 10
128 ST+ 21
129 RCL 09
130 *
131 +
132 XEQ 02
133 ST+ 22
134 DSE 17
135 GTO 01
136 RCL 10
137 RCL 22
138 *
139 3.6
140 FC? 02
141 ST/ Y
142 X^2
143 /
144 STO 06
145 END

( 216 bytes / SIZE 023 )

STACK	INPUT	OUTPUT
X	/	Integral

Example1: Evaluate I = §₁²§_x^{x^2}§_x+y^{x^2.y} [ Ln(1+x.y.z)] / sqrt ( x² + y² + z² ) dx dy dz

01 LBL "S"
02 STO T
03 ST* T
04 X<>Y
05 X^2
06 ST+ T
07 X<> L
08 *
09 X<>Y
10 X^2
11 ST+ Z
12 X<> L
13 *
14 LN1+X
15 X<>Y
16 SQRT
17 /
18 RTN
19 LBL "U"
20 X^2
21 LASTX
22 RTN
23 LBL "W"
24 STO Z
25 ST* Z
26 X<>Y
27 ST* Z
28 +
29 RTN
30 END

CF 02 ( triple integral )

S ASTO 00

U ASTO 04
V ASTO 05

1 STO 01
2 STO 02

-With n = 1    1 STO 03   XEQ "DTI"   >>>>   I = 1.226398672                     ---Execution time = 48s---
-With n = 2    2 STO 03        R/S          >>>>   I = 1.226795831                     ---Execution time = 5m26s---
-With n = 4    4 STO 03        R/S          >>>>   I = 1.226799708

Example2: Evaluate I = §₁²§_x^{x^2} [ Ln(1+x.y)] / sqrt ( x² + y² ) dx dy

01 LBL "S"
02 STO Z
03 ST* Z
04 X<>Y
05 X^2
06 ST+ Z
07 X<> L
08 *
09 LN1+X
10 X<>Y
11 SQRT
12 /
13 RTN
14 LBL "U"
15 X^2
16 LASTX
17 RTN
18 END

SF 02 ( double integral )

S ASTO 00

U ASTO 04

1 STO 01
2 STO 02

-With n = 1    1 STO 03   XEQ "DTI"   >>>>   I = 0.456387227                            ---Execution time = 15s---
-With n = 2    2 STO 03        R/S          >>>>   I = 0.456373589                            ---Execution time = 49s---
-With n = 4    4 STO 03        R/S          >>>>   I = 0.456373361                            ---Execution time = 184s---

2°) Gauss-Legendre 4-point Formula

-A better accuracy may of course be obtained by Gauss-Legendre 4-point formula:

§_-1¹ 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 thru R05 are to be initialized before executing "DTI4" )

                                      • R01 = a                                    • R04 = uv name       R06 = I      R07 to R24: temp
                                      • R02 = b                                    • R05 = tw name
                                      • R03 = n = Nb of subintervals

Flags: F02 CF 02 = Triple Integral
SF 02 = Double Integral

-Line 40 is a three-byte GTO 01

01 LBL "DTI4"
02 3
03 STO Y
04 4.8
05 SQRT
06 ST+ Z
07 -
08 28
09 ST/ Z
10 /
11 SQRT
12 STO 07
13 X<>Y
14 SQRT
15 STO 08
16 .5
17 STO Y
18 43.2
19 SQRT
20 1/X
21 ST+ Z
22 -
23 STO 24
24 /
25 STO 09
26 RCL 02
27 RCL 01
28 STO 20
29 -
30 RCL 03

31 STO 17
32 /
33 STO 22
34 2
35 STO 10
36 /
37 ST+ 20
38 CLX
39 STO 16
40 GTO 01
41 LBL 02
42 STO 06
43 XEQ IND 04
44 STO 13
45 -
46 RCL 03
47 STO 19
48 /
49 STO 14
50 RCL 10
51 /
52 ST+ 13
53 CLX
54 STO 18
55 LBL 03
56 RCL 13
57 RCL 14
58 RCL 08
59 *
60 -

61 XEQ 04
62 ST+ 18
63 RCL 13
64 RCL 14
65 RCL 07
66 *
67 -
68 XEQ 04
69 RCL 09
70 *
71 ST+ 18
72 RCL 13
73 RCL 14
74 RCL 07
75 *
76 +
77 XEQ 04
78 RCL 09
79 *
80 ST+ 18
81 RCL 13
82 RCL 14
83 ST+ 13
84 RCL 08
85 *
86 +
87 XEQ 04
88 ST+ 18
89 DSE 19
90 GTO 03

91 RCL 14
92 RCL 18
93 *
94 RTN
95 LBL 04
96 STO 15
97 RCL 06
98 FS? 02
99 GTO IND 00
100 XEQ IND 05
101 STO 11
102 -
103 RCL 03
104 STO 21
105 /
106 STO 12
107 RCL 10
108 /
109 ST+ 11
110 CLX
111 STO 23
112 LBL 05
113 RCL 11
114 RCL 12
115 RCL 08
116 *
117 -
118 RCL 15
119 RCL 06
120 XEQ IND 00

121 ST+ 23
122 RCL 11
123 RCL 12
124 RCL 07
125 *
126 -
127 RCL 15
128 RCL 06
129 XEQ IND 00
130 RCL 09
131 *
132 ST+ 23
133 RCL 11
134 RCL 12
135 RCL 07
136 *
137 +
138 RCL 15
139 RCL 06
140 XEQ IND 00
141 RCL 09
142 *
143 ST+ 23
144 RCL 11
145 RCL 12
146 ST+ 11
147 RCL 08
148 *
149 +

150 RCL 15
151 RCL 06
152 XEQ IND 00
153 ST+ 23
154 DSE 21
155 GTO 05
156 RCL 12
157 RCL 23
158 *
159 RTN
160 LBL 01
161 RCL 20
162 RCL 22
163 RCL 08
164 *
165 -
166 XEQ 02
167 ST+ 16
168 RCL 20
169 RCL 22
170 RCL 07
171 *
172 -
173 XEQ 02
174 RCL 09
175 *
176 ST+ 16
177 RCL 20
178 RCL 22

179 RCL 07
180 *
181 +
182 XEQ 02
183 RCL 09
184 *
185 ST+ 16
186 RCL 20
187 RCL 22
188 ST+ 20
189 RCL 08
190 *
191 +
192 XEQ 02
193 ST+ 16
194 DSE 17
195 GTO 01
196 RCL 16
197 RCL 22
198 *
199 RCL 24
200 2
201 /
202 FC? 02
203 ST* Y
204 X^2
205 *
206 STO 06
207 END

( 302 bytes / SIZE 025 )

STACK	INPUT	OUTPUT
X	/	Integral

Example1: Evaluate I = §₁²§_x^{x^2}§_x+y^{x^2.y} [ Ln(1+x.y.z)] / sqrt ( x² + y² + z² ) dx dy dz

CF 02 ( triple integral )

S ASTO 00

U ASTO 04
V ASTO 05

1 STO 01
2 STO 02

-With n = 1    1 STO 03   XEQ "DTI"   >>>>   I = 1.226803750                    ---Execution time = 1m51s---
-With n = 2    2 STO 03        R/S          >>>>   I = 1.226799889
-With n = 4    4 STO 03        R/S          >>>>   I = 1.226799707

Example2: Evaluate I = §₁²§_x^{x^2} [ Ln(1+x.y)] / sqrt ( x² + y² ) dx dy

01 LBL "S"
02 STO Z
03 ST* Z
04 X<>Y
05 X^2
06 ST+ Z
07 X<> L
08 *
09 LN1+X
10 X<>Y
11 SQRT
12 /
13 RTN
14 LBL "U"
15 X^2
16 LASTX
17 RTN
18 END

SF 02 ( double integral )

S ASTO 00

U ASTO 04

1 STO 01
2 STO 02

-With n = 1    1 STO 03   XEQ "DTI"   >>>>   I = 0.456373416                          ---Execution time = 26s---
-With n = 2    2 STO 03        R/S          >>>>   I = 0.456373357                          ---Execution time = 90s---
-With n = 4    4 STO 03        R/S          >>>>   I = 0.456373358

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++" ...............

hp41programs

Triple Integrals for the HP-41