The Museum of HP Calculators

# Romberg's method for the HP-41

This program is by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

1°) Romberg-Programs
a) Program#1  ( improved )
b) Program#2  ( new )
2°) Arc Length of a Curve
3°) Area of a Surface of Revolution
4°) Area of a 3-dimensional Surface ( new )
5°) Integrals
a) Program#1
b) Program#2  ( new )
c) Program#3  ( new )
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  In  tends to I  as  n tends to infinity  and  that the  "errors"  I - In  are nearly proportional to 1/nk
If errors were exactly proportional to 1/nk :      I = In + A/nk  = I2n + A/(2n)k ,      whence  I = ( 2k I2n - In ) / ( 2k - 1 )
When errors are only nearly proportional to 1/nk , the errors in this last formula are proportional to 1/nk+2

- "ROM" calculates   I1 , I2 , I4 , I8 ,  .....  ( contimually doubling n )   The formula above is applied to each pair of consecutive In
the results are        J1 , J2 , J4 ,  ....   and "ROM"  applies recursively the previous formula
producing               K1,K2,  ....
....,  ....   until two rounded consecutive approximations are equal.  Thus, accuracy is determined by the display setting.

-This program calls a subroutine which calculates In  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:  2k                 R24:  counter
R23:  2k , 22k , ...   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
44  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 , ....... , 2m  and with n = 1 , 2 , 4 , 8 , ....... , 2m  are equal.
-This trick generally avoids an unuseful extra-computation with  n = 2m+1
but unfortunately, it may sometimes cause a premature stop.

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

b) Program#2

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

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 = §ab  ( 1 + (y')2 )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:

•  R00:  f name     R04:  L                 Registers R00 thru R03 are to be initialized before executing "LNG"
•  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 "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  §ab  y ( 1 + (y')2 )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:

•  R00:  f name     R04:  S                 Registers R00 thru R03 are to be initialized before executing "SRV"
•  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  §ab §cd   ( 1 + fx2 + fy2  )1/2 dx dy

where  fx = f/x  and   fy = 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:

•  R00:  f name     •  R04:    c              Registers R00 thru R05 are to be initialized before executing "SKS"
•  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.

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                     a three-byte GTO
133  RCL 07
134  ST+ 09
135  DSE 17
136  GTO 01                     a three-byte GTO
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 ( -x2-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   §ab f(x) dx

a) Program#1

-The following program uses the very simple "midpoint" formula:   §ab 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:

•  R00:  f name     R04:    §ab f(x) dx              Registers R00 thru R03 are to be initialized before executing "IG"
•  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 =  §12  xx  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:  22 , 42 , 82 ,  ....
R12:  (b-a)/2n   R16:  2n                  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                                     the accuracy is controlled by the display setting
75  X#Y?                                    the iterations stop when 2 rounded successive approximations are equal.
76  GTO 01
77  RCL IND 18
78  END

( 114 bytes / SIZE 021+? )

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

Example:    Evaluate   §13  x ( 1+x3 )1/2  dx

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

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:  3n                  R16:  counter and 32 , 92 , 272 ,  ....
R11:  a                R14:  sums              R17:  counter
R12:  (b-a)/3n     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                                     the accuracy is controlled by the display setting
76  X#Y?                                    the iterations stop when 2 rounded successive approximations are equal.
77  GTO 01
78  RCL IND 17
79  END

( 117 bytes / SIZE 020+? )

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

Example:    Evaluate   §13  x ( 1+x3 )1/2  dx

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

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 loose the previous function evaluations,
but step trippling often leads to long execution times.

6°) Double Integrals    §ab §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:

•  R00:  f name     •  R04:    u name              Registers R00 thru R05 are to be initialized before executing "DIN"
•  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    §23 §xx^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   §ab  §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:

•  R00:  f name     •  R04:    u name              Registers R00 thru R05 are to be initialized before executing "TIN"
•  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    §01 §0x §-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
---- n = 2   ------------    0.236284830
---- n = 4   ------------    0.246481100
---- n = 8   ------------    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