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.




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




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







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





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




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













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

Click here to see example






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




* 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



      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


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


      PRINT 20, N, SUM


     +            ' SUBINTERVALS IS', F10.5)





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

* The integrand *




      REAL F, X


      F = X**2 + 1




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.