**Contents:**

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.

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.

For sufficiently small h, the term in f''' should become negligable. This form is the

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.

#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-01The minimum error is for

We are seeking quadrature formulae of the type

where the

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

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

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.