hp41programs

Romberg Romberg's method for the HP-41

Overview

1°) Romberg-Programs

a) Program#1
b) Program#2

2°) Arc Length of a Curve
3°) Area of a Surface of Revolution
4°) Area of a 3-dimensional Surface
5°) Integrals

   a) Program#1
   b) Program#2
   c) Program#3

6°) Double Integrals
7°) Triple Integrals

1°) Romberg-Programs

a) Program#1

-The Romberg's method is an "extrapolation to the limit" algorithm which can be applied to many problems:

-Suppose that a sequence I_n tends to I as n tends to infinity and that the "errors" I - I_n are nearly proportional to 1/n^k
If errors were exactly proportional to 1/n^k : I = I_n + A/n^k = I_2n + A/(2n)^k , whence I = ( 2^k I_2n - I_n ) / ( 2^k - 1 )
When errors are only nearly proportional to 1/n^k , the errors in this last formula are proportional to 1/n^k+2

- "ROM" calculates   I₁ , I₂ , I₄ , I₈ , ..... ( contimually doubling n )   The formula above is applied to each pair of consecutive I_n
         the results are        J₁ , J₂ , J₄ , ....   and "ROM" applies recursively the previous formula
              producing               K₁,K₂, ....
                    and so on               ...., ....   until two rounded consecutive approximations are equal. Thus, accuracy is determined by the display setting.

-This program calls a subroutine which calculates I_n for a given n-value. n is stored in register R03 and R03 must be unchanged by the subroutine.
-The name of this subroutine ( global label, 6 characters or less ) is to be stored into register R20
-The integer k ( the "order" of the method ) is to be stored into register R21

Data Registers:

R03: n ( change line 03 and line 10 if you store n in another register )
• R20: the name of the subroutine ( for instance an integrator ... )
• R21: k = the order of the method used by the subroutine

            ( k = 2 for the trapezoidal rule,
              k = 4 for Simpson's rule,
              k = 6 for the 3-point Gauss-Legendre formula ...

>>>> However, even Gaussian integration can behave like a first order method with singular integrals ! )

R22: 2^k R24: counter
R23: 2^k , 2^2k , ... R25 , R26: .... successive approximations.

Flags: /
Subroutines: It depends on the problem you have to solved

01 LBL "ROM"
02 1
03 STO 03
04 25
05 STO 24
06 XEQ IND 20
07 STO 25
08 LBL 01
09 2

10 ST* 03
11 RCL 21
12 Y^X
13 STO 22
14 STO 23
15 RCL 24
16 INT
17 E3
18 /

19 25
20 +
21 STO 24
22 XEQ IND 20
23 LBL 02
24 ENTER^
25 ENTER^
26 X<> IND 24
27 -

28 RCL 23
29 RCL 22
30 ST* 23
31 SIGN
32 -
33 /
34 +
35 ISG 24
36 GTO 02

37 STO IND 24
38 VIEW X
39 RND
40 X<>Y
41 RND
42 X#Y?
43 GTO 01
44 RCL IND 24
45 END

( 76 bytes SIZE 026+ ? )

STACK	INPUTS	OUTPUTS
X	/	Limit

N.B:

-The subroutines called by "ROM" are supposed to keep registers R20 R21 ... unchanged.
-Otherwise, replace for instance these registers by R30 R31 ... in the "ROM" listing, and replace line 04 and line 19 by 35.
-The Lagrange interpolation formula can also be used to obtain the same results ( cf the last example at the bottom of this page )

-In fact, "ROM" stops when the rounded approximations obtained with n = 2 , 4 , 8 , ....... , 2^m and with n = 1 , 2 , 4 , 8 , ....... , 2^m are equal.
-This trick generally avoids an unuseful extra-computation with n = 2^m+1
but unfortunately, it may sometimescause a premature stop.

-The variant hereafter really compares 2 successive approximations obtained with n = 1 , 2 , 4 , 8 , ....... , 2^m and n = 1 , 2 , 4 , 8 , ....... , 2^m+1

b) Program#2

Data Registers:

R22: 1 , 1/2 , 1/4 , 1/8 , ...... R24 , R25: .... successive approximations.
R23: counter

Flags: /
Subroutines: It depends on the problem you have to solved

01 LBL "ROM2"
02 1
03 STO 03
04 24
05 STO 23
06 XEQ IND 20
07 STO 24
08 LBL 01
09 2
10 ST* 03

11 SIGN
12 STO 22
13 RCL 23
14 INT
15 E3
16 /
17 24
18 +
19 STO 23
20 XEQ IND 20

21 LBL 02
22 ENTER^
23 X<> IND 23
24 STO Z
25 -
26 2
27 ST/ 22
28 CLX
29 RCL 22
30 RCL 21

31 Y^X
32 1
33 -
34 /
35 -
36 ISG 23
37 GTO 02
38 STO IND 23
39 VIEW X
40 RND

41 X<>Y
42 RND
43 X#Y?
44 GTO 01
45 RCL IND 23
46 END

( 77 bytes / SIZE 025+? )

STACK	INPUTS	OUTPUTS
X	/	Limit

N.B:

-The termination criterion is of course more rigorous, but the execution time is longer.
-Moreover, there is now a risk of infinite loop, especially if the display setting is FIX 9 or SCI 9.
-The following modification is a partial remedy:

Add LASTX 9 / + after line 39
and ENTER^ after line 37

-Now let's apply the Romberg method to solve a few problems:
( all the following formulas are second-order methods: R21 = 2 )

2°) Arc length of a curve

-The arc length of the curve y = f(x) ( a < x < b ) is given by L = §_a^b ( 1 + (y')² )^1/2 dx
"LNG" doesn't use this formula and avoids the calculation of dy/dx . It simply applies Pythagoras' theorem.
-This is a second order method: k = 2 STO 21

Data Registers: ( Registers R00 thru R03 are to be initialized before executing "LNG" )

   • R00: f name     R04: L                             "ROM" initializes R03 automatically
   • R01: a             R05: counter
   • R02: b             R06: (b-a)/n
   • R03: n             R07- R08: temp.

Flags: /
Subroutine: a program which takes x in X-register and calculates f(x) in X-register.

01 LBL "LNG"
02 RCL 02
03 RCL 01
04 STO 07
05 -
06 RCL 03
07 STO 05

08 /
09 STO 06
10 ST+ 07
11 RCL 01
12 XEQ IND 00
13 STO 08
14 CLX

15 STO 04
16 LBL 01
17 RCL 07
18 XEQ IND 00
19 ENTER^
20 X<> 08
21 -

22 X^2
23 RCL 06
24 ST+ 07
25 X^2
26 +
27 SQRT
28 ST+ 04

29 DSE 05
30 GTO 01
31 RCL 04
32 END

( 48 bytes / SIZE 009 )

STACK	INPUTS	OUTPUTS
X	/	Arc Length

Example:

-Calculate the arc length of the curve y = ln x 1 < x < 3

    LBL "T"
    LN
    RTN

"LNG"   ASTO 20    2    STO 21
"T"       ASTO 00    1    STO 01   3   STO 02
FIX 9    XEQ   "ROM"     the following approximations are displayed:

    2.300459681
    2.301931038
    2.301986835
    2.301987537
    2.301987533     ( 76 seconds )    The exact result is 2.301987535 ( rounded to 9 decimals )

Notes:

1- If you execute "ROM" in FIX 5 the program gives 2.30199
    At this step, if you want a better accuracy, you don't have to start all over again: simply modify the display format ( for instance FIX 9 ) and press XEQ 01
2-If the equation is given in polar coordinates, replace lines 21-24 by    ST+ Z - X^2 X<>Y 2   /   RCL 06 ST+ 07 *    and add    ENTER^   after line 19
   ( this is a second order method too: k = 2 STO 21 )
3-This method can be generalized to parametric equations ...

3°) Area of a surface of revolution

-The rotation of the curve y = f(x) ( a < x < b ) around x-axis generates a surface of revolution given by S = 2.pi §_a^b y ( 1 + (y')² )^1/2 dx
"SRV" doesn't use this formula and avoids the calculation of dy/dx : the area of a truncated cone is used instead.

-It is also a second order method : k = 2 STO 21

Data Registers: ( Registers R00 thru R03 are to be initialized before executing "SRV" )

   • R00: f name     R04: S
   • R01: a             R05: counter                   "ROM" initializes R03 automatically
   • R02: b             R06: (b-a)/n
   • R03: n             R07 R08 temp.

Flags: /
Subroutine: a program which takes x in X-register and calculates f(x) in X-register.

01 LBL "SRV"
02 RCL 02
03 RCL 01
04 STO 07
05 -
06 RCL 03
07 STO 05
08 /

09 STO 06
10 ST+ 07
11 RCL 01
12 XEQ IND 00
13 STO 08
14 CLX
15 STO 04
16 LBL 01

17 RCL 07
18 XEQ IND 00
19 ENTER^
20 ENTER^
21 X<> 08
22 ST+ Z
23 -
24 X^2

25 RCL 06
26 ST+ 07
27 X^2
28 +
29 SQRT
30 *
31 ST+ 04
32 DSE 05

33 GTO 01
34 PI
35 ST* 04
36 RCL 04
37 END

( 55 bytes / SIZE 009 )

STACK	INPUTS	OUTPUTS
X	/	Area

Example:

-Evaluate the area of the surface of revolution generated by the rotation of the curve y = sin x ( 0 < x < pi ) around the x-axis.

LBL "T"
SIN
RTN

"SRV"   ASTO 20    2    STO 21
"T"      ASTO 00    0    STO 01    PI    STO 02    set the HP-41 in RAD mode ( XEQ RAD )
FIX 9   XEQ   "ROM"     the following approximations are displayed:

15.59985804
14.26493243
14.42990383
14.42356623
14.42359884
14.42359950 ( 177 seconds )

4°) Area of a 3-Dimensional Surface

-This routine computes the area of a 3D-surface defined by: z = f(x,y) a < x < b , c < y < d

-The result could be obtained by the double integral §_a^b §_c^d ( 1 + f_x² + f_y² )^1/2 dx dy

where f_x = ¶f/¶x and f_y = ¶f/¶y are the partial derivatives with respect to x and y respectively.

-But "SKS" avoids the calculus of the partial derivatives:
-The intervals [a,b] and [c,d] are divided into n parts, and the approximate area is the sum of the areas of triangles obtained by Heron's formula.
-This is a second order method.

Data Registers: ( Registers R00 thru R05 are to be initialized before executing "SKS" )

   • R00: f name     • R04:    c
   • R01: a             • R05:    d                          "ROM" initializes R03 automatically
   • R02: b                R06:   Area
   • R03: n                R07 thru R18 temp.

Flags: /
Subroutine: A program which takes x in X-register and y in Y-register and calculates f(x,y) in X-register.

-Lines 132-136 are three-byte GTOs

01 LBL "SKS"
02 RCL 02
03 RCL 01
04 STO 09
05 -
06 STO 07
07 RCL 05
08 RCL 04
09 -
10 RCL 03
11 ST/ 07
12 STO 17
13 /
14 STO 08
15 CLX
16 STO 06
17 LBL 01
18 RCL 03
19 STO 18
20 RCL 04
21 STO 10
22 RCL 09
23 XEQ IND 00
24 STO 11
25 RCL 10
26 RCL 09
27 RCL 07
28 +

29 XEQ IND 00
30 STO 12
31 LBL 02
32 RCL 08
33 ST+ 10
34 RCL 10
35 RCL 09
36 RCL 07
37 +
38 XEQ IND 00
39 STO 14
40 RCL 10
41 RCL 09
42 XEQ IND 00
43 STO 13
44 RCL 12
45 -
46 X^2
47 RCL 07
48 X^2
49 +
50 RCL 08
51 X^2
52 +
53 SQRT
54 STO 15
55 RCL 12
56 RCL 11

57 -
58 X^2
59 RCL 07
60 X^2
61 +
62 SQRT
63 STO 16
64 +
65 RCL 13
66 RCL 11
67 -
68 X^2
69 RCL 08
70 X^2
71 +
72 SQRT
73 STO Z
74 +
75 2
76 /
77 ENTER^
78 ENTER^
79 R^
80 -
81 *
82 X<>Y
83 RCL 15
84 -

85 *
86 X<>Y
87 RCL 16
88 -
89 *
90 SQRT
91 ST+ 06
92 RCL 14
93 RCL 12
94 -
95 X^2
96 RCL 08
97 X^2
98 +
99 SQRT
100 STO 16
101 RCL 14
102 STO 12
103 RCL 13
104 STO 11
105 -
106 X^2
107 RCL 07
108 X^2
109 +
110 SQRT
111 STO Z
112 +

113 RCL 15
114 +
115 2
116 /
117 ENTER^
118 ENTER^
119 ST- 15
120 R^
121 -
122 *
123 RCL 15
124 *
125 RCL 16
126 R^
127 -
128 *
129 SQRT
130 ST+ 06
131 DSE 18
132 GTO 02
133 RCL 07
134 ST+ 09
135 DSE 17
136 GTO 01
137 RCL 06
138 END

( 170 bytes / SIZE 019 )

STACK	INPUTS	OUTPUTS
X	/	Area

Example: Calculate the area of the surface defined by z = exp ( -x²-y ) 0 < x < 2 , 0 < y < 3

LBL "Z"
X^2
+
CHS
E^X
RTN

"Z"   ASTO 00
   0    STO 01 STO 04
   2    STO 02
   3    STO 05

"SKS" ASTO 20
-This is a second order method: 2 STO 21

-For instance FIX 3 XEQ "ROM" the following approximations are displayed:

6.262
6.284
6.284 ( in 6mn46s )

-The last approximation is actually 6.283584 in X-register
whereas the exact result is 6.2833787 rounded to 7 decimals.

5°) Integrals §_a^b f(x) dx

a) Program#1

-The following program uses the very simple "midpoint" formula: §_a^b f(x) dx = (b-a).f((a+b)/2) It is a second order method ( k = 2 STO 21 )

-It does not have the refinements of the PPC ROM program "IG" ( cf PPC-ROM user's manual page 220-225 )
therefore, the convergence may be disappointing for singular integrals,
but it's a good example of how the Romberg's method can produce superb acceleration.
-Furthermore, it's probably one of the shortest integration program for the HP-41.

Data Registers: ( Registers R00 thru R03 are to be initialized before executing "IG" )

   • R00: f name     R04:    §_a^b f(x) dx
   • R01: a             R05: counter                   "ROM" initializes R03 automatically
   • R02: b             R06: (b-a)/n
   • R03: n             R07   temp.

Flags: /
Subroutine: a program which takes x in X-register and calculates f(x) in X-register.

01 LBL "IG"
02 RCL 02
03 RCL 01
04 STO 07
05 -

06 RCL 03
07 STO 05
08 /
09 STO 06
10 2

11 /
12 ST+ 07
13 CLX
14 STO 04
15 LBL 01

16 RCL 07
17 XEQ IND 00
18 ST+ 04
19 RCL 06
20 ST+ 07

21 DSE 05
22 GTO 01
23 ST* 04
24 RCL 04
25 END

( 40 bytes / SIZE 008 )

STACK	INPUTS	OUTPUTS
X	/	Integral

Example:

-Evaluate I = §₁² x^x dx to 7 decimals

LBL "T"
ENTER^
Y^X
RTN

"IG"     ASTO 20   2   STO 21
"T"      ASTO 00   1   STO 01 2 STO 02
FIX 7   XEQ "ROM"    the following values are displayed:

2.0438805
2.0503885
2.0504461
2.0504462 ( in 40 seconds )

-In fact, the result is correct to almost 9 decimals: I = 2.050446234 ( the last digit should be a 5 )

b) Program#2

-This routine uses the Romberg method with the trapezoidal rule.
-So, f(a) and f(b) must be finite.
-The previous sums are used the next iterations.
-The Romberg extrapolation is performed inside the program itself: it doesn't use "ROM"

Data Registers: R00 thru R09 are unused ( Register R10 is to be initialized before executing "IG2" )

   • R10: f name     R13-R14: temp
      R11: a              R15: counter           R17: 2² , 4² , 8² , ....
      R12: (b-a)/2ⁿ   R16: 2ⁿ                  R18: counter                           R19 , R20 , .... = successive approximations

Flags: /
Subroutine: a program which calculates f(x) assuming x is in X-register upon entry.

01 LBL "IG2"
02 STO 12
03 X<>Y
04 STO 11
05 X<>Y
06 XEQ IND 10
07 STO 14
08 RCL 11
09 ST- 12
10 XEQ IND 10
11 RCL 14
12 +
13 2
14 /
15 STO 14
16 RCL 12

17 *
18 STO 19
19 19
20 STO 18
21 SIGN
22 STO 16
23 LBL 01
24 RCL 16
25 STO 15
26 RCL 11
27 RCL 12
28 2
29 /
30 +
31 LBL 02
32 STO 13

33 XEQ IND 10
34 ST+ 14
35 RCL 12
36 RCL 13
37 +
38 DSE 15
39 GTO 02
40 2
41 ST/ 12
42 ST* 16
43 X^2
44 STO 17
45 RCL 18
46 INT
47 E3
48 /

49 19
50 +
51 STO 18
52 RCL 12
53 RCL 14
54 *
55 LBL 03
56 ENTER^
57 ENTER^
58 X<> IND 18
59 STO 13
60 -
61 RCL 17
62 4
63 ST* 17
64 SIGN

65 -
66 /
67 +
68 ISG 18
69 GTO 03
70 STO IND 18
71 VIEW X
72 RND
73 RCL 13
74 RND
75 X#Y?
76 GTO 01
77 RCL IND 18
78 END

( 114 bytes / SIZE 021+? )

STACK	INPUTS	OUTPUTS
Y	a	/
X	b	§_a^b f(x) dx

Example: Evaluate §₁³ x ( 1+x³ )^1/2 dx

LBL "T"      X^2           SIGN     *                       "T"   ASTO 10
ENTER^     *               +            RTN
ENTER^    ENTER^    SQRT

-The accuracy is controlled by the display setting ( line 74 )

FIX 7    1   ENTER^
               3   XEQ "IG2"   >>>>   13.7629072
                                                     13.7691196
                                                     13.7693295
                                                     13.7693320
                                                     13.7693320     ( in 36 seconds )

-The last value is in fact 13.76933202 all the decimals are correct.

c) Program#3

-Like "IG" , "IG3" uses the midpoint formula,
but in order to use the previous evaluations of the function,
the number of steps is trippled with each iteration.

Data Registers: R00 thru R09 are unused ( Register R10 is to be initialized before executing "IG3" )

   • R10: f name       R13: 3ⁿ                  R16: counter and 3² , 9² , 27² , ....
      R11: a               R14: sums              R17: counter
      R12: (b-a)/3ⁿ     R15: x                    R18, R19 , .... = successive approximations

Flags: /
Subroutine: a program which calculates f(x) assuming x is in X-register upon entry.

01 LBL "IG3"
02 STO 12
03 X<>Y
04 STO 11
05 ST- 12
06 +
07 2
08 /
09 XEQ IND 10
10 STO 14
11 RCL 12
12 *
13 STO 18
14 18
15 STO 17
16 SIGN

17 STO 13
18 LBL 01
19 RCL 13
20 STO 16
21 3
22 ST/ 12
23 ST* 13
24 RCL 11
25 RCL 12
26 2
27 /
28 +
29 LBL 02
30 STO 15
31 XEQ IND 10
32 ST+ 14

33 RCL 12
34 ST+ X
35 ST+ 15
36 RCL 15
37 XEQ IND 10
38 ST+ 14
39 RCL 12
40 RCL 15
41 +
42 DSE 16
43 GTO 02
44 9
45 STO 16
46 RCL 17
47 INT
48 E3

49 /
50 18
51 +
52 STO 17
53 RCL 12
54 RCL 14
55 *
56 LBL 03
57 ENTER^
58 ENTER^
59 X<> IND 17
60 STO 15
61 -
62 RCL 16
63 9
64 ST* 16

65 SIGN
66 -
67 /
68 +
69 ISG 17
70 GTO 03
71 STO IND 17
72 VIEW X
73 RND
74 RCL 15
75 RND
76 X#Y?
77 GTO 01
78 RCL IND 17
79 END

( 117 bytes / SIZE 020+? )

STACK	INPUTS	OUTPUTS
Y	a	/
X	b	§_a^b f(x) dx

Example: Evaluate §₁³ x ( 1+x³ )^1/2 dx

LBL "T"      X^2           SIGN     *                       "T"   ASTO 10
ENTER^     *               +            RTN
ENTER^    ENTER^    SQRT

-The accuracy is controlled by the display setting ( lines 73-75 )

FIX 7    1   ENTER^
               3   XEQ "IG3"   >>>>   13.7718432
                                                     13.7693525
                                                     13.7693320
                                                     13.7693320     ( in 64 seconds )

-The last value is 13.76933201 ( the last decimal should be a "2" )
-Actually, the next to last result was exact: 13.76933202 ( in R15 and in L-register )

-Unlike "IG" , "IG3" doesn't lose the previous function evaluations, but step trippling often leads to long execution times.

6°) Double Integrals §_a^b §_u(x)^v(x) f(x,y) dx dy

-The mid-point formula is used in both x and y-axis.
( a second order method : k = 2 STO 21 )

Data Registers: ( Registers R00 thru R05 are to be initialized before executing "DIN" )

   • R00: f name     • R04:    u name
   • R01: a             • R05:    v name                    "ROM" initializes R03 automatically
   • R02: b                R06 = §§
   • R03: n                R07 thru R13: temp.

Flags: /
Subroutines:

   A program which takes x in X-register and y in Y-register and calculates f(x,y) in X-register.
   A program which takes x in X-register and calculates u(x) in X-register.
   A program which takes x in X-register and calculates v(x) in X-register.

01 LBL "DIN"
02 RCL 02
03 RCL 01
04 STO 09
05 -
06 RCL 03
07 STO 07
08 /
09 STO 08
10 2

11 /
12 ST+ 09
13 CLX
14 STO 06
15 LBL 01
16 RCL 09
17 XEQ IND 04
18 STO 10
19 RCL 09
20 XEQ IND 05

21 RCL 10
22 -
23 RCL 03
24 STO 11
25 /
26 STO 12
27 2
28 /
29 ST+ 10
30 CLX

31 STO 13
32 LBL 02
33 RCL 10
34 RCL 09
35 XEQ IND 00
36 ST+ 13
37 RCL 12
38 ST+ 10
39 DSE 11
40 GTO 02

41 RCL 13
42 *
43 ST+ 06
44 RCL 08
45 ST+ 09
46 DSE 07
47 GTO 01
48 ST* 06
49 RCL 06
50 END

( 72 bytes / SIZE 014 )

STACK	INPUTS	OUTPUTS
X	/	Integral

Example:

-Calculate §₂³ §_x^{x^2} ( 1 + xy )^1/2 dx dy to 7 decimals

LBL "T"
*
1
+
SQRT
RTN
LBL "U"
RTN
LBL "V"
X^2
RTN

"DIN"   ASTO 20    2    STO 21
"T"      ASTO 00    2    STO 01 3 STO 02
"U"     ASTO 04   "V" ASTO 05
FIX 7   XEQ "ROM"    the HP-41 displays:

13.7800635
13.7747491
13.7746576
13.7746565 ( 263 seconds )

7°) Triple Integrals §_a^b §_u(x)^v(x) §_w(x,y)^t(x,y) f(x,y,z) dx dy dz

-The mid-point formula is again used ( in x, y and z-axis ).
( a second order method : k = 2 STO 21 )

Data Registers: ( Registers R00 thru R07 are to be initialized before executing "TIN" )

       • R00: f name     • R04:    u name
       • R01: a             • R05:    v name                 "ROM" initializes R03 automatically
       • R02: b             • R06:   w name
       • R03: n             • R07    t name                    R08 = §§§ ,    R09 thru R19 temp.

Flags: /
Subroutines:

A program which takes x in X-register, y in Y-register and z in Z-register and calculates f(x,y,z) in X-register.
A program which takes x in X-register and calculates u(x) in X-register.
A program which takes x in X-register and calculates v(x) in X-register.
A program which takes x in X-register and y in Y-register and calculates w(x,y) in X-register.
A program which takes x in X-register and y in Y-register and calculates t(x,y) in X-register.

01 LBL "TIN"
02 RCL 02
03 RCL 01
04 STO 11
05 -
06 RCL 03
07 STO 09
08 /
09 STO 10
10 2
11 /
12 ST+ 11
13 CLX
14 STO 08
15 LBL 01
16 RCL 11

17 XEQ IND 04
18 STO 12
19 RCL 11
20 XEQ IND 05
21 RCL 12
22 -
23 RCL 03
24 STO 13
25 /
26 STO 14
27 2
28 /
29 ST+ 12
30 CLX
31 STO 15
32 LBL 02

33 RCL 12
34 RCL 11
35 XEQ IND 06
36 STO 16
37 RCL 12
38 RCL 11
39 XEQ IND 07
40 RCL 16
41 -
42 RCL 03
43 STO 17
44 /
45 STO 18
46 2
47 /
48 ST+ 16

49 CLX
50 STO 19
51 LBL 03
52 RCL 16
53 RCL 12
54 RCL 11
55 XEQ IND 00
56 ST+ 19
57 RCL 18
58 ST+ 16
59 DSE 17
60 GTO 03
61 RCL 19
62 *
63 ST+ 15
64 RCL 14

65 ST+ 12
66 DSE 13
67 GTO 02
68 RCL 15
69 *
70 ST+ 08
71 RCL 10
72 ST+ 11
73 DSE 09
74 GTO 01
75 ST* 08
76 RCL 08
77 END

( 114 bytes / SIZE 020 )

STACK	INPUTS	OUTPUTS
X	/	Integral

Example: Evaluate §₀¹ §₀^x §_-x-y^-y 1 / ( 1 + x + y + z ) dx dy dz to 5 places

LBL "P"
+
+
1
+
1/X
RTN
LBL "U"
CLX
RTN
LBL "V"
RTN
LBL "W"
+
CHS
RTN
LBL "T"
X<>Y
CHS
RTN

"TIN"   ASTO 20      2    STO 21
   "P"     ASTO 00      0    STO 01     1 STO 02
   "U"     ASTO 04     V ASTO 05
   "W"    ASTO 06     T   ASTO 07
FIX 5    XEQ "ROM"    the calculator displays:

0.24838
0.24998
0.25000 ( 9 minutes ) ( exact value = 1/4 )

The last approximation is actually 0.249999724 and was obtained from the following results:

with n = 1   "TIN" yields    0.200000000
with n = 2   "TIN" yields    0.236284830
with n = 4   "TIN" yields    0.246481100
with n = 8   "TIN" yields    0.249114231

-If you use these 4 numbers as explained in "Lagrange Interpolation Formula for the HP-41" (example2) you'll get 0.249999724 too

-With n = 16 , the execution time increases very much , but you could perhaps XEQ "TIN" with n = 10 ( STO 03 ): the result is 0.249432635
and then, the Lagrange polynomial gives a slightly better accuracy: with these 5 numbers, you'll find: L(0) = 0.249999997

-"ROM" has the advantage to be an automatic process, but execution time may become a problem because n is continually doubled.
- In the Lagrange approach, you have the opportunity to choose the n-values by yourself in a possibly more economical way.

References:

[1]   "Numerical Analysis" - F. Scheid - Mc Graw Hill   ISBN 07-055197-9
[2]   "PPC-ROM user's manual" ( page 220 ..... )
[3]   "Extended Functions made easy" - Keith Jarett ( synthetix )
[4]   "Numerical Recipes in C++" - Press , Vetterling , Teukolsky , Flannery - Cambridge University Press - ISBN 0-521-75033-4