Infinite Integrals by Gaussian Quadrature

[Eddie W. Shore]

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

Source

HP-45 Applications Handbook Hewlett Packard Company. 1974.


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.

Code:

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))

Code:

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 ∞

Code:

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 ∞

Code:

LBL^T FX
CHS
E↑X
RTN

A? 0
Result: 0.50 (it turns out 0.5 is the exact answer)