2. Solving Integrals

Numerical integration evaluates definite integrals of the form:

where a and b are given and f(x) is either an analytical function or is given by a table of values. Geometrically, J is the area under the curve of f(x) between a and b. When f(x) is such that there is a corresponding differentiable function F(x); that is, where:

F'(x) = f(x)

then J can be evaluated using:

J = F(b) - F(a)

There are many natural phenomena and man-made systems which cannot be represented in this convenient analytical form. Then, numerical integration methods are used to approximate the integral of f(x). The simplest methods subdivide the interval between a and b into n equal widths using:

h = (b - a)/n

A reasonably accurate approximate integral is obtained using the trapezoidal rule:

J = h·[0.5·f(a) + f(x₁) + f(x₂) + . . . + f(x_n-1) + 0.5·f(b)]

Coding the trapezoidal rule algorithm

An algorithm that computes the numerical solution to an integral equation must be given the following information:

1. The initial value of the independent variable t.
2. The number of iterations or segments to achieve a desired accuracy.
3. The interval length from the initial value of t to the final value of t.

It must then return a value of the integrand at the end of a computation interval.

A pseudocode representation of the process is shown here:

    begin

        t := a
        h := (b - a)/n
        J := f(a)/2

        if n > 1 then

            i := 1

            loop

                exit when i = n

                t := t + h
                J := J + f(t)
                i := i + 1

            end loop

        end if

        J := J + f(b) / 2
    
        J := J * h

    end

Example

The C function integrate( double, double, int, double (f)(double t) ) in the example program finds the definite integral of a specified function using the trapezoidal rule for numerical integration. In the example, the function to be integrated is:

    y(t) = (1 - e-t)·cos(2·pi·t) 

/*----------------------------------------------------------------------------*\
**              Put your mathematical functions to plot here.                 **
\*----------------------------------------------------------------------------*/

int calculate( POINT *, POINT *, int );
double integrate( double, double, int, double (f)(double) );
double Y( double );

/*----------------------------------------------------------------------------*
**  "run" performs calculations and plots results.  
**
**  This is a shell in which the function calls for a specific calculation
**  are entered.
**
**  returns: nothing
*/
void run()
{
    POINT curve1[110], curve2[110];
    int i;
            
    // first, define a plot region 
    // this must be done before calculating points to be plotted
    np.plot_region( -0.1, 2.0, 2.0, -2.0 );

    i = calculate( curve1, curve2, 110 );

    // plot the results
    np.draw_axes();
    np.red_curve( curve1, i );
    np.blue_curve( curve2, i );
}

/*----------------------------------------------------------------------------*
**  "calculate" is a shell for the calculation algorithms.
**
**  returns: number of items in POINT array
*/
int calculate( POINT *a0, POINT *a1, int points )
{
    int i;
    double amps  = 0.0; 
    double t     = 0.0;
    double delta = 0.02;

    a0[0] = np.plot_point( 0, 0 );
    a1[0] = np.plot_point( 0, 0 );

    for( i = 1; i < points; i++ )
    {
        // next iteration
        amps += integrate( t, delta, 2, Y );
        t    += delta;
	
        // append results to curves
        a0[i] = np.plot_point( t, amps );
        a1[i] = np.plot_point( t, Y( t ) );
    }
    return i;
}

/*----------------------------------------------------------------------------* 
**  "integrate" finds the definite integral of a specified function
**  using trapezoidal integration.
**
**  returns: nothing
*/
double integrate( double t_beg,                 // lower limit of integral
                  double delta,                 // period for integration
                  int    n,                     // number of segments > 0 
                  double (f)(double t) )        // integrand function
{
    double area;
    double t = t_beg;                           // initialize time
    double dt = delta / n;                      // initialize step size
    int    i;                                   // loop counter

    area = f( t_beg ) / 2.0;                    // initial value of integral 

    if( n > 1 )                                 // compute the n - 1 terms
    {                                           //   in the summation
        for( i = 1; i < n; i++ )
        {
            t    += dt;                         // increment time
            area += f( t );                     // add next term to sum
        }
    }

    area += f( t_beg + delta ) / 2.0;           // add last term
    area *= dt;                                 // scale by step size

    return area;
}

const double TWO_PI = 4 * acos( 0 );

/*----------------------------------------------------------------------------*
**  "Y" is the integrand function for this problem.
**
**  returns: value corresponding to argument
*/
double Y( double t )
{
    return ( ( 1 - exp( -t ) ) * cos( TWO_PI * t ) );
}