Swiss Micros DM41X: Infinite Integrals by Gaussian Quadrature
Introduction
The program INFGAUS calculates the integral:
∫ e^-x * f(x) dx from x = a to x = ∞
f(x) is the subroutine FX. The subroutine starts with the x value on the stack and ends with the RTN command.
The HP 45 algorithm, which is incorporated into the program INFGAUS estimates the integral by the sum
e^-a * Σ(( w_i * f(z_i + a)) for i=1 to 3)
where
w_1 = 0.71109390099
z_1 = 0.4157745568
w_2 = 0.2785177336
z_2 = 2.29428036
w_3 = 0.0103892565
z_3 = 6.289945083
DM41X Program: INFGAUS
The code can work is for the entire HP 41C/DM41 family. The calculator is set to Fix 2 mode. Registers 01 through 08 are needed.
01 LBL^T INFGAUS
02 .4157745568
03 STO 01
04 2.29428036
05 STO 02
06 6.289945083
07 STO 03
08 .7110930099
09 STO 04
10 .278517736
11 STO 05
12 .0103892565
13 STO 06
14 FIX 2
15 ^T A?
16 PROMPT
17 STO 07
18 RCL 01
19 +
20 XEQ^T FX
21 RCL 04
22 *
23 STO 08
24 RCL 07
25 RCL 02
26 +
27 XEQ^T FX
28 RCL 05
29 *
30 ST+ 08
31 RCL 03
32 RCL 07
33 +
34 XEQ^T FX
35 RCL 06
36 *
37 RCL 08
38 +
39 RCL 07
40 CHS
41 E↑X
42 *
43 STO 08
44 END
Examples
Example 1
∫ e^-x * x^3.99 dx from x = 0 to ∞ (calculate Γ(4.99))
LBL^T FX
3.99
Y↑X
RTN
A? 0
Result: 23.64
Example 2
∫ e^-x * x^2 ÷ (x - 1) dx from x = 2 to ∞
LBL^T FX
X↑2
LASTx
1
-
/
RTN
A? 2
Result: 0.62
Example 3
∫ (e^-x)^2 dx from x = 0 to ∞
= ∫ (e^-x) * (e^-x) dx from x = 0 to ∞
LBL^T FX
CHS
E↑X
RTN
A? 0
Result: 0.50 (it turns out 0.5 is the exact answer)
Source
HP-45 Applications Handbook Hewlett Packard Company. 1974.
Remember, my regular posting schedule for January and February 2023 is on Saturdays only.
Eddie
All original content copyright, © 2011-2023. Edward Shore. Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited. This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.