(11C) Gaussian integration

[Thomas Klemm]

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:

1. 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.
2. 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
3. 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

[Jeff_Kearns]

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

[Namir]

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

[Thomas Klemm]

(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

[Namir]

Thank you!!! The implementation combines simplicity and effectiveness. Very clever!!

Namir

[Jeff_Kearns]

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

[Dieter]

(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

[Paul Dale]

(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

[walter b]

It's a pity both Valentin and Karl are no longer posting on this forum.

d:-(