# 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  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,  ....
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 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 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 , ....... , 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:                                           ( 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  §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:                                         ( 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  §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:                                              ( 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 ( -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:                                           ( Registers R00 thru R03 are to be initialized before executing "IG" )

•  R00:  f name     R04:    §ab 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 =  §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  75  X#Y?  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

-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:  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     76  X#Y?    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

-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    §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:                                                 ( 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    §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:                                                 ( 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    §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
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