2. Finite Difference Techniques

To numerically solve PDEs, it is necessary to discretize the continuous differential equations. To this end, the space-time domain on which the equations are solved is replaced with a discrete grid of points, the continuous functions entering the PDEs are replaced with grid functions that take values on the points of the grid and the derivative operators are replaced by finite difference operators. For the four-dimensional wave equation, the spatial domain [ x_min ; x_max ] × [ y_min ; y_max ] × [ z_min ; z_max ] is discretized as N_x × N_y × N_z points and the time domain [ t_min ; t_max ] is discretized as n time levels, for a total of N_x × N_y × N_z × n points. N_x, N_y, N_z and n are defined in terms of dx, dy, dz and dt as (summation convention suppressed):

In general, in finite difference techniques dx, dy and dz are chosen equal, and the ratio between dxⁱ and dt, called the Courant factor &lambda, is kept constant, so that dt alone is sufficient to specify the discretization level along all axes.

The finite difference operators can be determined by the following general procedure, which will be exemplified for operator D_x which corresponds to the first derivative with respect to x. As explained in these lecture notes by Matt Choptuik, for any grid function u defined at grid point ( i, j, k), it is possible to expand u( i+1, j, k) and u( i-1, j, k) as a Taylor series around ( i, j, k) as:

Solving for D_x from these two equations gives the operator corresponding to the first derivative with respect to x as:

This formula is correct to second order in dx. In a similar manner, the operator corresponding to the second derivative with respect to x, also accurate to second order in dx, can be determined to be:

It should be noted that nesting operators also provides a correct method for determining higher-order operators, although in the case of second derivatives nesting the above expression for D_x gives a 5-point stencil rather than a 3-point stencil. However, nesting operators can be used to provide working formulas for mixed partial derivative operators.

There is a complication worth mentioning regarding finite difference operators. Although the continuous wave equation is valid in the entire spatial domain, operator expressions need not be. For instance, the above expression for D_x obviously cannot be valid for points with i = 1, as it references i - 1. When such cases occur it is necessary to divide the spatial domain into subdomains, and deduce the correct operator formulas for each domain. Particular care should be given to this caveat, as incorrect use of the finite difference operators is a prime cause of segfaults.

For more information on finite difference techniques, see the lecture notes by Matt Choptuik, given at the Taller de Verano 1999 summer school.

Some of the material presented on this web page is based upon work supported in part by the Princeton Applied and Computational Mathematics Program and by the National Science Foundation under Grant No. 0745779.