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)

[Albert Chan]

(11-23-2022 01:55 PM)Eddie W. Shore Wrote:  The program INFGAUS calculates the integral:

∫ e^-x * f(x) dx from x = a to x = ∞

...

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.

The handbook only listed the shifts and weights ... How does the numbers derived ?

I was expecting Gaussian Quadrature of z integral. But with 3 points, results are terrible.
Note: in XCas, Int(...) == quote(int(...))

XCas> Int(e^-x * f(x), x, a, inf) (x=-ln(y))

∫(f(-ln(y)), y, 0, exp(-a))

XCas> Int(f(-ln(y)), y, 0, exp(-a)) (y = (z+1)*exp(-a)/2)

exp(-a) * ∫(f(a - ln((1+z)/2))/2, z, -1, 1)

[Albert Chan]

2nd attempt, points and weights are closer.
We assume a=0, then later put back a.

Let x=(1-y)/(1+y), y=[-1 ..1] map to x = [inf .. 0]
Trivia: inverse has same shape, y=(1-x)/(1+x)

XCas> Int(e^-x * f(x), x, 0, inf) (x=(1-y)/(1+y))

\(\displaystyle \int _{-1}^{1}
\frac{2 \exp \left(- \frac{1-y}{1+y}\right)\;f\left(\frac{1-y}{1+y}\right)}
{\left(1+y\right)^{2}}\,dy\)

Add back a, we have:

\(\displaystyle \int _{a}^{\infty} e^{-x}\,f(x)\,dx = e^{-a}\;\int _{-1}^{1}
\frac{2 \exp \left(- \frac{1-y}{1+y}\right)\;f\left(\frac{1-y}{1+y} + a\right)}
{\left(1+y\right)^{2}}\,dy\)

With 3 points Gaussian quadrature, we have:

e^-a * Σ(( w_i * f(z_i + a)) for i=1 to 3)

where

w_1 = 0.310738836
z_1 = 0.1270166538

w_2 = 0.6540078954
z_2 = 1.0

w_3 = 0.008329973286
z_3 = 7.872983346

Numbers strongly depends on how infinite intervals is mapped to [-1, 1]
For OP examples, above numbers are still not as good. (but, perhaps good for others?)

[Albert Chan]

I figured out how the numbers are derived. The trick is *no* transformation!
Again, we assume a=0, then later put back a (if needed)

Assume f is a quintic polynomial, we solved for integral exact solution.

XCas> P(x) := horner([A,B,C,D,E,F],x) // quintic polynomial
XCas> R := int(e^-x * P(x), x=0 .. inf)

120*A + 24*B + 6*C + 2*D + E+F

XCas> M := w1*P(z1) + w2*P(z2) + w3*P(z3) - R

If weights/points correctly picked, M=0, for any P coefficients.
For simplicity, we use identity matrix for P coefficients.

XCas> M := normal([
M(A=1,B=0,C=0,D=0,E=0,F=0),
M(A=0,B=1,C=0,D=0,E=0,F=0),
M(A=0,B=0,C=1,D=0,E=0,F=0),
M(A=0,B=0,C=0,D=1,E=0,F=0),
M(A=0,B=0,C=0,D=0,E=1,F=0),
M(A=0,B=0,C=0,D=0,E=0,F=1)
])

XCas> transpose(M)

\(\left(\begin{array}{c}-120+w_{1} z_{1}^{5}+w_{2} z_{2}^{5}+w_{3} z_{3}^{5}\\-24+w_{1} z_{1}^{4}+w_{2} z_{2}^{4}+w_{3} z_{3}^{4}\\-6+w_{1} z_{1}^{3}+w_{2} z_{2}^{3}+w_{3} z_{3}^{3}\\-2+w_{1} z_{1}^{2}+w_{2} z_{2}^{2}+w_{3} z_{3}^{2}\\-1+w_{1} z_{1}+w_{2} z_{2}+w_{3} z_{3}\\-1+w_{1}+w_{2}+w_{3}\end{array}\right)\)

Last equation, we have sum of weights = 1, matching OP
6 equations, 6 unknown, using previous post numbers as guess.

XCas> fsolve(M=0,[w1,z1, w2,z2, w3,z3] = [0.311,0.127, 0.654,1.00, 0.00833,7.87], [w1,z1, w2,z2, w3,z3])

[0.278517733569, 2.29428036028,
 0.711093009929, 0.415774556783,
 0.0103892565016, 6.28994508294]

Weights/Points now match OP (in different orders, but that's OK)