Contents:

Functions on a Numerical Lattice

The theory of functions underlies virtually all scientific and engineering disciplines. A natural part of our introductory study of computational physics is to see how functions are represented numerically, and to understand which parts of the mathematical ``baggage'' we know about functions can be adapted for numerical applications. Here we look at the basic operations of numerical differentiation and integration (quadrature) of functions.

If a function f is known analytically, it is almost always possible, with enough patience and care (and possibly with the help of an algebraic software package, such as Maple) to differentiate it analytically. Integrating it, on the other hand, generally proves more challenging. However, sometimes the analytic expression for f(x) is too cumbersome to be useful; alternatively, f(x) may be the result of a numerical calculation, with no analytic formula at all. In such cases, f(x) is known only at particular values of the abcissa x--it will be defined on a numerical lattice.

Let's define a numerical lattice on which the function f(x) is defined, i.e., . To start with, we will use equally spaced lattice points, i.e., The position ``0'' will refer to the x position used in deriving our formulae.

f(x) on an equally-spaced lattice.

The basic analytical tool allowing the use of functions on a lattice is the Taylor Series expansion:


valid for small values of delta. On the lattice, this translates into

The basic mathematical operations on the lattice follow by interpolation of functions between the lattice points using simple mathematical expressions, such as low-order polynomials.

Differentiation

We seek formulae to compute the derivative of a function on a uniformly-spaced lattice, knowing the values of the function at the grid points. The starting point is the Taylor series expansion. For instance, subtract the two forms of eq. (2.1) at +h and -h to get the following form for the first derivative at x = 0:

For sufficiently small h, the term in f''' should become negligable. This form is the 3-point first derivative form.

If the function were a polynomial of order 2, the 3-point first derivative would be exact (since f''' = 0). Therefore, the geometrical interpretation of this formula is that the function f(x) is locally approximated by a polynomial of order 2--that is, a parabola.

A local approximation of f(x) by a straight line yields the less accurate ``forward'' and ``backward'' difference formulae (easily derived from eq. (2.1)):


Local approximations using higher-order polynomials yield more accurate, but also more complex, formulae. Some of these are listed in the accompanying
table of derivative forms. You may want to confirm these formulae starting from the Taylor series [eq. (2.1) and (2.2)].

Higher-order derivatives are derived similarly. For instance, adding the two forms of eq. (2.1) and subtracting twice f(0) yields the 3-point second derivative form:


Forms based on higher-order polynomials, and derived starting from eq. (2.1) and (2.3), are also listed in the table of derivative forms.

Simple Examples

The following program computes the sine function on a lattice, stores the results in an array, then computes numerically the first derivative of the function. The results are compared to the exact form. The reader should plot these data.
#include < math.h >

#define N_POINTS 300	              /* Array size */

void main()
{
    int i;
    float angle[N_POINTS], f[N_POINTS], f1_3[N_POINTS];
    float d_angle;
                                      /* compute the data     */

    d_angle = 2.0*M_PI / (float)(N_POINTS-1);

    for (i = 0; i < N_POINTS; i = i + 1 ) {
        angle[i] = i*d_angle;
        f[i] = sin(angle[i]);
    }
                                      /* 2-point derivative   */
                                      /* at beginning and end */
                                      /* of lattice           */

    f1_3[0] = (f[1] - f[0]) / d_angle;
    f1_3[N_POINTS-1] = (f[N_POINTS-1] - f[N_POINTS-2]) / d_angle;

                                      /* 3-point derivative   */
                                      /* everywhere else      */

    for (i = 1; i < N_POINTS-1; i = i + 1) {
        f1_3[i] = (f[i+1] - f[i-1]) / (2*d_angle);
    }

                                      /* output results and error */

    for (i = 0; i < N_POINTS; i = i + 1) {
        printf("%f %f %f %f\n ", angle[i], f[i], f1_3[i],
				 fabs(f1_3[i]-cos(angle[i])) );
    }

}

The error is plotted below.

Note that the error is relatively small. Note also that there is a systematic error introduced by the 3-point form which causes oscillations in the error function of same wavelength as those of the original (sine) function. In addition, there are random fluctuations due to the finite precision in the arithmetic. The latter could be diminished by using double type variables systematically, but not the former.

One might think that using the higher-order derivative forms would always be better; this is general the case, but with an important caveat! These forms may involve large cancellations in the numerator as h becomes small, thereby producing bad results due to round-off. The following code illustrates this point. It computes the derivative of sin(x) at 45 degrees.

#include < math.h >

void main()
{
    float x = M_PI/4.0, f1_3, f1_5, f1_e;
    float d_angle;
    int i;

    d_angle = 5.0;

    printf("\n\nTable of errors in 3- and 5-point derivatives\n\n");
    printf("   angle           3-point            5-point \n");

    for (i = 1; i < 10; i = i + 1) {

        d_angle = d_angle / 10.0;
        f1_3 = (sin(x+d_angle) - sin(x-d_angle)) / (2.0*d_angle);
        f1_5 = (sin(x-2*d_angle) - 8*sin(x-d_angle) 
		          + 8*sin(x+d_angle) - sin(x+2*d_angle))
		 / (12.0*d_angle);
        f1_e = cos(x);

        printf(" %e   %+12.8e   %+12.8e \n", d_angle, f1_3 - f1_e, f1_5 - f1_e);

    }

}

Running this code produces the following output:

 
Table of errors in 3- and 5-point derivatives 

   angle           3-point            5-point 
 5.000000e-01   -2.90966630e-02   -1.42991543e-03 
 5.000000e-02   -2.94446945e-04   +0.00000000e+00 
 5.000000e-03   -3.63588333e-06   -6.55651093e-07 
 5.000000e-04   +3.30805779e-05   +4.71472740e-05 
 5.000000e-05   +1.17421150e-04   +1.17421150e-04 
 4.999999e-06   +9.60350037e-04   +9.60350037e-04 
 4.999999e-07   -3.27571034e-02   -4.68060970e-02 
 4.999999e-08   +1.35830283e-01   +1.35830283e-01 
 5.000000e-09   -7.07106769e-01   -7.07106769e-01 

The minimum error is for d_angle = 0.005 for the 3-point formula, and 0.05 for the 5-point formula. The error becomes larger for smaller intervals due to the arithmetical errors. Note the use of single precision in the code to make the effect more evident!

Quadrature

Quadrature in its simplest form determines the 1-dimensional definite integral of f(x) over some finite range (a,b), where a < b. Eventually, we will consider more complicated cases, in which the limits may be infinite, the integrand may not be particularly well-behaved, or the number of dimensions may be greater than one.

We are seeking quadrature formulae of the type


where the fi are the values of the function at the lattice points and the Wi are the ``weights'' that insure that the finite sum over the lattice will give the value of the integral.

Consider a uniformly-spaced lattice of lattice spacing h; look at the two intervals (-h,0) and (0,h) and use (2.1):


This is equivalent to describing f(x) by straight-line segments in the two intervals. We can also use the three-point formula for f'(0) and f''(0) in eq. (2.1) to get

and substitute into the integral to obtain

Working formulae are obtained by summing these formulae over the domain of integration:


Note that the number of lattice points N = (b-a)/h must be even to use the Simpson formula.

You can improve the accuracy of these procedures by using higher-order polynomials. Using cubic and quartic polynomials yield, respectively, Simpson's 3/8 and Bode's rules:


For the simple and well-behaved functions we will consider in this course, these forms will suffice.