The Museum of HP Calculators

HP Articles Forum

Romberg Integration on the wp34s

Posted by Les Wright on 4 Dec 2011, 7:43 p.m.

Romberg's Method of numerical integration entails iteratively applying a simple quadrature scheme such as the trapezoidal rule or midpoint (rectangle) rule by adding more function evaluations at each pass through the algorithm and accelerating convergence using the century-old method of Richardson extrapolation. The underlying mathematics are not terribly complicated, and I can recommend the Wikipedia article on the subject. (However, one should keep in mind that that article's treatment builds the algorithm around the trapezoidal rule, whereas what follows uses the rectangle rule. But the principles are the same.)

An adapted version of Romberg's Method underlies the numerical integration abilities of the HP34C, HP15C, HP33s, HP35s, and the INTEG routine of the HP41C Advantage Pac. As of version 2.2, it is NOT used by the wp34s, which instead uses a very fast and good 10/21 Gauss-Kronrod scheme. However, some quadrature aficionados may prefer a small and relatively fast alternative on the wp34s that performs more like the built-in integrator on the above calculators. I can take no credit for the underlying programming, but I am glad to provide the following port.

The following is a nearly keystroke-by-keystroke translation of the HP41C PPC ROM Romberg integrator, IG. I have long been fascinated by the economy of this routine since discovering it in the PPC ROM manual on the MOHPC DVD set. The earliest HP implementation of Romberg's Method was in the 34C by William Kahan. In keeping with Kahan's work, IG just doesn't apply Romberg directly to the function. There are first two transformations of variable--the first a linear one to standardize the interval of integration to [-1,1], the second a nonlinear one so that points are sampled at uneven intervals and problems with some periodic functions are avoided. It is at this point that I put in an unapologetic plug for the MOHPC DVD--the PPC ROM manual and the August 1980 HP Journal article describing the 34C integration capacities are, for me, in themselves well worth the cost.

Without ado, here is the routine. It is in proper syntax and can be simply cut, pasted, and assembled for use in the emulator and transferred to the calculator:

```001 LBL'IG'
002 STO 17
003 x[<->] Y
004 -
005 4
006 /
007 STO 16
008 STO- 17
009 STO- 17
010 0
011 STO 15
012 STO 11
013 STO 18
014 SF 09
015 LBL 01
016 1
017 ENTER[^]
018 2
019 STO 14
020 RCL 11
021 +/-
022 y[^x]
023 STO[times] 14
024 1
025 -
026 LBL 02
027 STO 12
028 x[^2]
029 -
030 STO 13
031 2
032 +
033 RCL[times] 12
034 RCL[times] 16
035 RCL+ 17
036 XEQ[alpha]
037 RCL[times] 13
038 STO+ 15
039 1
040 RCL 12
041 RCL+ 14
042 x<? Y
043 GTO 02
044 RCL 11
045 STO 13
046 1
047 8
048 STO 12
049 1
050 STO+ 11
051 RCL 15
052 RCL 16
053 1
054 .
055 5
056 [times]
057 [times]
058 RCL[times] 14
059 LBL 03
060 R[^]
061 4
062 [times]
063 ENTER[^]
064 DSE Y
065 x[<->] Z
066 ENTER[^]
067 x[<->][->]12
068 STO- Y
069 ROUND
070 x[<->] Z
071 /
072 RCL+[->]12
073 INC 12
074 DSE 13
075 GTO 03
076 STO[->]12
077 FS? 10
078 VIEW X
079 FS?C 09
080 GTO 01
081 ROUND
082 x[!=]? Y
083 GTO 01
084 RCL L
085 RTN
```

The program uses the 4-stack and R11 through R[->]12. As we can see in the Wikipedia article (and the PPC ROM article for those who have it), Romberg's Method produces a lower triangular matrix with the best estimate after a given iteration located in the rightmost cell of the last row. In IG, only the most recent row is saved in R18 to R[->]12. The routine continues until the rightmost elements of the two most recent rows agree to the desired number of digits. (Note that in the Wikipedia article a different convergence criterion is used--i.e., the last two elements of the last row agree. This is less generally stringent.)

Since comparison of rounded values determines when IG stops, the setting of FIX, SCI, ENG or ALL determines accuracy, as it does in old HP calcs (except for the 42S, where an actual number for accuracy was entered). In general, SCI nn gives more predictable results--we know, for example, that SCI 04 will give us five significant digits, whereas with FIX nn that can vary according to the number of leading zeros. ALL nn gives the same results as SCI 11 for any nn.

As in the original IG, if flag 10 is set interim estimates are displayed until convergence is achieved. I like keeping this set so I can watch the convergence happen. I can tell things are done when the little "RCL" annunciator stops flashing. But this port differs from the original IG in that the global alpha label of the routine to be integrated is stored in the alpha register, not R10. IG as it is presently written does not integrate routines labeled by numeric or local alpha labels (unless a globally labeled wrapper is provided). This behaviour makes it similar to the Advantage Pac's INTEG.

Here are instructions for use:

1. Enter the function to be integrated. It must be labeled by a global alpha label, end with RTN, not use R11 through to about R30 (rarely IG will need more for the last row of the Romberg matrix), and return the function value to the X register.

2. Enter the alpha name of this function in the alpha register.

3. Set flag 10 if you want to see interim results.

4. Set the desired accuracy through the display settings--SCI nn is recommended to give nn+1 significant digits.

5. Put the limits of integration on the stack: lower [ENTER] upper.

6. Execute the routine ([XEQ] [ENTER] IG [ENTER]) and watch the magic unfold.

Here is an example given in the 1980 Kahan article: integrate 2*w^2/((w-1)*(w+1)) - w/ln(w) from w = 0 to w = 1.

I enter the following routine:

```001 LBL'KAH'
002 FILL
003 1
004 -
005 x[<->] Y
006 1
007 +
008 [times]
009 /
010 ENTER[^]
011 +
012 x[<->] Y
013 LN
014 1/x
015 -
016 [times]
017 RTN
```

`SF 10 0 [ENTER] 1 [f] [alpha] KAH [ENTER] [SCI] 11 [XEQ] [ENTER] IG [ENTER]`
eventually gives 3.64899739786E-2. This is indeed the correct result to 12 digits. It took about 80 seconds on my wp34s, but I should note my calc is clocked down to half speed due to battery drain. When the batteries were fresher, this integral took at most 45 seconds. Of course, it is much faster if the desired digits are set lower.

Of note, SHOW tells us that the internal 16-digit result is 3.648997397857771e-2. According to Mathematica, the integral does indeed have an exact result--it is 2 - Euler's gamma - ln(4)--which is 3.648997397857652e-2 to 16 digits. This is remarkably close agreement.

I close with a couple of notes. First of all, there is nothing particularly sacred about using R11 and up, and that can be changed. Second, some users may choose to modify the code so that no global alpha labelling of the function routine is necessary. Finally, in contrast with to the output of the 34C, 35s, and 15C, we get no error estimate explicitly with this routine. Personally, I have been using R[->]12 - R18, the difference between the final result and the most recent unextrapolated rectangle estimate (stored as the first element in the the most recent row for the Romberg tableau) as a rough upper limit on the absolute error.

Reference: