Chapter 2

                       NUMERICAL INTEGRATION

                              Summary

               Numerical integration functions support the
               evaluation of integrals and the numerical
               solution of systems of differential equations.
               Areas covered by these functions are:

                     o Numerical Integrals
                     o Differential Equations


 Notes on Contents

     The numerical integration functions are used to evaluate integrals and to
 solve systems of ordinary differential equations.

 o  Numerical Integrals:

    fintg  -------- integrate a function over a definite interval.
    chintg  ------- compute a Tchebycheff expansion to evaluate an
                    integral as a function of its upper limit.
    fchb  --------- evaluate the Tchebycheff expansion generated by
                    chintg.

 o  Differential Equations:

    deqsy  -------- solve a system of first order differential
                    equations.


 General Technical Comments

     The function 'fintg' uses an extrapolation technique to adjust step size.
 The efficiency of this step size adjustment can be enhanced by dividing the
 range of integration into a number of subregions and employing a separate
 function call for each subregion. This technique is particularly important
 when the integrand is singular at some points. The function chintg develops
 a Tchebycheff series to evaluate the integral as a function of its upper
 limit. This capability is useful when the function represented by the
 integral must be evaluated at many points.

     The integration of a system of first order differential equations uses
 the Bulirsch-Stoer technique combining a modified midpoint integration step
 with repeated Richardson extrapolation. This yields highly accurate and
 reliable solutions. Differential equations of higher order can readily be
 represented as systems of first order equations.

     Functions in this library receive pointers to user defined functions,
 with a specified set of call parameters. Use external variables to deliver
 additional parameters to these integrand functions.


-------------------------------------------------------------------------------
                         FUNCTION SYNOPSES
-------------------------------------------------------------------------------

 Numerical Integrals:
-------------------------------------------------------------------------------

fintg

     Compute the integral of func(x) over a definite interval.

     #include <math.h>
     double fintg(double a,double b,int n,double te,double (*func)())
       a = lower limit of integration
       b = upper limit of integration
       n = initial step number (step size = (b-a)/n)
       te = convergence criteria threshold
            ( te specifies the maximum relative correction
              to the integral. )
       func = pointer to integrand evaluation function
               ( called by: y = (*func)(x). )
      return value: I = value of integral
                    ( I = HUGE > convergence failure )

            I = Intg(a to b){ func(x) dx } .

     ------------------------------------------------------------------

chintg

     Compute a Tchebycheff expansion supporting evaluation of an integral
     as a function of its upper limit.

     #include <math.h>
     double chintg(double a[],int m,double (*func)())
       a = m+1 dimensional array of expansion coefficients
       m = maximum degree of expansion
       func = pointer to integrand evaluation function
               ( called by: y = (*func)(x). )
      return value: h = maximum absolute value among last 3 coefficients

         The integral evaluation is given by:

          I(x) = Intg(y= -1 to x){ func(y)dy }=
              Sum(i=0 to m){ a[i]*Ti(x) } .

     ----------------------------------------------------------------------

fchb

     Evaluate a Tchebycheff series, such as the chintg function's output.

     double fchb(double x,double a[],int m)
       x = value of independent variable ( -1.0 <= x <= 1.0 )
       a = m+1 dimensional array of series coefficients
       m = maximum degree of polynomial in the series
      return value: S = value of the series

          S = Sum(i=0 to m){ a[i]*Ti(x)} , where Ti(x) is
          the ith Tchebycheff polynomial.

-------------------------------------------------------------------------------

 Differential Equations:
-------------------------------------------------------------------------------

deqsy

     Solve a system of first order differential equations,
     dy[k]/dx = f[k](x,y), using an extrapolation technique.

     int deqsy(double y[],int n,double a,double b,int nd,double te, \
               int (*fsys)())
       y = n-vector of system variables (dy[k]/dx = f[k](x,y))
       n = dimension of system
       a = initial value of x
       b = final value of x (solution target point)
       nd = initial step number (step size = (b-a)/nd)
       te = convergence criteria threshold
             ( te is the maximum relative correction to the components
               of y. ( |dr[k]| <= te*|y[k]| for all k )
       fsys = pointer to function which evaluates the derivatives
               of system variables
               ( (*fsys)(double x,double *y,double *dr) returns the
                 values of dy[k]/dx in the n dimensional array dr )
      return value: m = number of extrapolations
                        m<0 -> convergence failure