14 Chapter Vibration and Diffusion in OneDimensional Media
In this chapter, we study the solutions, both analytical and numerical, to the two most important equations of one-dimensional continuum dynamics. The heat equation models the diffusion of thermal energy in a body; here, we analyze 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. Again, we restrict our attention here to the case of waves in a one-dimensional medium, e.g., a string, or a bar, or a column of air. The two- and three-dimensional versions of these fundamental equations will be analyzed in the later Chapters 17 and 18. As we saw in Section 12.1, the basic solution strategy is inspired by our eigenvaluebased methods for solving linear systems of ordinary differential equations. Substituting the appropriate exponential or trigonometric ansatz will effectively reduce the partial differential equation to a one-dimensional boundary value problem. The linear superposition principle implies that general solution can then be expressed as a infinite series in the resulting eigenfunction solutions. In both cases considered here, the eigenfunctions of the one-dimensional boundary value problem are trigonometric, and so the solution to the partial differential equation takes the form of a time-dependent Fourier series. Although we cannot, in general, analytically sum the Fourier series to produce a simpler formula for the solution, there are a number of useful observations that can be gleaned from it. In the case of the heat equation, the solutions decay exponentially fast to thermal equilibrium, at a rate governed by the smallest positive eigenvalue of the associated boundary value problem. The higher order Fourier modes damp out very rapidly, and so the heat equation can be used to automatically smooth and denoise signals and images. It also implies that the heat equation cannot be run backwards in time -- determining the initial temperature profile from a later measurement is an ill-posed problem. The response to a concentrated unit impulse leads to the fundamental solution, which can then be used to construct integral representations of the solution to the inhomogeneous heat equation. We will also explain how to exploit the symmetry properties of the differential equation in order to construct new solutions from known solutions. In the case of the wave equation, each Fourier mode vibrates with its natural frequency. In a stable situation, the full solution is a linear combination of these fundamental vibrational modes, while an instability induces an extra linearly growing mode. For one-dimensional media, the natural frequencies are integral multiples of a single lowest frequency, and hence the solution is periodic, which, in particular, explains the tonal qualities 2/25/07 752
c 2006 Peter J. Olver
of string and wind instruments. The one-dimensional wave equation admits an alternative explicit solution formula, due to d'Alembert, which points out the role of characteristics in signal propagation and the behavior of solutions. The analytic and series solution methods serve to shed complementary lights on the physical phenomena of waves and vibration in continuous media. Following our analytical study, we introduce several basic numerical solution methods for both the heat and the wave equations. 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 Chapter 10. In preparation for our treatment of partial differential equations in higher dimensions, the final section introduces a general framework for dynamics that incorporates both discrete dynamics, modeled by linear systems of ordinary differential equations, and continuum dynamics, modeled by the heat and wave equations, their multidimensional counterparts, as well as the equilibrium boundary value problems governing bars, beams, membranes and solid bodies. Common features, based on eigenvalues, are discussed.
14.1. The Diffusion and Heat Equations.
Let us begin with a physical derivation of the heat equation from first principles of thermodynamics. The reader solely interested in the mathematical developments can skip ahead to the following subsection. However, physical insight can often play an critical role in understanding the underlying mathematics, and is neglected at one's peril. We consider a bar -- meaning a thin, heat-conducting body of length . "Thin" means that we can regard the bar as a one-dimensional continuum with no significant transverse temperature variation. We use 0 x to denote the position along the bar. Our goal is to find the temperature u(t, x) of the bar at position x and time t. The dynamical equations governing the temperature are based on three fundamental physical laws. The first law is that, in the absence of external sources, thermal energy can only enter the bar through its ends. In physical terms, we are assuming that the bar is fully insulated along its length. Let (t, x) to denote the thermal energy in the bar at position x and time t. Consider a small section of the bar lying between x and x+x. The total amount of heat
x+x
energy contained in this section is obtained by integrating (summing):
x
(t, y) dy.
Further, let w(t, x) denote the heat flux , i.e., the rate of flow of thermal energy along the bar. We use the convention that w(t, x) > 0 means that the energy is moving to the right, while w(t, x) < 0 if it moves to the left. The first law implies that the rate of change in the thermal energy in any section of the bar is equal to the total heat flux, namely the amount of the heat passing through its ends. Therefore, in view of our sign convention on 2/25/07 753
c 2006 Peter J. Olver
the flux, t
x+x x
(t, y) dy = - w(t, x + x) + w(t, x),
the latter two terms denoting the respective flux of heat into the section of the bar at its right and left ends. Assuming sufficient regularity of the integrand, we are permitted to bring the derivative inside the integral. Thus, dividing both sides of the resulting equation by x, x+x w(t, x + x) - w(t, x) 1 (t, y) dy = - . x x t x In the limit as the length x 0, the right hand side of this equation converges to minus the x derivative of w(t, x), while, by the Fundamental Theorem of Calculus, the left hand side converges to the integrand /t at the point x; the net result is the fundamental differential equation w = - (14.1) t x relating thermal energy and heat flux w. A partial differential equation of this particular form is known as a conservation law , and, in this instance, formulates the law of conservation of thermal energy. See Exercise for details. The second physical law is a constitutive assumption, based on experimental evidence. In most physical materials, thermal energy is found to be proportional to temperature, (t, x) = (x) u(t, x). The factor (x) = (x) (x) > 0 is the product of the density of the material and its specific heat , which is the amount of heat energy required to raise the temperature of a unit mass of the material by one unit. Note that we are assuming the bar is not changing in time, and so physical quantities such as density and specific heat depend only on position x. We also assume, perhaps with less physical justification, that the material properties do not depend upon the temperature; otherwise, we would be led to a much more difficult nonlinear diffusion equation. The third physical law relates the heat flux to the temperature. Physical experiments in a wide variety of materials indicate that the heat energy moves from hot to cold at a rate that is in direct proportion to the rate of change -- meaning the derivative -- of the temperature. The resulting linear constitutive relation w(t, x) = - (x) u x (14.3) (14.2)
is known as Fourier's Law of Cooling. The proportionality factor (x) > 0 is called the thermal conductivity of the bar at position x. A good heat conductor, e.g., silver, will have high conductivity, while a poor conductor, e.g., glass, will have low conductivity. u (t, x) > 0 the The minus sign tells us that heat energy moves from hot to cold; if x 2/25/07 754
c 2006 Peter J. Olver
temperature is increasing from left to right, and so the heat energy moves back to the left, with consequent flux w(t, x) < 0. Combining the three laws (14.1), (14.2) and (14.3) produces the basic partial differential equation u (x) , 0 < x < , (14.4) (x) u = t x x governing the diffusion of heat in a non-uniform bar. The resulting linear diffusion equation is used to model a variety of diffusive processes, including heat flow, chemical diffusion, population dispersion, and the spread of infectious diseases. If, in addition, we allow external heat sources h(t, x), then the linear diffusion equation acquires an inhomogeneous term: u (x) + h(t, x), 0 < x < . (14.5) (x) u = t x x In order to uniquely prescribe the solution u(t, x) to the diffusion equation, we need to specify the initial temperature distribution u(t0 , x) = f (x), 0 x , (14.6)
along the bar at an initial time t0 . In addition, we must impose suitable boundary conditions at the two ends of the bar. As with the equilibrium equations discussed in Chapter 11, there are three common physical types. The first is a Dirichlet boundary condition, where an end of the bar is held at prescribed temperature. Thus, the boundary condition u(t, 0) = (t) (14.7)
fixes the temperature at the left hand end of the bar. Alternatively, the Neumann boundary condition u (14.8) (t, 0) = (t) x u (t, 0) at the left hand end. In particular, the prescribes the heat flux w(t, 0) = - (0) x homogeneous Neumann condition with (t) 0 corresponds to an insulated end, where no heat can flow in or out. Each end of the bar should have one or the other of these boundary conditions. For example, a bar with both ends having prescribed temperatures is governed by the pair of Dirichlet boundary conditions u(t, 0) = (t), u(t, ) = (t), (14.9)
whereas a bar with two insulated ends requires two homogeneous Neumann boundary conditions u u (14.10) (t, 0) = 0, (t, ) = 0. x x The mixed case, with one end fixed and the other insulated, is similarly formulated. Finally, the periodic boundary conditions u(t, 0) = u(t, ), 2/25/07 755 u u (t, 0) = (t, ), x x
c 2006
(14.11)
Peter J. Olver
correspond to a circular ring of length . As before, we are assuming the heat is only allowed to flow around the ring -- insulation prevents any radiation of heat from one side of the ring to the other. The Heat Equation In this book, we will retain the term "heat equation" to refer to the homogeneous case, in which the bar is made of a uniform material, and so its density , conductivity , and specific heat are all positive constants. Under these assumptions, the homogeneous diffusion equation (14.4) reduces to the heat equation 2u u = t x2 for the temperature u(t, x) in the bar at time t and position x. The constant = = (14.13) (14.12)
is called the thermal diffusivity of the bar, and incorporates all of its relevant physical properties. The solution u(t, x) will be uniquely prescribed once we specify initial conditions (14.6) and a suitable pair of boundary conditions at the ends of the bar. As we learned in Section 12.1, the elementary, separable solutions to the heat equation are based on the exponential ansatz u(t, x) = e- t v(x), (14.14)
where v(x) is a time-independent function. Substituting the solution formula (14.14) into (14.12) and canceling the common exponential factors, we find that v(x) must solve the ordinary differential equation - v = v. In other words, v is an eigenfunction with eigenvalue , for the second derivative operator K = - D2 . Once we determine the eigenvalues and eigenfunctions, we will be able to reconstruct the solution u(t, x) as a linear combination, or, rather, infinite series in the corresponding separable eigenfunction solutions. Let us consider the simplest case of a uniform bar held at zero temperature at each end. For simplicity, we take the initial time to be t0 = 0, and so the initial and boundary conditions are u(t, 0) = 0, u(t, ) = 0, t 0, (14.15) u(0, x) = f (x), 0 < x < . According to the general prescription, we need to solve the eigenvalue problem d2 v + v = 0, dx2 v(0) = 0, v() = 0. (14.16)
As noted in Exercise , positive definiteness of the underlying differential operator K = - D2 when subject to Dirichlet boundary conditions implies that we need only look for positive eigenvalues: > 0. In Exercises , the skeptical reader is asked to check 2/25/07 756
c 2006 Peter J. Olver
explicitly that if 0 or is complex, then the boundary value problem (14.16) admits only the trivial solution v(x) 0. Setting = 2 with > 0, the general solution to the differential equation is a trigonometric function v(x) = a cos x + b sin x, where a, b are constants whose values will be fixed by the boundary conditions. The boundary condition at x = 0 requires a = 0. The second boundary condition requires v() = b sin = 0. Therefore, assuming b = 0, as otherwise the solution is trivial, must be an integer multiple of , and so 2 3 = , , , ... . We conclude that the eigenvalues and eigenfunctions of the boundary value problem (14.16) are n x n 2 (14.17) , n = 1, 2, 3, . . . . , vn (x) = sin n = The corresponding separable solutions (14.14) to the heat equation with the given boundary conditions are un (t, x) = exp - n2 2 t 2 sin n x , n = 1, 2, 3, . . . . (14.18)
Each represents a trigonometrically oscillating temperature profile that maintains its form while decaying to zero at an exponential rate. The first of these, u1 (t, x) = exp - 2 t 2 sin x ,
experiences the slowest decay. The higher "frequency" modes un (t, x), n 2, all go to zero at a faster rate, with those having a highly oscillatory temperature profile, where n 0, effectively disappearing almost instantaneously. Thus, small scale temperature fluctuations tend to rapidly cancel out through diffusion of heat energy. Linear superposition is used to assemble the general series solution
u(t, x) =
n=1
bn un (t, x) =
n=1
bn exp
-
n2 2 t 2
sin
n x
(14.19)
as a combination of the separable solutions. Assuming that the series converges, the initial temperature profile is n x = f (x). (14.20) u(0, x) = bn sin n=1 This has the form of a Fourier sine series (12.44) on the interval [ 0, ]. By orthogonality of the eigenfunctions -- which is a direct consequence of the self-adjointness of the underlying 2/25/07 757
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 14.1.
A Solution to the Heat Equation.
boundary value problem (14.16) -- the coefficients are determined by the inner product formulae (12.45), and so bn = 2
f (x) sin
0
n x dx,
n = 1, 2, 3, . . . .
(14.21)
The resulting solution (14.19) describes the Fourier sine series for the temperature u(t, x) of the bar at each later time t 0. It can be rigorously proved that, for quite general initial conditions, the Fourier series does indeed converge to a solution to the initial-boundary value problem, [183]. Example 14.1. Consider the initial temperature profile 1 0 x 5, - x, 1 7 u(0, x) = f (x) = x - 2, 5 5 x 10 , 7 1 - x, 10 x 1, b1 = .0448 . . . , b5 = - .0081 . . . ,
(14.22)
on a bar of length 1, plotted in the first graph in Figure 14.1. Using (14.21), the first few Fourier coefficients of f (x) are computed as b6 = .0066 . . . ,
b2 = - .096 . . . ,
b3 = - .0145 . . . , b7 = .0052 . . . ,
b4 = 0, b8 = 0,
... .
Setting = 1, the resulting Fourier series solution to the heat equation is u(t, x) =
n=1
bn un (t, x) =
n=1
2
b n e- n
2
2 t
sin n x sin 2 x - .0145 e- 9
2
= .0448 e-
t
In Figure 14.1, 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, 2/25/07 758
c 2006 Peter J. Olver
sin x - .096 e- 4
2
t
t
sin 3 x - .
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. Smoothing and Long Time Behavior The fact that we can write the solution to an initial-boundary value problem in the form of an infinite series is progress of a sort. However, because it cannot be summed in closed form, this "solution" is much less satisfying than a direct, explicit formula. Nevertheless, there are important qualitative and quantitative features of the solution that can be easily gleaned from such series expansions. If the initial data f (x) is piecewise continuous, then its Fourier coefficients are uniformly bounded; indeed, for any n 1, | bn | 2
f (x) sin
0
2 n x dx
0
| f (x) | dx M.
(14.23)
This property holds even for quite irregular data; for instance, the Fourier coefficients (12.59) of the delta function are also uniformly bounded. Under these conditions, each term in the series solution (14.19) is bounded by an exponentially decaying function bn exp - n2 2 t 2 sin n x M exp - n2 2 t 2 .
This means that, as soon as t > 0, most of the high frequency terms, n 0, will be extremely small. Only the first few terms will be at all noticeable, and so the solution essentially degenerates into a finite sum over the first few Fourier modes. As time increases, more and more of the Fourier modes will become negligible, and the sum further degenerates into progressively fewer significant terms. Eventually, as t , all of the Fourier modes will decay to zero. Therefore, the solution will converge exponentially fast to a zero temperature profile: u(t, x) 0 as t , representing the bar in its final uniform thermal equilibrium. The fact that its equilibrium temperature is zero is the result of holding both ends of the bar fixed at zero temperature, and any initial heat energy will eventually be dissipated away through the ends. The last term to disappear is the one with the slowest decay, namely u(t, x) b1 exp - 2 t 2 sin x , where b1 = 1
f (x) sin x dx.
0
(14.24)
Generically, b1 = 0, and the solution approaches thermal equilibrium exponentially fast with rate equal to the smallest eigenvalue, 1 = 2 /2 , which is proportional to the thermal diffusivity divided by the square of the length of the bar. The longer the bar, or the smaller the diffusivity, the longer it takes for the effect of holding the ends at zero temperature to propagate along the entire bar. Also, again provided b1 = 0, the asymptotic shape of the temperature profile is a small sine arc, just as we observed in Example 14.1. In exceptional situations, namely when b1 = 0, the solution decays even faster, at a rate 2/25/07 759
c 2006 Peter J. Olver
7 6 5 4 3 2 1 0.2 7 6 5 4 3 2 1 0.2 0.4 0.6 0.8 1 0.4 0.6 0.8 1
7 6 5 4 3 2 1 0.2 7 6 5 4 3 2 1 0.2 0.4 0.6 0.8 1 0.4 0.6 0.8 1
7 6 5 4 3 2 1 0.2 7 6 5 4 3 2 1 0.2 0.4 0.6 0.8 1 0.4 0.6 0.8 1
Figure 14.2.
Denoising a Signal with the Heat Equation.
equal to the eigenvalue k = k 2 2 /2 corresponding to the first nonzero term, bk = 0, in the series; its asymptotic shape now oscillates k times over the interval. The heat equation's smoothing effect on irregular initial data by fast damping of the high frequency modes underlies its effectiveness for smoothing out and denoising signals. We take the initial data u(0, x) = f (x) to be a noisy signal, and then evolve the heat equation forward to a prescribed time t > 0. The resulting function g(x) = u(t , x) will be a smoothed version of the original signal f (x) in which most of the high frequency noise has been eliminated. Of course, if we run the heat flow for too long, all of the low frequency features will be also be smoothed out and the result will be a uniform, constant signal. Thus, the choice of stopping time t is crucial to the success of this method. Figure 14.2 shows the effect running the heat equation , with = 1, to times t = 0., .00001, .00005, .0001, .001, .01 on the same signal from Figure 13.5. Observe how quickly the noise is removed. By the final time, the overall smoothing effect of the heat flow has caused significant degradation (blurring) of the original signal. The heat equation approach to denoising has the advantage that no Fourier coefficients need be explicitly computed, nor does one need to reconstruct the smoothed signal from its remaining Fourier coefficients. The final section discusses some numerical methods that can be used to solve the heat equation directly. Another, closely related observation is that, for any fixed time t > 0 after the initial moment, the coefficients in the Fourier series (14.19) decay exponentially fast as n . According to the discussion at the end of Section 12.3, this implies that the solution u(t, x) is a very smooth, infinitely differentiable function of x at each positive time t, no matter how unsmooth the initial temperature profile. We have discovered the basic smoothing property of heat flow. Theorem 14.2. If u(t, x) is a solution to the heat equation with piecewise continuous
To be honest, we are using periodic boundary conditions in the figures, although the Dirichlet version leads to similar results
2/25/07
760
c 2006
Peter J. Olver
initial data f (x) = u(0, x), or, more generally, initial data satisfying (14.23), then, for any t > 0, the solution u(t, x) is an infinitely differentiable function of x. After even a very short amount of time, the heat equation smoothes out most, and, eventually, all of the fluctuations in the initial temperature profile. As a consequence, it becomes impossible to reconstruct the initial temperature u(0, x) = f (x) by measuring the temperature distribution h(x) = u(t, x) at a later time t > 0. Diffusion is irreversible -- we cannot run the heat equation backwards in time! Indeed, if the initial data u(0, x) = f (x) is not smooth, there is no function u(t, x) for t < 0 that could possibly yield such an initial distribution because all corners and singularities are smoothed out by the diffusion process as t goes forward! Or, to put it another way, the Fourier coefficients (14.21) of any purported solution will be exponentially growing when t < 0, and so high frequency noise will completely overwhelm the solution. For this reason, the backwards heat equation is said to be ill-posed . On the other hand, the unsmoothing effect of the backwards heat equation does have potential benefits. For example, in image processing, diffusion will gradually blur an image. Image enhancement is the reverse process, and requires running the heat flow backwards in some stable manner. One option is to restrict to the backwards evolution to the first few Fourier modes, which prevents the small scale fluctuations from overwhelming the computation. Similar issues arise in the reconstruction of subterranean profiles from seismic data, a problem of great concern in the oil and gas industry. In forensics, determining the time of death based on the current temperature of a corpse also requires running the equations governing the dissipation of body heat backwards in time. For these and other applications, a key issue in contemporary research is how to cleverly circumventing the ill-posedness of the backwards heat flow. Remark : The irreversibility of the heat equation points out a crucial distinction between partial differential equations and ordinary differential equations. Ordinary differential equations are always reversible -- unlike the heat equation, existence, uniqueness and continuous dependence properties of solutions are all equally valid in reverse time (although the detailed qualitative and quantitative properties of solutions can very well depend upon whether time is running forwards or backwards). The irreversibility of partial differential equations modeling the diffusive processes in our universe may well be why, in our experience, Time's Arrow points exclusively to the future. The Heated Ring Let us next consider the periodic boundary value problem modeling heat flow in an insulated circular ring. Let us fix the length of the ring to be = 2 , with - < x < representing "angular" coordinate around the ring. For simplicity, we also choose units in which the thermal diffusivity is = 1. Thus, we seek to solve the heat equation 2u u = , - < x < , t x2 subject to periodic boundary conditions u(t, - ) = u(t, ), 2/25/07 t > 0, (14.25)
u u (t, - ) = (t, ), x x 761
t 0.
c 2006
(14.26)
Peter J. Olver
The initial temperature distribution is u(0, x) = f (x), - < x < . (14.27)
The resulting temperature u(t, x) will be a periodic function in x of period 2 . Substituting the separable solution ansatz u(t, x) = e- t v(x) into the heat equation and the boundary conditions leads to the periodic eigenvalue problem d2 v + v = 0, dx2 v(- ) = v(), v(- ) = v(). (14.28)
As we know, in this case the eigenvalues are n = n2 where n = 0, 1, 2, . . . is a non-negative integer, and the corresponding eigenfunction solutions are the trigonometric functions vn (x) = cos n x, v n (x) = sin n x, n = 0, 1, 2, . . . .
Note that 0 = 0 is a simple eigenvalue, with constant eigenfunction cos 0x = 1 -- the sine solution sin 0x 0 is trivial -- while the positive eigenvalues are, in fact, double, each possessing two linearly independent eigenfunctions. The corresponding separable solutions to the heated ring equation are un (t, x) = e- n
2
t
cos n x,
un (t, x) = e- n
2
t
sin n x,
n = 0, 1, 2, 3, . . . .
The resulting infinite series solution is
u(t, x) =
1 2
a0 +
n=1
an e- n t cos n x + bn e- n
2
2
t
sin n x .
(14.29)
The initial conditions require
u(0, x) =
1 2
a0 +
n=1
an cos n x + bn sin n x = f (x),
(14.30)
which is precisely the Fourier series of the initial temperature profile f (x). Consequently, an = 1
f (x) cos n x dx,
-
bn =
1
f (x) sin n x dx,
-
(14.31)
are the usual Fourier coefficients of f (x). As in the Dirichlet problem, after the initial instant, the high frequency terms in 2 the series (14.29) become extremely small, since e- n t 1 for n 0. Therefore, as soon as t > 0, the solution essentially degenerates into a finite sum over the first few Fourier modes. Moreover, as t , all of the Fourier modes will decay to zero with the exception of the constant one, with null eigenvalue 0 = 0. Therefore, the solution will converge exponentially fast to a constant temperature profile: u(t, x) - 2/25/07 1 1 a0 = 2 2 762
f (x) dx,
- c 2006 Peter J. Olver
which equals the average of the initial temperature profile. Physically, we observe that the heat energy is redistributed so that the ring achieves a uniform constant temperature and is in thermal equilibrium. Indeed, the total heat
H(t) =
-
u(t, x) dx = constant
(14.32)
is conserved, meaning constant, for all time; the proof of this fact is left as an Exercise . (On the other hand, the Dirichlet boundary value problem does not conserve hear energy.) Prior to equilibrium, only the lowest frequency Fourier modes will still be noticeable, and so the solution will asymptotically look like u(t, x) where 1 a1 = r1 cos 1 = 2
1 2
a0 + e- t (a1 cos x + b1 sin x) =
1 2
a0 + r1 e- t cos(x + 1 ),
(14.33)
f (x) cos x dx,
-
1 b1 = r1 sin 1 = 2
f (x) sin x dx.
-
Thus, for most initial data, the solution approaches thermal equilibrium exponentially fast, at a unit rate. The exceptions are when r1 = a2 + b2 = 0, for which the rate of 1 1 2 convergence is even faster, namely at a rate e- k t where k is the smallest integer such that rk = a2 + b2 = 0. k k Inhomogeneous Boundary Conditions So far, we have concentrated our attention on homogeneous boundary conditions. There is a simple trick that will convert a boundary value problem with inhomogeneous but constant Dirichlet boundary conditions, u(t, 0) = , u(t, ) = , t 0, (14.34)
into a homogeneous Dirichlet problem. We begin by solving for the equilibrium temperature profile, which is the affine function that satisfies the two boundary conditions, namely u (x) = + The difference - x. (14.35)
- x (14.36) measures the deviation of the solution from equilibrium. It clearly satisfies the homogeneous boundary conditions at both ends: u(t, x) = u(t, x) - u (x) = u(t, x) - - u(t, 0) = 0 = u(t, ). Moreover, by linearity, since both u(t, x) and u (x) are solutions to the heat equation, so is u(t, x). The initial data must be similarly adapted: u(0, x) = f (x) = f (x) - u (x) = f (x) - - 2/25/07 763 - x.
c 2006 Peter J. Olver
Solving the resulting homogeneous initial value problem, we write u(t, x) in Fourier series form (14.19), where the Fourier coefficients are computed from the modified initial data f (x). The solution to the inhomogeneous boundary value problem thus has the series form - u(t, x) = + x + where 2 bn =
0
bn exp
n=1
n2 2 t - 2
sin
n x ,
(14.37)
f (x) sin
n x dx,
n = 1, 2, 3, . . . .
(14.38)
Since, for any reasonable initial data, u(t, 0) will decay to zero at an exponential rate as t , the actual temperature profile (14.37) will asymptotically decay to the equilibrium profile, - x u(t, x) - u (x) = + at the same exponentially fast rate. This method does not apply when the boundary conditions are time-dependent: u(t, 0) = (t), u(t, ) = (t). Attempting to mimic the preceding technique, we discover that the deviation (t) - (t) u(t, x) = u(t, x) - u (t, x), where u (t, x) = (t) + x, (14.39) does satisfy the homogeneous boundary conditions, but now solves an inhomogeneous version of the heat equation: u 2u u = - h(t, x), where h(t, x) = (t, x). t x2 t Solution techniques in this case will be discussed below. (14.40)
14.2. Symmetry and the Maximum Principle.
So far we have relied on the method of separation of variables to construct explicit solutions to partial differential equations. A second useful solution technique relies on exploiting inherent symmetry properties of the differential equation. Unlike separation of variables , symmetry methods can be also successfully applied to produce solutions to a broad range of nonlinear partial differential equations; examples can be found in Chapter 22. While we do not have the space or required mathematical tools to develop the full apparatus of symmetry techniques, we can introduce the important concept of a similarity solution, applied in the particular context of the heat equation. In general, by a symmetry of an equation, we mean a transformation, either linear (as in Section 7.2), affine (as in Section 7.3), or even nonlinear, that takes solutions to solutions. Thus, if we have a symmetry, and know one solution, then we can construct a
This is not quite fair: separation of variables can be applied to some special nonlinear partial differential equations such as HamiltonJacobi equations, [ 130 ].
2/25/07
764
c 2006
Peter J. Olver
second solution by applying the symmetry. And, possibly, a third solution by applying the symmetry yet again. And so on. If we know lots of symmetries, then we can produce lots and lots of solutions by this simple device. Remark : General symmetry techniques are founded on the theory of Lie groups, named after the influential nineteenth century Norwegian mathematician Sophus Lie (pronounced "Lee"). Lie's theory provides an algorithm for completely determining all the symmetries of a given differential equation, but this is beyond the scope of this introductory text. However, direct inspection and/or physical intuition will often detect the most important symmetries without appealing to such a sophisticated theory. Modern applications of Lie's symmetry methods to partial differential equations arising in physics and engineering can be traced back to the influential book of G. Birkhoff, [18], on hydrodynamics. A complete and comprehensive treatment of symmetry methods can be found in the first author's book [143], and, at a more introductory level, in the recent books by Hydon, [106], and Cantwell, [39], with particular emphasis on fluid mechanics. The heat equation serves as an excellent testing ground for the general symmetry methodology, as it admits a rich variety of symmetry transformations that take solutions to solutions. The simplest are the translations. Moving the space and time coordinates by a fixed amount, t - t - a, x - x - b, (14.41) where a, b are constants, changes the function u(t, x) into the translated function U (t, x) = u(t - a, x - b). (14.42)
A simple application of the chain rule proves that the partial derivatives of U with respect to t and x agree with the corresponding partial derivatives of u, so U u U u 2U 2u = , = , = , t t x x x2 x2 and so on. In particular, the function U (t, x) is a solution to the heat equation Ut = Uxx whenever u(t, x) also solves ut = uxx . Physically, the translation symmetries formalize the property that the heat equation models a homogeneous medium, and hence the solution does not depend on the choice of reference point or origin of our coordinate system. As a consequence, each solution to the heat equation will produce an infinite family of translated solutions. For example, starting with the separable solution u(t, x) = e- t sin x, we immediately produce the additional solutions u(t, x) = e- (t-a) sin (x - b), valid for any choice of constants a, b. Warning: Typically, the symmetries of a differential equation do not respect initial or boundary conditions. For instance, if u(t, x) is defined for t 0 and in the domain 0 x , then its translated version U (t, x) is defined for t a and in the translated domain b x + b, and so will solve a translated initial-boundary value problem. 2/25/07 765
c 2006 Peter J. Olver
A second, even more important class of symmetries are the scaling invariances. We already know that if u(t, x) is a solution, so is any scalar multiple c u(t, x); this is a simple consequence of linearity of the heat equation. We can also add an arbitrary constant to the temperature, noting that U (t, x) = c u(t, x) + k (14.43) is a solution for any choice of constants c, k. Physically, the transformation (14.43) amounts to a change in the scale for measuring temperature. For instance, if u is measured degrees 9 Celsius, and we set c = 5 and k = 32, then U = 9 u + 32 will be measured in degrees 5 Fahrenheit. Thus, reassuringly, the physical processes described by the heat equation do not depend upon our choice of thermometer. More interestingly, suppose we rescale the space and time variables: t - t, x - x, (14.44)
where , > 0 are positive constants. The effect of such a scaling transformation is to convert u(t, x) into a rescaled function U (t, x) = u( t, x). (14.45)
The derivatives of U are related to those of u according to the following formulae, which are direct consequences of the multi-variable chain rule: U u = , t t U u = , x x 2U 2u = 2 2 . x2 x
Therefore, if u satisfies the heat equation ut = uxx , then U satisfies the rescaled heat equation Ut = ut = uxx = 2 Uxx , which we rewrite as Ut = Uxx , where = . 2 (14.46)
Thus, the net effect of scaling space and time is merely to rescale the diffusion coefficient in the heat equation. Physically, the scaling symmetry (14.44) corresponds to a change in the physical units used to measure time and distance. For instance, to change from seconds to minutes, set = 60, and from meters to yards, set = 1.0936. The net effect (14.46) on the diffusion coefficient is a reflection of its physical units, namely distance2 /time. In particular, if we choose = 1 , = 1,
then the rescaled diffusion coefficient becomes = 1. This observation has the following important consequence. If U (t, x) solves the heat equation for a unit diffusivity, = 1, then u(t, x) = U ( t, x) (14.47) 2/25/07 766
c 2006 Peter J. Olver
solves the heat equation for the diffusivity . Thus, the only effect of the diffusion coefficient is to speed up or slow down time! A body with diffusivity = 2 will cool down twice as fast as a body (of the same shape subject to the same boundary conditions and initial conditions) with diffusivity = 1. Note that this particular rescaling has not altered the space coordinates, and so U (t, x) is defined on the same domain as u(t, x). On the other hand, if we set = 2 , then the rescaled diffusion coefficient is exactly the same as the original: = . Thus, the transformation t - 2 t, x - x, (14.48)
does not alter the equation, and hence defines a scaling symmetry, also known as a similarity transformation, for the heat equation. Combining (14.48) with the linear rescaling u c u, we make the elementary, but important observation that if u(t, x) is any solution to the heat equation, then so is the function U (t, x) = c u( 2 t, x), (14.49)
2
for the same diffusion coefficient . For example, rescaling the solution u(t, x) = e- t cos x 2 2 leads to the solution U (t, x) = c e- t cos x. Warning: As in the case of translations, rescaling space by a factor = 1 will alter the domain of definition of the solution. If u(t, x) is defined for 0 x , then U (t, x) is defined for 0 x /. Suppose that we have solved the heat equation for the temperature u(t, x) on a bar of length 1, subject to certain initial and boundary conditions. We are then given a bar composed of the same material of length 2. Since the diffusivity coefficient has not changed, 1 and we can directly construct the new solution U (t, x) by rescaling. Setting = 2 will 1 serve to double the length. If we also rescale time by a factor = 2 = 4 , then the 1 rescaled function U (t, x) = u 1 t, 2 x will be a solution of the heat equation on the longer 4 bar with the same diffusivity constant. The net effect is that the rescaled solution will be evolving four times as slowly as the original. Thus, it effectively takes a bar that is twice the length four times as long to cool down. The Maximum Principle The Maximum Principle is the mathematical formulation of the thermodynamical law that heat cannot, in the absence of external sources, achieve a value at a point that strictly larger than its surroundings. We will help in the proof to formulate it for the more general case when the heat equation is subject to external forcing. Theorem 14.3. Let > 0. Suppose u(t, x) is a solution to the forced heat equation u 2u = + F (t, x), t x2 on the rectangular domain R = { a x b, 0 t c. }
Suppose F (t, x) 0 for all (t, x) R. Then the global maximum of u(t, x) on the domain R occurs either at t = 0 or at x = a or x = b. 2/25/07 767
c 2006 Peter J. Olver
In other words, if we are no;t adding in any external heat, the maximum temperature ion the bar occurs either at the initial time, or at one of the endpoints. Proof : First let us prove the result under the assumption that F (t, x) < 0, and hence u 2u < t x2 (14.50)
everywhere in R. Suppose u(t, x) has a (local) maximum at a point in the interior of R. Then, by multivariable calculus (see also Theorem 19.43), its gradient must vanish, so ut = ux = 0, and its Hessian matrix must be positive definite, which, in particular, requires that uxx 0. But these contradict the inequality (14.50). Thus, the solution cannot have a local interior maximum. If the maximum were to occur on the top of the rectangle, at a point (t, x) where t = c, then we would necessarily have ut 0 there, as otherwise u(t, x) would be decreasing as a function of t and hence (c, x) could not be a local maximum on R, and also uxx 0, again contradicting (14.50). To generalize to the case when F (t, x) 0 -- which includes the heat equation when F (t, x) 0, requires a little trick. We set v(t, x) = u(t, x) + x2 , Then, where > 0.
v 2v 2v = - 2 + F (t, x) = F (t, x), t x2 x + F (t, x) = F (t, x) - 2 < 0
where everywhere in R. Thus, by the previous paragraph, the maximum of v occurs on either the bottom or sides of the rectangle. Now we let 0 and conclude the same for u. More precisely, let u(t, x) M on t = 0 or x = a or x = b. Then v(t, x) M + max{ a2 , b2 } and hence, on all of R, u(t, x) v(t, x) M + max{ a2 , b2 }. letting 0 proves that u(t, x) M everywhere, which completes the proof. Q.E.D. Thus, solution u(t, x) to the heat equation (14.12) can only achieve a maximum on the boundary of its domain of definition. In other words, a solution cannot have a maximum at any point (t, x) that lies in the interior of the domain and after the initial time t > 0. Physically, if the temperature in a fully insulated bar starts out everywhere above freezing, then, in the absence of external heat sources, it can never dip below freezing at any later time. Let us apply the Maximum Principle to prove uniqueness of solutions to the heat equation. Theorem 14.4. There is at most one solution to the boundary value problem. 2/25/07 768
c 2006 Peter J. Olver
Proof : Suppose u and u are any two solutions. Then their difference v = u - u solves the homogeneous boundary value problem. Thus, by the maximu principle v(t, x) 0 at all points ofR. But - v = u - u also solves the homogeneous boundary value problem, and hence - v 0 too. This implies v(t, x) 0 and hence u u everywhbere, proving uniqueness. Q.E.D. Remark : Existence follows from the Fourier series solution -- assuming the initial and boundary data and the forcing function are sufficiently nice.
14.3. The Fundamental Solution.
One disadvantage of the Fourier series solution to the heat equation is that it is not nearly as explicit as one might desire for either practical applications, numerical computations, or even further theoretical investigations and developments. An alternative, general approach is based on the idea of the fundamental solution, which derives its inspiration from the Green's function method for solving boundary value problems. For the heat equation, the fundamental solution measures the effect of a concentrated heat source. Let us restrict our attention to homogeneous boundary conditions. The idea is to analyze the case when the initial data u(0, x) = y (x) = (x - y) is a delta function, which we can interpret as a highly concentrated unit heat source, e.g., a soldering iron or laser beam, that is instantaneously applied at a position y along the bar. The heat will diffuse away from its initial concentration, and the resulting fundamental solution is denoted by u(t, x) = F (t, x; y), with F (0, x; y) = (x - y). (14.51)
For each fixed y, the fundamental solution F (t, x; y), considered as a function of t > 0 and x, must satisfy the differential equation, so F 2F = , (14.52) t x2 as well as the specified homogeneous boundary conditions. Once we have determined the fundamental solution, we can then use linear superposition to reconstruct the general solution to the initial-boundary value problem. Namely, we first write the initial data
u(0, x) = f (x) =
0
(x - y) f (y) dy
(14.53)
as a superposition of delta functions, as in (11.36). Linearity implies that the solution is then the corresponding superposition of the responses to those concentrated delta profiles:
u(t, x) =
0
F (t, x; y) f (y) dy.
(14.54)
Assuming that we can differentiate under the integral sign, the fact that F (t, x; y) satisfies the differential equation and the homogeneous boundary conditions for each fixed y immediately implies that the integral (14.54) is also a solution with the correct initial and (homogeneous) boundary conditions. 2/25/07 769
c 2006 Peter J. Olver
Unfortunately, most boundary value problems do not have fundamental solutions that can be written down in closed form. An important exception is the case of an infinitely long homogeneous bar, which requires solving the heat equation u 2u (14.55) = , for - < x < , t > 0. t x2 For simplicity, we have chosen units in which the thermal diffusivity is = 1. The solution u(t, x) is defined for all x R, and has initial conditions u(0, x) = f (x) for - < x < . In order to specify the solution uniquely, we shall require that the temperature be squareintegrable at all times, so that
-
| u(t, x) |2 dx <
for all
t 0.
(14.56)
Roughly speaking, we are requiring that the temperature be small at large distances, which are the relevant boundary conditions for this situation. On an infinite interval, the Fourier series solution to the heat equation becomes a Fourier integral. We write the initial temperature distribution as a superposition 1 f (x) = 2
e i k x f (k) dk,
-
of complex exponentials e i k x , where f (k) is the Fourier transform (13.83) of f (x). The corresponding separable solutions to the heat equation are u(t, x) = e- k
2
t
e i k x = e- k
2
t
cos kx + i sin kx ,
(14.57)
where the frequency variable k is allowed to assume any real value. We invoke linear superposition to combine these complex solutions into a Fourier integral 1 u(t, x) = 2
e- k
-
2
t
e i k x f (k) dk
(14.58)
that forms the solution to the initial value problem for the heat equation. In particular, to recover the fundamental solution, we take the initial temperature profile to be a delta function y (x) = (x - y) concentrated at x = y. According to (13.113), its Fourier transform is e- i ky . y (k) = 2 Plugging this into (14.58), and then referring to our table of Fourier transforms, we find the following explicit formula for the fundamental solution: F (t, x; y) = 2/25/07 1 2
e- k
-
2
t
2 1 e i k(x-y) dk = e-(x-y) /(4 t) . 2 t
(14.59)
Peter J. Olver
770
c 2006
1.4 1.2 1 0.8 0.6 0.4 0.2 -6 -4 -2 2 4 6 -6 -4 -2
1.4 1.2 1 0.8 0.6 0.4 0.2 2 4 6
1.4 1.2 1 0.8 0.6 0.4 0.2 -6 -4 -2 2 4 6 -6 -4 -2
1.4 1.2 1 0.8 0.6 0.4 0.2 2 4 6
Figure 14.3.
The Fundamental Solution to the Heat Equation.
As you can verify, for each fixed y, the function F (t, x; y) is indeed a solution to the heat equation for all t > 0. In addition, lim+ F (t, x; y) =
0, ,
x = y, x = y. (14.60)
t0
Furthermore, its integral
F (t, x; y) dx = 1,
-
which represents the total heat energy, is constant -- in accordance with the law of conservation of energy; see Exercise . Therefore, as t 0+ , the fundamental solution satisfies the original limiting definition (11.2930) of the delta function, and so F (0, x; y) = y (x) has the desired initial temperature profile. In Figure 14.3 we graph F (t, x; 0) at times t = .05, .1, 1., 10.. It starts life as a delta spike concentrated at the origin, and then immediately smoothes out into a tall and narrow bell-shaped curve, centered at x = 0. As time increases, the solution shrinks and widens, decaying everywhere to zero. Its maximal amplitude is proportional to t-1/2 , while its overall width is proportional to t1/2 . The total heat energy (14.60), which is the area under the graph, remains fixed while gradually spreading out over the entire real line. Remark : In probability, these exponentially bell-shaped curves are known as normal or Gaussian distributions. The width of the bell curve corresponds to the standard deviation. For this reason, the fundamental solution to the heat equation sometimes referred to as a "Gaussian filter". Remark : One of the non-physical artifacts of the heat equation is that the heat energy propagates with infinite speed. Indeed, the effect of any initial concentration of heat energy will immediately be felt along the entire length of an infinite bar, because, at any t > 0, the fundamental solution is nonzero for all x. (The graphs in Figure 14.3 are a little misleading because they fail to show the extremely small, but still positive, exponentially decreasing tails.) This effect, while more or less negligible at large distances, is nevertheless 2/25/07 771
c 2006 Peter J. Olver
in clear violation of physical intuition -- not to mention relativity that postulates that signals cannot propagate faster than the speed of light. Despite this non-physical property, the heat equation remains an extremely accurate model for heat propagation and similar diffusive phenomena. With the fundamental solution in hand, we can then adapt the linear superposition formula (14.54) to reconstruct the general solution u(t, x) = 1 2 t
e- (x-y)
-
2
/(4 t)
f (y) dy
(14.61)
to our initial value problem (14.55). Comparing with (13.127), we see that the solutions are obtained by convolution, u(t, x) = g(t, x) f (x), where g(t, x) = F (t, x; 0) = e- x
2
/(4 t)
,
of the initial data with a one-parameter family of progressively wider and shorter Gaussian filters. Since u(t, x) solves the heat equation, we conclude that Gaussian filter convolution has the same smoothing effect on the initial signal f (x). Indeed, the convolution integral (14.61) serves to replace each initial value f (x) by a weighted average of nearby values, the weight being determined by the Gaussian distribution. Such a weighted averaging has the effect of smoothing out high frequency variations in the signal, and, consequently, the Gaussian convolution formula (14.61) provides an effective method for denoising signals and images. In fact, for practical reasons, the graphs displayed earlier in Figure 14.2 were computed by using a standard numerical integration routine to evaluate the convolution integral (14.61), rather than by a numerical solution scheme for the heat equation. Example 14.5. An infinite bar is initially heated to unit temperature along a finite interval. This corresponds to an initial temperature profile u(0, x) = f (x) = (x - a) - (x - b) = 1, 0, a < x < b, otherwise.
The corresponding solution to the heat equation is obtained by the integral formula (14.61), producing u(t, x) = where 1 2 t
b
e- (x-y)
a
2
/(4 t)
dy =
1 erf 2
x 0
x-a 2 t e- z dz
2
- erf
x-b 2 t
,
(14.62)
2 erf x =
(14.63)
is known as the error function due to its applications in probability and statistics, [67]. A graph appears in Figure 14.4. The error function integral cannot be written in terms of elementary functions. Nevertheless, its importance in various applications means that its properties have been well studied, and its values tabulated, [60]. In particular, it has asymptotic values lim erf x = 1, lim erf x = -1. (14.64)
x x -
2/25/07
772
c 2006
Peter J. Olver
1
0.5
-2
-1 -0.5
1
2
-1
Figure 14.4.
1 0.8 0.6 0.4 0.2
The Error Function.
1 0.8 0.6 0.4 0.2 1 0.8 0.6 0.4 0.2
-10
-5 1 0.8 0.6 0.4 0.2
5
10
-10
-5 1 0.8 0.6 0.4 0.2
5
10
-10
-5 1 0.8 0.6 0.4 0.2
5
10
-10
-5
5
10
-10
-5
5
10
-10
-5
5
10
Figure 14.5.
Error Function Solution to the Heat Equation.
A graph of the heat equation solution (14.62) when a = -5, b = 5, at successive times t = 0., .1, 1, 5, 30, 300, is displayed in Figure 14.5. Note the initial smoothing or blurring of the sharp interface, followed by a gradual decay to thermal equilibrium. The Forced Heat Equation The fundamental solution can be also used to solve the inhomogeneous heat equation ut = uxx + h(t, x), (14.65)
that models a bar under an external heat source h(t, x), that might depend upon both position and time. We begin by solving the particular case ut = uxx + (t - s) (x - y), (14.66)
whose inhomogeneity represents a heat source of unit magnitude that is concentrated at a position 0 < y < and applied instantaneously at a single time t = s > 0. Physically, we apply a soldering iron or laser beam to a single spot on the bar for a brief moment. Let 2/25/07 773
c 2006 Peter J. Olver
us also impose homogeneous initial conditions u(0, x) = 0 (14.67)
as well as homogeneous boundary conditions of one of our standard types. The resulting solution u(t, x) = G(t, x; s, y) (14.68) will be referred to as the general fundamental solution to the heat equation. Since a heat source which is applied at time s will only affect the solution at later times t s, we expect that G(t, x; s, y) = 0 for all t < s. (14.69) Indeed, since u(t, x) solves the unforced heat equation at all times t < s subject to homogeneous boundary conditions and has zero initial temperature, this follows immediately from the uniqueness of the solution to the initial-boundary value problem. Once we know the general fundamental solution (14.68), we are able to solve the problem for a general external heat source (14.65) by appealing to linearity. We first write the forcing as a superposition
0
h(t, x) =
0
h(s, y) (t - s) (x - y) dy ds
(14.70)
of concentrated instantaneous heat sources. Linearity allows us to conclude that the solution is given by the self-same superposition formula
t
u(t, x) =
0 0
h(s, y) G(t, x; s, y) dy ds.
(14.71)
The fact that we only need to integrate over times 0 s t follows from (14.69). Remark : If we have a nonzero initial condition, u(0, x) = f (x), then we appeal to linear superposition to write the solution
t
u(t, x) =
0
F (t, x; y) f (y) dy +
0 0
h(s, y) G(t, x; s, y) dy ds
(14.72)
as a combination of (a) the solution with no external heat source, but nonzero initial conditions, plus (b) the solution with homogeneous initial conditions but nonzero heat source. Let us solve the forced heat equation in the case of a infinite bar, so - < x < . We begin by computing the general fundamental solution to (14.66), (14.67). As before, we take the Fourier transform of both sides of the partial differential equation with respect to x. In view of (13.113), (13.117), we find u 1 e- i k y (t - s). + k2 u = t 2 2/25/07 774
c 2006
(14.73)
Peter J. Olver
which is an inhomogeneous first order ordinary differential equation for the Fourier transform u(t, k) of u(t, x). Assuming s > 0, by (14.69), the initial condition is u(0, k) = 0. (14.74)
We solve the initial value problem (14.7374) by the usual method, [27]. Multiplying the 2 differential equation by the integrating factor ek t yields t ek t u
2 2 1 ek t- i k y (t - s), = 2
Integrating both sides from 0 to t and using the initial condition, we find
2 1 ek (t-s)- i k y (t - s), u(t, k) = 2
where (t) is the usual step function (11.42). Finally, we apply the inverse Fourier transform formula (13.86), and then (14.59), to deduce that u(t, x) = G(t, x; s, y) = = (t - s) 2 2 (t - s)
e i k (y-x)+k
-
2
(t-s)
dk = (t - s)F (t - s, x - y) .
(x - y)2 exp - 4(t - s) (t - s)
Thus, the general fundamental solution is obtained by translating the fundamental solution F (t, x; y) for the initial value problem to a starting time of t = s instead of t = 0. Thus, an initial condition has the same aftereffect on the temperature as an instantaneous applied heat source of the same magnitude. Finally, the superposition principle (14.71) produces the solution t h(s, y) (x - y)2 u(t, x) = exp - dy ds. (14.75) 4 (t - s) (t - s) 0 - 2 to the heat equation with source term on an infinite bar. The Root Cellar Problem As a final example, we discuss a problem that involves analysis of the heat equation on a semi-infinite interval. The question is: how deep should you dig a root cellar? In the prerefrigeration era, a root cellar was used to keep food cool in the summer, but not freeze in the winter. We assume that the temperature in the earth only depends on the depth and the time of year. Let u(t, x) denote the deviation in the temperature in the earth, from its annual mean, at depth x > 0 and time t. We shall assume that the temperature at the earth's surface, x = 0, fluctuates in a periodic manner; specifically, we set u(t, 0) = a cos t, where the oscillatory frequency = 2/25/07 2 = 2.0 10-7 sec-1 365.25 days 775
c 2006
(14.76)
(14.77)
Peter J. Olver
refers to yearly temperature variations. In this model, we shall ignore daily temperature fluctuations as their effect is not significant below a very thin surface layer. At large depth the temperature is assumed to be unvarying: u(t, x) - 0 as x - , (14.78)
where 0 refers to the mean temperature. Thus, we must solve the heat equation on a semi-infinite bar 0 < x < , with timedependent boundary conditions (14.76), (14.78) at the ends. The analysis will be simplified a little if we replace the cosine by a complex exponential, and so look for a complex solution with boundary conditions u(t, 0) = a e i t , Let us try a separable solution of the form u(t, x) = v(x) e i t . Substituting this expression into the heat equation ut = uxx leads to i v(x) e i t = v (x) e i t . Canceling the common exponential factors, we conclude that v(x) should solve the boundary value problem v (x) = i v, v(0) = a,
x x
lim u(t, x) = 0.
(14.79)
(14.80)
lim v(x) = 0.
The solutions to the ordinary differential equation are v2 (x) = e- i / x = e- /2 (1+ i ) x . v1 (x) = e i / x = e /2 (1+ i ) x , The first solution is exponentially growing as x , and so not appropriate to our problem. The solution to the boundary value problem must therefore be a multiple, - /2 (1+ i ) x v(x) = a e of the exponentially decaying solution. Substituting back into (14.80), we find the (complex) solution to the root cellar problem to be - x /2 i t- /2 x e u(t, x) = a e . (14.81) The corresponding real solution is obtained by taking the real part, u(t, x) = a e
-x
/2
cos
t-
x 2
.
(14.82)
The first term in (14.82) is exponentially decaying as a function of the depth. Thus, the further down one goes, the less noticeable the effect of the surface temperature fluctuations. The second term is periodic with the same annual frequency . The interesting feature is 2/25/07 776
c 2006 Peter J. Olver
the phase lag in the response. The temperature at depth x is out of phase with respect to the surface temperature fluctuations, having an overall phase lag = x 2
that depends linearly on depth. In particular, a cellar built at a depth where is an odd multiple of will be completely out of phase, being hottest in the winter, and coldest in the summer. Thus, the (shallowest) ideal depth at which to build a root cellar would take = , corresponding to a depth of x= 2 .
For typical soils in the earth, 10-6 meters2 sec-1 , and hence, by (14.77), x 9.9 meters. However, at this depth, the relative amplitude of the oscillations is e- x /2 = e- = .04 and hence there is only a 4% temperature fluctuation. In Minnesota, the temperature varies, roughly, from - 40 C to + 40 C, and hence our 10 meter deep root cellar would experience only a 3.2 C annual temperature deviation from the winter, when it is the warmest, to the summer, where it is the coldest. Building the cellar twice as deep would lead to a temperature fluctuation of .2%, now in phase with the surface variations, which means that the cellar is, for all practical purposes, at constant temperature year round.
14.4. The Wave Equation.
The second important class of dynamical partial differential equations are those modeling vibrations of continuous media. As we saw in Chapter 9, Newton's Law implies that the free vibrations of a discrete mechanical system are governed by a second order system of ordinary differential equations of the form M d2 u = - K u, dt2
in which M is the positive definite, diagonal mass matrix, while K = A A = AT CA is the positive definite (or semi-definite in the case of an unstable system) stiffness matrix. The corresponding dynamical equations describing the small vibrations of continuous media take an entirely analogous form 2u = - K[ u ]. t2 (14.83)
In this framework, describes the density, while K = L L is the same self-adjoint differential operator, with appropriate boundary conditions, that appears in the equilibrium 2/25/07 777
c 2006 Peter J. Olver
equations (11.89). For one-dimensional media, such as a vibrating bar or string, we are led to a partial differential equation in the particular form (x) 2u = 2 t x (x) u x , 0 < x < , (14.84)
where (x) is the density of the bar or string at position x, while (x) > 0 denotes its stiffness or tension. The second order partial differential equation (14.84) models the dynamics of vibrations and waves in a broad range of applications, including elastic vibrations of a bar, sound vibrations in a column of air, e.g., inside a wind instrument, and also transverse vibrations of a string, e.g., a violin string. (However, bending vibrations of a beam lead to a fourth order partial differential equation; see Exercise .) The wave equation is also used to model small amplitude water waves, electromagnetic waves, including light, radio and microwaves, gravitational waves, and many others. A detailed derivation of the model from first principles in the case of a vibrating string can be found in [183]. To specify the solution, we must impose suitable boundary conditions. The usual suspects -- Dirichlet, Neumann, mixed, and periodic boundary conditions -- continue to play a central role, and have immediate physical interpretations. Tying down an end of the string imposes a Dirichlet condition u(t, 0) = , while a free end is prescribed by a homogeneous Neumann boundary condition ux (t, 0) = 0. Periodic boundary conditions, as in (14.11), correspond to the vibrations of a circular ring. As with all second order Newtonian systems of ordinary differential equations, the solution to the full boundary value problem for the second order partial differential equation will then be uniquely specified by its initial displacement and initial velocity: u(0, x) = f (x), u (0, x) = g(x). t (14.85)
The simplest situation occurs when the medium is homogeneous, and so both its density and stiffness are constant. Then the general vibration equation (14.84) reduces to the one-dimensional wave equation 2u 2u = c2 2 . t2 x The constant c= > 0 (14.87) (14.86)
is known as the wave speed , for reasons that will soon become apparent. The method for solving such second order systems is motivated by our solution in the discrete case discussed in Section 9.5. To keep matters simple, we shall concentrate on the wave equation (14.86), although the method is easily extended to the general homogeneous Newtonian system (14.84). We will try a separable solution with trigonometric time dependence: u(t, x) = cos( t) v(x), (14.88) 2/25/07 778
c 2006 Peter J. Olver
in which both the frequency and profile v(x) are to be determined. Differentiating (14.88), we find 2u = - 2 cos( t) v(x), 2 t 2u = cos( t) v (x). 2 x
Substituting into the wave equation (14.86) and canceling the common cosine factors, we deduce that v(x) must satisfy the ordinary differential equation c2 d2 v + 2 v = 0. 2 dx (14.89)
Thus 2 = can be viewed as an eigenvalue with eigenfunction v(x) for the second order differential operator K = - c2 D2 . Ignoring the boundary conditions for the moment, if x x > 0, the solutions are the trigonometric functions cos c , sin c . Thus, (14.88) leads to the pair of explicit solutions cos t cos x , c cos t sin x , c
to the wave equation. Now, the computation will work just as well with a sine function, which yields two additional solutions sin t cos x , c sin t sin x . c
Each of these four solutions represents a spatially periodic standing wave form of period 2 c/, that is vibrating with frequency . Observe that the smaller scale waves vibrate faster. On the other hand, if = 0, then (14.89) has the solution v = x + , leading to the solutions u(t, x) = 1, and u(t, x) = x. (14.90) The first is a constant, nonvibrating solution, while the second is also constant in time, but will typically not satisfy the boundary conditions and so can be discarded. As we learned in Chapter 9, the existence of a zero eigenvalue corresponds to an unstable mode in the physical system, in which the displacement grows linearly in time. In the present situation, these correspond to the two additional solutions u(t, x) = t, and u(t, x) = x t, (14.91)
of the wave equation. Again, the second solution will typically not satisfy the homogeneous boundary conditions, and can usually be ignored. Such null eigenfunction modes will only arise in unstable situations. The boundary conditions will distinguish the particular eigenvalues and natural frequencies of vibration. Consider first the case of a string of length with two fixed ends, and thus subject to homogeneous Dirichlet boundary conditions u(t, 0) = 0 = u(t, ). 2/25/07 779
c 2006 Peter J. Olver
This forms a positive definite boundary value problem, and so there is no unstable mode. Indeed, the eigenfunctions of the boundary value problem (14.89) with Dirichlet boundary conditions v(0) = 0 = v() were found in (14.17): n x n c , n = , n = 1, 2, 3, . . . . Therefore, we can write the general solution as a Fourier sine series vn (x) = sin
u(t, x) =
n=1
bn cos rn cos
n=1
n x n c t n x n c t sin + dn sin sin n c t + n n x sin . (14.92)
=
The solution is thus a linear combination of the natural Fourier modes vibrating with frequencies n c n n = = , n = 1, 2, 3, . . . . (14.93) The longer the length of the string, or the higher its density , the slower the vibrations; whereas increasing its stiffness or tension speeds them up -- in exact accordance with our physical intuition. The Fourier coefficients bn and dn in (14.92) will be uniquely determined by the initial conditions (14.85). Differentiating the series term by term, we discover that we must represent the initial displacement and velocity as Fourier sine series
u(0, x) =
n=1
bn sin
n x = f (x), 2
u n x n c (0, x) = sin = g(x). dn t n=1 n x dx, n = 1, 2, 3, . . . . (14.94)
Therefore, bn = f (x) sin
0
are the Fourier sine coefficients (12.84) of the initial displacement f (x), while dn = 2 n c
g(x) sin
0
n x dx,
n = 1, 2, 3, . . . .
(14.95)
are rescaled versions of the Fourier sine coefficients of the initial velocity g(x). Example 14.6. A string of unit length is held taut in the center and then released. Our goal is to describe the ensuing vibrations. Let us assume the physical units are chosen so that c2 = 1, and so we are asked to solve the initial-boundary value problem utt = uxx , u(0, x) = f (x), ut (0, x) = 0, u(t, 0) = u(t, 1) = 0. (14.96)
To be specific, we assume that the center of the string has been displaced by half a unit, and so the initial displacement is x, 0 x 1, 2 f (x) = 1 1 - x, x 1. 2 2/25/07 780
c 2006 Peter J. Olver
0.4 0.2
0.4 0.2
0.4 0.2
0.2 -0.2 -0.4
0.4
0.6
0.8
1 -0.2 -0.4
0.2
0.4
0.6
0.8
1 -0.2 -0.4
0.2
0.4
0.6
0.8
1
0.4 0.2
0.4 0.2
0.4 0.2
0.2 -0.2 -0.4
0.4
0.6
0.8
1 -0.2 -0.4
0.2
0.4
0.6
0.8
1 -0.2 -0.4
0.2
0.4
0.6
0.8
1
Figure 14.6.
Plucked String Solution of the Wave Equation.
The vibrational frequencies n = n are the integral multiples of , and so the natural modes of vibration are cos n t sin n x
and
sin n t sin n x
for
n = 1, 2, . . . .
Consequently, the general solution to the boundary value problem is u(t, x) =
n=1
bn cos n t sin n x + dn sin n t sin n x ,
1/2
where
1
bn = 2
f (x) sin n x dx =
0
4
0
x sin n x dx =
4 (- 1)k , (2 k + 1)2 2
n = 2 k + 1, n = 2 k,
0,
are the Fourier sine coefficients of the initial displacement, while dn = 0 are the Fourier sine coefficients of the initial velocity. Therefore, the solution is the Fourier sine series
u(t, x) = 4
k=0
(- 1)k
cos(2 k + 1) t sin(2 k + 1) x , (2 k + 1)2 2
(14.97)
whose graph at times t = 0, .2, .4, .6, .8, 1. is depicted in Figure 14.6. At time t = 1, the original displacement is reproduced exactly, but upside down. The subsequent dynamics proceeds as before, but in mirror image form. The original displacement reappears at time t = 2, after which time the motion is periodically repeated. Interestingly, at times tk = .5, 1.5, 2.5, . . ., the displacement is identically zero: u(tk , x) 0, although the velocity ut (tk , x) 0. The solution appears to be piecewise affine, i.e., its graph is a collection of straight lines. This fact will be verified in Exercise , where you are asked to construct an exact analytical formula for this solution. Unlike the heat equation, the wave equation does not smooth out discontinuities and corners in the initial data. While the series form (14.92) of the solution is not entirely satisfying, we can still use it to deduce important qualitative properties. First of all, since each term is periodic in t with 2/25/07 781
c 2006 Peter J. Olver
period 2 /c, the entire solution is time periodic with that period: u(t + 2 /c, x) = u(t, x). In fact, after half the period, at time t = /c, the solution reduces to u ,x c n x = =- (-1) bn sin n=1
n
bn sin
n=1
n ( - x) = - u(0, - x) = - f ( - x).
In general, u t+ , x c = - u(t, - x), u t+ 2 , x c = u(t, x). (14.98)
Therefore, the initial wave form is reproduced, first as an upside down mirror image of itself at time t = /c, and then in its original form at time t = 2 /c. This has the important consequence that vibrations of (homogeneous) one-dimensional media are purely periodic phenomena! There is no quasi-periodicity because the fundamental frequencies are all integer multiples of each other. Remark : The preceding analysis has important musical consequences. To the human ear, sonic vibrations that are integral multiples of a single frequency are harmonic, whereas those that admit quasi-periodic vibrations, with irrationally related frequencies, sound percussive. This is why most tonal instruments rely on vibrations in one dimension, be it a violin string, a column of air in a wind instrument (flute, clarinet, trumpet or saxophone), a xylophone bar or a triangle. On the other hand, most percussion instruments rely on the vibrations of two-dimensional media, e.g., drums and cymbals, or three-dimensional solid bodies, e.g., blocks. As we shall see in Chapters 17 and 18, the frequency ratios of the latter are typically irrational. A bar with both ends left free, and so subject to the Neumann boundary conditions u u (t, 0) = 0 = (t, ), x x (14.99)
will have a slightly different behavior, owing to the instability of the underlying equilibrium equations. The eigenfunctions of (14.89) with Neumann boundary conditions v (0) = 0 = v () are now vn (x) = cos n x with n = n c , n = 0, 1, 2, 3, . . . .
The resulting solution takes the form of a Fourier cosine series
u(t, x) = a0 + c0 t +
n=1
an cos
n x n c t n x n c t cos + cn sin cos
.
(14.100)
In accordance with (14.90), the first two terms come from the null eigenfunction v0 (x) = 1 with 0 = 0. The bar vibrates with the same fundamental frequencies (14.93) as in the fixed end case, but there is now an additional unstable mode c0 t that is no longer periodic, but grows linearly in time. 2/25/07 782
c 2006 Peter J. Olver
Substituting (14.100) into the initial conditions (14.85), we find the Fourier coefficients are prescribed, as before, by the initial displacement and velocity, an = 2
f (x) cos
0
n x dx,
cn =
2 n c
g(x) cos
0
n x dx,
n = 1, 2, 3, . . . .
The order zero coefficients , a0 = 1
f (x) dx,
0
c0 =
1
g(x) dx,
0
are equal to the average initial displacement and average initial velocity of the bar. In particular, when c0 = 0 there is no net initial velocity, and the unstable mode is not excited. In this case, the solution is time-periodic, oscillating around the position given by the average initial displacement. On the other hand, if c0 = 0, the bar will move off with constant average speed c0 , all the while vibrating at the same fundamental frequencies. Similar considerations apply to the periodic boundary value problem for the wave equation on a circular ring. The details are left as Exercise for the reader. Forcing and Resonance In Section 9.6, we learned that periodically forcing an undamped mechanical structure (or a resistanceless electrical circuit) at a frequency that is distinct from its natural vibrational frequencies leads, in general, to a quasi-periodic response. The solution is a sum of the unforced vibrations superimposed with an additional vibrational mode at the forcing frequency. However, if forced at one of its natural frequencies, the system may go into a catastrophic resonance. The same type of quasi-periodic/resonant response is also observed in the partial differential equations governing periodic vibrations of continuous media. To keep the analysis as simple as possible, we restrict our attention to the forced wave equation for a homogeneous bar 2u 2u = c2 2 + F (t, x), (14.101) t2 x subject to specified homogeneous boundary conditions. The external forcing function F (t, x) may depend upon both time t and position x. We will be particularly interested in a periodically varying external force of the form F (t, x) = cos( t) h(x), (14.102)
where the function h(x) is fixed. As always, the solution to an inhomogeneous linear equation can be written as a combination, u(t, x) = u (t, x) + z(t, x) (14.103)
Note that, we have not included the usual (14.100).
1 2
factor in the constant terms in the Fourier series
2/25/07
783
c 2006
Peter J. Olver
of a particular solution u (t, x) plus the general solution z(t, x) to the homogeneous equation, namely 2z 2z = c2 2 . (14.104) t2 x The boundary and initial conditions will serve to uniquely prescribe the solution u(t, x), but there is some flexibility in its two constituents (14.103). For instance, we may ask that the particular solution u satisfy the homogeneous boundary conditions along with zero (homogeneous) initial conditions, and thus represents the pure response of the system to the forcing. The homogeneous solution z(t, x) will then reflect the effect of the initial and boundary conditions unadulterated by the external forcing. The final solution will equal the sum of the two individual responses. In the case of periodic forcing (14.102), we look for a particular solution u (t, x) = cos( t) v (x) (14.105)
that vibrates at the forcing frequency. Substituting the ansatz (14.105) into the equation (14.101), and canceling the common cosine factors, we discover that v (x) must satisfy the boundary value problem prescribed by
- c2 v - 2 v = h(x),
(14.106)
supplemented by the relevant homogeneous boundary conditions -- Dirichlet, Neumann, mixed, or periodic. At this juncture, there are two possibilities. If the unforced, homogeneous boundary value problem has only the trivial solution v 0, then there is a solution to the forced boundary value problem for any form of the forcing function h(x). On the other hand, the homogeneous boundary value problem has a nontrivial solution v(x) when 2 = is an eigenvalue, and so is a natural frequency of vibration to the homogeneous problem; the solution v(x) is the corresponding eigenfunction appearing in the solution series (14.92). In this case, according to the Fredholm alternative, Theorem 5.55, the boundary value problem (14.106) has a solution if and only if the forcing function h(x) is orthogonal to the eigenfunction(s): h , v = 0. (14.107) See Example 11.3 and Exercise 11.5.7 for details. If we force in a resonant manner -- meaning that (14.107) is not satisfied -- then the solution will be a resonantly growing vibration u (t, x) = t sin( t) v (x). In a real-world situation, such large resonant (or near resonant) vibrations will either cause a catastrophic breakdown, e.g., the bar breaks or the string snaps, or will send the system into a different, nonlinear regime that helps mitigate the resonant effects, but is no longer modeled by the simple linear wave equation. Example 14.7. As a specific example, consider the initial-boundary value problem modeling the forced vibrations of a uniform bar of unit length and fixed at both ends: utt = c2 uxx + cos( t) h(x), u(t, 0) = 0 = u(t, 1), 2/25/07 u(0, x) = f (x), 784 ut (0, x) = g(x).
c 2006
(14.108)
Peter J. Olver
The particular solution will have the nonresonant (14.105) form provided there exists a solution v (x) to the boundary value problem
c2 v + 2 v = - h(x),
v (0) = 0 = v (1).
(14.109)
The natural frequencies and associated eigenfunctions are n = n c , vn (x) = sin n x, n = 1, 2, 3, . . . .
The boundary value problem (14.109) will have a solution, and hence the forcing is not resonant, provided either = n is not a natural frequency, or = n , but
1
0 = h , vn =
h(x) sin n x dx
0
(14.110)
is orthogonal to the associated eigenfunction. Otherwise, the forcing profile will induce a resonant response. For example, under periodic forcing of frequency with trigonometric profile h(x) sin k x, the particular solution to (14.109) is v (x) = 2 sin k x , - k 2 2 c2 so that u (t, x) = cos t sin k x , 2 - k 2 2 c2 (14.111)
which is a valid solution as long as = k = k c. Note that we may allow the forcing frequency = n to coincide with any other resonant forcing frequency, n = k, because the sine profiles are mutually orthogonal and so the nonresonance condition (14.110) holds. On the other hand, if = k = k c, then the particular solution u (t, x) = t sin k c t sin k x , 2k c (14.112)
is resonant, and grows linearly in time. To obtain the full solution to the initial-boundary value problem, we write u = u + z where z(t, x) must satisfy ztt - c2 zxx = 0, along with the modified initial conditions z(0, x) = f (x) - sin k x , 2 - k 2 2 c2 u (0, x) = g(x), x z(t, 0) = 0 = z(t, 1),
stemming from the fact that the particular solution (14.111) has non-trivial initial data. (In the resonant case (14.112), there is no extra term in the initial data.) Note that, the closer is to the resonant frequency, the larger the modification of the initial data, and hence the larger the response of the system to the periodic forcing. As before, the solution z(t, x) to the homogeneous equation can be written as a Fourier sine series (14.92). The final formulae are left to the reader to write out. 2/25/07 785
c 2006 Peter J. Olver
14.5. d'Alembert's Solution.
The one-dimensional wave equation is distinguished by admitting an explicit solution formula, originally discovered by the eighteenth century French mathematician Jean d'Alembert, that entirely avoids the complicated Fourier series form. D'Alembert's solution provides additional valuable insight into the behavior of the solutions. Unfortunately, unlike series methods that have a very broad range of applicability, d'Alembert's method only succeeds in this one very special situation: the homogeneous wave equation in a single space variable. The starting point is to write the wave equation (14.86) in the suggestive form
2 2 u = (t - c2 x ) u = utt - c2 uxx = 0.
(14.113)
2 2 Here = t - c2 x is a common mathematical notation for the wave operator , while t , x are convenient shorthands for the partial derivative operators with respect to t and x. We note that is a linear, second order partial differential operator. In analogy with the elementary polynomial factorization
t2 - c2 x2 = (t - c x)(t + c x), we shall factor the wave operator into a product of two first order partial differential operators: 2 2 = t - c2 x = (t - c x ) (t + c x ). (14.114) Now, if the second factor annihilates the function u(t, x), meaning (t + c x ) u = ut + c ux = 0, then u is automatically a solution to the wave equation, since u = (t - c x ) (t + c x ) u = (t - c x ) 0 = 0. In other words, every solution to the simpler first order partial differential equation (14.115) is a solution to the wave equation (14.86). (The converse is, of course, not true.) It is relatively easy to solve linear first order partial differential equations. Proposition 14.8. Every solution u(t, x) to the partial differential equation u u +c =0 t x has the form u(t, x) = p(x - c t), where p() is an arbitrary function of the characteristic variable = x - c t.
See Chapter 22 for more details, including extensions to first order nonlinear partial differential equations.
(14.115)
(14.116) (14.117)
2/25/07
786
c 2006
Peter J. Olver
1
1
1
0.5
0.5
0.5
0.2 -0.5
0.4
0.6
0.8
1
1.2
1.4 -0.5
0.2
0.4
0.6
0.8
1
1.2
1.4 -0.5
0.2
0.4
0.6
0.8
1
1.2
1.4
-1
-1
-1
Figure 14.7.
Traveling Wave.
Proof : We adopt a linear change of variables to rewrite the solution u(t, x) = p(t, x - c t) = p(t, ) as a function of the characteristic variable and the time t. Applying the chain rule, we express the derivatives of u in terms of the derivatives of p as follows: u p p = -c , t t and hence u p = , x
u p p p p u +c = -c +c = . t x t t
Therefore, u is a solution to (14.116) if and only if p(t, ) is a solution to the very simple partial differential equation p = 0. t This clearly implies that p = p() does not depend on the variable t, and hence u = p() = p(x - c t) is of the desired form. Q.E.D. Therefore, any function of the characteristic variable, e.g., 2 + 1, or cos , or e , will produce a corresponding solution, (x - c t)2 + 1, or cos(x - c t), or ex-c t , to the first order partial differential equation (14.116), and hence a solution to the wave equation (14.86). The functions of the form u(t, x) = p(x - c t) are known as traveling waves. At t = 0 the wave has the initial profile u(0, x) = p(x). As t progresses, the wave moves to the right, with constant speed c > 0 and unchanged in form; see Figure 14.7. For this reason, (14.116) is sometimes referred to as the one-way or unidirectional wave equation. Proposition 14.8 tells us that every traveling wave with wave speed c is a solution to the full wave equation (14.86). But keep in mind that such solutions do not necessarily respect the boundary conditions, which, when present, will affect their ultimate behavior.
More rigorously, one should also assume that, at each time t, the domain of definition of p() is a connected interval. A similar technical restriction should be imposed upon the solutions in the statement of Proposition 14.8. See Exercise for a detailed example.
2/25/07
787
c 2006
Peter J. Olver
Now, since c is constant, the factorization (14.114) can be written equally well in the reverse order: 2 2 = t - c2 x = (t + c x ) (t - c x ). (14.118) The same argument tells us that any solution to the alternative first order partial differential equation u u -c = 0, (14.119) t x also provides a solution to the wave equation. This is also a unidirectional wave equation, but with the opposite wave speed - c. Applying Proposition 14.8, now with c replaced by - c, we conclude that the general solution to (14.119) has the form u(t, x) = q(x + c t) (14.120) where q() is an arbitrary differentiable function of the alternate characteristic variable = x + c t. The solutions (14.120) represent traveling waves moving to the left with constant speed c > 0 and unchanged in form. The wave equation (14.113) is bidirectional and admits both left and right traveling wave solutions. Linearity of the wave equation implies that the sum of solutions is again a solution. In this way, we can produce solutions which are superpositions of left and right traveling waves. The remarkable fact is that every solution to the wave equation can be so represented. Theorem 14.9. Every solution to the wave equation (14.86) can be written as a combination u(t, x) = p() + q() = p(x - c t) + q(x + c t) (14.121) of right and left traveling waves, each depending on its characteristic variable = x - c t, = x + c t. (14.122)
Proof : The key is to use a linear changes of variables to rewrite the wave equation entirely in terms of the characteristic variables. We set u(t, x) = w(x - c t, x + c t) = w(, ), whereby w(, ) = u - + , 2c 2 .
Then, using the chain rule to compute the partial derivatives, u =c t and hence 2u = c2 t2 Therefore u= 2/25/07 2w 2w 2w -2 + 2 2 , 2u 2w 2w 2w = +2 + . x2 2 2 w w - , u w w = + . x
2u 2w 2u - c2 2 = - 4 c2 . t2 x 788
c 2006 Peter J. Olver
We conclude that u(t, x) solves the wave equation second order partial differential equation 2w = 0,
u = 0 if and only if w(, ) solves the w
which we write in the form
= 0.
This partial differential equation can be integrated once with respect to , resulting in w = r(), where r is an arbitrary function of the characteristic variable . Integrating both sides of the latter partial differential equation with respect to , we find w(, ) = p() + q(), where q () = r(),
while p() represents the integration "constant". Replacing the characteristic variables by their formulae in terms of t and x completes the proof. Q.E.D. Remark : As noted above, we have been a little cavalier with our specification of the domain of definition of the functions and the differentiability assumptions required. Sorting out the precise technical details is left to the meticulous reader. Remark : As we know, the general solution to a second order ordinary differential equation depends on two arbitrary constants. Here we observe that the general solution to a second order partial differential equation depends on two arbitrary functions -- in this case p() and q(). This observation serves as a useful rule of thumb, but should not be interpreted too literally. Let us now see how this new form of solution to the wave equation can be used to effectively solve initial value problems. The simplest case is that of a bar or string of infinite length, in which case we have a pure initial value problem 2u 2u u = c2 2 , u(0, x) = f (x), (0, x) = g(x), for - < x < . t2 x t Substituting the solution formula (14.121) into the initial conditions, we find u (0, x) = - c p (x) + c q (x) = g(x). t To solve this pair of linear equations for p and q, we differentiate the first equation: u(0, x) = p(x) + q(x) = f (x), p (x) + q (x) = f (x). Subtracting the second equation divided by c, we find 2 p (x) = f (x) - Therefore, p(x) = 2/25/07 1 1 f (x) - 2 2c 789
0 c 2006 Peter J. Olver
1 g(x). c
x
g(z) dz + a,
1
1
1
0.5
0.5
0.5
0.5 -0.5
1
1.5
2 -0.5
0.5
1
1.5
2 -0.5
0.5
1
1.5
2
-1
-1
-1
1
1
1
0.5
0.5
0.5
0.5 -0.5
1
1.5
2 -0.5
0.5
1
1.5
2 -0.5
0.5
1
1.5
2
-1
-1
-1
Figure 14.8.
Interaction of Waves.
where a is an integration constant. The first equation then yields 1 1 q(x) = f (x) - p(x) = f (x) + 2 2c
x 0
g(z) dz - a.
Substituting these two expressions back into (14.121), we find 1 f () + f () + u(t, x)= p() + q() = 2 2c = 1 f () + f () + 2 2c
-
+
0 0
g(z) dz
g(z) dz,
where , are the characteristic variables (14.122). In this fashion, we have derived d'Alembert's solution to the wave equation on the entire line - < x < . Theorem 14.10. The solution to the initial value problem 2u 2u = c2 2 , t2 x is given by u(0, x) = f (x), u (0, x) = g(x), t
x+c t
- < x < . g(z) dz.
(14.123)
f (x - c t) + f (x + c t) 1 u(t, x) = + 2 2c
(14.124)
x-c t
. Let us investigate the implications of d'Alembert's solution formula (14.124). First, suppose there is no initial velocity, so g(x) 0, and the motion is purely the result of the initial displacement u(0, x) = f (x). In this case, the solution (14.124) reduces to u(t, x) =
1 2
f (x - c t) +
1 2
f (x + c t).
The effect is that the initial displacement f (x) splits into two waves, one traveling to the right and one traveling to the left, each of constant speed c, and each of exactly the same 2/25/07 790
c 2006 Peter J. Olver
shape as the initial displacement f (x) but only half as tall. For example, if the initial displacement is a localized pulse centered at the origin, say u(0, x) = e- x , then the solution u(t, x) =
1 2
2
u (0, x) = 0, t
2
e- (x-c t) +
1 2
e- (x+c t)
2
consists of two half size pulses running away from the origin in opposite directions with equal speed c. If we take two separated pulses, say u(0, x) = e- x + 2 e- (x-1) , centered at x = 0 and x = 1, then the solution u(t, x) =
1 2
2 2
u (0, x) = 0, x
e- (x-c t) + e- (x-1-c t) +
2
2
1 2
e- (x+c t) + e- (x-1+c t)
2
2
will consist of four pulses, two moving to the right and two to the left, all with the same speed, as pictured in Figure 14.8. An important observation is that when a right-moving pulse collides with a left-moving pulse, they emerge from the collision unchanged -- a consequence of the inherent linearity of the wave equation. The first picture shows the initial displacement. In the second and third pictures, the two localized bumps have each split into two copies moving in opposite directions. In the fourth and fifth, the larger right moving bump is in the process of interacting with the smaller left moving bump. Finally, in the last picture the interaction is complete, and the two left moving bumps and two right moving bumps travel in tandem with no further collisions. Remark : If the initial displacement has bounded support, and so f (x) = 0 if x < a or x > b for some a < b, then after a finite time the right and left-moving waves will be completely separated, and the observer will see two exact half size replicas running away, with speed c, in opposite directions. If the displacement is not localized, then the left and right traveling waves will never fully disengage, and one might be hard pressed (just as in our earlier discussion of quasi-periodic phenomena) in recognizing that a complicated solution pattern is, in reality, just the superposition of two simple traveling waves. For example, using a trigonometric identity, any separable solution, e.g., cos c t sin x =
1 2
cos(x - c t) +
1 2
cos(x + c t),
can be rewritten in d'Alembert form. The interpretation of the two solution formulas is quite different. Most observers will see a standing sinusoidal wave, vibrating with frequency c, as represented on the left side of the equation. However, the right hand side says that this is the same as the difference of a right and left traveling cosine wave. The interactions of their peaks and troughs reproduces the standing wave. Thus, the same solution can be interpreted in seemingly incompatible ways! 2/25/07 791
c 2006 Peter J. Olver
5
4
3
2
1
0.5
1
1.5
2
2.5
3
-1
Figure 14.9.
Characteristic Lines for the Wave Equation.
are known as the characteristics of the wave equation. The two characteristics going through a point on the x axis, where the initial data is prescribed, are illustrated in Figure 14.9. The reader should note that, in this figure, the t axis is horizontal, while the x axis is vertical. In general, signals propagate along characteristics. More specifically, if we start out with an initial displacement concentrated very close to a point x = a, t = 0, then the solution will be concentrated along the two characteristic lines emanating from the point, namely x-c t = a and x+c t = a. In the limit, a unit impulse or delta function displacement at x = a, corresponding to the initial condition u(0, x) = (x - a), will result in a solution u(t, x) =
1 2
In general, the lines of slope c in the (t, x)plane, where the characteristic variables are constant, = x - c t = a, = x + c t = b, (14.125)
u (0, x) = 0, t
1 2
(14.126)
(x - c t - a) +
(x + c t - a)
(14.127)
consisting of two half-strength delta spikes traveling away from the starting position concentrated on the two characteristic lines. Let us return to the general initial value problem (14.123). Suppose now that there is no initial displacement, u(0, x) = f (x) 0, but rather a concentrated initial velocity, say a delta function u (0, x) = a (x) = (x - a). t Physically, this would correspond to striking the string with a highly concentrated blow at the point x = a. The d'Alembert solution (14.124) to this "hammer blow" problem is 1 , x+c t x - c t < a < x + c t, 1 2c (14.128) a (z) dz = u(t, x) = 2 c x-c t 0, otherwise, 2/25/07 792
c 2006 Peter J. Olver
5
(0, x + c t)
4
3
2
(t, x)
1
0.5
1
1.5
2
2.5
3
(0, x - c t)
-1
Figure 14.10.
1.2 1 0.8 0.6 0.4 0.2 0.2 -0.2 0.4 0.6 0.8 1 1.2 1.4 -0.2 1.2 1 0.8 0.6 0.4 0.2 0.2
Domain of Dependence.
1.2 1 0.8 0.6 0.4 0.2 0.4 0.6 0.8 1 1.2 1.4 -0.2 0.2 0.4 0.6 0.8 1 1.2 1.4
Figure 14.11.
Concentrated Initial Velocity for Wave Equation.
and consists of a constant displacement, of magnitude 1/(2 c), between the two characteristic lines x - c t = a = x + c t based at x = a, t = 0 -- the shaded region of Figure 14.9. This region is known as the domain of influence of the point (a, 0) since, in general, the value of the initial data at that point will only affect the solution values in the region. Vice versa, the value of the solution u(t, x) a point with t > 0 depends only on that part of the initial data that lies within the domain of dependence of the point (t, x), which is the triangle whose sides have slope c extending back to the x-axis; see Figure 14.10. The particular solution, which is plotted in Figure 14.11, has two jump discontinuities between the undisturbed state and the displaced state, each propagating along its characteristic line with speed c, but in opposite directions. Note that, unlike a concentrated initial displacement, where the signal remains concentrated and each point along the bar is temporarily displaced, eventually returning to its undisturbed state, a concentrated initial velocity has a lasting effect, and the bar remains permanently deformed by an amount 1/(2 c) . Solutions on Bounded Intervals So far, we have been using d'Alembert formula to solve the wave equation on an infinite interval. The formula can still be used on bounded intervals, but in a suitably modified format so as to respect the boundary conditions. The easiest to deal with is the 2/25/07 793
c 2006 Peter J. Olver
1 0.5 -4 -2 -0.5 -1 2 4
Figure 14.12.
Odd Periodic Extension of a Concentrated Pulse.
periodic problem on 0 x , with boundary conditions u(t, 0) = u(t, ), ux (t, 0) = ux (t, ). (14.129) If we extend the initial displacement f (x) and velocity g(x) to be periodic functions of period , so f (x+) = f (x) and g(x+) = g(x) for all x R, then the resulting d'Alembert solution (14.124) will also be periodic in x, so u(t, x + ) = u(t, x). In particular, it satisfies the boundary conditions (14.129) and so coincides with the desired solution. Details can be found in Exercises . Next, suppose we have fixed (Dirichlet) boundary conditions u(t, 0) = 0, u(t, ) = 0. (14.130)
The resulting solution can be written as a Fourier sine series (14.92), and hence is both odd and 2 periodic in x. Therefore, to write the solution in d'Alembert form, we extend the initial displacement f (x) and velocity g(x) to be odd, periodic functions of period 2 : f (- x) = - f (x), f (x + 2 ) = f (x), g(- x) = - g(x), g(x + 2 ) = g(x).
This will ensure that the d'Alembert solution also remains odd and periodic. As a result, it satisfies the boundary conditions (14.130) for all t. Keep in mind that, while the solution u(t, x) is defined for all x, the only physically relevant values occur on the interval 0 x . Nevertheless, the effects of displacements in the unphysical regime will eventually be "felt" as the propagating waves pass through the physical interval. For example, consider an initial displacement which is concentrated near x = a for some 0 < a < . Its odd, periodic extension consists of two sets of replicas: those of the same form occurring at positions a 2 , a 4 , . . . , and their mirror images at the intermediate positions - a, - a 2 , - a 4 , . . . ; Figure 14.12 shows a representative example. The resulting solution begins with each of the pulses, both positive and negative, splitting into two half-size replicas that propagate with speed c in opposite directions. As the individual pulses meet, they interact as they pass unaltered through each other. The process repeats periodically, with an infinite row of half-size pulses moving to the right kaleidoscopically interacting with an infinite row moving to the left. However, only the part of this solution that lies on 0 x is actually realized on the physical bar. The net effect is as if we were forced to view the complete solution as it passes by a window of length that blocks out all other regions of the real axis. What the viewer effectively sees assumes a somewhat different interpretation. To wit, the original pulse at position 0 < a < splits up into two half size replicas that move off in opposite directions. As each half-size pulse reaches an end of the bar, it meets a mirror image pulse that has been propagating in the opposite direction from the non-physical regime. The pulse appears to be reflected at the end of the interval, and changes into an upside down 2/25/07 794
c 2006 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
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
Figure 14.13.
Solution to Wave Equation with Fixed Ends.
mirror image moving in the opposite direction. The original positive pulse has moved off the end of the bar just as its mirror image has moved into the physical regime. (A common physical illustration is a pulse propagating down a jump rope that is held fixed at its end; the reflected pulse returns upside down.) A similar reflection occurs as the other half-size pulse hits the other end of the physical interval, after which the solution consists of two upside down half-size pulses moving back towards each other. At time t = /c they recombine at the point - a to instantaneously form a full-sized, but upside-down mirror image of the original disturbance -- in accordance with (14.98). The recombined pulse in turn splits apart into two upside down half-size pulses that, when each collides with the end, reflects and returns to its original upright form. At time t = 2 /c, the pulses recombine to exactly reproduce the original displacement. The process then repeats, and the solution is periodic in time with period 2 /c. In Figure 14.13, the first picture displays the initial displacement. In the second, it has splits into left and right moving, half-size clones. In the third picture, the left moving bump is in the process of colliding with the left end of the bar. In the fourth picture, it has emerged from the collision, and is now upside down, reflected, and moving to the right. Meanwhile, the right moving pulse is starting to collide with the right end. In the fifth picture, both pulses have completed their collisions and are now moving back towards each other, where, in the last picture, they recombine into an upside-down mirror image of the original pulse. The process then repeats itself, in mirror image, finally recombining to the original pulse, at which point the entire process starts over. The Neumann (free) boundary value problem u (t, 0) = 0, x u (t, ) = 0, x (14.131)
is handled similarly. Since the solution has the form of a Fourier cosine series in x, we extend the initial conditions to be even, 2 periodic functions f (- x) = f (x), 2/25/07 f (x + 2 ) = f (x), 795 g(- x) = g(x), g(x + 2 ) = g(x).
c 2006 Peter J. Olver
The resulting d'Alembert solution is also even and 2 periodic in x, and hence satisfies the boundary conditions. In this case, when a pulse hits one of the ends, its reflection remains upright, but becomes a mirror image of the original; a familiar physical illustration is a water wave that reflects off a solid wall. Further details are left to the reader in Exercise In summary, we have now learned two different ways to solve the one-dimensional wave equation. The first, based on Fourier analysis, emphasizes the vibrational or wave character of the solutions, while the second, based on the d'Alembert formula, emphasizes their particle aspects, where individual wave packets collide with each other, or reflect at the boundary, all the while maintaining their overall form. Some solutions look like vibrating waves, while others seem much more like interacting particles. But, like the blind men describing the elephant, these are merely two facets of the same solution. The Fourier series formula shows how every particle-like solution can be decomposed into its constituent vibrational modes, while the d'Alembert formula demonstrates how vibrating waves combine as moving particles. The coexistence of particle and wave features is reminiscent of the long running historical debate over the nature of light. In the beginning, Newton and his disciples advocated a particle basis, in the form of photons. However, until the beginning of the twentieth century, most physicists advocated a wave or vibrational viewpoint. Einstein's explanation of the photoelectric effect in 1905 served to resurrect the particle interpretation. Only with the establishment of quantum mechanics was the debate resolved -- light, and, indeed, all subatomic particles are both, manifesting particle and wave features, depending upon the experiment and the physical situation. But the theoretical evidence for wave-particle duality already existed in the competing solution formulae of the classical wave equation!
14.6. Numerical Methods.
As you know, most differential equations are much too complicated to be solved analytically. Thus, to obtain quantitative results, one is forced to construct a sufficiently accurate numerical approximation to the solution. Even in cases, such as the heat and wave equations, where explicit solution formulas (either closed form or infinite series) exist, numerical methods still can be profitably employed. Moreover, justification of the numerical algorithm is facilitated by the ability compare it with an exact solution. Moreover, the lessons learned in the design of numerical algorithms for "solved" problems prove to be of immense value when one is confronted with more complicated problems for which solution formulas no longer exist. In this final section, we present some of the most basic numerical solution techniques for the heat and wave equations. 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., [28, 35, 109]. Finite Differences Numerical solution methods for differential equations can be partitioned into two principal classes. (In this oversimplified presentation, we are ignoring more specialized methods of less general applicability.) The first category, already introduced in Section 11.6, are the 2/25/07 796
c 2006 Peter J. Olver
One-Sided Difference Figure 14.14.
Central Difference Finite Difference Approximations.
finite element methods. Finite elements are designed for the differential equations describing equilibrium configurations, since they rely on minimizing a functional. A competitive alternative is to directly approximate the derivatives appearing in the differential equation, which requires knowing appropriate numerical differentiation formulae. 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 (14.132)
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 14.14. 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 1 u(x + h) = u(x) + u (x) h + 2 u () h2 . (14.133) We have used the Cauchy form (C.8) 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 2/25/07 797
1 2
u () h.
c 2006
(14.134)
Peter J. Olver
Since the error is proportional to h, we say that the finite difference quotient (14.134) 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). (14.135) 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 14.11. Let u(x) = sin x. Let us try to approximate u (1) = cos 1 = 0.5403023 . . . by computing finite difference quotients sin(1 + h) - sin 1 . h The result for different values of h is listed in the following table. cos 1 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 h2 h3 u(x + h) = u(x) + u (x) h + u (x) + u (x) + O(h4 ), 2 6 h2 h3 u(x - h) = u(x) - u (x) h + u (x) - u (x) + O(h4 ), 2 6
(14.136)
where the error terms are proportional to h4 . Adding the two formulae together gives 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 ), (14.137) 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) = 2/25/07 798
c 2006 Peter J. Olver
Example 14.12. 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 (14.137): 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
1 Each reduction in step size by a factor of 10 reduces the size of the error by a factor of 1 100 and results in a gain of two new decimal digits of accuracy, confirming that the finite difference approximation is of second order. 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 14.11, 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 (14.135) 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 (14.136). To solve for u (x), we subtract the two formulae, and so h3 + O(h4 ). 3 Rearranging the terms, we are led to the well-known centered difference formula u(x + h) - u(x - h) = 2 u (x) h + u (x) u(x + h) - u(x - h) + O(h2 ), (14.138) 2h 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 u (x) =
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 .
2/25/07
799
c 2006
Peter J. Olver
x - h, u(x - h) and x + h, u(x + h) on the graph of u centered symmetrically about the point x. Figure 14.14 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 14.13. Return to the function u(x) = sin x considered in Example 14.11. The centered difference approximation to its derivative u (1) = cos 1 = 0.5403023 . . . is cos 1 The results are tabulated as follows:
h approximation error .1 0.53940225217 -0.00090005370 .01 0.54029330087 -0.00000900499 .001 0.54030221582 -0.00000009005 .0001 0.54030230497 -0.00000000090
sin(1 + h) - sin(1 - h) . 2h
As advertised, the results are much more accurate than the one-sided finite difference approximation used in Example 14.11 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. Numerical Algorithms for the Heat Equation Consider the heat equation 2u u = , t x2 0 < x < , t 0, (14.139)
representing a bar of length and thermal diffusivity > 0, which is assumed to be constant. To be concrete, we impose time-dependent Dirichlet boundary conditions u(t, 0) = (t), u(t, ) = (t), t 0, (14.140)
specifying the temperature at the ends of the bar, along with the initial conditions u(0, x) = f (x), 0 x , (14.141)
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 = 2/25/07 800 and 0 = t 0 < t1 < t2 < .
c 2006 Peter J. Olver
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, (14.142)
to denote the numerical approximation to the solution value at the indicated mesh point. 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 (14.137), and hence u(ti , xj+1 ) - 2 u(ti , xj ) + u(ti , xj-1 ) 2u (ti , xj ) + O(h2 ) 2 x h2 ui,j+1 - 2 ui,j + ui,j-1 + O(h2 ), 2 h
(14.143)
where the error in the approximation is proportional to h2 . Similarly, the one-sided finite difference approximation (14.135) is used for the time derivative, and so u(ti+1 , xj ) - u(ti , xj ) ui+1,j - ui,j u (ti , xj ) + O(k) + O(k), t k k (14.144)
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 (14.144) by the second order central difference approximation (14.138). However, this produces significant complications, and the resulting numerical scheme is not practical. Replacing the derivatives in the heat equation (14.145) by their finite difference approximations (14.143), (14.144), 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 = k . h2 i = 0, 1, 2, . . . , j = 1, . . . , n - 1, (14.145)
(14.146)
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 . 2/25/07 801
c 2006 Peter J. Olver
The initial condition (14.141) means that we should initialize our numerical data by sampling the initial temperature at the mesh points: u0,j = fj = f (xj ), j = 1, . . . , n - 1. i = 0, 1, 2, . . . . (14.147)
Similarly, the boundary conditions (14.140) require that ui,0 = i = (ti ), ui,n = i = (ti ), (14.148)
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 , fn = f () = u(0, ) = (0) = 0 .
The three equations (14.145148) completely prescribe the numerical approximation algorithm for solving the initial-boundary value problem (14.139141). 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
(14.149)
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 (14.140). Then (14.145) assumes the compact vectorial form u(i+1) = A u(i) + b(i) , (14.150) where 1 - 2 1 - 2 1 - 2 .. . A= . .. .. . . .. , i 0 0 . . = . . 0 i
b(i)
(14.151)
1 - 2
The coefficient matrix A is symmetric and tridiagonal. The contributions (14.148) 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 14.14. Let us fix the diffusivity = 1 and the bar length = 1. For illustrative purposes, we take a spatial step size of h = .1. In Figure 14.15 we compare two (slightly) different time step sizes on the same initial data as used in (14.22). 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 14.1. 2/25/07 802
c 2006 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 14.15.
Numerical Solutions for the Heat Equation Based on the Explicit Scheme.
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, . . . , (14.152) (14.153)
and so (14.150) reduces to a homogeneous, linear iterative system u(i+1) = A u(i) .
According to Proposition 10.11, 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. Example 14.14 is indicating that for mesh size h = .1, the time step 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 10.14, 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, which solves Exercise 8.2.48. It is a direct consequence of Exercise 8.2.47, which contains the explicit formulae for the eigenvectors. Lemma 14.15. 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 (14.154) k = a + 2 b cos , k = 1, . . . , n - 1. n In our particular case, a = 1 - 2 and b = , and hence the eigenvalues of the matrix A given by (14.151) are k = 1 - 2 + 2 cos 2/25/07 k , n k = 1, . . . , n - 1.
c 2006 Peter J. Olver
803
Since the cosine term ranges between -1 and +1, the eigenvalues satisfy 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 , h2 2 or k h2 . 2 (14.155)
With some additional analytical work, [109], it can be shown that this is sufficient to conclude that the numerical scheme (14.145148) converges to the true solution to the initial-boundary 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 (14.155) 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 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 (14.156)
to approximate the temporal derivative. Substituting (14.156) and the same approximation (14.143) 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 , (14.157) j = 1, . . . , n - 1, where the parameter = k/h2 is as above. The initial and boundary conditions also have the same form (14.147, 148). The system can be written in the matrix form A u(i+1) = u(i) + b(i+1) , (14.158)
where A is obtained from the matrix A in (14.151) 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 1.7, 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 (14.152), the system takes the form u(i+1) = A-1 u(i) , 2/25/07 804
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 14.16.
Numerical Solutions for the Heat Equation Based on the Implicit Scheme.
and the convergence is now governed by the eigenvalues of A-1 . Lemma 14.15 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 (14.158) is convergent for any choice of step sizes h, k, and hence unconditionally stable. Example 14.16. Consider the same initial-boundary value problem considered in Example 14.14. In Figure 14.16, 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 (14.159)
which can be obtained by averaging the explicit and implicit schemes (14.145, 157). We can write the iterative system in matrix form B u(i+1) = C u(i) + 2/25/07 805
1 2
b(i) + b(i+1) ,
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 14.17.
Numerical Solutions for the Heat Equation Based on the CrankNicolson Scheme.
where 1+ -1 2 1 1 -2 1+ -2 .. B= . -1 2 .. . . . , . .. .
1 1- 2 1 2 1- C= 1 2
1 2
. .
.. ..
Convergence is governed by the generalized eigenvalues of the tridiagonal matrix pair B, C, or, equivalently, the eigenvalues of the product B -1 C, cf. Exercise 9.5.33. According to Exercise , these are 1- 1+ k n k 1 - cos n 1 - cos
. . . . .. .
(14.160)
k =
,
k = 1, . . . , n - 1.
(14.161)
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 14.17 shows the result of running the scheme on the initial data (14.22). 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 2/25/07 806
c 2006 Peter J. Olver
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. 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 (14.124), 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 (14.162) = c2 , 0 < x < , t 0, t2 x2 modeling vibrations of a homogeneous bar of length with constant wave speed c > 0. We 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 0 x . (14.164) u(t, ) = (t), t 0. (14.163)
. n In order to discretize the wave equation, we replace the second order derivatives by their standard finite difference approximations (14.137), namely h= u(ti+1 , xj ) - 2 u(ti , xj ) + u(ti-1 , xj ) 2u (ti , xj ) + O(h2 ), 2 2 t k u(ti , xj+1 ) - 2 u(ti , xj ) + u(ti , xj-1 ) 2u (t , x ) + O(k 2 ), x2 i j h2
(14.165)
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 (14.165) into the partial differential equation (14.162), 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, (14.166)
for the numerical approximations ui,j u(ti , xj ) to the solution values at the mesh points. The positive parameter ck > 0 (14.167) = h 2/25/07 807
c 2006 Peter J. Olver
depends upon the wave speed and the ratio of space and time step sizes. The boundary conditions (14.163) require that ui,0 = i = (ti ), ui,n = i = (ti ), u(i+1) = B u(i) - u(i-1) + b(i) , where 2 (1 - 2 ) 2 2 2 (1 - 2 ) 2 .. . B= . .. 2 .. .. . . i = 0, 1, 2, . . . . (14.168)
This allows us to rewrite the system in matrix form (14.169)
2 ui,1 i ui,2 0 . . (i) , u = . , b(i) = . . . . 0 u 2 i,n-2 2 i uii,n-1 2 2 (1 - 2 ) (14.170) The entries of u(i) are, as in (14.149), the numerical approximations to the solution values at the interior nodes. Note that the system (14.169) 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 u(k, xj ) - u(0, xj ) u1,j - fj u (0, xj ) t k k to compute the required values u1,j = fj + k gj . gj = (14.171)
However, the approximation (14.171) 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 ) u k 2u u c2 k 2 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 - 2 fj + fj-1 ) , = f (xj ) + k g(xj ) + 2 2 h2 j+1 808
c 2006 Peter J. Olver
2/25/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 14.18.
Numerically Stable Waves.
where we can use the finite difference approximation (14.137) 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 B u(0) + k g + 2 b(0) , 1 2
2 fj+1 + (1 - 2 )fj +
1 2
2 fj-1 + k gj ,
(14.172)
(14.173)
we will have maintained the desired order k 2 (and h2 ) accuracy. Example 14.17. 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 14.18 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 14.19 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 (14.169) into a first order system. In analogy with u(i) Example 10.6, this is accomplished by introducing the vector z(i) = R 2n-2 . u(i-1) Then z(i+1) = C z(i) + c(i) , 2/25/07 where 809 C= B I -I O . (14.174)
Peter J. Olver
c 2006
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 14.19.
Numerically Unstable Waves.
Therefore, the stability of the method will be determined by the eigenvalues of the coeffiu , can be written out cient matrix C. The eigenvector equation C z = z, where z = v in its individual components: B u - v = u, ( B - 2 - 1) v = 0, or 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 14.15, in which a = 2(1 - 2 ) and b = 2 , and hence are Bv = + + 1 =2 1 - 2 + 2 cos k n , k = 1, . . . , n - 1. k < 1. n
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 (14.175)
Each pair of solutions to these n - 1 quadratic equations, namely = ak k a2 - 1 , k (14.176)
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 (14.176) are complex numbers of modulus 1, indicated stability (but not convergence) of the matrix C. Therefore, in view of (14.175), we should require that h ck < 1, or k< , (14.177) = h c 2/25/07 810
c 2006 Peter J. Olver
Figure 14.20.
The Courant Condition.
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 (14.177) 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 14.20, the wave speed is c = 1.25. The first figure has equal mesh spacing k = h, and does not satisfy the Courant condition (14.177), whereas the 1 second figure has k = 2 h, which does. Note how the characteristic lines starting at a given mesh point have progressed beyond the neighboring mesh points after one time step in the first case, but not in the second.
14.7. A General Framework for Dynamics.
According to Section 12.1, the one-dimensional heat and wave equations are specific instances of two broad classes of dynamical systems that include, in a common framework, both discrete dynamics modeled by systems of ordinary differential equations and continuum systems modeled by (systems of) partial differential equations. In preparation for their multi-dimensional generalizations, it will be useful to summarize the general mathematical framework, which can be regarded as the dynamical counterpart to the framework for equilibrium developed in Sections 7.5 and 15.4, which you are advised to review before going through this section. Readers more attuned to concrete examples, on the other hand, might prefer to skip this material entirely, referring back as necessary. In all situations, the starting point is a linear function that maps a vector space U to another vector space V . In mechanics, the elements of U represent displacements, while the elements of V represent strains (elongations). In 2/25/07 811
c 2006 Peter J. Olver
L : U - V
(14.178)
electromagnetism and gravitation, elements of U represent potentials and elements of V electric or magnetic or gravitational fields. In thermodynamics, elements of U represent temperature distributions, and elements of V temperature gradients. In fluid mechanics, U contains potential functions and V is fluid velocities. And so on. In the discrete, finite-dimensional setting, when U = R n and V = R m , the linear function L[ u ] = A u is represented by multiplication by an m n matrix A -- the incidence matrix and its generalizations. In the infinite-dimensional function space situation presented in Chapter 11 and earlier in this chapter, the linear map L is a differential operator; in the case of one-dimensional elastic bars, it is the first order derivative L = Dx , while for 2 beams it is a second order derivative L = Dx . In the multi-dimensional physics treated in subsequent chapters, the most important example is the gradient operator, L = , that maps scalar potentials to vector fields. The vector spaces U and V are further each endowed with inner products, that encapsulate the consitutive assumptions underlying the physical problem. Definition 7.53 explains how to construct the resulting adjoint map L : V - U (14.179)
which goess in the reverse direction. In finite dimensions, when U = R n and V = R m are both equipped with the Euclidean dot product, the adjoint corresponds to the simple transpose operation. Modifications under more general weighted inner products were discussed in Section 7.5. In infinite-dimensional contexts, the computation of the adjoint presented in Section 11.3 relied on an integration by parts argument, supplemented by suitable boundary conditions. In the next chapter, we will adapt this argument to determine adjoints of multi-dimensional differential operators. The crucial operator underlying the equilibrium and dynamical equations for a remarkably broad range of physics is the self-adjoint combination K = L L: U - U. (14.180)
According to Theorem 7.60, the operator K is positive semi-definite in all situations, and positive definite if and only if ker L = {0}. In the finite-dimensional context, K is represented by the symmetric positive (semi-)definite Gram matrix AT A when both U = R n and V = R m have the dot product, by the symmetric combination AT C A when V has a weighted inner product represented by the positive definite matrix C > 0, and the more general self-adjoint form M -1 AT C A when U also is endowed with a weighted inner product represented by M > 0. In one-dimensional bars, K is represented by a self-adjoint second order differential operator, whose form depends upon the underlying inner products., while in beams it becomes a fourth order differential operator. The definiteness of K (or lack thereof) depends upon the imposed boundary conditions. In higher dimensions, as we discuss in Section 15.4, K becomes the Laplacian or a more general elliptic partial differential operator. With this set-up, the basic equilibrium equation has the form K[ u ] = f, 2/25/07 812
c 2006
(14.181)
Peter J. Olver
where f represents an external forcing function. In finite dimensions, this is a linear system consisting of n equations in n unknowns with positive (semi-)definite coefficient matrix. In function space, it becomes a self-adjoint boundary value problem for the unknown function u. If K is positive definite, the solution is unique. (But rigorously proving the existence of a solution is not trivial, requiring serious analysis beyond the scope of this text.) Theorem 7.61 says that the solutions are characterized by a minimization principle, as the minimizers of the quadratic function(al) p[ u ] =
1 2
L[ u ]
2
- u,f ,
(14.182)
which typically represents the potential energy in the system. If K is only positive semidefinite, existence of a solution requires that the forcing function satisfy the Fredholm condition that it be orthogonal to the the unstable modes, that is the elements of ker K = ker L With the equilibrium operator K in hand, there are two basic types of dynamical systems of importance in physical models. Unforced diffusion processes are modeled by a dynamical system of the form ut = - K[ u ], where K = L L (14.183)
is the standard self-adjoint combination of a linear operator L. In the discrete case, K is a matrix and so this represents a first order system of ordinary differential equations, that has the form of a linear gradient flow (9.18), so named because it decreases the energy function q[ u ] = L[ u ] 2 , (14.184) which is (14.182) when f = 0, as rapidly as possible. In the continuous case, K is a differential operator, and (14.183) represents a partial differential equation for the timevarying function u = u(t, x). The solution to the general diffusion equation (14.183) mimics earlier our separation of variables method for the one-dimensional heat equation, as well as the original solution technique for linear systems of ordinary differential equations, as developed in Chapter 9. The separable solutions are of exponential form u(t) = e- t v, (14.185)
where v U is a fixed element of the domain space -- i.e., a vector in the discrete context, or a function v = v(x) that only depends on the spatial variables in the continuum versions. Since the operator K is linear and does not involve t differentiation, we find u = - e- t v, while K[ u ] = e- t K[ v ]. t Substituting back into (14.183) and canceling the common exponential factors, we are led to the eigenvalue problem K[ v ] = v. (14.186) Thus, v must be an eigenvector/eigenfunction for the linear operator K, with the corresponding eigenvalue. 2/25/07 813
c 2006 Peter J. Olver
Generalizing our earlier observations on the eigenvalues of positive definite matrices and the boundary value problems associated with the one-dimensional heat equation, let us establish the positivity of eigenvalues of such general self-adjoint, positive (semi-)definite linear operators. Theorem 14.18. All eigenvalues of the linear operator K = L L are real and nonnegative: 0. If, moreover, K is positive definite, or, equivalently, ker K = ker L = 0, then all eigenvalues are strictly positive: > 0. Proof : Suppose K[ u ] = u with u = 0. Then u
2
= u , u = K[ u ] , u = L L[ u ] , u = L[ u ] , L[ u ] = L[ u ]
2
0,
by the defining equation (7.74) of the adjoint operator. Since u 2 > 0, this immediately implies that 0. Furthermore, in the positive definite case ker L = {0}, and so L[ u ] = 0. Thus, L[ u ] 2 > 0, proving that > 0. Q.E.D. We index the eigenvalues in increasing order: 0 < 1 2 3 . (14.187)
where the eigenvalues are repeated according to their multiplicities, and 0 = 0 is an eigenvalue only in the positive semi-definite case. Each eigenvalue induces a separable eigensolution uk (t) = e- k t vk (14.188) to the diffusion equation (14.183). Those associated with strictly positive definite eigenvalues are exponentially decaying, at a rate equal to the eigenvalue k > 0, while any null eigenvalue modes correspond to constant solutions u0 (t) v0 where v0 is any null eigenvector/eigenfunction, i.e., element of ker K = ker L. The general solution is built up as a linear combination u(t) =
k
ck uk (t) =
k
ck e- k t vk
(14.189)
of the eigensolutions. Thus, all solutions decay exponentially fast as t to an equilibrium solution, which is 0 in the positive definite case, or, more generally, the null eigenspace component c0 v0 . The overall rate of decay is, generally, prescribed by the smallest positive eigenvalue 1 > 0. In the discrete version, the summation (14.189) has only finitely many terms, corresponding to the n eigenvalues of the matrix representing K. Moreover, thanks to the Spectral Theorem 8.26, the eigensolutions form a basis for the solution space to the diffusion equation. In infinite-dimensional function space, there are, in many instances, an infinite number of eigenvalues, with k as k . Completeness of the resulting eigensolutions is a more delicate issue. Often, as in the one-dimensional heat equation, the eigenfunctions are complete and the (Fourier) series (14.189) converges and represents the general solution. But in situations involving unbounded domains, like the hydrogen atom to be discussed in Section 18.7, there are additional separable solutions corresponding to the so-called continuous spectrum that are not represented in terms of eigenfunctions, 2/25/07 814
c 2006 Peter J. Olver
and the analysis is considerably more involved, requiring analogs of the Fourier transform. A full discussion of completeness and convergence of eigenfunction expansions must be relegated to an advanced course in analysis, [47, 151]. Assuming completeness, the eigensolution coefficients ck in (14.189) are prescribed by the initial conditions, which require ck vk = h,
k
(14.190)
where h = u(0) represents the initial data. To compute the coefficients ck in the eigenfunction expansion (14.190), we appeal, as in the case of ordinary Fourier series, to orthogonality of the eigenvectors/eigenfunctions vk . Orthogonality is proved by a straightforward adaptation of our earlier proof of part (b) of Theorem 8.20, guaranteeing the orthogonality of the eigenvectors of a symmetric matrix. Theorem 14.19. Two eigenvectors/eigenfunctions u, v of the self-adjoint linear operator K = L L that are associated with distinct eigenvalues = are orthogonal. Proof : Self-adjointness of the operator K means that K[ u ] , v = L[ u ] , L[ v ] = u , K[ v ] for any for any u, v. In particular, if K[ u ] = u, K[ v ] = v, are eigenfunctions, then u , v = K[ u ] , v = u , K[ v ] = u , v , which, assuming = , immediately implies orthogonality: u , v = 0. Q.E.D.
If an eigenvalue admits more than one independent eigenvector/eigenfunction, we can apply the GramSchmidt process to produce an orthogonal basis of its eigenspace V = ker(K - I ). In this way, the entire set of eigenvectors/eigenfunctions can be assumed to be mutually orthogonal: vj , vk = 0, j = k.
As a consequence, taking the inner product of both sides of (14.190) with the eigenfunction vk leads to the equation ck vk
2
= h , vk
and hence
ck =
h , vk . vk 2
(14.191)
In this manner we recover our standard orthogonality formula (5.7) for expressing elements of a vector space in terms of an orthogonal basis. The second important class of dynamical systems consists of second order (in time) vibration equations utt = - K[ u ]. (14.192)
This assumes that the eigenspace V is finite-dimensional -- which is assured whenever the eigenvalues k as k .
2/25/07
815
c 2006
Peter J. Olver
Newton's equations of motion in the absence of non-conservative frictional forces, the propagation of waves in fluids and solids, as well as electromagnetic waves, and many other related physical problems lead to such vibrational systems. In this case, the separable solutions are of trigonometric form u(t, x) = cos( t) v(x) or sin( t) v(x). (14.193)
Substituting this ansatz back into the vibration equation (14.192) results in the same eigenvalue problem (14.186) with eigenvalue = 2 equal to the square of the vibrational frequency. We conclude that the normal mode or separable eigensolutions take the form uk (t) = cos(k t) vk , uk (t) = sin(k t) vk , provided
2 k = k > 0
is a non-zero eigenvalue. In the stable, positive definite case, there are no zero eigenvalues, and so the general solution is built up as a (quasi-)periodic combination u(t) =
k
ck uk (t) + dk uk (t) =
k
rk cos(k t + k ) vk ,
(14.194)
of the eigenmodes. The initial conditions g = u(0) =
k
ck vk ,
h = ut (0) =
k
dk k vk ,
(14.195)
are used to specify the coefficients ck , dk , using the same orthogonality formula (14.191): ck = f , vk , vk 2 dk = f , vk k vk
2
.
(14.196)
In the unstable, positive semi-definite cases, the null eigensolutions have the form u0 (t) = v0 , u0 (t) = t v0 ,
where v0 ker K = ker L, and must be appended to the series solution. The unstable mode u0 (t) is excited if and only if the initial velocity is not orthogonal to the kernel element: h , v0 = 0. In classical mechanics, the diffusion and vibration equations and their variants (see, for instance, Exercises for versions with external forcing and Exercise for vibrations with frictional effects) are the most important classes of dynamical systems. In quantum mechanics, the basic dynamical system is known as the Schrdinger equation, first written o down by the the German physicist Erwin Schrdinger, one of the founders of quantum o mechanics. The abstract form of the Schrdinger equation is o i ut = K[ u ].
(14.197)
The solution is periodic if and only if the frequencies appearing in the sum are all integer multiples of a common frequency: k = nk for nk N.
2/25/07
816
c 2006
Peter J. Olver
Here i =
-1, while =
h 1.055 10-34 Joule seconds (14.198) 2 is Planck's constant, whose value governs the quantization of all physical quantities. At each time t, the solution u(t, x) to the Schrdinger equation represents the wave function o of the quantum system, and so is a complex-valued square integrable function of constant L2 norm: u = 1. (See Sections 12.5 and 13.3 for the basics of quantum mechanics and Hilbert space.) As usual, we interpret the wave fuction as a probability density on the possible quantum states, and so the Schrdinger equation governs the dynamical evolution o of quantum probabilities. The operator K = L L is known as the Hamiltonian for the quantum mechanical system governed by (14.197), and, typically represents the quantum energy operator. For physical systems such as atoms and nuclei, the relevant Hamiltonian operator is constructed from the classical energy through the rather mysterious process of "quantization". The interested reader should consult a basic text on quantum mechanics, e.g., [124, 129], for full details on both the physics and underlying mathematics. Proposition 14.20. If u(t) is a solution to the Schrdinger equation, its Hermitian o L norm u is constant.
2
Proof : Since the solution is complex-valued, we use the sesquilinearity of the underlying Hermitian inner product, as in (3.92), to compute d u dt
2
= ut , u + u , u t = - i K[ u ] , u + u, - i K[ u ] =- i K[ u ] , u + i u , K[ u ] = 0,
which vanishes since K is self-adjoint. Since its derivative vanishes everywhere, this implies that u 2 is constant. Q.E.D. As a result, if the initial data u(t0 ) = h is a quantum mechanical wave function, meaning that h = 1, then, at each time t, the solution to the Schrdinger equation also o has norm 1, and hence remains a wave function for all t. Apart from the extra factor of i , the Schrdinger equation looks like a diffusion o equation (14.183). (Warning: Despite this superficial similarity, their solutions have radically different behavior.) This inspires us to seek separable solutions with an exponential ansatz: u(t, x) = e t v(x). Substituting this expression into the Schrdinger equation (14.197) and canceling the como mon exponential factors reduces us to the usual eigenvalue problem K[ v ] = v, with eigenvalue = - i .
Let vk (x) denote the normalized eigenfunction associated with the k th eigenvalue k . Th corresponding separable solution of the Schrdinger equation is the complex wave functions o uk (t, x) = e i k t/ vk (x). 2/25/07 817
c 2006 Peter J. Olver
Observe that, in contrast to the exponentially decaying solutions to the diffusion equation, the eigensolutions to the Schrdinger equation are periodic, of frequencies proportional o to the eigenvalues: k = k / . (Along with constant solutions corresponding to the null eigenmodes, if any.) The general solution is a (quasi-)periodic series in the fundamental eigensolutions. The periodicity of the summands has the additional implication that, again unlike the diffusion equation, the Schrdinger equation can be run backwards in time. So, o we can figure out both the past and future behavior of a quantum system from its present configuration. Example 14.21. In a single space dimension, the simplest version of the Schrdinger o equation is based on the derivative operator L = Dx , for which, assuming appropriate 2 boundary conditions, the self-adjoint combination K = L L = - Dx . Thus, (14.197) reduces to the second order partial differential equation i ut = - uxx . Imposing the Dirichlet boundary conditions u(t, 0) = u(t, ) = 0, the Schrdinger equation governs the dynamics of a quantum particle confined to the o interval 0 < x < ; the boundary conditions imply that there is zero probability of the particle escaping from the interval. According to Section 14.1, the eigenfunctions of the Dirichlet eigenvalue problem vxx + v = 0, are vk (x) = k 2 sin x for k = 1, 2, . . . , with eigenvalue k = k2 2 , 2 v(0) = v() = 0, (14.199)
where the initial factor is ensures that vk has unit L2 norm, and hence is a bona fide wave function. The corresponding separable solutions or eigenmodes are uk (t, x) = 2 exp i k2 2 t 2 sin k x. (14.200)
The eigenvalues represent the energy levels of the particle, which can be observed from the spectral lines emitted by the system. For instance, when an electron jumps from one level to another, it conserves energy by emitting a photon, with energy equal to the difference in energy between the two quantum levels. These emitted photons define the observed electromagnetic spectral lines, hence the adoption of the physics term "spectrum" to describe the eigenvalues of the Hamiltonian operator K.
2/25/07
818
c 2006
Peter J. Olver