Double Exp Integration

Numerical Integration: Double Exponential Formulas for the HP41

Overview

1°) Tanh-Sinh Quadrature §_-1⁺¹ f(x) dx

1-a) With M-Code Routines SINH COSH TANH X#Y?? ( 2 programs )
1-a) Without M-Code Routine ( 2 programs )

2°) Tanh-Sinh Quadrature §_a^b f(x) dx

2-a) With M-Code Routines SINH COSH TANH X#Y?? ( 2 programs )
2-b) Without M-Code Routine ( 2 programs )

3°) Sinh-Sinh Quadrature §_-infinity^+infinity f(x) dx

3-a) With M-Code Routines SINH COSH X#Y??
3-b) Without M-Code Routine

4°) Exp-Sinh Quadrature §₀^+infinity f(x) dx

4-a) With M-Code Routines SINH COSH X#Y??
4-b) Without M-Code Routine

-The 1st versions of these programs employ M-Code routines which are listed in "A Few M-Code Routines for the HP-41"
-The 2nd versions avoid M-Code and calculate hyperbolic functions with the built-in E^X and E^X-1

1°) Tanh-Sinh Quadrature I = §_-1⁺¹ f(x) dx

a) With M-Code Routines SINH COSH TANH X#Y??

Formula: For a given "small" number h , I ~ SUM_{k=...-2,-1,0,1,2,...} w_k f(x_k)

where x_k = Tanh ( Sinh ( k.h ) ) ( The original formulas involve a factor PI / 2 that is omitted here )
and w_k = h Cosh ( k.h ) / Cosh² ( Sinh ( k.h )

Data Registers: • R00 = Function name ( Register R00 is to be initialized before executing "TSI" )

                                         R01 = I = R07                          R03 = k
                                         R02 = h , h/2 , h/4 , ....              R04 to R06: temp
Flags: /
Subroutine:   A program that takes x in X-register and returns f(x) in X-register

01 LBL "TSI"
02 STO 02
03 E99
04 STO 07
05 SIGN
06 STO 06
07 CLX
08 STO 03
09 XEQ IND 00
10 RCL 02
11 *
12 STO 01

13 GTO 01
14 LBL 00
15 2
16 ST/ 01
17 ST/ 02
18 STO 06
19 SIGN
20 CHS
21 STO 03
22 LBL 01
23 RCL 06
24 ST+ 03

25 RCL 02
26 RCL 03
27 *
28 SINH
29 STO 05
30 TANH
31 STO 04
32 XEQ IND 00
33 X<> 04
34 CHS
35 XEQ IND 00
36 RCL 04

37 +
38 RCL 02
39 ST* Y
40 RCL 03
41 *
42 COSH
43 *
44 RCL 05
45 COSH
46 X^2
47 /
48 RCL 01

49 +
50 STO 01
51 LASTX
52 X#Y?
53 GTO 01
54 VIEW X
55 X<> 07
56 X#Y??
57 GTO 00
58 RCL 07
59 END

( 87 bytes / SIZE 008 )

STACK	INPUT	OUTPUT
X	h	I

Example: I = §_-1⁺¹ sqrt(1-x^2) dx

-The subroutine may be

01 LBL "T"
02 X^2
03 SIGN
04 LASTX
05 -
06 SQRT
07 RTN
08 END

"T" ASTO 00 and if you choose h = 0.2

0.2 XEQ "TSI" >>>> 1.570796325 The HP-41 displays the successive approximations
1.570796326 ---Execution time = 109s---

-The exact value is PI / 2

Notes:

-Thogh it will work in many cases, the criterion lines 52-53 is obviously not very good:
-If you choose h = 1 and if f(-1) + f(1) = 0, the program would stop with a wrong result

-A better - though not perfect - criterion could be:

1-Add DSE 08 GTO 02 after line 53
2-Add 3 STO 08 LBL 02 after line 22

3 or another positive integer...

-If the function f has a singularity for x = +/-1 , 1 - x will not calculated with a good precision and the program may fail.

>>> In the following variant, the HP41 uses 1 - x = [ exp ( - sinh x ) } / cosh ( sinh x ) if x > 0

Data Registers: • R00 = Function name ( Register R00 is to be initialized before executing "TSI" )

                                         R01 = I = R07                          R03 = k
                                         R02 = h , h/2 , h/4 , ....              R04 to R09: temp
Flags: /
Subroutine:   A program that takes x in X-register ,
                                                   1 - x in Y-register ,
                                                   1 + x in Z-register and returns f(x) in X-register

01 LBL "TSI"
02 STO 02
03 E99
04 STO 07
05 SIGN
06 ENTER
07 STO 06
08 0
09 STO 03
10 XEQ IND 00
11 RCL 02
12 *
13 STO 01
14 GTO 01
15 LBL 00

16 2
17 ST/ 01
18 ST/ 02
19 STO 06
20 SIGN
21 CHS
22 STO 03
23 LBL 01
24 RCL 06
25 ST+ 03
26 RCL 02
27 RCL 03
28 *
29 SINH
30 STO 05

31 TANH
32 ENTER
33 STO 04
34 1
35 +
36 STO 08
37 RCL 05
38 CHS
39 E^X
40 RCL 05
41 COSH
42 STO 05
43 /
44 STO 09
45 R^

46 XEQ IND 00
47 X<> 04
48 CHS
49 RCL 08
50 RCL 09
51 X<> Z
52 XEQ IND 00
53 RCL 04
54 +
55 RCL 02
56 ST* Y
57 RCL 03
58 *
59 COSH
60 *

61 RCL 05
62 X^2
63 /
64 RCL 01
65 +
66 STO 01
67 LASTX
68 X#Y?
69 GTO 01
70 VIEW X
71 X<> 07
72 X#Y??
73 GTO 00
74 RCL 07
75 END

( 104 bytes / SIZE 010 )

STACK	INPUT	OUTPUT
X	h > 0	I

Example1: I = §_-1⁺¹ 1 / sqrt(1-x^2) dx

-Load this subroutine in main memory

01 LBL "T"
02 RDN
03 *
04 SQRT
05 1/X
06 RTN
07 END

"T" ASTO 00 and if you choose h = 0.2

0.2 XEQ "TSI" >>>> 3.141592655 The HP-41 displays the successive approximations
3.141592654 ---Execution time = 2m40s---

-The exact value is PI

Example2: I = §_-1⁺¹ 1 / sqrt(1-x^4) dx

-Since 1 - x^4 = ( 1 - x ) ( 1 + x ) ( 1 + x^2 ) we can load the following subroutine in main memory

01 LBL "T"
02 X^2
03 SIGN
04 LASTX
05 +
06 *
07 *
08 SQRT
09 1/X
10 RTN
11 END

"T" ASTO 00 and if you choose h = 0.2

0.2 XEQ "TSI" >>>> 2.622057557 The HP-41 displays the successive approximations
2.622057557 ---Execution time = 2m52s---

Note:

-The criterion of lines 68-69 may be improved like with the 1st version:

1-Add DSE 08 GTO 02 after line 69
2-Add 3 STO 08 LBL 02 after line 23

b) Without M-Code Routine

-The same results may be obtained without M-Code:

1st Version:

01 LBL "TSI"
02 STO 02
03 E99
04 STO 07
05 SIGN
06 STO 06
07 CLX
08 STO 03
09 XEQ IND 00
10 RCL 02
11 *
12 STO 01
13 GTO 01
14 LBL 00
15 2
16 ST/ 01

17 ST/ 02
18 STO 06
19 SIGN
20 CHS
21 STO 03
22 LBL 01
23 RCL 06
24 ST+ 03
25 RCL 02
26 RCL 03
27 *
28 E^X-1
29 LASTX
30 CHS
31 E^X-1
32 -

33 STO 05
34 E^X-1
35 STO Y
36 2
37 +
38 /
39 STO 04
40 XEQ IND 00
41 X<> 04
42 CHS
43 XEQ IND 00
44 RCL 04
45 +
46 RCL 02
47 ST* Y
48 RCL 03

49 *
50 E^X
51 ENTER
52 1/X
53 +
54 ST+ X
55 *
56 RCL 05
57 2
58 /
59 E^X
60 ENTER
61 1/X
62 +
63 X^2
64 /

65 RCL 01
66 +
67 STO 01
68 LASTX
69 X#Y?
70 GTO 01
71 VIEW X
72 X<> 07
73 RND
74 X<>Y
75 RND
76 X#Y?
77 GTO 00
78 RCL 07
79 END

( 104 bytes / SIZE 008 )

STACK	INPUT	OUTPUT
X	h	I

-Here, the precision is controlled by the display format
-FIX 9 could lead to unuseful loops and FIX 8 is often enough to get the best accuracy.

-Same remark as above about the criterion of lines 69-70

2nd Version:

01 LBL "TSI"
02 STO 02
03 E99
04 STO 07
05 SIGN
06 ENTER
07 STO 06
08 0
09 STO 03
10 XEQ IND 00
11 RCL 02
12 *
13 STO 01
14 GTO 01
15 LBL 00
16 2
17 ST/ 01
18 ST/ 02
19 STO 06
20 SIGN

21 CHS
22 STO 03
23 LBL 01
24 RCL 06
25 ST+ 03
26 RCL 02
27 RCL 03
28 *
29 E^X-1
30 LASTX
31 CHS
32 E^X-1
33 -
34 STO 05
35 2
36 /
37 X<> 05
38 E^X-1
39 STO Y
40 2

41 +
42 /
43 STO 04
44 RCL 05
45 CHS
46 E^X
47 RCL X
48 RCL 05
49 E^X
50 +
51 2
52 /
53 STO 05
54 /
55 STO 09
56 1
57 R^
58 +
59 STO 08
60 X<> Z

61 XEQ IND 00
62 X<> 04
63 CHS
64 RCL 08
65 RCL 09
66 X<> Z
67 XEQ IND 00
68 RCL 04
69 +
70 RCL 02
71 ST* Y
72 RCL 03
73 *
74 E^X
75 ENTER
76 1/X
77 +
78 *
79 RCL 05
80 X^2

81 ST+ X
82 /
83 RCL 01
84 +
85 STO 01
86 LASTX
87 X#Y?
88 GTO 01
89 VIEW X
90 X<> 07
91 RND
92 X<>Y
93 RND
94 X#Y?
95 GTO 00
96 RCL 07
97 END

( 104 bytes / SIZE 008 )

STACK	INPUT	OUTPUT
X	h > 0	I

-Here, the precision is also controlled by the display format
-FIX 9 could lead to unuseful loops and FIX 8 is often enough to get the best accuracy.

-Same remark as above about the criterion of lines 87-88

2°) Tanh-Sinh Quadrature I = §_a^b f(x) dx

a) With M-Code Routines SINH COSH TANH X#Y??

-We can use the programs above with the change of arguments x = t ( b - a ) / 2 + ( a + b ) / 2
-But your HP-41 can do the job for you.

Data Registers: • R00 = Function name ( Register R00 is to be initialized before executing "TSI2" )

                                         R01 = I = R07                          R03 = k                      R08 = ( a + b ) / 2
                                         R02 = h , h/2 , h/4 , ....              R04 to R06: temp       R09 = ( b - a ) / 2
Flags: /
Subroutine:   A program that takes x in X-register and returns f(x) in X-register

01 LBL "TSI2"
02 STO 02
03 RDN
04 STO 09
05 X<>Y
06 ST- 09
07 +
08 2
09 ST/ 09
10 /
11 STO 08
12 E99
13 STO 07
14 SIGN
15 STO 06
16 CLX

17 STO 03
18 RCL 08
19 XEQ IND 00
20 RCL 02
21 *
22 RCL 09
23 *
24 STO 01
25 GTO 01
26 LBL 00
27 2
28 ST/ 01
29 ST/ 02
30 STO 06
31 SIGN
32 CHS

33 STO 03
34 LBL 01
35 RCL 06
36 ST+ 03
37 RCL 02
38 RCL 03
39 *
40 SINH
41 STO 05
42 TANH
43 RCL 09
44 *
45 STO 04
46 RCL 08
47 +
48 XEQ IND 00

49 X<> 04
50 CHS
51 RCL 08
52 +
53 XEQ IND 00
54 RCL 04
55 +
56 RCL 02
57 ST* Y
58 RCL 03
59 *
60 COSH
61 *
62 RCL 09
63 *
64 RCL 05

65 COSH
66 X^2
67 /
68 RCL 01
69 +
70 STO 01
71 LASTX
72 X#Y?
73 GTO 01
74 VIEW X
75 X<> 07
76 X#Y??
77 GTO 00
78 RCL 07
79 END

( 110 bytes / SIZE 010 )

STACK	INPUT	OUTPUT
Z	a	/
Y	b	/
X	h	I

Example: I = §_-2³ exp ( -x^2 ) dx

01 LBL "T"
02 X^2
03 CHS
04 E^X
05 RTN
06 END

"T" ASTO 00 and if you choose h = 0.2

     -2 ENTER^
      3 ENTER^
    0.2 XEQ "TSI2" >>>>   1.768288737                          The HP-41 displays the successive approximations
                                            1.768288738                                                                                                       ---Execution time = 123s---

Notes:

-Like in the 1st paragraph, the termination criterion of lines 72-73 may be improved:

1-Add DSE 10 GTO 02 after line 73
2-Add 3 STO 10 LBL 02 after line 34

-If the function f has a singularity for x = a and/or x = b , x - a and/or b - x will not calculated with a good precision and the program may fail.

>>> In the following variant, the HP41 uses 1 - x = [ exp ( - sinh x ) } / cosh ( sinh x ) if x > 0

Data Registers: • R00 = Function name ( Register R00 is to be initialized before executing "TSI2" )

                                         R01 = I = R07                          R03 = k                     R10 = ( b - a ) / 2
                                         R02 = h , h/2 , h/4 , ....              R04 to R09: temp      R11 = ( a + b ) / 2
Flags: /
Subroutine:   A program that takes x in X-register ,
                                                   b - x in Y-register ,
                                                   x - a in Z-register and returns f(x) in X-register

01 LBL "TSI2"
02 STO 02
03 RDN
04 STO 10
05 X<>Y
06 ST- 10
07 +
08 2
09 ST/ 10
10 /
11 STO 11
12 E99
13 STO 07
14 SIGN
15 STO 06
16 CLX
17 STO 03
18 RCL 10
19 ENTER
20 R^

21 XEQ IND 00
22 RCL 02
23 *
24 RCL 10
25 *
26 STO 01
27 GTO 01
28 LBL 00
29 2
30 ST/ 01
31 ST/ 02
32 STO 06
33 SIGN
34 CHS
35 STO 03
36 LBL 01
37 RCL 06
38 ST+ 03
39 RCL 02
40 RCL 03

41 *
42 SINH
43 STO 05
44 TANH
45 RCL 10
46 *
47 STO 04
48 RCL 11
49 +
50 RCL 04
51 RCL 10
52 +
53 STO 08
54 RCL 05
55 CHS
56 E^X
57 RCL 05
58 COSH
59 STO 05
60 /

61 RCL 10
62 *
63 STO 09
64 R^
65 XEQ IND 00
66 X<> 04
67 CHS
68 RCL 11
69 +
70 RCL 08
71 RCL 09
72 X<> Z
73 XEQ IND 00
74 RCL 04
75 +
76 RCL 02
77 ST* Y
78 RCL 03
79 *
80 COSH

81 *
82 RCL 10
83 *
84 RCL 05
85 X^2
86 /
87 RCL 01
88 +
89 STO 01
90 LASTX
91 X#Y?
92 GTO 01
93 VIEW X
94 X<> 07
95 X#Y??
96 GTO 00
97 RCL 07
98 END

( 130 bytes / SIZE 012 )

STACK	INPUT	OUTPUT
Z	a	/
Y	b	/
X	h > 0	I

Example: I = §₂⁸ 1 / sqrt(-x^4+20.x^3-150.x^2+500.x-544) dx = §₂⁸ 1 / sqrt((8-x).(x-2).(x^2-10.x+34)) dx

-Load this subroutine in main memory

01 LBL "T"
02 STO 13
03 10
04 -
05 RCL 13
06 *
07 34
08 +
09 *
10 *
11 SQRT
12 1/X
13 RTN
14 END

"T" ASTO 00 and if you choose h = 0.2

     2    ENTER^
     8    ENTER^
    0.2 XEQ "TSI2" >>>>   0.874019185                          The HP-41 displays the successive approximations
                                            0.874019185                                                                                                        ---Execution time = 3m23s---

Note:

-The criterion of lines 68-69 may be improved like with the 1st version:

1-Add DSE 12 GTO 02 after line 92
2-Add 3 STO 12 LBL 02 after line 36

b) Without M-Code Routine

-The same results may be obtained without M-Code:

1st Version:

01 LBL "TSI2"
02 STO 02
03 RDN
04 STO 09
05 X<>Y
06 ST- 09
07 +
08 2
09 ST/ 09
10 /
11 STO 08
12 E99
13 STO 07
14 SIGN
15 STO 06
16 CLX
17 STO 03
18 RCL 08
19 XEQ IND 00
20 RCL 02

21 *
22 RCL 09
23 *
24 STO 01
25 GTO 01
26 LBL 00
27 2
28 ST/ 01
29 ST/ 02
30 STO 06
31 SIGN
32 CHS
33 STO 03
34 LBL 01
35 RCL 06
36 ST+ 03
37 RCL 02
38 RCL 03
39 *
40 E^X-1

41 LASTX
42 CHS
43 E^X-1
44 -
45 STO 05
46 E^X-1
47 STO Y
48 2
49 +
50 /
51 RCL 09
52 *
53 STO 04
54 RCL 08
55 +
56 XEQ IND 00
57 X<> 04
58 CHS
59 RCL 08
60 +

61 XEQ IND 00
62 RCL 04
63 +
64 RCL 02
65 ST* Y
66 RCL 03
67 *
68 E^X
69 ENTER
70 1/X
71 +
72 ST+ X
73 *
74 RCL 09
75 *
76 RCL 05
77 2
78 /
79 E^X
80 ENTER

81 1/X
82 +
83 X^2
84 /
85 RCL 01
86 +
87 STO 01
88 LASTX
89 X#Y?
90 GTO 01
91 VIEW X
92 X<> 07
93 RND
94 X<>Y
95 RND
96 X#Y?
97 GTO 00
98 RCL 07
99 END

( 127 bytes / SIZE 010 )

STACK	INPUT	OUTPUT
Z	a	/
Y	b	/
X	h	I

-Here, the precision is controlled by the display format
-FIX 9 could lead to unuseful loops and FIX 8 is often enough to get the best accuracy.

-Same remark as above about the criterion of lines 89-90

2nd Version:

01 LBL "TSI2"
02 STO 02
03 RDN
04 STO 10
05 X<>Y
06 ST- 10
07 +
08 2
09 ST/ 10
10 /
11 STO 11
12 E99
13 STO 07
14 SIGN
15 STO 06
16 CLX
17 STO 03
18 RCL 10
19 ENTER
20 R^
21 XEQ IND 00
22 RCL 02
23 *
24 RCL 10

25 *
26 STO 01
27 GTO 01
28 LBL 00
29 2
30 ST/ 01
31 ST/ 02
32 STO 06
33 SIGN
34 CHS
35 STO 03
36 LBL 01
37 RCL 06
38 ST+ 03
39 RCL 02
40 RCL 03
41 *
42 E^X-1
43 LASTX
44 CHS
45 E^X-1
46 -
47 STO 05
48 2

49 /
50 X<> 05
51 E^X-1
52 STO Y
53 2
54 +
55 /
56 RCL 10
57 *
58 STO 04
59 RCL 11
60 +
61 RCL 05
62 CHS
63 E^X
64 RCL X
65 RCL 05
66 E^X
67 +
68 2
69 /
70 STO 05
71 /
72 RCL 10

73 *
74 STO 09
75 RCL 04
76 RCL 10
77 +
78 STO 08
79 X<> Z
80 XEQ IND 00
81 X<> 04
82 CHS
83 RCL 11
84 +
85 RCL 08
86 RCL 09
87 X<> Z
88 XEQ IND 00
89 RCL 04
90 +
91 RCL 02
92 ST* Y
93 RCL 03
94 *
95 E^X
96 ENTER

97 1/X
98 +
99 *
100 RCL 10
101 *
102 RCL 05
103 X^2
104 ST+ X
105 /
106 RCL 01
107 +
108 STO 01
109 LASTX
110 X#Y?
111 GTO 01
112 VIEW X
113 X<> 07
114 RND
115 X<>Y
116 RND
117 X#Y?
118 GTO 00
119 RCL 07
120 END

( 120 bytes / SIZE 012 )

STACK	INPUT	OUTPUT
Z	a	/
Y	b	/
X	h > 0	I

-Here again, the precision is controlled by the display format
-FIX 9 could lead to unuseful loops and FIX 8 is often enough to get the best accuracy.

-Same remark as above about the criterion of lines 110-111

3°) Sinh-Sinh Quadrature I = §_-infinity^+infinity f(x) dx

a) With M-Code Routines SINH COSH X#Y??

-Here, we use the change of variable x = Sinh ( Sinh t ) where the factor PI / 2 has again been omitted

Data Registers: • R00 = Function name ( Register R00 is to be initialized before executing "SSI" )

01 LBL "SSI"
02 STO 02
03 E99
04 STO 07
05 SIGN
06 STO 06
07 CLX
08 STO 03
09 XEQ IND 00
10 RCL 02
11 *
12 STO 01

13 GTO 01
14 LBL 00
15 2
16 ST/ 01
17 ST/ 02
18 STO 06
19 SIGN
20 CHS
21 STO 03
22 LBL 01
23 RCL 06
24 ST+ 03

25 RCL 02
26 RCL 03
27 *
28 SINH
29 STO 05
30 SINH
31 STO 04
32 XEQ IND 00
33 X<> 04
34 CHS
35 XEQ IND 00
36 RCL 04

37 +
38 RCL 02
39 ST* Y
40 RCL 03
41 *
42 COSH
43 *
44 RCL 05
45 COSH
46 *
47 RCL 01
48 +

49 STO 01
50 LASTX
51 X#Y?
52 GTO 01
53 VIEW X
54 X<> 07
55 X#Y??
56 GTO 00
57 RCL 07
58 END

( 86 bytes / SIZE 008 )

STACK	INPUT	OUTPUT
X	h	I

Example: Evaluate I = §_-infinity^+infinity exp ( -x^2 ) dx

01 LBL "T"
02 X^2
03 CHS
04 E^X
05 RTN
06 END

"T" ASTO 00 and if you choose h = 0.2

     0.2   XEQ "SSI" >>>>     1.772418685                          The HP-41 displays the successive approximations
                                              1.772453849
                                              1.772453851                                                                                                        ---Execution time = 124s---

Notes:

-The exact result is sqrt(PI)
-As usual, the criterion of lines 51-52 may be improved like this:

1-Add DSE 08 GTO 02 after line 52
2-Add 3 STO 08 LBL 02 after line 22

b) Without M-Code Routine

-Without SINH , COSH and X#Y?? it yields:

01 LBL "SSI"
02 STO 02
03 E99
04 STO 07
05 SIGN
06 STO 06
07 CLX
08 STO 03
09 XEQ IND 00
10 RCL 02
11 *
12 STO 01
13 GTO 01
14 LBL 00
15 2
16 ST/ 01

17 ST/ 02
18 STO 06
19 SIGN
20 CHS
21 STO 03
22 LBL 01
23 RCL 06
24 ST+ 03
25 RCL 02
26 RCL 03
27 *
28 E^X-1
29 LASTX
30 CHS
31 E^X-1
32 -

33 2
34 /
35 STO 05
36 E^X-1
37 LASTX
38 CHS
39 E^X-1
40 -
41 2
42 /
43 STO 04
44 XEQ IND 00
45 X<> 04
46 CHS
47 XEQ IND 00
48 RCL 04

49 +
50 RCL 02
51 ST* Y
52 RCL 03
53 *
54 E^X
55 ENTER
56 1/X
57 +
58 *
59 RCL 05
60 E^X
61 ENTER
62 1/X
63 +
64 *

65 4
66 /
67 RCL 01
68 +
69 STO 01
70 LASTX
71 X#Y?
72 GTO 01
73 VIEW X
74 X<> 07
75 RND
76 X<>Y
77 RND
78 X#Y?
79 GTO 00
80 RCL 07
81 END

( 104 bytes / SIZE 008 )

STACK	INPUT	OUTPUT
X	h	I

-Here, the precision is controlled by the display format
-FIX 9 could lead to unuseful loops and FIX 8 is often enough to get the best accuracy.

-Same remarks as above about the criterion of lines 71-72

4°) Exp-Sinh Quadrature I = §₀^+infinity f(x) dx

a) With M-Code Routines SINH COSH X#Y??

-The change of variable is now: x = Exp ( Sinh t ) where the factor PI / 2 is again omitted

Data Registers: • R00 = Function name ( Register R00 is to be initialized before executing "ESI" )

01 LBL "ESI"
02 STO 02
03 E99
04 STO 07
05 CLX
06 STO 03
07 SIGN
08 STO 06
09 XEQ IND 00
10 RCL 02
11 *
12 STO 01

13 GTO 01
14 LBL 00
15 2
16 ST/ 01
17 ST/ 02
18 STO 06
19 SIGN
20 CHS
21 STO 03
22 LBL 01
23 RCL 06
24 ST+ 03

25 RCL 02
26 RCL 03
27 *
28 SINH
29 E^X
30 STO 04
31 STO 05
32 XEQ IND 00
33 X<> 04
34 1/X
35 XEQ IND 00
36 RCL 05

37 /
38 RCL 04
39 RCL 05
40 *
41 +
42 RCL 02
43 ST* Y
44 RCL 03
45 *
46 COSH
47 *
48 RCL 01

49 +
50 STO 01
51 LASTX
52 X#Y?
53 GTO 01
54 VIEW X
55 X<> 07
56 X#Y??
57 GTO 00
58 RCL 07
59 END

( 85 bytes / SIZE 008 )

STACK	INPUT	OUTPUT
X	h	I

Example: Evaluate I = §₀^+infinity exp ( -x^2 ) dx

01 LBL "T"
02 X^2
03 CHS
04 E^X
05 RTN
06 END

"T" ASTO 00 and if you choose h = 0.2

     0.2   XEQ "ESI" >>>>     0.886226890                          The HP-41 displays the successive approximations
                                              0.886226926
                                              0.886226926                                                                                                        ---Execution time = 3m48s---

Notes:

-The exact result is sqrt(PI) / 2
-As usual, the criterion of lines 52-53 may be improved like this:

1-Add DSE 08 GTO 02 after line 53
2-Add 3 STO 08 LBL 02 after line 22

b) Without M-Code Routine

-If you want to avoid M-Code:

01 LBL "ESI"
02 STO 02
03 E99
04 STO 07
05 CLX
06 STO 03
07 SIGN
08 STO 06
09 XEQ IND 00
10 RCL 02
11 *
12 STO 01
13 GTO 01
14 LBL 00
15 2

16 ST/ 01
17 ST/ 02
18 STO 06
19 SIGN
20 CHS
21 STO 03
22 LBL 01
23 RCL 06
24 ST+ 03
25 RCL 02
26 RCL 03
27 *
28 E^X-1
29 LASTX
30 CHS

31 E^X-1
32 -
33 2
34 /
35 E^X
36 STO 04
37 STO 05
38 XEQ IND 00
39 X<> 04
40 1/X
41 XEQ IND 00
42 RCL 05
43 /
44 RCL 04
45 RCL 05

46 *
47 +
48 RCL 02
49 ST* Y
50 RCL 03
51 *
52 E^X
53 ENTER
54 1/X
55 +
56 *
57 2
58 /
59 RCL 01
60 +

61 STO 01
62 LASTX
63 X#Y?
64 GTO 01
65 VIEW X
66 X<> 07
67 RND
68 X<>Y
69 RND
70 X#Y?
71 GTO 00
72 RCL 07
73 END

( 96 bytes / SIZE 008 )

STACK	INPUT	OUTPUT
X	h	I

-Here, the precision is controlled by the display format
-FIX 9 could lead to unuseful loops and FIX 8 is often enough to get the best accuracy.

-Same remarks as above about the criterion of lines 63-64

References:

[1] Hidetosi TAKAHASI and Masatake MORI - "Double Exponential Formulas for Numerical Integration"
[2] David H. Bailey - "Tanh-Sinh High-Precision Quadrature"

hp41programs

Numerical Integration: Double Exponential Formulas for the HP41