+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (11C) Gaussian integration (/thread-381.html)
(11C) Gaussian integration - Thomas Klemm - 01-09-2014 09:38 PM
Published in PPC Journal V7N6 page 10 (ROM Progress section), Jul/Aug 1980 by Valentin Albillo
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) = 6x^5\) 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)=\frac{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.
This article is a copy of the original thread: Numerical integration on the 11C
RE: Gaussian integration for the HP-11C - Jeff_Kearns - 01-11-2014 09:23 PM
This excellent program works just as well on the HP-34C if you:
Change GSB E to GSB 1
Change LBL E to LBL 1
And then enter your f(x) to be integrated into program memory, starting at 47 LBL 1.
The execution times are most impressive!
Example 1 returns 8,303,718,967 in 9 seconds (N=1).
Example 2 returns 1.605418622 in 9 seconds (N=1)
Too bad HP didn't use this method on the 34C instead of Romberg.
Would love to see it adapted to the HP-29C... ;-)
Jeff Kearns
RE: Gaussian integration for the HP-11C - Namir - 01-14-2014 01:40 PM
Hi Thomas,
I am curious about the algorithm for your HP-11C. Is it a Gaussian Quadrature??? If so, I don't see where the quadrature weights are stored! Do you have a link to the algorithm?
Cheers,
Namir
RE: Gaussian integration for the HP-11C - Thomas Klemm - 01-14-2014 02:27 PM
(01-14-2014 01:40 PM)Namir Wrote: Is it a Gaussian Quadrature???Sure.
Quote: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.Thus I assume that 3 points are used: 0, \(\pm\sqrt{\frac{3}{5}}\). The weights are \(\frac{8}{9}\) and \(\frac{5}{9}\).
Quote:If so, I don't see where the quadrature weights are stored!
10 .6
12 SQRT
13 *
22 5
23 *
34 8
35 *
Quote:Do you have a link to the algorithm?Gaussian quadrature
HTH
Thomas
RE: Gaussian integration for the HP-11C - Namir - 01-14-2014 09:52 PM
Thank you!!! The implementation combines simplicity and effectiveness. Very clever!!
Namir
RE: Gaussian integration for the HP-11C - Jeff_Kearns - 06-18-2014 12:12 AM
PPC Journal V7N6 reference
What a shame this didn't make it onto the PPC ROM! The routine is one of the best examples of Valentin's HP RPN programming mastery.
Jeff
RE: Gaussian integration for the HP-11C - Dieter - 06-19-2014 03:13 PM
(01-11-2014 09:23 PM)Jeff_Kearns Wrote: Too bad HP didn't use this method on the 34C instead of Romberg.
Definitely not. This method surely has its advantages, but - unlike the Romberg method - it is not adaptive. The modified Romberg algorithm in the 34C will reliably handle cases that will not work well with the Gaussian method. The WP34s originally used a (much more sophisticated) Gauss quadrature algorithm. This was changed to Romberg's method because of the mentioned reasons.
Quote:Would love to see it adapted to the HP-29C... ;-)
Well, that's easy. Simply change the labels (e.g. 0→2, A→0 and E→1) and replace register 0 with 4 and I with 0 – the 19C and 29C use R0 instead of I as their loop counter. Of course "DSE" becomes "DSZ" here.
Dieter
RE: Gaussian integration for the HP-11C - Paul Dale - 06-20-2014 01:46 AM
(06-19-2014 03:13 PM)Dieter Wrote: The modified Romberg algorithm in the 34C will reliably handle cases that will not work well with the Gaussian method. The WP34s originally used a (much more sophisticated) Gauss quadrature algorithm. This was changed to Romberg's method because of the mentioned reasons.
I would have liked to do an adaptive algorithm that used a Gauss-Kronrod quadrature for the sub steps. The GSL integration uses a variety of Gauss-Kronrod quadratures of increasing point count (but reusing the same points).
- Pauli
RE: Gaussian integration for the HP-11C - walter b - 06-20-2014 05:27 AM
It's a pity both Valentin and Karl are no longer posting on this forum.
d:-(