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 HP-34C (and subsequently HP-15C, HP-71B's Math ROM, etc) are Gaussian ones, were it not for the fact that they do require a certain number of full-precision real constants that would require a noticeable amount of ROM space, not available in the HP-34C 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 low-point variants. For instance, here's my HP-11C version of a Gaussian integration routine I originally wrote for the HP-41C, as published in PPC Journal V7N6 page 10 (ROM Progress section), Jul/Aug 1980, i.e. exactly 26 years ago:
Gaussian integration for the HP-11C
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 STO-2
06 STO 0 22 5 37 RCL 0
07 2 23 * 38 STO+1
08 / 24 STO-2 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 STO-2 47 LBL E
As you may see, it's a very small, 46-step program which uses just R0-R3 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) 5th, 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 3rd degree only (not 2nd, as stated in another post).
To use it to compute the integral of an arbitrary f(x) between x=a and x=b, using N-subinterval 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) = 6x5 between x=-6 and x=45.
As f(x) is a 5th 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:
1E-99, 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 8-digit 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.
|