Swiss Micros DM41X: Infinite Integrals by Gaussian Quadrature

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. 

Casio fx-9750GIII: Integrals with Infinite Limits

Casio fx-9750GIII:  Integrals with Infinite Limits

A Substitution to Get to Infinity... 

It is quite a challenge to calculate numerical integrals with infinite limits such as 

∫( f(x) dx, a, ∞)

∫( f(x) dx, -∞, a)

∫( f(x) dx, -∞, ∞)

A trick is to substitute x = tan Θ.  Then:

Θ = arctan x

dx = sec^2 Θ dΘ = cos^-2 Θ dΘ

Note that

lim t→∞ arctan t = π/2

lim t→-∞ arctan t = -π/2  

In this blog, assume that radian angle mode is used in all calculations.

Let's go over some examples and see how it works.  I used a Casio fx-9750GIII, however, this technique should work with all calculators with numerical integral calculations.  Furthermore, changing the integral will allow for Simpson's Rule or Trapezoid Rule approximation.

For the upper limit, I approximate π/2.

A = 1.5708  (π/2 to 4 places)

B = 1.570796  (π/2 to 6 places)

C = 1.57079633 (π/2 to 8 places)

Example 1:  

∫( e^(-x^2) dx, 0, ∞)  

transforms to

∫( e^(-tan^2 Θ)/cos^2 Θ dΘ, 0, ≈π/2)

Calculations (fx-9750GIII):

Upper Limit: A (see above),  Result:  0.8862269255

Upper Limit: B,  Result:  0.8862269255

Upper Limit: C,  Result:  0.8862269255

Example 2:  

∫(x^0.5 * e^(-x) dx, 0, ∞)

transforms to

∫((tan Θ)^0.5 * e^(-tan Θ))/cos^2 Θ dΘ, 0, π/2)

Calculations:

Upper limits A, B, C:  0.862269255

Coincidently, ∫( e^(-x^2) dx, 0, ∞)  = ∫(x^0.5 * e^(-x) dx, 0, ∞) = Γ(1.5) = √(π)/2

Example 3:

∫( e^x*(x^2 + 1) dx, -∞, 0)

transforms to

∫( e^(tan Θ) * (tan^2 Θ + 1)/cos^2 Θ dΘ, -π/2, 0)

=  ∫( e^(tan Θ) * 1/cos^2 Θ * 1/cos^2 Θ dΘ, -π/2, 0)

∫( e^(tan Θ)/cos^4 Θ dΘ, -π/2, 0)

For upper limits A, B, C, the answer returned is 3, which is the exact answer.  

For integrals like this all is needed is a four digit approximation of π/2 = 1.5708.

Let's keep going:

Example 4:

∫( (x^2 + 1)/(x^4 - 1) dx, 0, ∞) 

transforms to

∫( 1/cos^Θ * (tan^2 Θ + 1)/(tan^4 Θ - 1) dΘ, 0, π/2)

= ∫( 1/(sin^4 Θ - cos^4 Θ) dΘ, 0, π/2)

On the fx-9750GIII get MA Error.  

Example 5:

∫( 1/√(x^2 + 3*x + 1) dx, 0, ∞)

transforms to 

∫( 1/(cos^2 Θ * √(tan^2 Θ + 3 * tan Θ + 1)) dΘ, 0, π/2)

like example 4, I get the MA Error.

Observation:  the transformation works best if the integral involves some form of either of the following:

f(x) * e^(g(x))

or 

f(x) * e^(-g(x))

Gamma

The Gamma Function is an excellent candidate for this transformation.  

Γ(t) = ∫( x^(t-1) * e^-x dx, 0, ∞)

transforms to

∫( tan^(t-1) Θ * e^(-tan Θ)/cos^2 Θ dΘ, 0, π/2)

Source

Mier-Jedrzejuwicz, W.A.C. Ph.D.    Tips And Programs for the HP 32S   Synthetix Publication.  Berkeley, CA.  September 1988.  ISBN 0-937637-05-X

Download the book here:  https://literature.hpcalc.org/items/1756

Happy calculating,

Eddie

All original content copyright, © 2011-2022.  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. 

The functions e^x, e^-x, e^(-x^2), erf(x) and Taylor Series

The functions e^x, e^-x, e^(-x^2), erf(x) and Taylor Series

Accurate digits are highlighted in green.  Calculations are used with a TI 84 Plus CE.

e^x = 1 + x + x^2/2! + x^3/3! + x^4/4 + … = Σ(x^n/n!, from n = 0 to ∞)

x =	e^x	10 terms	25 terms	50 terms
1	2.718281828	2.718281801	2.718281828	2.718281828
3	20.08553692	20.07966518	20.08553692	20.08553692
5	148.4131591	146.380601	148.4131591	148.4131591
9.9	19930.37044	11869.50538	19930.07221	19930.37044

e^(-x) = 1 – x + x^2/2! – x^3/3! + x^4/4 - … = Σ( (-x)^n/n!, from n = 0 to ∞)

x =	e^(-x)	10 terms	25 terms	50 terms
1	0.3678794412	0.3678794643	0.3678794412	0.3678794412
3	0.0497870684	0.0533258929	0.0497870684	0.0497870684
5	0.006737947	0.8640390763	0.0067379439	0.006737947
9.9	5.017468206E-5	1207.799663	-0.1392914019	5.017463241E-5

I think you know where I’m going.

e^(-x^2) = 1 – x^2 + x^4/2! – x^6/3! + x^8/4! = Σ( (-x)^(2*n)/n!, from n = 0 to ∞)

x =	e^(-x^2)	10 terms	25 terms	50 terms
1	0.3678794412	0.3678794643	0.3678794412	0.3678794412
3	1.234098041E-4	442.2750223	-0.0118646275	1.234194001E-4
5	1.38879439E-11	18613495.8	-2834107793	85689.40174
9.9	2.72143414E-43	2.04347238E13	-3.10254183E24	7.951057508E34

(Something really goes bonkers as x increases and n increases)

Error Function

erf(x) = 2/√π * ∫(e^(-t^2) dt, 0, x)

= 2/√π * (x – x^3/3 + x^5/(5*2!) – x^7/(7*3!) + x^9/(9*4!) - ...)

= 2/√π * Σ( (-x^(2n+1)/((2n+1)*n!) from n = 0 to ∞ )

x =	erf(x)	10 terms	25 terms	50 terms
1	0.8427007929	0.8427007941	0.8427007929	0.8427007929
3	0.9999779095	68.58627744	0.9992050426	0.9999779095
5	1	4853382.901	-3070260210.4	4724.331354
9.9	1	1.076461715E13	-6.7395908E23	*overflows during calculation* (Result: 8.73442E33 from WolframAlpha) (erf(x) is practically 1 for x > 3)

Note: 9.9^(2*50+1) ≈ 3.623E100

Thoughts:

*  Taylor series are great when x is near its center point.  In the all the cases above, the center point is x = 0. 

*  The more simple the expression, the better range of accuracy with less terms. 

*  Before you recommend a Taylor Series to approximate f(x), check the accuracy and the range.  A cautionary tale. 

Eddie

This blog is property of Edward Shore, 2016.