INTRODUCTION

Recall from Section 13.1 that linear second-order PDEs in two independent variables are classified as elliptic, parabolic, and hyperbolic. Roughly, elliptic PDEs involve partial derivatives with respect to spatial variables only and as a consequence solutions of such equations are determined by boundary conditions alone. Parabolic and hyperbolic equations involve partial derivatives with respect to both spatial and time variables, and so solutions of such equations generally are determined from boundary and initial conditions. A solution of an elliptic PDE (such as Laplace’s equation) can describe a physical system whose state is in equilibrium (steady state), a solution of a parabolic PDE (such as the heat equation) can describe a diffusional state, whereas a hyperbolic PDE (such as the wave equation) can describe a vibrational state.

In this section we begin our discussion with approximation methods appropriate for elliptic equations. Our focus will be on the simplest but probably the most important PDE of the elliptic type: Laplace’s equation.

Difference Equation Replacement

Suppose that we are seeking a solution u(x, y) of Laplace’s equation

in a planar region R that is bounded by some curve C. See FIGURE 16.1.1. Analogous to (6) of Section 6.5, using the central differences

approximations for the second partial derivatives u_xx and u_yy can be obtained using the difference quotients

(1)

(2)

Now by adding (1) and (2) we obtain a five-point approximation to the Laplacian:

Hence we can replace Laplace’s equation by the difference equation

(3)

If we adopt the notation u(x, y) = u_ij and

then (3) becomes

(4)

To understand (4) a little better, suppose a rectangular grid consisting of horizontal lines spaced h units apart and vertical lines spaced h units apart is placed over the region R. The number h is called the mesh size. See FIGURE 16.1.2(a). The points P_ij = P(ih, jh), i and j integers, of intersection of the horizontal and vertical lines, are called mesh points or lattice points. A mesh point is an interior point if its four nearest neighboring mesh points are points of R. Points in R or on C that are not interior points are called boundary points. For example, in Figure 16.1.2(a) we have

P₂₀ = P(2h, 0), P₁₁ = P(h, h), P₂₁ = P(2h, h), P₂₂ = P(2h, 2h),

and so on. Of the points listed, P₂₁ and P₂₂ are interior points, whereas P₂₀ and P₁₁ are boundary points. In Figure 16.1.2(a) interior points are the dots shown in red and the boundary points are shown in black. Now from (4) we see that

(5)

and so, as shown in Figure 16.1.2(b), the value of u_ij at an interior mesh point of R is the average of the values of u at four neighboring mesh points. The neighboring points P_{i + 1, j}, P_{i, j + 1}, P_{i − 1, j}, and P_{i, j − 1} correspond, respectively, to the four points on a compass: E, N, W, and S.

Dirichlet Problem

Recall that in the Dirichlet problem for Laplace’s equation ∇²u = 0, the values of u(x, y) are prescribed on the boundary C of a region R. The basic idea is to find an approximate solution to Laplace’s equation at interior mesh points by replacing the partial differential equation at these points by the difference equation (4). Hence the approximate values of u at the mesh points—namely, the u_ij—are related to each other and, possibly, to known values of u if a mesh point lies on the boundary C. In this manner we obtain a system of linear algebraic equations that we solve for the unknown u_ij. The following example illustrates the method for a square region.

EXAMPLE 1 A Boundary-Value Problem Revisited

In Problem 16 of Exercises 13.5 you were asked to solve the boundary-value problem

utilizing the superposition principle. To apply the present numerical method, let us start with a mesh size of h = . As we see in FIGURE 16.1.3, that choice yields four interior points and eight boundary points. The numbers listed next to the boundary points are the exact values of u obtained from the specified condition along that boundary. For example, at P₃₁ = P(3h, h) = P(2, ) we have x = 2 and y = , and so the condition u(2, y) gives u(2, ) = (2 − ) = . Similarly, at P₁₃ = P(, 2), the condition u(x, 2) gives u(, 2) = . We now apply (4) at each interior point. For example, at P₁₁ we have i = 1 and j = 1, so (4) becomes

u₂₁ + u₁₂ + u₀₁ + u₁₀ − 4u₁₁ = 0.

Since u₀₁ = u(0, ) = 0 and u₁₀ = u(, 0) = 0, the foregoing equation becomes −4u₁₁ + u₂₁ + u₁₂ = 0. Repeating this, in turn, at P₂₁, P₁₂, and P₂₂, we get three additional equations:

(6)

Using a computer algebra system to solve the system, we find the approximate temperatures at the four interior points to be

≡

As in the discussion of ordinary differential equations, we expect that a smaller value of h will improve the accuracy of the approximation. However, using a smaller mesh size means, of course, that there are more interior mesh points, and correspondingly there is a much larger system of equations to be solved. For a square region whose length of side is L, a mesh size of h = L/n will yield a total of (n − 1)² interior mesh points. In Example 1, for n = 8, the mesh size is a reasonable h = = , but the number of interior points is (8 − 1)² = 49. Thus we have 49 equations in 49 unknowns. In the next example we use a mesh size of h = .

EXAMPLE 2 Example 1 with More Mesh Points

As we see in FIGURE 16.1.4, with n = 4, a mesh size h = = for the square in Example 1 gives 3² = 9 interior mesh points. Applying (4) at these points and using the indicated boundary conditions, we get nine equations in nine unknowns. So that you can verify the results, we give the system in an unsimplified form:

(7)

In this case, a CAS yields

After we simplify (7), it is interesting to note that the 9 × 9 matrix of coefficients is

(8)

This is an example of a sparse matrix in that a large percentage of the entries are zeros. The matrix (8) is also an example of a banded matrix. These kinds of matrices are characterized by the properties that the entries on the main diagonal and on diagonals (or bands) parallel to the main diagonal are all nonzero. The bands shown in red in (8) are separated by diagonals consisting of all zeros or not.

Gauss–Seidel Iteration

Problems requiring approximations to solutions of partial differential equations invariably lead to large systems of linear algebraic equations. It is not uncommon to have to solve systems involving hundreds of equations. Although a direct method of solution such as Gaussian elimination leaves unchanged the zero entries outside the bands in a matrix such as (8), it does fill in the positions between the bands with nonzeros. Since storing very large matrices uses up a large portion of computer memory, it is usual practice to solve a large system in an indirect manner. One popular indirect method is called Gauss–Seidel iteration. This method is named after the German mathematicians Johann Carl Friedrich Gauss (1777–1855) and Philipp Ludwig von Seidel (1821–1896).

We shall illustrate this method for the system in (6). For the sake of simplicity we replace the double-subscripted variables u₁₁, u₂₁, u₁₂, and u₂₂ by x₁, x₂, x₃, and x₄, respectively.

EXAMPLE 3 Gauss–Seidel Iteration

Step 1: Solve each equation for the variables on the main diagonal of the system. That is, in (6), solve the first equation for x₁, the second equation for x₂, and so on:

(9)

These equations can be obtained directly by using (5) rather than (4) at the interior points.

Step 2: Iterations. We start by making an initial guess for the values of x₁, x₂, x₃, and x₄. If this were simply a system of linear equations and we knew nothing about the solution, we could start with x₁ = 0, x₂ = 0, x₃ = 0, x₄ = 0. But since the solution of (9) represents approximations to a solution of a boundary-value problem, it would seem reasonable to use as the initial guess for the values of x₁ = u₁₁, x₂ = u₂₁, x₃ = u₁₂, and x₄ = u₂₂ the average of all the boundary conditions. In this case the average of the numbers at the eight boundary points shown in Figure 16.1.2 is approximately 0.4. Thus our initial guess is x₁ = 0.4, x₂ = 0.4, x₃ = 0.4, and x₄ = 0.4. Iterations of the Gauss–Seidel method use the x values as soon as they are computed. Note that the first equation in (9) depends only on x₂ and x₃; thus substituting x₂ = 0.4 and x₃ = 0.4 gives x₁ = 0.2. Since the second and third equations depend on x₁ and x₄, we use the newly calculated values x₁ = 0.2 and x₄ = 0.4 to obtain x₂ = 0.3722 and x₃ = 0.3167. The fourth equation depends on x₂ and x₃, so we use the new values x₂ = 0.3722 and x₃ = 0.3167 to get x₄ = 0.5611. In summary, the first iteration has given the values

Note how close these numbers are already to the actual values given at the end of Example 1.

The second iteration starts with substituting x₂ = 0.3722 and x₃ = 0.3167 into the first equation. This gives x₁ = 0.1722. From x₁ = 0.1722 and the last computed value of x₄ (namely, x₄ = 0.5611), the second and third equations give, in turn, x₂ = 0.4055 and x₃ = 0.3500. Using these two values, we find from the fourth equation that x₄ = 0.5678. At the end of the second iteration we have

The third through seventh iterations are summarized in Table 16.1.1.

TABLE 16.1.1 Iteration

To apply Gauss–Seidel iteration to a general system of n linear equations in n unknowns, the variable x_i must actually appear in the ith equation of the system. Moreover, after each equation is solved for x_i, i = 1, 2, …, n, the resulting system has the form X = AX + B, where all the entries on the main diagonal of A are zero.

REMARKS

(i) In the examples given in this section, the values of u_ij were determined using known values of u at boundary points. But what do we do if the region is such that boundary points do not coincide with the actual boundary C of the region R? In this case the required values can be obtained by interpolation.

(ii) It is sometimes possible to cut down the number of equations to solve by using symmetry. Consider the rectangular region defined by 0 ≤ x ≤ 2, 0 ≤ y ≤ 1, shown in FIGURE 16.1.5. The boundary conditions are u = 0 along the boundaries x = 0, x = 2, y = 1, and u = 100 along y = 0. The region is symmetric about the lines x = 1 and y = , and the interior points P₁₁ and P₃₁ are equidistant from the neighboring boundary points at which the specified values of u are the same. Consequently we assume that u₁₁ = u₃₁, and so the system of three equations in three unknowns reduces to two equations in two unknowns. See Problem 2 in Exercises 16.1.

(iii) It may not be noticeable on a computer, but convergence of Gauss–Seidel iteration (or Liebman’s method) may not be particularly fast. Also, in a more general setting, Gauss–Seidel iteration may not converge at all. For conditions that are sufficient to guarantee convergence of Gauss–Seidel iteration you are encouraged to consult texts on numerical analysis.

16.1 Exercises Answers to selected odd-numbered problems begin on page ANS-40.

In Problems 1–8, use a computer as a computational aid.

In Problems 1–4, use (4) to approximate the solution of Laplace’s equation at the interior points of the given region. Use symmetry when possible.

1. u(0, y) = 0, u(3, y) = y(2 − y), 0 < y < 2
   u(x, 0) = 0, u(x, 2) = x(3 − x), 0 < x < 3
   mesh size: h = 1
2. u(0, y) = 0, u(2, y) = 0, 0 < y < 1
   u(x, 0) = 100, u(x, 1) = 0, 0 < x < 2
   mesh size: h =
3. u(0, y) = 0, u(1, y) = 0, 0 < y < 1
   u(x, 0) = 0, u(x, 1) = sin πx, 0 < x < 1
   mesh size: h =
4. u(0, y) = 108y²(1 − y), u(1, y) = 0, 0 < y < 1
   u(x, 0) = 0, u(x, 1) = 0, 0 < x < 1
   mesh size: h =

In Problems 5 and 6, use (5) and Gauss–Seidel iteration to approximate the solution of Laplace’s equation at the interior points of a unit square. Use the mesh size h = . In Problem 5, the boundary conditions are given; in Problem 6, the values of u at boundary points are given in FIGURE 16.1.6.

5. u(0, y) = 0, u(1, y) = 100y, 0 < y < 1
   u(x, 0) = 0, u(x, 1) = 100x, 0 < x < 1

6. 

7. 1. In Problem 12 of Exercises 13.6, you solved a potential problem using a special form of Poisson’s equation . Show that the difference equation replacement for Poisson’s equation is

      u_{i + 1, j} + u_{i, j + 1} + u_{i − 1, j} + u_{i, j − 1} − 4u_ij = h²f(x, y).

   2. Use the result in part (a) to approximate the solution of the Poisson equation = − 2 at the interior points of the region in FIGURE 16.1.7. The mesh size is h = , u = 1 at every point along ABCD, and u = 0 at every point along DEFGA. Use symmetry and, if necessary, Gauss–Seidel iteration.

   

8. Use the result in part (a) of Problem 7 to approximate the solution of the Poisson equation = − 64 at the interior points of the region in FIGURE 16.1.8. The mesh is h = , and u = 0 at every point on the boundary of the region. If necessary, use Gauss–Seidel iteration.