|Integration speeds: JK's HP 34C review |
Message #9 Posted by Karl Schneider on 4 May 2010, 1:50 a.m.,
in response to message #2 by Gene Wright
(Gene Wright) did some additional comparisons of this integration program with the examples given in (an) HP-34C review from the PPC Journal V6N5P19-20.
Thanks for presenting the table. The author's review of the HP-34C was favorable, yet he chose six integrals that, in a sense, made the HP-34C's INTEG capability 'look bad' by comparison to 16-point Gaussian quadrature. The GQ results, the author points out, were at least as accurate and were obtained in less time. Furthermore, they were produced by an RPN routine programmed on the HP-34C -- not a built-in microcoded function. He speculates that the GQ approach "may not have been used by HP because of internal storage requirements."
Maybe so, but I notice that all the examples utilized very well-behaved integrands over narrow intervals, for which even simple methods will tend to produce good results. HP's objective in this inaugural effort was robustness and practicality -- to provide sufficiently-accurate real-world results in science and engineering without the user's having to specify the number of evaluations, or avoid endpoints for which the function was undefined, or risk wrong answers due to an inadequate sampling rate for a periodic function.
I do note that the GQ program uses all 21 storage registers and consumes 40 lines of programming. The user's program defining the integrand must not employ any register, because no unused ones are available. It also cannot exceed 30 lines, lest register R19 ('.9') be automatically deallocated, rendering the GQ program inoperative. (Score a point for the successor HP-15C and its manual memory allocation.)
The exact value for the fifth integral is 6 - 236*e^(-5) = 4.40984450821.
The exact value for the sixth integral is 3*pi/96 = 0.0981747704247.
The HP-34C GQ result and timing for the "Kahan integral" in the paper linked by Gerson:
f(x) = sqrt(x)/(x-1) - 1/ln(x)
0 <= x <= 1
At each endpoint, f(x) is undefined but does approach zero, as shown by limit analysis. From x = 0.000, f(x) rises quickly to a 'smooth' local maximum of 0.117680329841 at x = 0.00473796324935, then gradually and monotonically declines to a local minimum of ...
-0.117680329841 at x = 211.061155894 (which is 1 / 0.00473796324935), as this function has the property f(1/x) = -f(x) for all x (where defined).
1 (x + 1)
df/dx = ---------- - -----------------
x*(ln x)^2 2*sqrt(x)*(x-1)^2
The exact value of the integral is (2 - ln 4 - gamma), where gamma is the Euler-Mascheroni Constant of 0.577215664901532...
FIX 5 FIX 6 'Exact' Gaussian
0.03661813695 0.03649023200 0.036489973980 0.0364900935
19 seconds 268 seconds 45 seconds
The 16-point GQ outperforms the built-in Romberg method in this instance, too. Over this short interval of [0, 1], the 16 points seem to provide adequate coverage for high accuracy, and it does not evaluate at endpoints -- a definite 'plus' in this instance.
The highly-accurate GQ result might have been 'dumb luck', as no evaluation point (node) fell between zero and 0.004738.
I do wonder about its performance over large intervals with erratic integrands, but the Gauss - Kronrod method (used by the TI-82 and likely its successors) handles that situation by comparing the difference between estimations and taking more evaluations as necessary, much as HP's modified Romberg method does.
I may add a section to my HP Articles Archive contribution (#556), with a deeper analysis of this integral and computational methods.
Gerson said above:
... but even Simpson 3/8 with 19500 intervals performs badly on this one:
Anyway, this is a tough integral for numerical methods:
In reasonable amounts of time with endpoints of 1E-8 and 1 - 1E-8, the Casio fx-115MS obtains 0.0365 using Simpson's Rule at 512 points, but the successor fx-115ES obtains 0.03648997346 using Gauss-Kronrod (and, I might add, a perfect answer of 0.03648997398 using 1E-10 and 1 - 1E-10 as the endpoints). Then, after a period of inactivity, they shut themselves off and forgot the integrals I'd entered...
Gene said below:
Karl, what I'd rather do is make sure you get an HP 30b soon to try your hand at programming it for things like this.
Nah, I'd rather not program things that should be built into a calculator, even though the HP-30b is a business model. Besides, the HP-34C LED 'light show' should be played from time to time...
Edited: 11 May 2010, 11:48 a.m. after one or more responses were posted