Infinite Integrals by Gaussian Quadrature - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: HP-41C Software Library (/forum-11.html) +--- Thread: Infinite Integrals by Gaussian Quadrature (/thread-19191.html) Infinite Integrals by Gaussian Quadrature - Eddie W. Shore - 11-23-2022 01:55 PM 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) RE: Infinite Integrals by Gaussian Quadrature - Albert Chan - 12-14-2022 08:21 PM (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) RE: Infinite Integrals by Gaussian Quadrature - Albert Chan - 12-19-2022 02:22 PM 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?) RE: Infinite Integrals by Gaussian Quadrature - Albert Chan - 12-19-2022 08:11 PM 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)