Lecture AIMS Notes 2006
Peter J. Olver
11. Numerical Solution of the Heat and Wave Equations
In this part, we study numerical solution methodss for the two most important equations of one-dimensional continuum dynamics. The heat equation models the diffusion of thermal energy in a body; here, we treat the case of a one-dimensional bar. The wave equation describes vibrations and waves in continuous media, including sound waves, water waves, elastic waves, electromagnetic waves, and so on. For simplicity, we restrict our attention to the case of waves in a one-dimensional medium, e.g., a string, bar, or column of air. We begin with a general discussion of finite difference formulae for numerically approximating derivatives of functions. The basic finite difference scheme is obtained by replacing the derivatives in the equation by the appropriate numerical differentiation formulae. However, there is no guarantee that the resulting numerical scheme will accurately approximate the true solution, and further analysis is required to elicit bona fide, convergent numerical algorithms. In dynamical problems, the finite difference schemes replace the partial differential equation by an iterative linear matrix system, and the analysis of convergence relies on the methods covered in Section 7.1. We will only introduce the most basic algorithms, leaving more sophisticated variations and extensions to a more thorough treatment, which can be found in numerical analysis texts, e.g., [5, 7, 28].
11.1. Finite Differences.
In general, to approximate the derivative of a function at a point, say f (x) or f (x), one constructs a suitable combination of sampled function values at nearby points. The underlying formalism used to construct these approximation formulae is known as the calculus of finite differences. Its development has a long and influential history, dating back to Newton. The resulting finite difference numerical methods for solving differential equations have extremely broad applicability, and can, with proper care, be adapted to most problems that arise in mathematics and its many applications. The simplest finite difference approximation is the ordinary difference quotient u(x + h) - u(x) u (x), h 4/20/07 186
c 2006
(11.1)
Peter J. Olver
One-Sided Difference Figure 11.1.
Central Difference Finite Difference Approximations.
used to approximate the first derivative of the function u(x). Indeed, if u is differentiable at x, then u (x) is, by definition, the limit, as h 0 of the finite difference quotients. Geometrically, the difference quotient equals the slope of the secant line through the two points x, u(x) and x + h, u(x + h) on the graph of the function. For small h, this should be a reasonably good approximation to the slope of the tangent line, u (x), as illustrated in the first picture in Figure 11.1. How close an approximation is the difference quotient? To answer this question, we assume that u(x) is at least twice continuously differentiable, and examine the first order Taylor expansion (11.2) u(x + h) = u(x) + u (x) h + 1 u () h2 . 2 We have used the Cauchy formula for the remainder term, in which represents some point lying between x and x + h. The error or difference between the finite difference formula and the derivative being approximated is given by u(x + h) - u(x) - u (x) = h
1 2
u () h.
(11.3)
Since the error is proportional to h, we say that the finite difference quotient (11.3) is a first order approximation. When the precise formula for the error is not so important, we will write u(x + h) - u(x) + O(h). (11.4) u (x) = h The "big Oh" notation O(h) refers to a term that is proportional to h, or, more rigorously, bounded by a constant multiple of h as h 0. Example 11.1. Let u(x) = sin x. Let us try to approximate u (1) = cos 1 = 0.5403023 . . . by computing finite difference quotients cos 1 sin(1 + h) - sin 1 . h
The result for different values of h is listed in the following table. 4/20/07 187
c 2006 Peter J. Olver
h approximation error
1 0.067826 -0.472476
.1 0.497364 -0.042939
.01 0.536086 -0.004216
.001 0.539881 -0.000421
.0001 0.540260 -0.000042
1 We observe that reducing the step size by a factor of 10 reduces the size of the error by approximately the same factor. Thus, to obtain 10 decimal digits of accuracy, we anticipate needing a step size of about h = 10-11 . The fact that the error is more of less proportional to the step size confirms that we are dealing with a first order numerical approximation.
To approximate higher order derivatives, we need to evaluate the function at more than two points. In general, an approximation to the nth order derivative u(n) (x) requires at least n+1 distinct sample points. For simplicity, we shall only use equally spaced points, leaving the general case to the exercises. For example, let us try to approximate u (x) by sampling u at the particular points x, x + h and x - h. Which combination of the function values u(x - h), u(x), u(x + h) should be used? The answer to such a question can be found by consideration of the relevant Taylor expansions u(x + h) = u(x) + u (x) h + u (x) h2 h3 + u (x) + O(h4 ), 2 6 h2 h3 u(x - h) = u(x) - u (x) h + u (x) - u (x) + O(h4 ), 2 6 u(x + h) + u(x - h) = 2 u(x) + u (x) h2 + O(h4 ). Rearranging terms, we conclude that u(x + h) - 2 u(x) + u(x - h) + O(h2 ), (11.6) h2 The result is known as the centered finite difference approximation to the second derivative of a function. Since the error is proportional to h2 , this is a second order approximation. u (x) = Example 11.2. Let u(x) = ex , with u (x) = (4 x2 + 2) ex . Let us approximate u (1) = 6 e = 16.30969097 . . . by using the finite difference quotient (11.6): e(1+h) - 2 e + e(1-h) . 6e h2 The results are listed in the following table.
h approximation error 1 50.16158638 33.85189541 .1 16.48289823 0.17320726 .01 16.31141265 0.00172168 .001 16.30970819 0.00001722 .0001 16.30969115 0.00000018
2 2 2 2
(11.5)
where the error terms are proportional to h4 . Adding the two formulae together gives
1 Each reduction in step size by a factor of 10 reduces the size of the error by a factor of 1 and results in a gain of two new decimal digits of accuracy, confirming that the finite 100 difference approximation is of second order.
4/20/07
188
c 2006
Peter J. Olver
However, this prediction is not completely borne out in practice. If we take h = .00001 then the formula produces the approximation 16.3097002570, with an error of 0.0000092863 -- which is less accurate that the approximation with h = .0001. The problem is that round-off errors have now begun to affect the computation, and underscores the difficulty with numerical differentiation. Finite difference formulae involve dividing very small quantities, which can induce large numerical errors due to round-off. As a result, while they typically produce reasonably good approximations to the derivatives for moderately small step sizes, to achieve high accuracy, one must switch to a higher precision. In fact, a similar comment applied to the previous Example 11.1, and our expectations about the error were not, in fact, fully justified as you may have discovered if you tried an extremely small step size. Another way to improve the order of accuracy of finite difference approximations is to employ more sample points. For instance, if the first order approximation (11.4) to the first derivative based on the two points x and x + h is not sufficiently accurate, one can try combining the function values at three points x, x + h and x - h. To find the appropriate combination of u(x -h), u(x), u(x +h), we return to the Taylor expansions (11.5). To solve for u (x), we subtract the two formulae, and so u(x + h) - u(x - h) = 2 u (x) h + u (x) h3 + O(h4 ). 3
Rearranging the terms, we are led to the well-known centered difference formula u (x) = u(x + h) - u(x - h) + O(h2 ), 2h (11.7)
which is a second order approximation to the first derivative. Geometrically, the centered difference quotient represents the slope of the secant line through the two points x - h, u(x - h) and x + h, u(x + h) on the graph of u centered symmetrically about the point x. Figure 11.1 illustrates the two approximations; the advantages in accuracy in the centered difference version are graphically evident. Higher order approximations can be found by evaluating the function at yet more sample points, including, say, x + 2 h, x - 2 h, etc. Example 11.3. Return to the function u(x) = sin x considered in Example 11.1. The centered difference approximation to its derivative u (1) = cos 1 = 0.5403023 . . . is cos 1 The results are tabulated as follows: sin(1 + h) - sin(1 - h) . 2h
This next computation depends upon the computer's precision; here we used single precision in Matlab.
The terms O(h4 ) do not cancel, since they represent potentially different multiples of h4 .
4/20/07
189
c 2006
Peter J. Olver
h approximation error
.1 0.53940225217 -0.00090005370
.01 0.54029330087 -0.00000900499
.001 0.54030221582 -0.00000009005
.0001 0.54030230497 -0.00000000090
As advertised, the results are much more accurate than the one-sided finite difference approximation used in Example 11.1 at the same step size. Since it is a second order 1 approximation, each reduction in the step size by a factor of 10 results in two more decimal places of accuracy. Many additional finite difference approximations can be constructed by similar manipulations of Taylor expansions, but these few very basic ones will suffice for our subsequent purposes. In the following subsection, we apply the finite difference formulae to develop numerical solution schemes for the heat and wave equations.
11.2. Numerical Algorithms for the Heat Equation.
Consider the heat equation 2u u = , t x2 0<x< , t 0, (11.8)
representing a homogeneous diffusion process of, sqy, heat in bar of length and constant thermal diffusivity > 0. The solution u(t, x) represents the temperature in the bar at time t 0 and position 0 x . To be concrete, we will impose time-dependent Dirichlet boundary conditions u(t, 0) = (t), u(t, ) = (t), t 0, (11.9)
specifying the temperature at the ends of the bar, along with the initial conditions u(0, x) = f (x), 0x , (11.10)
specifying the bar's initial temperature distribution. In order to effect a numerical approximation to the solution to this initial-boundary value problem, we begin by introducing a rectangular mesh consisting of points (ti , xj ) with 0 = x0 < x1 < < xn = and 0 = t 0 < t1 < t2 < .
For simplicity, we maintain a uniform mesh spacing in both directions, with h = xj+1 - xj = n , k = ti+1 - ti ,
representing, respectively, the spatial mesh size and the time step size. It will be essential that we do not a priori require the two to be the same. We shall use the notation ui,j u(ti , xj ) where ti = i k, xj = j h, (11.11)
to denote the numerical approximation to the solution value at the indicated mesh point. 4/20/07 190
c 2006 Peter J. Olver
As a first attempt at designing a numerical method, we shall use the simplest finite difference approximations to the derivatives. The second order space derivative is approximated by (11.6), and hence u(ti , xj+1 ) - 2 u(ti , xj ) + u(ti , xj-1 ) 2u (ti , xj ) + O(h2 ) x2 h2 ui,j+1 - 2 ui,j + ui,j-1 + O(h2 ), h2
(11.12)
where the error in the approximation is proportional to h2 . Similarly, the one-sided finite difference approximation (11.4) is used for the time derivative, and so ui+1,j - ui,j u(ti+1 , xj ) - u(ti , xj ) u (11.13) (ti , xj ) + O(k) + O(k), t k k where the error is proportion to k. In practice, one should try to ensure that the approximations have similar orders of accuracy, which leads us to choose k h2 . Assuming h < 1, this requirement has the important consequence that the time steps must be much smaller than the space mesh size. Remark : At this stage, the reader might be tempted to replace (11.13) by the second order central difference approximation (11.7). However, this produces significant complications, and the resulting numerical scheme is not practical. Replacing the derivatives in the heat equation (11.14) by their finite difference approximations (11.12), (11.13), and rearranging terms, we end up with the linear system ui+1,j = ui,j+1 + (1 - 2 )ui,j + ui,j-1 , in which = i = 0, 1, 2, . . . , j = 1, . . . , n - 1, (11.14)
k . (11.15) h2 The resulting numerical scheme takes the form of an iterative linear system for the solution values ui,j u(ti , xj ), j = 1, . . . , n - 1, at each time step ti . The initial condition (11.10) means that we should initialize our numerical data by sampling the initial temperature at the mesh points: u0,j = fj = f (xj ), ui,0 = i = (ti ), j = 1, . . . , n - 1. (11.16)
Similarly, the boundary conditions (11.9) require that ui,n = i = (ti ), i = 0, 1, 2, . . . . (11.17)
For consistency, we should assume that the initial and boundary conditions agree at the corners of the domain: f0 = f (0) = u(0, 0) = (0) = 0 , 4/20/07 191 fn = f ( ) = u(0, ) = (0) = 0 .
c 2006 Peter J. Olver
0.2 0.1
0.2 0.1
0.2 0.1
0.2 -0.1 -0.2
0.4
0.6
0.8
1 -0.1 -0.2
0.2
0.4
0.6
0.8
1 -0.1 -0.2
0.2
0.4
0.6
0.8
1
0.2 0.1
0.2 0.1
0.2 0.1
0.2 -0.1 -0.2
0.4
0.6
0.8
1 -0.1 -0.2
0.2
0.4
0.6
0.8
1 -0.1 -0.2
0.2
0.4
0.6
0.8
1
Figure 11.2.
A Solution to the Heat Equation.
The three equations (11.1417) completely prescribe the numerical approximation algorithm for solving the initial-boundary value problem (11.810). Let us rewrite the scheme in a more transparent matrix form. First, let u(i) = ui,1 , ui,2 , . . . , ui,n-1
T
u(ti , x1 ), u(ti , x2 ), . . . , u(ti , xn-1 )
T
(11.18)
be the vector whose entries are the numerical approximations to the solution values at time ti at the interior nodes. We omit the boundary nodes x0 = 0, xn = , since those values are fixed by the boundary conditions (11.9). Then (11.14) assumes the compact vectorial form (11.19) u(i+1) = A u(i) + b(i) , where 1 - 2 A= 1 - 2 1 - 2 .. .. . . .. .. . 1 - 2 , i 0 0 . . = . . 0 i
b(i)
(11.20)
.
The coefficient matrix A is symmetric and tridiagonal. The contributions (11.17) of the boundary nodes appear in the vector b(i) . This numerical method is known as an explicit scheme since each iterate is computed directly without relying on solving an auxiliary equation -- unlike the implicit schemes to be discussed below. Example 11.4. Let us fix the diffusivity = 1 and the bar length the initial temperature profile 0 x 1, - x, 5 1 7 x - 2, x 10 , u(0, x) = f (x) = 5 5 7 1 - x, 10 x 1, 4/20/07 192
c 2006
= 1. Consider
(11.21)
Peter J. Olver
1
1
1
0.5
0.5
0.5
0.2 -0.5
0.4
0.6
0.8
1 -0.5
0.2
0.4
0.6
0.8
1 -0.5
0.2
0.4
0.6
0.8
1
-1
-1
-1
0.2 0.1
0.2 0.1
0.2 0.1
0.2 -0.1 -0.2
0.4
0.6
0.8
1 -0.1 -0.2
0.2
0.4
0.6
0.8
1 -0.1 -0.2
0.2
0.4
0.6
0.8
1
Figure 11.3.
Numerical Solutions for the Heat Equation Based on the Explicit Scheme.
on a bar of length 1, plotted in the first graph in Figure 11.2. The solution is plotted at the successive times t = ., .02, .04, . . . , .1. Observe that the corners in the initial data are immediately smoothed out. As time progresses, the solution decays, at an exponential rate of 2 9.87, to a uniform, zero temperature, which is the equilibrium temperature distribution for the homogeneous Dirichlet boundary conditions. As the solution decays to thermal equilibrium, it also assumes the progressively more symmetric shape of a single sine arc, of exponentially decreasing amplitude. In our numerical solution, we take the spatial step size h = .1. In Figure 11.3 we compare two (slightly) different time step sizes on the same initial data as used in (11.21). The first sequence uses the time step k = h2 = .01 and plots the solution at times t = 0., .02, .04. The solution is already starting to show signs of instability, and indeed soon thereafter becomes completely wild. The second sequence takes k = .005 and plots the solution at times t = 0., .025, .05. (Note that the two sequences of plots have different vertical scales.) Even though we are employing a rather coarse mesh, the numerical solution is not too far away from the true solution to the initial value problem, which can be found in Figure 11.2. In light of this calculation, we need to understand why our scheme sometimes gives reasonable answers but at other times utterly fails. To this end, let us specialize to homogeneous boundary conditions u(t, 0) = 0 = u(t, ), whereby i = i = 0 for all i = 0, 1, 2, 3, . . . , (11.22) (11.23)
and so (11.19) reduces to a homogeneous, linear iterative system u(i+1) = A u(i) .
According to Proposition 7.8, all solutions will converge to zero, u(i) 0 -- as they are supposed to (why?) -- if and only if A is a convergent matrix. But convergence depends upon the step sizes. 11.4 Example is indicating that for mesh size h = .1, the time step 4/20/07 193
c 2006 Peter J. Olver
k = .01 yields a non-convergent matrix, while k = .005 leads to a convergent matrix and a valid numerical scheme. As we learned in Theorem 7.11, the convergence property of a matrix is fixed by its spectral radius, i.e., its largest eigenvalue in magnitude. There is, in fact, an explicit formula for the eigenvalues of the particular tridiagonal matrix in our numerical scheme, which follows from the following general result. Lemma 11.5. The eigenvalues of an (n - 1) (n - 1) tridiagonal matrix all of whose diagonal entries are equal to a and all of whose sub- and super-diagonal entries are equal to b are k (11.24) , k = 1, . . . , n - 1. k = a + 2 b cos n Proof : The corresponding eigenvectors are vk = k 2k sin , sin , n n ... nk sin n
T
.
Indeed, the j th entry of the eigenvalue equation A vk = k vk reads a sin jk + b n sin (j - 1) k (j + 1) k + sin n n = a + 2 b cos k n sin jk , n
which follows from the trigonometric identity sin + sin = 2 cos + - sin . 2 2 Q.E .D.
In our particular case, a = 1 - 2 and b = , and hence the eigenvalues of the matrix A given by (11.20) are k , k = 1, . . . , n - 1. n Since the cosine term ranges between -1 and +1, the eigenvalues satisfy k = 1 - 2 + 2 cos 1 - 4 < k < 1. Thus, assuming that 0 < 1 guarantees that all | k | < 1, and hence A is a convergent 2 matrix. In this way, we have deduced the basic stability criterion = 1 k , 2 h 2 or k h2 . 2 (11.25)
With some additional analytical work, [28], it can be shown that this is sufficient to conclude that the numerical scheme (11.1417) converges to the true solution to the initialboundary value problem for the heat equation. Since not all choices of space and time steps lead to a convergent scheme, the numerical method is called conditionally stable. The convergence criterion (11.25) places a severe restriction on the time step size. For instance, if we have h = .01, and = 1, then we can only use a time step size k .00005, which is minuscule. It would take an inordinately 4/20/07 194
c 2006 Peter J. Olver
large number of time steps to compute the value of the solution at even a moderate times, e.g., t = 1. Moreover, owing to the limited accuracy of computers, the propagation of round-off errors might then cause a significant reduction in the overall accuracy of the final solution values. An unconditionally stable method -- one that does not restrict the time step -- can be constructed by using the backwards difference formula u(ti , xj ) - u(ti-1 , xj ) u + O(hk ) (ti , xj ) t k (11.26)
to approximate the temporal derivative. Substituting (11.26) and the same approximation (11.12) for uxx into the heat equation, and then replacing i by i + 1, leads to the iterative system i = 0, 1, 2, . . . , ui+1,j - ui+1,j+1 - 2 ui+1,j + ui+1,j-1 = ui,j , (11.27) j = 1, . . . , n - 1, where the parameter = k/h2 is as above. The initial and boundary conditions also have the same form (11.16), (11.17). The system can be written in the matrix form A u(i+1) = u(i) + b(i+1) , (11.28)
where A is obtained from the matrix A in (11.20) by replacing by - . This defines an implicit method since we have to solve a tridiagonal linear system at each step in order to compute the next iterate u(i+1) . However, as we learned in Section 4.5, tridiagonal systems can be solved very rapidly, and so speed does not become a significant issue in the practical implementation of this implicit scheme. Let us look at the convergence properties of the implicit scheme. For homogeneous Dirichlet boundary conditions (11.22), the system takes the form u(i+1) = A-1 u(i) , and the convergence is now governed by the eigenvalues of A-1 . Lemma 11.5 tells us that the eigenvalues of A are k = 1 + 2 - 2 cos k , n k = 1, . . . , n - 1.
As a result, its inverse A-1 has eigenvalues 1 = k 1 1 + 2 k 1 - cos n , k = 1, . . . , n - 1.
Since > 0, the latter are always less than 1 in absolute value, and so A is always a convergent matrix. The implicit scheme (11.28) is convergent for any choice of step sizes h, k, and hence unconditionally stable. 4/20/07 195
c 2006 Peter J. Olver
0.2 0.1 0.2 -0.1 -0.2 0.4 0.6 0.8 1
0.2 0.1 0.2 -0.1 -0.2 0.4 0.6 0.8 1
0.2 0.1 0.2 -0.1 -0.2 0.4 0.6 0.8 1
0.2 0.1 0.2 -0.1 -0.2 0.4 0.6 0.8 1
0.2 0.1 0.2 -0.1 -0.2 0.4 0.6 0.8 1
0.2 0.1 0.2 -0.1 -0.2 0.4 0.6 0.8 1
Figure 11.4.
Numerical Solutions for the Heat Equation Based on the Implicit Scheme.
Example 11.6. Consider the same initial-boundary value problem considered in Example 11.4. In Figure 11.4, we plot the numerical solutions obtained using the implicit scheme. The initial data is not displayed, but we graph the numerical solutions at times t = .2, .4, .6 with a mesh size of h = .1. On the top line, we use a time step of k = .01, while on the bottom k = .005. Unlike the explicit scheme, there is very little difference between the two -- both come much closer to the actual solution than the explicit scheme. Indeed, even significantly larger time steps give reasonable numerical approximations to the solution. Another popular numerical scheme is the CrankNicolson method ui+1,j - ui,j = ui+1,j+1 - 2 ui+1,j + ui+1,j-1 + ui,j+1 - 2 ui,j + ui,j-1 . 2 (11.29)
which can be obtained by averaging the explicit and implicit schemes (11.14, 27). We can write the iterative system in matrix form B u(i+1) = C u(i) + where -1 2 1+ -1 2 .. 1 . -2 .. . . . , . .. .
1 2
b(i) + b(i+1) ,
1 2 1 2
1+ -1 2 B=
1- 1 2 C=
1-
1 2
. .
.. ..
. . . . .. .
(11.30)
Convergence is governed by the generalized eigenvalues of the tridiagonal matrix pair B, C, or, equivalently, the eigenvalues of the product B -1 C, which are 1- k = 1+ 1 - cos k n k 1 - cos n 196
,
k = 1, . . . , n - 1.
(11.31)
4/20/07
c 2006
Peter J. Olver
0.2 0.1 0.2 -0.1 -0.2 0.4 0.6 0.8 1
0.2 0.1 0.2 -0.1 -0.2 0.4 0.6 0.8 1
0.2 0.1 0.2 -0.1 -0.2 0.4 0.6 0.8 1
0.2 0.1 0.2 -0.1 -0.2 0.4 0.6 0.8 1
0.2 0.1 0.2 -0.1 -0.2 0.4 0.6 0.8 1
0.2 0.1 0.2 -0.1 -0.2 0.4 0.6 0.8 1
Figure 11.5.
Numerical Solutions for the Heat Equation Based on the CrankNicolson Scheme.
Since > 0, all of the eigenvalues are strictly less than 1 in absolute value, and so the CrankNicolson scheme is also unconditionally stable. A detailed analysis will show that the errors are of the order of k 2 and h2 , and so it is reasonable to choose the time step to have the same order of magnitude as the space step, k h. This gives the CrankNicolson scheme one advantage over the previous two methods. However, applying it to the initial value problem considered earlier points out a significant weakness. Figure 11.5 shows the result of running the scheme on the initial data (11.21). The top row has space and time step sizes h = k = .1, and does a rather poor job replicating the solution. The second row uses h = k = .01, and performs better except near the corners where an annoying and incorrect local time oscillation persists as the solution decays. Indeed, since most of its eigenvalues are near -1, the CrankNicolson scheme does not do a good job of damping out the high frequency modes that arise from small scale features, including discontinuities and corners in the initial data. On the other hand, most of the eigenvalues of the fully implicit scheme are near zero, and it tends to handle the high frequency modes better, losing out to CrankNicolson when the data is smooth. Thus, a good strategy is to first evolve using the implicit scheme until the small scale noise is dissipated away, and then switch to CrankNicolson to use a much larger time step for final the large scale changes.
11.3. Numerical Solution Methods for the Wave Equation.
Let us now look at some numerical solution techniques for the wave equation. Although this is in a sense unnecessary, owing to the explicit d'Alembert solution formula, the experience we gain in designing workable schemes will serve us well in more complicated situations, including inhomogeneous media, and higher dimensional problems, when analytic solution formulas are no longer available. Consider the wave equation 2u 2u = c2 , t2 x2 0<x< , t 0, (11.32)
modeling vibrations of a homogeneous bar of length with constant wave speed c > 0. We 4/20/07 197
c 2006 Peter J. Olver
impose Dirichlet boundary conditions u(t, 0) = (t), and initial conditions u (0, x) = g(x), t We adopt the same uniformly spaced mesh u(0, x) = f (x), ti = i k, xj = j h, where 0x . (11.34) u(t, ) = (t), t 0. (11.33)
. n In order to discretize the wave equation, we replace the second order derivatives by their standard finite difference approximations (11.6), namely u(ti+1 , xj ) - 2 u(ti , xj ) + u(ti-1 , xj ) 2u (ti , xj ) + O(h2 ), 2 2 t k (11.35) u(ti , xj+1 ) - 2 u(ti , xj ) + u(ti , xj-1 ) 2u 2 (t , x ) + O(k ), x2 i j h2 Since the errors are of orders of k 2 and h2 , we anticipate to be able to choose the space and time step sizes of comparable magnitude: k h. Substituting the finite difference formulae (11.35) into the partial differential equation (11.32), and rearranging terms, we are led to the iterative system ui+1,j = 2 ui,j+1 + 2 (1 - 2 ) ui,j + 2 ui,j-1 - ui-1,j , i = 1, 2, . . . , j = 1, . . . , n - 1, (11.36)
h=
for the numerical approximations ui,j u(ti , xj ) to the solution values at the mesh points. The positive parameter ck > 0 (11.37) = h depends upon the wave speed and the ratio of space and time step sizes. The boundary conditions (11.33) require that ui,0 = i = (ti ), ui,n = i = (ti ), u(i+1) = B u(i) - u(i-1) + b(i) , where 2 2 (1 - 2 ) 2 2 (1 - 2 ) 2 .. .. B= . . 2 .. .. . .
2
i = 0, 1, 2, . . . .
(11.38)
This allows us to rewrite the system in matrix form (11.39)
, u(j)
2 2 (1 - 2 ) 198
2 u1,j j u2,j 0 . . = . , b(j) = . . . . 0 u n-2,j 2 j un-1,j (11.40)
c 2006 Peter J. Olver
4/20/07
1 0.75 0.5 0.25 0.2 -0.25 -0.5 -0.75 -1 0.4 0.6 0.8 1
1 0.75 0.5 0.25 0.2 -0.25 -0.5 -0.75 -1 0.4 0.6 0.8 1
1 0.75 0.5 0.25 0.2 -0.25 -0.5 -0.75 -1 0.4 0.6 0.8 1
1 0.75 0.5 0.25 0.2 -0.25 -0.5 -0.75 -1 0.4 0.6 0.8 1
1 0.75 0.5 0.25 0.2 -0.25 -0.5 -0.75 -1 0.4 0.6 0.8 1
1 0.75 0.5 0.25 0.2 -0.25 -0.5 -0.75 -1 0.4 0.6 0.8 1
Figure 11.6.
Numerically Stable Waves.
The entries of u(i) are, as in (11.18), the numerical approximations to the solution values at the interior nodes. Note that the system (11.39) is a second order iterative scheme, since computing the next iterate u(i+1) requires the value of the preceding two, u(i) and u(i-1) . The one difficulty is getting the method started. We know u(0) since u0,j = fj = f (xj ) is determined by the initial position. However, we also need to find u(1) with entries u1,j u(k, xj ) at time t1 = k in order launch the iteration, but the initial velocity ut (0, x) = g(x) prescribes the derivatives ut (0, xj ) = gj = g(xj ) at time t0 = 0 instead. One way to resolve this difficult would be to utilize the finite difference approximation gj = u1,j - gj u(k, xj ) - u(0, xj ) u (0, xj ) t k k u1,j = fj + k gj . However, the approximation (11.41) is only accurate to order k, whereas the rest of the scheme has error proportional to k 2 . Therefore, we would introduce an unacceptably large error at the initial step. To construct an initial approximation to u(1) with error on the order of k 2 , we need to analyze the local error in more detail. Note that, by Taylor's theorem, u(k, xj ) - u(0, xj ) k 2u u c2 k 2 u u (0, xj ) = (0, xj ) , (0, xj ) + (0, xj ) + k t 2 t2 t 2 x2 where the error is now of order k 2 , and, in the final equality, we have used the fact that u is a solution to the wave equation. Therefore, we find u(k, xj ) u(0, xj ) + k u c2 k 2 2 u (0, xj ) + (0, xj ) t 2 x2 c2 k 2 c2 k 2 f (xj ) fj + k gj + = f (xj ) + k g(xj ) + (f - 2 fj + fj-1 ) , 2 2 h2 j+1 199
c 2006 Peter J. Olver
(11.41)
to compute the required values
4/20/07
1 0.75 0.5 0.25 0.2 -0.25 -0.5 -0.75 -1 0.4 0.6 0.8 1
1 0.75 0.5 0.25 0.2 -0.25 -0.5 -0.75 -1 0.4 0.6 0.8 1
1 0.75 0.5 0.25 0.2 -0.25 -0.5 -0.75 -1 0.4 0.6 0.8 1
1 0.75 0.5 0.25 0.2 -0.25 -0.5 -0.75 -1 0.4 0.6 0.8 1
1 0.75 0.5 0.25 0.2 -0.25 -0.5 -0.75 -1 0.4 0.6 0.8 1
1 0.75 0.5 0.25 0.2 -0.25 -0.5 -0.75 -1 0.4 0.6 0.8 1
Figure 11.7.
Numerically Unstable Waves.
where we can use the finite difference approximation (11.6) for the second derivative of f (x) if no explicit formula is known. Therefore, when we initiate the scheme by setting u1,j = or, in matrix form, u(0) = f , u(1) =
1 2 1 2
2 fj+1 + (1 - 2 )fj +
1 2
2 fj-1 + k gj ,
(11.42)
B u(0) + k g + 1 b(0) , 2
(11.43)
we will have maintained the desired order k 2 (and h2 ) accuracy. Example 11.7. Consider the particular initial value problem utt = uxx , u(0, x) = e- 400 (x-.3) ,
2
ut (0, x) = 0,
0 x 1, t 0,
u(t, 0) = u(1, 0) = 0,
subject to homogeneous Dirichlet boundary conditions on the interval [ 0, 1 ]. The initial data is a fairly concentrated single hump centered at x = .3, and we expect it to split into two half sized humps, which then collide with the ends. Let us choose a space discretization 1 consisting of 90 equally spaced points, and so h = 90 = .0111 . . . . If we choose a time step of k = .01, whereby = .9, then we get reasonably accurate solution over a fairly long time range, as plotted in Figure 11.6 at times t = 0, .1, .2, . . . , .5. On the other hand, if we double the time step, setting k = .02, so = 1.8, then, as plotted in Figure 11.7 at times t = 0, .05, .1, .14, .16, .18, we observe an instability eventually creeping into the picture that eventually overwhelms the numerical solution. Thus, the numerical scheme appears to only be conditionally stable. The stability analysis of this numerical scheme proceeds as follows. We first need to recast the second order iterative system (11.39) into a first order system. In analogy with u(i) R 2n-2 . Example 7.4, this is accomplished by introducing the vector z(i) = u(i-1) Then B -I z(i+1) = C z(i) + c(i) , where C= . (11.44) I O 4/20/07 200
c 2006 Peter J. Olver
Therefore, the stability of the method will be determined by the eigenvalues of the coeffiu cient matrix C. The eigenvector equation C z = z, where z = , can be written out v in its individual components: B u - v = u, ( B - 2 - 1) v = 0, u = v.
Substituting the second equation into the first, we find 1 v. The latter equation implies that v is an eigenvector of B with + -1 the corresponding eigenvalue. The eigenvalues of the tridiagonal matrix B are governed by Lemma 11.5, in which a = 2(1 - 2 ) and b = 2 , and hence are or Bv = + + 1 =2 1 - 2 + 2 cos k n , k = 1, . . . , n - 1.
Multiplying both sides by leads to a quadratic equation for the eigenvalues, 2 - 2 ak + 1 = 0, where 1 - 2 2 < ak = 1 - 2 + 2 cos k < 1. n (11.45)
Each pair of solutions to these n - 1 quadratic equations, namely = ak k a2 - 1 , k (11.46)
yields two eigenvalues of the matrix C. If ak > 1, then one of the two eigenvalues will be larger than one in magnitude, which means that the linear iterative system has an exponentially growing mode, and so u(i) as i for almost all choices of initial data. This is clearly incompatible with the wave equation solution that we are trying to approximate, which is periodic and hence remains bounded. On the other hand, if | ak | < 1, then the eigenvalues (11.46) are complex numbers of modulus 1, indicated stability (but not convergence) of the matrix C. Therefore, in view of (11.45), we should require that h ck < 1, or k< , (11.47) = h c which places a restriction on the relative sizes of the time and space steps. We conclude that the numerical scheme is conditionally stable. The stability criterion (11.47) is known as the Courant condition, and can be assigned a simple geometric interpretation. Recall that the wave speed c is the slope of the characteristic lines for the wave equation. The Courant condition requires that the mesh slope, which is defined to be the ratio of the space step size to the time step size, namely h/k, must be strictly greater than the characteristic slope c. A signal starting at a mesh point (ti , xj ) will reach positions xj k/c at the next time ti+1 = ti + k, which are still between the mesh points xj-1 and xj+1 . Thus, characteristic lines that start at a mesh point are not allowed to reach beyond the neighboring mesh points at the next time step. For instance, in Figure 11.8, the wave speed is c = 1.25. The first figure has equal mesh spacing k = h, and does not satisfy the Courant condition (11.47), whereas the 4/20/07 201
c 2006 Peter J. Olver
Figure 11.8.
The Courant Condition.
second figure has k = 1 h, which does. Note how the characteristic lines starting at a 2 given mesh point have progressed beyond the neighboring mesh points after one time step in the first case, but not in the second.
4/20/07
202
c 2006
Peter J. Olver