Gaussian Integration for the HP-41
Overview
1°) Gauss-Legendre Formulae
a) 2 M-Code Routines: GX & GW
b) Gauss-Legendre 10-point Formula
2°) Gauss-Laguerre
a) 10-point Formula
b) 15-point Formula
3°) Gauss-Hermite
a) 20-point Formula
b) 30-point Formula
4°) Double Integrals
a) Gauss-Legendre Formulae
b) Gauss-Laguerre Formulae
c) Gauss-Hermite Formulae
5°) Triple Integrals
a) Gauss-Legendre Formulae
a-1) Triple Integrals only
a-2) Simple, Double & Triple
Integrals
b) Gauss-Laguerre Formulae
c) Gauss-Hermite Formulae
-This page completes "Numerical Integration" , "Double Integrals" and
"Triple Integrals".
1°) Gauss-Legendre Formulae
a) 2 M-Code Routines:
GX & GW
-These routines take an integer i between 1 and 5 in X-register and
replace it by the abscissa xi or the weight
wi of the 10-point Gauss-Legendre formula.
098 "X"
007 "G"
2A0 SETDEC
0F8 C=X
2FC RCR 13
27E C=C-1 MS
10E A=C ALL
04E C=0 ALL
35C PT=12
1BE A=A-1 MS
05B JNC+11d
050 LD@PT- 1
110 LD@PT- 4
210 LD@PT- 8
210 LD@PT- 8
1D0 LD@PT- 7
110 LD@PT- 4
0D0 LD@PT- 3
0D0 LD@PT- 3
250 LD@PT- 9
193 JNC+50d
1BE A=A-1 MS
063 JNC+12d
110 LD@PT- 4
0D0 LD@PT- 3
0D0 LD@PT- 3
0D0 LD@PT- 3
250 LD@PT- 9
150 LD@PT- 5
0D0 LD@PT- 3
250 LD@PT- 9
110 LD@PT- 4
050 LD@PT- 1
12B JNC+37d
1BE A=A-1 MS
063 JNC+12d
190 LD@PT- 6
1D0 LD@PT- 7
250 LD@PT- 9
110 LD@PT- 4
010 LD@PT- 0
250 LD@PT- 9
150 LD@PT- 5
190 LD@PT- 6
210 LD@PT- 8
0D0 LD@PT- 3
0C3 JNC+24d
1BE A=A-1 MS
063 JNC+12d
210 LD@PT- 8
190 LD@PT- 6
150 LD@PT- 5
010 LD@PT- 0
190 LD@PT- 6
0D0 LD@PT- 3
0D0 LD@PT- 3
190 LD@PT- 6
190 LD@PT- 6
1D0 LD@PT- 7
05B JNC+11d
250 LD@PT- 9
1D0 LD@PT- 7
0D0 LD@PT- 3
250 LD@PT- 9
010 LD@PT- 0
190 LD@PT- 6
150 LD@PT- 5
090 LD@PT- 2
210 LD@PT- 8
150 LD@PT- 5
266 C=C-1 S&X
0E8 X=C
3E0 RTN
097 "W"
007 "G"
2A0 SETDEC
0F8 C=X
2FC RCR 13
27E C=C-1 MS
10E A=C ALL
04E C=0 ALL
35C PT=12
266 C=C-1 S&X
1BE A=A-1 MS
063 JNC+12d
090 LD@PT- 2
250 LD@PT- 9
150 LD@PT- 5
150 LD@PT- 5
090 LD@PT- 2
110 LD@PT- 4
090 LD@PT- 2
090 LD@PT- 2
110 LD@PT- 4
1D0 LD@PT- 7
19B JNC+51d
1BE A=A-1 MS
063 JNC+12d
090 LD@PT- 2
190 LD@PT- 6
250 LD@PT- 9
090 LD@PT- 2
190 LD@PT- 6
190 LD@PT- 6
1D0 LD@PT- 7
050 LD@PT- 1
250 LD@PT- 9
0D0 LD@PT- 3
133 JNC+36d
1BE A=A-1 MS
063 JNC+12d
090 LD@PT- 2
050 LD@PT- 1
250 LD@PT- 9
010 LD@PT- 0
210 LD@PT- 8
190 LD@PT- 6
0D0 LD@PT- 3
190 LD@PT- 6
090 LD@PT- 2
150 LD@PT- 5
0CB JNC+25d
1BE A=A-1 MS
063 JNC+12d
050 LD@PT- 1
110 LD@PT- 4
250 LD@PT- 9
110 LD@PT- 4
150 LD@PT- 5
050 LD@PT- 1
0D0 LD@PT- 3
110 LD@PT- 4
250 LD@PT- 9
090 LD@PT- 2
063 JNC+12d
190 LD@PT- 6
190 LD@PT- 6
190 LD@PT- 6
1D0 LD@PT- 7
050 LD@PT- 1
0D0 LD@PT- 3
110 LD@PT- 4
110 LD@PT- 4
0D0 LD@PT- 3
050 LD@PT- 1
266 C=C-1 S&X
0E8 X=C
3E0 RTN
( 148 words )
"GX"
X = 1 -----> X = 0.1488743390
X = 2 -----> X = 0.4333953941
X = 3 -----> X = 0.6794095683
X = 4 -----> X = 0.8650633667
X = 5 -----> X = 0.9739065285
"GW"
X = 1 -----> X = 0.2955242247
X = 2 -----> X = 0.2692667193
X = 3 -----> X = 0.2190863625
X = 4 -----> X = 0.1494513492
X = 5 -----> X = 0.06667134431
-Only the 1st digit of the mantissa in X-input is taken into account.
-This digit is shifted in the MS field of CPU register A.
-So, similar routines may be written for Gauss-Legendre m-point formulae,
provided m is even and between 2 and 18.
-It would be hardly more difficult if m = 20 , 22 , ...
-Though there is no check for alpha data, there is no risk of crash.
-These routines are used in the focal programs "GL10" "DGL10"
"TGL10" & "SDTI"
Improvement:
-In these 4 programs, "GW" is always followed by a multiplication, so we may insert this operation in the M-code routine itself:
Replace the last 2 words of "GW" ( 0E8 X=C 3E0 RTN ) by
10E A=C ALL
0B8 C=Y
135 C=
060 A*C
0E8 X=C
3E0 RTN
and delete * in the focal programs.
-Likewise, "GX" is always used inside instructions of the type:
RCL aa
STO nn
RCL bb
RCL cc
GX
*
ST+ nn
-
They may be replaced by
RCL aa
RCL bb
RCL cc
GX
STO nn
X<>Y
if you replace the last 2 words of "GX" ( 0E8 X=C
3E0 RTN ) by
10E A=C ALL
0B8 C=Y
135 C=
060 A*C
070 N=C ALL
2BE C=-C
10E A=C ALL
078 C=Z
01D C=
060 A+C
0A8 Y=C
078 C=Z
10E A=C
0B0 C=N ALL
01D C=
060 A+C
0E8 X=C
3E0 RTN
-"GX" will then change the stack
Z = z
Y = y
X = i into
Y = z - y wi
X = z + y wi
-These improvements. will save both bytes & time:
-Especially with triple integrals, the execution time will be reduced
by 30 seconds.
b) Gauss-Legendre 10-point
Formula
-Gauss-Legendre formulas are defined by §-11
f(x).dx ~ w1.f(x1) + w2.f(x2)
+ ....... + wm.f(xm)
where wi and xi are choosen so that
the formula produces perfect accuracy when f(x) is a polynomial of
degree < 2m
-A linear change of variable maps an interval [a,b] to the standard
[-1,+1]
-The interval [a,b] may also be divided into n subintervals for a better
accuracy:
-Theoretically, the approximations tend to the exact result as n tends
to infinity.
-"GL10" evaluates §ab f(x) dx
with the 10-point formula.
Data Registers: • R00 = function name ( Registers R00 thru R03 are to be initialized before executing "GL10" )
• R01 = a
R04 = Integral
• R02 = b
• R03 = n
R05 thru R09: temp
Flags: /
Subroutine: A program that takes x in X-register
and returns f(x) in X-register
-If you don't want to use the M-Code routines GX & GW, replace:
lines 31-32 by DSE 08 RCL
IND 08
lines 21-22 by RCL IND 08
line 15 by 19.009
and store the following constants into R10 to R19 before executing "GL10"
R10 = 0.2955242247
R11 = 0.1488743390
R12 = 0.2692667193
R13 = 0.4333953941
R14 = 0.2190863625
R15 = 0.6794095683
R16 = 0.1494513492
R17 = 0.8650633667
R18 = 0.06667134431
R19 = 0.9739065285
| 01 LBL "GL10"
02 CLX 03 STO 04 04 RCL 02 05 RCL 01 06 STO 05 07 - 08 RCL 03 09 STO 06 |
10 ST+ X
11 / 12 STO 07 13 LBL 01 14 ST+ 05 15 5 16 STO 08 17 LBL 02 18 RCL 05 |
19 STO 09
20 RCL 07 21 RCL 08 22 GX 23 * 24 ST+ 09 25 - 26 XEQ IND 00 27 X<> 09 |
28 XEQ IND 00
29 RCL 09 30 + 31 RCL 08 32 GW 33 * 34 ST+ 04 35 DSE 08 36 GTO 02 |
37 RCL 07
38 ST+ 05 39 DSE 06 40 GTO 01 41 ST* 04 42 RCL 04 43 END |
( 67 bytes / SIZE 010 )
| STACK | INPUT | OUTPUT |
| X | / | I |
where I = §ab f(x) dx
Example: Evaluate §04 Sin x2 dx
LBL "T" any global
LBL maximum 6 characters
X^2
SIN
RTN
XEQ "RAD" "T" ASTO 00 0 STO 01 4 STO 02
-With n = 1 1 STO 03 XEQ "GL10" >>>> I = 0.748650151
-With n = 2 2 STO 03 R/S >>>> I = 0.747133892
-With n = 3 3 STO 03
R/S >>>>
I = 0.747133845 ( correct to
9D )
---Execution time = 29s----
Notes:
-If you've written GX & GW for, say, Gauss-Legendre 18-point formula, simply replace 5 ( line 15 ) by 9
-If you prefer the improved versions of "GX" & "GW" as suggested in the note paragraph 1°)a)
Delete line 33 and replace lines 18 to 25 by
RCL 05
RCL 07
RCL 08
GX
STO 09
X<>Y
2°) Gauss-Laguerre
a) 10-Point Formula
-If the upper or the lower limit of integration is infinite, one may
use a change of variable like x = tan u.
-But special types of Gaussian formulae have also been derived for
these integrals.
-Gauss-Laguerre m-point formulas are of the type:
§0+¥ exp(-x)
f(x).dx ~ w1.f(x1) + w2.f(x2)
+ ....... + wm.f(xm)
where wi and xi are choosen so that
the formula produces perfect accuracy when f(x) is a polynomial of
degree < 2m
-"GLA10" employs the 10-point formula.
Data Registers: • R00 = function name ( Registers R00 and R03 thru R22 are to be initialized before executing "GLA10" )
R01 = Integral R02: temp
• R03 = 3.084411158 E-01
• R04 = 1.377934705 E-01
• R05 = 4.011199292 E-01
• R06 = 7.294545495 E-01
• R07 = 2.180682876 E-01
• R08 = 1.808342902
• R09 = 6.208745610 E-02
• R10 = 3.401433698
• R11 = 9.501516975 E-03
• R12 = 5.552496140
• R13 = 7.530083886 E-04
• R14 = 8.330152747
• R15 = 2.825923350 E-05
• R16 = 11.84378584
• R17 = 4.249313985 E-07
• R18 = 16.27925783
• R19 = 1.839564824 E-09
• R20 = 21.99658581
• R21 = 9.911827220 E-13
• R22 = 29.92069701
Flags: /
Subroutine: A program that takes x in X-register
and returns f(x) in X-register
| 01 LBL "GLA10"
02 CLX 03 STO 01 04 22.002 05 STO 02 06 LBL 01 07 RCL IND 02 08 XEQ IND 00 09 DSE 02 10 RCL IND 02 11 * 12 ST+ 01 13 DSE 02 14 GTO 01 15 RCL 01 16 END |
( 38 bytes / SIZE 023 )
| STACK | INPUT | OUTPUT |
| X | / | I |
where I = §0+¥ exp(-x) f(x) dx
Example: Evaluate I = §0+¥ exp(-x) Ln ( 1 + x ) dx
LBL "T"
LN1+X
RTN
"T" ASTO 00 XEQ "GLA10" >>>> I ~ 0.596347721 ---Execution time = 10s---
Note:
-If you have to calculate §a+¥
exp(-x) f(x) dx , use the change of variable x = a + u
b) 15-Point Formula
-Just change a few lines and the coefficients to use the Gauss-Laguerre
15-point formula
Data Registers: • R00 = function name ( Registers R00 and R03 thru R32 are to be initialized before executing "GLA15" )
R01 = Integral R02: temp
• R03 = 2.182348859 E-01
• R04 = 0.09330781202
• R05 = 3.422101779 E-01
• R06 = 0.4926917403
• R07 = 2.630275779 E-01
• R08 = 1.215595412
• R09 = 1.264258181 E-01
• R10 = 2.269949526
• R11 = 4.020686492 E-02
• R12 = 3.667622722
• R13 = 8.563877804 E-03
• R14 = 5.425336627
• R15 = 1.212436147 E-03
• R16 = 7.565916227
• R17 = 1.116743923 E-04
• R18 = 10.12022857
• R19 = 6.459926762 E-06
• R20 = 13.13028248
• R21 = 2.226316907 E-07
• R22 = 16.65440771
• R23 = 4.227430385 E-09
• R24 = 20.77647890
• R25 = 3.921897267 E-11
• R26 = 25.62389423
• R27 = 1.456515264 E-13
• R28 = 31.40751917
• R29 = 1.483027051 E-16
• R30 = 38.53068331
• R31 = 1.600594906 E-20
• R32 = 48.02608557
Flags: /
Subroutine: A program that takes x in X-register
and returns f(x) in X-register
| 01 LBL "GLA15"
02 CLX 03 STO 01 04 32.002 05 STO 02 06 LBL 01 07 RCL IND 02 08 XEQ IND 00 09 DSE 02 10 RCL IND 02 11 * 12 ST+ 01 13 DSE 02 14 GTO 01 15 RCL 01 16 END |
( 38 bytes / SIZE 023 )
| STACK | INPUT | OUTPUT |
| X | / | I |
where I = §0+¥ exp(-x) f(x) dx
Example: Evaluate I = §0+¥ exp(-x) Ln ( 1 + x ) dx
LBL "T"
LN1+X
RTN
"T" ASTO 00 XEQ "GLA15" >>>>
I ~ 0.596347721
---Execution time = 12s---
Notes:
-Is it the best evaluation of this integral that we can get from these
formulae ?
-Let's try the change of variable x = u/2.1 the subroutine
becomes:
LBL "T"
2.1
/
ENTER^
LN1+X
X<>Y
1.1
*
E^X
*
2.1
/
RTN
XEQ "GLA15" >>>> I ~ 0.5963473625 whereas "GLA10" gives 0.596347379
-The exact result - rounded to 10D - is 0.5963473623
-This trick will work sometimes ... but not always !
-Similar routines may be written for any Gauss-Laguerre formula: simply
modify line 04 & the coefficients.
3°) Gauss-Hermite
a) 20-Point Formula
-Gauss-Hermite m-point formulas are: §-¥+¥
exp(-x2) f(x).dx ~ w1.f(x1)
+ w2.f(x2) + ....... + wm.f(xm)
where wi and xi are choosen so that
the formula produces perfect accuracy when f(x) is a polynomial of
degree < 2m
-The following programs assume that m = 2 p is even so that the formula may be re-written:
§-¥+¥ exp(-x2) f(x).dx ~ w1.[ f(x1) + f(-x1) ] + w2.[ f(x2) + f(-x2) ] + ....... + wp.[ f(xp) + f(-xp) ]
-"GH20" employs the 20-point formula.
Data Registers: • R00 = function name ( Registers R00 and R04 thru R23 are to be initialized before executing "GH20" )
R01 = Integral R02-R03: temp
• R04 = 4.622436696 E-01
• R05 = 0.2453407083
• R06 = 2.866755054 E-01
• R07 = 0.7374737285
• R08 = 1.090172060 E-01
• R09 = 1.234076215
• R10 = 2.481052089 E-02
• R11 = 1.738537712
• R12 = 3.243773342 E-03
• R13 = 2.254974002
• R14 = 2.283386360 E-04
• R15 = 2.788806058
• R16 = 7.802556479 E-06
• R17 = 3.347854567
• R18 = 1.086069371 E-07
• R19 = 3.944764040
• R20 = 4.399340992 E-10
• R21 = 4.603682450
• R22 = 2.229393646 E-13
• R23 = 5.387480890
Flags: /
Subroutine: A program that takes x in X-register
and returns f(x) in X-register
| 01 LBL "GH20"
02 23.003 03 STO 03 04 CLX 05 STO 01 06 LBL 01 07 RCL IND 03 08 STO 02 09 XEQ IND 00 10 X<> 02 11 CHS 12 XEQ IND 00 13 RCL 02 14 + 15 DSE 03 16 RCL IND 03 17 * 18 ST+ 01 19 DSE 03 20 GTO 01 21 RCL 01 22 END |
( 45 bytes / SIZE 024 )
| STACK | INPUT | OUTPUT |
| X | / | I |
where I = §-¥+¥ exp(-x2) f(x).dx
Example: Evaluate I = §-¥+¥ exp(-x2) Ln ( 1 + x + x2 ) dx
LBL "T"
ENTER^
X^2
+
LN1+X
RTN
"T" ASTO 00 XEQ "GH20" >>>>
I ~ 0.451490093
---Execution time = 15s---
b) 30-Point Formula
-Lines 02 and 01 are the unique differences between "GH20" & "GH30"
Data Registers: • R00 = function name ( Registers R00 and R04 thru R33 are to be initialized before executing "GH30" )
R01 = Integral R02-R03: temp
• R04 = 3.863948895 E-01
• R05 = 0.2011285765
• R06 = 2.801309308 E-01
• R07 = 0.6039210586
• R08 = 1.467358475 E-01
• R09 = 1.008338271
• R10 = 5.514417687 E-02
• R11 = 1.415527800
• R12 = 1.470382970 E-02
• R13 = 1.826741144
• R14 = 2.737922473 E-03
• R15 = 2.243391468
• R16 = 3.483101243 E-04
• R17 = 2.667132125
• R18 = 2.938725229 E-05
• R19 = 3.099970530
• R20 = 1.579094887 E-06
• R21 = 3.544443873
• R22 = 5.108522451 E-08
• R23 = 4.003908604
• R24 = 9.178580424 E-10
• R25 = 4.483055357
• R26 = 8.106186297 E-12
• R27 = 4.988918969
• R28 = 2.878607081 E-14
• R29 = 5.533147152
• R30 = 2.810333603 E-17
• R31 = 6.138279220
• R32 = 2.908254700 E-21
• R33 = 6.863345294
Flags: /
Subroutine: A program that takes x in X-register
and returns f(x) in X-register
| 01 LBL "GH30"
02 33.003 03 STO 03 04 CLX 05 STO 01 06 LBL 01 07 RCL IND 03 08 STO 02 09 XEQ IND 00 10 X<> 02 11 CHS 12 XEQ IND 00 13 RCL 02 14 + 15 DSE 03 16 RCL IND 03 17 * 18 ST+ 01 19 DSE 03 20 GTO 01 21 RCL 01 22 END |
( 45 bytes / SIZE 024 )
| STACK | INPUT | OUTPUT |
| X | / | I |
where I = §-¥+¥ exp(-x2) f(x).dx
Example: Evaluate I = §-¥+¥ exp(-x2) Ln ( 1 + x + x2 ) dx
LBL "T"
ENTER^
X^2
+
LN1+X
RTN
"T" ASTO 00 XEQ "GH30" >>>> I ~ 0.451471191 ---Execution time = 22s---
Notes:
-With the change of argument x = u / 1.5 , "GH20" gives 0.451469515 and "GH30" returns 0.451469593 which is correct to 9D.
-Similar routines may be written for any Gauss-Hermite formula, provided
the number of points is even: simply modify line 02 & the coefficients.
-If the number of points is odd, f(0) must be computed without adding
f(-0) before multiplying by the corresponding weight.
-So, the program would be slightly longer.
4°) Double Integrals
a) Gauss-Legendre Formulae
-"DGL10" uses the Gauss-Legendre 10-point formula to evaluate:
§ab
§u(x)v(x)
f(x;y) dx dy
Data Registers: • R00 = f name ( Registers R00 thru R05 are to be initialized before executing "DGL10" )
• R01 = a
• R04 = u name
R07 thru R18: temp
• R02 = b
• R05 = v name
• R03 = n
R06 = Integral
Flags: /
Subroutines:
-A program that takes x in X-register & y in Y-register and returns
f(x,y) in X-register
-A program that takes x in X-register and returns u(x)
in X-register
-A program that takes x in X-register and returns v(x)
in X-register
-If you don't want to use the M-Code routines GX & GW, replace:
lines 78-79 by DSE 07 RCL IND 07 & lines
31-32 by DSE 11 RCL IND 11
lines 66-67 by RCL IND 07 & lines 21-22
by RCL IND 11
lines 15 & 60 by 28.018
and store the following constants into R19 to R28:
R19 = 0.2955242247
R20 = 0.1488743390
R21 = 0.2692667193
R22 = 0.4333953941
R23 = 0.2190863625
R24 = 0.6794095683
R25 = 0.1494513492
R26 = 0.8650633667
R27 = 0.06667134431
R28 = 0.9739065285
| 01 LBL "DGL10"
02 CLX 03 STO 06 04 RCL 02 05 RCL 01 06 STO 08 07 - 08 RCL 03 09 STO 18 10 ST+ X 11 / 12 STO 10 13 LBL 01 14 ST+ 08 15 5 16 STO 11 17 LBL 02 18 RCL 08 19 STO 12 |
20 RCL 10
21 RCL 11 22 GX 23 * 24 ST+ 12 25 - 26 XEQ 03 27 X<> 12 28 XEQ 03 29 RCL 12 30 + 31 RCL 11 32 GW 33 * 34 ST+ 06 35 DSE 11 36 GTO 02 37 RCL 10 38 ST+ 08 |
39 DSE 18
40 GTO 01 41 ST* 06 42 GTO 06 43 LBL 03 44 STO 13 45 XEQ IND 04 46 STO 14 47 CLX 48 STO 17 49 RCL 13 50 XEQ IND 05 51 RCL 14 52 - 53 RCL 03 54 STO 16 55 ST+ X 56 / 57 STO 15 |
58 LBL 04
59 ST+ 14 60 5 61 STO 07 62 LBL 05 63 RCL 14 64 STO 09 65 RCL 15 66 RCL 07 67 GX 68 * 69 ST+ 09 70 - 71 RCL 13 72 XEQ IND 00 73 X<> 09 74 RCL 13 75 XEQ IND 00 76 RCL 09 |
77 +
78 RCL 07 79 GW 80 * 81 ST+ 17 82 DSE 07 83 GTO 05 84 RCL 15 85 ST+ 14 86 DSE 16 87 GTO 04 88 RCL 17 89 * 90 RTN 91 LBL 06 92 RCL 06 93 END |
( 141 bytes / SIZE 019 )
| STACK | INPUTS | OUTPUTS |
| X | / | I |
where I = §ab §u(x)v(x) f(x;y) dx dy
Example: Evaluate I = §12
§xx^2
(1+x4 y4)1/2 dx dy.
01 LBL "T"
( any global labels, maximum of 6 characters )
02 *
03 X^2
04 X^2
05 SIGN
06 LAST X
07 +
08 SQRT
09 RTN
10 LBL "U"
11 RTN
12 LBL "V"
13 X^2
14 RTN
"T" ASTO 00
1 STO 01
2 STO 02
"U" ASTO 04
"V" ASTO 05
n = 1 1 STO 03
XEQ "DGL10" >>>> 15.46686245
in X-register & R06
---Execution time = 1m33s---
n = 2 2 STO 03
R/S >>>>
15.46686247 in X-register & R06
---Execution time = 5m44s---
Notes:
-The exact result ( rounded to 9 decimals ) is already obtained with
n = 1
-With n = 2 , roundoff-errors are slightly larger...
-If you've written GX & GW for, say, Gauss-Legendre 18-point formula, simply replace 5 ( lines 15 & 60 ) by 9
-If you prefer the improved versions of "GX" & "GW" as suggested in the note paragraph 1°)a)
Delete line 80 & replace lines 63 to 70 by
RCL 14
RCL 15
RCL 07
GX
STO 09
X<>Y
and delete line 33 & replace lines 18 to 25 by
RCL 08
RCL 10
RCL 11
GX
STO 12
X<>Y
b) Gauss-Laguerre Formulae
-The Gauss-Laguerre 10-point & 15-point formulae may be used to
compute I = §0+¥
§0+¥ exp(-x-y)
f(x,y) dx dy
Data Registers: • R00 = function name ( Registers R00 and R05 thru R24 are to be initialized before executing "DGLA10" )
R01 = Integral R02 to R04: temp
• R05 = 3.084411158 E-01
• R06 = 1.377934705 E-01
• R07 = 4.011199292 E-01
• R08 = 7.294545495 E-01
• R09 = 2.180682876 E-01
• R10 = 1.808342902
• R11 = 6.208745610 E-02
• R12 = 3.401433698
• R13 = 9.501516975 E-03
• R14 = 5.552496140
• R15 = 7.530083886 E-04
• R16 = 8.330152747
• R17 = 2.825923350 E-05
• R18 = 11.84378584
• R19 = 4.249313985 E-07
• R20 = 16.27925783
• R21 = 1.839564824 E-09
• R22 = 21.99658581
• R23 = 9.911827220 E-13
• R24 = 29.92069701
Flags: /
Subroutine: A program that takes x in X-register
& y in Y-register and returns f(x,y) in X-register
| 01 LBL 'DGLA10"
02 CLX 03 STO 01 04 24.004 05 STO 03 06 STO 04 07 LBL 01 08 CLX 09 STO 02 10 LBL 02 11 RCL IND 04 12 RCL IND 03 13 XEQ IND 00 14 DSE 04 15 RCL IND 04 16 * 17 ST+ 02 18 DSE 04 19 GTO 02 20 RCL 02 21 DSE 03 22 RCL IND 03 23 * 24 ST+ 01 25 20 26 ST+ 04 27 DSE 03 28 GTO 01 29 RCL 01 30 END |
( 61 bytes / SIZE 025 )
| STACK | INPUT | OUTPUT |
| X | / | I |
where I = §0+¥ §0+¥ exp(-x-y) f(x,y) dx dy
Example: Evaluate I = §0+¥ §0+¥ exp(-x-y) Sin ( x + y ) dx dy
LBL "T"
+
SIN
RTN
XEQ "RAD" "T" ASTO 00
XEQ "DGLA10" >>>> I ~ 0.500000715
---Execution time = 1m45s---
Note:
-A quite similar program will use the Gauss-Laguerre 15-point formula:
only lines 04 and 25 must be modified
Data Registers: • R00 = function name ( Registers R00 and R05 thru R34 are to be initialized before executing "DGLA15" )
R01 = Integral R02 to R04: temp
• R05 = 2.182348859 E-01
• R06 = 0.09330781202
• R07 = 3.422101779 E-01
• R08 = 0.4926917403
• R09 = 2.630275779 E-01
• R10 = 1.215595412
• R11 = 1.264258181 E-01
• R12 = 2.269949526
• R13 = 4.020686492 E-02
• R14 = 3.667622722
• R15 = 8.563877804 E-03
• R16 = 5.425336627
• R17 = 1.212436147 E-03
• R18 = 7.565916227
• R19 = 1.116743923 E-04
• R20 = 10.12022857
• R21 = 6.459926762 E-06
• R22 = 13.13028248
• R23 = 2.226316907 E-07
• R24 = 16.65440771
• R25 = 4.227430385 E-09
• R26 = 20.77647890
• R27 = 3.921897267 E-11
• R28 = 25.62389423
• R29 = 1.456515264 E-13
• R30 = 31.40751917
• R31 = 1.483027051 E-16
• R32 = 38.53068331
• R33 = 1.600594906 E-20
• R34 = 48.02608557
Flags: /
Subroutine: A program that takes x in X-register
& y in Y-register and returns f(x,y) in X-register
| 01 LBL 'DGLA15"
02 CLX 03 STO 01 04 34.004 05 STO 03 06 STO 04 07 LBL 01 08 CLX 09 STO 02 10 LBL 02 11 RCL IND 04 12 RCL IND 03 13 XEQ IND 00 14 DSE 04 15 RCL IND 04 16 * 17 ST+ 02 18 DSE 04 19 GTO 02 20 RCL 02 21 DSE 03 22 RCL IND 03 23 * 24 ST+ 01 25 30 26 ST+ 04 27 DSE 03 28 GTO 01 29 RCL 01 30 END |
( 61 bytes / SIZE 025 )
| STACK | INPUT | OUTPUT |
| X | / | I |
where I = §0+¥ §0+¥ exp(-x-y) f(x,y) dx dy
Example: Evaluate again I = §0+¥ §0+¥ exp(-x-y) Sin ( x + y ) dx dy
LBL "T"
+
SIN
RTN
XEQ "RAD" "T" ASTO 00
XEQ "DGLA15" >>>> I ~ 0.5000000001
---Execution time = 3m54s---
c) Gauss-Hermite Formulae
-The Gauss-Hermite 20-point & 30-point formulae may also be used
to compute I = §-¥+¥
§-¥+¥
exp(-x2-y2) f(x,y) dx dy
Data Registers: • R00 = function name ( Registers R00 and R06 thru R25 are to be initialized before executing "DGH20" )
R01 = Integral R02 to R05: temp
• R06 = 4.622436696 E-01
• R07 = 0.2453407083
• R08 = 2.866755054 E-01
• R09 = 0.7374737285
• R10 = 1.090172060 E-01
• R11 = 1.234076215
• R12 = 2.481052089 E-02
• R13 = 1.738537712
• R14 = 3.243773342 E-03
• R15 = 2.254974002
• R16 = 2.283386360 E-04
• R17 = 2.788806058
• R18 = 7.802556479 E-06
• R19 = 3.347854567
• R20 = 1.086069371 E-07
• R21 = 3.944764040
• R22 = 4.399340992 E-10
• R23 = 4.603682450
• R24 = 2.229393646 E-13
• R25 = 5.387480890
Flags: /
Subroutine: A program that takes x in X-register
& y in Y-register and returns f(x,y) in X-register
| 01 LBL "DGH20"
02 CLX 03 STO 01 04 25.005 05 STO 04 06 STO 05 07 LBL 01 08 CLX 09 STO 02 10 LBL 02 |
11 RCL IND 05
12 RCL IND 04 13 XEQ IND 00 14 STO 03 15 RCL IND 05 16 RCL IND 04 17 CHS 18 XEQ IND 00 19 ST+ 03 20 RCL IND 05 |
21 CHS
22 RCL IND 04 23 XEQ IND 00 24 ST+ 03 25 RCL IND 05 26 CHS 27 RCL IND 04 28 CHS 29 XEQ IND 00 30 RCL 03 |
31 +
32 DSE 05 33 RCL IND 05 34 * 35 ST+ 02 36 DSE 05 37 GTO 02 38 RCL 02 39 DSE 04 40 RCL IND 04 |
41 *
42 ST+ 01 43 20 44 ST+ 05 45 DSE 04 46 GTO 01 47 RCL 01 48 END |
( 89 bytes / SIZE 026 )
| STACK | INPUT | OUTPUT |
| X | / | I |
where I = §-¥+¥ §-¥+¥ exp(-x2-y2) f(x,y) dx dy
Example: Evaluate I = §-¥+¥ §-¥+¥ exp(-x2-y2) Cos ( x + y ) dx dy
LBL "T"
+
COS
RTN
XEQ "RAD" "T" ASTO 00
XEQ "DGH20" >>>> I ~ 1.905472266
---Execution time = 5m26s---
-The exact result is PI / sqrt(e) = 1.905472265 rounded to 9 decimals
-And with the 30-point Gauss-Hermite formula:
Data Registers: • R00 = function name ( Registers R00 and R06 thru R25 are to be initialized before executing "DGH30" )
R01 = Integral R02 to R05: temp
• R06 = 3.863948895 E-01
• R07 = 0.2011285765
• R08 = 2.801309308 E-01
• R09 = 0.6039210586
• R10 = 1.467358475 E-01
• R11 = 1.008338271
• R12 = 5.514417687 E-02
• R13 = 1.415527800
• R14 = 1.470382970 E-02
• R15 = 1.826741144
• R16 = 2.737922473 E-03
• R17 = 2.243391468
• R18 = 3.483101243 E-04
• R19 = 2.667132125
• R20 = 2.938725229 E-05
• R21 = 3.099970530
• R22 = 1.579094887 E-06
• R23 = 3.544443873
• R24 = 5.108522451 E-08
• R25 = 4.003908604
• R26 = 9.178580424 E-10
• R27 = 4.483055357
• R28 = 8.106186297 E-12
• R29 = 4.988918969
• R30 = 2.878607081 E-14
• R31 = 5.533147152
• R32 = 2.810333603 E-17
• R33 = 6.138279220
• R34 = 2.908254700 E-21
• R35 = 6.863345294
Flags: /
Subroutine: A program that takes x in X-register
& y in Y-register and returns f(x,y) in X-register
| 01 LBL "DGH30"
02 CLX 03 STO 01 04 35.005 05 STO 04 06 STO 05 07 LBL 01 08 CLX 09 STO 02 10 LBL 02 |
11 RCL IND 05
12 RCL IND 04 13 XEQ IND 00 14 STO 03 15 RCL IND 05 16 RCL IND 04 17 CHS 18 XEQ IND 00 19 ST+ 03 20 RCL IND 05 |
21 CHS
22 RCL IND 04 23 XEQ IND 00 24 ST+ 03 25 RCL IND 05 26 CHS 27 RCL IND 04 28 CHS 29 XEQ IND 00 30 RCL 03 |
31 +
32 DSE 05 33 RCL IND 05 34 * 35 ST+ 02 36 DSE 05 37 GTO 02 38 RCL 02 39 DSE 04 40 RCL IND 04 |
41 *
42 ST+ 01 43 30 44 ST+ 05 45 DSE 04 46 GTO 01 47 RCL 01 48 END |
( 89 bytes / SIZE 036 )
| STACK | INPUT | OUTPUT |
| X | / | I |
where I = §-¥+¥ §-¥+¥ exp(-x2-y2) f(x,y) dx dy
Example: Evaluate again I = §-¥+¥ §-¥+¥ exp(-x2-y2) Cos ( x + y ) dx dy
LBL "T"
+
COS
RTN
XEQ "RAD" "T" ASTO 00
XEQ "DGH30" >>>> I ~ 1.905472265
---Execution time = 12m14s---
-"DGH30" gives the exact result = PI / sqrt(e) = 1.905472265
rounded to 9 decimals.
5°) Triple Integrals
a) Gauss-Legendre Formulae
a-1)
Triple Integrals Only
-"TGL10" uses the Gauss-Legendre 10-point formula to evaluate:
§ab
§u(x)v(x)
§w(x,y)t(x,y) f(x,y,z) dx dy dz
Data Registers: • R00 = f name ( Registers R00 thru R07 are to be initialized before executing "TGL10" )
• R01 = a
• R04 = u name
• R06 = w name
• R02 = b
• R05 = v name
• R07 = t name
• R03 = n
R08 = Integral
R09 to R27: temp
Flags: /
Subroutines:
-A program that takes x , y , z in registers X , Y , Z and returns
f(x,y,z) in X-register
-A program that takes x in X-register and returns u(x)
in X-register
-A program that takes x in X-register and returns v(x)
in X-register
-A program that takes x , y in registers X & Y and returns
w(x,y) in X-register
-A program that takes x , y in registers X & Y and returns
t(x,y) in X-register
-If you don't want to use the M-Code routines GX & GW, replace:
lines 131-132 by DSE 15 RCL IND 15 , lines
79-80 by DSE 14 RCL IND 14 , lines 32-33 by
DSE 13 RCL IND 13
lines 117-118 by RCL IND 15 , lines 69-70
by RCL IND 14 , lines 22-23 by RCL IND 13
lines 16-63-111 by 37.027
and store the following constants into R28 to R37:
R28 = 0.2955242247
R29 = 0.1488743390
R30 = 0.2692667193
R31 = 0.4333953941
R32 = 0.2190863625
R33 = 0.6794095683
R34 = 0.1494513492
R35 = 0.8650633667
R36 = 0.06667134431
R37 = 0.9739065285
| 01 LBL "TGL10"
02 LBL 00 03 CLX 04 STO 08 05 RCL 02 06 RCL 01 07 STO 16 08 - 09 RCL 03 10 STO 19 11 ST+ X 12 / 13 STO 09 14 LBL 01 15 ST+ 16 16 5 17 STO 13 18 LBL 02 19 RCL 16 20 STO 10 21 RCL 09 22 RCL 13 23 GX 24 * 25 ST+ 10 26 - 27 XEQ 03 28 X<> 10 29 XEQ 03 |
30 RCL 10
31 + 32 RCL 13 33 GW 34 * 35 ST+ 08 36 DSE 13 37 GTO 02 38 RCL 09 39 ST+ 16 40 DSE 19 41 GTO 01 42 ST* 08 43 RCL 08 44 RTN 45 GTO 00 46 LBL 03 47 STO 22 48 XEQ IND 04 49 STO 17 50 CLX 51 STO 23 52 RCL 22 53 XEQ IND 05 54 RCL 17 55 - 56 RCL 03 57 STO 20 58 ST+ X |
59 /
60 STO 11 61 LBL 04 62 ST+ 17 63 5 64 STO 14 65 LBL 05 66 RCL 17 67 STO 12 68 RCL 11 69 RCL 14 70 GX 71 * 72 ST+ 12 73 - 74 XEQ 06 75 X<> 12 76 XEQ 06 77 RCL 12 78 + 79 RCL 14 80 GW 81 * 82 ST+ 23 83 DSE 14 84 GTO 05 85 RCL 11 86 ST+ 17 87 DSE 20 |
88 GTO 04
89 RCL 23 90 * 91 RTN 92 LBL 06 93 STO 24 94 RCL 22 95 XEQ IND 06 96 STO 18 97 CLX 98 STO 25 99 RCL 24 100 RCL 22 101 XEQ IND 07 102 RCL 18 103 - 104 RCL 03 105 STO 21 106 ST+ X 107 / 108 STO 26 109 LBL 07 110 ST+ 18 111 5 112 STO 15 113 LBL 08 114 RCL 18 115 STO 27 116 RCL 26 |
117 RCL 15
118 GX 119 * 120 ST+ 27 121 - 122 RCL 24 123 RCL 22 124 XEQ IND 00 125 X<> 27 126 RCL 24 127 RCL 22 128 XEQ IND 00 129 RCL 27 130 + 131 RCL 15 132 GW 133 * 134 ST+ 25 135 DSE 15 136 GTO 08 137 RCL 26 138 ST+ 18 139 DSE 21 140 GTO 07 141 RCL 25 142 * 143 RTN 144 END |
( 236 bytes / SIZE 028 )
| STACK | INPUTS | OUTPUTS |
| X | / | I |
where I = §ab §u(x)v(x) §w(x,y)t(x,y) f(x,y,z) dx dy dz
Example: Evaluate I = §12
§xx^2
§x+yx.y ( x2+y2+z2
)1/2 dx dy dz
01 LBL "S"
02 X^2
03 X<>Y
04 X^2
05 +
06 X<>Y
07 X^2
08 +
09 SQRT
10 RTN
11 LBL "U"
12 RTN
13 LBL "V"
14 X^2
15 RTN
16 LBL "W"
17 +
18 RTN
19 LBL "T"
20 *
21 RTN
"S" 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 "TGL10" >>>> 0.813302277
in X-register & R08
---Execution time = 15mn58s---
n = 2 2 STO 03
R/S >>>>
0.813302283 in X-register & R08
n = 4 4 STO 03
R/S >>>>
0.813302278 in X-register & R08
Notes:
-Due to roundoff-errors, it's difficult to decide which result is best
! may be the first one ?
-If you've written GX & GW for, say, Gauss-Legendre 18-point formula,
simply replace 5 ( lines 16-63-111 ) by 9
-If you prefer the improved versions of "GX" & "GW" as suggested in the note paragraph 1°)a)
Delete line 133 & replace lines 114 to 121 by
RCL 18
RCL 26
RCL 15
GX
STO 27
X<>Y
delete line 81 & replace lines 66 to 73 by
RCL 17
RCL 11
RCL 14
GX
STO 12
X<>Y
and delete line 34 & replace lines 19 to 26 by
RCL 16
RCL 09
RCL 13
GX
STO 10
X<>Y
a-2)
Simple, Double & Triple Integrals
-We can combine "GL10" , "DGL10" & "TGL10" to save memory room:
Data Registers: • R00 = f name ( Registers R00 thru R03 or R05 or R07 are to be initialized before executing "SDTI" )
• R01 = a
• R04 = u name
• R06 = w name
• R02 = b
• R05 = v name
• R07 = t name
• R03 = n
R08 = Integral
R09 to R27: temp
Flags: F01-F02
SF01
= Simple Integral
CF 01 & SF 02 = Double Inegral
CF 01 & CF 02 = Triple Integral
Subroutines: Programs that compute f , u , v , w , t assuming x , y , z are in registers X , Y , Z upon entry.
-If you don't want to use the M-Code routines GX & GW, replace:
RCL nn GW by DSE nn
RCL IND nn ( lines 136-137 , 105-106 , 51-52 )
RCL nn GX by RCL IND
nn ( lines 126-127 , 91-92 , 41-42 )
line 04 by 38.028
and store the following constants into R29 to R38:
R29 = 0.2955242247
R30 = 0.1488743390
R31 = 0.2692667193
R32 = 0.4333953941
R33 = 0.2190863625
R34 = 0.6794095683
R35 = 0.1494513492
R36 = 0.8650633667
R37 = 0.06667134431
R38 = 0.9739065285
-Line 15 is a three-byte GTO 01 ( not really important since this instruction
is executed only once )
| 01 LBL "SDTI"
02 CLX 03 STO 08 04 5 05 STO 12 06 RCL 02 07 RCL 01 08 STO 09 09 - 10 RCL 03 11 STO 10 12 ST+ X 13 / 14 STO 11 15 GTO 01 16 LBL 03 17 FS? 01 18 GTO IND 00 19 STO 15 20 XEQ IND 04 21 STO 16 22 CLX 23 STO 17 24 RCL 15 25 XEQ IND 05 26 RCL 16 27 - 28 RCL 03 29 STO 18 30 ST+ X |
31 /
32 STO 19 33 LBL 04 34 ST+ 16 35 RCL 12 36 STO 20 37 LBL 05 38 RCL 16 39 STO 21 40 RCL 19 41 RCL 20 42 GX 43 * 44 ST+ 21 45 - 46 XEQ 06 47 X<> 21 48 XEQ 06 49 RCL 21 50 + 51 RCL 20 52 GW 53 * 54 ST+ 17 55 DSE 20 56 GTO 05 57 RCL 19 58 ST+ 16 59 DSE 18 60 GTO 04 |
61 RCL 17
62 * 63 RTN 64 LBL 06 65 STO 22 66 RCL 15 67 FS? 02 68 GTO IND 00 69 XEQ IND 06 70 STO 23 71 CLX 72 STO 28 73 RCL 22 74 RCL 15 75 XEQ IND 07 76 RCL 23 77 - 78 RCL 03 79 STO 24 80 ST+ X 81 / 82 STO 25 83 LBL 07 84 ST+ 23 85 RCL 12 86 STO 26 87 LBL 08 88 RCL 23 89 STO 27 90 RCL 25 |
91 RCL 26
92 GX 93 * 94 ST+ 27 95 - 96 RCL 22 97 RCL 15 98 XEQ IND 00 99 X<> 27 100 RCL 22 101 RCL 15 102 XEQ IND 00 103 RCL 27 104 + 105 RCL 26 106 GW 107 * 108 ST+ 28 109 DSE 26 110 GTO 08 111 RCL 25 112 ST+ 23 113 DSE 24 114 GTO 07 115 RCL 28 116 * 117 RTN 118 LBL 01 119 ST+ 09 120 RCL 12 |
121 STO 13
122 LBL 02 123 RCL 09 124 STO 14 125 RCL 11 126 RCL 13 127 GX 128 * 129 ST+ 14 130 - 131 XEQ 03 132 X<> 14 133 XEQ 03 134 RCL 14 135 + 136 RCL 13 137 GW 138 * 139 ST+ 08 140 DSE 13 141 GTO 02 142 RCL 11 143 ST+ 09 144 DSE 10 145 GTO 01 146 ST* 08 147 RCL 08 148 END |
( 246 bytes / SIZE 015 or 023 or 029 )
| STACK | INPUT | OUTPUT |
| X | / | I |
where
I = §ab f(x) dx
if SF 01
or
I = §ab
§u(x)v(x)
f(x;y) dx dy
if CF 01 & SF 02
or I =
§ab
§u(x)v(x)
§w(x,y)t(x,y) f(x,y,z) dx dy dz
if CF 01 & CF 02
Notes:
-The same examples give the same results. I is in X-register and in R08.
-Registers R04-R05-R06-R07 are unused to calculate simple integrals.
-Registers R06-R07 are unused to calculate double integrals.
-Instead of using flags F01 & F02 to determine the type of integrals,
one may use flags F02 & F03:
CF 02 CF 03 for simple integrals, SF 02 for double
integrals, CF 02 SF 03 for triple integrals.
-In this case, replace line 17 by FC? 02 FS? 03
FS? 30
-If you prefer the improved versions of "GX" & "GW" as suggested in the note paragraph 1°)a)
Delete line 138 & replace lines 123 to 130 by
RCL 09
RCL 11
RCL 13
GX
STO 14
X<>Y
delete line 107 & replace lines 88 to 95 by
RCL 23
RCL 25
RCL 26
GX
STO 27
X<>Y
and delete line 53 & replace lines 38 to 45 by
RCL 16
RCL 19
RCL 20
GX
STO 21
X<>Y
-The execution time ( for n = 1 ) becomes 15m28s instead of 15m58s.
b) Gauss-Laguerre Formulae
-We now try to compute I = §0+¥
§0+¥ §0+¥
exp(-x-y-z) f(x,y,z) dx dy dz
Data Registers: • R00 = function name ( Registers R00 and R07 thru R26 are to be initialized before executing "TGLA10" )
R01 = Integral R02 to R06: temp
• R07 = 3.084411158 E-01
• R08 = 1.377934705 E-01
• R09 = 4.011199292 E-01
• R10 = 7.294545495 E-01
• R11 = 2.180682876 E-01
• R12 = 1.808342902
• R13 = 6.208745610 E-02
• R14 = 3.401433698
• R15 = 9.501516975 E-03
• R16 = 5.552496140
• R17 = 7.530083886 E-04
• R18 = 8.330152747
• R19 = 2.825923350 E-05
• R20 = 11.84378584
• R21 = 4.249313985 E-07
• R22 = 16.27925783
• R23 = 1.839564824 E-09
• R24 = 21.99658581
• R25 = 9.911827220 E-13
• R26 = 29.92069701
Flags: /
Subroutine: A program that takes x in X-register
, y in Y-register , z in Z-register and returns f(x,y,z) in X-register
| 01 LBL "TGLA10"
02 CLX 03 STO 01 04 26.006 05 STO 04 06 STO 05 07 STO 06 08 LBL 01 09 CLX |
10 STO 02
11 LBL 02 12 CLX 13 STO 03 14 LBL 03 15 RCL IND 06 16 RCL IND 05 17 RCL IND 04 18 XEQ IND 00 |
19 DSE 06
20 RCL IND 06 21 * 22 ST+ 03 23 DSE 06 24 GTO 03 25 RCL 03 26 DSE 05 27 RCL IND 05 |
28 *
29 ST+ 02 30 20 31 ST+ 06 32 DSE 05 33 GTO 02 34 ST+ 05 35 RCL 02 36 DSE 04 |
37 RCL IND 04
38 * 39 ST+ 01 40 DSE 04 41 GTO 01 42 RCL 01 43 END |
( 80 bytes / SIZE 027 )
| STACK | INPUT | OUTPUT |
| X | / | I |
where I = §0+¥ §0+¥ §0+¥ exp(-x-y-z) f(x,y,z) dx dy dz
Example: Evaluate I = §0+¥§0+¥ §0+¥ exp(-x-y-z) Sin ( x + y + z ) dx dy dz
LBL "T"
+
+
SIN
RTN
XEQ "RAD" "T" ASTO 00
XEQ "TGLA10" >>>> I ~ 0.250000765
---Execution time = 18m51s---
Note:
-And with the Gauss-Laguerre 15-point formula.
Data Registers: • R00 = function name ( Registers R00 and R07 thru R36 are to be initialized before executing "TGLA15" )
R01 = Integral R02 to R06: temp
• R07 = 2.182348859 E-01
• R08 = 0.09330781202
• R09 = 3.422101779 E-01
• R10 = 0.4926917403
• R11 = 2.630275779 E-01
• R12 = 1.215595412
• R13 = 1.264258181 E-01
• R14 = 2.269949526
• R15 = 4.020686492 E-02
• R16 = 3.667622722
• R17 = 8.563877804 E-03
• R18 = 5.425336627
• R19 = 1.212436147 E-03
• R20 = 7.565916227
• R21 = 1.116743923 E-04
• R22 = 10.12022857
• R23 = 6.459926762 E-06
• R24 = 13.13028248
• R25 = 2.226316907 E-07
• R26 = 16.65440771
• R27 = 4.227430385 E-09
• R28 = 20.77647890
• R29 = 3.921897267 E-11
• R30 = 25.62389423
• R31 = 1.456515264 E-13
• R32 = 31.40751917
• R33 = 1.483027051 E-16
• R34 = 38.53068331
• R35 = 1.600594906 E-20
• R36 = 48.02608557
Flags: /
Subroutine: A program that takes x in X-register
, y in Y-register , z in Z-register and returns f(x,y,z) in X-register
| 01 LBL "TGLA15"
02 CLX 03 STO 01 04 36.006 05 STO 04 06 STO 05 07 STO 06 08 LBL 01 09 CLX |
10 STO 02
11 LBL 02 12 CLX 13 STO 03 14 LBL 03 15 RCL IND 06 16 RCL IND 05 17 RCL IND 04 18 XEQ IND 00 |
19 DSE 06
20 RCL IND 06 21 * 22 ST+ 03 23 DSE 06 24 GTO 03 25 RCL 03 26 DSE 05 27 RCL IND 05 |
28 *
29 ST+ 02 30 30 31 ST+ 06 32 DSE 05 33 GTO 02 34 ST+ 05 35 RCL 02 36 DSE 04 |
37 RCL IND 04
38 * 39 ST+ 01 40 DSE 04 41 GTO 01 42 RCL 01 43 END |
( 80 bytes / SIZE 027 )
| STACK | INPUT | OUTPUT |
| X | / | I |
where I = §0+¥ §0+¥ §0+¥ exp(-x-y-z) f(x,y,z) dx dy dz
Example: Evaluate I = §0+¥§0+¥ §0+¥ exp(-x-y-z) Sin ( x + y + z ) dx dy dz
LBL "T"
+
+
SIN
RTN
XEQ "RAD" "T" ASTO 00 XEQ "TGLA15" >>>> I ~ 0.2500000000 ---Execution time = 1h02m33s---
Note:
-This is the exact result: no roundoff-error at all !
c) Gauss-Hermite Formulae
-We try to evaluate I = §-¥+¥
§-¥+¥
§-¥+¥
exp(-x2-y2-z2) f(x,y,z) dx dy dz
with the 20-point formula
Data Registers: • R00 = function name ( Registers R00 and R08 thru R27 are to be initialized before executing "TGH20" )
R01 = Integral R02 to R07: temp
• R08 = 4.622436696 E-01
• R09 = 0.2453407083
• R10 = 2.866755054 E-01
• R11 = 0.7374737285
• R12 = 1.090172060 E-01
• R13 = 1.234076215
• R14 = 2.481052089 E-02
• R15 = 1.738537712
• R16 = 3.243773342 E-03
• R17 = 2.254974002
• R18 = 2.283386360 E-04
• R19 = 2.788806058
• R20 = 7.802556479 E-06
• R21 = 3.347854567
• R22 = 1.086069371 E-07
• R23 = 3.944764040
• R24 = 4.399340992 E-10
• R25 = 4.603682450
• R26 = 2.229393646 E-13
• R27 = 5.387480890
Flags: /
Subroutine: A program that takes x in X-register
, y in Y-register & z in Z-register and returns f(x,y,z) in X-register
-Lines 82 and 90 are three-byte GTO's
| 01 LBL "TGH20"
02 CLX 03 STO 01 04 27.007 05 STO 05 06 STO 06 07 STO 07 08 LBL 01 09 CLX 10 STO 02 11 LBL 02 12 CLX 13 STO 03 14 LBL 03 15 RCL IND 07 16 RCL IND 06 17 RCL IND 05 18 XEQ IND 00 19 STO 04 |
20 RCL IND 07
21 RCL IND 06 22 RCL IND 05 23 CHS 24 XEQ IND 00 25 ST+ 04 26 RCL IND 07 27 RCL IND 06 28 CHS 29 RCL IND 05 30 XEQ IND 00 31 ST+ 04 32 RCL IND 07 33 RCL IND 06 34 CHS 35 RCL IND 05 36 CHS 37 XEQ IND 00 38 ST+ 04 |
39 RCL IND 07
40 CHS 41 RCL IND 06 42 RCL IND 05 43 XEQ IND 00 44 ST+ 04 45 RCL IND 07 46 CHS 47 RCL IND 06 48 RCL IND 05 49 CHS 50 XEQ IND 00 51 ST+ 04 52 RCL IND 07 53 CHS 54 RCL IND 06 55 CHS 56 RCL IND 05 57 XEQ IND 00 |
58 ST+ 04
59 RCL IND 07 60 CHS 61 RCL IND 06 62 CHS 63 RCL IND 05 64 CHS 65 XEQ IND 00 66 RCL 04 67 + 68 DSE 07 69 RCL IND 07 70 * 71 ST+ 03 72 DSE 07 73 GTO 03 74 RCL 03 75 DSE 06 76 RCL IND 06 |
77 *
78 ST+ 02 79 20 80 ST+ 07 81 DSE 06 82 GTO 02 83 ST+ 06 84 RCL 02 85 DSE 05 86 RCL IND 05 87 * 88 ST+ 01 89 DSE 05 90 GTO 01 91 RCL 01 92 END |
( 165 bytes / SIZE 028 )
| STACK | INPUT | OUTPUT |
| X | / | I |
where I = §-¥+¥ §-¥+¥ §-¥+¥ exp(-x2-y2-z2) f(x,y,z) dx dy dz
Example: Evaluate I = I = §-¥+¥ §-¥+¥ §-¥+¥ exp(-x2-y2-z2) Cos ( x + y + z ) dx dy dz
LBL "T"
+
+
COS
RTN
XEQ "RAD" "T" ASTO 00 XEQ "TGH20" >>>> I ~ 2.630291901
-The exact result is I = ( PI / sqrt(e) )3/2 = 2.630291900 ( rounded to 9 decimals )
Notes:
-I've used V41 to get this result.
-With a real HP-41, the execution time is probably of the order of
2 hours.
-The Gauss-Hermite 30-point formula leads to a similar program:
-Just change lines 01-04 and 79
Data Registers: • R00 = function name ( Registers R00 and R08 thru R37 are to be initialized before executing "TGH30" )
R01 = Integral R02 to R07: temp
• R08 = 3.863948895 E-01
• R09 = 0.2011285765
• R10 = 2.801309308 E-01
• R11 = 0.6039210586
• R12 = 1.467358475 E-01
• R13 = 1.008338271
• R14 = 5.514417687 E-02
• R15 = 1.415527800
• R16 = 1.470382970 E-02
• R17 = 1.826741144
• R18 = 2.737922473 E-03
• R19 = 2.243391468
• R20 = 3.483101243 E-04
• R21 = 2.667132125
• R22 = 2.938725229 E-05
• R23 = 3.099970530
• R24 = 1.579094887 E-06
• R25 = 3.544443873
• R26 = 5.108522451 E-08
• R27 = 4.003908604
• R28 = 9.178580424 E-10
• R29 = 4.483055357
• R30 = 8.106186297 E-12
• R31 = 4.988918969
• R32 = 2.878607081 E-14
• R33 = 5.533147152
• R34 = 2.810333603 E-17
• R35 = 6.138279220
• R36 = 2.908254700 E-21
• R37 = 6.863345294
Flags: /
Subroutine: A program that takes x in X-register
, y in Y-register & z in Z-register and returns f(x,y,z) in X-register
-Lines 82 and 90 are three-byte GTO's
| 01 LBL "TGH30"
02 CLX 03 STO 01 04 37.007 05 STO 05 06 STO 06 07 STO 07 08 LBL 01 09 CLX 10 STO 02 11 LBL 02 12 CLX 13 STO 03 14 LBL 03 15 RCL IND 07 16 RCL IND 06 17 RCL IND 05 18 XEQ IND 00 19 STO 04 |
20 RCL IND 07
21 RCL IND 06 22 RCL IND 05 23 CHS 24 XEQ IND 00 25 ST+ 04 26 RCL IND 07 27 RCL IND 06 28 CHS 29 RCL IND 05 30 XEQ IND 00 31 ST+ 04 32 RCL IND 07 33 RCL IND 06 34 CHS 35 RCL IND 05 36 CHS 37 XEQ IND 00 38 ST+ 04 |
39 RCL IND 07
40 CHS 41 RCL IND 06 42 RCL IND 05 43 XEQ IND 00 44 ST+ 04 45 RCL IND 07 46 CHS 47 RCL IND 06 48 RCL IND 05 49 CHS 50 XEQ IND 00 51 ST+ 04 52 RCL IND 07 53 CHS 54 RCL IND 06 55 CHS 56 RCL IND 05 57 XEQ IND 00 |
58 ST+ 04
59 RCL IND 07 60 CHS 61 RCL IND 06 62 CHS 63 RCL IND 05 64 CHS 65 XEQ IND 00 66 RCL 04 67 + 68 DSE 07 69 RCL IND 07 70 * 71 ST+ 03 72 DSE 07 73 GTO 03 74 RCL 03 75 DSE 06 76 RCL IND 06 |
77 *
78 ST+ 02 79 30 80 ST+ 07 81 DSE 06 82 GTO 02 83 ST+ 06 84 RCL 02 85 DSE 05 86 RCL IND 05 87 * 88 ST+ 01 89 DSE 05 90 GTO 01 91 RCL 01 92 END |
( 165 bytes / SIZE 038 )
| STACK | INPUT | OUTPUT |
| X | / | I |
where I = §-¥+¥ §-¥+¥ §-¥+¥ exp(-x2-y2-z2) f(x,y,z) dx dy dz
Example: Evaluate I = I = §-¥+¥ §-¥+¥ §-¥+¥ exp(-x2-y2-z2) Cos ( x + y + z ) dx dy dz
LBL "T"
+
+
COS
RTN
XEQ "RAD" "T" ASTO 00 XEQ "TGH20" >>>> I ~ 2.630291899
-The exact result is I = ( PI / sqrt(e) )3/2 = 2.630291900
( rounded to 9 decimals )
-"TGH20" & "TGH30" are equivalent for this integral.
Notes:
-I've again used V41 to get this result.
-With a real HP-41, the execution time is probably more than 6 hours.
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 C++" - Cambridge University Press - ISBN 0-521-75033-4