FYI: Update to my SOLVE/INTEG article (#556)

[Karl Schneider]

All --

For those who are interested, I've made some corrections to my article (#556) in the MoHPC Articles Forum. The changes pertain to INTEG input-function uncertainty and estimated error for the pre-Saturn models. This particular material about INTEG (for all RPN-based models), as well as a list of primary references, was added in early August.

I'll be able to answer any questions only through Saturday; then it'll have to wait until the week after next.

http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/articles.cgi?read=556

-- KS

Edited: 15 Sept 2006, 12:03 a.m.

[Kiyoshi Akima]

I know this question gets asked periodically, but I have yet to see the definitive answer so I'll risk the wrath of the long-timers and ask it again.

Does anybody know the precise algorithm used by the internal integrator on the 34C and 15C (and/or later models)?

I know it's a Romberg with a non-uniform transformation. That much is documented in Kahan's 1980 article and expanded upon in the PPC ROM manual. However, the PPC ROM routine (which was intended to be similar but not identical to the HP rountine) does not give the same results as the 34C/15C or the 41C AdvantagePac (beyond what can be accounted for by the difference in precision). Hugh posted a routine which apparently gives the same results as the 34C/15C with seven function evaluations, but differs with fifteen or more.

Also, the 34C/15C can produce different values for the same number of function evaluations. I don't mean the 34C and the 15C differ, I mean on the same machine you can get two different values for the same set of function evaluations. Try integrating exp(-x^2) from 0 to 1 in FIX 3 and FIX 4.

Any help would be appreciated.

[Karl Schneider]

Hi, Koyoshi --

You posted:

Quote:

---

Also, the 34C/15C can produce different values for the same number of function evaluations. I don't mean the 34C and the 15C differ, I mean on the same machine you can get two different values for the same set of function evaluations. Try integrating exp(-x^2) from 0 to 1 in FIX 3 and FIX 4.

---

I got the following results:

FIX     Answer     Estimated error  Elapsed Time     

 3    0.746825137      0.0005          23.32 sec

 4    0.746823672      0.00005         24.96 sec

 9    0.746824133      0.0000000005    88    sec

How do you know that the sets of samples were identical? There was at least a half-second difference in time of execution, although it's difficult to ascribe that to more samples taken.

The estimated errors are as indicated by the formula given in the HP-34C and HP-15C manuals, which I quoted in my article. Only five digits (including the leading zero) are needed to display the value of the integrand under FIX 4, so no "roundoff" was performed.

-- KS

[Kiyoshi Akima]

Quote:

---

How do you know that the sets of samples were identical? There was at least a half-second difference in time of execution, although it's difficult to ascribe that to more samples taken.

---

I put code like [1] [STO + 0] [Rdown] inside the function being integrated. This gives me the number of function evaluations (assuming you know what was in R0 at the beginning). The set of points at which the function is evaluated is known from the non-uniform transformation. For a given number of evals, the set of points is always the same (before the transformation).

I've noticed the slight time difference, though I didn't measure it to 0.01 seconds :-) What it's doing in that time, I don't know.

The best I can figure is that in FIX 3 it was satisfied with the initial result and presented that. In FIX 4 it wasn't satisfied, but rather than doubling the number of points it first tried some other tweak, which did satisfy it. In FIX 5 the tweak wasn't enough so it does go ahead and double the number of points.

But just what that "tweak" on top of the Romberg extrapolation is, I don't know. And if that tweak only takes half a second, why doesn't it always apply it?

BTW, have you ever looked at all ten digits of a typical error estimate? Where is that other non-zero digit coming from, and why isn't it always in the same place relative to the first one?