11-23-2022, 01:55 PM
Post: #1
 Eddie W. Shore Senior Member Posts: 1,383 Joined: Dec 2013
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)
12-14-2022, 08:21 PM
Post: #2
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
RE: Infinite Integrals by Gaussian Quadrature
(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)
12-19-2022, 02:22 PM
Post: #3
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
RE: Infinite Integrals by Gaussian Quadrature
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$$

$$\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?)
12-19-2022, 08:11 PM
Post: #4
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
RE: Infinite Integrals by Gaussian Quadrature
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)
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)