Re: Numerical integration on the 11C Message #4 Posted by Valentin Albillo on 17 June 2006, 11:57 p.m., in response to message #1 by Gerson W. Barbosa
Hi, Gerson:
Nice program, but Simpson's method is to numerical integration what bubble sort is to sorting, i.e., vastly inefficient. The best methods, and the ones who would've been included as the basis of INTEG in the HP34C (and subsequently HP15C, HP71B's Math ROM, etc) are Gaussian ones, were it not for the fact that they do require a certain number of fullprecision real constants that would require a noticeable amount of ROM space, not available in the HP34C project, among other considerations, as discussed in Mr. Kahan's wonderful article on INTEG, originally published in the HP Journal eons ago.
However, Gaussian methods are still way better than Simpson's rule even for lowpoint variants. For instance, here's my HP11C version of a Gaussian integration routine I originally wrote for the HP41C, as published in PPC Journal V7N6 page 10 (ROM Progress section), Jul/Aug 1980, i.e. exactly 26 years ago:
Gaussian integration for the HP11C
01 LBL A 17 LBL 0 32 RCL 1
02 STO 1 18 RCL 1 33 GSB E
03  19 RCL 3 34 8
04 RCL I 20 + 35 *
05 / 21 GSB E 36 STO2
06 STO 0 22 5 37 RCL 0
07 2 23 * 38 STO+1
08 / 24 STO2 39 DSE
09 STO+1 25 RCL 1 40 GTO 0
10 .6 26 RCL 3 41 RCL 2
12 SQRT 27  42 *
13 * 28 GSB E 43 18
14 STO 3 29 5 45 /
15 CLX 30 * 46 RTN
16 STO 2 31 STO2 47 LBL E
As you may see, it's a very small, 46step program which uses just R0R3 for scratch and RI as a decrementing counter for the number of subintervals. It delivers exact results, even using just 1 subinterval, for f(x) being a polynomial of degrees up to (and including) 5^{th}, while evaluating f(x) just 3 times per subinterval. That's about twice as precise as Simpson's rule, which delivers exact results for polynomials up to 3^{rd} degree only (not 2^{nd}, as stated in another post).
To use it to compute the integral of an arbitrary f(x) between x=a and x=b, using Nsubinterval Gaussian integration, just do the following:
 Enter your f(x) to be integrated into program memory, starting at
47 LBL E
ending it either with a RTN instruction or with the end of program memory.
 Store N, the number of subintervals you want to use, in Register I. N must be an integer number equal or greater than 1. The larger N, the more precise the result will be and the longer it will take to run.
N, STO I
 Enter the limits of integration, a and b, into the stack and call the integration routine:
a, ENTER, b, GSB A
The computation will proceed and the result will be displayed upon termination. Let's see a couple of examples:
1. Compute the integral of f(x) = 6x^{5} between x=6 and x=45.
As f(x) is a 5^{th} degree polynomial, we expect an exact result, save for minor rounding errors in the very last place. Let's compute it:
 Define f(x):
LBL E, 5, y^x, 6, *, RTN
 Store the number of subintervals, just one will do:
1, STO I
 Enter the limits and compute the integral:
6, CHS, ENTER, 45, GSB A
 The result is returned within 6 seconds:
> 8,303,718,967
The exact integral is 8,303,718,969, so our
result is exact to 10 digits within 2 ulps, despite using
just one subinterval and despite the large interval of integration.
Now for your own example:
2. Compute the integral of f(x)=Sin(x)/x between x=0 and x=2
 Set RAD mode and define f(x):
LBL E, SIN, LASTX, /, RTN
 Store the number of subintervals, let's try just 1:
1, STO I
 Enter the limits and compute the integral:
1E99, ENTER, 2, GSB A
 The result is returned within 5 seconds:
> 1.605418622
Testing with 2,3, and 4 subintervals we get (the exact integral being 1.605412977):
N subintervals Computed integral Time

1 1.605418622 5 sec.
2 1.605413059 11 sec.
3 1.605412984 16 sec.
4 1.605412978 22 sec.
so even using just 2 subintervals does provide 8digit accuracy, and using 4 nails down the result to 10 digits save for a single unit in the last place.
That's all. I think you'll agree this simple Gaussian method beats Simpson's and Trapezoidal Rules hands down. :)
Best regards from V.
