15

The Numerical Solution of Elliptic and Parabolic Partial Differential Equations

DAVID YOUNG

PROFESSOR OF MATHEMATICS AND DIRECTOR OF THE COMPUTATION CENTER
THE UNIVERSITY OF TEXAS

15.1    Introduction

Many important problems in science and engineering can be formulated in terms of linear second-order partial differential equations. For instance, a number of problems in elasticity,⁴⁴ steady-state heat conduction,^8,39 and fluid mechanics^8,39 involve the solution of Laplace’s equation in a two-dimensional region subject to certain auxiliary conditions at the boundary of the region. Sometimes the function itself is prescribed on the boundary; in other cases, the normal derivative or even a linear combination of the function and its normal derivative may be prescribed. Some problems—for instance, those involving steady-state heat conduction with variable conductivity—involve linear second-order elliptic partial differential equations other than Laplace’s. Other problems, such as those arising in the study of the neutron-flux distribution in a nuclear reactor, may involve the simultaneous solution of several equations.

Except in very special cases, it is not practical to obtain numerical solutions to such problems by analytic methods. Instead, one is nearly always forced to resort to approximate numerical methods. Of such methods, probably the most popular and the one that is most frequently used is the method of finite differences.⁸ A network of regularly spaced straight lines, each parallel to one of the coordinate axes, is superimposed over the region. Instead of seeking the value of the solution of the differential equation at every point in the continuous region, one seeks to determine approximate values of the solution at the mesh points, the intersections of the network lines with each other and with the boundary of the region. The partial differential equation is replaced by a partial difference equation. Although more elaborate techniques may be used, one method of obtaining the difference equation at a given interior mesh point is to replace the partial derivatives in the differential equation by partial difference quotients in the standard way, thus obtaining a linear algebraic equation corresponding to each interior mesh point. In the case of a boundary-value problem, since the values at the boundary mesh points are prescribed, the solution of the difference equation involves the solution of a system of N linear algebraic equations with N unknowns, N being the number of interior mesh points.

While it is easy to show that under rather general conditions the system of equations has a unique solution, the actual numerical determination of the solution presents a serious computational problem. Since very little is known about the accuracy of the solution of the difference equation as a solution of the differential equation, it is desirable to make the distance between adjacent mesh points, here denoted as the mesh size, small in order to be reasonably sure of obtaining satisfactory accuracy. Thus the number of interior mesh points and hence the number of equations will generally be large; in some cases it may be in the thousands. Direct methods such as the use of determinants, as in Cramer’s rule, and the Gauss elimination method appear out of the question. Since the matrix of the system is very sparse, each equation involving the unknown values only at the corresponding mesh point and at a few neighboring points, the use of iterative methods appears appropriate. One selects an arbitrary initial approximation to each of the unknowns and successively modifies the approximate values according to a given rule of procedure. If the iterative method is properly chosen, the approximate values for any given point will converge to the value of the exact solution of the difference equation at that point. Of course, what is important from the computational standpoint is not only whether the method converges but also how fast it converges.

The Gauss-Seidel method is probably one of the simplest iterative methods and is very easy to use on a high-speed computer; however, although convergence is guaranteed for a wide class of problems, in many cases the convergence is exceedingly slow. For instance, for Laplace’s equation the number of complete iterations—i.e., the number of times that the value of the function at each point is modified—necessary to achieve a specified degree of convergence is approximately proportional to h⁻², where h is the mesh size. As h decreases, the number of iterations increases very rapidly. By a simple modification of the Gauss-Seidel method, the number of iterations can be substantially reduced. The modified method, which was called the “extrapolated Liebmann method” by Frankel¹⁷ and the “successive-overrelaxation method” by Young,⁵⁷ involves a parameter known as the relaxation factor. If this relaxation factor is properly chosen, then in the case involving Laplace’s equation discussed above the number of iterations varies approximately as h⁻¹, thus affording a very substantial improvement over the Gauss-Seidel method.

Another iterative method that is designed to give rapid convergence is the alternating-direction implicit method of Peaceman and Rachford.³⁵ Whereas both the Gauss-Seidel method and the successive-overrelaxation method are “point” iterative procedures in the sense that at each step the value at one single point is modified, in the Peaceman-Rachford method the values for all points of a line of mesh points, either a row or a column, are modified simultaneously. The modifications are carried out alternately first for all rows and then for all columns. Although in order to perform a row or a column iteration it is necessary to solve a linear system with as many equations and unknowns as there are mesh points in the row or column, nevertheless since the matrix of the system is usually tridiagonal, the solution of the system presents very little practical difficulty. If a set of iteration parameters is properly chosen, then for the case of Laplace’s equation on a rectangle the number of complete iterations required to achieve convergence varies approximately as |log h|⁻¹. The theoretical improvement over the successive-overrelaxation method in this case is considerable. Moreover, although the theoretical advantage in convergence has so far been proved only in the special case mentioned above, nevertheless actual numerical experiences indicate a much greater generality. At the present time, the successive-overrelaxation method and the Peaceman-Rachford method are each finding considerable usage for problems involving elliptic partial differential equations while other methods are being used only occasionally. Quite recently, successive-line overrelaxation and even two-line over-relaxation has been used by Varga and others.¹¹

Parabolic partial differential equations arising in scientific and engineering problems are often of the form

U_t = L(U)

where L(U) is a second-order elliptic partial differential operator with either one or two independent “space” variables. Typical examples are the equation

U_t = U_xx

involving one space variable, and the equation

U_t = U_xx + U_yy

involving two space variables. Let R denote the region of the space variable or variables and let S denote the boundary of R. One frequently seeks to solve an initial-value problem in which the values of U in R are specified for t = t₀ and the values of U on S are specified for all t > t₀. As in the case of problems involving elliptic equations, the method of finite differences is frequently used to obtain approximate numerical solutions. In addition to the network that is superimposed over the region R, a time increment k is introduced. One seeks approximate numerical values of U for all mesh points in R and for t = t₀ + k, t₀ + 2k, …, etc., the values of U for t = t₀ being already given. Here, in contrast to the situation for elliptic equations, it is possible to develop an explicit procedure for determining the values of U for any given value of t in terms of those values that have already been obtained for previous values of t. Unfortunately, however, this procedure, which is sometimes known as the “forward-difference method,” may be numerically unstable unless k is chosen so small relative to the space-mesh size that the computational labor involved in proceeding from t₀ to another value of t, say t₁, in steps of length k may be prohibitive even with a fast computer. A larger value of k could not be used because the excessive growth of rounding errors that occurs with an unstable process would destroy all accuracy within a few time steps.

For problems involving one space dimension, this difficulty may be avoided by use of the Crank-Nicolson method.¹⁰ The procedure is implicit in that the process of proceeding from one time step to another involves the solution of a system of linear algebraic equations, there being as many equations as there are interior space-mesh points. Since the matrix of the system is tridiagonal, however, there is little difficulty involved. The extension of the Crank-Nicolson method to problems involving two space dimensions would involve solving an elliptic boundary-value problem for each time step. Since this procedure would be very slow indeed, the alternating-direction method of Peaceman, Rachford, and Douglas is recommended. This method involves computing approximate values of U successively at t₀ + k/2, t₀ + k, t₀ + 3k/2, etc., using a procedure that is alternately implicit on rows and explicit on columns and vice versa. Since each step involves the solution of a linear system with a tridiagonal matrix, the computing time is much less than would be required if the Crank-Nicolson method were used. The local accuracy of the method is of the same order of magnitude as for the Crank-Nicolson method. The alternating-direction method is also numerically stable.

Prior to the advent of high-speed computing machines, numerical solutions were obtained by manual methods. For boundary-value problems involving equations of elliptic type, the difference equations were solved by “relaxation methods”; see, for instance, the books of Southwell⁴⁵ and Allen.² The relaxation method is basically an iterative method in which, instead of following a fixed procedure, the computer is allowed considerable freedom in modifying the approximate values. A skilled and experienced computer could often obtain a numerical solution to a complicated problem in a remarkably short time. Nevertheless, the amount of drudgery involved with manual methods places a limit on the accuracy, on the speed with which a solution can be obtained, and on the number of related cases that it is practical to solve in a given study.

The availability of modern high-speed machines and the development of efficient and effective numerical procedures have made practical the solution of a number of problems. In Sec. 15.9, the solutions of three such problems are described; these include a flow problem, a nuclear-reactor problem, and a thermal-ignition problem. It would be interesting to try to estimate how many millions of hours of tedious hand computing have been avoided by the use of these programs.

One difficulty in the use of large computers for solving partial differential equations is the amount of effort required in the preparation of a program for the machine. This is of less relative importance when a great many cases are to be solved, but if only a few cases are required, then it may not always be easy to decide whether or not the problem should be done on a computer. In any case, there is a considerable time lag between the formulation of the problem and its solution on the computer.

Attempting to remedy this situation, a group of computing installations, primarily on the West Coast, established a collaborative project known as the SPADE project. The objective of the SPADE project was the development of a set of programs for the IBM 704 and 709 computers for solving a class of problems involving systems of elliptic and parabolic partial differential equations. It would be necessary for the program user to specify the boundary (for instance, by giving the equations of a finite number of boundary segments), the mesh size, formulas for the determination of the coefficients of the differential equations, and the iterative method to be used. The program of SPADE would then determine the coordinates of the mesh points, determine the coefficients of the difference equations, and solve the difference equations by using the prescribed iterative procedure. Since the programs were to be available as closed subroutines that might be used in various combinations, a considerable degree of flexibility would be afforded. By means of the SPADE programs, the amount of labor and the time involved in preparing a problem for solution on a high-speed computer would be greatly reduced. Actually, in fact, the SPADE project has recently been abandoned. Its failure can be attributed to lack of sufficient effort and does not imply that the program could not be successfully carried out.

15.2    Boundary-value Problems and the Method of Finite Differences

Let R be a bounded plane region with boundary S and let f be a function defined and continuous on S. We consider the problem of finding a function U that is continuous in R + S, is twice differentiable in R, satisfies in R the linear second-order partial differential equation

and satisfies on S the condition

Here A, C, D, E, F, and G are analytic functions of the independent variables x and y in R + S and satisfy the conditions A > 0, C > 0, and F ≤ 0. Because of the conditions on A and C, Eq. (15.1) is said to be elliptic.^8,39 More generally, an equation of the form

is said to be elliptic in R + S if B² − AC < 0 in R + S. We remark that by introducing new variables ξ and η satisfying the conditions

one obtains an equation of the form

where α > 0, γ > 0 in R′ + S′, the region in the ξη plane corresponding to R + S. Equation (15.4) is similar to Eq. (15.1) in so far as the terms involving the second-order derivatives are concerned. Much of our subsequent discussion will apply to Eq. (15.4) as well as to Eq. (15.1) and will also apply to certain cases in which, instead of prescribing U on S, one prescribes on S the normal derivative ∂U/∂n or a linear combination of U and ∂U/∂n.

A special case of the above is the solution of Laplace’s equation

in the circle x² + y² < 1 with values prescribed on the boundary of the circle. Such a problem corresponds to the solution of a steady-state heat-conduction problem for an infinite right-circular cylinder with axis normal to the xy plane and with external temperature held constant along lines parallel to the axis of the cylinder. The variation of the temperature on the circumference of any cross section of the cylinder is prescribed. If, instead of prescribing the temperature everywhere, one assumes that part of the surface is insulated, then the corresponding condition for the appropriate part of the circle x² + y² = 1 is ∂U/∂n = 0.

In order to apply the method of finite differences to obtain an approximate numerical solution of the problem defined by Eqs. (15.1) and (15.2), the first step is to construct a network. We choose an arbitrary but fixed point () and a positive number h that we shall call the mesh size. We choose one mesh size for simplicity, although we could have different mesh sizes in the two coordinate directions. We let ∑_h denote the set of points (x, y) such that both and are integers. Such points are called mesh points. Two mesh points (x,y) and (x′, y′) are adjacent if

(x − x′)² + (y − y′)² = h²

We let R_h denote the set of all points that are in both R and Σ_h. Points of R_h are known as interior mesh points. We determine the set S_h of boundary mesh points as follows: For each point (x, y) of R_h, we consider the four adjacent mesh points. Let (x′, y′) denote one such point. If no point of S lies on the line segment joining (x,y) to (x′,y′), then no point of S_h lies on that segment. Otherwise, if at least one point of S lies on the line segment, then the nearest such point to (x,y), which may be (x′,y′), is said to belong to S_h. The set S_h consists of all points found by considering each point of R_h and the corresponding four line segments to the four adjacent points.

As an example, consider the ellipse

with and h = 1; see Fig. 15.1. Points of R_h and S_h are designated by solid black circles and by open circles, respectively. Other mesh points are designated by open squares.

Fig. 15.1 Interior and boundary mesh points.

Once we have determined a set of mesh points, our objective will be to determine approximate values of U at these points. To this end, we next construct a difference-equation representation of Eq. (15.1). First we consider a typical regular point (x₀,y₀) of R_h, that is, a point of R_h such that the four adjacent points of Σ_h belong either to R_h or to S_h. We seek to determine a difference equation corresponding to Eq. (15.1) for the point (x₀,y₀).

Fig. 15.2 Typical neighborhood configuration for a regular interior mesh point.

Let us consider the configuration shown in Fig. 15.2; this consists of (x₀,y₀) and the four adjacent mesh points, which for convenience are labeled as shown. We first seek an expression for U_xx. We adopt the notation that U_i = U(x_i,y_i), i = 0, 1, 2, 3, 4, and (U_x)_i = U_x(x_i,y_i), etc. Assuming that U (x,y) has partial derivatives of the fourth order in a sufficiently large neighborhood about (x₀,y₀), we have, by Taylor’s theorem,

with

where ξ is some value of x such that

x₀ < ξ < x₁

Similarly,

By combining Eqs. (15.6) and (15.7), we obtain

By using the same method, we can also obtain

By neglecting the remainder terms, i.e., the terms involving brackets, in Eqs. (15.8) to (15.11) and by substituting in Eq. (15.1) for each regular point (x,y) in R_h, we obtain the difference equation

where A, C, D, E, F, and G are evaluated at (x,y). If we write Eq. (15.12) in the form

where the α_i, i = 0, 1, 2, 3, 4, are functions of x and y, we observe that, because A > 0, C > 0, and F ≤ 0, if the α_i are positive then

Fig. 15.3 Typical neighborhood configuration for an irregular interior mesh point.

Even for irregular interior mesh points, i.e., points of R_h that are not regular interior mesh points, we can write the difference equation in the form (15.13). Thus if, as in Fig. 15.3, x₁ − x₀ = s₁h and x₃ − x₀ = s₃h, where 0 < s₁ ≤ 1, 0 < s₃ ≤ 1, with the aid of Taylor’s theorem we can easily derive the following formulas:

Similar expressions can be obtained for (U_y)₀ and (U_yy)₀. If the resulting expressions are substituted into Eq. (15.1), then one obtains

where

It is evident that every α_i appearing in Eqs. (15.13) and (15.17) will be positive provided h is chosen so small that

h < min (M₁, M₂)

where

For M₁ and M₂, the minima are taken over all points of R + S. That M₁ and M₂ are positive is easily seen by nothing that the functions A, C, D, and E are continuous in the closed region R + S and that A > 0, C > 0. It is evident that D and E are bounded and that A and C have positive minima.

If Eq. (15.1) is self-adjoint, i.e., if it can be written in the form

then one can obtain a difference equation with the following properties: The local accuracy is as good as that for Eqs. (15.13) and (15.17); the coefficients are positive; and if R_h contains only regular points, then the difference equation is symmetric. That is to say, the coefficient of U(x′,y′) in the difference equation for the point (x,y) is the same as the coefficient of U(x,y) in the equation for the point (x′,y′). In such a case, the matrix of the linear system corresponding to the difference equation would be symmetric, a situation that has certain theoretical and, perhaps, practical advantages.

In order for Eq. (15.1) to be self-adjoint, it is necessary and sufficient that we have

If Eqs. (15.20) hold, then the differential equation can readily be put in the self-adjoint form. In deriving the difference equation for the neighborhood configuration shown in Fig. 15.3, we use finite-difference expressions of the form

and derive a difference equation of the form (15.17) with

We remark that, if there are irregular interior mesh points, then the matrix of the linear system that one obtains will not in general be symmetric. So far, experience indicates that this is probably not a serious disadvantage. We further remark that there are many other possible difference equations as well as other methods for deriving them, but space does not permit us to discuss them here.

Here and subsequently we assume that our difference equation can be written in the form (15.17), where the α_i are positive and where

We further assume that the region R_h is connected; i.e., any two points of R_h can be joined by a broken line consisting of line segments connecting pairs of adjacent points of R_h.

Under the foregoing assumptions, it is not difficult to show that Eq. (15.17) has a unique solution. We shall give here a well-known proof. It is clearly sufficient to show that the determinant of the matrix of the related system of linear equations does not vanish. This will be true provided that the homogeneous system, which we obtain by letting the function G and the values of U on the boundary of the region vanish, has only the trivial solution that vanishes throughout R_h. Suppose that the homogeneous system had a nontrivial solution V. Then for some point of R_h we would have V ≠ 0. We can assume that V > 0 at this point; otherwise, we could consider the function −V, which would also satisfy the homogeneous system. Let M > 0 denote the maximum value of V in R_h, and let (x₀,y₀) be the interior mesh point for which V(x₀,y₀) = M. By Eqs. (15.17) and (15.22), and the fact that the α_i are positive, it follows that V must assume the value M at each of the four adjacent points. By repeating this argument and using the connectedness of R_h, one can show that V(x,y) = M for all (x,y) in R_h and in S_h. This contradicts the assumption that V vanishes on S_h and thus proves that V vanishes identically in R_h + S_h. Consequently, the homogeneous system admits only the trivial solution. It then follows that Eq. (15.17) has a unique solution.

Before considering, as we shall in the next three sections, the formidable computational problem of finding the unique solution that we now know exists, we shall discuss the accuracy of the exact solution of the difference equation as a solution of the differential equation. It is shown in the well-known paper of Courant, Friedrichs, and Lewy⁹ that, as the mesh size h approaches zero, the solutions of the corresponding difference equations approach the solution of the differential equation; nevertheless, very little is known about the accuracy of the solution of the difference equation for a given value of h.

Besides the doubtful methods of intuition and guessing, one frequently used procedure is to solve the difference equation for some mesh size h and then again for the mesh size h/2. The difference between the two solutions so obtained is taken as a measure of the error. Indeed, if one assumes that the error is precisely proportional to h², which assumption is fairly well satisfied in certain favorable cases (see, for instance, Ref. 59), then one can do even better by using a process that was introduced by L. F. Richardson³⁶ and that is referred to as “extrapolation to zero grid size.” If u_h denotes the solution of the difference equation with mesh size h and if U denotes the exact solution of the differential equation, then under the given assumptions we have

Although the theoretical basis for using this procedure is not very substantial, in some cases one obtains a remarkable improvement in accuracy by using it.

Some empirical studies on the accuracy of difference-equation solutions for Laplace’s equation indicate that extrapolation to zero grid size might work rather well provided that the boundary has no corners with interior angle greater than 180° and provided that the boundary values are sufficiently smooth. The effect of curved boundaries, which would involve irregular mesh configurations, could be very unfavorable.

Gerschgorin²¹ obtained a bound for the maximum error in terms of the local error. In the case of Laplace’s equation for a region containing only regular interior points, Gerschgorin showed that, for any point (x,y) in R_h,

where r is the radius of any circle containing the given region and where

M₄ = max (X₄, Y₄)

with

X₄ = max |U_xxxx|Y₄ = max |U_yyyy|

the maxima for X₄ and Y₄ being taken in R + S. Since the solution U is not known, much less its fourth partial derivatives, this estimate has limited usefulness, particularly since in some important cases the fourth partial derivatives are not even bounded. In some cases, however, one can estimate U_xxyy by using the difference quotient

and then estimate U_xxxx and U_yyyy by the relations

U_xxyy = −U_xxxx = −U_yyyy

which are derived from Laplace’s equation. When the largest derivative thus found is substituted in Gerschgorin’s error formula, we obtain an estimate of the accuracy.

Although a great deal of effort has been and is being spent in an attempt to find a satisfactory method of error estimation (see, for instance, papers by Rosenbloom,³⁸ Wasow,^55,56 Walsh and Young,^53,54 and Laasonen^27,28), a great deal of research remains to be done in this direction.

15.3    Point Iterative Methods

Although the existence of a unique solution of the difference equation (15.17) under rather general conditions has already been established, nevertheless the problem of actually computing the solution presents a serious practical difficulty. In view of the uncertainty about the accuracy of the finite-difference solution, one is inclined to wish to use a small value of the mesh size h. This, of course, implies that the number of mesh points will be large. In many practical cases, several thousand mesh points are considered. While the use of direct methods—such as Cramer’s rule, which involves determinants, and the Gauss elimination method—is almost out of the question, the fact that the matrix of the system is sparse suggests the use of iterative methods. We note from Eq. (15.17) that at most five elements in any row of the matrix are nonzero. The sparseness of the matrix is of considerable advantage for the use of iterative methods but does not help appreciably for direct methods.

The use of an iterative method involves first the selection of a function u⁽⁰⁾(x,y), which is arbitrary except that the closer it agrees with u(x,y) the less computing effort will be required. One then modifies the values of u⁽⁰⁾(x,y) at the interior mesh points according to a prescribed procedure in such a way that the approximate values so obtained converge to u(x,y). In this section, we consider “point” iterative methods for which at each stage the approximate value at one single point is modified. In the next section we shall consider iteration procedures such that at each step several values are modified simultaneously.

For convenience, we let u(x,y) denote the exact solution of Eq. (15.17). If we solve Eq. (15.17) for u(x,y), we get

where

and where, for i = 1, 2, 3, 4,

We note that every β_i is positive and that

We remark that, in solving a problem on a high-speed computer, one would normally evaluate all of the coefficients β_i and then store them in an auxiliary memory such as a magnetic drum or magnetic tape memory. Large blocks of β_i could then be transferred to the high-speed memory when needed. Of course, if many of the coefficients are identical, as for instance in the case of Laplace’s equation where for every regular point each β_i equals ¼, then some other procedure would be more efficient.

Perhaps the simplest point iterative method is the Gauss-Seidel method, which is also known as the method of “successive displacements.”²⁰ The approximate values of u are improved at points of R_h according to an arbitrary but fixed ordering, using Eq. (15.23) and using improved values as soon as available. Thus, if we consider the ordering in which (x,y) follows (x′, y′) provided that either y > y′ or that y = y′ and x > x′, then for n = 1, 2, 3, … the successive approximations are determined by

A complete iteration consists in improving the approximate values once at all points of R_h. Having traversed the points of R_h, one starts over again at the “first” point and repeats the process until d_n < , where is the prescribed tolerance and where

d_n = max |u⁽ⁿ⁾(x,y) − u⁽ⁿ⁻¹⁾(x,y)|

the maximum being taken over all points of R_h.

We remark that in a similar method known as the “Jacobi method,” or the “method of simultaneous displacements,”²⁰ one does not use improved values until after a complete iteration. The improvement formula for this method is

This method requires approximately twice as many iterations to attain a specified convergence as does the Gauss-Seidel method.

As an example of the use of the Gauss-Seidel method, let us consider the problem of solving Laplace’s equation

U_xx + U_yy = 0

for the unit square 0 ≤ x ≤ 1, 0 ≤ y ≤ 1. Let us choose M = h⁻¹ = 3. For boundary conditions, we let U = 1,000 on the side of the square where y = 1, and we let U = 0 elsewhere on the boundary. Evidently, each of the four points of R_h is a regular interior mesh point. By Eqs. (15.17), (15.18), and (15.23), the difference equation for each point of R_h becomes

Table 15.1 illustrates the computational procedure. First, the initial approximation u⁽⁰⁾(x,y) to the solution of the difference equation was set equal to zero at the four interior mesh points. The ordering of the points was as follows: (⅓,⅓), (⅔,⅓), (⅓, ⅔) (⅔,⅔).

Since the values of u⁽⁶⁾ are identical with those of u⁽⁷⁾, the iteration process was terminated after seven iterations. The final values represent the exact solution of the difference equation, as can be verified by substituting in Eq. (15.26). For convenience in studying the behavior of the method, the values of d⁽ⁿ⁾ = u⁽ⁿ⁾ − u⁽ⁿ⁻¹⁾ are given in Table 15.1. The ratios d⁽ⁿ⁾/d⁽ⁿ⁻¹⁾ give an indication of how rapidly the errors are decreasing. In this case it can be seen that the ratios are approximately ¼ hence, the convergence is quite rapid. The general formula for the limiting ratio λ, however, is

See, for instance, the paper by Young.⁵⁷ We shall refer to λ as the spectral radius† of the linear transformation defined by the Gauss-Seidel method.

Table 15.1 An Example of the Use of the Gauss-Seidel Method

That the convergence of the Gauss-Seidel method is very slow can be seen from considering the case M = 20. Since λ = 0.975 in this case, the error in the approximate solution is reduced by only about 2.5 per cent in each iteration.

The number N of iterations needed to reduce the error to a specified factor ρ of its original value is given approximately by

Hence,

By this argument, it can be seen that the number of iterations varies as M² = h⁻². For a more rigorous discussion, the reader might consult Ref. 57.

By a simple modification of Eq. (15.25), a substantial improvement can be made in the rate of convergence. The modified formula is

Here ω is a parameter that we shall refer to as the relaxation factor. If ω = 1, then the procedure reduces to the Gauss-Seidel method. The proper choice of the relaxation factor determines the rate of convergence of the method.

The method defined by Eq. (15.28) is referred to as the “successive overrelaxation method” by Young;⁵⁷ when applied to Laplace’s equation, it is called the “extrapolated Liebmann method” by Frankel.¹⁷ We remark that the ordinary “Liebmann method” is just the Gauss-Seidel method as applied to the difference-equation analogue of Laplace’s equation; see the papers by Liebmann³³ and Shortley and Weller.⁴²

The use of the successive-overrelaxation method will be illustrated by applying it to the same numerical example as was used in connection with the Gauss-Seidel method. Table 15.2 shows the results obtained using 1.1 as relaxation factor. It should be noted that, since the values of u⁽³⁾ and u⁽⁴⁾ are identical, the process is terminated after four iterations. The final values differ from the exact solution of Eq. (15.28) by one unit. This discrepancy could have been avoided if more significant figures had been retained throughout the iteration process. The erratic behavior of the ratios d⁽ⁿ⁾/d⁽ⁿ⁻¹) should be noted and contrasted to the relatively smooth behavior for the Gauss-Seidel method.

Table 15.2 An Example of the Use of the Successive-overrelaxation Method

As shown in Ref. 57, the optimum value of ω is given by

where λ is the spectral radius of the Gauss-Seidel method. In the example treated above, we have λ = ¼, and upon substituting this value in Eq. (15.29), we get ω_b = 1.072. It is not difficult to show that the convergence of the successive-overrelaxation method with the best value of ω is approximately as rapid as it would be if the ratio d_n+1/d_n tended to λ*, where

Here λ* is the spectral radius of the successive-overrelaxation method. In the case M = 20, we get ω_b = 1.73 and hence λ* = 0.73, as compared with λ = 0.975 for the Gauss-Seidel method. Thus, on the average, one complete iteration of the successive-overrelaxation method reduces the error by approximately 27 per cent. The number N of iterations required to reduce the error to a specified fraction ρ of its initial value is given approximately by

Hence

from which it follows that the number of iterations varies approximately as M = h⁻¹. This represents a very substantial improvement over the Gauss-Seidel method and one that has been well verified in actual practice; see, for instance, Ref. 59.

In order to achieve the maximum possible improvement in the speed of convergence, it is necessary to select the proper value of ω. By Eq. (15.29), this choice depends on an estimate of λ, the spectral radius of the Gauss-Seidel method. For the case of Laplace’s equation and for a rectangle with sides a = Rh and b = Sh, where R and S are integers, λ is given by

For other problems involving Laplace’s equation but for regions other than rectangles, one may compute λ for a rectangle that contains the given region and then use the value of ω computed from Eq. (15.29). The value of ω so obtained will be somewhat larger than the true optimum value. It has been shown theoretically and verified experimentally by Young^57,59 that an overestimation of the optimum value of ω does not result in a serious reduction in the rate of convergence. A slightly better estimate of the optimum value could probably be obtained by the use of a rectangle of approximately the same area and proportions.

The problem of choosing a suitable relaxation factor is in general more difficult for problems involving differential equations other than Laplace’s. Sometimes special properties of the equation can be used. It may be possible to show that the appropriate value of λ will be close to that for Laplace’s equation. Alternatively, we can estimate λ by performing a number of iterations with ω = 1 and then estimating the average of the ratios d⁽ⁿ⁾(x,y)/d⁽ⁿ⁻¹⁾(x, y). Finally, and this is certainly a reasonable approach if a large number of cases are to be treated, one can simply try different values of ω for different cases and determine the best value by experiment.

The improvement in convergence rate that is afforded by the successive-overrelaxation method has been shown by Young⁵⁷ to hold under certain conditions. These conditions are satisfied by symmetric five-point difference equations of the form (15.17), provided that the ordering of the points is properly chosen. Recent theoretical work of Garabedian¹⁹ and Kahan,²⁶ however, indicates that the improvement in convergence rate can be realized under much more general conditions. For instance, it was observed by Young⁵⁹ in certain numerical cases that improvement in convergence rate occurs when a nine-point difference equation is used to represent Laplace’s equation. Garabedian¹⁹ provided theoretical justification for this. Kahan’s results indicate that the improvement in convergence rate will hold provided that the coefficients α_i; of the difference equation when written in the form (15.17), though possibly with more terms, are positive, and provided that the matrix of the corresponding linear system is symmetric and positive definite. The ordering according to which the approximate values are modified can be arbitrary.

15.4    Peaceman-Rachford Iterative Method

Whereas the Gauss-Seidel and the successive-overrelaxation iterative methods involve at each step the modification of the approximate solution at one single point of R_h, the method that we shall consider in this section involves the simultaneous modification of the values at all points of a line of points. Such a line of points may be either a row of points or a column of points. The method was developed by Peaceman and Rachford³⁵ and is frequently known as the Peaceman-Rachford method. Each iteration consists of two half iterations: first, all the values on all the rows are modified; second, all the values on all the columns are modified. To illustrate, let us consider the five-point difference-equation representation of Laplace’s equation for a region containing only regular interior mesh points. Our difference equation (15.13) becomes

As before, one selects an arbitrary function u⁽⁰⁾(x,y) defined at all points of R_h. In general, having determined u⁽ⁿ⁾(x,y), we then determine u^(n+½)(x,y) by a row iteration and u⁽ⁿ⁺¹⁾(x,y) by a column iteration according to the following formulas:

Here r is a parameter that may depend on n but must be the same for both halves of any iteration.

In order to determine u^(n+½) (x,y) from the first of the above equations, it is necessary to solve, for each value of y, a system of M simultaneous linear equations, where M is the number of interior mesh points on the row. It is easy to see that the problem of doing this is the same as that of solving a system of M linear algebraic equations of the form

where the T_i are unknown and where the A_i, B_i, C_i, and D_i are given coefficients. Since the matrix of such a system is tridiagonal, i.e., has no nonzero elements except on the main diagonal and on the diagonals adjacent to the main diagonal, there is a very simple algorithm for solving the equations. In this algorithm, which is actually just the Gauss elimination method, the number of arithmetic operations is proportional to M, as opposed to approximately ⅓M³ operations as required with a general matrix. The algorithm appears to have been first used for solving elliptic partial differential equations by L. H. Thomas.⁴⁷ It was used first in connection with parabolic partial differential equations by Bruce, Peaceman, Rachford, and Rice.⁶

The following formulas are used in the determination of the T_i:

The numerical procedure defined by Eqs. (15.33) has been observed to be very satisfactory in a large number of cases, no large growth of rounding errors having been noted.

Peaceman and Rachford³⁵ and Wachspress and Habetler⁵² have given methods for choosing iteration parameters. We first describe a procedure that is very similar to that given in Ref. 35. Let M − 1 denote the number of interior mesh points in the row or column of greatest length and let

We choose a set of m different values of r that will be used in a cyclic order: r₁, r₂, …, r_m, r₁, r₂, …. First, m is chosen as the smallest positive integer such that

The r_k are then determined by

for k = 1, 2, …, m. After each set of m iterations, the error will be reduced by a factor

As an example, let us consider the case of the unit square with . Since M = 20 we have, by Eq. (15.34),

From inequality (15.35), we find that m = 3; hence by Eq. (15.36) we have

Thus with three double sweeps, the error will theoretically be reduced by a factor of 0.16. This may be compared with a factor of reduction of 0.73 for the same problem after one iteration of the successive-overrelaxation method and with 0.39 after three iterations. Because of the added complication involved in using the Peaceman-Rachford method and the fact that each iteration involves a double sweep, it would appear that the successive-overrelaxation method is slightly better for the particular problem under discussion. Actually, however, the reduction factor for the Peaceman-Rachford method is somewhat less than that indicated above. Moreover, for smaller values of h the Peaceman-Rachford method is relatively better, the number of iterations being approximately proportional to |log h|⁻¹ instead of to h⁻¹ as for the successive-overrelaxation method; see the discussion given in Ref. 35.

An alternative procedure for choosing the r_k, based on more detailed analysis, has been given by Wachspress and Habetler.⁵² The method described below represents a slight variant of their procedure. Again letting M − 1 denote the largest number of interior mesh points in any row or column, one first computes a number z from the relation

z = σ^1/m−1

where

The r_k are then determined by

According to the analysis given in Ref. 51, the factor of reduction of the error that is achieved after m double sweeps is approximately

Since the average factor of reduction per double sweep is S_m = (P_m)^1/m, it seems reasonable to choose m so as to minimize S_m. Because of the complicated formula for S_m as a function of m, it is probably simplest to determine the optimum value of m by computing S_m for a number of different values of m.

For example, in the example considered above we have

By computing S_m for various values of m, we obtain

S₆ = 0.3152   S₈ = 0.3060   S₉ = 0.3056   S₁₀ = 0.3064

Other values of m give larger values of S_m. Thus the optimum value for S_m is taken as 0.3056, which is assumed for m = 9. The corresponding values of r_k are

The average factor of reduction per double sweep thus obtained, namely, 0.3056, compares favorably with 0.73 for the successive-over-relaxation method even though each iteration of the latter method involves only a single sweep.

As mentioned above, the number of iterations required to achieve a specified accuracy by using the Peaceman-Rachford method is approximately proportional to |log h|⁻¹ as compared with h⁻¹ for the successive-overrelaxation method. Before one is tempted to discard the latter method entirely, however, there are several factors to be considered. In the first place, the result stated above has been proved only for the five-point difference-equation representation of Laplace’s equation and for a rectangle. As a matter of fact, Birkhoff and Varga⁵ showed that in the case of Laplace’s equation the methods used in the theoretical analysis of the Peaceman-Rachford method are not applicable for regions other than rectangles. On the other hand, numerical experiences to date^52,60 indicate that the method works well for considerably more general problems. In the second place, for problems involving a moderately small mesh size the actual machine time required for using the successive-overrelaxation method may be less, either because the total number of sweeps is less or because the time required per sweep is less, even though the total number of sweeps might be slightly greater. Finally, the successive-overrelaxation method is much simpler both because of the arithmetic operations involved and because it is not necessary with the successive-overrelaxation method to obtain the approximate values first by rows and then by columns. This may be a considerable advantage if magnetic-tape storage must be used.

The theoretical analysis underlying the Peaceman-Rachford method has not yet been extended to include differential equations that are essentially different from Laplace’s. Nevertheless, the method can be formally applied. One approach would be first to estimate a and b and then determine the r_k as described above. After some experience with several cases involving a certain differential equation, improved estimates can probably be obtained.

15.5    Other Iterative Methods for Solving Elliptic Equations

It is probable that at the present time one would choose either the successive-overrelaxation method or the Peaceman-Rachford method when solving an elliptic partial differential equation by finite-difference methods; nevertheless, it is well to be aware of the existence of other methods. One of these other methods may actually be best suited for a given problem. Moreover, it may well happen that such a method, or a slight modification thereof, may turn out to be generally superior to either of the two methods that are now more frequently used.

A method of L. F. Richardson,³⁶ which dates back to 1910, can probably best be described in terms of its application to a linear system

where u is an unknown vector, or column matrix, d is a given column matrix, and A is a given square matrix that is positive definite and symmetric. The difference-equation representation described in Sec. 15.2 of a self-adjoint elliptic equation will have a positive-definite symmetric matrix if every interior mesh point is regular. Starting with an arbitrary initial approximation u⁽⁰⁾ to the exact solution, one computes u⁽¹⁾, u⁽²⁾, … according to the formula

Here the β_n are constants, and their selection determines the rapidity of convergence. Normally, one chooses a set of m different values of β and employs them in a cyclic order, β₁, β₂, …, β_m, β₁, β₂, …. If we let e⁽ⁿ⁾ denote the error after n iterations as defined by

e⁽ⁿ⁾ = u⁽ⁿ⁾ − u

then, since Au + d = 0, we have, by Eq. (15.39),

As described by Shortley⁴¹ and by Young,⁵⁸ one chooses the coefficients β_n so as to minimize the maximum absolute value of the polynomial

over the range , where a and b are positive lower and upper bounds, respectively, for the eigenvalues of A. The polynomial with the smallest maximum absolute value is found by using methods involving Chebyshev polynomials. It can be shown that, when Richardson’s method is applied to Laplace’s equation, the number of iterations varies as h⁻¹, as the mesh size h tends to zero. The number of iterations is thus of the same order of magnitude as for the successive-overrelaxation method; however, since the latter method converges at least twice as fast as Richardson’s method and since it is simpler and requires less storage, its use is recommended in preference to Richardson’s method. Moreover, it has been found by Young and Warlick⁶³ that with Richard-son’s method there is serious difficulty in controlling rounding errors even when great care is taken in selecting the order in which the β_n are to be used. This difficulty may be at least partially overcome by the use, in place of Eq. (15.39), of a formula involving u⁽ⁿ⁻¹⁾ as well as u⁽ⁿ⁾ and u⁽ⁿ⁺¹⁾, as described in Ref. 46. This formula is based on the use of a three-term recurrence relationship for Chebyshev polynomials. Even with this improvement, however, Richardson’s method appears inferior to the successive-overrelaxation method for solving elliptic difference equations, particularly since the convergence is slower and since the method is more complicated and requires considerably more storage. [Quite recently, though, Golub and Varga (in “Chebyshev Semi-iterative Methods, Successive Overrelaxation Iterative Methods, and Second Order Richardson Iterative Methods,” Report 1028, Case Institute of Technology, Cleveland, Ohio, February, 1960) have modified the method so that it converges more rapidly than the successive-overrelaxation method. The modified method closely resembles the successive-overrelaxation method.]

Sheldon⁴⁰ developed a method that combined the ideas used in Aitken’s “to-and-fro” method,¹ the successive-overrelaxation method, and Richardson’s method. The to-and-fro method as applied to the successive-overrelaxation procedure involves improving the approximate values according to a given ordering on the first half of an iteration and then proceeding in exactly the reverse order for the second half of the iteration. This procedure has the feature of making the eigenvalues of the associated linear transformation real, whereas in the ordinary successive-overrelaxation method the eigenvalues are complex. By having real eigenvalues, we can then combine successive iterants according to a procedure based on the use of Chebyshev polynomials, as in Richardson’s method, in order to increase the rate of convergence still further. Unfortunately, however, since no definite results are available concerning the spectral radius of the successive-overrelaxation method under the to-and-fro ordering, it is difficult to estimate the effectiveness of the method. Until this is possible with some degree of confidence, it appears doubtful that Sheldon’s method will receive general usage, particularly in view of the complications and extra storage requirements involved.

Arms and Gates,³ with Zondek, and Friedman¹⁸ considered the use of “line relaxation,” by which at each step, instead of modifying the approximate value of u at a single point, one modifies simultaneously the values on an entire line of points. The procedure for modifying the values in a line of points is the same as that described in the preceding section. With line relaxation, however, one treats the lines—which are usually either rows or columns—successively, new values being used as soon as available. By the introduction of overrelaxation, it is possible to obtain for the Dirichlet problem a rate of convergence that is faster by a factor of approximately than the (point) successive-overrelaxation method. Moreover, as shown recently by Cuthill and Varga,¹¹ in many cases it is possible to perform the line relaxation in the same number of arithmetic operations per iteration as required by the successive-overrelaxation method. While the line-relaxation method is somewhat more complicated, nevertheless, if a number of cases are to be solved for the same problem, then line relaxation should be considered as a means of reducing machine time by nearly 30 per cent.

The method of Douglas and Rachford¹⁴ is rather similar to the Peaceman-Rachford method described in the previous section. It has the advantage of applying to three-dimensional problems as well as to two-dimensional problems. On the other hand, for two-dimensional problems in which a direct comparison is possible, the Peaceman-Rachford method converges considerably more rapidly. It is thus not surprising that the Douglas-Rachford method has so far not been used as much. This situation may change, however, with the development of larger and faster computers capable of solving three-dimensional problems.

The discussions of the various iterative methods considered in this section have necessarily been brief. The objective has been to provide sufficient information concerning each method so that the reader can decide whether the method might be useful to his particular problem. More of the details involved in using and in estimating the effectiveness of the methods can be found in the appropriate references.

15.6    Parabolic Equations—Forward-difference Method

The parabolic partial differential equations that we shall consider in this section and in the next two sections are of the form

where L(U) is a second-order elliptic partial differential expression with either one or two “space” variables. With one space variable, L(U) would have the form

whereas with two space variables it would have the form

Here A and C are functions of the space variables such that A and C are positive in the domain of the space variables under consideration. The functions ϕ and ψ are assumed to be analytic in their respective arguments.

A typical problem involving Eq. (15.40) would be the following: Let R be a region in the domain of the space variables and let S be the boundary of R. We wish to find a function U(x,t) that satisfies Eq. (15.40) for t > t₀ and for x in R such that U equals a prescribed function g on S for t > t₀. Here x = x in the case of one space variable and x = (x,y) in the case of two space variables.

The problem just formulated is known as an initial-value problem. Since, as in the case of boundary-value problems for elliptic equations, analytic methods are seldom effective as a means for obtaining numerical solutions, approximate methods are usually employed. As before, the method of finite differences is probably the best such method. In order to apply the method of finite differences, one superimposes a network over the region R + S in the same manner as described in Sec. 15.2. In addition, a time-mesh size k is selected. For each discrete value of t, one desires to determine approximate values of U for all mesh points in R.

Fig. 15.4 Forward-difference method: typical neighborhood configuration.

The forward-difference method (see Fig. 15.4) for determining approximate values of U involves the replacement of the time derivative in Eq. (15.40) by the forward-difference quotient

The differential expression L(U), on the right-hand side of Eq. (15.40), is replaced by a difference expression L*(U), which is derived by means of the methods of Sec. 15.2. By letting u(x,t) denote the solution of the difference equation and solving for u(x, t + k), we obtain

In the case of the differential equation

for example, the difference equation would be

where h is the mesh size for the variable x and where

Whereas for elliptic equations the value of the solution at any one point depends on all the boundary values, in the case of the parabolic equations considered here the value of U(x,t₁) depends only on the initial values and on the values of g for t₀ ≤ t ≤ t₁. This means that it should be possible to solve the difference equation by a kind of “marching” process rather than being required, as in the case of elliptic equations, to determine the solution at all points in order to be able to evaluate it at any one point. In the case of the forward-difference method, the marching process is particularly easy to carry out. By knowing U(x,t₀), one can determine L*(U) for t₀ and thus compute U(x, t₀ + k). In the same way one can compute U(x, t₀ + 2k), U(x, t₀ + 3k), etc. Unfortunately, as we shall presently see by considering a specific example, the usefulness of this very simple procedure is very limited because of the problem of stability. Without attempting to formulate a precise mathematical definition of the term, we shall simply state that a method that is not stable is unsuitable for numerical computation because of the rapid growth of rounding errors that may, after a few steps, cause all accuracy to be lost.

In order to obtain stability with the forward-difference method, one is frequently forced to use a value of k that is extremely small relative to the space-mesh size. In such a case, the number of time steps necessary to reach a desired value t₁ starting from t = t₀ may be prohibitively large.

To illustrate, let us consider a special case of Eq. (15.44), namely, the diffusion equation

in the region 0 < x < 1, t > 0, subject to the boundary conditions

The interval 0 ≤ x ≤ 1 is divided into subintervals of length h = M⁻¹, where M is an integer. The difference equation (15.45) becomes

u(x, t + k) = u(x,t) + r[u(x + h, t) + u(x − h, t) − 2u(x,t)]

By letting u_i,j = u(ih,jk), f_i = f(ih), g_1,j = g₁(jk), etc., we have

u_i,j+1 = r(u_i+1,j + u_i−1,j) + (1 − 2r)u_i,j

for i = 1, 2, …, M − 1; j = 0, 1, 2, …. The boundary conditions are

Now, letting f(x) ≡ 1, g₁(t) ≡ g₂(t) ≡ 0, and M = 4, for the case r = ½ we get

The numerical solution for j = 0, 1, 2, …, 5 is given in Table 15.3(I).

Table 15.3 Examples of the Use of the Forward-difference Method

Ignoring for the moment the quantities in parentheses, one sees that even though the values obtained are not smooth, they at least appear to be reasonable. It is not difficult to believe that, as h and k approach zero, with the ratio r = k/h² being held fixed, the numbers obtained by this process would approach the exact solution of the differential equation. To investigate the stability of the process, we introduce an error of amount in u_2,2. Because of linearity, we can study the effect of the error separately as indicated in the table. We note that the numbers in the parentheses decrease in magnitude as j increases. This is characteristic of a stable process.

Let us now consider the case where r = 1. The difference equation becomes

u_i,j+1 = u_i+1,j + u_i−1,j − u_i,j

The results obtained by using this formula are given in Table 15.3(II). Since it is easy to show that the exact solution of the differential equation is always bounded between zero and unity, it is clear that the computed values are very much in error and that the error increases rapidly with j. We are not surprised to learn that the method is not convergent. Moreover, a stability analysis similar to that performed previously for the case r = ½ indicates that the process is highly unstable.

If the initial values were defined by f(x) ≡ sin πx, then for r = 1 the forward-difference method would be convergent though not stable. One would obtain the values indicated in Table 15.3(III). Here the values are at least respectable until one reaches j = 5, when, because of rounding errors, they begin to oscillate. If more decimal places had been retained in the calculations, then the occurrence of the oscillations would have been delayed. Finally, if exact values had been used throughout—i.e., if, for instance, had been used in place of 0.707—then the oscillations would never have occurred. Of course, this latter consideration is largely academic except that it does illustrate the distinction between the concepts of stability and convergence.

Detailed discussions of convergence and stability for the forward-difference method and other methods of solving parabolic equations have been given by John,²³ Douglas,¹³ O’Brien, Hyman, and Kaplan,³⁴ Lax and Richtmyer,³⁰ Leutert,³² Todd,⁴⁸ Juncosa and Young,²⁴ and others. The forward-difference method as applied to the problem specified by Eqs. (15.47) and (15.48), with g₁(t) ≡ g₂(t) ≡ 0 and with f(x) piecewise continuous, is convergent and stable provided that

For the corresponding differential equation with two space variables,

U_t = U_xx + U_yy

the condition on r for stability is

where the mesh size in each space-coordinate direction is taken equal to h. In both cases, the use of a larger value of r would make the process unstable. Consequently, a time-mesh size that is of the order of the square of the space-mesh size must be used. Thus the total computational effort involved in proceeding from t₀ to t₁, with space-mesh size h and with time-mesh size k = rh², is proportional to h⁻³.

Thus, while the forward-difference method is very simple and straightforward for use on a high-speed computer, stability considerations usually require the use of such a small time increment k that the method is generally not practical. In the next two sections, we shall consider methods that, though less simple and though requiring considerably more arithmetic operations per time step, nevertheless more than compensate for this by making possible the use of a much larger time-step size.

15.7    The Crank-Nicolson Method

The difference equation used in the Crank-Nicolson method (see Fig. 15.5) for solving Eq. (15.40) is the same as the forward-difference method except that the differential expression L(U) is replaced by

½[L*(U) + L**(U)]

instead of by L*(U). Here, as before, L*(U) is the difference-expression representation of L(U) for the current value of t. The difference-expression representation of L(U) for t + k is represented by L**(U). The difference equation thus obtained can be written in the form

Fig. 15.5 Crank-Nicolson method: typical neighborhood configuration.

In the case of the differential equation (15.44), this becomes

Evidently both Eq. (15.49) and Eq. (15.50) are implicit equations in the sense that values of u for t + k appear on both sides of the equations. Moreover, since in Eq. (15.50) the quantities u(x, t + k), u(x − h, t + k), and u(x + h, t + k) appear, it is not possible to solve for any one of them explicitly. In order to proceed from one time step to another, it is thus necessary to solve a system of equations. If the function ϕ appearing in Eq. (15.50) is nonlinear in u, then the equations are nonlinear. For the problem defined by Eqs. (15.47) and (15.48), however, by using the notation u_i,j = u(ih,jk), etc., we get

One method of determining, for each j, the quantities u_1,j+1, u_2,j+1, u_3,j+1, … would be to use the successive-overrelaxation method. It is not difficult to show that the spectral radius of the Gauss-Seidel method for the system satisfies the inequality

from which the optimum value of ω can be determined by

This procedure was used by Young and Ehrlich.⁶¹ Since each iteration of the successive-overrelaxation method on the average multiplies the magnitude of the error by ω_b − 1, it is clear that, as r increases, the rapidity of the convergence decreases. If k is chosen proportional to h, so that r increases as h⁻¹, then the number of iterations required varies as h^−½. Since the number of mesh points involved in proceeding from t₀ to t₁, with space-mesh size h and with time-mesh size proportional to h, is proportional to h⁻², it follows that the total computational effort involved is proportional to , as compared with h⁻³ for the forward-difference method.

Another method for solving Eq. (15.51) is to use the special procedure described in Sec. 15.4 for solving a linear system with a tridiagonal matrix. The specific formulas are

As an example, let us consider the numerical example given in the preceding section with r = 1. Table 15.4 shows the calculations involved in proceeding from j = 0 to j = 1. One can easily verify that the values obtained satisfy the Crank-Nicolson difference equation (15.51). Furthermore, as mentioned earlier, the process is stable and is otherwise satisfactory for use on a high-speed computer. Since the number of calculations per point is clearly independent of the number of points involved, it follows that the total computational effort involved in proceeding from t₀ to t₁ is proportional to h⁻², as compared with with the successive-overrelaxation method and as compared with h⁻³ with the forward-difference method.

Table 15.4 Example of the Computations Involved in Using the Crank-Nicolson Method for f(x) ≡ 1, r = 1

In the above analysis we have assumed that it would be possible to let k be proportional to h, rather than to h² as with the forward-difference method. For the problem specified by Eqs. (15.47) and (15.42), it is not difficult to show that the Crank-Nicolson method is stable for any positive h and k; see for instance Ref. 34. For the case where g₁(t) ≡ g₂(t) ≡ 0 and where f(x) is piecewise continuous, Juncosa and Young²⁵ have shown that, as h approaches zero, the solutions obtained by the Crank-Nicolson method converge to the true solution of the differential equation provided that k = O(h/|log h|). The latter condition implies that k approaches zero slightly faster than h. Actually, it is probable that the result could be proved for k = O(h); some numerical experiments of Young and Ehrlich⁶¹ suggest that this is the case. These experiments, together with theoretical results of Juncosa and Young,^24,25 suggest that the order of accuracy of the solution of the Crank-Nicolson difference equation is about the same as the order of the accuracy of the solution of the forward-difference method.

Let us now consider the differential equation (15.44), for which the Crank-Nicolson difference equation is given by Eq. (15.50). If ϕ is nonlinear in u, then one is faced with the problem of solving a system of nonlinear equations. One possibility would be to use the successive-overrelaxation method, wherein on any iteration the value of u(x, t + k) from the previous iteration would be used to evaluate ϕ[x, t + k, u(x, t + k)]. For the cases ϕ ≡ e^u and ϕ ≡ e^−1/u, the experience of Young and Ehrlich⁶¹ indicated no significant loss in convergence rate because of the nonlinear term. Probably a better procedure would be to iterate a fixed number of times by using the special procedure for solving a linear system with a tridiagonal matrix. The term ϕ[x, t + k, u(x, t + k)] would be evaluated from the previous iteration and treated as a known quantity. If, as seems reasonable, the use of a fixed number of iterations would suffice for all mesh sizes, then the total computational labor would still be proportional to h⁻² as for the linear case.

In applying the Crank-Nicolson method to problems involving two space variables, one would have to solve, at each time step, an elliptic boundary-value problem. For the case of the differential equation

U_t = U_xx + U_yy

however, we can get a fairly accurate estimate of the optimum relaxation factor to be used with the successive-overrelaxation method. It can be shown that the spectral norm λ of the Gauss-Seidel method satisfies the inequality

Consequently, if one chooses k proportional to h, then the number of iterations would be proportional to h^−½, as compared with h⁻¹, as would be required if one were solving the equation

U_xx + U_yy = 0

In spite of this saving, it would certainly seem desirable to have a method that would permit one to proceed from t to t + k without the necessity of solving a system of equations and without the limitations of the forward-difference method.

15.8    The Alternating-direction Method for Parabolic Equations Involving Two Space Variables

Peaceman and Rachford³⁵ and Douglas¹² presented a method for solving the parabolic equation

The process of passing from one time step to the next is done in two parts, one implicit on rows and the other implicit on columns. The specific formulas are

where r = k/2h².

The formula (15.41) is very similar to those involved in performing an iteration by using the Peaceman-Rachford method for elliptic equations. As a matter of fact, the alternating-direction method for parabolic equations and the Peaceman-Rachford method for elliptic equations are very closely related, the latter method being suggested by the former.^12,35 The computations for both methods involve at each stage the solution of linear systems with tridiagonal matrices. After writing each system in the form (15.32), one uses Eqs. (15.33), first to obtain u(x, y, t + k/2) for all (x,y) and then to obtain u(x, y, t + k).

It is shown in Ref. 12 that the local approximation of the difference equations (15.53) to the differential equation (15.52) is of the same order in h and k as that of the Crank-Nicolson difference equation. It is further shown that, for the case of the rectangle, the method is stable for any positive h and k. On the other hand, either part of Eqs. (15.53) when used alone would not, in general, be stable. Although stability has so far been proved only for the case of the differential equation (15.52) for the rectangle, it may well hold for more general differential equations and for more general regions. Certainly, the method can be formally applied in more general situations.

It is clear that the method, where applicable, affords a substantial improvement over the Crank-Nicolson method. The amount of computation required per point is independent of h and k, while for the Crank-Nicolson method, using the successive-overrelaxation method, the effort per point is proportional to h^−½ as stated in the previous section. Whenever the alternating-direction method is stable and whenever the time-mesh size k can be taken proportional to h without loss of accuracy, the method would appear to be superior. It may be very difficult, however, to estimate stability and accuracy in advance.

The method of Douglas and Rachford¹⁴ is similar to the alternating-direction method; moreover, it can be extended to apply to three-dimensional problems. As previously stated in Sec. 15.5, the Douglas-Rachford method also leads to an effective iterative method for solving elliptic equations.

15.9    Illustrative Examples

Prior to the advent of modern high-speed computing machines, manual methods were used for solving elliptic and parabolic partial differential equations. Southwell⁴⁵ and his colleagues solved a number of problems by using relaxation methods, wherein the human computer was free to use ingenuity and intuition in order to speed up the convergence. Problems were solved in the following areas: torsion and shear stress for shafts with various cross sections; electrostatics; magnetism; conformal mapping; steady-state heat conduction; film lubrication; gas flow; and hydrodynamics, including flows where the boundaries are not initially known. Allen² also described the solution of a number of problems by relaxation methods. He considered membrane-vibration problems as well as plate-bending, compression, and stretching problems involving the biharmonic equation. Shortley, Weller, Darby, and Gamble⁴³ described the solution of problems in electrostatics and torsion. Instead of relaxation methods, systematic iteration methods were used, including the Gauss-Seidel method, block iteration, and an extrapolation procedure.

Many problems that involve the three space variables x, y, and z were treated as two-dimensional problems by assuming that the functions involved are independent of z. The following examples are given of the mathematical formulations of typical problems:

a. The torsion problem for a shaft of uniform cross section involves the solution of the Poisson equation

ψ_xx + ψ_yy = −2

in a plane region subject to the requirement that ψ = 0 on the boundary of the region.

b. The steady-state temperature distribution in a homogeneous material of uniform cross section involves the solution of Laplace’s equation

T_xx + T_yy = 0

in a region R. Frequently, the temperature distribution on the boundary is prescribed.

c. The problem of determining the current in a plane conducting sheet with resistivity ρ involves the solution of the equation

(ρψ_x)_x + (ρψ_y)_y = 0

where the current function ψ assumes prescribed values on the boundary. The currents in the x and y directions are ψ_y and − ψ_x, respectively.

With the availability of high-speed computers, problems of the type just described could be handled very easily, once the necessary computer program had been prepared. A program for eliminating or reducing the effort required to prepare such programs is discussed in the next section. Our discussion for the remainder of this section will be devoted to describing work on three problems for which approximate numerical solutions have been obtained by use of large computers.

a. Cavity-flow Problem.   An attempt to determine the axially symmetric flow past a disk was described by Young, Gates, Arms, and Eliezer.⁶² In order to obtain a cavity of finite length behind the plate, a second artificial plate was introduced downstream, as shown in Fig. 15.6.

Fig. 15.6 Cavity-flow problem.

This is the so-called Riabouchinsky model.³⁷ In the interior of the flow, the Stokes’ stream function ψ satisfies the equation

where x is the distance along the axis of symmetry and r is the distance normal to the axis of symmetry. The stream function vanishes on the axis of symmetry, on the plates, and on the cavity boundary, the last of these being unknown. In addition, the stream function satisfies the conditions

and, on the cavity boundary,

where ψ_n denotes the inward normal derivative of ψ and where λ is a constant to be determined. The last condition, which states that the velocity along the cavity boundary must be constant, is derived from the Bernoulli equation and from the requirement that the pressure in the cavity be constant.

In order to apply the method of finite differences, artificial boundaries DC and DE for the flow were introduced as shown in Fig. 15.7. The limiting values of ψ, as x² + r² → ∞, were assumed on these artificial boundaries. By symmetry, only one-quarter of the flow was considered, the symmetry condition ψ_x = 0 being imposed on BC. A network was superimposed over the region. The spacing of the mesh points was made to vary so that the spacing became finer as one approached the separation point A. The basic procedure was as follows: The curve AB and initial values for ψ at all interior mesh points were assumed. By means of the successive-overrelaxation method, values of ψ were determined at all interior mesh points. For this determination, only the condition ψ = 0 was imposed on AB. After the values of ψ were determined, the next step was to compute approximate values for r⁻¹ψ_n. If these values were constant along AB, then the problem would be considered solved, with λ being selected as the constant value. Otherwise, a new curve AB would be assumed and the process repeated.

Fig. 15.7 Mathematical formulation of cavity-flow problem.

The attempt to solve the problem on the NORC computer at the U.S. Naval Proving Ground, Dahlgren, was only partially successful, largely because of difficulties in the neighborhood of the separation point A. Although a fine mesh was used close to this point, it was necessary to use an analytic approximation in the immediate neighborhood of A based on the known two-dimensional solution. It was not possible to get the computed velocities constant to an accuracy of better than ±1.5 per cent of the average value. The fact that no systematic method was known for modifying the boundary added to the difficulties. Estimates are needed for the errors introduced by the imposition of the artificial boundaries, by the use of finite-difference methods, and by the use of approximate analytic methods at the stagnation point.

As for the numerical techniques employed, it should be mentioned that, once the curve AB was specified, the mesh points and the distances to the boundary were automatically determined by the machine. There were approximately 1,150 interior mesh points. Unfortunately, the successive-overrelaxation method, while considerably faster than the Gauss-Seidel method, was much less effective than expected. In fact, since the use of the estimated best value of ω led to a divergent process, it was necessary to use a much smaller value of ω with the accompanying loss in convergence speed. It is believed that the difficulty was caused by the choice of the difference-equation representations of the differential equation for points connecting regions of a fine mesh and of a coarser mesh. When the value 1.4, which appeared to be optimal, was used for ω, approximately 100 iterations were required to reduce the residuals to 10⁻⁵, starting with reasonable initial estimates. Since the NORC performed one iteration in approximately 5 sec, it was possible to treat one trial cavity boundary in less than 10 min.

b. Nuclear-reactor Problems.   An important problem arising in the design of nuclear reactors is the determination of the neutron flux inside the reactor. Figure 15.8 shows a cross section of a very highly idealized reactor. The regions R₁ and R₂ are diffusion regions. The C_i, i = 1, 2, 3, 4, represent control-rod regions. For simplicity, we shall consider the two-group case discussed by Varga.⁵⁰ Here it is assumed that each neutron belongs either to the fast group or to the slow group. The functions ϕ(x,y) and ψ(x,y) represent, respectively, the slow-neutron flux and the fast-neutron flux. The functions ϕ and ψ satisfy the equations

Fig. 15.8 Nuclear-reactor problem.

The coefficients D_s, D_f, σ_s, σ_f, S_s^f and S₂₅ are assumed to be regionwise constant. The constant η is to be determined. The last term in the first equation represents the degeneration of fast neutrons to slow neutrons; the last term in the second equation represents the creation of fast neutrons by the fission of slow neutrons.

On the outer boundary, both ϕ and ψ vanish. Across interfaces between diffusion regions, the functions ϕ, ψ, D_s∇ϕ and D_f∇ψ are required to be continuous. On interfaces between diffusion regions and control-rod regions, both ψ and D_f∇ψ must be continuous. Although ϕ is not defined in a control-rod region, the condition

must be satisfied on the boundary of each such region.

In order to apply the method of finite differences, Varga⁵⁰ superimposed a network over the entire region, requiring that all interfaces lie on mesh lines. The difference equations are derived by using discrete representations of certain integrals obtained from the differential equations by means of Green’s theorem. Such a procedure allows for the incorporation of the interface conditions in a natural way. The values of ϕ and ψ at each mesh point are expressed as linear combinations of the values at the four neighboring mesh points.

Having determined two difference equations analogous to Eqs. (15.54) and (15.55), one carries out the following procedure. Starting with arbitrary initial approximations ϕ⁽⁰⁾, ψ⁽⁰⁾, and η₀ for ϕ, ψ, and η, respectively, one then solves the first equation for ϕ⁽¹⁾. Here the first equation is treated as an elliptic boundary-value problem; the successive-overrelaxation method of iteration is used. Next, by means of the value of ϕ⁽¹⁾ just obtained, the second equation is solved for ψ⁽¹⁾. A new value of η is determined from

where, in general, the inner product (f, g) of two functions denotes the sum of the products fg taken over all mesh points. The entire process is repeated starting with ϕ⁽¹⁾, and ψ⁽¹⁾, and η₁, is continued until convergence is achieved.

The entire process of determining η_i+1 from η_i is called an outer iteration. Each outer iteration consists of two sets of inner iterations, the inner iterations being merely iterations of the successive-overrelaxation method.

Two similar codes have been prepared at the Westinghouse Bettis Plant for the IBM 704 computer. The first, which is known as the QED code and is described in Ref. 50, is restricted to two groups. The second, which is known as the PDQ code and is described in Ref. 4, is a substantially improved version of the QED code and handles problems involving up to four groups. In both programs, the coefficients of the difference equation are computed automatically. Next, upper and lower bounds for the optimum values of ω to be used for the n sets of inner iterations are computed, n being the number of groups. With the value of ω chosen for each group, the iteration process begins. The PDQ program is capable of handling from 1,250 to 6,500 mesh points, depending on the size of the computer-core storage. It is reported in Ref. 4 that a typical two-group, 2,500-point problem can be solved in less than one hour.

Wachspress⁵¹ described the CURE program for solving similar problems. The CURE program, which has also been prepared for the IBM 704 computer, can handle problems involving up to three groups. As many as 644 mesh points can be taken. The Peaceman-Rachford method of iteration is used. Because of the relatively few mesh points involved, however, it is by no means clear that the use of the Peaceman-Rachford method is preferable to the use of the successive-overrelaxation method. On the other hand, even though the basic theory for the Peaceman-Rachford method has not been shown to hold except for Laplace’s equation and a rectangle, experiences in which the CURE program was used indicate that the method converges as rapidly as predicted, even for much more general cases.

c. Thermal-ignition Problem.   Hicks, Kelso, and Davis²² described the mathematical formulation of a thermal-ignition problem. An explosive material is placed in contact with a hot gas. Soon the temperature begins to rise very rapidly. The objective is to determine the time required for the maximum temperature to reach a certain value, known as the ignition temperature. A highly simplified version of the problem involves the solution of the differential equation

U_t = U_xx + e^−1/U0 < x t > 0

in the region 0 < x, subject to the initial condition

U(x, 0) = U₀0 < x

and to the boundary conditions

Here units have been chosen to make the coefficients of U_xx and of e^−1/U in the differential equation unity. Thus U represents a “reduced temperature.” The constants a, b, and U₀ represent, respectively, the reduced heat-transfer coefficient, the reduced gas temperature, and the reduced initial temperature.

In order that the problem could be attacked numerically, the second boundary condition was replaced by

U_x(A,t) = 0t > 0

where A is an appropriately chosen number. Young and Ehrlich⁶¹ computed an approximate solution on the ORDVAC computer at the Aberdeen Proving Ground, Maryland, by using the Crank-Nicolson method as described in Sec. 15.7. In proceeding from t to t + k, however, the function e^−1/U was always evaluated at t. In this way, it was sufficient to solve a linear system with a tridiagonal matrix. Greater accuracy would undoubtedly have been obtained if the average value of e^−1/U for t and t + k had been used. This improvement could probably have been accomplished with one or two iterations of the procedure that was actually used.

Other examples of elliptic and parabolic partial differential equations include a weather-prediction problem described by Charney and Phillips⁷ and a thermal-melting problem described by Ehrlich.¹⁵ In the latter problem, the material was assumed to melt upon reaching a given temperature. The molten material then remained in the system but with different thermal properties. One of the objectives of the computation was to determine the rate at which the melting occurred. In proceeding from one value of t to the next, it was necessary to determine the amount of melting as well as the changes in temperature.

With the availability of machines having very large storage capacities, it is expected that the solution of problems involving three space variables will be possible. This should greatly improve the accuracy obtainable and hence the value of the computations for scientific and engineering purposes.

15.10    The SPADE Project for the Development of a Computer Program for Solving Elliptic and Parabolic Equations

One of the principal obstacles encountered when one attempts to solve a partial differential equation by using a high-speed computer is the time and effort required to prepare an appropriate computer program. Each of the programs described in the preceding section required several months to prepare. It would thus have been very useful to have had available a comprehensive program, or a set of programs, capable of handling, with a minimum amount of input information, any one of a large class of problems.

The preparation of such a “library” of computer programs for the IBM 704 and 709 computers was the objective of the SPADE Committee.† The programs of the library were to be capable of performing the arithmetic, logical, and bookkeeping operations involved in the solution of the class of partial differential equations by finite-difference methods. While the user still would be fully responsible for the numerical analysis, including choice of mesh size and difference equation, as well as the stability, convergence, and accuracy, the programs of SPADE would relieve him of much of the programming drudgery. The user would specify the region, mesh size, boundary conditions, differential equations, iterative method, etc. He also would prepare a small coordinating program to connect those routines of the library that are involved. Alternatively, there would be complete program “packages” available for handling frequently used classes of problems.

The problems included within the scope of the SPADE programs would involve systems of partial differential equations. There would be two independent space variables and, possibly, one independent “time-like” variable. Each equation would have the form

where the U⁽¹⁾, U⁽²⁾, … are the dependent variables, t is the timelike variable, and the L_i are differential expressions independent of t, such that the system

L_i(U⁽¹⁾, U⁽²⁾, …, U^(m)) = 0i = 1, 2, …, m

is elliptic. The coefficients a_i may be either zero or one. It is assumed that there exist suitable difference representations for the L_i of the form

Here U_i^(k) designates the value of U^(k) at the mesh point P_i, and p is the number of neighboring points of the point P₀ that is under consideration. The coefficients α_{k, j}, and β^(k) may depend on U_j⁽¹⁾, U_j⁽²⁾, …, U_j^(m), j = 1, 2, …, p, as well as on P₀. It may be well to emphasize here that the routines of SPADE are designed to handle the operations that are needed to solve problems of the kind now discussed, on the assumption that the usual methods apply. The responsibility for determining whether or not this is the case, however, would rest entirely with the user of the library.

The boundary conditions are of the form

where the coefficients α_k^(v), β_k^(v), and C^(v) may depend on U⁽¹⁾, U⁽²⁾, …, U^(m), as well as on the independent variables. In addition to the boundary of the region, there may be one or more interfaces in the region, each with appropriate continuity conditions for the dependent variables and for certain of the derivatives normal to such interfaces.

The SPADE programs have been divided into three classes, or “phases.” In Phase I, the coordinates and other appropriate items of information are computed for each mesh point. The input to Phase I consists of specifications of each boundary and interface as a series of paths. A path may be specified either as the arc of a conic section the equation of which is prescribed as input data, or as a series of points. In the latter case, the machine treats consecutive points as though they were joined by line segments. In addition, the mesh sizes are specified. While each mesh line is required to be parallel to a coordinate axis, the spacing need not be uniform. Furthermore, subdivisions of the mesh by factors of powers of two in certain areas are allowed. At the completion of Phase I, the information for each boundary and interior mesh point is stored on magnetic tape as a point record.

In Phase II, each point record is augmented to include the coefficients involved in the difference equation at the point. The routines to be prepared first require that algebraic formulas for the coefficients be provided in terms of the computed distances from the point to its neighbors and in terms of the coefficients of the differential equations. The latter would be available as subroutines; hence, the numerical coefficients could be computed automatically. More sophisticated routines might actually allow the machine to derive the difference equation automatically.

In Phase III, the difference equations would be solved by an iterative process. In the case of a time-dependent problem, the result would be merely the completion of a single time step followed by a return to Phase II for a recalculation of the coefficients. The same would be true in the nonlinear case in which the coefficients, as functions of the dependent variables, would be recomputed every few iterations. For the time-independent linear case, however, it is clear that, since the system of linear algebraic equations is completely specified by the information on magnetic tape, practically any iterative method could in principle be used, including the successive-overrelaxation method, the Peaceman-Rachford method, or other methods mentioned in Sec. 15.5.

REFERENCES

1.   Aitken, A. C., On the Iterative Solution of a System of Linear Equations, Proc. Roy. Soc. Edinburgh, Sect. A, vol. 63, pp. 52–60, 1950.

2.   Allen, D. N. de G., “Relaxation Methods,” McGraw-Hill Book Company, Inc., New York, 1954.

3.   Arms, R. J., and L. D. Gates, Jr., “The Computation of an Axially Symmetric Free Boundary Problem on NORC, Part II,” U.S. Naval Proving Ground Report No. 1533, Dahlgren, Va., Apr. 5, 1957.

4.   Bilodeau, G. C., W. R. Cadwell, J. P. Dorsey, J. G. Fairey, R. S. Varga, “PDQ—An IBM-704 Code to Solve the Two-dimensional Few Group Neutron Diffusion Equations,” Report WAPD-TM-70, Bettis Plant, Westinghouse Electric Corporation, Pittsburgh, Pa., August, 1957. Copies available from the Office of Technical Services, U.S. Department of Commerce, Washington, D.C.

5.   Birkhoff, Garrett, and Richard S. Varga, Implicit Alternating Direction Methods, Trans. Amer. Math. Soc., vol. 92, pp. 13–24, 1959.

6.   Bruce, G. H., D. W. Peaceman, H. H. Rachford, Jr., and J. D. Rice, Calculation of Unsteady-state Gas Flow through Porous Media, Petroleum Trans., A.I.M.E., vol. 198, pp. 79–92, 1953.

7.   Charney, J. G., and N. A. Phillips, Numerical Integration of the Quasi-geostrophic Equations for Barotropic and Simple Baroclinic Flows, J. Meteorol., vol. 10, pp. 71–99, 1953.

8.   Courant, Richard, Hyperbolic Partial Differential Equations and Applications, chap. 5 in “Modern Mathematics for the Engineer,” First Series, edited by E. F. Beckenbach, McGraw-Hill Book Company, Inc., New York, 1956.

9.   ———, K. Friedrichs, and H. Lewy, Über die partiellen Differenzengleichungen der mathematischen Physik, Math. Ann., vol. 100, pp. 32–74, 1928.

10.   Crank, J., and P. Nicolson, A Practical Method for Numerical Evaluation of Solutions of Partial Differential Equations of the Heat-conduction Type, Proc. Cambridge Philos. Soc., vol. 43, pp. 50–67, 1947.

11.   Cuthill, Elizabeth H., and Richard S. Varga, A Method of Normalized Block Iteration, J. Assoc. Comput. Mach., vol. 6, pp. 236–244, 1959.

12.   Douglas, Jim, Jr., On the Numerical Integration of ∂²u/∂x² + ∂²u/∂y² = ∂u/∂t by Implicit Methods, J. Soc. Indust. Appl. Math., vol. 3, pp. 42–65, 1955.

13.   ———, On the Relation between Stability and Convergence in the Numerical Solution of Linear Parabolic and Hyperbolic Differential Equations, J. Soc. Indust. Appl. Math., vol. 4, pp. 20–37, 1956.

14.   ——— and H. Rachford, On the Numerical Solution of Heat Conduction Problems in Two and Three Space Variables, Trans. Amer. Math. Soc., vol. 82, pp. 421–439, 1956.

15.   Ehrlich, L. W., A Numerical Method of Solving a Heat Flow Problem with Moving Boundary, J. Assoc. Comput. Mach., vol. 5, pp. 161–176, 1958.

16.   Forsythe, George E., What Are Relaxation Methods? chap. 17 in “Modern Mathematics for the Engineer,” First Series, edited by E. F. Beckenbach, McGraw-Hill Book Company, Inc., New York, 1956.

17.   Frankel, S., Convergence Rates of Iterative Treatments of Partial Differential Equations, Math. Tables Aids Comput., vol. 4, pp. 65–75, 1950.

18.   Friedman, Bernard, “The Iterative Solution of Elliptic Difference Equations,” A.E.C. Research and Development Report NYO-7698, Institute of Mathematical Sciences, New York University, June 1, 1957.

19.   Garabedian, P. R., Estimation of the Relaxation Factor for Small Mesh Size, Math. Tables Aids Comput., vol. 10, pp. 183–185, 1956.

20.   Geiringer, H., On the Solution of Systems of Linear Equations by Certain Iteration Methods, pp. 365–393 in “Reissner Anniversary Volume, Contributions to Applied Mechanics,” J. W. Edwards, Publisher, Inc., Ann Arbor, Mich., 1949.

21.   Gerschgorin, S., Fehlerabschätzung für das Differenzenverfahren zur Lösung partieller Differentialgleichungen, Z. Angew. Math. Mech., vol. 10, pp. 373–382, 1930.

22.   Hicks, Bruce L., J. W. Kelso, and Julian Davis, “Mathematical Theory of the Ignition Process Considered as a Thermal Reaction,” Ballistic Research Laboratories Report No. 756, Aberdeen Proving Ground, Md., 1957.

23.   John, F., On Integration of Parabolic Equations by Difference Methods, Comm. Pure Appl. Math., vol. 5, pp. 155–211, 1952.

24.   Juncosa, M. L., and David Young, On the Order of Convergence of Solutions of a Difference Equation to a Solution of the Diffusion Equation, J. Soc. Indust. Appl. Math., vol. 1, pp. 111–135, 1953.

25.   ——— and ———, On the Crank-Nicolson Procedure for Solving Parabolic Partial Differential Equations, Proc. Cambridge Philos. Soc., vol. 53, part 2, pp. 448–461, 1957.

26.   Kahan, W., “The Rate of Convergence of the Extrapolated Gauss-Seidel Iteration,” abstract of paper presented at the Conference on Matrix Computations, Wayne State University, Detroit, Mich., Sept. 4, 1957.

27.   Laasonen, P., On the Degree of Convergence of Discrete Approximations for the Solutions of the Dirichlet Problem, Ann. Acad. Sci. Fenn., Ser. A. I., vol. 246, pp. 1–19, 1957.

28.   ———, On the Solution of Poisson’s Difference Equation, J. Assoc. Comput. Mach., vol. 5, pp. 370–382, 1958.

29.   Lanczos, C., Solution of Systems of Linear Equations by Minimized Iterations, J. Res. Nat. Bur. Standards, vol. 49, pp. 33–53, 1952.

30.   Lax, P. D., and R. D. Richtmyer, Survey of Stability of Linear Finite Difference Equations, Comm. Pure Appl. Math., vol. 9, pp. 267–293, 1956.

31.   Lehmer, Derrick H., High-speed Computing Devices and Their Applications, chap. 19 in “Modern Mathematics for the Engineer,” First Series, edited by E. F. Beckenbach, McGraw-Hill Book Company, Inc., New York, 1956.

32.   Leutert, W. W., On the Convergence of Unstable Approximate Solutions of the Heat Equation to the Exact Solution, J. Math. Phys., vol. 30, pp. 245–251, 1952.

33.   Liebmann, H., Die angenährte Ermittelung harmonischer Funktionen und konformer Abbildung, Sitzungsberichte der mathematische-naturwissenschaftlichen Klasse der bayerischen Akademie der Wissenschaften zu München, 1918, pp. 385–416.

34.   O’Brien, G. G., M. A. Hyman, and S. Kaplan, A Study of the Numerical Solution of Partial Differential Equations, J. Math. Phys., vol. 29, pp. 223–251, 1951.

35.   Peaceman, D. W., and H. H. Rachford, Jr., The Numerical Solution of Parabolic and Elliptic Differential Equations, J. Soc. Indust. Appl. Math., vol. 3, pp. 28–41, 1955.

36.   Richardson, L. F., The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving Differential Equations, with Application to the Stresses in a Masonry Dam, Philos. Trans. Roy. Soc. London, Ser. A., vol. 210, pp. 307–357, 1910.

37.   Riabouchinsky, D., On Steady Fluid Motions with Free Surface, Proc. London Math. Soc., vol. 19, pp. 206–215, 1920–1921.

38.   Rosenbloom, P. C., “The Difference Equation Method for Solving the Dirichlet Problem,” Nat. Bur. Standards Applied Mathematics Series, vol. 18, 1952

39.   Schiffer, Menahem M., Boundary-value Problems in Elliptic Partial Differential Equations, chap. 6 in “Modern Mathematics for the Engineer,” First Series, edited by E. F. Beckenbach, McGraw-Hill Book Company, Inc., New York, 1956.

40.   Sheldon, J., On the Numerical Solution of Elliptic Difference Equations, Math. Tables Aids Comput., vol. 9, pp. 101–112, 1955.

41.   Shortley, G., Use of Tschebyscheff-polynomial Operators in the Numerical Solution of Boundary-value Problems, J. Appl. Phys., vol. 24, pp. 392–396, 1953.

42.   ——— and R. Weller, The Numerical Solution of Laplace’s Equation, J. Appl. Phys., vol. 9, pp. 334–344, 1938.

43.   ———, ———, Paul Darby, and Edward H. Gamble, Numerical Solution of Axisymmetrical Problems with Applications to Electrostatics and Torsion, J. Appl. Phys., vol. 18, pp. 116–129, 1947.

44.   Sokolnikoff, Ivan S., The Elastostatic Boundary-value Problems, chap. 7 in “Modern Mathematics for the Engineer,” First Series, edited by E. F. Beckenbach, McGraw-Hill Book Company, Inc., New York, 1956.

45.   Southwell, R. V., “Relaxation Methods in Theoretical Physics,” Oxford University Press, New York, 1946.

46.   Stiefel, E., Recent Developments in Relaxation Techniques, in “Proceedings of the International Congress of Mathematicians,” vol. 1, Amsterdam, 1954.

47.   Thomas, L. H., Elliptic Problems in Linear Difference Equations over a Network, unpublished manuscript, Watson Scientific Computing Laboratory.

48.   Todd, John, A Direct Approach to the Problem of Stability in the Numerical Solution of Partial Differential Equations, Comm. Pure Appl. Math., vol. 9, pp. 597–612, 1956.

49.   Tompkins, Charles B., Methods of Steep Descent, chap. 18 in “Modern Mathematics for the Engineer,” First Series, edited by E. F. Beckenbach, McGraw-Hill Book Company, Inc., New York, 1956.

50.   Varga, Richard S., “Numerical Solution of the Two-group Diffusion Equation in x-y Geometry,” Report WAPD-159, Bettis Plant, Westinghouse Electric Corp., Pittsburgh, Pa., August, 1956.

51.   Wachspress, E. L., “CURE: A Generalized Two-space-dimension Multigroup Coding for the IBM-704,” Report KAPL-1724, Knolls Atomic Power Laboratory, General Electric Company, Schenectady, N.Y., Apr. 30, 1957.

52.   ——— and G. J. Habetler, An Alternating Direction-implicit Iteration Technique, J. Soc. Indust. Appl. Math., vol. 8, pp. 403–424, 1960.

53.   Walsh, J. L., and David Young, On the Accuracy of the Numerical Solution of the Dirichlet Problem by Finite Differences, J. Research Nat. Bur. Standards, vol. 51, pp. 343–363, 1953.

54.   ——— and ———, On the Degree of Convergence of Solutions of Difference Equations to the Solution of the Dirichlet Problem, J. Math. Phys., vol. 33, pp. 80–93, 1954.

55.   Wasow, W., On the Truncation Error in the Solution of Laplace’s Equation by Finite Differences, J. Res. Nat. Bur. Standards, vol. 48, pp. 345–348, 1952.

56.   ———, The Accuracy of Difference Approximations to Plane Dirichlet Problems with Piecewise Analytic Boundary Values, Quart. Appl. Math., vol. 15, pp. 53–63, 1957.

57.   Young, David, Iterative Methods for Solving Partial Difference Equations of Elliptic Type, Trans. Amer. Math. Soc., vol. 76, pp. 92–111, 1954.

58.   ———, On Richardson’s Method for Solving Linear Systems with Positive Definite Matrices, J. Math. Phys., vol. 32, pp. 254–255, 1954.

59.   ———, Ordvac Solutions of the Dirichlet Problem, J. Assoc. Comput. Mach., vol. 2, pp. 137–161, 1955.

60.   ——— and Louis Ehrlich, Some Numerical Studies of Iterative Methods for Solving Elliptic Difference Equations, pp. 143–162 in “Boundary Problems in Differential Equations,” University of Wisconsin Press, Madison, Wis., 1960.

61.   ——— and ———, “On the Numerical Solution of Linear and Non-linear Parabolic Equations on the Ordvac,” Interim Technical Report No. 18, Office of Ordnance Research Contract DA-36-034-ORD-1486, University of Maryland, February, 1956.

62.   ———, L. D. Gates, Jr., R. J. Arms, and D. F. Eliezer, “The Computation of an Axially Symmetric Free Boundary Problem on NORC,” U.S. Naval Proving Ground Report No. 1413, Dahlgren, Va., Dec. 16, 1955.

63.   ——— and Charles H. Warlick, “On the Use of Richardson’s Method for the Numerical Solution of Laplace’s Equation on the ORDVAC,” Ballistic Research Laboratories Memorandum Report No. 707, Aberdeen Proving Ground, Md., July, 1953.

† This quantity is sometimes referred to as the spectral norm, for instance in Ref. 57. The spectral radius is a norm in the mathematical sense for linear transformations associated with symmetric matrices, but not for general linear transformations.

† The SPADE Committee was made up of representatives of several computing installations, primarily on the West Coast. The organizations represented included the RAND Corporation, Space Technology Laboratories, Lockheed Aircraft Corporation, North American Aviation, AiResearch, General Electric Company, International Business Machines Corporation, California Research Corporation, and The University of Texas.