4 Chapter Separation of Variables
There are three paradigmatic linear second order partial differential equations that have collectively driven the development of the entire subject. The first two we have already encountered: The wave equation describes vibrations and waves in continuous media, including sound waves, water waves, elastic waves, electromagnetic waves, and so on. The heat equation models diffusion processes, including thermal energy in solids, solutes in liquids, and biological populations. Third, and most important of all, is the Laplace equation and its inhomogeneous counterpart, the Poisson equation, that govern equilibrium mechanics. These two equations arise in an astonishing variety of mathematical and physical contexts, ranging through elasticity and solid mechanics, fluid mechanics, electromagnetism, potential theory, thermomechanics, geometry, probability, number theory, and many more. The solutions to the Laplace equation are known as harmonic functions, and the discovery of their many remarkable properties forms one of the most celebrated chapters in the history of mathematics. All three equations, along with their multi-dimensional kin, appear repeatedly throughout this text. The aim of the current chapter is to develop the method of separation of variables for solving these key partial differential equations in their two independent variable incarnations. For the wave and heat equations, the variables are time, t, and a single space coordinate, x, leading to initial-boundary value problems modeling the dynamical behavior of a one-dimensional medium. For the Laplace and Poisson equations, both variables represent space coordinates, x and y, and the associated boundary value problems model the equilibrium configuration of a planar body, e.g., the deformations of a membrane. Separation of variables seeks special solutions that can be written as the product of functions of the individual variables, thereby reducing the partial differential equation to a pair of ordinary differential equations. More general solutions can then be expressed as infinite series in the appropriate separable solutions. For the two-variable equations considered here, this results in a Fourier series representation of the solution. In the case of the wave equation, separation of variables serves to focus attention on the vibrational character of the solution, whereas the earlier d'Alembert approach emphasizes its particle-like aspects. Unfortunately, for the Laplace equation, separation of variables only applies to boundary value problems in very special geometries, e.g., rectangles and disks. Further development of the separation of variables method for solving partial differential equations in three or more variables can be found in Chapters 12 and 13. In the final section, we take the opportunity to summarize the fundamental tripartite classification of planar second order partial differential equations: hyperbolic, such as the wave equation; parabolic, such as the heat equation; and elliptic, such as the Laplace and 1/19/12 100
c 2012 Peter J. Olver
Poisson equations. Each category enjoys a number of distinctive properties and features, both analytical and numerical, and, in effect, forms its own mathematical subdiscipline.
4.1. The Diffusion and Heat Equations.
Let us begin with a brief physical derivation of the heat equation from first principles. We consider a bar -- meaning a thin, heat-conducting body. "Thin" means that we can regard the bar as a one-dimensional continuum with no significant transverse temperature variation. We will assume that the bar is fully insulated along its length, and so heat can only enter (or leave) through its uninsulated endpoints. We use t to represent time, and a x b to denote spatial position along the bar, which occupies the interval [ a, b ]. 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 principles. First is the Law of Conservation of Heat Energy. Recalling the general Definition 2.7, this particular conservation law takes the form w + = 0, t x (4.1)
in which (t, x) represents the thermal energy density at time t and position x, while w(t, x) denotes the heat flux , i.e., the rate of flow of thermal energy along the bar. Our sign convention is that w(t, x) > 0 at points where heat energy flows in the direction of increasing x (left to right). The integrated form (2.48) of the conservation law, namely d dt
b a
(t, x) dx = w(t, a) - w(t, b),
(4.2)
states that the rate of change in the thermal energy within the bar is equal to the total heat flux passing through its uninsulated ends. The signs of the boundary terms confirm that heat flux into the bar results in an increase in temperature. The second ingredient is a constitutive assumption concerning the bar's material properties. It has been observed that, under reasonable conditions, thermal energy is proportional to temperature: (t, x) = (x) u(t, x). (4.3) 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 degree. Note that we are assuming the medium 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 its material properties do not depend upon the temperature; otherwise, we would be forced to deal with a much thornier nonlinear diffusion equation. The third physical principle relates heat flux and temperature. Physical experiments show that the heat energy moves from hot to cold at a rate that is in direct proportion to 1/19/12 101
c 2012 Peter J. Olver
the temperature gradient which, in the one-dimensional case, means its derivative u/x. The resulting relation u (4.4) w(t, x) = - (x) x is known as Fourier's Law of Cooling. The proportionality factor (x) > 0 is the thermal conductivity of the bar at position x, and the minus sign reflects the everyday observation that heat energy moves from hot to cold. A good heat conductor, e.g., silver, will have high conductivity, while a poor conductor, e.g., glass, will have low conductivity. Combining the three laws (4.1, 3, 4) produces the linear diffusion equation (x) u = t x (x) u x , a < x < b, (4.5)
governing the thermodynamics of a one-dimensional medium. It is also used to model a wide variety of diffusive processes, including chemical diffusion, diffusion of contaminants in liquids and gases, population dispersion, and the spread of infectious diseases. If there is an external heat source along the length of the bar, then the diffusion equation acquires an additional inhomogeneous term: (x) u = t x (x) u x + h(t, x), a < x < b. (4.6)
In order to uniquely prescribe the solution u(t, x), we need to specify an initial temperature distribution u(t0 , x) = f (x), a x b. (4.7) In addition, we must impose a suitable boundary condition at each end of the bar. There are three common types. The first is a Dirichlet boundary condition, where the end is held at prescribed temperature. For example, u(t, a) = (t)
(4.8)
fixes the temperature (possibly time-varying) at the left end. Alternatively, the Neumann boundary condition u (4.9) (t, a) = (t) x prescribes the heat flux w(t, a) = - (a)ux (t, a) there. In particular, a homogeneous Neumann condition, ux (t, a) 0, models an insulated end that prevents heat energy flowing in or out. The Robin boundary condition u (4.10) (t, a) + k u(t, a) = (t), x with k > 0, models the heat exchange resulting from the end of the bar being placed in a reservoir at temperature (t). Each end of the bar is required to satisfy one 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, a) = (t), 1/19/12 102 u(t, b) = (t),
c 2012
(4.11)
Peter J. Olver
whereas a bar with two insulated ends requires two homogeneous Neumann boundary conditions u u (4.12) (t, a) = 0, (t, b) = 0. x x Mixed boundary conditions, with one end at a fixed temperature and the other insulated, are similarly formulated, e.g., u(t, a) = (t), Finally, the periodic boundary conditions u(t, a) = u(t, b), u u (t, a) = (t, b), x x (4.14) u (t, b) = 0. x (4.13)
correspond to a circular ring obtained by joining the two ends of the bar. As before, we are assuming the heat is only allowed to flow around the ring -- insulation prevents the radiation of heat from one side of the ring affecting the other side. The Heat Equation In this book, we will retain the term "heat equation" to refer to the case in which the bar is composed of a uniform material, and so its density , conductivity , and specific heat are all positive constants. We also exclude external heat sources, meaning that the bar remains insulated along its entire length. Under these assumptions, the general diffusion equation (4.5) reduces to the homogeneous heat equation 2u u = t x2 for the temperature u(t, x) at time t and position x. The constant = = (4.16) (4.15)
is called the thermal diffusivity, and incorporates all of the bar's relevant physical properties. The solution u(t, x) will be uniquely prescribed once we specify initial conditions (4.7) and a suitable boundary condition at both of its endpoints. As we learned in Section 3.1, the separable solutions to the heat equation are based on the exponential ansatz u(t, x) = e- t v(x), (4.17) where v(x) depends only on the spatial variable. Functions of this form, which "separate" into a product of a function of t times a function of x, are known as separable solutions.
Anticipating the eventual signs of the eigenvalues, and to facilitate later discussions, we now include a minus sign in the exponential term.
1/19/12
103
c 2012
Peter J. Olver
Substituting (4.17) into (4.15) and canceling the common exponential factors, we find that v(x) must solve the second order linear ordinary differential equation d2 v = v. dx2 Each nontrivial solution v(x) 0 is an eigenfunction, with eigenvalue , for the linear differential operator L[ v ] = - v (x). With the separable eigensolutions (4.17) in hand, we will then be able to reconstruct the desired solution u(t, x) as a linear combination, or, rather, infinite series thereof. Let us concentrate on the simplest case: a uniform, insulated bar of length that is held at zero temperature at both ends. We specify its initial temperature f (x) at time t0 = 0, and so the relevant initial and boundary conditions are - u(t, 0) = 0, u(t, ) = 0, u(0, x) = f (x), t 0, 0 x . (4.18)
The eigensolutions (4.17) are found by solving the Dirichlet boundary value problem d2 v (4.19) + v = 0, v(0) = 0, v() = 0. dx2 Repeating the analysis of Section 3.1, we find that if is either complex, or real and 0, then the only solution to the boundary value problem (4.19) is the trivial solution v(x) 0. Hence, all the eigenvalues must necessarily be real and positive. In fact, the reality and positivity of the eigenvalues does need not be explicitly checked, but, rather, follows from very general properties of positive definite boundary value problems, of which (4.19) is a particular case. See Section 10.5 for the underlying theory. When > 0, the general solution to the differential equation is a trigonometric function v(x) = a cos x + b sin x, where = / , and a and b are arbitrary constants. The first boundary condition requires v(0) = a = 0. Using this to eliminate the cosine term, the second boundary condition requires v() = b sin = 0. Therefore, since b = 0 -- as otherwise the solution is trivial and does not qualify as an eigenfunction -- must be an integer multiple of , and so 2 3 , , , ... . We conclude that the eigenvalues and eigenfunctions of the boundary value problem (4.19) are n x n 2 (4.20) , vn (x) = sin , n = 1, 2, 3, . . . . n = The corresponding eigensolutions (4.17) are = un (t, x) = exp - 1/19/12 n2 2 t 2 sin n x , n = 1, 2, 3, . . . .
c 2012
(4.21)
Peter J. Olver
104
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 4.1.
A Solution to the Heat Equation.
Each represents a trigonometrically oscillating temperature profile that maintains its form while decaying to zero at an exponentially fast rate. To solve the general initial value problem, we assemble the eigensolutions into an infinite series,
u(t, x) =
n=1
bn un (t, x) =
n=1
bn exp
n2 2 t - 2
sin
n x ,
(4.22)
whose coefficients bn are to be fixed by the initial conditions. Indeed, assuming that the series converges, the initial temperature profile is
u(0, x) =
n=1
bn sin
n x = f (x).
(4.23)
This has the form of a Fourier sine series (3.52) on the interval [ 0, ]. Thus, the coefficients are determined by the Fourier formulae (3.53), and so 2 n x (4.24) f (x) sin dx, n = 1, 2, 3, . . . . 0 As we will later prove, the resulting formula (4.22) describes the Fourier sine series for the temperature u(t, x) of the bar at each later time t 0. bn = Example 4.1. Consider the initial temperature profile 1 0 x 5, - x, 1 7 x 10 , x - 2, u(0, x) = f (x) = 5 5 7 1 - x, 10 x 1, b1 .045, b5 - .0081, 1/19/12 b2 - .096, b6 .0066, 105 b3 - .0145, b7 .0052,
(4.25)
on a bar of length 1, plotted in the first graph in Figure 4.1. Using (4.24), the first few Fourier coefficients of f (x) are computed (by either exact or numerical integration) to be b4 = 0, b8 = 0, ... .
Peter J. Olver
c 2012
The resulting Fourier series solution to the heat equation is
u(t, x) =
n=1
bn un (t, x) =
n=1 - 2 t
bn e- n
2
2 t
sin n x sin 2 x - .0145 e- 9
2
In Figure 4.1, the solution, for = 1, is plotted at the successive times t = 0., .02, .04, . . . , .1. Observe that the corners in the initial profile are immediately smoothed out. As time pro2 gresses, the solution decays, at a fast exponential rate of e- t e- 9.87 t , to a uniform, zero temperature, which is the equilibrium temperature distribution for the homogeneous Dirichlet boundary conditions. As the solution decays to thermal equilibrium, the higher Fourier modes rapidly disappear, and the solution assumes the progressively more symmetric shape of a single sine arc, of rapidly 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 (4.22) is progress of a sort. However, because we are unable to sum the series 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 integrable (e.g., piecewise continuous), then its Fourier coefficients are uniformly bounded; indeed, for any n 1, 2 | bn |
0
.045 e
sin x - .096 e- 4
2
t
t
sin 3 x - .
n x 2 f (x) sin dx
0
| f (x) | dx M.
(4.26)
This property holds even for quite irregular data. Under these conditions, each term in the series solution (4.22) 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, whereby any initial heat energy is eventually dissipated away through the ends. The small scale temperature fluctuations tend to rapidly cancel out through diffusion of heat energy, and the last term to disappear is the one with the slowest decay, namely u(t, x) b1 exp - 1/19/12 2 t 2 sin x , 106 where b1 = 1
f (x) sin x dx.
0 c 2012
(4.27)
Peter J. Olver
For generic initial data, the coefficient b1 = 0, and the solution approaches thermal equilibrium at an exponential rate prescribed by 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 its entire length. Also, again provided b1 = 0, the asymptotic shape of the temperature profile is a small, exponentially decaying sine arc, just as we observed in Example 4.1. In exceptional situations, namely when b1 = 0, the solution decays even faster, at a rate equal to the eigenvalue k = k 2 2 /2 corresponding to the first nonzero term, bk = 0, in the Fourier series; its asymptotic shape now oscillates k times over the interval. Another, closely related observation is that, for any fixed time t > 0 after the initial moment, the coefficients in the Fourier sine series (4.22) decay exponentially fast as n . According to the discussion at the end of Section 3.3, this implies that the Fourier series converges to an 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, which we state for a general initial time t0 . Theorem 4.2. If u(t, x) is a solution to the heat equation with piecewise continuous initial data f (x) = u(t0 , x), or, more generally, initial data satisfying (4.26), then, for any t > t0 , the solution u(t, x) is an infinitely differentiable function of x. In other words, the heat equation instantaneously smoothes out any discontinuities and corners in the initial temperature profile by fast damping of the high frequency modes. The heat equation's effect on irregular initial data 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 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 4.2 shows the effect running the heat equation , with = 1, to times t = 0., .00001, .00005, .0001, .001, .01 on a signal that has been contaminated by random noise. 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. Basic numerical solution schemes for the heat equation are to be discussed in Chapter 5. An important theoretical consequence of the smoothing property is that diffusion is a one-way process -- one cannot run time backwards and accurately infer what a temperature distribution looked like in the past. In particular, if the initial data u(0, x) = f (x) is not
To avoid artifacts at the ends of the interval, we are, in fact, using the periodic boundary conditions in the plots, as discussed below. Away from the ends, running the equation with Dirichlet boundary conditions leads to almost identical results.
1/19/12
107
c 2012
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 4.2.
Denoising a Signal with the Heat Equation.
smooth, then the value of u(t, x) for any t < 0 cannot be defined, because, if u(t0 , x) were defined and integrable at some t0 < 0, then, by Theorem 4.2, u(t, x) would be smooth at all subsequent times t > t0 , including t = 0, in contradiction to our assumption. Moreover, for most initial data, the Fourier coefficients in the solution formula (4.22) are, at any t < 0, exponentially growing as n , indicating that high frequency noise has completely overwhelmed the solution, thereby precluding any kind of convergence of the Fourier series. Mathematically, we can reverse future and past by changing t to - t. In the differential equation, this merely reverses the sign of the time derivative term; the x derivatives are unaffected. Thus, by the above reasoning, the backwards heat equation 2u u = - , t x2 with a negative diffusion coefficient - < 0,
is an ill-posed problem in the sense that, for most initial data, the solution is not defined in forwards time t > 0 (although it is well-posed if we run t backwards). The same holds for more general diffusion processes, e.g., (4.5). If, as in all physically relevant cases, the coefficient of uxx is everywhere positive, then the initial value problem is well-posed for t > 0, but ill-posed for t < 0. On the other hand, if the coefficient is everywhere negative the reverse holds. A coefficient that changes signs would cause the differential equation to be ill-posed in both directions. While theoretically undesirable, the unsmoothing effect of the backwards heat equation does have potential benefits in certain contexts. For example, in image processing, diffusion will gradually blur an image by damping out the high frequency modes. Image enhancement is the reverse process, and can be based on running the heat flow backwards in some stable manner. 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. One option would be to restrict the backwards evolution to the first few Fourier modes, which prevents the small scale fluctuations from overwhelming the computation. Ill-posed problems also arise in the reconstruction of subterranean profiles 1/19/12 108
c 2012 Peter J. Olver
from seismic data, a crucial problem in the oil and gas industry. These and other applications are driving contemporary research into how to cleverly circumvent the ill-posedness of backwards diffusion processes. Remark : The irreversibility of the heat equation, like our earlier discussion of the irreversibility of nonlinear transport in the presence of shock waves, highlights a crucial distinction between partial differential equations and ordinary differential equations. Ordinary differential equations are always reversible -- the existence, uniqueness and continuous dependence properties of solutions are all equally valid in reverse time (although their detailed qualitative and quantitative properties will, of course, depend upon whether time is running forwards or backwards). The irreversibility and ill-posedness of partial differential equations modeling the diffusive processes in our universe may explain why Time's Arrow points exclusively to the future. The Heated Ring Redux Let us next consider the periodic boundary value problem modeling heat flow in an insulated circular ring. We fix the length of the ring to be = 2 , with - x representing the "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 u 2u = , t x2 subject to periodic boundary conditions u(t, - ) = u(t, ), u u (t, - ) = (t, ), x x t 0, (4.29) - < x < , t > 0, (4.28)
that ensure continuity of the solution when the angular coordinate switches from - to . The initial temperature distribution is u(0, x) = f (x), - < x . (4.30)
The resulting temperature u(t, x) will be a periodic function in x of period 2 . Substituting the separable solution ansatz (3.15) into the heat equation and the boundary conditions results in the periodic eigenvalue problem d2 v + v = 0, dx2 v(- ) = v(), v (- ) = v (). (4.31)
As we already noted in Section 3.1, the eigenvalues of this particular boundary value problem are n = n2 where n = 0, 1, 2, . . . is a non-negative integer; the corresponding eigenfunctions 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 1/19/12 109
c 2012 Peter J. Olver
possessing two linearly independent eigenfunctions. The corresponding eigensolutions to the heated ring equation (4.2829) are un (t, x) = e-n t cos n x,
2
un (t, x) = e-n t sin n x,
2
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 t sin n x ,
2
2
(4.32)
with as yet unspecified coefficients an , bn . The initial conditions require
u(0, x) =
1 2
a0 +
n=1
an cos n x + bn sin n x = f (x),
(4.33)
which is precisely the complete Fourier series (3.34) of the initial temperature profile f (x). Consequently, an = 1
f (x) cos n x dx,
-
bn =
1
f (x) sin n x dx,
-
(4.34)
are its usual Fourier coefficients (3.35). As in the Dirichlet problem, after the initial instant, the high frequency terms in the 2 series (4.32) become extremely small, since e-n t 1 for n 0. Therefore, as soon as t > 0, the solution becomes instantaneously smooth, and quickly degenerates into what is in essence 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 mode, associated with the null eigenvalue 0 = 0. Consequently, the solution will converge, at an exponential rate, to a constant temperature profile: u(t, x) -
1 2
a0 =
1 2
f (x) dx,
-
which equals the average of the initial temperature profile. In physical terms, since the insulation prevents any heat energy from escaping the ring, it rapidly redistributes itself so that the ring achieves a uniform constant temperature -- its eventual equilibrium state. Prior to attaining equilibrium, only the very lowest frequency Fourier modes will still be noticeable, and so the solution will asymptotically look like u(t, x) where a1 = r1 cos 1 = 1 2
1 2
a0 + e- t (a1 cos x + b1 sin x) =
1 2
a0 + r1 e- t cos(x + 1 ),
(4.35)
f (x) cos x dx,
-
b1 = r1 sin 1 =
1 2
f (x) sin x dx.
-
Thus, for most initial data, the solution approaches thermal equilibrium at an exponential rate of e- t . The exceptions are when a1 = b1 = 0, for which the rate of convergence is 1/19/12 110
c 2012 Peter J. Olver
even faster, namely at a rate e- k t where k is the smallest integer such that at least one of the k th order Fourier coefficients ak , bk is nonzero. In fact, once we are convinced that the bar must tend to thermal equilibrium as t , we can predict the final temperature without knowing the explicit solution formula. Our derivation in Section 4.1 implies that the heat equation has the form of a conservation law (4.1), with the conserved density being the temperature u(t, x). As in (4.2), the integrated form of the conservation law reads d dt
-
2
u(t, x) dx =
u 2u (t, x) dx = (t, x) dx 2 - t - x u u (t, ) - (t, - ) = 0, = x x
where the flux terms cancel thanks to the periodic boundary conditions (4.29). Physically, any flux out of one end of the circular bar is immediately fed into the other, abutting end, and so there is no net loss of heat energy. We conclude that, for the periodic boundary value problem , the total heat
H(t) =
-
u(t, x) dx = constant
(4.36)
remains constant for all time. In general, a system is in equilibrium if it remains unaltered as time progresses. Thus, any equilibrium configuration has the form u = u (x), and hence satisfies u /t = 0. If, in addition, u (x) is an equilibrium solution to the periodic heat equation (4.2831), then it must satisfy u 2 u u u (4.37) =0= , u (- ) = u (), (- ) = (). t x2 x x In other words, u is a solution to the periodic boundary value problem (4.31) for the null eigenvalue = 0. The null eigenfunctions (including the zero solution) are the possible equilibrium solutions. In particular, for the periodic boundary value problem, the null eigenfunctions are constant, and therefore solutions to the periodic heat equation will tend to a constant equilibrium temperature. Now, once we know the solution tends to a constant, u(t, x) a as t , then its total heat tends to
H(t) =
-
u(t, x) dx -
a dx = 2 a
-
as
t - .
On the other hand, as we just demonstrated, the total heat is constant, so
H(t) = H(0) =
-
u(0, x) dx =
-
f (x) dx.
In contrast, this does not hold for the Dirichlet boundary value problem. In that case, the heat energy steadily goes to 0 due to the out-flux of heat through the ends of the bar. See Exercise for further details.
1/19/12
111
c 2012
Peter J. Olver
Combining these two, we conclude that
f (x) dx = 2 a,
-
and so the equilibrium temperature a =
1 2
f (x) dx
-
equals the average of the initial temperature. This reconfirms our earlier result, but avoids having to rely on the explicit series solution formula. 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 2u (4.38) = , u(t, 0) = , u(t, ) = , t 0, t x2 into a homogeneous Dirichlet problem. We begin by solving for the equilibrium temperature profile. As in (4.37), the equilibrium does not depend on t and hence satisfies the boundary value problem 2 u u =0= , u (0) = , u () = . t x2 Solving the ordinary differential equation, u (x) = a + b x, where the constants a, b are fixed by the boundary conditions; we conclude that the equilibrium solution is a straight line connecting the boundary values: u (x) = + The difference - x. (4.39)
- x (4.40) 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: - x f (x). (4.41) Solving the resulting homogeneous initial-boundary value problem, we write u(t, x) in Fourier series form (4.22), where the Fourier coefficients are specified by the modified initial data f (x) in (4.41). The solution to the inhomogeneous boundary value problem thus has the series form u(0, x) = u(t, x) - u (x) = f (x) - - - x + u(t, x) = + 1/19/12
bn exp
n=1
-
n2 2 t 2
sin
n x ,
(4.42)
Peter J. Olver
112
c 2012
where 2 bn =
0
f (x) sin
n x dx,
n = 1, 2, 3, . . . .
(4.43)
Since u(t, 0) decays to zero at an exponential rate as t , the actual temperature profile (4.42) will asymptotically decay to the equilibrium profile, - x u(t, x) - u (x) = + at the same exponentially fast rate, governed by the first eigenvalue 1 = 2 /2 -- unless b1 = 0, in which case the decay rate is even faster. This method does not work as well 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, (4.44) does satisfy the homogeneous boundary conditions, but now solves an inhomogeneous version of the heat equation: 2u u (t) - (t) u = + h(t, x), where h(t, x) = - (t, x) = - (t) - x. (4.45) t x2 t Solution techniques for the latter partial differential equation will be discussed below. 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 inside 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 = 2.0 10-7 sec-1 365.25 days (4.47) (4.46)
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 - ,
(4.48)
In this case, u (t, x) is not an equilibrium solution. Indeed, we do not expect the bar to go to equilibrium if the temperature of its endpoints it constantly changing.
1/19/12
113
c 2012
Peter J. Olver
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 (4.46, 48) 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.
(4.49)
(4.50)
lim v(x) = 0.
The solutions to the ordinary differential equation are i / x = e /2 (1+ i ) x , v2 (x) = e- i / x = e- /2 (1+ i ) x . v1 (x) = e The first solution is exponentially growing as x , and so not germane to our problem. The solution to the boundary value problem must therefore be a multiple, v(x) = a e- /2 (1+ i ) x of the exponentially decaying solution. Substituting back into (4.50), we find the (complex) solution to the root cellar problem to be u(t, x) = a e- x /2 e i t- /2 x . (4.51) The corresponding real solution is obtained by taking the real part, u(t, x) = a e- x /2 cos t - x . 2
(4.52)
The first factor in (4.52) 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 factor is periodic with the same annual frequency . The interesting feature is 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 1/19/12 114
c 2012 Peter J. Olver
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 (4.47), 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 Minneapolis, 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.
4.2. The Wave Equation.
Let us return to the one-dimensional wave equation 2u 2u = c2 , t2 x2 (4.53)
with constant wave speed c, used to model the vibrations of bars and strings. In Chapter 2, we learned how to solve the wave equation by the method of d'Alembert. Unfortunately, d'Alembert's approach does not extend to other equations of interest to us, and so alternative solution techniques, particularly those based on Fourier methods, are worth developing. Indeed, the resulting series solutions provide valuable insight into wave dynamics on bounded intervals. Separation of Variables and Fourier Series Solutions One of the oldest -- and still one of the most widely used -- techniques for constructing explicit analytical solutions to a wide range of linear partial differential equations is the method of separation of variables. We have, in fact, already employed a simplified version of the method when constructing the eigensolutions to the heat equation as an exponential function of t times a function of x. In general, the separation of variables method seeks solutions to the partial differential equation which can be written as the product of functions of the individual independent variables. For the wave equation, we seek solutions u(t, x) = w(t) v(x) (4.54) that can be written as the product of a function of t alone times a function of x alone. When the method succeeds (which is not guaranteed in advance), both factors are found as solutions to certain ordinary differential equations. 1/19/12 115
c 2012 Peter J. Olver
Let us see whether such an expression can possibly solve the wave equation. First of all, 2u = w (t) v(x), t2 2u = w(t) v (x), x2
where the primes indicate ordinary derivatives. Substituting these expressions into the wave equation (4.53), we find w (t) v(x) = c2 w(t) v (x). Dividing both sides by w(t) v(x) (which we assume is not identically zero, as otherwise the solution would be trivial) yields w (t) v (x) = c2 , w(t) v(x) which effectively "separates" the t and x variables on each side of the equation, hence the name "separation of variables". Now, how could a function of t alone be equal to a function of x alone? A moment's reflection should convince the reader that this can happen if and only if the two functions are constant , so
w (t) 2 v (x) =c = , w(t) v(x)
(4.55)
where we use to indicate the common separation constant. Thus, the individual factors w(t) and v(x) must satisfy ordinary differential equations d2 w - w = 0, dt2 d2 v - 2 v = 0, 2 dx c
as promised. We already know how to solve both of these ordinary differential equations by elementary techniques. There are three different cases, depending on the sign of the separation constant . As a result, each value of leads to 4 independent separable solutions to the wave equation, as listed in the following table.
Technical detail: one should assume that the underlying domain is connected for this to be valid as stated. However, in practice, this technicality can safely be ignored.
1/19/12
116
c 2012
Peter J. Olver
Separable Solutions to the Wave Equation
w(t)
v(x) x x cos , sin c c 1, x e- x/c , e x/c
u(t, x) = w(t) v(x) cos t cos x x , cos t sin , c c x x sin t cos , sin t sin c c 1, x, t, t x e- (t+x/c) , e (t-x/c) , e- (t-x/c) , e (t+x/c)
= - 2 < 0 =0 = 2 > 0
cos t, sin t 1, t e- t , e t
So far, we have not taken the boundary conditions into account. 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, ). Substituting the separable ansatz (4.55), we find that v(x) must satisfy d2 v - 2 v = 0, 2 dx c v(0) = 0 = v(). (4.56)
The complete system of solutions to this boundary value problem were found in (4.20): n x n c 2 , n = - , n = 1, 2, 3, . . . . Hence, referring to the table, the corresponding separable solutions are vn (x) = sin n x n c t n x n c t (4.57) sin , un (t, x) = sin sin . We will now employ these solutions to construct a candidate series solution to the wave equation subject to the prescribed boundary conditions: un (t, x) = cos
u(t, x) =
n=1
bn cos
n x n c t n x n c t sin + dn sin sin
.
(4.58)
The solution is thus a linear combination of the natural Fourier modes vibrating with frequencies n c n n = = , n = 1, 2, 3, . . . . (4.59) Observe that, 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. 1/19/12 117
c 2012 Peter J. Olver
The Fourier coefficients bn and dn in (4.58) will be uniquely determined by the initial conditions u (0, x) = g(x), 0 < x < . u(0, x) = f (x), t 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),
n c u n x (0, x) = dn sin = g(x). t n=1 n x dx, n = 1, 2, 3, . . . . (4.60)
Therefore, 2 bn = f (x) sin
0
are the Fourier sine coefficients (3.83) of the initial displacement f (x), while dn = 2 n c
g(x) sin
0
n x dx,
n = 1, 2, 3, . . . .
(4.61)
are rescaled versions of the Fourier sine coefficients of the initial velocity g(x). Example 4.3. A string of unit length is held taut in the center and then released. Our task 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. (4.62)
To be specific, we assume that the center of the string has been displaced by half a unit, and so the initial displacement is f (x) = x, 1 - x,
1 0 x 2, 1 x 1. 2
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, cos(2 k + 1) t sin(2 k + 1) x , (2 k + 1)2
c 2012
while dn = 0. Therefore, the solution is the Fourier sine series 4 u(t, x) = 2 1/19/12 (- 1)k
k=0
(4.63)
Peter J. Olver
118
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 4.3.
Plucked String Solution of the Wave Equation.
whose graph at times t = 0, .2, .4, .6, .8, 1. is depicted in Figure 4.3. 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 can, in fact, be proved as a consequence of the d'Alembert formula; see Exercise . Observe that, unlike the heat equation, the wave equation does not smooth out discontinuities and corners in the initial data. While the series form (4.58) of the solution is perhaps less satisfying than a d'Alembertstyle formula, we can still use it to deduce important qualitative properties of the solutions. First of all, since each term is periodic in t with period 2 /c, the entire solution is time periodic with that period: u(t + 2 /c, x) = u(t, x). In fact, after half a period, the solution reduces to u ,x c
=
n=1
(-1)n bn sin
n x =-
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). (4.64)
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 inherently periodic, because the fundamental frequencies (4.59) are all integer multiples of the lowest one: n = n 1 . Remark : The preceding analysis has important musical consequences. To the human ear, sonic vibrations that are integral multiples of a single frequency, and thus periodic 1/19/12 119
c 2012 Peter J. Olver
in time, sound harmonic, whereas those with irrationally related frequencies, and hence experiencing non-periodic vibrations, 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 12 and 13, the frequency ratios of the latter are irrationally related, and hence their motion is not periodic. For some reason, our appreciation of music is psychologically attuned to the differences between rationally related/periodic and irrationally related/quasi-periodic vibrations. Consider next a string with both ends left free, and so subject to the Neumann boundary conditions u u (t, 0) = 0 = (t, ). (4.65) x x The solutions of (4.56) satisfying v (0) = 0 = v () are now n c n x with n = , n = 0, 1, 2, 3, . . . . The resulting solution takes the form of a Fourier cosine series vn (x) = cos
u(t, x) = a0 + c0 t +
n=1
an cos
n x n c t n x n c t cos + cn sin cos
.
(4.66)
The first two terms come from the null eigenfunction v0 (x) = 1 with 0 = 0. The string vibrates with the same fundamental frequencies (4.59) 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. In general, the presence of null eigenfunctions implies that the wave equation admits unstable modes. Substituting (4.66) into the initial conditions u (0, x) = g(x), 0 < x < , t we find the Fourier coefficients are prescribed, as before, by the initial displacement and velocity, u(0, x) = f (x), 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
Note that we have not included the usual (4.66).
1 2
factor in the constant terms in the Fourier series
1/19/12
120
c 2012
Peter J. Olver
are equal to the average initial displacement and average initial velocity of the string. 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 string 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. The d'Alembert Formula for Bounded Intervals In Theorem 2.14, we derived the explicit d'Alembert formula 1 f (x - c t) + f (x + c t) + u(t, x) = 2 2c
x+c t
g(z) dz,
x-c t
(4.67)
for solving the basic initial value problem for the wave equation on an infinite interval: 2u 2u u = c2 , u(0, x) = f (x), (0, x) = - g(x), < x < . 2 2 t x t In this section we explain how to adapt the formula in order to solve initial-boundary value problems on bounded intervals, thereby effectively summing the Fourier series solution. The easiest case to deal with is the periodic problem on 0 x , with boundary conditions u(t, 0) = u(t, ), ux (t, 0) = ux (t, ). (4.68) 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 (4.67) will also be periodic in x, so u(t, x + ) = u(t, x). In particular, it satisfies the boundary conditions (4.68) 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. (4.69)
The resulting solution can be written as a Fourier sine series (4.58), and hence is both odd and 2 periodic in x. Therefore, to write the solution in d'Alembert form (4.67), 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, cf. Exercise , it satisfies the homogeneous Dirichlet boundary conditions (4.69) 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 = y for some 0 < y < . Its odd, periodic extension consists of two sets of replicas: those of 1/19/12 121
c 2012 Peter J. Olver
Figure 4.4.
Odd Periodic Extension of a Concentrated Pulse.
the same form occurring at positions y 2 , y 4 , . . . , and their upside down mirror images at the intermediate positions - y, - y 2 , - y 4 , . . . ; Figure 4.4 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. When a left and right moving pulse meet, they emerge from the interaction unaltered. 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 observed on the physical string. The effect is as if one were watching the full solution as it passes by a window of length . Such an observer will interpret what they see a bit differently. To wit, the original pulse starting at position 0 < y < splits up into two half size replicas that move off in opposite directions. As each half-size pulse reaches an end of the string, it meets a mirror image pulse that has been propagating in the opposite direction from the non-physical regime. The pulse is reflected at the end of the interval, and changes into an upside down mirror image moving in the opposite direction. The original positive pulse has moved off the end of the string just as its mirror image has moved into the physical regime. (A common physical realization 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 - y to instantaneously form a full-sized, but upside-down mirror image of the original disturbance -- in accordance with (4.64). 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 4.5, the first picture displays the initial displacement. In the second, it has split 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 string. 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 1/19/12 122 u (t, ) = 0, x
c 2012
(4.70)
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 4.5.
Solution to Wave Equation with Fixed Ends.
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), f (x + 2 ) = f (x), g(- x) = g(x), g(x + 2 ) = g(x).
The resulting d'Alembert solution (4.67) is also even and 2 periodic in x, and hence satisfies the boundary conditions, cf. Exercise (b). 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 studied two very different ways to solve the one-dimensional wave equation. The first, based on the d'Alembert formula, emphasizes their particle-like aspects, where individual wave packets collide with each other, or reflect at the boundary, all the while maintaining their overall form, while the second, based on Fourier analysis, emphasizes the vibrational or wave-like character of the solutions. Some solutions look like vibrating waves, while others appear 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 solutions combine into moving wave packets. The coexistence of particle and wave features is reminiscent of the long-running historical debate over the nature of light. Newton and his disciples proposed a particle-based theory, anticipating the modern concept of photons. However, until the beginning of the twentieth century, most physicists advocated a wave-like 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 both particle and wave features, depending upon the experiment and the physical situation. But a theoretical basis for the perplexing wave-particle duality could have already been found in the competing solution formulae for the classical wave equation! 1/19/12 123
c 2012 Peter J. Olver
Figure 4.6.
Planar Domain.
4.3. The Planar Laplace and Poisson Equations.
The two-dimensional Laplace equation is the second order linear partial differential equation 2u 2u + 2 = 0, (4.71) x2 y named in honor of the influential eighteenth century French mathematician PierreSimon Laplace. It, along with its higher dimensional versions, is arguably the most important differential equation in all of mathematics. A real-valued solution u(x, y) to the Laplace equation is known as a harmonic function. The space of harmonic functions can thus be identified as the kernel of the second order linear partial differential operator = 2 2 + 2, x2 y (4.72)
known as the Laplace operator , or Laplacian for short. The inhomogeneous or forced version, namely 2u 2u - 2 = f (x, y) (4.73) - [ u ] = - x2 y is known as Poisson's equation, named after SimonDenis Poisson, who was taught by e Laplace. The mathematical and physical reasons for including the minus signs will gradually become clear. Besides their theoretical importance, the Laplace and Poisson equations arise as the basic equilibrium equations in a remarkable variety of physical systems. For example, we may interpret u(x, y) as the displacement of a membrane, e.g., a drum skin; the inhomogeneity f (x, y) in the Poisson equation represents an external forcing over the surface of the membrane. Another example is in the thermal equilibrium of flat plates; here u(x, y) represents the temperature and f (x, y) an external heat source. In fluid mechanics, u(x, y) represents the potential function whose gradient v = u is the velocity vector of 1/19/12 124
c 2012 Peter J. Olver
h(x, y)
Figure 4.7. Dirichlet Boundary Conditions.
a steady planar fluid flow. Similar considerations apply to two-dimensional electrostatic and gravitational potentials. The dynamical counterparts to the Laplace equation are the two-dimensional versions of the heat and wave equations, to be analyzed in Chapter 12. Since both the Laplace and Poisson equations describe equilibrium configurations, they almost always appear the context of boundary value problems. We seek a solution u(x, y) to the partial differential equation defined at points (x, y) belonging to a bounded, open domain R 2 . The solution is required to satisfy suitable conditions on the boundary of the domain, denoted , which will consist of one or more simple, closed curves, as illustrated in Figure 4.6. As in one-dimensional boundary value problems, there are three especially important types of boundary conditions. The first are the fixed or Dirichlet boundary conditions, which specify the value of the function u on the boundary: u(x, y) = h(x, y) for (x, y) . (4.74)
Under mild regularity conditions on the domain , the boundary values h and the forcing function f , the Dirichlet conditions (4.74) serve to uniquely specify the solution u(x, y) to the Laplace or the Poisson equation. Physically, in the case of a free or forced membrane, the Dirichlet boundary conditions correspond to gluing the edge of the membrane to a wire at height h(x, y) over each boundary point (x, y) , as illustrated in Figure 4.7. Similarly, in the modeling of thermal equilibrium, a Dirichlet boundary condition represents the imposition of a prescribed temperature distribution, represented by the function h, along the boundary of the plate. The second important class are the Neumann boundary conditions u = u n = k(x, y) n
on
,
(4.75)
Extending existence and uniqueness results to irregular domains, e.g., those with fractal boundaries, is an active area of contemporary research. Interestingly, lack of regularity at sharp cusps underlies the electromagnetic phenomenon known as St. Elmo's fire, [ 144 ].
1/19/12
125
c 2012
Peter J. Olver
in which the normal derivative of the solution u on the boundary is prescribed. In general, n denotes the unit outwards normal to the boundary , i.e., the vector of unit length, n = 1, that is orthogonal to the tangent to the boundary and pointing away from the domain. For example, in thermomechanics, a Neumann boundary condition specifies the heat flux out of the plate through its boundary. The "no-flux" or homogeneous Neumann boundary conditions, where k(x, y) 0, correspond to a fully insulated boundary. In the case of a membrane, homogeneous Neumann boundary conditions correspond to a free, unattached edge of the drum. In fluid mechanics, the Neumann conditions prescribe the fluid flux through the boundary; in particular, homogeneous Neumann boundary conditions correspond to a solid boundary that the fluid cannot penetrate. There are also Robin boundary conditions u + a(x, y) u = k(x, y) on , n where a(x, y) > 0, that model insulated plates in heat reservoirs, or membranes attached to springs. Finally, one can mix the previous kinds of boundary conditions, imposing Dirichlet conditions on part of the boundary, and Neumann conditions on the complementary part. A typical mixed boundary value problem has the form u (4.76) = k on N, n with the boundary = D N being the disjoint union of a "Dirichlet segment", denoted by D, and a "Neumann segment" N . For example, if u represents the equilibrium temperature in a plate, then the Dirichlet segment of the boundary is where the temperature is fixed, while the Neumann segment is insulated, or, more generally, has prescribed heat flux. Similarly, when modeling the displacement of a membrane, the Dirichlet segment is where the edge of the drum is attached to a support, while the homogeneous Neumann segment is left hanging free. - u = f in , u = h on D, Separation of Variables Our first approach to solving the Laplace equation u = 2u 2u + 2 =0 x2 y (4.77)
will be based on the method of separation of variables. As in (4.54), we seek solutions that can be written as a product u(x, y) = v(x) w(y) (4.78) of a function of x alone times a function of y alone. We compute 2u = v (x) w(y), 2 x and so u = 1/19/12 2u = v(x) w (y), 2 y
2u 2u + 2 = v (x) w(y) + v(x) w (y) = 0. x2 y 126
c 2012 Peter J. Olver
We then separate the variables by placing all of the terms involving x on one side of the equation and all the terms involving y on the other; this is accomplished by dividing by v(x) w(y) and then writing the resulting equation in the separated form v (x) w (y) =- = . v(x) w(y) (4.79)
As we argued in (4.55), the only way a function of x alone can be equal to a function of y alone is if both functions are equal to a common separation constant . Thus, the factors v(x) and w(y) must satisfy the elementary ordinary differential equations v - v = 0, w + w = 0.
As before, the solution formulas depend on the sign of the separation constant . We list the resulting collection of separable harmonic functions in the following table. Separable Solutions to Laplace's Equation = - 2 < 0 =0 = >0
2
v(x) cos x, sin x 1, x e
- x
w(y) e- y , e y , 1, y
x
u(x, y) = v(x) w(y) e y cos x, e y sin x, e- y cos x, e- y sin x 1, x, y, x y e x cos y, e x sin y, e- x cos y, e- x sin y
, e
cos y, sin y
Since Laplace's equation is a homogeneous linear system, any linear combination of solutions is also a solution. So, we can build more general solutions as finite linear combinations, or, provided we pay proper attention to convergence issues, infinite series in the separable solutions. Our goal is to solve boundary value problems, and so we must ensure that the resulting combination satisfies the boundary conditions. But this is not such an easy task, unless the underlying domain has a rather special geometry. In fact, the only bounded domains on which we can explicitly solve boundary value problems using the preceding separable solutions are rectangles. So, we will concentrate on boundary value problems for Laplace's equation u = 0 on a rectangle R = { 0 < x < a, 0 < y < b }. (4.80)
To make progress, we will only allow nonzero boundary values on one of the four sides of the rectangle. To illustrate, we will focus on the following Dirichlet boundary conditions: u(x, 0) = f (x), u(x, b) = 0, u(0, y) = 0, u(a, y) = 0. (4.81)
Once we know how to solve this type of problem, we can employ linear superposition to solve the general Dirichlet boundary value problem on a rectangle; see Exercise for details. 1/19/12 127
c 2012 Peter J. Olver
Other boundary conditions can be treated in a similar fashion -- with the proviso that the condition on each side of the rectangle is either entirely Dirichlet or entirely Neumann. To solve the boundary value problem (4.8081), the first step is to narrow down the separable solutions to only those that respect the three homogeneous boundary conditions. The separable function u(x, y) = v(x) w(y) will vanish on the top, right and left sides of the rectangle provided v(0) = v(a) = 0, and w(b) = 0.
1 where sinh z = 2 (ez - e- z ) is the usual hyperbolic sine function. However, the second and third cases cannot satisfy the second boundary condition v(a) = 0, and so we discard them. The first case leads to the condition
Referring to the preceding table, the first condition v(0) = 0 requires sin x, = - 2 < 0, v(x) = x, = 0, sinh x, = 2 > 0,
v(a) = sin a = 0,
and hence
a = , 2 , 3 , . . . ,
an integral multiple of . The corresponding separation constants and solutions (up to constant multiple) are n = - 2 = - n2 2 , a2 vn (x) = sin n x , a n = 1, 2, 3, . . . . (4.82)
Note: We have merely recomputed the known eigenvalues and eigenfunctions of the familiar boundary value problem v - v = 0, v(0) = v(a) = 0. Since = - 2 < 0, we have w(y) = c1 e y + c2 e- y for constants c1 , c2 . The third boundary condition w(b) = 0 then requires that, up to constant multiple, wn (y) = sinh (b - y) = sinh We conclude that the harmonic functions un (x, y) = sin n x n (b - y) sinh , a a n = 1, 2, 3, . . . , (4.84) n (b - y) . a (4.83)
provide a complete list of separable solutions that satisfy the three homogeneous boundary conditions. It remains to analyze the inhomogeneous boundary condition along the bottom edge of the rectangle. To this end, let us try a linear superposition of the relevant separable solutions in the form of an infinite series:
u(x, y) =
n=1
cn un (x, y) =
n=1
cn sin
n (b - y) n x sinh , a a
c 2012 Peter J. Olver
1/19/12
128
whose coefficients c1 , c2 , . . . are to be prescribed by the remaining boundary condition. At the bottom edge, y = 0, we find
u(x, 0) =
n=1
cn sinh
n x n b sin = f (x), a a
0 x a,
(4.85)
which takes the form of a Fourier sine series for the function f (x). Let 2 bn = a
a
f (x) sin
0
n x dx a
(4.86)
be its Fourier sine coefficients, whence cn = bn / sinh(n b/a). We thus anticipate that the solution to the boundary value problem can be expressed as the infinite series
u(x, y) =
n=1
bn sin
n x n (b - y) sinh a a . n b sinh a
(4.87)
Does this series actually converge to the solution to the boundary value problem? Fourier analysis says that, under very mild conditions on the boundary function f (x), the answer is "yes". Suppose that its Fourier coefficients are uniformly bounded, | bn | M
a
for all
n 1,
(4.88)
which, according to (4.26) is true whenever f (x) is piecewise continuous or, more generally, integrable:
0
| f (x) | dx < . In this case, as you are asked to prove in Exercise , the n (b - y) a - 0 n b sinh a
coefficients of the Fourier sine series (4.87) go to zero exponentially fast: bn sinh
as
n -
for all
0 < y b,
(4.89)
and so, at each point inside the rectangle, the series can be well approximated by partial summation. Theorem 3.31 tells us that, for each 0 < y b, the solution u(x, y) is an infinitely differentiable function of x. Moreover, by term-wise differentiation of the series with respect to y, and use of Proposition 3.28, we also establish that the solution is infinitely differentiable with respect to y; see Exercise . (In fact, as we shall see, solutions to the Laplace equation are always analytic functions inside their domain of definition -- even when their boundary values are rather rough.) Since the individual terms all satisfy the Laplace equation, we conclude that the series (4.87) is indeed a classical solution to the boundary value problem. Example 4.4. A membrane is stretched over a wire in the shape of a unit square 1/19/12 129
c 2012 Peter J. Olver
Figure 4.8.
Square Membrane on a Wire. precise boundary conditions are y = 0, y = 0, y = 1, 0 y 1,
with one side bent in half, as graphed in Figure 4.8. The 0 x 1, x, 2 1 1 - x, x 1, 2 u(x, y) = 0, 0 x 1, 0, x = 0, 0, x = 1, f (x) = 4 = 2 x, 1 - x,
1 0 x 2, 1 2
The Fourier sine series of the inhomogeneous boundary function is readily computed:
0 y 1.
x 1, 4 = 2
sin 3 x sin 5 x sin x - + - 9 25
(-1)j
j =0
sin(2 j + 1)x . (2 j + 1)2
Specializing (4.87) when a = b = 1, we conclude that the solution to the boundary value problem can be expressed as a Fourier series 4 u(x, y) = 2
(-1)j
j =0
sin(2 j + 1) x sinh(2 j + 1) (1 - y) . (2 j + 1)2 sinh(2 j + 1)
In Figure 4.8 we plot the sum of the first 10 terms in the series. This gives a reasonably good approximation to the actual solution, except when we are very close to the raised corner of the boundary wire -- which is the point of maximal displacement of the membrane. Polar Coordinates The method of separation of variables can be successfully exploited in certain other very special geometries. One particularly important case is a circular disk. To be specific, 1/19/12 130
c 2012 Peter J. Olver
let us take the disk to have radius 1 and centered at the origin. Consider the Dirichlet boundary value problem u = 0, x2 + y 2 < 1, and u = h, x2 + y 2 = 1, (4.90)
so that the function u(x, y) satisfies the Laplace equation on the unit disk and the specified Dirichlet boundary conditions on the unit circle. For example, u(x, y) might represent the displacement of a circular drum that is attached to a wire of height h(x, y) = h(cos , sin ) h(), - < , (4.91)
at each point (x, y) = (cos , sin ) on its edge. The rectangular separable solutions are not particularly helpful in this situation, and so we look for solutions that are better adapted to a circular geometry. This inspires us to adopt polar coordinates x = r cos , y = r sin , or r= x2 + y 2 , = tan-1 y , x (4.92)
and write the solution u(r, ) as a function thereof. Warning: We will usually retain the same symbol, e.g., u, when rewriting a function in a different coordinate system. This is the convention of tensor analysis, physics, and differential geometry, [2], that treats the function (scalar field) as an intrinsic object, which is concretely realized through its formula in any chosen coordinate system. For instance, if u(x, y) = x2 + 2 y in rectangular coordinates, then its expression in polar coordinates is u(r, ) = (r cos )2 + 2 r sin , not r 2 + 2 . This convention avoids the inconvenience of introducing new symbols when changing coordinates. We need to relate derivatives with respect to x and y to those with respect to r and . Performing a standard chain rule computation based on (4.92), we find = cos + sin , r x y = - r sin + r cos , x y sin = cos - x r r cos = sin + y r r , .
so
(4.93)
Applying the squares of the latter differential operators to u(r, ), we find, after a calculation in which many of the terms cancel, u = 2 u 1 u 1 2u 2u 2u + 2 = + + 2 = 0. x2 y r 2 r r r 2 (4.94)
The boundary conditions are imposed on the unit circle r = 1, and so, by (4.91), take the form u(1, ) = h(). (4.95) Keep in mind that, in order to be single-valued functions of x, y, the solution u(r, ) and its boundary values h() must both be 2 periodic functions of the angular coordinate: u(r, + 2 ) = u(r, ), 1/19/12 131 h( + 2 ) = h().
c 2012
(4.96)
Peter J. Olver
Polar separation of variables is based on the ansatz u(r, ) = v(r) w() (4.97)
that assumes that the solution is a product of functions of the individual variables. Substituting (4.97) into the polar form (4.94) of Laplace's equation yields v (r) w() + 1 1 v (r) w() + 2 v(r) w () = 0. r r
We now separate variables by moving all the terms involving r onto one side of the equation and all the terms involving onto the other. This is accomplished by first multiplying the equation by r 2 /v(r) w(), and then moving the final term to the right hand side: w () r 2 v (r) + r v (r) =- = . v(r) w() As in the rectangular case, a function of r can equal a function of if and only if both are equal to a common separation constant, which we call . The partial differential equation thus splits into a pair of ordinary differential equations r 2 v + r v - v = 0, w + w = 0, (4.98)
that will prescribe the separable solution (4.97). Observe that both have the form of an eigenfunction equation in which the separation constant plays the role of the eigenvalue. We are, as always, only interested in nonzero solutions. We have already solved the eigenvalue problem for w(). According to (4.96), w( + 2 ) = w() must be a 2 periodic function. Therefore, by our earlier discussion, this periodic boundary value problem has the nonzero eigenfunctions 1, sin n , cos n , for n = 1, 2, . . . . (4.99)
corresponding to the eigenvalues (separation constants) n = n2 , where n = 0, 1, 2, . . .. Fixing the value of , the equation for the radial component, r 2 v + r v - n2 v = 0, (4.100)
is no longer a constant coefficient ordinary differential equation. But, fortunately, it has the form of a second order Euler ordinary differential equation, [23, 108], and hence can be readily solved by substituting the power ansatz v(r) = r k . (See also Exercise .) Note that v (r) = k r k-1 , v (r) = k (k - 1) r k-2 , and hence, substituting into the differential equation, r 2 v + r v - n2 v = k (k - 1) + k - n2 r k = (k 2 - n2 )r k . Thus, r k is a solution if and only if k 2 - n2 = 0, 1/19/12 and hence 132 k = n.
c 2012 Peter J. Olver
For n = 0, we have found the two linearly independent solutions: v1 (r) = r n , v2 (r) = r -n , n = 1, 2, . . . . (4.101)
When n = 0, the power ansatz only yields the constant solution. But in this case, the equation r 2 v + r v = 0 is effectively first order and linear in v , and hence readily integrated. This provides the two independent solutions v1 (r) = 1, v2 (r) = log r, n = 0. (4.102)
Combining (4.99) and (4.101102), we produce a complete list of separable polar coordinate solutions to the Laplace equation: 1, log r, r n cos n , r -n cos n , r n sin n , r -n sin n , n = 1, 2, 3, . . . . (4.103)
Now, the solutions in the top row of (4.103) are continuous (in fact analytic) at the origin, where r = 0, whereas the solutions in the bottom row have singularities as r 0. The latter are not of use in the present situation since we require the solution remain bounded and smooth -- even at the center of the disk. Thus, we should only use the non-singular solutions to concoct a candidate series solution a u(r, ) = 0 + 2
an r n cos n + bn r n sin n .
n=1
(4.104)
The coefficients an , bn will be prescribed by the boundary conditions (4.95). Substituting r = 1, we find a0 + an cos n + bn sin n = h(). u(1, ) = 2 n=1 We recognize this as a standard Fourier series (3.29) (with replacing x) for the 2 periodic function h(). Therefore, an = 1
h() cos n d,
-
bn =
1
h() sin n d,
-
(4.105)
are precisely its Fourier coefficients, cf. (3.35). In this manner, we have produced a series solution (4.104) to the boundary value problem (4.9495). Remark : Introducing the complex variable z = x + i y = r e i = r cos + i r sin allows us to write z n = r n e i n = r n cos n + i r n sin n . Therefore, the non-singular separable solutions are the harmonic polynomials r n cos n = Re z n , 1/19/12 133 r n sin n = Im z n .
c 2012
(4.106) (4.107)
(4.108)
Peter J. Olver
The first few are listed in the following table: n 0 1 2 3 4 Re z n 1 x 2 x - y2 x3 - 3 x y 2 x4 - 4 x2 y 2 + y 4 Im z n 0 y 2xy 2 3 x y - y3 4 x3 y - 4 x y 3
The general formula is obtained by using the Binomial Formula to compute z n = (x + i y)n n n-3 n n-2 x ( i y)3 + + ( i y)n x ( i y)2 + 3 2 n n-3 3 n n-2 2 x y + , x y -i = xn + i n xn-1 y - 3 2 = xn + n xn-1 ( i y) + where n k = n! k ! (n - k) ! (4.109)
are the usual binomial coefficients. Separating the real and imaginary terms, we find the explicit formulae r n cos n = Re z n = xn - n n-4 4 n n-2 2 x y + , x y + 4 2 n n-5 5 n n-3 3 x y + , x y + r n sin n = Im z n = n xn-1 y - 5 3
(4.110)
for the two independent harmonic polynomials of degree n. Example 4.5. Consider the Dirichlet boundary value problem on the unit disk with u(1, ) = for - < < . (4.111)
The boundary data can be interpreted as a wire in the shape of a single turn of a spiral helix sitting over the unit circle. The wire has a single jump discontinuity, of magnitude 2 , at the boundary point (-1, 0). The required Fourier series h() = 2 sin - sin 2 sin 3 sin 4 + - + 2 3 4
was already computed in Example 3.3. Therefore, invoking our solution formula (4.104 105), r 2 sin 2 r 3 sin 3 r 4 sin 4 u(r, ) = 2 r sin - + - + (4.112) 2 3 4 1/19/12 134
c 2012 Peter J. Olver
Figure 4.9.
Membrane Attached to a Helical Wire.
is the desired solution, and is plotted in Figure 4.9. In fact, this series can be explicitly summed. In view of (4.108), u = 2 Im where = tan-1 z- z3 z4 z2 + - + = 2 Im log(1 + z) = 2 ph(1 + z) = 2 , 2 3 4 y 1+x (4.113)
(4.114)
is the angle that the line passing through the two points (x, y) and (-1, 0) makes with the x-axis, as sketched in Figure 4.10. You should try to convince yourself that, on the unit circle, 2 = has the correct boundary values. Observe that, even though the boundary values are discontinuous, the solution is an analytic function inside the disk. In fact, unlike the rectangular series (4.87), the general polar series solution formula (4.104) can, in fact, be summed in closed form! If we substitute the explicit Fourier formulae (4.105) into (4.104) -- remembering to change the integration variable to, say, to avoid a notational conflict -- we find a u(r, ) = 0 + an r n cos n + bn r n sin n 2 n=1 = 1/19/12 1 2
h() d
-
(4.115) 135
c 2012 Peter J. Olver
(x, y)
Figure 4.10.
Geometrical Construction of the Solution.
+
n=1
r n cos n 1 + 2 1 + 2
h() cos n d +
-
r n sin n
h() sin n d
-
= =
1 1
h()
-
r n cos n cos n + sin n sin n
n=1 n=1
d
h()
-
r n cos n ( - ) d.
We next show how to sum the final series. Using (4.107), we can write it as the real part of a geometric series: 1 + 2
r cos n = Re
n=1
n
1 + zn 2 n=1
= Re
= Re
(1 + z)(1 - z) 2 | 1 - z |2
Re (1 + z - z - | z |2 ) 1 - | z |2 1 - r2 = = = . 2 | 1 - z |2 2 | 1 - z |2 2 (1 + r 2 - 2 r cos )
z 1 + 2 1-z
= Re
1+z 2 (1 - z)
Substituting back into (4.115) establishes the important Poisson Integral Formula for the solution to the boundary value problem. Theorem 4.6. The solution to the Laplace equation in the unit disk subject to Dirichlet boundary conditions u(1, ) = h() is 1 u(r, ) = 2
h()
-
1 - r2 d. 1 + r 2 - 2 r cos( - )
(4.116)
Example 4.7. A uniform metal disk of unit radius has half of its circular boundary held at 1 , while the other half is held at 0 . Our task is to find the equilibrium temperature u(x, y). In other words, we seek the solution to the Dirichlet boundary value problem u = 0, 1/19/12 x + y < 1,
2 2
u= 136
1, 0,
x2 + y 2 = 1, y > 0, x2 + y 2 = 1, y < 0.
c 2012
(4.117)
Peter J. Olver
Figure 4.11.
Equilibrium Temperature of a Disk.
In polar coordinates, the boundary data is a (periodic) step function h() = 1, 0 < < , 0, - < < 0.
Therefore, according to the Poisson formula (4.116), the solution is given by 2 1 - 1 tan-1 1 - r , 0 < < , 2 r sin 1 - r2 1 1 d = , = 0, , u(r, ) = 2 2 0 1 + r 2 - 2 r cos( - ) 2 - 1 tan-1 1 - r , - < < 0. 2 r sin (4.118) 1 1 -1 where we use the usual branch - 2 < tan t < 2 of the inverse tangent. Reverting to rectangular coordinates, the equilibrium temperature has the explicit formula 2 2 1 - 1 tan-1 1 - x - y , x2 + y 2 < 1, y > 0, 2y 1 u(x, y) = (4.119) , x2 + y 2 < 1, y = 0, 2 2 2 - 1 tan-1 1 - x - y , x2 + y 2 < 1, y < 0. 2y The result is depicted in Figure 4.11. Averaging, the Maximum Principle, and Analyticity Let us investigate some important consequences of the Poisson integral formula (4.116).
The detailed derivation of the final formula is left to the reader as Exercise .
1/19/12
137
c 2012
Peter J. Olver
First, setting r = 0 yields the formula u(0, ) = 1 2
h() d.
-
(4.120)
The left hand side is the value of u at the origin -- the center of the disk -- and so independent of ; the right hand side is the average of its boundary values around the unit circle. This formula is a particular instance of an important general fact. Theorem 4.8. Let u(x, y) be harmonic inside a disk of radius a centered at a point (x0 , y0 ) with piecewise continuous (or, more generally, integrable) boundary values on the circle C = { (x - x0 )2 + (y - y0 )2 = a2 }. Then its value at the center of the disk is equal to the average of its values on the boundary circle: u(x0 , y0 ) = 1 2 a u ds =
C
1 2
-
u(x0 + a cos , y0 + a sin ) d.
(4.121)
Proof : We use the scaling and translation symmetries of the Laplace equation, cf. Exercises , to map the disk of radius a centered at (x0 , y0 ) to the unit disk centered at the origin. Specifically, we set U (x, y) = u(x0 + a x, y0 + a y). (4.122)
An easy chain rule computation proves that U (x, y) also satisfies the Laplace equation on the unit disk x2 + y 2 < 1, with boundary values h() = U (cos , sin ) = u(x0 + a cos , y0 + a sin ). Therefore, by (4.120) , U (0, 0) = 1 2
h() d =
-
1 2
U (cos , sin ) d.
-
Replacing U by its formula (4.122) produces the desired result.
Q.E.D.
An important consequence of the integral formula (4.121) is the Maximum Principle for harmonic functions. Proposition 4.9. If u is a nonconstant harmonic function defined on a domain , then u does not have a local maximum or local minimum at any interior point of . Proof : The average of a continuous real function lies strictly between its maximum and minimum values -- except in the trivial case when the function is constant. Since u is harmonic, it is continuous inside . So Theorem 4.8 implies that the value of u at a point (x, y) lies strictly between its maximal and minimal values on any sufficiently small circle centered at (x, y). This clearly excludes the possibility of u having a local maximum or minimum at (x, y). Q.E.D. Consequently, on any bounded domain, a harmonic function u achieves its maximum and minimum values only at boundary points. 1/19/12 138
c 2012 Peter J. Olver
Theorem 4.10. Let u(x, y) be a harmonic function defined on a bounded domain that is continuous on . Let m = min { u(x, y) | (x, y) } , m u(x, y) M, M = max { u(x, y) | (x, y) }, for all (x, y) .
be, respectively, its maximum and minimum values on the boundary. Then
Physically, if we interpret u(x, y) as the vertical displacement of a membrane stretched over a wire, then Theorem 4.10 says that, in the absence of external forcing, the membrane cannot have any internal bumps -- its highest and lowest points are necessarily on the boundary of the domain. This reconfirms our physical intuition: the restoring force exerted by the stretched membrane will serve to flatten any bump, and hence a membrane with a local maximum or minimum cannot be in equilibrium. A similar interpretation holds for heat conduction. A body in thermal equilibrium can achieve its maximum and minimum temperature only on the boundary of the domain. Indeed, heat energy would flow away from any internal maximum, or towards any local minimum, and so if the body contained a local maximum or minimum on its interior, it could not be in thermal equilibrium. The Maximum Principle immediately implies the uniqueness of solutions to the Dirichlet boundary value problem for both the Laplace and Poisson equations: Theorem 4.11. If u and u both satisfy the same Poisson equation - u = - u = f within a bounded domain , and u = u on , then u u throughout . Proof : By linearity, the difference v = u - u satisfies the homogeneous boundary value problem v = 0 in and v = 0 on . Our assumption implies that the maximum and minimum boundary values of v are both 0 = m = M . Theorem 4.10 implies that 0 v(x, y) 0 at all (x, y) . Thus v 0, and hence u u everywhere in . Q.E.D. Remark : The existence of solutions on general domains is a more challenging issue, and we refer the interested reader to more advanced texts, e.g., [35, 45, 52, 75], for precise statements of theorems and proofs. nth Finally, let us discuss the analyticity of harmonic functions. In view of (4.108), the order term in the polar series solution (4.104), namely, an r n cos n + bn r n sin n = an Re z n + bn Im z n = Re (an - i bn ) z n , is, in fact, a homogeneous polynomial in (x, y) of degree n. This means that, when written in rectangular coordinates x and y, (4.104) is, in fact, a power series for the harmonic function u(x, y). It is well known -- see [7] -- that any convergent power series converges to an analytic function -- in this case u(x, y). Moreover, the power series must, in fact, be the Taylor series for u(x, y) based at the origin, and so its coefficients are multiples of the derivatives of u at x = y = 0. Details are worked out in Exercise . We can adapt this argument to prove analyticity of all solutions to the Laplace equation. Note especially the contrast with the wave equation, which has many non-analytic solutions. 1/19/12 139
c 2012 Peter J. Olver
Theorem 4.12. A harmonic function is analytic at every point in the interior of its domain of definition. Proof : Let u(x, y) be a solution to the Laplace equation on the open domain R 2 . Let x0 = (x0 , y0 ) , and choose a > 0 such that the closed disk of radius a centered at x0 is entirely contained within : Da (x0 ) = { x - x0 a } , where is the usual Euclidean norm. Then the function U (x, y) defined by (4.122) is harmonic on the unit disk, with well-defined boundary values. Thus, by the preceding remarks, U (x, y) is analytic at every point inside the unit disk, and hence so is u(x, y) = U x - x0 y - y 0 , a a
at every point (x, y) in the interior of the disk Da (x0 ). Since x0 was arbitrary, this establishes the analyticity of u throughout the domain. Q.E.D. This concludes our discussion of the method of separation of variables and its consequences for the planar Laplace equation. The method can be used in a few other special coordinate systems. See [92, 95, 97] for a complete account, including the fascinating connections with the underlying symmetry properties.
4.4. Classification of Linear Partial Differential Equations.
We have, at last, been introduced to the three cardinal linear, second order partial differential equations for functions of two variables. The homogeneous versions of the trinity are (a) The wave equation: utt - c2 uxx = 0, hyperbolic, (b) The heat equation: ut - uxx = 0, parabolic, (c) Laplace's equation: uxx + uyy = 0, elliptic. The last column indicates the equation's type, in accordance with the standard taxonomy of partial differential equations; the explanation will appear momentarily. The wave, heat and Laplace equations are the prototypical representatives of these three fundamental genres of partial differential equations. Each genre has distinctive analytical features, physical manifestations, and even numerical solution schemes. Equations governing vibrations, such as the wave equation, are typically hyperbolic. Equations modeling diffusion, such as the heat equation, are parabolic. Hyperbolic and parabolic equations both typically represent dynamical processes, and so one of the independent variables is identified as time. On the other hand, equations modeling equilibrium phenomena, including the Laplace and Poisson equations, are usually elliptic, and only involve spatial variables. Elliptic partial differential equations are associated with boundary value problems, whereas parabolic and hyperbolic equations involve initial and initial-boundary value problems. While this tripartite classification into hyperbolic, parabolic, and elliptic equations initially appears in the bivariate context, the terminology, underlying properties, and associated physical models carry over to second order partial differential equations in higher 1/19/12 140
c 2012 Peter J. Olver
dimensions. Most of the partial differential equations arising in applications fall into one of these three categories, and it is fair to say that the field of partial differential equations splits into three distinct subfields. Or, rather four subfields, the last containing all the equations, including higher order equations, that do not fit into the preceding categorization. (One important example appears in Section 9.5.) The full classification of real, linear, second order partial differential equations for a scalar-valued function u(x, y) depending on two variables proceeds as follows. The most general such equation has the form L[ u ] = A uxx + B uxy + C uyy + D ux + E uy + F u = G, (4.123)
where the coefficients A, B, C, D, E, F are all allowed to be functions of (x, y), as is the inhomogeneity or forcing function G(x, y). The equation is homogeneous if and only if G 0. We assume that at least one of the leading coefficients A, B, C is not identically zero, as otherwise the equation degenerates to a first order equation. The key quantity that determines the type of such a partial differential equation is its discriminant = B 2 - 4 A C. (4.124) This should (and for good reason) remind the reader of the discriminant of the quadratic equation Q(, ) = A 2 + B + C 2 + D + E + F = 0. (4.125) Its solutions trace out a plane curve -- a conic section. In the nondegenerate cases, the discriminant (4.124) fixes its geometrical type: (a) a hyperbola when > 0, (b) a parabola when = 0, or (c) an ellipse when < 0. This classification provides the underlying rationale for the choice of terminology used to classify second order partial differential equations. Definition 4.13. At a point (x, y), the linear, second order partial differential equation (4.123) is called (a) hyperbolic (x, y) > 0, (b) parabolic (c) elliptic (d) singular if and only if (x, y) = 0, but A2 + B 2 + C 2 = 0, (x, y) < 0, A = B = C = 0.
In particular: The wave equation uxx - uyy = 0 has discriminant = 4, and is hyperbolic. The heat equation uxx - uy = 0 has discriminant = 0, and is parabolic. The Poisson equation uxx + uyy = - f has discriminant = -4, and is elliptic.
For dynamical equations, we identify y with the time variable t.
1/19/12
141
c 2012
Peter J. Olver
Example 4.14. When the coefficients A, B, C vary, the type of the partial differential equation may not remain fixed over the entire domain. Equations that change type are less common, as well as being much harder to analyze and solve, both analytically and numerically. One example arising in the theory of supersonic aerodynamics is the Tricomi equation 2u 2u - 2 = 0. (4.126) y x2 y Comparing with (4.123), we find that A = y, B = 0, C = -1, The discriminant in this particular case is = B 2 - 4 A C = 4 y, and hence the equation is hyperbolic when y > 0, elliptic when y < 0, and parabolic on the transition line y = 0. In the physical model, the hyperbolic region corresponds to subsonic flow, while the supersonic regions are of elliptic type. The transitional parabolic boundary represents the shock line between the sub- and super-sonic regions -- the familiar sonic boom as an airplane crosses the sound barrier. Remark : The classification into hyperbolic, parabolic, elliptic, and singular types carries over as stated to quasi-linear second order equations, whose coefficients A, . . . , G are allowed to depend on u and its first order derivatives, ux , uy . Now the type of the equation can vary with both the point in the domain and the particular solution being considered. Even more generally, for a fully nonlinear second order partial differential equation H(x, y, u, ux, uy , uxx , uxy , uyy ) = 0, one defines its discriminant to be = H uxy
2
while
D = E = F = G = 0.
(4.127)
-4
H H . uxx uyy
(4.128)
Its sign determines the type of the equation as above -- again depending on the point in the domain and the solution under consideration. Characteristics In Chapter 2, we learned how the characteristics guide the behavior of solutions to partial differential equations that govern wave phenomena. Characteristics play a similarly fundamental role in the study of more general hyperbolic partial differential equations. Indeed, they can be used to distinguish among the three classes of second order partial differential equations. Definition 4.15. The graph of the function y = y(x) is called a characteristic curve for the second order linear partial differential equation (4.123) if A(x, y) 1/19/12 dy dx
2
- B(x, y) 142
dy + C(x, y) = 0. dx
c 2012
(4.129)
Peter J. Olver
Alternatively, if the curve is given by the graph of x = x(y), then the characteristic equation (4.129) becomes dx dx 2 A(x, y) - B(x, y) = 0. (4.130) + C(x, y) dy dy For example, consider the hyperbolic wave equation utt - c2 uxx = 0. According to (4.129), a characteristic curve x(t) satisfies dx dt
2
- c2 = 0,
which implies that
dx = c. dt
We conclude that, in accordance with out previous usage, the characteristic curves are the straight lines of slope c, and there are two characteristic curves passing through each point of the (t, x)plane. On the other hand, the elliptic Laplace equation uxx + uyy = 0 has no (real) characteristic curves, since the characteristic equation (4.129) reduces to dy dx Finally, for the parabolic heat equation uxx - ut = 0, the characteristic curve equation is simply dt dx
2 2
+ 1 = 0.
= 0,
(since the first derivative term plays no role), and so there is only one characteristic curve passing through each point, namely the vertical line t = a. We note that the characteristic curve equation (4.129) is a quadratic equation for dy/dx. The number of real solutions to the equation depends on its discriminant = B 2 - 4 A C: In the hyperbolic case, > 0, and there are two real characteristic curves passing through each point; in the parabolic case, = 0, and there is just one real characteristic curve passing through each point; in the elliptic case, < 0, and there are no real characteristic curves. In this manner, elliptic, parabolic, and hyperbolic partial differential equations are distinguished by the number of (real) characteristic curves passing through a point -- namely, zero, one, and two, respectively. With further analysis, it can be shown that, as with the wave equation, signals and disturbances propagate along characteristic curves. Thus, hyperbolic equations share many qualitative properties in common with the wave equation, with signals moving in two different directions. For example, light rays move along characteristic curves, and are thereby subject to the optical phenomena of refraction and focusing. Similarly, since the 1/19/12 143
c 2012 Peter J. Olver
characteristic curves for the parabolic heat equation are the vertical lines t = a, this indicates that the effect of a disturbance at a point (t, x) = (a, b) is simultaneously felt along the entire contemporaneous vertical line t = a. Indeed, Section 9.1 analytically justifies the infinite speed of propagation of disturbances in the heat equation. As a result, the effect of an initial concentrated heat source is immediately felt all along the bar; see Section 9.1 for an amplification of this counterintuitive fact. Elliptic equations have no characteristics, and as a consequence, do not admit propagating signals; the effect of a localized disturbance, is felt throughout the domain. For example, an external force that is concentrated near a single point induces a displacement throughout the entire membrane. Remark : First order partial differential equations are not covered by the preceding classification, but are generally viewed as hyperbolic owing to the behavior of their solutions along their characteristic curves.
1/19/12
144
c 2012
Peter J. Olver