Numerical Integration and Differentiation

 

Given the existence of programs like Mathematica, etc, why bother with integration and differentiation?  Why “numerical”?  The answer is that in science, one frequently encounters functions that are “crazy”, i.e., they do not behave in a regular and predictable manner and thus they are difficult to deal with analytically.  Most importantly, however, one needs frequently to integrate (or differentiate) numerical data that is the results of experiments.

 

Differentiation

 

Starting with the basic definition of the problem given in the figure below:

 

 

We can use the following simple difference formulas to compute the various derivatives.

 

 

One does have to careful in applying these formulas to take account of roundoff and truncation error.  This is especially true for the higher derivatives.  Generally speaking the symmetric derivatives are preferred over the forward or backward derivatives.  Another useful approximation to the derivative is the “5 point formula”, i.e.,

which is equivalent to using a 4th degree polynomial.  The following table indicates the errors associated with various methods of determining the derivative.

 

Some improved formulas exist for the higher derivatives

 

 

In some circumstances, like those depicted below, one may wish to smooth the derivatives  to get a more realistic picture of what is happening.  That is the case with the analysis of projectile motion.

 

 

Computer Methods

 

Excel

 

The following example shows the use of numerical derivatives in Excel and how they naturally occur in solving problems in analytical chemistry.

Click here to see example

 

 

 

 

Mathematica

 

In Mathematica, the function for numerical derivative is ND.  ND[f,x,x0] is the numerical derivative df/dx at x=x0.  ND[f,{x,n},x0] is the nth derivative

 

The following example shows how derivatives are taken of List data.

Click here to see example

 

FORTRAN90

 

 

Numerical Integration

 

There are several methods of numerical integration of varying accuracy and ease of use.  The most commonly used methods are the simplest, the trapezoidal rule and Simpson’s rule.  The following handwritten notes summarize some essential features of these methods.

 

 

 

 

 

 

Gauss Quadrature (unequally spaced points).

 

The Cadillac of numerical integration methods is that of Gauss quadrature.  Formally the value of the integral I is approximated as I = c1f(x1) + c2f(x2) + … cnf(xn)  where the ciare asset of predetermined numerical coefficients.  The following table shows these coefficients.

 

 

Computer Methods

 

Excel

 

The following examples shows the use of the trapezoidal rule and Simpson’s rule in a spreadsheet  The analytic value of the integral is 2.0000.

Click here to see example

 

 

 

 

 

 

 

 

 

 

Mathematica

 

The relevant Mathematica function for numerical integration is NIntegrate.  The following example shows the use of this Mathematica function.

Click here to see example

 

 

 

FORTRAN90

 

Thje following program, taken from Nyhoff (Fig 6-7.f) shows a FORTRAN90 implementation of the trapezoidal ruile.

 

      PROGRAM AREA

************************************************************************

* Program to approximate the integral of a function over the interval  *

* [A,B] using the trapezoidal method.  Identifiers used are:           *

*     A, B : the endpoints of the interval of integration              *

*     N    : the number of subintervals used                           *

*     I    : counter                                                   *

*     DELX : the length of the subintervals                            *

*     X    : a point of subdivision                                    *

*     Y    : the value of the function at X                            *

*     SUM  : the approximating sum                                     *

*     F    : the integrand                                             *

*                                                                      *

* Input:  A, B, and N                                                  *

* Output: Approximation to integral of F on [A, B]                     *

************************************************************************

 

      REAL F, A, B, DELX, X, Y, SUM

      INTEGER N, I

 

      PRINT *, 'ENTER THE INTERVAL ENDPOINTS AND THE # OF SUBINTERVALS'

      READ *, A, B, N

 

* Calculate subinterval length

* and initialize the approximating SUM and X

 

      DELX = (B - A) / REAL(N)

      X = A

      SUM = 0.0

 

* Now calculate and display the sum

 

      DO 10 I = 1, N - 1

         X = X + DELX

         Y = F(X)

         SUM = SUM + Y

 10   CONTINUE

      SUM = DELX * ((F(A) + F(B)) / 2.0 + SUM)

 

      PRINT 20, N, SUM

 20   FORMAT (1X, 'APPROXIMATE VALUE USING', I4,

     +            ' SUBINTERVALS IS', F10.5)

 

      END

 

 

** F(X) *********

* The integrand *

*****************

 

      FUNCTION F(X)

      REAL F, X

 

      F = X**2 + 1

 

      END

 

The following program implements Simpson’s rule in FORTRAN.

 

 

The following program implements Gaussian quadrature in FORTRAN.

 

 

Another technique of numerical integration, which we discuss in another lesson, is that of Monte Carlo integration.  We present the relevant FORTRAN code here for completeness.