(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