Integrate on HP-41

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 step-size 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

1. SIZE 030 is the recommended minimum. A few integrals may require a larger size.
2. 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.
3. If a printer is connected, the approximations will be printed.
4. 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 R10-R29 or use flags F09 and F10.
5. Alpha-store the global label name from step 4 (six or less alpha characters) in R10.
6. 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.
7. 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 X-register 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 Euler-MacLaurin 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 shut-off 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 sub-intervals or divisions which decreases the error, without actually increasing the number of sub-intervals. This process is called extrapolation to the limit.

Romberg integration is iterative, automatic, and non-adaptive. 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 non-adaptive 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 non-uniform by a non-linear substitution. Such as substitution was implemented in the HP-34C 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

1. Approximate

   $$
   \int_0^1 \frac{4}{x^2+1}\,\rm{dx} \nonumber
   $$

2. SIZE 030 minimum.
3. Select a display setting of SCI 4.
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.

5. 01 LBL "FX1"
   02 X^2
   03 1
   04 +
   05 4
   06 X<>Y
   07 /
   08 RTN

6. 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.
7. 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
  • X-Functions module on the HP-41cv
• 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      "F-NAME ?"       ;   store constants (b-a)/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 089-090)
21      LBL 01          ; calculate u0 and the step size 2^(k-1)
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      X<Y?
53      GTO 02
54      RCL 11          ; calculate M(k,0)
55      STO 13
56      18
57      STO 12
58       E
59      ST+ 11
60      RCL 15
61      RCL 16
62      1.5
63      *
64      *
65      RCL 14
66      *
67      LBL 03          ; calculate M(k,j)
68      R^
69      4
70      *
71      ENTER^
72      DSE Y
73      X<> Z
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