SPHINTEG

Integrals over a Sphere & its Surface for the HP-41

Overview

1°) Integrals over a Sphere

1-a) Program#1
1-b) Program#2

2°) Integrals on the Surface of a Sphere

2-a) Program#1
2-b) Program#2

3°) Integrals over the Unit Hypersphere

3-a) Program#1 ( 5th degree of accuracy )
3-b) Program#2 ( 7th degree of accuracy )

4°) Integrals on the Surface of the Unit Hypersphere

4-a) Program#1 ( 5th degree of accuracy )
4-b) Program#2 ( 7th degree of accuracy )

-We could employ usual formulae - like Gaussian integration - in several directions
-But specific formulae for spheres and hyperspheres require much less evaluations of the function.

1°) Integrals over a Sphere

a) Program#1

-We want to estimate I = §§§_{x^2+y^2+z^2<=R^2} f(x,y,z) dx dy dz

Data Registers: • R00 = function name ( Register R00 is to be initialized before executing "I3SPH" )

R01 = Integral R02 = R
Flags: /
Subroutine: A program that takes x , y , z in registers X , Y , Z respectively and returns f(x,y,z) in X-register.

01 LBL "I3SPH"
02 STO 02
03 CLST
04 XEQ IND 00
05 4
06 *
07 STO 01
08 CLST
09 RCL 02
10 XEQ IND 00

11 ST+ 01
12 CLST
13 RCL 02
14 CHS
15 XEQ IND 00
16 ST+ 01
17 CLST
18 RCL 02
19 X<>Y
20 XEQ IND 00

21 ST+ 01
22 CLST
23 RCL 02
24 CHS
25 X<>Y
26 XEQ IND 00
27 ST+ 01
28 RCL 02
29 0
30 ENTER

31 XEQ IND 00
32 ST+ 01
33 RCL 02
34 CHS
35 0
36 ENTER
37 XEQ IND 00
38 ST+ 01
39 RCL 02
40 3

41 Y^X
42 .4
43 *
44 PI
45 *
46 3
47 /
48 ST* 01
49 RCL 01
50 END

( 75 bytes / SIZE 003 )

STACK	INPUT	OUTPUT
X	R	Integral

where R = radius of the sphere

Example: I = §§§_{x^2+y^2+z^2<=R^2} Ln ( 16 + x + y² + z³ ) dx dy dz with R = 2

-This short routine calculates Ln ( 16 + x + y² + z³ )

01 LBL "T"
02 X<>Y
03 X^2
04 +
05 X<>Y
06 3
07 Y^X
08 +
09 16
10 +
11 LN
12 RTN

"T" ASTO 00

2 XEQ "I3SPH" >>>> I = 93.3891

b) Program#2

-The same integral I = §§§_{x^2+y^2+z^2<=R^2} f(x,y,z) dx dy dz is evaluated with a better accuracy

Data Registers: • R00 = function name ( Register R00 is to be initialized before executing "I3SPH" )

R01 = Integral R02 = R R03-R04: temp
Flags: /
Subroutine: A program that takes x , y , z in registers X , Y , Z respectively and returns f(x,y,z) in X-register.

01 LBL "I3SPH"
02 STO 02
03 CLST
04 XEQ IND 00
05 .4156003483
06 *
07 STO 01
08 RCL 02
09 .8326956271
10 *
11 STO 03
12 0
13 ENTER
14 XEQ IND 00
15 STO 04
16 RCL 03
17 CHS
18 0
19 ENTER
20 XEQ IND 00
21 ST+ 04
22 0
23 RCL 03
24 0
25 XEQ IND 00
26 ST+ 04
27 0
28 RCL 03
29 CHS
30 0
31 XEQ IND 00

32 ST+ 04
33 CLST
34 RCL 03
35 XEQ IND 00
36 ST+ 04
37 CLST
38 RCL 03
39 CHS
40 XEQ IND 00
41 RCL 04
42 +
43 .1994483079
44 *
45 ST+ 01
46 0
47 RCL 02
48 .7476506947
49 *
50 STO 03
51 ENTER
52 XEQ IND 00
53 STO 04
54 0
55 RCL 03
56 ENTER
57 CHS
58 XEQ IND 00
59 ST+ 04
60 0
61 RCL 03
62 CHS

63 RCL 03
64 XEQ IND 00
65 ST+ 04
66 0
67 RCL 03
68 CHS
69 ENTER
70 XEQ IND 00
71 ST+ 04
72 RCL 03
73 0
74 RCL 03
75 XEQ IND 00
76 ST+ 04
77 RCL 03
78 0
79 RCL 03
80 CHS
81 XEQ IND 00
82 ST+ 04
83 RCL 03
84 CHS
85 0
86 RCL 03
87 XEQ IND 00
88 ST+ 04
89 RCL 03
90 CHS
91 0
92 RCL Y
93 XEQ IND 00

94 ST+ 04
95 RCL 03
96 RCL 03
97 0
98 XEQ IND 00
99 ST+ 04
100 RCL 03
101 ENTER
102 CHS
103 0
104 XEQ IND 00
105 ST+ 04
106 RCL 03
107 CHS
108 RCL 03
109 0
110 XEQ IND 00
111 ST+ 04
112 RCL 03
113 CHS
114 ENTER
115 ENTER
116 CLX
117 XEQ IND 00
118 RCL 04
119 +
120 38.06761012 E-3
121 *
122 ST+ 01
123 RCL 02
124 .4294549988

125 *
126 STO 03
127 ENTER
128 ENTER
129 XEQ IND 00
130 STO 04
131 RCL 03
132 ENTER
133 ENTER
134 CHS
135 XEQ IND 00
136 ST+ 04
137 RCL 03
138 ENTER
139 CHS
140 RCL 03
141 XEQ IND 00
142 ST+ 04
143 RCL 03
144 ENTER
145 CHS
146 ENTER
147 XEQ IND 00
148 ST+ 04
149 RCL 03
150 CHS
151 RCL 03
152 ENTER
153 XEQ IND 00
154 ST+ 04
155 RCL 03

156 CHS
157 RCL 03
158 RCL Y
159 XEQ IND 00
160 ST+ 04
161 RCL 03
162 CHS
163 STO Y
164 RCL 03
165 XEQ IND 00
166 ST+ 04
167 RCL 03
168 CHS
169 ENTER
170 ENTER
171 XEQ IND 00
172 RCL 04
173 +
174 .264961086
175 *
176 RCL 01
177 +
178 RCL 02
179 3
180 Y^X
181 *
182 STO 01
183 END

( 317 bytes / SIZE 005 )

STACK	INPUT	OUTPUT
X	R	Integral

where R = radius of the sphere

Example: I = §§§_{x^2+y^2+z^2<=R^2} Ln ( 16 + x + y² + z³ ) dx dy dz with R = 2

With the same subroutine,

"T" ASTO 00

2 XEQ "I3SPH" >>>> I = 94.2545 ( instead of 93.3891 with the 1st program )

-A substantial improvement since the exact result - rounded to 4 decimals - is 94.2541

2°) Integrals on the Surface of a Sphere

a) Program#1

-Now we want to estimate I = §§_{x^2+y^2+z^2=R^2} f(x,y,z) ds where ds is the element of surface on a sphere

Data Registers: • R00 = function name ( Register R00 is to be initialized before executing "I2SPH" )

R01 = Integral R02 = R R03-R04: temp
Flags: /
Subroutine: A program that takes x , y , z in registers X , Y , Z respectively and returns f(x,y,z) in X-register.

01 LBL "I2SPH"
02 STO 02
03 3
04 SQRT
05 /
06 STO 03
07 ENTER
08 ENTER
09 XEQ IND 00
10 STO 01
11 RCL 03
12 ENTER
13 ENTER
14 CHS
15 XEQ IND 00
16 ST+ 01
17 RCL 03
18 ENTER
19 CHS
20 RCL 03
21 XEQ IND 00
22 ST+ 01
23 RCL 03
24 ENTER
25 CHS
26 ENTER
27 XEQ IND 00
28 ST+ 01
29 RCL 03
30 CHS

31 RCL 03
32 ENTER
33 XEQ IND 00
34 ST+ 01
35 RCL 03
36 CHS
37 RCL 03
38 ENTER
39 CHS
40 XEQ IND 00
41 ST+ 01
42 RCL 03
43 CHS
44 STO Y
45 RCL 03
46 XEQ IND 00
47 ST+ 01
48 RCL 03
49 CHS
50 ENTER
51 ENTER
52 XEQ IND 00
53 ST+ 01
54 27
55 ST* 01
56 RCL 02
57 2
58 SQRT
59 /
60 ENTER

61 STO 03
62 0
63 XEQ IND 00
64 STO 04
65 RCL 03
66 ENTER
67 CHS
68 0
69 XEQ IND 00
70 ST+ 04
71 RCL 03
72 CHS
73 RCL 03
74 0
75 XEQ IND 00
76 ST+ 04
77 RCL 03
78 CHS
79 STO Y
80 0
81 XEQ IND 00
82 ST+ 04
83 RCL 03
84 0
85 RCL 03
86 XEQ IND 00
87 ST+ 04
88 RCL 03
89 0
90 RCL 03

91 CHS
92 XEQ IND 00
93 ST+ 04
94 RCL 03
95 CHS
96 0
97 RCL 03
98 XEQ IND 00
99 ST+ 04
100 RCL 03
101 CHS
102 0
103 RCL Y
104 XEQ IND 00
105 ST+ 04
106 CLST
107 RCL 03
108 ENTER
109 XEQ IND 00
110 ST+ 04
111 CLST
112 RCL 03
113 ENTER
114 CHS
115 XEQ IND 00
116 ST+ 04
117 CLST
118 RCL 03
119 CHS
120 RCL 03

121 XEQ IND 00
122 ST+ 04
123 CLST
124 RCL 03
125 CHS
126 ENTER
127 XEQ IND 00
128 RCL 04
129 +
130 32
131 *
132 ST+ 01
133 CLST
134 RCL 02
135 XEQ IND 00
136 STO 04
137 CLST
138 RCL 02
139 CHS
140 XEQ IND 00
141 ST+ 04
142 CLST
143 RCL 02
144 X<>Y
145 XEQ IND 00
146 ST+ 04
147 CLST
148 RCL 02
149 CHS
150 X<>Y

151 XEQ IND 00
152 ST+ 04
153 RCL 02
154 0
155 ENTER
156 XEQ IND 00
157 ST+ 04
158 RCL 02
159 CHS
160 0
161 ENTER
162 XEQ IND 00
163 RCL 04
164 +
165 40
166 *
167 ST+ 01
168 RCL 02
169 X^2
170 PI
171 *
172 210
173 /
174 ST* 01
175 RCL 01
176 END

( 245 bytes / SIZE 005 )

STACK	INPUT	OUTPUT
X	R	Integral

where R = radius of the sphere

Example: I = §§_{x^2+y^2+z^2=R^2} Ln ( 16 + x + y² + z³ ) dx dy dz with R = 2

With the same subroutine,

"T" ASTO 00

2 XEQ "I2SPH" >>>> I = 142.1916

b) Program#2

-The same integral is now estimated with a better precision

Data Registers: • R00 = function name ( Register R00 is to be initialized before executing "I2SPH" )

R01 = Integral R02 = R R03-R04: temp
Flags: /
Subroutine: A program that takes x , y , z in registers X , Y , Z respectively and returns f(x,y,z) in X-register.

01 LBL "I2SPH"
02 STO 02
03 0
04 ENTER
05 XEQ IND 00
06 STO 01
07 RCL 02
08 CHS
09 0
10 ENTER
11 XEQ IND 00
12 ST+ 01
13 0
14 RCL 02
15 0
16 XEQ IND 00
17 ST+ 01
18 0
19 RCL 02
20 CHS
21 0
22 XEQ IND 00
23 ST+ 01
24 CLST
25 RCL 02
26 XEQ IND 00
27 ST+ 01
28 CLST
29 RCL 02
30 CHS
31 XEQ IND 00
32 RCL 01
33 +
34 2.652142441
35 *
36 STO 01
37 RCL 02
38 2
39 SQRT
40 /
41 ENTER

42 STO 03
43 0
44 XEQ IND 00
45 STO 04
46 RCL 03
47 CHS
48 RCL 03
49 0
50 XEQ IND 00
51 ST+ 04
52 RCL 03
53 ENTER
54 CHS
55 0
56 XEQ IND 00
57 ST+ 04
58 RCL 03
59 CHS
60 STO Y
61 0
62 XEQ IND 00
63 ST+ 04
64 RCL 03
65 0
66 RCL 03
67 XEQ IND 00
68 ST+ 04
69 RCL 03
70 CHS
71 0
72 RCL 03
73 XEQ IND 00
74 ST+ 04
75 RCL 03
76 0
77 RCL 03
78 CHS
79 XEQ IND 00
80 ST+ 04
81 RCL 03
82 CHS

83 0
84 RCL Y
85 XEQ IND 00
86 ST+ 04
87 CLST
88 RCL 03
89 ENTER
90 XEQ IND 00
91 ST+ 04
92 CLST
93 RCL 03
94 CHS
95 RCL 03
96 XEQ IND 00
97 ST+ 04
98 CLST
99 RCL 03
100 ENTER
101 CHS
102 XEQ IND 00
103 ST+ 04
104 CLST
105 RCL 03
106 CHS
107 ENTER
108 XEQ IND 00
109 RCL 04
110 +
111 1.993014763
112 *
113 ST+ 01
114 RCL 02
115 .3879073041
116 *
117 ENTER
118 STO 03
119 RCL 02
120 .8360955967
121 *
122 STO 04
123 XEQ IND 00

124 STO 05
125 RCL 03
126 CHS
127 RCL 03
128 RCL 04
129 XEQ IND 00
130 ST+ 05
131 RCL 03
132 RCL 03
133 RCL 04
134 CHS
135 XEQ IND 00
136 ST+ 05
137 RCL 03
138 CHS
139 RCL 03
140 RCL 04
141 CHS
142 XEQ IND 00
143 ST+ 05
144 RCL 03
145 ENTER
146 CHS
147 RCL 04
148 XEQ IND 00
149 ST+ 05
150 RCL 03
151 CHS
152 STO Y
153 RCL 04
154 XEQ IND 00
155 ST+ 05
156 RCL 03
157 ENTER
158 CHS
159 RCL 04
160 CHS
161 XEQ IND 00
162 ST+ 05
163 RCL 03
164 CHS

165 STO Y
166 RCL 04
167 CHS
168 XEQ IND 00
169 ST+ 05
170 RCL 03
171 RCL 04
172 RCL 03
173 XEQ IND 00
174 ST+ 05
175 RCL 03
176 CHS
177 RCL 04
178 RCL 03
179 XEQ IND 00
180 ST+ 05
181 RCL 03
182 RCL 04
183 RCL 03
184 CHS
185 XEQ IND 00
186 ST+ 05
187 RCL 03
188 CHS
189 RCL 04
190 RCL 03
191 CHS
192 XEQ IND 00
193 ST+ 05
194 RCL 03
195 RCL 04
196 CHS
197 RCL 03
198 XEQ IND 00
199 ST+ 05
200 RCL 03
201 CHS
202 RCL 04
203 CHS
204 RCL 03
205 XEQ IND 00

206 ST+ 05
207 RCL 03
208 RCL 04
209 CHS
210 RCL 03
211 CHS
212 XEQ IND 00
213 ST+ 05
214 RCL 03
215 CHS
216 RCL 04
217 CHS
218 RCL Y
219 XEQ IND 00
220 ST+ 05
221 RCL 04
222 RCL 03
223 ENTER
224 XEQ IND 00
225 ST+ 05
226 RCL 04
227 CHS
228 RCL 03
229 ENTER
230 XEQ IND 00
231 ST+ 05
232 RCL 04
233 RCL 03
234 ENTER
235 CHS
236 XEQ IND 00
237 ST+ 05
238 RCL 04
239 CHS
240 RCL 03
241 ENTER
242 CHS
243 XEQ IND 00
244 ST+ 05
245 RCL 04
246 RCL 03

247 CHS
248 RCL 03
249 XEQ IND 00
250 ST+ 05
251 RCL 04
252 CHS
253 RCL 03
254 CHS
255 RCL 03
256 XEQ IND 00
257 ST+ 05
258 RCL 04
259 RCL 03
260 CHS
261 ENTER
262 XEQ IND 00
263 ST+ 05
264 RCL 04
265 CHS
266 RCL 03
267 CHS
268 ENTER
269 XEQ IND 00
270 RCL 05
271 +
272 2.507123675
273 *
274 ST+ 01
275 RCL 02
276 X^2
277 PI
278 *
279 25
280 /
281 ST* 01
282 RCL 01
283 END

( 430 bytes / SIZE 006 )

STACK	INPUT	OUTPUT
X	R	Integral

where R = radius of the sphere

Example: I = §§_{x^2+y^2+z^2=R^2} Ln ( 16 + x + y² + z³ ) dx dy dz with R = 2

With the same subroutine,

"T" ASTO 00

2 XEQ "I2SPH" >>>> I = 142.2271 ( instead of 142.1916 with the 1st program )

-The exact value is I = 142.2311 rounded to 4 decimals.

3°) Integrals over the Unit HyperSphere

a) Program#1 ( 5th degree of accuracy )

"IHS5" evaluates I = §_Sn f(x₁,x₂,.....,x_n) dx₁ dx₂ ...... dx_n where (S_n): x₁² + x₂² + ........ + x_n² <= 1

Data Registers: • R00 = n > 1 ( Register R00 is to be initialized before executing "IHS5" )

                                         R01 = x₁   R02 = x₂   .....................   Rnn = x_n
Flags: /
Subroutine: A program LBL "T" that takes x₁ , x₂ , ............ , x_n in registers R01 R02 ...... Rnn   respectively
                       and returns f(x₁,x₂,.....,x_n) in X-register.

01 LBL "IHS5"
02 RCL 00
03 E3
04 /
05 ISG X
06 CLRGX
07 3
08 RCL 00
09 4
10 +
11 /
12 SQRT
13 STO 01
14 STO 02
15 CLA
16 X<>Y
17 STO N
18 .001
19 -
20 STO M
21 ISG N
22 LBL 01

23 XEQ "T"
24 ST+ O
25 RCL IND N
26 CHS
27 STO IND N
28 XEQ "T"
29 ST+ O
30 CLX
31 X<> IND N
32 CHS
33 ISG N
34 X<0?
35 GTO 00
36 STO IND N
37 GTO 01
38 LBL 00
39 RCL IND M
40 CHS
41 X>0?
42 GTO 00
43 STO IND M
44 RCL M

45 INT
46 1
47 +
48 RCL N
49 FRC
50 +
51 STO N
52 X<>Y
53 CHS
54 STO IND Y
55 GTO 01
56 LBL 00
57 CLX
58 X<> IND M
59 CHS
60 ISG M
61 X<0?
62 GTO 00
63 STO IND M
64 RCL M
65 INT
66 1

67 +
68 RCL N
69 FRC
70 +
71 STO N
72 X<>Y
73 STO IND Y
74 GTO 01
75 LBL 00
76 RCL 00
77 STO M
78 4
79 +
80 2
81 /
82 ST* O
83 CLX
84 STO N
85 X<>Y
86 LBL 03
87 STO IND M
88 XEQ "T"

89 ST+ N
90 RCL IND M
91 CHS
92 STO IND M
93 XEQ "T"
94 ST+ N
95 CLX
96 X<> IND M
97 CHS
98 DSE M
99 GTO 03
100 16
101 RCL 00
102 X^2
103 -
104 RCL N
105 *
106 ST+ O
107 XEQ "T"
108 RCL 00
109 3
110 -

111 RCL 00
112 *
113 10
114 -
115 RCL 00
116 *
117 36
118 +
119 *
120 ST+ O
121 RCL 00
122 ENTER
123 STO Z
124 2
125 ST/ Z
126 MOD
127 1
128 +
129 X<>Y
130 INT
131 STO Z
132 LBL 04

133 CLX
134 PI
135 ST* Y
136 X<> L
137 ST/ Y
138 SIGN
139 ST- L
140 DSE Z
141 GTO 04
142 CLX
143 RCL O
144 *
145 RCL 00
146 2
147 +
148 18
149 *
150 /
151 CLA
152 END

( 234 bytes / SIZE 006 )

STACK	INPUT	OUTPUT
X	/	I

where I = §_Sn f(x₁,x₂,.....,x_n) dx₁ dx₂ ...... dx_n with (S_n): x₁² + x₂² + ........ + x_n² <= 1

Example1: I = §§§§_S4 ( x² + y² + z² + u² ) dx dy dz du

01 LBL "T"
02 RCL 01
03 X^2
04 RCL 02
05 X^2
06 +
07 RCL 03
08 X^2
09 +
10 RCL 04
11 X^2
12 +
13 RTN

4 STO 00 XEQ "IHS5" >>>> I = 3.289868134 ---Execution time = 27s---

-The exact result is PI^2 / 3 , so all the decimals are correct ( rounded to 9 D )
-The formula gives perfect accuracy if f is a polynomial of degree < 6

Example2: I = §§§§§§_S6 Ln ( PI² + x + y² + z³ + u⁴ + v⁵ + w⁶ ) dx dy dz du dv dw

01 LBL "T"
02 RCL 01
03 RCL 02
04 X^2
05 +
06 RCL 03
07 3
08 Y^X
09 +
10 RCL 04
11 X^2
12 X^2
13 +
14 RCL 05
15 5
16 Y^X
17 +
18 RCL 06
19 6
20 Y^X
21 +
22 PI
23 X^2
24 +
25 LN
26 RTN

6 STO 00 XEQ "IHS5" >>>> I = 11.9174 ---Execution time = 2mn15s---

Notes:

-The formula is exact if f is a polynomial of degree < 6

-The subroutine must begin by LBL "T" otherwise change lines 23-28-88-93-107
-The function f is evaluated ( 2 n² + 1 ) times

-Registers R01 ..... Rnn are cleared by this program

-Lines 121 to 142 may be replaced by XEQ "HS" ( cf "Hyperspheres for the HP-41' )
-They return 0 if n > 187, so don't use this program in this case.
-They actually calculate PI^n/2/(n/2)! so the program may be simplified if you have an M-Code Gamma function.

b) Program#2 ( 7th degree of accuracy )

"IHS7" also evaluates I = §_Sn f(x₁,x₂,.....,x_n) dx₁ dx₂ ...... dx_n where (S_n): x₁² + x₂² + ........ + x_n² <= 1 but with a better precision

Data Registers: • R00 = n > 2 ( Register R00 is to be initialized before executing "IHS7" )

                                         R01 = x₁   R02 = x₂   .....................   Rnn = x_n   Rn+1: temp
Flags: /
Subroutine: A program LBL "T" that takes x₁ , x₂ , ............ , x_n in registers R01 R02 ...... Rnn   respectively
                       and returns f(x₁,x₂,.....,x_n) in X-register.

01 LBL "IHS7"
02 RCL 00
03 1
04 +
05 0
06 STO IND Y
07 RCL 00
08 E3
09 /
10 ISG X
11 CLRGX
12 3
13 RCL 00
14 6
15 +
16 /
17 SQRT
18 STO 01
19 STO 02
20 STO 03
21 CLA
22 CLX
23 2
24 +
25 STO O
26 1.001
27 -
28 STO N
29 LASTX
30 -
31 STO M
32 LBL 01
33 XEQ "T"
34 RCL 00
35 1
36 +
37 X<>Y
38 ST+ IND Y
39 RCL IND M
40 CHS
41 STO IND M
42 X<0?
43 GTO 01
44 RCL IND N
45 CHS
46 STO IND N
47 X<0?
48 GTO 01
49 RCL IND O
50 CHS

51 STO IND O
52 X<0?
53 GTO 01
54 CLX
55 X<> IND O
56 ISG O
57 X=0?
58 GTO 00
59 STO IND O
60 GTO 01
61 LBL 00
62 CLX
63 X<> IND N
64 ISG N
65 X=0?
66 GTO 00
67 LBL 02
68 STO IND N
69 RCL N
70 INT
71 RCL O
72 FRC
73 +
74 ISG X
75 STO O
76 X<>Y
77 STO IND Y
78 GTO 01
79 LBL 00
80 CLX
81 X<> IND M
82 ISG M
83 X=0?
84 GTO 00
85 STO IND M
86 RCL M
87 INT
88 RCL N
89 FRC
90 +
91 ISG X
92 STO N
93 X<>Y
94 GTO 02
95 LBL 00
96 STO 01
97 STO 02
98 CLX
99 X<> O
100 FRC

101 2
102 +
103 STO N
104 1.001
105 -
106 STO M
107 LBL 03
108 XEQ "T"
109 ST+ O
110 RCL IND M
111 CHS
112 STO IND M
113 X<0?
114 GTO 03
115 RCL IND N
116 CHS
117 STO IND N
118 X<0?
119 GTO 03
120 CLX
121 X<> IND N
122 ISG N
123 X=0?
124 GTO 00
125 STO IND N
126 GTO 03
127 LBL 00
128 CLX
129 X<> IND M
130 ISG M
131 X=0?
132 GTO 00
133 STO IND M
134 RCL M
135 INT
136 RCL N
137 FRC
138 +
139 ISG X
140 STO N
141 X<>Y
142 STO IND Y
143 GTO 03
144 LBL 00
145 5
146 RCL 00
147 -
148 4
149 /
150 RCL O

151 *
152 RCL 00
153 1
154 +
155 RDN
156 RCL IND T
157 8
158 /
159 +
160 STO O
161 X<>Y
162 STO 01
163 CLX
164 X<> N
165 FRC
166 ISG X
167 STO M
168 LBL 04
169 XEQ "T"
170 ST+ N
171 RCL IND M
172 CHS
173 STO IND M
174 X<0?
175 GTO 04
176 CLX
177 X<> IND M
178 ISG M
179 X=0?
180 GTO 00
181 STO IND M
182 GTO 04
183 LBL 00
184 CLX
185 X<> N
186 6
187 RCL 00
188 -
189 1
190 STO 01
191 LASTX
192 X^2
193 -
194 *
195 36
196 +
197 *
198 4
199 /
200 RCL 00

201 ST- M
202 3
203 +
204 /
205 ST+ O
206 LBL 05
207 XEQ "T"
208 ST+ N
209 RCL IND M
210 CHS
211 STO IND M
212 X<0?
213 GTO 05
214 CLX
215 X<> IND M
216 ISG M
217 X=0?
218 GTO 00
219 STO IND M
220 GTO 05
221 LBL 00
222 RCL N
223 81
224 *
225 RCL 00
226 3
227 +
228 RCL 00
229 6
230 +
231 X^2
232 *
233 /
234 ST+ O
235 XEQ "T"
236 RCL 00
237 45
238 *
239 324
240 +
241 RCL 00
242 *
243 216
244 +
245 RCL 00
246 6
247 +
248 X^2
249 STO M
250 /

251 RCL 00
252 6
253 -
254 X^2
255 29
256 +
257 RCL 00
258 *
259 6
260 /
261 -
262 *
263 ST+ O
264 RCL 00
265 RCL 00
266 RCL 00
267 2
268 ST/ Z
269 MOD
270 1
271 +
272 X<>Y
273 INT
274 STO Z
275 LBL 06
276 CLX
277 PI
278 ST* Y
279 X<> L
280 ST/ Y
281 SIGN
282 ST- L
283 DSE Z
284 GTO 06
285 X<> M
286 *
287 RCL 00
288 2
289 +
290 RCL 00
291 4
292 +
293 *
294 27
295 *
296 /
297 RCL O
298 *
299 CLA
300 END

( 443 bytes / SIZE n+2 )

STACK	INPUT	OUTPUT
X	/	I

where I = §_Sn f(x₁,x₂,.....,x_n) dx₁ dx₂ ...... dx_n with (S_n): x₁² + x₂² + ........ + x_n² <= 1

Example1: I = §§§§§§_S6 Ln ( PI + x² + y² + z² + u² + v² + w² ) dx dy dz du dv dw

01 LBL "T"
02 RCL 01
03 X^2
04 RCL 02
05 X^2
06 +
07 RCL 03
08 X^2
09 +
10 RCL 04
11 X^2
12 +
13 RCL 05
14 X^2
15 +
16 RCL 06
17 X^2
18 +
19 PI
20 +
21 LN
22 RTN

6 STO 00 XEQ "IHS7" >>>> I = 7.015497950 ---Execution time = 5m47s---

-The exact result is I = 7.015376313...

Example2: I = §§§§§§_S6 Ln ( PI² + x + y² + z³ + u⁴ + v⁵ + w⁶ ) dx dy dz du dv dw

01 LBL "T"
02 RCL 01
03 RCL 02
04 X^2
05 +
06 RCL 03
07 3
08 Y^X
09 +
10 RCL 04
11 X^2
12 X^2
13 +
14 RCL 05
15 5
16 Y^X
17 +
18 RCL 06
19 6
20 Y^X
21 +
22 PI
23 X^2
24 +
25 LN
26 RTN

6 STO 00 XEQ "IHS7" >>>> I = 11.91901135 ---Execution time = 7mn47s---

Notes:

-The formula is exact if f is a polynomial of degree < 8
-The subroutine must begin by LBL "T"
-Registers R01 ..... Rnn are cleared by this program

>>> n must be at least 3

-The function f is evaluated ( 4 n³ - 6 n² + 14 n + 3 ) / 3 times

n	3	4	5	6	7	8	9	10
eval	33	73	141	245	393	593	853	1181

4°) Integrals on the Surface of the Unit HyperSphere

a) Program#1 ( 5th degree of accuracy )

"ISHS5" evaluates I = §_dSn f(x₁,x₂,.....,x_n) ds where (dS_n): x₁² + x₂² + ........ + x_n² = 1

-After solving the equations (1) to (4) given in reference [3], we find the following weights: w

for the 2n points ( x , 0 , ........... , 0 ) x² = 1 w = A ( 4 - n ) / ( 2 n² + 4 n )
for the 2n(2n-1) points ( y , y , 0 , ..... , 0 ) y² = 1/2 w = A / ( n² + 2 n )

where A = Surface of the hypersphere

Data Registers: • R00 = n > 1 ( Register R00 is to be initialized before executing "ISHS5" )

01 LBL "ISHS5"
02 RCL 00
03 E3
04 /
05 ISG X
06 CLRGX
07 .5
08 SQRT
09 STO 01
10 STO 02
11 CLA
12 SIGN
13 +
14 STO N
15 1.001
16 -
17 STO M
18 LBL 01
19 XEQ "T"
20 ST+ O
21 RCL IND M
22 CHS
23 STO IND M

24 X<0?
25 GTO 01
26 RCL IND N
27 CHS
28 STO IND N
29 X<0?
30 GTO 01
31 CLX
32 X<> IND N
33 ISG N
34 X=0?
35 GTO 00
36 STO IND N
37 GTO 01
38 LBL 00
39 CLX
40 X<> IND M
41 ISG M
42 X=0?
43 GTO 00
44 STO IND M
45 RCL M
46 INT

47 RCL N
48 FRC
49 +
50 ISG X
51 STO N
52 X<>Y
53 STO IND Y
54 GTO 01
55 LBL 00
56 1
57 STO 01
58 CLX
59 X<> N
60 FRC
61 ISG X
62 STO M
63 LBL 02
64 XEQ "T"
65 ST+ N
66 RCL IND M
67 CHS
68 STO IND M
69 X<0?

70 GTO 02
71 CLX
72 X<> IND M
73 ISG M
74 X=0?
75 GTO 00
76 STO IND M
77 GTO 02
78 LBL 00
79 X<> N
80 4
81 RCL 00
82 -
83 *
84 2
85 /
86 ST+ O
87 RCL 00
88 RCL 00
89 RCL 00
90 2
91 ST/ Z
92 MOD

93 1
94 +
95 X<>Y
96 INT
97 STO Z
98 LBL 03
99 CLX
100 PI
101 ST* Y
102 X<> L
103 ST/ Y
104 SIGN
105 ST- L
106 DSE Z
107 GTO 03
108 X<> O
109 *
110 2
111 RCL 00
112 +
113 /
114 CLA
115 END

( 184 bytes / SIZE n+1 )

STACK	INPUT	OUTPUT
X	/	I

where I = §_dSn f(x₁,x₂,.....,x_n) ds with (dS_n): x₁² + x₂² + ........ + x_n² = 1

Example: Calculate I = §§§_dS4 Ln ( Pi + x² y + z t ) ds

01 LBL "T"
02 RCL 01
03 X^2
04 RCL 02
05 *
06 RCL 03
07 RCL 04
08 *
09 +
10 PI
11 +
12 LN
13 RTN

4 STO 00

XEQ "ISHS5" >>>> I = 22.53289243 ---Execution time = 35s---

Notes:

-The formula is exact if f is a polynomial of degree < 6
-The subroutine must begin by LBL "T"
-Registers R01 thru Rnn are cleared.
-There are 2 n² evaluations of the function

b) Program#2 ( 7th degree of accuracy )

"ISHS7" also evaluates I = §_dSn f(x₁,x₂,.....,x_n) ds where (dS_n): x₁² + x₂² + ........ + x_n² = 1 but with a better precision.

-It generalizes the program listed in paragraph 2°)a)

-After solving the equations (1) to (7) given in reference [3], we find the following weights: w

for the 2n points                     ( x , 0 , ............... , 0 )       x² = 1        w = A ( n² - 9 n + 38 ) / 4 / ( n³ + 6 n² + 8 n )
for the 2n(2n-1) points            ( y , y , 0 , ........... , 0 )       y² = 1/2     w = A ( - 2 n + 10 ) / ( n³ + 6 n² + 8 n )
for the 4n(n-1)(n-2)/3 points ( z , z , z , 0 , ....... , 0 )       z² = 1/3      w = 27 A / ( n³ + 6 n² + 8 n )

where A = Surface of the hypersphere

Data Registers: • R00 = n > 2 ( Register R00 is to be initialized before executing "ISHS7" )

01 LBL "ISHS7"
02 RCL 00
03 1
04 +
05 0
06 STO IND Y
07 RCL 00
08 E3
09 /
10 ISG X
11 CLRGX
12 3
13 1/X
14 SQRT
15 STO 01
16 STO 02
17 STO 03
18 CLA
19 CLX
20 2
21 +
22 STO O
23 1.001
24 -
25 STO N
26 LASTX
27 -
28 STO M
29 LBL 01
30 XEQ "T"
31 RCL 00
32 1

33 +
34 X<>Y
35 ST+ IND Y
36 RCL IND M
37 CHS
38 STO IND M
39 X<0?
40 GTO 01
41 RCL IND N
42 CHS
43 STO IND N
44 X<0?
45 GTO 01
46 RCL IND O
47 CHS
48 STO IND O
49 X<0?
50 GTO 01
51 CLX
52 X<> IND O
53 ISG O
54 X=0?
55 GTO 00
56 STO IND O
57 GTO 01
58 LBL 00
59 CLX
60 X<> IND N
61 ISG N
62 X=0?
63 GTO 00
64 LBL 02

65 STO IND N
66 RCL N
67 INT
68 RCL O
69 FRC
70 +
71 ISG X
72 STO O
73 X<>Y
74 STO IND Y
75 GTO 01
76 LBL 00
77 CLX
78 X<> IND M
79 ISG M
80 X=0?
81 GTO 00
82 STO IND M
83 RCL M
84 INT
85 RCL N
86 FRC
87 +
88 ISG X
89 STO N
90 X<>Y
91 GTO 02
92 LBL 00
93 2
94 1/X
95 SQRT
96 STO 01

97 STO 02
98 CLX
99 X<> O
100 FRC
101 2
102 +
103 STO N
104 1.001
105 -
106 STO M
107 LBL 03
108 XEQ "T"
109 ST+ O
110 RCL IND M
111 CHS
112 STO IND M
113 X<0?
114 GTO 03
115 RCL IND N
116 CHS
117 STO IND N
118 X<0?
119 GTO 03
120 CLX
121 X<> IND N
122 ISG N
123 X=0?
124 GTO 00
125 STO IND N
126 GTO 03
127 LBL 00
128 CLX

129 X<> IND M
130 ISG M
131 X=0?
132 GTO 00
133 STO IND M
134 RCL M
135 INT
136 RCL N
137 FRC
138 +
139 ISG X
140 STO N
141 X<>Y
142 STO IND Y
143 GTO 03
144 LBL 00
145 40
146 RCL 00
147 8
148 *
149 -
150 ST* O
151 RCL 00
152 1
153 +
154 RCL IND X
155 13.5
156 *
157 ST+ O
158 1
159 STO 01
160 CLX

161 X<> N
162 FRC
163 ISG X
164 STO M
165 LBL 04
166 XEQ "T"
167 ST+ N
168 RCL IND M
169 CHS
170 STO IND M
171 X<0?
172 GTO 04
173 CLX
174 X<> IND M
175 ISG M
176 X=0?
177 GTO 00
178 STO IND M
179 GTO 04
180 LBL 00
181 X<> N
182 RCL 00
183 9
184 -
185 RCL 00
186 *
187 38
188 +
189 *
190 ST+ O
191 RCL 00
192 RCL 00

193 RCL 00
194 2
195 ST/ Z
196 MOD
197 1
198 +
199 X<>Y
200 INT
201 STO Z
202 LBL 05
203 CLX
204 PI
205 ST* Y
206 X<> L
207 ST/ Y
208 SIGN
209 ST- L
210 DSE Z
211 GTO 05
212 X<> O
213 *
214 2
215 RCL 00
216 ST+ Y
217 4
218 +
219 *
220 4
221 *
222 /
223 CLA
224 END

( 346 bytes / SIZE n+2 )

STACK	INPUT	OUTPUT
X	/	I

where I = §_dSn f(x₁,x₂,.....,x_n) ds with (dS_n): x₁² + x₂² + ........ + x_n² = 1

Example: Calculate I = §§§_dS4 Ln ( Pi + x² y + z t ) ds

01 LBL "T"
02 RCL 01
03 X^2
04 RCL 02
05 *
06 RCL 03
07 RCL 04
08 *
09 +
10 PI
11 +
12 LN
13 RTN

4 STO 00

XEQ "ISHS7" >>>> I = 22.53840629 ---Execution time = 80s---

Notes:

-The formula is exact if f is a polynomial of degree < 8
-The subroutine must begin by LBL "T"
-Registers R01 thru Rnn are cleared.
-There are ( 4 n³ - 6 n² + 8 n ) / 3 evaluations of the function

n	3	4	5	6	7	8	9	10
eval	26	64	130	232	378	576	834	1160

References:

[1] Abramowitz and Stegun - "Handbook of Mathematical Functions" - Dover Publications - ISBN 0-486-61272-4
[2] Christopher A. Feuchter - "Numerical Integration Over a Sphere"
[3] Preston C. Hammer and Arthur H. Stroud - "Numerical Evaluation of Multiple Integrals II"
[4] L. N. Dobrodeev - A cubature formula of the seventh degree of accuracy for a hypersphere and a hypercube.

hp41programs

Integrals over a Sphere & its Surface for the HP-41