Re: Long Live the 33C!!! Message #5 Posted by Les Wright on 24 Dec 2006, 3:27 a.m., in response to message #1 by Les Wright
Without further ado, here is a 2point Gaussian quadrature routine for the 33C, which is basically Valentin Abillo's excellent 3point routine with a few steps saved since one of the function computations per subinterval is removed to make room for more steps from step 39 on for the function to be integrated:
01 STO 1
02 
03 RCL 7
04 /
05 STO 0
06 2
07 /
08 ST+ 1
09 3
10 SQRT
11 1/X
12 *
13 STO 3
14 CLX
15 STO 2
16 RCL 1
17 RCL 3
18 +
19 GSB 39
20 ST 2
21 RCL 1
22 RCL 3
23 
24 GSB 39
25 ST 2
26 RCL 0
27 ST+ 1
28 1
29 ST 7
30 RCL 7
31 X!=0?
32 GTO 16
33 RCL 0
34 RCL 2
35 *
36 2
37 /
38 RTN
To use, put the function to be integrated from line 39 on, making sure the result is returned to the X register, and end that segment with RTN. Then store the number of desired subintervals in R7, and the limits of integration on the stack as lower ENTER upper. Then GSB 01 and voila!
This is not nearly as effective as Valentin's routine. On a given interval of say width h, both threepoint Gaussian quadrature and the 5point NewtonCotes formula (Bode's Rule) have roughly the same theoretical error, though of course the Gaussian version saves two function evaluations on the first subinterval and one on every subinterval thereafter (since with Bode's Rule the connecting endpoints need be computed only once). So if you divide up the area of integration into 4 subintervals, Bode's Rule requires 5 + 3*4 = 17 function evaluations, whereas threepoint Gaussian quadrature requires 12. It is easy to see that if we divide up the area of integration into 6 subintervals for a 3point Gaussian integration, the accuracy should be much better than Bode's Rule at the cost of only one extra evaluation of the function. I would guess that it is also superior Simpson's rule computed over nine subintervals (19 function evalutions), but since the subinterval size is changing this is a little harder for me to demonstrate to myself.
With 2point Gaussian quadrature, the theoretical error on an interval of width h is of the order h^5/4320 (times some value of the fourth derivative of the function being integrated), vs. 3point NewtonCotes quadrature (Simpson's Rule), which is of order h^5/2880. So, at least in theory, perceptibly better. However, the savings in function evaluation isn't as impressiveSimpson's rule computes three times on the first subinterval, and twice on every subinterval thereafter, whereas repeated application of 2point Gaussian quadrature is simply twice per subinterval. So there is a saving of only one function evalution. However, I would guess, on balance, this routine is preferable to Simpson's rule as you get more accuracy with one fewer function evaluations, and one need not compute the function at the endpoints.
As a quick comparison, I computed the integral of sin(x)/x from 0 to 2 using the HP41 Math Pac's INTG routine (Simpson's Rule) with n = 12 (13 function evaluations), Valentin's original routine with 4 sub intervals (12 function evaluations), and the above routine with 6 subintervals (12 function evaluations). This is what I get:
Simpson's: 1.605413995
Abillo's: 1.605412978
Les's: 1.605412298
The actual 10 digit answer is 1.605412977, so Valentin's routine gives a full 10 digits within 1 ULP. Both mine and Simpson's get it to 7 digits within 1 ULP after rounding, but mine is a wee bit closer when you take it out to 8 digits7 ULP vs. 10 ULP. This is not earthshattering, but it is gratifying since it got a little closer with one less evaluation of the function.
I doubt that this routine has any practical implications given the superior integrators in HP calculators starting from the 34C and right up to the otherwise much maligned 33S and of course delightful 42S. But I do like the challenge of the small program, and I positively love the 33C, despite its modest powers. It was a great pleasure for me to adapt Valentin's excellent work to a more modest end, and I hope that the pleasure can be shared.
Les
Edited: 24 Dec 2006, 4:03 a.m.
