This routine uses the Romberg algorithm to calculate a numerical approximation of the definite integral of a function. The routine is iterative in that increasingly accurate approximations are calculated until two rounded consecutive approximations are equal. The routine is automatic in that no stepsize information has to be provided. The desired accuracy of the final approximation is determined by the display setting. The consecutive approximations will be shown. The code originates from the PPC ROM Module [1] released by the international group People Programming Computers on December 1981.\(\)
Instructions
 SIZE 030 is the recommended minimum. A few integrals may require a larger size.
 Select display mode. The display setting will control the accuracy of the final approximation. In general a display mode of SCO n will return a value correctly routed to n+1 significant digits. Larger values of n will cause the program to run longer so it is best to select the minimum value of n that is acceptable. The use of SCO or ENG display modes are generally preferable to the FIX mode.
 If a printer is connected, the approximations will be printed.
 Specify the integrand. The integrand represented by the function ƒ(x) must be programmed as a subroutine in program memory which starts with a global label name and ends with a RTN or END instruction. This label name should be six or less characters and will be stored in R10. The input x and the output ƒ(x) are both assumed to be in the X register. Since global label search begins at the bottom of memory, when the ƒ(x) subroutine is in RAM, it should be located at the bottom of RAM. The ƒ(x) program should not use registers R10R29 or use flags F09 and F10.
 Alphastore the global label name from step 4 (six or less alpha characters) in R10.
 Enter limits of integration. The lower and upper limits of integration, a and b, respectively, are to be keyed in as “a ENTER b“, so that a is keyed into the Y register and b is keyed in the X register.
 Execute “IG” program. Key XEQ “IG”. Program execution will commence, and if F10 is set, the consecutive approximations will be displayed. If a printer is connected the approximations will be printed. The final approximation will be left in the Xregister when the program ends.
Method
The method underlying the program is due to Romberg [2] and is essentially an application of Richardson’s extrapolation procedure to the EulerMacLaurin sum formula. Romberg was first to describe the method in recursive form. Commencing with improved midpoint rule estimates, the continued application of extrapolation to the limit produces a lower triangular matrix. Although the columns of the matrix converge to the solution, the diagonal elements converge asymptotically faster than any geometric series, or, super linearly. Program shutoff occurs when two rounded consecutive diagonal elements are equal. The diagonal elements M(k,k) are the values displayed.
Assuming that the number of divisions of the interval of integration is increased by improved midpoint rule estimates, one would assume that convergence could occur when we made the number of intervals high enough. However, at some point, round off errors eventually dominate and our effective accuracy decreases. The Romberg method allows the simulation of a high number of subintervals or divisions which decreases the error, without actually increasing the number of subintervals. This process is called extrapolation to the limit.
Romberg integration is iterative, automatic, and nonadaptive. It is iterative in that it produces increasingly accurate estimates of the solution until the convergence criterion is satisfied. It is automatic in that the number of function evaluations depends upon the behavior of that function over the interval of integration. It is nonadaptive in that the function evaluations occur at a fixed set of points, independent of the function.
As the Romberg method successively halves the interval of integration to produce improved midpoint rule estimates, it uses all previously computed functional evaluations at each stage. The retention of all previous calculated function evaluations is a significant aspect of the Romberg algorithm. As the function is evaluated at the center of each interval, the end points of the intervals are not used as sample points. Hence, the endpoints of the interval of integration, a and b are also not used as simple points. This allows certain improper integrals to be approximated. For example, LN(x) can be integrated over the interval (0,1] even though LN(0) is undefined.
A refinement was implemented in the basic Romberg scheme. If uniformly spaced sample points are taken, periodic integrands may sometimes cause a problem due to resonance phenomena. In this program the sampling has been made nonuniform by a nonlinear substitution. Such as substitution was implemented in the HP34C calculator integration key routine.
A discussion of the theory of Romberg Integration is given in [3] and [4]. Reference [5] provided the starting point for “IG”, and users are urged to consult that work.
Example

Approximate
$$
\int_0^1 \frac{4}{x^2+1}\,\rm{dx} \nonumber
$$  SIZE 030 minimum.
 Select a display setting of SCI 4.
 Key the integrand of the above integral as a function in program memory starting with a global label and ending with RTN. Assume x is in the X register and leave ƒ(x)in the X register.

01 LBL "FX1" 02 X^2 03 1 04 + 05 4 06 X<>Y 07 / 08 RTN
 Enter the limits of integration into the stack as 0 ENTER^ 1. XEQ “IG” and enter the global label name “FX1”. The following sequence of numbers will be displayed: 3.2000+00, 3.1409+00, 3.1413+00, 3.1416+00, 3.1416+00.
 This example takes about 49 seconds to run. The true answer is π, and the displayed result is accurate to the five significant digits. Switching to FIX 9, we see the last value returned is 3.141592651, but because we were in SCI 4 mode we should not expect more than 4 decimal places of accuracy.
For more examples refer to [1].
Listing
Available through the repository

Requires
 XFunctions module on the HP41cv

Available as
 source code
 raw code for the V41 emulator
 bar code for the HP Wand (HP82153A)

Size
 145 bytes including END, 2 tracks.

Note
 For line by line comments refer to [1].
01 LBL "IG" ; init by 02 "FNAME ?" ; store constants (ba)/4 in R16 03 AON ; store constants (b+a)/4 in R17 04 STOP 05 ASTO 10 06 AOFF 07 SF 10 08 STO 17 09 X<>Y 10  11 4 12 / 13 STO 16 14 ST 17 15 ST 17 16 . ; initialize Sk, k and M(k,k) for k=0 17 STO 15 18 STO 11 19 STO 18 20 SF 09 ; for at least two iterations (see 089090) 21 LBL 01 ; calculate u0 and the step size 2^(k1) 22 E 23 2 24 STO 14 25 RCL 11 26 CHS 27 Y^X 28 ST* 14 29 E 30  31 LBL 02 32 STO 12 33 X^2 ; calculate xi 34  35 STO 13 36 2 37 + 38 RCL 12 39 * 40 RCL 16 41 * 42 RCL 17 43 + 44 XEQ IND 10 ; calculate f(xi) 45 RCL 13 46 * 47 ST+ 15 ; calculate Sk 48 E 49 RCL 12 50 RCL 14 51 + 52 XZ 74 ENTER^ 75 X<> IND 12 76 ST Y 77 RND 78 X<> Z 79 / 80 RCL IND 12 81 + 82 ISG 12 83 STOP 84 DSE 13 85 GTO 03 86 STO IND 12 87 FS? 10 88 VIEW X 89 FS?C 09 90 GTO 01 91 RND ; end when two consecutive rounded approximations 92 X#Y? ; are equal 93 GTO 01 94 LASTX ; recall the final approximation and halt. 95 RTN 96 END
References
[1]  IG – Integrate John Kennedy, December 1981 PPC User Manual, pg. 220226 

[2]  Vereinfachte Numerische Integration W. Romberg, 1955 Norske videnskab. Selskab., Forh. (trondheim 28, No. 7, 1955) 

[3]  New aspects in numerical quadrature F.L. Bauer, H. Rutishauser and E.L. Stiefel, 1963 Experimental Arithmentic, High Speed Computing and Mathematics, pp. 199218 American Mathematical Society, Providence, Rhode Island 

[4]  An introduction to numerical mathematics E.L. Stiefel, 1963 Academic Press, New York 

[5]  Handheld Calculator Evaluates Integrals William H. Kahan, August 1980 HP Journal V31N8, pg. 2332 