INTRODUCTION

In this section we approximate a solution of the one-dimensional wave equation using the finite difference method used in the preceding two sections. The one-dimensional wave equation is the prototype hyperbolic partial differential equation.

Difference Equation Replacement

Suppose u(x, t) represents a solution of the one-dimensional wave equation

(1)

Using two central differences,

,

we replace equation (1) by

(2)

We solve (2) for u(x, t + k), which is u_{i, j + 1}. If λ = ck/h, then (2) yields

(3)

for i = 1, 2, …, n − 1 and j = 1, 2, …, m − 1.

In the case when the wave equation (1) is a model for the vertical displacements u(x, t) of a vibrating string, typical boundary conditions are u(0, t) = 0, u(a, t) = 0, t > 0, and initial conditions are u(x, 0) = f(x), ∂u/∂t|_t=0 = g(x), 0 < x < a. The functions f and g can be interpreted as the initial position and the initial velocity of the string. The numerical method based on equation (3), like the first method considered in Section 16.2, is an explicit finite difference method. As before, we use the difference equation (3) to approximate the solution u(x, t) of (1), using the boundary and initial conditions, over a rectangular region in the xt-plane defined by the inequalities 0 ≤ x ≤ a, 0 ≤ t ≤ T, where T is some specified value of time. If n and m are positive integers and

the vertical and horizontal grid lines on this region are defined by

x_i = ih, i = 0, 1, 2, …, n     and     t_j = jk,     j = 0, 1, 2, …, m.

As shown in FIGURE 16.3.1, (3) enables us to obtain the approximation u_{i, j+1} on the (j + 1)st time line from the values indicated on the jth and (j − 1)st time lines. Moreover, we use

and .

There is one minor problem in getting started. You can see from (3) that for j = 1 we need to know the values of u_{i, 1} (that is, the estimates of u on the first time line) in order to find u_{i, 2}. But from Figure 16.3.1, with j = 0, we see that the values of u_{i, 1} on the first time line depend on the values of u_{i, 0} on the zeroth time line and on the values of u_{i, −1}. To compute these latter values we make use of the initial-velocity condition u_t(x, 0) = g(x). At t = 0, it follows from (5) of Section 6.5 that

(4)

In order to make sense of the term u(x_i, − k) = u_{i, −1} in (4), we have to imagine u(x, t) extended backward in time. It follows from (4) that

This last result suggests that we define

(5)

in the iteration of (3). By substituting (5) into (3) when j = 0, we get the special case

(6)

EXAMPLE 1 Using the Finite Difference Method

Approximate the solution of the boundary-value problem

using (3) with n = 5 and m = 20.

SOLUTION

We make the identifications c = 2, a = 1, and T = 1. With n = 5 and m = 20 we get h = = 0.2, k = = 0.05, and λ = 0.5. Thus, with g(x) = 0, equations (6) and (3) become, respectively,

(7)

(8)

For i = 1, 2, 3, 4, equation (7) yields the following values for the u_{i, 1} on the first time line:

(9)

Note that the results given in (9) were obtained from the initial condition u(x, 0) = sin πx. For example, u₂₀ = sin(0.2π), and so on. Now j = 1 in (8) gives

u_{i, 2} = 0.25u_{i+1, 1} + 1.5u_{i, 1} + 0.25u_{i−1, 1} − u_{i, 0},

and so for i = 1, 2, 3, 4 we get

u₁₂ = 0.25u₂₁ + 1.5u₁₁ + 0.25u₀₁ − u₁₀

u₂₂ = 0.25u₃₁ + 1.5u₂₁ + 0.25u₁₁ − u₂₀

u₃₂ = 0.25u₄₁ + 1.5u₃₁ + 0.25u₂₁ − u₃₀

u₄₂ = 0.25u₅₁ + 1.5u₄₁ + 0.25u₃₁ − u₄₀.

Using the boundary conditions, the initial conditions, and the data obtained in (9), we get from these equations the approximations for u on the second time line. These last results and an abbreviation of the remaining calculations are summarized in Table 16.3.1 on page 827.

TABLE 16.3.1 Explicit Difference Equation Approximation with h = 0.2, k = 0.05, λ = 0.5

It is readily verified that the exact solution of the BVP in Example 1 is u(x, t) = sin πx cos 2πt. With this function we can compare the exact results with the approximations. For example, some selected comparisons are given in Table 16.3.2. As you can see in the table, the approximations are in the same “ball park” as the exact values, but the accuracy is not particularly impressive. We can, however, obtain more accurate results. The accuracy of this algorithm varies with the choice of λ. Of course, λ is determined by the choice of the integers n and m, which in turn determine the values of the step sizes h and k. It can be proved that the best accuracy is always obtained from this method when the ratio λ = kc/h is equal to one—in other words, when the step in the time direction is k = h/c. For example, the choice n = 8 and m = 16 yields h = , k = , and λ = 1. The sample values listed in Table 16.3.3 clearly show the improved accuracy.

Stability

We note in conclusion that this explicit finite difference method for the wave equation is stable when λ ≤ 1 and unstable when λ > 1.

16.3 Exercises Answers to selected odd-numbered problems begin on page ANS-43.

In Problems 1, 3, 5, and 6, use a computer as a computational aid.

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

   

   when
   (a) c = 1, a = 1, T = 1, f(x) = x(1 − x); n = 4 and m = 10
   (b) c = 1, a = 2, T = 1, f(x) = ; n = 5 and m = 10
   (c)

   

   n = 10 and m = 25.

2. Consider the boundary-value problem

   

   1. Use the methods of Chapter 13 to verify that the solution of the problem is u(x, t) = sin πx cos πt.
   2. Use the method of this section to approximate the solution of the problem without the aid of a computer program. Use n = 4 and m = 5.
   3. Compute the absolute error at each interior grid point.
3. Approximate the solution of the boundary-value problem in Problem 2 using a computer program with
   (a) n = 5, m = 10
   (b) n = 5, m = 20.
4. Given the boundary-value problem

   

   use h = k = in equation (6) to compute the values of u_{i, 1} by hand.

5. It was shown in Section 13.2 that the equation of a vibrating string is

   

   where T is the constant magnitude of the tension in the string and ρ is its mass per unit length. Suppose a string of length 60 centimeters is secured to the x-axis at its ends and is released from rest from the initial displacement

   

   Use the difference equation (3) in this section to approximate the solution of the boundary-value problem when h = 10, k = 5 and where ρ = 0.0225 g/cm, T = 1.4 × 10⁷ dynes. Use m = 50.

6. Repeat Problem 5 using

   

   and h = 10, k = 2.5. Use m = 50.