(11C) Gaussian integration

01092014, 09:38 PM
(This post was last modified: 06152017 01:17 PM by Gene.)
Post: #1




(11C) Gaussian integration
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 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) 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 Nsubinterval Gaussian integration, just do the following:
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: 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. This article is a copy of the original thread: Numerical integration on the 11C 

01112014, 09:23 PM
(This post was last modified: 01162014 09:27 PM by Jeff_Kearns.)
Post: #2




RE: Gaussian integration for the HP11C
This excellent program works just as well on the HP34C 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 HP29C... ;) Jeff Kearns 

01142014, 01:40 PM
Post: #3




RE: Gaussian integration for the HP11C
Hi Thomas,
I am curious about the algorithm for your HP11C. 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 

01142014, 02:27 PM
Post: #4




RE: Gaussian integration for the HP11C
(01142014 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 

01142014, 09:52 PM
(This post was last modified: 01142014 11:43 PM by Namir.)
Post: #5




RE: Gaussian integration for the HP11C
Thank you!!! The implementation combines simplicity and effectiveness. Very clever!!
Namir 

06182014, 12:12 AM
Post: #6




RE: Gaussian integration for the HP11C
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 

06192014, 03:13 PM
Post: #7




RE: Gaussian integration for the HP11C
(01112014 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 HP29C... ;) 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 

06202014, 01:46 AM
Post: #8




RE: Gaussian integration for the HP11C
(06192014 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 GaussKronrod quadrature for the sub steps. The GSL integration uses a variety of GaussKronrod quadratures of increasing point count (but reusing the same points).  Pauli 

06202014, 05:27 AM
Post: #9




RE: Gaussian integration for the HP11C
It's a pity both Valentin and Karl are no longer posting on this forum.
d:( 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)