INTRODUCTION

The basic idea in the following discussion is the same as in Section 16.1; we approximate a solution of a PDE—this time a parabolic PDE—by replacing the equation with a finite difference equation. But unlike the preceding section, we shall consider two finite-difference approximation methods for parabolic partial differential equations: one called an explicit method and the other an implicit method. For the sake of definiteness, we consider only the one-dimensional heat equation.

Difference Equation Replacement

To approximate the solution u(x, t) of the one-dimensional heat equation

(1)

we again replace the derivatives by difference quotients. Using the central difference approximation (2) of Section 16.1,

and the forward difference approximation (3) of Section 6.5,

equation (1) becomes

(2)

If we let λ = ck/h² and

then, after simplifying, (2) is

(3)

In the case of the heat equation (1), typical boundary conditions are u(0, t) = u₁, u(a, t) = u₂, t > 0, and an initial condition is u(x, 0) = f(x), 0 < x < a. The function f can be interpreted as the initial temperature distribution in a homogeneous rod extending from x = 0 to x = a; u₁ and u₂ can be interpreted as constant temperatures at the endpoints of the rod. Although we shall not prove it, the boundary-value problem consisting of (1) and these two boundary conditions and one initial condition has a unique solution when f is continuous on the closed interval [0, a]. This latter condition will be assumed, and so we replace the initial condition by u(x, 0) = f(x), 0 ≤ x ≤ a. Moreover, instead of working with the semi-infinite region in the xt-plane defined by the inequalities 0 ≤ x ≤ a, t ≥ 0, we use a rectangular region defined by 0 ≤ x ≤ a, 0 ≤ t ≤ T, where T is some specified value of time. Over this region we place a rectangular grid consisting of vertical lines h units apart and horizontal lines k units apart. See FIGURE 16.2.1. If we choose two positive integers n and m and define

then the vertical and horizontal grid lines are defined by

As illustrated in FIGURE 16.2.2, the plan here is to use formula (3) to estimate the values of the solution u(x, t) at the points on the (j + 1)st time line using only values from the jth time line. For example, the values on the first time line (j = 1) depend on the initial condition u_{i, 0} = u(x_i, 0) = f(x_i) given on the zeroth time line (j = 0). This kind of numerical procedure is called an explicit finite difference method.

EXAMPLE 1 Using the Finite Difference Method

Consider the boundary-value problem

First we identify c = 1, a = 1, and T = 0.5. If we choose, say, n = 5 and m = 50, then h = = 0.2, k = = 0.01, λ = 0.25,

Thus (3) becomes

u_{i, j + 1} = 0.25(u_{i + 1, j} + 2u_ij + u_{i + 1, j}).

By setting j = 0 in this formula, we get a formula for the approximations to the temperature u on the first time line:

u_{i, 1} = 0.25(u_{i + 1, 0} + 2u_{i, 0} + u_{i + 1, 0}).

If we then let i = 1, …, 4 in the last equation, we obtain, in turn,

u₁₁ = 0.25(u₂₀ + 2u₁₀ + u₀₀)

u₂₁ = 0.25(u₃₀ + 2u₂₀ + u₁₀)

u₃₁ = 0.25(u₄₀ + 2u₃₀ + u₂₀)

u₄₁ = 0.25(u₅₀ + 2u₄₀ + u₃₀).

The first equation in this list is interpreted as

From the initial condition u(x, 0) = sin πx, the last line becomes

u₁₁ = 0.25(0.951056516 + 2(0.587785252) + 0) = 0.531656755.

This number represents an approximation to the temperature u(0.2, 0.01).

Since it would require a rather large table of over 200 entries to summarize all the approximations over the rectangular grid determined by h and k, we give only selected values in Table 16.2.1.

TABLE 16.2.1 Explicit Difference Equation Approximation with h = 0.2, k = 0.01, λ = 0.25

You should verify, using the methods of Chapter 13, that an exact solution of the boundary-value problem in Example 1 is given by u(x, t) = sin πx. Using this solution, we compare in Table 16.2.2 a sample of exact values with their corresponding approximations.

Stability

These approximations are comparable to the exact values and accurate enough for some purposes. But there is a problem with the foregoing method. Recall that a numerical method is unstable if round-off errors or any other errors grow too rapidly as the computations proceed. The numerical procedure in Example 1 can exhibit this kind of behavior. It can be proved that the procedure is stable if λ is less than or equal to 0.5, but unstable otherwise. To obtain λ = 0.25 ≤ 0.5 in Example 1, we had to choose the value k = 0.01; the necessity of using very small step sizes in the time direction is the principal fault of this method. You are urged to work Problem 12 in Exercises 16.2 and witness the predictable instability when λ = 1.

Crank–Nicolson Method

There are implicit finite difference methods for solving parabolic partial differential equations. These methods require that we solve a system of equations to determine the approximate values of u on the (j + 1)st time line. However, implicit methods do not suffer from instability problems.

The finite difference algorithm introduced by the English mathematicians John Crank (1916–2006) and Phyllis Nicolson (1917–1968) in 1947 is used mostly for solving the heat equation or similar parabolic PDEs. The algorithm consists of replacing the second partial derivative in by an average of two central difference quotients, one evaluated at t and the other at t + k:

(4)

If we again define λ = ck/h², then after rearranging terms we can write (4) as

(5)

where α = 2(1 + 1/λ) and β = 2(1 − 1/λ), j = 0, 1, …, m − 1, and i = 1, 2, …, n − 1.

For each choice of j, the difference equation (5) for i = 1, 2, …, n − 1 gives n − 1 equations in n − 1 unknowns u_{i, j + 1}. Because of the prescribed boundary conditions, the values of u_{i, j + 1} are known for i = 0 and for i = n. For example, in the case n = 4, the system of equations for determining the approximate values of u on the (j + 1)st time line is

−u_{0, j + 1} + αu_{1, j + 1} − u_{2, j + 1} = u_{2, j} − βu_{1, j} + u_{0, j}

−u_{1, j + 1} + αu_{2, j + 1} − u_{3, j + 1} = u_{3, j} − βu_{2, j} + u_{1, j}

−u_{2, j + 1} + αu_{3, j + 1} − u_{4, j + 1} = u_{4, j} − βu_{3, j} + u_{2, j}

or (6)

where .

In general, if we use the difference equation (5) to determine values of u on the (j + 1)st time line, we need to solve a linear system AX = B, where the coefficient matrix A is a tridiagonal matrix,

,

and the entries of the column matrix B are

EXAMPLE 2 Using the Crank–Nicolson Method

Use the Crank–Nicolson method to approximate the solution of the boundary-value problem

using n = 8 and m = 30.

SOLUTION

From the identifications a = 2, T = 0.3, h = = 0.25, k = = 0.01, and c = 0.25, we get λ = 0.04. With the aid of a computer we get the results in Table 16.2.3. As in Example 1, the entries in this table represent only a selected number from the 210 approximations over the rectangular grid determined by h and k.

TABLE 16.2.3 Crank–Nicolson Method with h = 0.25, k = 0.01, λ = 0.25

Like Example 1, the boundary-value problem in Example 2 also possesses an exact solution given by u(x, t) = sin πx. The sample comparisons listed in Table 16.2.4 show that the absolute errors are of the order 10⁻² or 10⁻³. Smaller errors can be obtained by decreasing either h or k.

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

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

1. Use the difference equation (3) to approximate the solution of the boundary-value problem

   

   Use n = 8 and m = 40.

2. Using the Fourier series solution obtained in Problem 3 of Exercises 13.3 with L = 2, one can sum the first 20 terms to estimate the values for u(0.25, 0.1), u(1, 0.5), and u(1.5, 0.8) for the solution u(x, t) of Problem 1 above. A student wrote a computer program to do this and obtained the results u(0.25, 0.1) = 0.3794, u(1, 0.5) = 0.1854, and u(1.5, 0.8) = 0.0623. Assume these results are accurate for all digits given. Compare these values with the approximations obtained in Problem 1 above. Find the absolute errors in each case.
3. Solve Problem 1 by the Crank–Nicolson method with n = 8 and m = 40. Use the values for u(0.25, 0.1), u(1, 0.5), and u(1.5, 0.8) given in Problem 2 to compute the absolute errors.
4. Repeat Problem 1 using n = 8 and m = 20. Use the values for u(0.25, 0.1), u(1, 0.5), and u(1.5, 0.8) given in Problem 2 to compute the absolute errors. Why are the approximations so inaccurate in this case?
5. Solve Problem 1 by the Crank–Nicolson method with n = 8 and m = 20. Use the values for u(0.25, 0.1), u(1, 0.5), and u(1.5, 0.8) given in Problem 2 to compute the absolute errors. Compare the absolute errors with those obtained in Problem 4.
6. It was shown in Section 13.2 that if a rod of length L is made of a material with thermal conductivity K, specific heat γ, and density ρ, the temperature u(x, t) satisfies the partial differential equation

   

   Consider the boundary-value problem consisting of the foregoing equation and conditions

   

   Use the difference equation (3) in this section with n = 10 and m = 10 to approximate the solution of the boundary-value problem when

   1. L = 20, K = 0.15, ρ = 8.0, γ = 0.11, f(x) = 30
   2. L = 50, K = 0.15, ρ = 8.0, γ = 0.11, f(x) = 30
   3. L = 20, K = 1.10, ρ = 2.7, γ = 0.22, f(x) = 0.5x(20 − x)
   4. L = 100, K = 1.04, ρ = 10.6, γ = 0.06,

   

7. Solve Problem 6 by the Crank–Nicolson method with n = 10 and m = 10.
8. Repeat Problem 6 if the endpoint temperatures are u(0, t) = 0, u(L, t) = 20, 0 ≤ t ≤ 10.
9. Solve Problem 8 by the Crank–Nicolson method.
10. Consider the boundary-value problem in Example 2. Assume that n = 4.
    1. Find the new value of λ.
    2. Use the Crank–Nicolson difference equation (5) to find the system of equations for u₁₁, u₂₁, and u₃₁, that is, the approximate values of u on the first time line. [Hint: Set j = 0 in (5), and let i take on the values 1, 2, 3.]
    3. Solve the system of three equations without the aid of a computer program. Compare your results with the corresponding entries in Table 16.2.3.
11. Consider a rod whose length is L = 20 for which K = 1.05, ρ = 10.6, and γ = 0.056. Suppose

    

    1. Use the method outlined in Section 13.6 to find the steady-state solution ψ(x).
    2. Use the Crank–Nicolson method to approximate the temperatures u(x, t) for 0 ≤ t ≤ T_max. Select T_max large enough to allow the temperatures to approach the steady-state values. Compare the approximations for t = T_max with the values of ψ(x) found in part (a).
12. Use the difference equation (3) to approximate the solution of the boundary-value problem

    

    Use n = 5 and m = 25.