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
Quote:
(Gene Wright) did some additional comparisons of this integration program with the examples given in (an) HP34C review from the PPC Journal V6N5P1920.
Gene 
Thanks for presenting the table. The author's review of the HP34C was favorable, yet he chose six integrals that, in a sense, made the HP34C's INTEG capability 'look bad' by comparison to 16point 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 HP34C  not a builtin 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 wellbehaved 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 sufficientlyaccurate realworld 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 HP15C 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 HP34C GQ result and timing for the "Kahan integral" in the paper linked by Gerson:
f(x) = sqrt(x)/(x1)  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).
The derivative:
1 (x + 1)
df/dx =   
x*(ln x)^2 2*sqrt(x)*(x1)^2
The exact value of the integral is (2  ln 4  gamma), where gamma is the EulerMascheroni 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 16point GQ outperforms the builtin 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 highlyaccurate 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 TI82 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:
Quote:
... but even Simpson 3/8 with 19500 intervals performs badly on this one:
3.6498295709E02
Anyway, this is a tough integral for numerical methods:
In reasonable amounts of time with endpoints of 1E8 and 1  1E8, the Casio fx115MS obtains 0.0365 using Simpson's Rule at 512 points, but the successor fx115ES obtains 0.03648997346 using GaussKronrod (and, I might add, a perfect answer of 0.03648997398 using 1E10 and 1  1E10 as the endpoints). Then, after a period of inactivity, they shut themselves off and forgot the integrals I'd entered...
Gene said below:
Quote:
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 HP30b is a business model. Besides, the HP34C LED 'light show' should be played from time to time...
 KS
Edited: 11 May 2010, 11:48 a.m. after one or more responses were posted
