Home / Math and Physics Files / Math / ODEs / Peter Olver Notes

hw

PDF · 67 pages · 587.5 KB

Chapter 14 of Peter J. Olver's PDE/applied mathematics text (dated 2012), kept in the ODEs folder of Phil's files. It derives the heat equation from conservation of energy and Fourier's law, covers boundary conditions, Fourier series solutions, the wave equation with d'Alembert's formula, and finite difference numerical methods. The file name 'hw.pdf' suggests it may be used for homework.

AI-written summary; may contain errors.

Extracted text (machine-read; may contain errors)

Chapter 14 VibrationandDiﬀusion inOne–DimensionalMedia In this chapter, we study the solutions, both analyticaland numerical, to the two most important equations of one-dimensional continuum dynamic s. Theheat equation models the diﬀusion of thermal energy in a body; here, we analyze the case of a one-dimensional bar. The wave equation describes vibrations and waves in continuous media, includ ing sound waves, water waves, elastic waves, electromagnetic w aves, and so on. Again, we restrict our attention here to the case of waves in a one-dime nsional medium, e.g., a string, or a bar, or a column of air. The two- and three-dimensional ve rsions of these fundamental equations will be analyzed in the later Chapters 17 and 18. As we saw in Section 12.1, the basic solution strategy is insp ired by our eigenvalue- based methods for solving linear systems of ordinary diﬀere ntial equations. Substituting the appropriate exponential or trigonometric ansatz will e ﬀectively reduce the partial dif- ferential equation to a one-dimensional boundary value pro blem. The linear superposition principle implies that general solution can then be express ed as a inﬁnite series in the resulting eigenfunction solutions. In both cases consider ed here, the eigenfunctions of the one-dimensional boundary value problem are trigonometric , and so the solution to the partial diﬀerential equation takes the form of a time-depen dent Fourier series. Although we cannot, in general, analytically sum the Fourier series t o produce a simpler formula for the solution, there are a number of useful observations that can be gleaned from it. Inthecaseoftheheatequation, thesolutionsdecay exponen tiallyfast tothermalequi- librium, at a rate governed by the smallest positive eigenva lue of the associated boundary value problem. The higher order Fourier modes damp out very r apidly, and so the heat equation can be used to automatically smooth and denoise sig nals and images. It also implies that the heat equation cannot be run backwards in tim e — determining the initial temperature proﬁle from a later measurement is an ill-posed problem. The response to a concentrated unit impulse leads to the fundamental soluti on, which can then be used to construct integral representations of the solution to th e inhomogeneous heat equation. We will also explain how to exploit the symmetry properties o f the diﬀerential equation in order to construct new solutions from known solutions. In the case of the wave equation, each Fourier mode vibrates w ith its natural fre- quency. In a stable situation, the full solution is a linear c ombination of these fundamen- tal vibrational modes, while an instability induces an extr a linearly growing mode. For one-dimensional media, the natural frequencies are integr al multiples of a single lowest fre- quency, and hence the solution is periodic, which, in partic ular, explains the tonal qualities 12/11/12 739 c/ci∇cleco√y∇t2012 Peter J. Olver of string and wind instruments. The one-dimensional wave eq uation admits an alternative explicit solution formula, due to d’Alembert, which points out the role of characteristics in signal propagation and the behavior of solutions. The analy tic and series solution methods serve to shed complementary lights on the physical phenomen a of waves and vibration in continuous media. Following our analytical study, we introduce several basic numerical solution methods for both the heat and the wave equations. We begin with a gener al discussion of ﬁnite diﬀerence formulae for numerically approximating derivat ives of functions. The basic ﬁnite diﬀerence scheme is obtained by replacing the derivatives in the equation by t he appropri- ate numerical diﬀerentiation formulae. However, there is n o guarantee that the resulting numerical scheme will accurately approximate the true solu tion, and further analysis is required to elicit bona ﬁde, convergent numerical algorith ms. In dynamical problems, the ﬁnite diﬀerence schemes replace the partial diﬀerential eq uation by an iterative linear ma- trix system, and the analysis of convergence relies on the me thods covered in Chapter 10. In preparation for our treatment of partial diﬀerential equ ations in higher dimen- sions, the ﬁnal section introduces a general framework for d ynamics that incorporates both discrete dynamics, modeled by linear systems of ordina ry diﬀerential equations, and continuum dynamics, modeled by the heat and wave equations, their multidimensional counterparts, as well as the equilibrium boundary value pro blems governing bars, beams, membranes and solid bodies. Common features, based on eigen values, are discussed. 14.1. The Diﬀusion and Heat Equations. Let us begin with a physical derivation of the heat equation f rom ﬁrst principles of thermodynamics. The reader solely interested in the mathem atical developments can skip ahead to the following subsection. However, physical insig ht can often play an critical role in understanding the underlying mathematics, and is neglec ted at one’s peril. We consider a bar — meaning a thin, heat-conducting body of le ngthℓ. “Thin” means that we can regard the bar as a one-dimensional continuum wit h no signiﬁcant transverse temperature variation. We use 0 ≤x≤ℓto denote the position along the bar. Our goal is to ﬁnd the temperature u(t,x) of the bar at position xand time t. The dynamical equations governing the temperature are based on three fund amental physical laws. The ﬁrst law is that, in the absence of external sources, ther mal energy can only enter the bar through its ends. In physical terms, we are assuming t hat the bar is fully insulated along its length. Let ε(t,x) to denote the thermal energy in the bar at position xand time t. Consider asmallsectionofthebarlyingbetween xandx+∆x. Thetotalamountofheat energy contained in this section is obtained by integrating (summing):/integraldisplayx+∆x xε(t,y)dy. Further, let w(t,x) denote the heat ﬂux , i.e., the rate of ﬂow of thermal energy along the bar. We use the convention that w(t,x)>0 means that the energy is moving to the right, whilew(t,x)<0 if it moves to the left. The ﬁrst law implies that the rate of c hange in the thermal energy in any section of the bar is equal to the tot al heat ﬂux, namely the amount of the heat passing through its ends. Therefore, in vi ew of our sign convention on 12/11/12 740 c/ci∇cleco√y∇t2012 Peter J. Olver the ﬂux, ∂ ∂t/integraldisplayx+∆x xε(t,y)dy=−w(t,x+∆x)+w(t,x), the latter two terms denoting the respective ﬂux of heat intothe section of the bar at its right and left ends. Assuming suﬃcient regularity of the int egrand, we are permitted to bring the derivative inside the integral. Thus, dividing bo th sides of the resulting equation by ∆x, 1 ∆x/integraldisplayx+∆x x∂ε ∂t(t,y)dy=−w(t,x+∆x)−w(t,x) ∆x. In the limit as the length ∆ x→0, the right hand side of this equation converges to minus thexderivative of w(t,x), while, by the Fundamental Theorem of Calculus, the left ha nd side converges to the integrand ∂ε/∂tat the point x; the net result is the fundamental diﬀerential equation ∂ε ∂t=−∂w ∂x(14.1) relating thermal energy εand heat ﬂux w. A partial diﬀerential equation of this partic- ular form is known as a conservation law , and, in this instance, formulates the law of conservation of thermal energy. See Exercise for details. The second physical law is a constitutive assumption , based on experimental evidence. In most physical materials, thermal energy is found to be pro portional to temperature, ε(t,x) =σ(x)u(t,x). (14.2) The factor σ(x) =ρ(x)χ(x)>0 is the product of the densityρof the material and its speciﬁc heat χ, which is the amount of heat energy required to raise the temperature of a unit mas s of the material by one unit. Note that we are assuming the bar is not changing in time, and s o physical quantities such as density and speciﬁc heat depend only on position x. We also assume, perhaps with less physical justiﬁcation, that the material properties do not depend upon the temperature; otherwise, we would be led to a much more diﬃcult nonlinear di ﬀusion equation. The third physical law relates the heat ﬂux to the temperatur e. Physical experiments in a wide variety of materials indicate that the heat energy m oves from hot to cold at a rate that is in direct proportion to the rate of change — meani ng the derivative — of the temperature. The resulting linear constitutive relation w(t,x) =−κ(x)∂u ∂x(14.3) is known as Fourier’s Law of Cooling . The proportionality factor κ(x)>0 is called the thermal conductivity of the bar at position x. A good heat conductor, e.g., silver, will have high conductivity, while a poor conductor, e.g., glass , will have low conductivity. The minus sign tells us that heat energy moves from hot to cold ; if∂u ∂x(t,x)>0 the 12/11/12 741 c/ci∇cleco√y∇t2012 Peter J. Olver temperature is increasing from left to right, and so the heat energy moves back to the left, with consequent ﬂux w(t,x)<0. Combining the three laws (14.1), (14.2) and (14.2) produces the basic partial diﬀer- ential equation ∂ ∂t/bracketleftbig σ(x)u/bracketrightbig =∂ ∂x/parenleftbigg κ(x)∂u ∂x/parenrightbigg , 0< x < ℓ, (14.4) governingthediﬀusionofheatinanon-uniform bar. Theresu ltinglinear diﬀusion equation is used to model a variety of diﬀusive processes, including h eat ﬂow, chemical diﬀusion, population dispersion, and the spread of infectious diseas es. If, in addition, we allow external heat sources h(t,x), then the linear diﬀusion equation acquires an inhomogene ous term: ∂ ∂t/bracketleftbig σ(x)u/bracketrightbig =∂ ∂x/parenleftbigg κ(x)∂u ∂x/parenrightbigg +h(t,x), 0< x < ℓ. (14.5) In order to uniquely prescribe the solution u(t,x) to the diﬀusion equation, we need to specify the initial temperature distribution u(t0,x) =f(x), 0≤x≤ℓ, (14.6) along the bar at an initial time t0. In addition, we must impose suitable boundary condi- tionsatthetwoends ofthebar. Aswiththeequilibriumequat ionsdiscussed inChapter 11, there are three common physical types. The ﬁrst is a Dirichlet boundary condition , where an end of the bar is held at prescribed temperature. Thus, the boundary condition u(t,0) =α(t) (14 .7) ﬁxesthetemperature attheleft handend ofthebar. Alternat ively,the Neumann boundary condition ∂u ∂x(t,0) =ξ(t) (14.8) prescribes the heat ﬂux w(t,0) =−κ(0)∂u ∂x(t,0) at the left hand end. In particular, the homogeneous Neumann condition with ξ(t)≡0 corresponds to an insulated end, where no heat can ﬂow in or out. Each end of the bar should have one or t he other of these boundary conditions. For example, a bar with both ends havin g prescribed temperatures is governed by the pair of Dirichlet boundary conditions u(t,0) =α(t), u (t,ℓ) =β(t), (14.9) whereas a bar with two insulated ends requires two homogeneo us Neumann boundary conditions∂u ∂x(t,0) = 0,∂u ∂x(t,ℓ) = 0. (14.10) Themixedcase, withoneendﬁxedandtheotherinsulated, iss imilarlyformulated. Finally, theperiodic boundary conditions u(t,0) =u(t,ℓ),∂u ∂x(t,0) =∂u ∂x(t,ℓ), (14.11) 12/11/12 742 c/ci∇cleco√y∇t2012 Peter J. Olver correspond to a circular ringof length ℓ. As before, we are assuming the heat is only allowed to ﬂow around the ring — insulation prevents any radi ation of heat from one side of the ring to the other. The Heat Equation In this book, we will retain the term “heat equation” to refer to the homogeneous case, in which the bar is made of a uniform material, and so its densityρ, conductivity κ, and speciﬁc heat χare all positive constants. Under these assumptions, the ho mogeneous diﬀusion equation (14.4) reduces to the heat equation ∂u ∂t=γ∂2u ∂x2(14.12) for the temperature u(t,x) in the bar at time tand position x. The constant γ=κ σ=κ ρχ(14.13) iscalledthe thermal diﬀusivity ofthe bar, andincorporates allofitsrelevant physical pro p- erties. The solution u(t,x) will be uniquely prescribed once we specify initial condit ions (14.6) and a suitable pair of boundary conditions at the ends of the bar. As we learned in Section 12.1, the elementary, separable sol utions to the heat equation are based on the exponential ansatz u(t,x) =e−λtv(x), (14.14) wherev(x) is a time-independent function. Substituting the solutio n formula (14.14) into (14.12) and canceling the common exponential factors, we ﬁn d thatv(x) must solve the ordinary diﬀerential equation −γv′′=λv. In other words, vis aneigenfunction witheigenvalue λ, for the second derivative operator K=−γD2. Once we determine the eigenvalues and eigenfunctions, we w ill be able to reconstruct the solution u(t,x) as a linear combination, or, rather, inﬁnite series in the corresponding separable eigenfunction solutions. Let us consider the simplest case of a uniform bar held at zero temperature at each end. For simplicity, we take the initial time to be t0= 0, and so the initial and boundary conditions are u(t,0) = 0, u(t,ℓ) = 0, t ≥0, u(0,x) =f(x), 0< x < ℓ.(14.15) According to the general prescription, we need to solve the e igenvalue problem γd2v dx2+λv= 0, v (0) = 0, v (ℓ) = 0. (14.16) As noted in Exercise , positive deﬁniteness of the underlying diﬀerential opera torK= −γD2when subject to Dirichlet boundary conditions implies that we need only look for positive eigenvalues: λ >0. In Exercises –, the skeptical reader is asked to check 12/11/12 743 c/ci∇cleco√y∇t2012 Peter J. Olver explicitly that if λ≤0 orλis complex, then the boundary value problem (14.16) admits only the trivial solution v(x)≡0. Settingλ=γω2withω >0, the general solution to the diﬀerential equation is a trigonometric function v(x) =acosωx+bsinωx, wherea,bare constants whose values will be ﬁxed by the boundary condi tions. The boundary condition at x= 0 requires a= 0. The second boundary condition requires v(ℓ) =bsinωℓ= 0. Therefore, assuming b/ne}ationslash= 0, as otherwise the solution is trivial, ωℓmust be an integer multiple of π, and so ω=π ℓ,2π ℓ,3π ℓ, ... . Weconclude thattheeigenvaluesandeigenfunctionsoftheb oundaryvalueproblem (14.16) are λn=γ/parenleftignπ ℓ/parenrightig2 , vn(x) = sinnπx ℓ, n = 1,2,3,.... (14.17) Thecorrespondingseparablesolutions(14.14)totheheate quationwiththegivenboundary conditions are un(t,x) = exp/parenleftbigg −γn2π2t ℓ2/parenrightbigg sinnπx ℓ, n = 1,2,3,... . (14.18) Each represents a trigonometrically oscillating temperat ure proﬁle that maintains its form while decaying to zero at an exponential rate. The ﬁrst of the se, u1(t,x) = exp/parenleftbigg −γπ2t ℓ2/parenrightbigg sinπx ℓ, experiences the slowest decay. The higher “frequency” mode sun(t,x),n≥2, all go to zero at a faster rate, with those having a highly oscillatory temperature proﬁle, where n≫0, eﬀectively disappearing almost instantaneously. Thus, small scale temperature ﬂuctuations tend to rapidly cancel out through diﬀusion of h eat energy. Linear superposition is used to assemble the general series solution u(t,x) =∞/summationdisplay n=1bnun(t,x) =∞/summationdisplay n=1bnexp/parenleftbigg −γn2π2t ℓ2/parenrightbigg sinnπx ℓ(14.19) as a combination of the separable solutions. Assuming that t he series converges, the initial temperature proﬁle is u(0,x) =∞/summationdisplay n=1bnsinnπx ℓ=f(x). (14.20) This has the form of a Fourier sine series (12.44) on the inter val [0,ℓ]. By orthogonality of the eigenfunctions — which is a direct consequence of the sel f-adjointness of the underlying 12/11/12 744 c/ci∇cleco√y∇t2012 Peter J. Olver 0.2 0.4 0.6 0.8 1.0 -0.3-0.2-0.10.10.20.3 0.2 0.4 0.6 0.8 1.0 -0.3-0.2-0.10.10.20.3 0.2 0.4 0.6 0.8 1.0 -0.3-0.2-0.10.10.20.3 0.2 0.4 0.6 0.8 1.0 -0.3-0.2-0.10.10.20.3 0.2 0.4 0.6 0.8 1.0 -0.3-0.2-0.10.10.20.3 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 Figure 14.1. A Solution to the Heat Equation. boundary value problem (14.16) — the coeﬃcients are determi ned by the inner product formulae (12.45), and so bn=2 ℓ/integraldisplayℓ 0f(x) sinnπx ℓdx, n = 1,2,3,... . (14.21) The resulting solution (14.19) describes the Fourier sine s eries for the temperature u(t,x) of the bar at each later time t≥0. It can be rigorously proved that, for quitegeneral initia l conditions, the Fourier series does indeed converge to a sol ution to the initial-boundary value problem, [ 187]. Example 14.1. Consider the initial temperature proﬁle u(0,x) =f(x) =  −x, 0≤x≤1 5, x−2 5,1 5≤x≤7 10, 1−x,7 10≤x≤1,(14.22) on a bar of length 1, plotted in the ﬁrst graph in Figure 14.1. U sing (14.21), the ﬁrst few Fourier coeﬃcients of f(x) are computed as b1=.0448..., b2=−.096..., b3=−.0145..., b4= 0, b5=−.0081..., b6=.0066..., b7=.0052..., b8= 0,... . Settingγ= 1, the resulting Fourier series solution to the heat equati on is u(t,x) =∞/summationdisplay n=1bnun(t,x) =∞/summationdisplay n=1bne−n2π2tsinnπx =.0448e−π2tsinπx−.096e−4π2tsin2πx−.0145e−9π2tsin3πx− ···. In Figure 14.1, the solution is plotted at the successive tim est=.,.02,.04,...,.1. Observe that the corners in the initial data are immediately smoothe d out. As time progresses, the solution decays, at an exponential rate of π2≈9.87, to a uniform, zero temperature, which is the equilibrium temperature distribution for the h omogeneous Dirichlet boundary 12/11/12 745 c/ci∇cleco√y∇t2012 Peter J. Olver conditions. Asthesolutiondecays to thermal equilibrium, it also assumes the progressively more symmetric shape of a single sine arc, of exponentially d ecreasing amplitude. Smoothing and Long Time Behavior The fact that we can write the solution to an initial-boundar y value problem in the form of an inﬁnite series is progress of a sort. However, beca use it cannot be summed in closed form, this “solution” is much less satisfying than a d irect, explicit formula. Never- theless, there are important qualitative and quantitative features of the solution that can be easily gleaned from such series expansions. If the initial data f(x) is piecewise continuous, then its Fourier coeﬃcients are u ni- formly bounded; indeed, for any n≥1, |bn| ≤2 ℓ/integraldisplayℓ 0/vextendsingle/vextendsingle/vextendsinglef(x) sinnπx ℓ/vextendsingle/vextendsingle/vextendsingledx≤2 ℓ/integraldisplayℓ 0|f(x)|dx≡M. (14.23) This property holds even for quite irregular data; for insta nce, the Fourier coeﬃcients (12.59) of the delta function are also uniformly bounded. Un der these conditions, each term in the series solution (14.19) is bounded by an exponent ially decaying function /vextendsingle/vextendsingle/vextendsingle/vextendsinglebnexp/parenleftbigg −γn2π2 ℓ2t/parenrightbigg sinnπx ℓ/vextendsingle/vextendsingle/vextendsingle/vextendsingle≤Mexp/parenleftbigg −γn2π2 ℓ2t/parenrightbigg . This means that, as soon as t >0, most of the high frequency terms, n≫0, will be extremely small. Only the ﬁrst few terms will be at all notice able, and so the solution essentially degenerates into a ﬁnite sum over the ﬁrst few Fo urier modes. As time in- creases, more and more of the Fourier modes will become negli gible, and the sum further degenerates into progressively fewer signiﬁcant terms. Ev entually, as t→ ∞,allof the Fourier modes will decay to zero. Therefore, the solution wi ll converge exponentially fast to a zero temperature proﬁle: u(t,x)→0 ast→ ∞, representing the bar in its ﬁnal uni- form thermal equilibrium. The fact that its equilibrium tem perature is zero is the result of holding both ends of the bar ﬁxed at zero temperature, and a ny initial heat energy will eventually be dissipated away through the ends. The last ter m to disappear is the one with the slowest decay, namely u(t,x)≈b1exp/parenleftbigg −γπ2 ℓ2t/parenrightbigg sinπx ℓ,where b1=1 π/integraldisplayπ 0f(x)sinxdx.(14.24) Generically, b1/ne}ationslash= 0, and the solution approaches thermal equilibrium expone ntially fast with rate equal to the smallest eigenvalue, λ1=γπ2/ℓ2, which is proportional to the thermal diﬀusivity divided by the square of the length of the bar. The longer the bar, or the smaller the diﬀusivity, the longer it takes for the eﬀe ct of holding the ends at zero temperaturetopropagatealongtheentirebar. Also, againp rovidedb1/ne}ationslash= 0, theasymptotic shape of the temperature proﬁle is a small sine arc, just as we observed in Example 14.1. In exceptional situations, namely when b1= 0, the solution decays even faster, at a rate 12/11/12 746 c/ci∇cleco√y∇t2012 Peter J. Olver 0.2 0.4 0.6 0.8 11234567 0.2 0.4 0.6 0.8 11234567 0.2 0.4 0.6 0.8 11234567 0.2 0.4 0.6 0.8 11234567 0.2 0.4 0.6 0.8 11234567 0.2 0.4 0.6 0.8 11234567 Figure 14.2. Denoising a Signal with the Heat Equation. equal to the eigenvalue λk=γk2π2/ℓ2corresponding to the ﬁrst nonzero term, bk/ne}ationslash= 0, in the series; its asymptotic shape now oscillates ktimes over the interval. The heat equation’s smoothing eﬀect on irregular initial da ta by fast damping of the high frequency modes underlies its eﬀectiveness for smooth ing 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 ﬂow fo r too long, all of the low frequency features will be also be smoothed out and the re sult will be a uniform, constant signal. Thus, the choice of stopping time t⋆is crucial to the success of this method. Figure 14.2 shows the eﬀect running the heat equatio n†, withγ= 1, to times t= 0.,.00001,.00005,.0001,.001,.01 on the same signal from Figure 13.5. Observe how quickly the noise is removed. By the ﬁnal time, the overall sm oothing eﬀect of the heat ﬂow has caused signiﬁcant degradation (blurring) of the ori ginal signal. The heat equation approach to denoising has the advantage that no Fourier coeﬃ cients need be explicitly computed, nor doesone need toreconstruct thesmoothedsign al from itsremaining Fourier coeﬃcients. The ﬁnal section discusses some numerical meth ods that can be used to solve the heat equation directly. Another, closely related observation is that, for any ﬁxed t imet >0 after the initial moment, the coeﬃcients in the Fourier series (14.19) decay e xponentially fast as n→ ∞. According to the discussion at the end of Section 12.3, this i mplies that the solution u(t,x) is a very smooth, inﬁnitely diﬀerentiable function of xat each positive time t,no matter how unsmooth the initial temperature proﬁle . We have discovered the basic smoothing property of heat ﬂow. Theorem 14.2. Ifu(t,x)isasolutiontotheheatequationwithpiecewisecontinuous †To be honest, we are using periodic boundary conditions in the ﬁgures , although the Dirichlet version leads to similar results 12/11/12 747 c/ci∇cleco√y∇t2012 Peter J. Olver initial data f(x) =u(0,x), or, more generally, initial data satisfying (14.23), then, for any t >0, the solution u(t,x)is an inﬁnitely diﬀerentiable function of x. After even a very short amount of time, the heat equation smoo thes out most, and, eventually, all of the ﬂuctuations in the initial temperatu re proﬁle. As a consequence, it becomes impossible to reconstruct the initial temperature u(0,x) =f(x) by measuring the temperaturedistribution h(x) =u(t,x)atalatertime t >0.Diﬀusion is irreversible —we cannot run the heat equation backwards in time! Indeed, if th e initial data u(0,x) =f(x) is not smooth, there is nofunction u(t,x) fort <0 that could possibly yield such an initial distribution because all corners and singularitie s are smoothed out by the diﬀusion process as tgoes forward! Or, to put it another way, the Fourier coeﬃcien ts (14.21) of any purported solution will be exponentially growing when t <0, and so high frequency noise will completely overwhelm the solution. For this reason, th e backwards heat equation is said to be ill-posed. On the other hand, the unsmoothing eﬀect of the backwards hea t equation does have potentialbeneﬁts. Forexample, inimageprocessing, diﬀus ionwillgraduallybluranimage. Image enhancement is the reverse process, and requires runn ing the heat ﬂow backwards in some stable manner. One option is to restrict to the backwa rds evolution to the ﬁrst few Fourier modes, which prevents the small scale ﬂuctuations f rom overwhelming the com- putation. Similar issues arise in the reconstruction of sub terranean proﬁles from seismic data, a problem of great concern in the oil and gas industry. I n forensics, determining the time of death based on the current temperature of a corpse also requires running the equations governing the dissipation of body heat backwards in time. For these and other applications, a key issue in contemporary research is how to cleverly circumventing the ill-posedness of the backwards heat ﬂow. Remark: The irreversibility of the heat equation points out a cruci al distinction be- tween partial diﬀerential equations and ordinary diﬀerent ial equations. Ordinary diﬀer- ential equations are always reversible — unlike the heat equ ation, existence, uniqueness and continuous dependence properties of solutions are all e qually valid in reverse time (although the detailed qualitative and quantitative prope rties of solutions can very well depend upon whether timeis running forwards or backwards). The irreversibility of partial diﬀerential equations modeling the diﬀusive processes in o ur universe may well be why, in our experience, Time’s Arrow points exclusively to the futu re. The Heated Ring Let us next consider the periodic boundary value problem mod eling heat ﬂow in an insulated circular ring. Let us ﬁx the length of the ring to be ℓ= 2π, with−π < x < π representing “angular” coordinate around the ring. For sim plicity, we also choose units in which the thermal diﬀusivity is γ= 1. Thus, we seek to solve the heat equation ∂u ∂t=∂2u ∂x2,−π < x < π, t > 0, (14.25) subject to periodic boundary conditions u(t,−π) =u(t,π),∂u ∂x(t,−π) =∂u ∂x(t,π), t ≥0.(14.26) 12/11/12 748 c/ci∇cleco√y∇t2012 Peter J. Olver The initial temperature distribution is u(0,x) =f(x),−π < x < π. (14.27) The resulting temperature u(t,x) will be a periodic function in xof period 2 π. Substituting the separable solution ansatz u(t,x) =e−λtv(x) into the heat equation and the boundary conditions leads to the periodic eigenvalu e problem d2v dx2+λv= 0, v (−π) =v(π), v′(−π) =v′(π). (14.28) As we know, in this case the eigenvalues are λn=n2wheren= 0,1,2,...is a non-negative integer, and the corresponding eigenfunction solutions ar e the trigonometric functions vn(x) = cosnx, /tildewidevn(x) = sinnx, n = 0,1,2,... . Note that λ0= 0 is a simple eigenvalue, with constant eigenfunction cos0 x= 1 — the sine solution sin0 x≡0 is trivial — while the positive eigenvalues are, in fact, do uble, each possessing two linearly independent eigenfunctions. The c orresponding separable solutions to the heated ring equation are un(t,x) =e−n2tcosnx, /tildewideun(t,x) =e−n2tsinnx, n = 0,1,2,3,... . The resulting inﬁnite series solution is u(t,x) =1 2a0+∞/summationdisplay n=1/parenleftbig ane−n2tcosnx+bne−n2tsinnx/parenrightbig . (14.29) The initial conditions require u(0,x) =1 2a0+∞/summationdisplay n=1/parenleftbig ancosnx+bnsinnx/parenrightbig =f(x), (14.30) which is precisely the Fourier series of the initial tempera ture proﬁle f(x). Consequently, an=1 π/integraldisplayπ −πf(x)cosnxdx, bn=1 π/integraldisplayπ −πf(x)sinnxdx, (14.31) are the usual Fourier coeﬃcients of f(x). As in the Dirichlet problem, after the initial instant, the h igh frequency terms in the series (14.29) become extremely small, since e−n2t≪1 forn≫0. Therefore, as soon as t >0, the solution essentially degenerates into a ﬁnite sum ove r the ﬁrst few Fourier modes. Moreover, as t→ ∞,allof the Fourier modes will decay to zero with the exception of the constant one, with null eigenvalue λ0= 0. Therefore, the solution will converge exponentially fast to a constant temperature proﬁ le: u(t,x)−→1 2a0=1 2π/integraldisplayπ −πf(x)dx, 12/11/12 749 c/ci∇cleco√y∇t2012 Peter J. Olver which equals the averageof the initial temperature proﬁle. Physically, we observe t hat the heat energy is redistributed so that the ring achieves a unif orm constant temperature and is in thermal equilibrium. Indeed, the total heat H(t) =/integraldisplayπ −πu(t,x)dx= constant (14 .32) is conserved, meaning constant, for all time; the proof of th is fact is left as an Exercise . (On the other hand, the Dirichlet boundary value problem doe s not conserve hear energy.) Prior to equilibrium, only the lowest frequency Fourier mod es will still be noticeable, and so the solution will asymptotically look like u(t,x)≈1 2a0+e−t(a1cosx+b1sinx) =1 2a0+r1e−tcos(x+δ1), (14.33) where a1=r1cosδ1=1 2π/integraldisplayπ −πf(x)cosxdx, b1=r1sinδ1=1 2π/integraldisplayπ −πf(x)sinxdx. Thus, for most initial data, the solution approaches therma l equilibrium exponentially fast, at a unit rate. The exceptions are when r1=/radicalbig a2 1+b2 1= 0, for which the rate of convergence is even faster, namely at a rate e−k2twherekis the smallest integer such that rk=/radicalbig a2 k+b2 k/ne}ationslash= 0. Inhomogeneous Boundary Conditions So far, we have concentrated our attention on homogeneous bo undary conditions. There is a simple trick that will convert a boundary value pro blem with inhomogeneous but constant Dirichlet boundary conditions, u(t,0) =α, u (t,ℓ) =β, t ≥0, (14.34) into a homogeneous Dirichlet problem. We begin by solving fo r the equilibrium tempera- ture proﬁle, which is the aﬃne function that satisﬁes the two boundary conditions, namely u⋆(x) =α+β−α ℓx. (14.35) The diﬀerence /tildewideu(t,x) =u(t,x)−u⋆(x) =u(t,x)−α−β−α ℓx (14.36) measures the deviation of the solution from equilibrium. It clearly satisﬁes the homoge- neous boundary conditions at both ends: /tildewideu(t,0) = 0 = /tildewideu(t,ℓ). Moreover, by linearity, since both u(t,x) andu⋆(x) are solutions to the heat equation, so is/tildewideu(t,x). The initial data must be similarly adapted: /tildewideu(0,x) =/tildewidef(x) =f(x)−u⋆(x) =f(x)−α−β−α ℓx. 12/11/12 750 c/ci∇cleco√y∇t2012 Peter J. Olver Solving the resulting homogeneous initial value problem, w e write/tildewideu(t,x) in Fourier series form (14.19), where the Fourier coeﬃcients are computed fro m the modiﬁed initial data /tildewidef(x). The solution to the inhomogeneous boundary value problem thus has the series form u(t,x) =α+β−α ℓx+∞/summationdisplay n=1/tildewidebnexp/parenleftbigg −γn2π2 ℓ2t/parenrightbigg sinnπx ℓ, (14.37) where /tildewidebn=2 ℓ/integraldisplayℓ 0/tildewidef(x) sinnπx ℓdx, n = 1,2,3,... . (14.38) Since, for any reasonable initial data, /tildewideu(t,0) will decay to zero at an exponential rate as t→ ∞, the actual temperature proﬁle (14.37) will asymptoticall y decay to the equilibrium proﬁle, u(t,x)−→u⋆(x) =α+β−α ℓx at the same exponentially fast rate. This method does not apply when the boundary conditions are t ime-dependent: u(t,0) =α(t), u(t,ℓ) =β(t). Attempting to mimic the preceding technique, we discover that the deviation /tildewideu(t,x) =u(t,x)−u⋆(t,x),where u⋆(t,x) =α(t)+β(t)−α(t) ℓx,(14.39) does satisfy the homogeneous boundary conditions, but now s olves an inhomogeneous version of the heat equation: ∂/tildewideu ∂t=∂2/tildewideu ∂x2−h(t,x),where h(t,x) =∂u⋆ ∂t(t,x). (14.40) Solution techniques in this case will be discussed below. 14.2. Symmetry and the Maximum Principle. So far we have relied on the method of separation of variables to construct explicit solutions to partial diﬀerential equations. A second usefu l solution technique relies on exploiting inherent symmetry properties of the diﬀerentia l equation. Unlike separation of variables†, symmetry methods can be also successfully applied to produ ce solutions to a broad range of nonlinear partial diﬀerential equations ; examples can be found in Chapter 22. While we do not have the space or required mathema tical tools to develop the full apparatus of symmetry techniques, we can introduce the important concept of a similarity solution , applied in the particular context of the heat equation. In general, by a symmetry of an equation, we mean a transformation, either linear (as in Section 7.2), aﬃne (as in Section 7.3), or even nonline ar, that takes solutions to solutions. Thus, if we have a symmetry, and know one solution , then we can construct a †This is not quite fair: separation of variables can be applied to some spec ial nonlinear partial diﬀerential equations such as Hamilton–Jacobi equations, [ 131]. 12/11/12 751 c/ci∇cleco√y∇t2012 Peter J. Olver second solution by applying the symmetry. And, possibly, a t hird solution by applying the symmetry yet again. And so on. If we know lots of symmetries, t hen we can produce lots and lots of solutions by this simple device. Remark: General symmetry techniques are founded on the theory of Li e groups, named after the inﬂuential nineteenth century Norwegian ma thematician Sophus Lie (pro- nounced “Lee”). Lie’s theory provides an algorithm for comp letely determining all the symmetries of a given diﬀerential equation, but this is beyo nd the scope of this introduc- tory text. However, direct inspection and/or physical intu ition will often detect the most important symmetries without appealing to such a sophistic ated theory. Modern appli- cations of Lie’s symmetry methods to partial diﬀerential eq uations arising in physics and engineering can be traced back to the inﬂuential book of G. Bi rkhoﬀ, [17], on hydrody- namics. A complete and comprehensive treatment of symmetry methods can be found in the ﬁrst author’s book [ 146], and, at a more introductory level, in the recent books by Hydon, [ 106], and Cantwell, [ 39], with particular emphasis on ﬂuid mechanics. The heat equation serves as an excellent testing ground for t he general symmetry methodology, as it admits a rich variety of symmetry transfo rmations that take solutions to solutions. The simplest are the translations. Moving the space and time coordinates by a ﬁxed amount, t/ma√sto−→t−a, x /ma√sto−→x−b, (14.41) wherea,bare constants, changes the function u(t,x) into the translated function U(t,x) =u(t−a,x−b). (14.42) A simple application of the chain rule proves that the partia l derivatives of Uwith respect totandxagree with the corresponding partial derivatives of u, so ∂U ∂t=∂u ∂t,∂U ∂x=∂u ∂x,∂2U ∂x2=∂2u ∂x2, and so on. In particular, the function U(t,x) is a solution to the heat equation Ut=γUxx whenever u(t,x) also solves ut=γuxx. Physically, the translation symmetries formalize theproperty thattheheatequationmodelsahomogeneous med ium, andhence thesolution does not depend on the choice of reference point or origin of o ur coordinate system. As a consequence, each solution to the heat equation will pro duce an inﬁnite family of translated solutions. For example, starting with the sep arable solution u(t,x) =e−γtsinx, we immediately produce the additional solutions u(t,x) =e−γ(t−a)sinπ(x−b), valid for any choice of constants a,b. Warning : Typically, the symmetries of a diﬀerential equation do not respect initial or boundary conditions. For instance, if u(t,x) is deﬁned for t≥0 and in the domain 0≤x≤ℓ, then its translated version U(t,x) is deﬁned for t≥aand in the translated domainb≤x≤ℓ+b, and so will solve a translated initial-boundary value prob lem. 12/11/12 752 c/ci∇cleco√y∇t2012 Peter J. Olver A second, even more important class of symmetries are the sca ling invariances. We already know that if u(t,x) is a solution, so is any scalar multiple cu(t,x); this is a simple consequence of linearity of the heat equation. We can also ad d an arbitrary constant to the temperature, noting that U(t,x) =cu(t,x)+k (14.43) isasolutionforanychoiceofconstants c,k. Physically,thetransformation(14.43)amounts to a change in the scale for measuring temperature. For insta nce, ifuis measured degrees Celsius, and we set c=9 5andk= 32, then U=9 5u+ 32 will be measured in degrees Fahrenheit. Thus, reassuringly, the physical processes de scribed by the heat equation do not depend upon our choice of thermometer. More interestingly, suppose we rescale the space and time va riables: t/ma√sto−→αt, x /ma√sto−→βx, (14.44) whereα,β >0 are positive constants. The eﬀect of such a scaling transfo rmation is to convertu(t,x) into a rescaled function U(t,x) =u(αt,βx). (14.45) The derivatives of Uare related to those of uaccording to the following formulae, which are direct consequences of the multi-variable chain rule: ∂U ∂t=α∂u ∂t,∂U ∂x=β∂u ∂x,∂2U ∂x2=β2∂2u ∂x2. Therefore, if usatisﬁes the heat equation ut=γuxx, thenUsatisﬁes the rescaled heat equation Ut=αut=αγuxx=αγ β2Uxx, which we rewrite as Ut= ΓUxx,where Γ =γα β2. (14.46) Thus, the net eﬀect of scaling space and time is merely to resc ale the diﬀusion coeﬃcient in the heat equation. Physically, the scaling symmetry (14.44 ) corresponds to a change in the physical units used to measure time and distance. For instan ce, to change from seconds to minutes, set α= 60, and from meters to yards, set β= 1.0936. The net eﬀect (14.46) on the diﬀusion coeﬃcient is a reﬂection of its physical units, namely distance2/time. In particular, if we choose α=1 γ, β = 1, then the rescaled diﬀusion coeﬃcient becomes Γ = 1. This obse rvation has the following important consequence. If U(t,x) solves the heat equation for a unit diﬀusivity, Γ = 1, then u(t,x) =U(γt,x) (14 .47) 12/11/12 753 c/ci∇cleco√y∇t2012 Peter J. Olver solvestheheatequationforthediﬀusivity γ. Thus, theonlyeﬀectofthediﬀusioncoeﬃcient γis to speed up or slow down time! A body with diﬀusivity γ= 2 will cool down twice as fast as a body (of the same shape subject to the same boundar y conditions and initial conditions) with diﬀusivity γ= 1. Note that this particular rescaling has not altered the space coordinates, and so U(t,x) is deﬁned on the same domain as u(t,x). On the other hand, if we set α=β2, then the rescaled diﬀusion coeﬃcient is exactly the same as the original: Γ = γ. Thus, the transformation t/ma√sto−→β2t, x /ma√sto−→βx, (14.48) does not alter the equation, and hence deﬁnes a scaling symmetry , also known as a simi- larity transformation , for the heat equation. Combining (14.48) with the linear re scaling u/ma√sto→cu, we make the elementary, but important observation that if u(t,x) is any solution to the heat equation, then so is the function U(t,x) =cu(β2t,βx), (14.49) forthesamediﬀusioncoeﬃcient γ. Forexample,rescalingthesolution u(t,x) =e−γt2cosx leads to the solution U(t,x) =ce−γβ2t2cosβx. Warning : As in the case of translations, rescaling space by a factor β/ne}ationslash= 1 will alter the domain of deﬁnition of the solution. If u(t,x) is deﬁned for 0 ≤x≤ℓ, thenU(t,x) is deﬁned for 0 ≤x≤ℓ/β. Suppose that we have solved the heat equation for the tempera tureu(t,x) on a bar of length 1, subject to certain initial and boundary conditi ons. We are then given a bar composedofthesamematerialoflength2. Sincethediﬀusivi tycoeﬃcient hasnotchanged, and we can directly construct the new solution U(t,x) by rescaling. Setting β=1 2will serve to double the length. If we also rescale time by a factor α=β2=1 4, then the rescaled function U(t,x) =u/parenleftbig1 4t,1 2x/parenrightbig will be a solution of the heat equation on the longer bar with the same diﬀusivity constant. The net eﬀect is that t he rescaled solution will be evolving four times as slowly as the original. Thus, it eﬀect ively takes a bar that is twice the length four times as long to cool down. The Maximum Principle TheMaximum Principle is the mathematical formulation of the thermodynamical law that heat cannot, in the absence of external sources, achiev e a value at a point that strictly larger than its surroundings. We will help in the proof to formulate it for the more general c ase when the heat equation is subject to external forcing. Theorem 14.3. Letγ >0. Suppose u(t,x)is a solution to the forced heat equation ∂u ∂t=γ∂2u ∂x2+F(t,x),on the rectangular domain R={a≤x≤b,0≤t≤c.} Suppose F(t,x)≤0for all(t,x)∈R. Then the global maximum of u(t,x)on the domain Roccurs either at t= 0or atx=aorx=b. 12/11/12 754 c/ci∇cleco√y∇t2012 Peter J. Olver In other words, if we are not adding in any external heat, the m aximum temperature ion the bar occurs either at the initial time, or at one of the e ndpoints. Proof: First let us prove the result under the assumption that F(t,x)<0, and hence ∂u ∂t< γ∂2u ∂x2(14.50) everywhere in R. Suppose u(t,x) has a (local) maximum at a point in the interior of R. Then, by multivariable calculus (see also Theorem 19.43), its gradient must vanish, sout=ux= 0, and its Hessian matrix must be positive deﬁnite, which, i n particular, requires that uxx≤0. But these contradict the inequality (14.50). Thus, the so lution cannot have a local interior maximum. If the maximum were to o ccur on the top of the rectangle, at a point ( t,x) wheret=c, then we would necessarily have ut≥0 there, as otherwise u(t,x) would be decreasing as a function of tand hence ( c,x) could not be a local maximum on R, and also uxx≤0, again contradicting (14.50). To generalize to the case when F(t,x)≤0 — which includes the heat equation when F(t,x)≡0, requires a little trick. We set v(t,x) =u(t,x)+εx2,where ε >0. Then, ∂v ∂t=γ∂2v ∂x2−2γε+F(t,x) =γ∂2v ∂x∂+/tildewideF(t,x), where /tildewideF(t,x) =F(t,x)−2γε <0 everywhere in R. Thus, by the previous paragraph, the maximum of voccurs on either the bottom or sides of the rectangle. Now we let ε→0 and conclude the same for u. More precisely, let u(t,x)≤Mont= 0 orx=aorx=b. Then v(t,x)≤M+εmax{a2,b2} and hence, on all of R, u(t,x)≤v(t,x)≤M+εmax{a2,b2}. lettingε→0 proves that u(t,x)≤Meverywhere, which completes the proof. Q.E.D. Thus, solution u(t,x) to the heat equation (14.12) can only achieve a maximum on th e boundary of its domain of deﬁnition. In other words, a soluti on cannot have a maximum at any point ( t,x) that lies in the interior of the domain and after the initial timet >0. Physically, if the temperature in a fully insulated bar star ts out everywhere above freezing, then, in the absence of external heat sources, it can never di p below freezing at any later time. Let us apply the Maximum Principle to prove uniqueness of sol utions to the heat equation. Theorem 14.4. There is at most one solution to the boundary value problem. 12/11/12 755 c/ci∇cleco√y∇t2012 Peter J. Olver Proof: Suppose uand/tildewideuare any two solutions. Then their diﬀerence v=u−/tildewideusolves the homogeneous boundary value problem. Thus, by the maximu principle v(t,x)≤0 at all points of R. But−v=/tildewideu−ualso solves the homogeneous boundary value problem, and hence −v≤0 too. This implies v(t,x)≡0 and hence u≡/tildewideueverywhbere, proving uniqueness. Q.E.D. Remark: Existence follows from the Fourier series solution — assum ing the initialand boundary data and the forcing function are suﬃciently nice. 14.3. The Fundamental Solution. One disadvantage of the Fourier series solution to the heat e quation is that it is not nearly as explicit as one might desire for either practical a pplications, numerical computa- tions, or even further theoretical investigations and deve lopments. An alternative, general approach is based on the idea of the fundamental solution , which derives its inspiration from the Green’s function method for solving boundary value problems. For the heat equation, the fundamental solution measures the eﬀect of a c oncentrated heat source. Let us restrict our attention to homogeneous boundary condi tions. The idea is to analyze the case when the initial data u(0,x) =δy(x) =δ(x−y) is a delta function, which we can interpret as a highly concentrated unit heat source, e .g., a soldering iron or laser beam, that is instantaneously applied at a position yalong the bar. The heat will diﬀuse away from its initial concentration, and the resulting fundamental solution is denoted by u(t,x) =F(t,x;y),with F(0,x;y)=δ(x−y). (14.51) For each ﬁxed y, the fundamental solution F(t,x;y), considered as a function of t >0 and x, must satisfy the diﬀerential equation, so ∂F ∂t=γ∂2F ∂x2, (14.52) as well as the speciﬁed homogeneous boundary conditions. Once we have determined the fundamental solution, we can the n use linear superpo- sition to reconstruct the general solution to the initial-b oundary value problem. Namely, we ﬁrst write the initial data u(0,x) =f(x) =/integraldisplayℓ 0δ(x−y)f(y)dy (14.53) as a superposition of delta functions, as in (11.37). Linear ity implies that the solution is then the corresponding superposition of the responses to th ose concentrated delta proﬁles: u(t,x) =/integraldisplayℓ 0F(t,x;y)f(y)dy. (14.54) Assuming that we can diﬀerentiate under the integral sign, t he fact that F(t,x;y) satis- ﬁes the diﬀerential equation and the homogeneous boundary c onditions for each ﬁxed y immediately implies that the integral (14.54) is also a solu tion with the correct initial and (homogeneous) boundary conditions. 12/11/12 756 c/ci∇cleco√y∇t2012 Peter J. Olver Unfortunately, most boundary value problems do not have fun damental solutions that can be written down in closed form. An important exception is the case of an inﬁnitely long homogeneous bar, which requires solving the heat equat ion ∂u ∂t=∂2u ∂x2,for −∞< x <∞, t > 0. (14.55) For simplicity, we have chosen units in which the thermal diﬀ usivity is γ= 1. The solution u(t,x) is deﬁned for all x∈R, and has initial conditions u(0,x) =f(x) for −∞< x <∞. In order to specify the solution uniquely, we shall require t hat the temperature be square- integrable at all times, so that /integraldisplay∞ −∞|u(t,x)|2dx <∞ for all t≥0. (14.56) Roughly speaking, we are requiring that the temperature be s mall at largedistances, which are the relevant boundary conditions for this situation. On an inﬁnite interval, the Fourier series solution to the he at equation becomes a Fourier integral. We write the initial temperature distrib ution as a superposition f(x) =1√ 2π/integraldisplay∞ −∞eikx/hatwidef(k)dk, of complex exponentials eikx, where /hatwidef(k) is the Fourier transform (13.84) of f(x). The corresponding separable solutions to the heat equation are u(t,x) =e−k2teikx=e−k2t/parenleftbig coskx+ i sinkx/parenrightbig , (14.57) where the frequency variable kis allowed to assume any real value. We invoke linear superposition to combine these complex solutions into a Fou rier integral u(t,x) =1√ 2π/integraldisplay∞ −∞e−k2teikx/hatwidef(k)dk (14.58) that forms the solution to the initial value problem for the h eat equation. In particular, to recover the fundamental solution, we take the initial temperature proﬁle to be a delta function δy(x) =δ(x−y) concentrated at x=y. According to (13.114), its Fourier transform is /hatwideδy(k) =e−iky √ 2π. Plugging this into (14.58), and then referring to our table o f Fourier transforms, we ﬁnd the following explicit formula for the fundamental solutio n: F(t,x;y) =1 2π/integraldisplay∞ −∞e−k2teik(x−y)dk=1 2√ πte−(x−y)2/(4t).(14.59) 12/11/12 757 c/ci∇cleco√y∇t2012 Peter J. Olver -6 -4 -2 2 4 60.20.40.60.811.2 -6 -4 -2 2 4 60.20.40.60.811.2 -6 -4 -2 2 4 60.20.40.60.811.2 -6 -4 -2 2 4 60.20.40.60.811.2 Figure 14.3. The Fundamental Solution to the Heat Equation. As you can verify, for each ﬁxed y, the function F(t,x;y) is indeed a solution to the heat equation for all t >0. In addition, lim t→0+F(t,x;y) =/braceleftbigg0, x/ne}ationslash=y, ∞, x=y. Furthermore, its integral/integraldisplay∞ −∞F(t,x;y)dx= 1, (14.60) which represents the total heat energy, is constant — in acco rdance with the law of conser- vation of energy; see Exercise . Therefore, as t→0+, the fundamental solution satisﬁes the original limiting deﬁnition (11.29–30) of the delta fun ction, and so F(0,x;y) =δy(x) has the desired initial temperature proﬁle. In Figure 14.3 w e graph F(t,x;0) at times t=.05,.1,1.,10.. It starts life as a delta spike concentrated at the origin, a nd then im- mediately smoothes out into a tall and narrow bell-shaped cu rve, centered at x= 0. As time increases, the solution shrinks and widens, decaying e verywhere to zero. Its maximal amplitude is proportional to t−1/2, while its overall width is proportional to t1/2. The total heat energy (14.60), which is the area under the graph, remains ﬁxed while gradually spreading out over the entire real line. Remark: In probability, these exponentially bell-shaped curves a re known as normal orGaussian distributions . The width of the bell curve corresponds to the standard devi- ation. For this reason, the fundamental solution to the heat equat ion sometimes referred to as a “Gaussian ﬁlter”. Remark: One of the non-physical artifacts of the heat equation is th at the heat energy propagateswithinﬁnitespeed. Indeed, the eﬀect ofany init ialconcentrationof heat energy will immediately be felt along the entire length of an inﬁnit e bar, because, at any t >0, the fundamental solution is nonzero for all x. (The graphs in Figure 14.3 are a little misleading because they fail to show the extremely small, bu t still positive, exponentially decreasing tails.) Thiseﬀect, while moreor less negligibl eat largedistances, isnevertheless in clear violation of physical intuition — not to mention rel ativity that postulates that signals cannot propagate faster than the speed of light. Des pite this non-physical property, 12/11/12 758 c/ci∇cleco√y∇t2012 Peter J. Olver the heat equation remains an extremely accurate model for he at propagation and similar diﬀusive phenomena. With the fundamental solution in hand, we can then adapt the l inear superposition formula (14.54) to reconstruct the general solution u(t,x) =1 2√ πt/integraldisplay∞ −∞e−(x−y)2/(4t)f(y)dy (14.61) to our initial value problem (14.55). Comparing with (13.12 8), we see that the solutions are obtained by convolution, u(t,x) =g(t,x)∗f(x),where g(t,x) =F(t,x;0) =e−x2/(4t), of the initial data with a one-parameter family of progressi vely wider and shorter Gaussian ﬁlters. Since u(t,x) solves the heat equation, we conclude that Gaussian ﬁlter c onvolution has the same smoothing eﬀect on the initial signal f(x). Indeed, the convolution integral (14.61) serves to replace each initial value f(x) by a weighted average of nearby values, the weight being determined by the Gaussian distribution. S uch a weighted averaging has the eﬀect of smoothing out high frequency variations in the s ignal, and, consequently, the Gaussian convolution formula (14.61) provides an eﬀective method for denoising signals and images. In fact, for practical reasons, the graphs displ ayed earlier in Figure 14.2 were computed by using a standard numerical integration routine to evaluate the convolution integral (14.61), rather than by a numerical solution schem e for the heat equation. Example 14.5. An inﬁnite bar is initially heated to unit temperature along a ﬁnite interval. This corresponds to an initial temperature proﬁl e u(0,x) =f(x) =σ(x−a)−σ(x−b) =/braceleftbigg1, a < x < b, 0,otherwise . Thecorresponding solutiontotheheatequationisobtained bytheintegralformula(14.61), producing u(t,x) =1 2√ πt/integraldisplayb ae−(x−y)2/(4t)dy=1 2/bracketleftbigg erf/parenleftbiggx−a 2√ t/parenrightbigg −erf/parenleftbiggx−b 2√ t/parenrightbigg/bracketrightbigg ,(14.62) where erfx=2√π/integraldisplayx 0e−z2dz (14.63) is known as the error function due to its applications in probability and statistics, [ 66]. A graph appears in Figure 14.4. The error function integral c annot be written in terms of elementary functions. Nevertheless, its importance in v arious applications means that its properties have been well studied, and its values tabula ted, [145]. In particular, it has asymptotic values lim x→∞erfx= 1, lim x→−∞erfx=−1.(14.64) A graph of the heat equation solution (14.62) when a=−5,b= 5, at successive times t= 0.,.1,1,5,30,300, is displayed in Figure 14.5. Note the initial smoothing or blurring of the sharp interface, followed by a gradual decay to therma l equilibrium. 12/11/12 759 c/ci∇cleco√y∇t2012 Peter J. Olver -2 -1 1 2 -1-0.50.51 Figure 14.4. The Error Function. -10 -5 5 100.20.40.60.81 -10 -5 5 100.20.40.60.81 -10 -5 5 100.20.40.60.81 -10 -5 5 100.20.40.60.81 -10 -5 5 100.20.40.60.81 -10 -5 5 100.20.40.60.81 Figure 14.5. Error Function Solution to the Heat Equation. The Forced Heat Equation The fundamental solution can be also used to solve the inhomo geneous heat equation ut=uxx+h(t,x), (14.65) that models a bar under an external heat source h(t,x), that might depend upon both position and time. We begin by solving the particular case ut=uxx+δ(t−s)δ(x−y), (14.66) whose inhomogeneity represents a heat source of unit magnit ude that is concentrated at a position 0 < y < ℓ and applied instantaneously at a single time t=s >0. Physically, we apply a soldering iron or laser beam to a single spot on the bar for a brief moment. Let us also impose homogeneous initial conditions u(0,x) = 0 (14 .67) as well as homogeneous boundary conditions of one of our stan dard types. The resulting solution u(t,x) =G(t,x;s,y) (14 .68) 12/11/12 760 c/ci∇cleco√y∇t2012 Peter J. Olver will be referred to as the general fundamental solution to the heat equation. Since a heat source which is applied at time swill only aﬀect the solution at later times t≥s, we expect that G(t,x;s,y)= 0 for all t < s. (14.69) Indeed, since u(t,x) solves the unforced heat equation at all times t < ssubject to homo- geneous boundary conditions and has zero initial temperatu re, this follows immediately from the uniqueness of the solution to the initial-boundary value problem. Once we know the general fundamental solution (14.68), we ar e able to solve the problem for a general external heat source (14.65) by appeal ing to linearity. We ﬁrst write the forcing as a superposition h(t,x) =/integraldisplay∞ 0/integraldisplayℓ 0h(s,y)δ(t−s)δ(x−y)dyds (14.70) of concentrated instantaneous heat sources. Linearity all ows us to conclude that the solu- tion is given by the self-same superposition formula u(t,x) =/integraldisplayt 0/integraldisplayℓ 0h(s,y)G(t,x;s,y)dyds. (14.71) The fact that we only need to integrate over times 0 ≤s≤tfollows from (14.69). Remark: If we have a nonzero initial condition, u(0,x) =f(x), then we appeal to linear superposition to write the solution u(t,x) =/integraldisplayℓ 0F(t,x;y)f(y)dy+/integraldisplayt 0/integraldisplayℓ 0h(s,y)G(t,x;s,y)dyds (14.72) as a combination of ( a) the solution with no external heat source, but nonzero init ial conditions, plus ( b) the solution with homogeneous initial conditions but nonz ero heat source. Let us solve the forced heat equation in the case of a inﬁnite b ar, so−∞< x <∞. We begin by computing the general fundamental solution to (1 4.66), (14.67). As before, we take the Fourier transform of both sides of the partial diﬀ erential equation with respect tox. In view of (13.114), (13.118), we ﬁnd ∂/hatwideu ∂t+k2/hatwideu=1√ 2πe−ikyδ(t−s). (14.73) which is an inhomogeneous ﬁrst order ordinary diﬀerential e quation for the Fourier trans- form/hatwideu(t,k) ofu(t,x). Assuming s >0, by (14.69), the initial condition is /hatwideu(0,k) = 0. (14.74) We solve the initial value problem (14.73–74) by the usual me thod, [26]. Multiplying the diﬀerential equation by the integrating factor ek2tyields ∂ ∂t/parenleftig ek2t/hatwideu/parenrightig =1√ 2πek2t−ikyδ(t−s), 12/11/12 761 c/ci∇cleco√y∇t2012 Peter J. Olver Integrating both sides from 0 to tand using the initial condition, we ﬁnd /hatwideu(t,k) =1√ 2πek2(t−s)−ikyσ(t−s), whereσ(t) is the usual step function (11.43). Finally, we apply the in verse Fourier trans- form formula (13.87), and then (14.59), to deduce that u(t,x) =G(t,x;s,y) =σ(t−s) 2π/integraldisplay∞ −∞eik(y−x)+k2(t−s)dk =σ(t−s) 2/radicalbig π(t−s)exp/bracketleftbigg −(x−y)2 4(t−s)/bracketrightbigg =σ(t−s)F(t−s,x−y). Thus, thegeneral fundamental solutionisobtainedby trans latingthefundamental solution F(t,x;y) for the initial value problem to a starting time of t=sinstead of t= 0. Thus, an initial condition has the same aftereﬀect on the temperatur e as an instantaneous applied heat source of the same magnitude. Finally, the superpositi on principle (14.71) produces the solution u(t,x) =/integraldisplayt 0/integraldisplay∞ −∞h(s,y) 2/radicalbig π(t−s)exp/bracketleftbigg −(x−y)2 4(t−s)/bracketrightbigg dyds. (14.75) to the heat equation with source term on an inﬁnite bar. The Root Cellar Problem As a ﬁnal example, we discuss a problem that involves analysi s of the heat equation on a semi-inﬁnite interval. The question is: how deep should you dig a root cellar? In the prerefrigeration era, a root cellar was used to keep food coo l in the summer, but not freeze in the winter. We assume that the temperature in the earth onl y 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, ﬂuctuates in a periodic manner; speciﬁcally, we set u(t,0) =acosωt, (14.76) where the oscillatory frequency ω=2π 365.25 days= 2.0×10−7sec−1(14.77) refers to yearly temperature variations. In this model, we s hall ignore daily temperature ﬂuctuations as their eﬀect is not signiﬁcant below a very thi n surface layer. At large depth the temperature is assumed to be unvarying: u(t,x)−→0 as x−→ ∞, (14.78) where 0 refers to the mean temperature. Thus, we must solve the heat equation on a semi-inﬁnite bar 0 < x <∞, with time- dependent boundary conditions(14.76),(14.78)at the ends . The analysiswill besimpliﬁed 12/11/12 762 c/ci∇cleco√y∇t2012 Peter J. Olver a littleifwereplace the cosineby a complex exponential, an dso lookfor a complex solution with boundary conditions u(t,0) =aeiωt, lim x→∞u(t,x) = 0. (14.79) Let us try a separable solution of the form u(t,x) =v(x)eiωt. (14.80) Substituting this expression into the heat equation ut=γuxxleads to iωv(x)eiωt=γv′′(x)eiωt. Canceling the common exponential factors, we conclude that v(x) should solve the bound- ary value problem γv′′(x) = iωv, v (0) =a, lim x→∞v(x) = 0. The solutions to the ordinary diﬀerential equation are v1(x) =e√ iω/γ x=e√ ω/2γ(1+i)x, v2(x) =e−√ iω/γ x=e−√ ω/2γ(1+i)x. The ﬁrst solution is exponentially growing as x→ ∞, and so not appropriate to our problem. The solution to the boundary value problem must the refore be a multiple, v(x) =ae−√ ω/2γ(1+i)x of the exponentially decaying solution. Substituting back into (14.80), we ﬁnd the (com- plex) solution to the root cellar problem to be u(t,x) =ae−x√ ω/2γeiω/parenleftbig t−√ ω/2γ/parenrightbig x. (14.81) The corresponding real solution is obtained by taking the re al part, u(t,x) =ae−x√ ω/2γcos/parenleftbigg ωt−/radicalbiggω 2γx/parenrightbigg . (14.82) The ﬁrst term in (14.82) is exponentially decaying as a funct ion of the depth. Thus, the further downonegoes, thelessnoticeabletheeﬀect ofthesu rface temperatureﬂuctuations. The second term is periodic with the same annual frequency ω. The interesting feature is the phase lag in the response. The temperature at depth xis out of phase with respect to the surface temperature ﬂuctuations, having an overall pha se lag δ=/radicalbiggω 2γx that depends linearly on depth. In particular, a cellar buil t at a depth where δis an odd multiple of πwill be completely out of phase, being hottest in the winter, and coldest in 12/11/12 763 c/ci∇cleco√y∇t2012 Peter J. Olver the summer. Thus, the (shallowest) ideal depth at which to bu ild a root cellar would take δ=π, corresponding to a depth of x=π/radicalbigg 2γ ω. For typical soils in the earth, γ≈10−6meters2sec−1, and hence, by (14.77), x≈ 9.9 meters. However, at this depth, the relative amplitude of t he oscillations is e−x√ ω/2γ=e−π=.04 and hence there is only a 4% temperature ﬂuctuation. In Minne sota, 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 th e warmest, to the summer, where it is the coldest. Building the cellar twice as deep would lead to a temperature ﬂuctuation of .2%, now in phase with the surface variations, which means that the cellar is, for all practical purposes, at cons tant temperature year round. 14.4. The Wave Equation. The second important class of dynamical partial diﬀerentia l equations are those mod- eling vibrations of continuous media. As we saw in Chapter 9, Newton’s Law implies that the free vibrations of a discrete mechanical system are gove rned by a second order system of ordinary diﬀerential equations of the form Md2u dt2=−Ku, in which Mis the positive deﬁnite, diagonal mass matrix, while K=A∗A=ATCAis the positive deﬁnite (or semi-deﬁnite in the case of an unstable system) stiﬀness matrix. The corresponding dynamical equations describing the smal l vibrations of continuous media take an entirely analogous form ρ∂2u ∂t2=−K[u]. (14.83) In this framework, ρdescribes the density, while K=L∗◦Lis the same self-adjoint dif- ferential operator, with appropriate boundary conditions , that appears in the equilibrium equations (11.90). For one-dimensional media, such as a vib rating bar or string, we are led to a partial diﬀerential equation in the particular form ρ(x)∂2u ∂t2=∂ ∂x/parenleftbigg κ(x)∂u ∂x/parenrightbigg , 0< x < ℓ, (14.84) whereρ(x) is the density of the bar or string at position x, whileκ(x)>0 denotes its stiﬀ- ness or tension. The second order partial diﬀerential equat ion (14.84)models the dynamics of vibrations and waves in a broad range of applications, inc luding elastic vibrations of a bar, sound vibrations in a column of air, e.g., inside a wind i nstrument, and also transverse 12/11/12 764 c/ci∇cleco√y∇t2012 Peter J. Olver vibrations of a string, e.g., a violin string. (However, ben ding vibrations of a beam lead to a fourth order partial diﬀerential equation; see Exercis e.) The wave equation is also used to model small amplitude water waves, electromagnetic waves, including light, radio and microwaves, gravitational waves, and many others. A det ailed derivation of the model from ﬁrst principles in the case of a vibrating string can be f ound in [ 187]. To specify the solution, we must impose suitable boundary co nditions. The usual suspects — Dirichlet, Neumann, mixed, and periodic boundar y conditions — continue to play a central role, and have immediate physical interpreta tions. Tying down an end of the string imposes a Dirichlet condition u(t,0) =α, while a free end is prescribed by a homogeneous Neumann boundary condition ux(t,0) = 0. Periodic boundary conditions, as in (14.11), correspond to the vibrations of a circular rin g. As with all second order New- tonian systems of ordinary diﬀerential equations, the solu tion to the full boundary value problem for the second order partial diﬀerential equation w ill then be uniquely speciﬁed by its initial displacement and initial velocity: u(0,x) =f(x),∂u ∂t(0,x) =g(x). (14.85) The simplest situation occurs when the medium is homogeneou s, and so both its density and stiﬀness are constant. Then the general vibrati on equation (14.84) reduces to the one-dimensional wave equation ∂2u ∂t2=c2∂2u ∂x2. (14.86) The constant c=/radicalbiggκ ρ>0 (14 .87) is known as the wave speed , for reasons that will soon become apparent. The method for solving such second order systems is motivate d by our solution in the discrete case discussed in Section 9.5. To keep matters simp le, we shall concentrate on the wave equation (14.86), although the method is easily ext ended to the general homoge- neous Newtonian system (14.84). We will try a separable solu tion with trigonometric time dependence: u(t,x) = cos(ωt)v(x), (14.88) in which both the frequency ωand proﬁle v(x) are to be determined. Diﬀerentiating (14.88), we ﬁnd ∂2u ∂t2=−ω2cos(ωt)v(x),∂2u ∂x2= cos(ωt)v′′(x). Substituting into the wave equation (14.86) and canceling t he common cosine factors, we deduce that v(x) must satisfy the ordinary diﬀerential equation c2d2v dx2+ω2v= 0. (14.89) Thusω2=λcan be viewed as an eigenvalue with eigenfunction v(x) for the second order diﬀerential operator K=−c2D2. Ignoring the boundary conditions for the moment, if 12/11/12 765 c/ci∇cleco√y∇t2012 Peter J. Olver ω >0, the solutions are the trigonometric functions cosωx c, sinωx c. Thus, (14.88) leads to the pair of explicit solutions cosωtcosωx c, cosωtsinωx c, to the wave equation. Now, the computation will work just as w ell with a sine function, which yields two additional solutions sinωtcosωx c, sinωtsinωx c. Each of these four solutions represents a spatially periodi c standing wave form of period 2πc/ω, that is vibrating with frequency ω. Observe that the smaller scale waves vibrate faster. On the other hand, if ω= 0, then (14.89) has the solution v=αx+β, leading to the solutions u(t,x) = 1, and u(t,x) =x. (14.90) The ﬁrst is a constant, nonvibrating solution, while the sec ond is also constant in time, but will typically not satisfy the boundary conditions and s o can be discarded. As we learned in Chapter 9, the existence of a zero eigenvalue corr esponds to an unstable mode in the physical system, in which the displacement grows line arly in time. In the present situation, these correspond to the two additional solution s u(t,x) =t, and u(t,x) =xt, (14.91) of the waveequation. Again, thesecond solution willtypica llynot satisfy the homogeneous boundary conditions, and can usually be ignored. Such null e igenfunction modes will only arise in unstable situations. The boundary conditions will distinguish the particular ei genvalues and natural fre- quencies of vibration. Consider ﬁrst the case of a string of l engthℓwith two ﬁxed ends, and thus subject to homogeneous Dirichlet boundary conditi ons u(t,0) = 0 = u(t,ℓ). This forms a positive deﬁnite boundary value problem, and so there is no unstable mode. Indeed, the eigenfunctions of the boundary value problem (1 4.89) with Dirichlet boundary conditions v(0) = 0 = v(ℓ) were found in (14.17): vn(x) = sinnπx ℓ, ωn=nπc ℓ, n = 1,2,3,... . Therefore, we can write the general solution as a Fourier sin e series u(t,x) =∞/summationdisplay n=1/bracketleftbigg bncosnπct ℓsinnπx ℓ+dnsinnπct ℓsinnπx ℓ/bracketrightbigg =∞/summationdisplay n=1rncos/parenleftbiggnπct ℓ+δn/parenrightbigg sinnπx ℓ.(14.92) 12/11/12 766 c/ci∇cleco√y∇t2012 Peter J. Olver The solution is thus a linear combination of the natural Four ier modes vibrating with frequencies ωn=nπc ℓ=nπ ℓ/radicalbiggκ ρ, n = 1,2,3,... . (14.93) The longer the length ℓof the string, or the higher its density ρ, the slower the vibrations; whereas increasing its stiﬀness or tension κspeeds them up — in exact accordance with our physical intuition. The Fourier coeﬃcients bnanddnin (14.92) will be uniquely determined by the initial conditions (14.85). Diﬀerentiating the series ter m by term, we discover that we must represent the initial displacement and velocity as Fou rier sine series u(0,x) =∞/summationdisplay n=1bnsinnπx ℓ=f(x),∂u ∂t(0,x) =∞/summationdisplay n=1dnnπc ℓsinnπx ℓ=g(x). Therefore, bn=2 ℓ/integraldisplayℓ 0f(x) sinnπx ℓdx, n = 1,2,3,... . (14.94) are the Fourier sine coeﬃcients (12.84) of the initial displ acement f(x), while dn=2 nπc/integraldisplayℓ 0g(x) sinnπx ℓdx, n = 1,2,3,... . (14.95) are rescaled versions of the Fourier sine coeﬃcients of the i nitial velocity g(x). Example 14.6. A string of unit length is held taut in the center and then rele ased. Our goal is to describe the ensuing vibrations. Let us assume the physical units are chosen so thatc2= 1, and so we are asked to solve the initial-boundary value pr oblem utt=uxx, u(0,x) =f(x), ut(0,x) = 0, u(t,0) =u(t,1) = 0.(14.96) To be speciﬁc, we assume that the center of the string has been displaced by half a unit, and so the initial displacement is f(x) =/braceleftigg x, 0≤x≤1 2, 1−x,1 2≤x≤1. The vibrational frequencies ωn=nπare the integral multiples of π, and so the natural modes of vibration are cosnπtsinnπx and sin nπtsinnπx forn= 1,2,... . Consequently, the general solution to the boundary value pr oblem is u(t,x) =∞/summationdisplay n=1/bracketleftbig bncosnπtsinnπx+dnsinnπtsinnπx/bracketrightbig , where bn= 2/integraldisplay1 0f(x) sinnπxdx=  4/integraldisplay1/2 0xsinnπxdx=4(−1)k (2k+1)2π2, n= 2k+1, 0, n= 2k, 12/11/12 767 c/ci∇cleco√y∇t2012 Peter J. Olver 0.2 0.4 0.6 0.8 1 -0.4-0.20.20.4 0.2 0.4 0.6 0.8 1 -0.4-0.20.20.4 0.2 0.4 0.6 0.8 1 -0.4-0.20.20.4 0.2 0.4 0.6 0.8 1 -0.4-0.20.20.4 0.2 0.4 0.6 0.8 1 -0.4-0.20.20.4 0.2 0.4 0.6 0.8 1 -0.4-0.20.20.4 Figure 14.6. Plucked String Solution of the Wave Equation. are the Fourier sine coeﬃcients of the initial displacement , whiledn= 0 are the Fourier sine coeﬃcients of the initial velocity. Therefore, the sol ution is the Fourier sine series u(t,x) = 4∞/summationdisplay k=0(−1)kcos(2k+1)πtsin(2k+1)πx (2k+1)2π2, (14.97) whose graph at times t= 0,.2,.4,.6,.8,1.is depicted in Figure 14.6. At time t= 1, the original displacement is reproduced exactly, but upside do wn. The subsequent dynamics proceeds as before, but in mirror image form. The original di splacement reappears at timet= 2, after which time the motion is periodically repeated. In terestingly, at times tk=.5,1.5,2.5,..., the displacement is identically zero: u(tk,x)≡0, although the velocity ut(tk,x)/ne}ationslash≡0. The solution appears to be piecewise aﬃne, i.e., its graph is a collection of straight lines. This fact will be veriﬁed in Exercise , where you are asked to construct an exact analytical formula for this solution. Unlike the he at equation, the wave equation doesnotsmooth out discontinuities and corners in the initial data. Whiletheseries form (14.92)ofthesolutionisnot entirely satisfying, wecanstilluse it todeduce important qualitativeproperties. First ofall, s ince each term isperiodicin twith period 2 ℓ/c, the entire solution is time periodic with that period: u(t+2ℓ/c,x) =u(t,x). In fact, after half the period, at time t=ℓ/c, the solution reduces to u/parenleftbiggℓ c,x/parenrightbigg =∞/summationdisplay n=1(−1)nbnsinnπx ℓ=−∞/summationdisplay n=1bnsinnπ(ℓ−x) ℓ=−u(0,ℓ−x) =−f(ℓ−x). In general, u/parenleftbigg t+ℓ c, x/parenrightbigg =−u(t,ℓ−x), u/parenleftbigg t+2ℓ c, x/parenrightbigg =u(t,x). (14.98) Therefore, theinitialwaveform isreproduced, ﬁrst asanup side downmirror imageofitself at timet=ℓ/c, and then in its original form at time t= 2ℓ/c. This has the important consequence that vibrations of (homogeneous) one-dimensi onal media are purely periodic 12/11/12 768 c/ci∇cleco√y∇t2012 Peter J. Olver phenomena! There is no quasi-periodicity because the funda mental frequencies are all integer multiples of each other. Remark: The preceding analysis has important musical consequence s. To the human ear, sonic vibrations that are integral multiples of a singl e frequency are harmonic, whereas those that admit quasi-periodic vibrations, with irration ally related frequencies, sound percussive. This is why most tonal instruments rely on vibra tions in one dimension, be it aviolinstring, acolumnofairinawindinstrument (ﬂute, cl arinet, trumpet orsaxophone), a xylophone bar or a triangle. On the other hand, most percuss ion instruments rely on the vibrations of two-dimensional media, e.g., drums and cy mbals, or three-dimensional solid bodies, e.g., blocks. As we shall see in Chapters 17 and 18, the frequency ratios of the latter are typically irrational. A bar with both ends left free, and so subject to the Neumann bo undary conditions ∂u ∂x(t,0) = 0 =∂u ∂x(t,ℓ), (14.99) willhave a slightlydiﬀerent behavior, owing to theinstabi lityofthe underlying equilibrium equations. The eigenfunctions of (14.89) with Neumann boun dary conditions v′(0) = 0 = v′(ℓ) are now vn(x) = cosnπx ℓwith ωn=nπc ℓ, n= 0,1,2,3,... . The resulting solution takes the form of a Fourier cosine ser ies u(t,x) =a0+c0t+∞/summationdisplay n=1/parenleftbigg ancosnπct ℓcosnπx ℓ+cnsinnπct ℓcosnπx ℓ/parenrightbigg .(14.100) In accordance with (14.90), the ﬁrst two terms come from the n ull eigenfunction v0(x) = 1 withω0= 0. The bar vibrates with the same fundamental frequencies ( 14.93) as in the ﬁxed end case, but there is now an additional unstable mode c0tthat is no longer periodic, but grows linearly in time. Substituting (14.100)intotheinitialconditions(14.85) ,weﬁnd theFouriercoeﬃcients are prescribed, as before, by the initial displacement and v elocity, an=2 ℓ/integraldisplayℓ 0f(x) cosnπx ℓdx, cn=2 nπc/integraldisplayℓ 0g(x) cosnπx ℓdx, n = 1,2,3,... . The order zero coeﬃcients†, a0=1 ℓ/integraldisplayℓ 0f(x)dx, c0=1 ℓ/integraldisplayℓ 0g(x)dx, †Note that, we have not included the usual1 2factor in the constant terms in the Fourier series (14.100). 12/11/12 769 c/ci∇cleco√y∇t2012 Peter J. Olver are equal to the average initial displacement and average in itial velocity of the bar. In particular, when c0= 0 there is no net initial velocity, and the unstable mode is n ot excited. In this case, the solution is time-periodic, oscil lating around the position given by the average initial displacement. On the other hand, if c0/ne}ationslash= 0, the bar will move oﬀ with constant average speed c0, all the while vibrating at the same fundamental frequencie s. Similar considerations apply to the periodic boundary valu e problem for the wave equation on a circular ring. The details are left as Exercise for the reader. Forcing and Resonance In Section 9.6, we learned that periodically forcing an unda mped mechanical struc- ture (or a resistanceless electrical circuit) at a frequenc y that is distinct from its natural vibrational frequencies leads, in general, to a quasi-peri odic response. The solution is a sum of the unforced vibrations superimposed with an additio nal vibrational mode at the forcing frequency. However, if forced at one of its natural f requencies, the system may go into a catastrophic resonance. The same type of quasi-periodic/resonant response is also o bserved in the partial diﬀerential equations governing periodic vibrations of co ntinuous media. To keep the analysis as simple as possible, we restrict our attention to the forced wave equation for a homogeneous bar ∂2u ∂t2=c2∂2u ∂x2+F(t,x), (14.101) subject to speciﬁed homogeneous boundary conditions. The e xternal forcing function F(t,x) may depend upon both time tand position x. We will be particularly interested in a periodically varying external force of the form F(t,x) = cos(ωt)h(x), (14.102) where the function h(x) is ﬁxed. As always, the solution to an inhomogeneous linear equation can be written as a combination, u(t,x) =u⋆(t,x)+z(t,x) (14 .103) of a particular solution u⋆(t,x) plus the general solution z(t,x) to the homogeneous equa- tion, namely ∂2z ∂t2=c2∂2z ∂x2. (14.104) The boundary and initial conditions will serve to uniquely p rescribe the solution u(t,x), but there is some ﬂexibility in its two constituents (14.103 ). For instance, we may ask that the particular solution u⋆satisfy the homogeneous boundary conditions along with zer o (homogeneous) initial conditions, and thus represents the pure response of the system to the forcing. The homogeneous solution z(t,x) will then reﬂect the eﬀect of the initial and boundary conditions unadulterated by the external forcing . The ﬁnal solution will equal the sum of the two individual responses. In the case of periodic forcing (14.102), we look for a partic ular solution u⋆(t,x) = cos(ωt)v⋆(x) (14 .105) 12/11/12 770 c/ci∇cleco√y∇t2012 Peter J. Olver that vibrates at the forcing frequency. Substituting the an satz (14.105) into the equa- tion (14.101), and canceling the common cosine factors, we d iscover that v⋆(x) must satisfy the boundary value problem prescribed by −c2v′′ ⋆−ω2v⋆=h(x), (14.106) supplemented by the relevant homogeneous boundary conditi ons — Dirichlet, Neumann, mixed, or periodic. At this juncture, there are two possibilities. If the unforc ed, homogeneous boundary value problem has only the trivial solution v≡0, then there is a solution to the forced boundary value problem for any form of the forcing function h(x). On the other hand, the homogeneous boundary value problem has a nontrivial soluti onv(x) whenω2=λis an eigenvalue, and so ωis a natural frequency of vibration to the homogeneous probl em; the solution v(x) is the corresponding eigenfunction appearing in the solut ion series (14.92). In this case, according to the Fredholm alternative, Theore m 5.55, the boundary value problem (14.106) has a solution if and only if the forcing fun ctionh(x) is orthogonal to the eigenfunction(s): /an}b∇acketle{th;v/an}b∇acket∇i}ht= 0. (14.107) See Example 11.3 and Exercise for details. If we force in a resonant manner — meaning that (14.107) is not satisﬁed — then the solution will be a res onantly growing vibration u⋆(t,x) =tsin(ωt)v⋆(x). In a real-world situation, such large resonant (or near reso nant) vibrations will either cause a catastrophic breakdown, e.g., the bar breaks or the string snaps, or will send the system into a diﬀerent, nonlinear regime that helps mitigate the re sonant eﬀects, but is no longer modeled by the simple linear wave equation. Example 14.7. As a speciﬁc example, consider the initial-boundary value p roblem modeling the forced vibrations of a uniform bar of unit lengt h and ﬁxed at both ends: utt=c2uxx+cos(ωt)h(x), u(t,0) = 0 = u(t,1), u(0,x) =f(x), ut(0,x) =g(x).(14.108) The particular solution will have the nonresonant form (14. 105) provided there exists a solution v⋆(x) to the boundary value problem c2v′′ ⋆+ω2v⋆=−h(x), v⋆(0) = 0 = v⋆(1). (14.109) The natural frequencies and associated eigenfunctions are ωn=ncπ, vn(x) = sinnπx, n = 1,2,3,... . The boundary value problem (14.109) will have a solution, an d hence the forcing is not resonant, provided either ω/ne}ationslash=ωnis not a natural frequency, or ω=ωn, but 0 =/an}b∇acketle{th;vn/an}b∇acket∇i}ht=/integraldisplay1 0h(x)sinnπxdx (14.110) 12/11/12 771 c/ci∇cleco√y∇t2012 Peter J. Olver is orthogonal to the associated eigenfunction. Otherwise, the forcing proﬁle will induce a resonant response. For example, under periodic forcing of frequency ωwith trigonometric proﬁle h(x)≡ sinkπx, the particular solution to (14.109) is v⋆(x) =sinkπx ω2−k2π2c2,so that u⋆(t,x) =cosωtsinkπx ω2−k2π2c2, (14.111) which is a valid solution as long as ω/ne}ationslash=ωk=kπc. Note that we may allow the forcing frequency ω=ωnto coincide with any other resonant forcing frequency, n/ne}ationslash=k, because the sine proﬁles are mutually orthogonal and so the nonreson ance condition (14.110) holds. On the other hand, if ω=ωk=kπc, then the particular solution u⋆(t,x) =tsinkπctsinkπx 2kπc, (14.112) is resonant, and grows linearly in time. To obtain the full solution to the initial-boundary value pr oblem, we write u=u⋆+z wherez(t,x) must satisfy ztt−c2zxx= 0, z (t,0) = 0 = z(t,1), along with the modiﬁed initial conditions z(0,x) =f(x)−sinkπx ω2−k2π2c2,∂u ∂x(0,x) =g(x), stemming from the fact that the particular solution (14.111 ) has non-trivial initial data. (In the resonant case (14.112), there is no extra term in the i nitial data.) Note that, the closerωis to the resonant frequency, the larger the modiﬁcation of t he initial data, and hence the larger the response of the system to the periodic fo rcing. As before, the solution z(t,x) to the homogeneous equation can be written as a Fourier sine series (14.92). The ﬁnal formulae are left to the reader to write out. 14.5. d’Alembert’s Solution. The one-dimensional wave equation is distinguished by admi tting an explicit solu- tion formula, originally discovered by the eighteenth cent ury French mathematician Jean d’Alembert, that entirely avoids the complicated Fourier s eries form. D’Alembert’s solu- tion provides additional valuable insight into the behavio r of the solutions. Unfortunately, unlike series methods that have a very broad range of applica bility, d’Alembert’s method only succeeds in this one very special situation: the homoge neous wave equation in a single space variable. The starting point is to write the wave equation (14.86) in th e suggestive form /squareu= (∂2 t−c2∂2 x)u=utt−c2uxx= 0. (14.113) Here/square=∂2 t−c2∂2 xis a common mathematical notation for the wave operator , while ∂t,∂xare convenient shorthands for the partial derivative opera tors with respect to tand 12/11/12 772 c/ci∇cleco√y∇t2012 Peter J. Olver x. We note that /squareis a linear, second order partial diﬀerential operator. In a nalogy with the elementary polynomial factorization t2−c2x2= (t−cx)(t+cx), we shall factor the wave operator into a product of two ﬁrst or der partial diﬀerential operators: /square=∂2 t−c2∂2 x= (∂t−c∂x)(∂t+c∂x). (14.114) Now, if the second factor annihilates the function u(t,x), meaning (∂t+c∂x)u=ut+cux= 0, (14.115) thenuis automatically a solution to the wave equation, since /squareu= (∂t−c∂x)(∂t+c∂x)u= (∂t−c∂x)0 = 0. Inotherwords, everysolutiontothesimplerﬁrstorderpart ialdiﬀerentialequation(14.115) is a solution to the wave equation (14.86). (The converse is, of course, not true.) It is relatively easy to solve linear†ﬁrst order partial diﬀerential equations. Proposition 14.8. Every solution u(t,x)to the partial diﬀerential equation ∂u ∂t+c∂u ∂x= 0 (14 .116) has the form u(t,x) =p(x−ct), (14.117) wherep(ξ)is an arbitrary function of the characteristic variable ξ=x−ct. Proof: We adopt a linear change of variables to rewrite the solutio n u(t,x) =p(t,x−ct) =p(t,ξ) as a function of the characteristic variable ξand the time t. Applying the chain rule, we express the derivatives of uin terms of the derivatives of pas follows: ∂u ∂t=∂p ∂t−c∂p ∂ξ,∂u ∂x=∂p ∂ξ, and hence ∂u ∂t+c∂u ∂x=∂p ∂t−c∂p ∂ξ+c∂p ∂ξ=∂p ∂t. Therefore, uis a solution to (14.116) if and only if p(t,ξ) is a solution to the very simple partial diﬀerential equation ∂p ∂t= 0. †See Chapter 22 for more details, including extensions to ﬁrst order n onlinear partial diﬀer- ential equations. 12/11/12 773 c/ci∇cleco√y∇t2012 Peter J. Olver Figure 14.7. Traveling Wave. This clearly‡implies that p=p(ξ) does not depend on the variable t, and hence u=p(ξ) =p(x−ct) is of the desired form. Q.E.D. Therefore, anyfunction of the characteristic variable, e.g., ξ2+1, or cos ξ, oreξ, will produce a corresponding solution, ( x−ct)2+1, or cos( x−ct), orex−ct, to the ﬁrst order partial diﬀerential equation (14.116), and hence a solutio n to the wave equation (14.86). The functions of the form u(t,x) =p(x−ct) are known as traveling waves . Att= 0 the wave has the initial proﬁle u(0,x) =p(x). Astprogresses, the wave moves to the right, withconstant speed c >0andunchanged inform; seeFigure14.7. Forthisreason, (14 .116) is sometimes referred to as the one-way orunidirectional wave equation . Proposition 14.8 tells us that every traveling wave with wave speed cis a solution to the full wave equation (14.86). But keep in mind that such solutions do not necessar ily respect the boundary conditions, which, when present, will aﬀect their ultimate behavior. Now, since cis constant, the factorization (14.114) can be written equa lly well in the reverse order: /square=∂2 t−c2∂2 x= (∂t+c∂x)(∂t−c∂x). (14.118) The same argument tells us that any solution to the alternati ve ﬁrst order partial diﬀer- ential equation ∂u ∂t−c∂u ∂x= 0, (14.119) also provides a solution to the wave equation. This is also a u nidirectional wave equation, but with the opposite wave speed −c. Applying Proposition 14.8, now with creplaced by −c, we conclude that the general solution to (14.119) has the fo rm u(t,x) =q(x+ct) (14 .120) whereq(η) is an arbitrary diﬀerentiable function of the alternate ch aracteristic variable η=x+ct. The solutions (14.120) represent traveling waves moving t o theleftwith constant speed c >0 and unchanged in form. The wave equation (14.113) is bidirectional and admits both left and right traveling wave solutions. Linearity of the wave equation implies that the sum of solutions is again a solution. In this way, we can produce solutions which are sup erpositions of left and right traveling waves. The remarkable fact is that everysolution to the wave equation can be so represented. ‡More rigorously, one should also assume that, at each time t, the domain of deﬁnition of p(ξ) is a connected interval. A similar technical restriction should be imposed upon the solutions in the statement of Proposition 14.8. See Exercise for a detailed example. 12/11/12 774 c/ci∇cleco√y∇t2012 Peter J. Olver Theorem 14.9. Every solution to the wave equation (14.86)can be written as a combination u(t,x) =p(ξ)+q(η) =p(x−ct)+q(x+ct) (14 .121) of right and left traveling waves, each depending on its char acteristic variable ξ=x−ct, η =x+ct. (14.122) Proof: The key is to use a linear changes of variables to rewrite the wave equation entirely in terms of the characteristic variables. We set u(t,x) =w(x−ct,x+ct) =w(ξ,η),whereby w(ξ,η) =u/parenleftbiggη−ξ 2c,ξ+η 2/parenrightbigg . Then, using the chain rule to compute the partial derivative s, ∂u ∂t=c/parenleftbigg∂w ∂ξ−∂w ∂η/parenrightbigg ,∂u ∂x=∂w ∂ξ+∂w ∂η. and hence ∂2u ∂t2=c2/parenleftbigg∂2w ∂ξ2−2∂2w ∂ξ∂η+∂2w ∂η2/parenrightbigg ,∂2u ∂x2=∂2w ∂ξ2+2∂2w ∂ξ∂η+∂2w ∂η2. Therefore /squareu=∂2u ∂t2−c2∂2u ∂x2=−4c2∂2w ∂ξ∂η. We conclude that u(t,x) solves the wave equation /squareu= 0 if and only if w(ξ,η) solves the second order partial diﬀerential equation ∂2w ∂ξ∂η= 0,which we write in the form∂ ∂ξ/parenleftbigg∂w ∂η/parenrightbigg = 0. This partial diﬀerential equation can be integrated once wi th respect to ξ, resulting in ∂w ∂η=r(η), whereris an arbitrary function of the characteristic variable η. Integrating both sides of the latter partial diﬀerential equation with respect to η, we ﬁnd w(ξ,η) =p(ξ)+q(η),where q′(η) =r(η), whilep(ξ) represents the integration “constant”. Replacing the cha racteristic variables by their formulae in terms of tandxcompletes the proof. Q.E.D. Remark: As noted above, we have been a little cavalier with our speci ﬁcation of the domainofdeﬁnitionofthefunctionsandthediﬀerentiabili tyassumptionsrequired. Sorting out the precise technical details is left to the meticulous r eader. 12/11/12 775 c/ci∇cleco√y∇t2012 Peter J. Olver Remark: As we know, the general solution to a second order ordinary diﬀerential equation depends on two arbitrary constants. Here we observ e that the general solution to a second order partialdiﬀerential equation depends on two arbitrary functions — i n this casep(ξ) andq(η). This observation serves as a useful rule of thumb, but shou ld not be interpreted too literally. Let us now see how this new form of solution to the wave equatio n can be used to eﬀectively solve initial value problems. The simplest case is that of a bar or string of inﬁnite length, in which case we have a pure initial value pro blem ∂2u ∂t2=c2∂2u ∂x2, u(0,x) =f(x),∂u ∂t(0,x) =g(x),for−∞< x <∞. Substituting the solution formula (14.121) into the initia l conditions, we ﬁnd u(0,x) =p(x)+q(x) =f(x),∂u ∂t(0,x) =−cp′(x)+cq′(x) =g(x). To solve this pair of linear equations for pandq, we diﬀerentiate the ﬁrst equation: p′(x)+q′(x) =f′(x). Subtracting the second equation divided by c, we ﬁnd 2p′(x) =f′(x)−1 cg(x). Therefore, p(x) =1 2f(x)−1 2c/integraldisplayx 0g(z)dz+a, whereais an integration constant. The ﬁrst equation then yields q(x) =f(x)−p(x) =1 2f(x)+1 2c/integraldisplayx 0g(z)dz−a. Substituting these two expressions back into (14.121), we ﬁ nd u(t,x)=p(ξ)+q(η) =f(ξ)+f(η) 2+1 2c/bracketleftigg −/integraldisplayξ 0+/integraldisplayη 0/bracketrightigg g(z)dz =f(ξ)+f(η) 2+1 2c/integraldisplayη ξg(z)dz, whereξ,ηare the characteristic variables (14.122). In this fashion , we have derived d’Alembert’s solution to the wave equation on the entire line −∞< x <∞. Theorem 14.10. The solution to the initial value problem ∂2u ∂t2=c2∂2u ∂x2, u(0,x) =f(x),∂u ∂t(0,x) =g(x),−∞< x <∞.(14.123) 12/11/12 776 c/ci∇cleco√y∇t2012 Peter J. Olver Figure 14.8. Interaction of Waves. is given by u(t,x) =f(x−ct)+f(x+ct) 2+1 2c/integraldisplayx+ct x−ctg(z)dz. (14.124) . Let us investigate the implications of d’Alembert’s soluti on formula (14.124). First, suppose there is no initial velocity, so g(x)≡0, and the motion is purely the result of the initial displacement u(0,x) =f(x). In this case, the solution (14.124) reduces to u(t,x) =1 2f(x−ct)+1 2f(x+ct). The eﬀect is that the initial displacement f(x) splits into two waves, one traveling to the right and one traveling to the left, each of constant speed c, and each of exactly the same shape as the initial displacement f(x) but only half as tall. For example, if the initial displacement is a localized pulse centered at the origin, sa y u(0,x) =e−x2,∂u ∂t(0,x) = 0, then the solution u(t,x) =1 2e−(x−ct)2+1 2e−(x+ct)2 consists of two half size pulses running away from the origin in opposite directions with equal speed c. If we take two separated pulses, say u(0,x) =e−x2+2e−(x−1)2,∂u ∂t(0,x) = 0, centered at x= 0 and x= 1, then the solution u(t,x) =1 2e−(x−ct)2+e−(x−1−ct)2+1 2e−(x+ct)2+e−(x−1+ct)2 will consist of four pulses, two moving to the right and two to the left, all with the same speed, as pictured in Figure 14.8. Animportantobservationisthatwhenaright-movingpulsec ollideswithaleft-moving pulse, they emerge from the collision unchanged — a conseque nce of the inherent linearity of the wave equation. The ﬁrst picture shows the initial disp lacement. In the second and third pictures, the two localized bumps have each split into two copies moving in opposite directions. In the fourth and ﬁfth, the larger right moving b ump is in the process of interacting with the smaller left moving bump. Finally, in t he last picture the interaction is complete, and the two left moving bumps and two right movin g bumps travel in tandem with no further collisions. 12/11/12 777 c/ci∇cleco√y∇t2012 Peter J. Olver Figure 14.9. Characteristic Lines for the Wave Equation. Remark: If the initial displacement has bounded support, and so f(x) = 0 ifx < a orx > bfor some a < b, then after a ﬁnite time the right and left-moving waves will be completely separated, and the observer will see two exact ha lf size replicas running away, with speed c, in opposite directions. If the displacement is not localiz ed, then the left and right traveling waves will never fully disengage, and one mi ght be hard pressed (just as in our earlier discussion of quasi-periodic phenomena) in r ecognizing that a complicated solution pattern is, in reality, just the superposition of t wo simple traveling waves. For example, using a trigonometric identity, any separable sol ution, e.g., cosctsinx=1 2cos(x−ct)+1 2cos(x+ct), can be rewritten in d’Alembert form. The interpretation of t he two solution formulas is quitediﬀerent. Mostobserverswillseeastandingsinusoid alwave, vibratingwithfrequency c, as represented on the left side of the equation. However, th e right hand side says that this is the same as the diﬀerence of a right and left traveling cosine wave. The interactions of their peaks and troughs reproduces the standing wave. Thu s, the same solution can be interpreted in seemingly incompatible ways! In general, the lines of slope ±cin the (t,x)–plane, where the characteristic variables are constant, ξ=x−ct=a, η =x+ct=b, (14.125) are known as the characteristics of the wave equation. The two characteristics going through a point on the xaxis, where the initial data is prescribed, are illustrated in Figure 14.9. The reader should note that, in this ﬁgure, the taxis is horizontal, while the xaxis is vertical. In general, signals propagate along characteristics. More speciﬁcally, if we start out with an initial displacement concentrated very close to a po intx=a,t= 0, then the solution will be concentrated along the two characteristic lines emanating from the point, namelyx−ct=aandx+ct=a. Inthelimit,aunitimpulseordeltafunctiondisplacement atx=a, corresponding to the initial condition u(0,x) =δ(x−a),∂u ∂t(0,x) = 0, (14.126) will result in a solution u(t,x) =1 2δ(x−ct−a)+1 2δ(x+ct−a) (14 .127) 12/11/12 778 c/ci∇cleco√y∇t2012 Peter J. Olver (t,x)(0,x+ct) (0,x−ct) Figure 14.10. Domain of Dependence. consisting of two half-strength delta spikes traveling awa y from the starting position con- centrated on the two characteristic lines. Let us return to the general initial value problem (14.123). Suppose now that there is no initial displacement, u(0,x) =f(x)≡0, but rather a concentrated initial velocity, say a delta function ∂u ∂t(0,x) =δa(x) =δ(x−a). Physically, this would correspond to striking the string wi th a highly concentrated blow at the point x=a. The d’Alembert solution (14.124) to this “hammer blow” pro blem is u(t,x) =1 2c/integraldisplayx+ct x−ctδa(z)dz=  1 2c, x−ct < a < x +ct, 0,otherwise ,(14.128) and consists of a constant displacement, of magnitude 1 /(2c), between the two character- istic lines x−ct=a=x+ctbased at x=a,t= 0 — the shaded region of Figure 14.9. This region is known as the domain of inﬂuence of the point ( a,0) since, in general, the value of the initial data at that point will only aﬀect the sol ution values in the region. Vice versa, the value of the solution u(t,x) a point with t >0 depends only on that part of the initial data that lies within the domain of dependence of the point ( t,x), which is the triangle whose sides have slope ±cextending back to the x-axis; see Figure 14.10. The particular solution, which is plotted in Figure 14.11, h as two jump discontinuities between the undisturbed state and the displaced state, each propagating along its char- acteristic line with speed c, but in opposite directions. Note that, unlike a concentrat ed initial displacement, where the signal remains concentrat ed and each point along the bar is temporarily displaced, eventually returning to its undist urbed state, a concentrated initial velocity has a lasting eﬀect, and the bar remains permanentl y deformed by an amount 1/(2c). 12/11/12 779 c/ci∇cleco√y∇t2012 Peter J. Olver 0.2 0.4 0.6 0.8 11.2 1.4 -0.20.20.40.60.811.2 0.2 0.4 0.6 0.8 11.2 1.4 -0.20.20.40.60.811.2 0.2 0.4 0.6 0.8 11.2 1.4 -0.20.20.40.60.811.2 Figure 14.11. Concentrated Initial Velocity for Wave Equation. Figure 14.12. Odd Periodic Extension of a Concentrated Pulse. Solutions on Bounded Intervals So far, we have been using d’Alembert formula to solve the wav e equation on an inﬁnite interval. The formula can still be used on bounded in tervals, but in a suitably modiﬁed format so as to respect the boundary conditions. The easiest to deal with is the periodic problem on 0 ≤x≤ℓ, with boundary conditions u(t,0) =u(t,ℓ), ux(t,0) =ux(t,ℓ). (14.129) If we extend the initial displacement f(x) and velocity g(x) to be periodic functions of periodℓ, sof(x+ℓ) =f(x)andg(x+ℓ) =g(x)for allx∈R, then the resulting d’Alembert solution (14.124) will also be periodic in x, sou(t,x+ℓ) =u(t,x). In particular, it satisﬁes the boundary conditions (14.129) and so coincides with the d esired solution. Details can be found in Exercises –. Next, suppose we have ﬁxed (Dirichlet) boundary conditions u(t,0) = 0, u (t,ℓ) = 0. (14.130) The resulting solution can be written as a Fourier sine serie s (14.92), and hence is both odd and 2 ℓperiodic in x. Therefore, to write the solution in d’Alembert form, we ext end 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 o dd and periodic. As a result, it satisﬁes the boundary conditions (14.130) for all t. Keep in mind that, while the solution u(t,x)isdeﬁnedforall x,theonlyphysicallyrelevantvaluesoccurontheinterval0 ≤x≤ℓ. Nevertheless, the eﬀects of displacements in the unphysica l regime will eventually be “felt” as the propagating waves pass through the physical interval . For example, consider an initial displacement which is conc entrated near x=afor some 0< a < ℓ . Its odd, periodic extension consists of two sets of replica s: those of the same form occurring at positions a±2ℓ,a±4ℓ, ..., and their mirror images at the intermediate positions −a,−a±2ℓ,−a±4ℓ, ...; Figure 14.12 shows a representative 12/11/12 780 c/ci∇cleco√y∇t2012 Peter J. Olver 0.2 0.4 0.6 0.8 1 -1-0.50.51 0.2 0.4 0.6 0.8 1 -1-0.50.51 0.2 0.4 0.6 0.8 1 -1-0.50.51 0.2 0.4 0.6 0.8 1 -1-0.50.51 0.2 0.4 0.6 0.8 1 -1-0.50.51 0.2 0.4 0.6 0.8 1 -1-0.50.51 Figure 14.13. Solution to Wave Equation with Fixed Ends. example. The resulting solution begins with each of the puls es, both positive and negative, splitting into two half-size replicas that propagate with s peedcin opposite directions. As the individual pulses meet, they interact as they pass unalt ered through each other. The process repeats periodically, with an inﬁnite row of half-s ize pulses moving to the right kaleidoscopically interacting with an inﬁnite row moving t o the left. However, only the part of this solution that lies on 0 ≤x≤ℓis actually realized on the physical bar. The net eﬀect is as if we were forced to view t he complete solution as it passes by a window of length ℓthat blocks out all other regions of the real axis. What the viewer eﬀectively sees assumes a somewhat diﬀerent interpr etation. To wit, the original pulse at position 0 < a < ℓ splits up into two half size replicas that move oﬀ in opposite directions. As each half-size pulse reaches an end of the bar , it meets a mirror image pulse that has been propagating in the opposite direction from the non-physical regime. The pulse appears to be reﬂected at the end of the interval, and ch anges into an upside down mirror image moving in the opposite direction. The original positive pulse has moved oﬀ the end of the bar just as its mirror image has moved into the physical regime. (A common physical illustration is a pulse propagating down a j ump rope that is held ﬁxed at its end; the reﬂected pulse returns upside down.) A similar r eﬂection 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=ℓ/cthey recombine at the point ℓ−ato instantaneously form a full-sized, but upside-down mirr or image of the original disturbance — in accordance with (14.9 8). The recombined pulse in turn splits apart into two upside down half-size pulses th at, when each collides with the end, reﬂects and returns to its original upright form. At timet= 2ℓ/c, the pulses recombine to exactly reproduce the original displacement. The process then repeats, and the solution is periodic in time with period 2 ℓ/c. In Figure 14.13, the ﬁrst picture displays the initial displ acement. In the second, it has splits into left and right moving, half-size clones. In t he third picture, the left moving bump is in the process of colliding with the left end of the bar . In the fourth picture, it has emerged from the collision, and is now upside down, reﬂ ected, and moving to the 12/11/12 781 c/ci∇cleco√y∇t2012 Peter J. Olver right. Meanwhile, the right moving pulse is starting to coll ide with the right end. In the ﬁfth picture, both pulses have completed their collisions a nd are now moving back towards each other, where, in the last picture, they recombine into a n upside-down mirror image of the original pulse. The process then repeats itself, in mi rror image, ﬁnally recombining to the original pulse, at which point the entire process star ts over. The Neumann (free) boundary value problem ∂u ∂x(t,0) = 0,∂u ∂x(t,ℓ) = 0, (14.131) is handled similarly. Since the solution has the form of a Fou rier 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 is also even and 2 ℓperiodic in x, and hence satisﬁes the boundary conditions. In this case, when a pulse hits one of th e ends, its reﬂection remains upright, but becomes a mirror image of the original; a famili ar physical illustration is a water wave that reﬂects oﬀ a solid wall. Further details are l eft to the reader in Exercise In summary, we have now learned two diﬀerent ways to solve the one-dimensional wave equation. The ﬁrst, based on Fourier analysis, emphasi zes the vibrational or wave character of the solutions, while the second, based on the d’ Alembert formula, emphasizes their particle aspects, where individual wave packets coll ide with each other, or reﬂect at the boundary, all the while maintaining their overall for m. Some solutions look like vibrating waves, while others seem much more like interacti ng particles. But, like the blind men describing the elephant, these are merely two face ts of the samesolution. The Fourier series formula shows how every particle-like solut ion can be decomposed into its constituent vibrational modes, while the d’Alembert formu la demonstrates how vibrating waves combine as moving particles. The coexistence of particle and wave features is reminiscen t of the long running histor- ical debate over the nature of light. In the beginning, Newto n and his disciples advocated a particle basis, in the form of photons. However, until the b eginning of the twentieth cen- tury, most physicists advocated a wave or vibrational viewp oint. Einstein’s explanation of the photoelectric eﬀect in 1905 served to resurrect the part icle interpretation. Only with the establishment of quantum mechanics was the debate resol ved — light, and, indeed, all subatomic particles are both, manifesting particle and wave features, depending upon the experiment and the physical situation. But the theoreti cal evidence for wave-particle duality already existed in the competing solution formulae of the classical wave equation! 14.6. Numerical Methods. As you know, most diﬀerential equations are much too complic ated to be solved an- alytically. Thus, to obtain quantitative results, one is fo rced to construct a suﬃciently accurate numerical approximation to the solution. Even in c ases, such as the heat and wave equations, where explicit solution formulas (either c losed form or inﬁnite series) exist, numerical methods still can be proﬁtably employed. Moreove r, justiﬁcation of the numeri- cal algorithm is facilitated by the ability compare it with a n exact solution. Moreover, the 12/11/12 782 c/ci∇cleco√y∇t2012 Peter J. Olver lessons learned in the design of numerical algorithms for “s olved” problems prove to be of immense value when one is confronted with more complicated p roblems for which solution formulas no longer exist. In this ﬁnal section, we present some of the most basic numeri cal solution techniques for the heat and wave equations. We will only introduce the mo st basic algorithms, leaving more sophisticated variations and extensions to a more thor ough treatment, which can be found in numerical analysis texts, e.g., [ 27,35,109]. Finite Diﬀerences Numerical solutionmethods for diﬀerential equationscan b epartitionedinto two prin- cipalclasses. (Inthisoversimpliﬁedpresentation, weare ignoringmorespecializedmethods of less general applicability.) The ﬁrst category, already introduced in Section 11.6, are the ﬁnite element methods . Finite elements are designed for the diﬀerential equation s describ- ing equilibrium conﬁgurations, since they rely on minimizi ng a functional. A competitive alternativeis to directly approximate the derivatives app earing in the diﬀerential equation, which requires knowing appropriate numerical diﬀerentiat ion formulae. In general, to approximate the derivative of a function at a p oint, say f′(x) orf′′(x), one constructs a suitable combination of sampled function v alues at nearby points. The underlying formalism used to construct these approximatio n formulae is known as the calculus of ﬁnite diﬀerences . Its development has a long and inﬂuential history, dating back to Newton. The resulting ﬁnite diﬀerence numerical methods for solving diﬀerential equations have extremely broad applicability, and can, wit h proper care, be adapted to most problems that arise in mathematics and its many applica tions. The simplest ﬁnite diﬀerence approximation is the ordinary diﬀerence quotient u(x+h)−u(x) h≈u′(x), (14.132) used to approximate the ﬁrst derivative of the function u(x). Indeed, if uis diﬀerentiable atx, thenu′(x) is, by deﬁnition, the limit, as h→0 of the ﬁnite diﬀerence quotients. Geometrically, the diﬀerence quotient equals the slope of t he secant line through the two points/parenleftbig x,u(x)/parenrightbig and/parenleftbig x+h,u(x+h)/parenrightbig on the graph of the function. For small h, this should be a reasonably good approximation to the slope of the tangent line, u′(x), as illustrated in the ﬁrst picture in Figure 14.14. How close an approximation is the diﬀerence quotient? To ans wer this question, we assume that u(x) is at least twice continuously diﬀerentiable, and examine the ﬁrst order Taylor expansion u(x+h) =u(x)+u′(x)h+1 2u′′(ξ)h2. (14.133) We have used the Cauchy form (C.8) for the remainder term, in w hichξrepresents some point lying between xandx+h. Theerroror diﬀerence between the ﬁnite diﬀerence formula and the derivative being approximated is given by u(x+h)−u(x) h−u′(x) =1 2u′′(ξ)h. (14.134) 12/11/12 783 c/ci∇cleco√y∇t2012 Peter J. Olver One-Sided Diﬀerence Central Diﬀerence Figure 14.14. Finite Diﬀerence Approximations. Since the error is proportional to h, we say that the ﬁnite diﬀerence quotient (14.134) is a ﬁrst order approximation. When the precise formula for the error is not so important, we will write u′(x) =u(x+h)−u(x) h+O(h). (14.135) The “big Oh” notation O( h) refers to a term that is proportional to h, or, more rigorously, bounded by a constant multiple of hash→0. Example 14.11. Letu(x) = sinx. Let us try to approximate u′(1) = cos1 = 0.5403023...by computing ﬁnite diﬀerence quotients cos1≈sin(1+h)−sin1 h. The result for diﬀerent values of his listed in the following table. h 1 .1 .01 .001 .0001 approximation 0.067826 0 .497364 0 .536086 0 .539881 0 .540260 error −0.472476 −0.042939 −0.004216 −0.000421 −0.000042 We observe that reducing the step size by a factor of1 10reduces the size of the error by approximatelythesamefactor. Thus, toobtain10decimaldi gitsofaccuracy, weanticipate needing a step size of about h= 10−11. The fact that the error is more of less proportional to the step size conﬁrms that we are dealing with a ﬁrst order n umerical approximation. To approximate higher order derivatives, we need to evaluat e the function at more than two points. In general, an approximation to the nthorder derivative u(n)(x) requires at leastn+1 distinct sample points. For simplicity, we shall only use equallyspaced points, leaving the general case to the exercises. For example, let us try to approximate u′′(x) by sampling uat the particular points x, x+handx−h. Which combination of the function values u(x−h),u(x),u(x+h) should be used? The answer to such a question can be found by consider ation of the relevant 12/11/12 784 c/ci∇cleco√y∇t2012 Peter J. Olver Taylor expansions u(x+h) =u(x)+u′(x)h+u′′(x)h2 2+u′′′(x)h3 6+O(h4), u(x−h) =u(x)−u′(x)h+u′′(x)h2 2−u′′′(x)h3 6+O(h4),(14.136) where the error terms are proportional to h4. Adding the two formulae together gives u(x+h)+u(x−h) = 2u(x)+u′′(x)h2+O(h4). Rearranging terms, we conclude that u′′(x) =u(x+h)−2u(x)+u(x−h) h2+O(h2), (14.137) The result is known as the centered ﬁnite diﬀerence approximation to the second derivative of a function. Since the error is proportional to h2, this is a second order approximation. Example 14.12. Letu(x) =ex2, withu′′(x) = (4x2+2)ex2. Let us approximate u′′(1) = 6e= 16.30969097 ...by using the ﬁnite diﬀerence quotient (14.137): 6e≈e(1+h)2−2e+e(1−h)2 h2. The results are listed in the following table. h 1 .1 .01 .001 .0001 approximation 50.16158638 16 .48289823 16 .31141265 16 .30970819 16 .30969115 error 33.85189541 0 .17320726 0 .00172168 0 .00001722 0 .00000018 Each reduction in step size by a factor of1 10reduces the size of the error by a factor of 1 100and results in a gain of two new decimal digits of accuracy, co nﬁrming that the ﬁnite diﬀerence approximation is of second order. However, thispredictionisnotcompletelyborneoutinprac tice. Ifwetake†h=.00001 thentheformulaproducestheapproximation16 .3097002570,withanerrorof0 .0000092863 — which is lessaccurate that the approximation with h=.0001. The problem is that round-oﬀ errors have now begun to aﬀect the computation, and underscores the diﬃculty with numerical diﬀerentiation. Finite diﬀerence formulae involve dividing very small quan- tities, which can induce large numerical errors due to round -oﬀ. As a result, while they typically produce reasonably good approximations to the de rivatives for moderately small step sizes, to achieve high accuracy, one must switch to a hig her precision. In fact, a similar comment applied to the previous Example 14.11, and o ur expectations about the error were not, in fact, fully justiﬁed as you may have discov ered if you tried an extremely small step size. †This next computation depends upon the computer’s precision; here we used single precision inMatlab . 12/11/12 785 c/ci∇cleco√y∇t2012 Peter J. Olver Another way to improve the order of accuracy of ﬁnite diﬀeren ce approximations is to employ more sample points. For instance, if the ﬁrst order ap proximation (14.135) to the ﬁrst derivative based on the two points xandx+his not suﬃciently accurate, one can try combining the function values at three points x,x+handx−h. To ﬁnd the appropriate combination of u(x−h),u(x),u(x+h), we return to the Taylor expansions (14.136). To solve for u′(x), we subtract‡the two formulae, and so u(x+h)−u(x−h) = 2u′(x)h+u′′′(x)h3 3+O(h4). Rearranging the terms, we are led to the well-known centered diﬀerence formula u′(x) =u(x+h)−u(x−h) 2h+O(h2), (14.138) which is a second order approximation to the ﬁrst derivative . Geometrically, the cen- tered diﬀerence quotient represents the slope of the secant line through the two points/parenleftbig x−h,u(x−h)/parenrightbig and/parenleftbig x+h,u(x+h)/parenrightbig on the graph of ucentered symmetrically about the point x. Figure 14.14 illustrates the two approximations; the adva ntages in accuracy in thecentered diﬀerence versionaregraphicallyevident. Hi gherorder approximationscanbe found by evaluating the function at yet more sample points, i ncluding, say, x+2h,x−2h, etc. Example 14.13. Return to the function u(x) = sinxconsidered in Example 14.11. The centered diﬀerence approximation to its derivative u′(1) = cos1 = 0 .5403023 ...is cos1≈sin(1+h)−sin(1−h) 2h. The results are tabulated as follows: h .1 .01 .001 .0001 approximation 0.53940225217 0 .54029330087 0 .54030221582 0 .54030230497 error −0.00090005370 −0.00000900499 −0.00000009005 −0.00000000090 As advertised, the results are much more accurate than the on e-sided ﬁnite diﬀerence approximation used in Example 14.11 at the same step size. Si nce it is a second order approximation, each reduction in the step size by a factor of1 10results in two more decimal places of accuracy. Many additionalﬁnite diﬀerence approximationscanbe cons tructed by similar manip- ulations of Taylor expansions, but these few very basic ones will suﬃce for our subsequent purposes. In the following subsection, we apply the ﬁnite di ﬀerence formulae to develop numerical solution schemes for the heat and wave equations. ‡The terms O( h4) donotcancel, since they represent potentially diﬀerent multiples of h4. 12/11/12 786 c/ci∇cleco√y∇t2012 Peter J. Olver Numerical Algorithms for the Heat Equation Consider the heat equation ∂u ∂t=γ∂2u ∂x2,0< x < ℓ, t ≥0, (14.139) representing a bar of length ℓand thermal diﬀusivity γ >0, which is assumed to be constant. To be concrete, we impose time-dependent Dirichl et boundary conditions u(t,0) =α(t), u (t,ℓ) =β(t), t ≥0, (14.140) specifying the temperature at the ends of the bar, along with the initial conditions u(0,x) =f(x), 0≤x≤ℓ, (14.141) specifying the bar’s initial temperature distribution. In order to eﬀect a numerical approx- imation to the solution to this initial-boundary value prob lem, we begin by introducing a rectangular mesh consisting of points ( ti,xj) with 0 =x0< x1<···< xn=ℓand 0 = t0< t1< t2<···. For simplicity, we maintain a uniform mesh spacing in both di rections, with h=xj+1−xj=ℓ n, k =ti+1−ti, representing, respectively, the spatial mesh size and the t ime step size. It will be essential that we do nota priori require the two to be the same. We shall use the notati on ui,j≈u(ti,xj) where ti=ik, xj=jh, (14.142) to denote the numerical approximation to the solution value at the indicated mesh point. As a ﬁrst attempt at designing a numerical method, we shall us e the simplest ﬁnite diﬀerence approximations to the derivatives. The second or der space derivative is approx- imated by (14.137), and hence ∂2u ∂x2(ti,xj)≈u(ti,xj+1)−2u(ti,xj)+u(ti,xj−1) h2+O(h2) ≈ui,j+1−2ui,j+ui,j−1 h2+O(h2),(14.143) where the error in the approximation is proportional to h2. Similarly, the one-sided ﬁnite diﬀerence approximation (14.135) is used for the time deriv ative, and so ∂u ∂t(ti,xj)≈u(ti+1,xj)−u(ti,xj) k+O(k)≈ui+1,j−ui,j k+O(k),(14.144) where the error is proportion to k. In practice, one should try to ensure that the approxi- mations have similar orders of accuracy, which leads us to ch oose k≈h2. Assuming h <1, this requirement has the important consequence that the t ime steps must bemuchsmaller than the space mesh size. 12/11/12 787 c/ci∇cleco√y∇t2012 Peter J. Olver Remark: At this stage, the reader might be tempted to replace (14.14 4) by the sec- ond order central diﬀerence approximation (14.138). Howev er, this produces signiﬁcant complications, and the resulting numerical scheme is not pr actical. Replacing the derivatives in the heat equation (14.145) by t heir ﬁnite diﬀerence ap- proximations (14.143), (14.144), and rearranging terms, w e end up with the linear system ui+1,j=µui,j+1+(1−2µ)ui,j+µui,j−1,i= 0,1,2,... , j= 1,...,n−1,(14.145) in which µ=γk h2. (14.146) The resulting numerical scheme takes the form of an iterativ elinear system for the solution valuesui,j≈u(ti,xj),j= 1,...,n−1, at each time step ti. The initial condition (14.141) means that we should initial ize our numerical data by sampling the initial temperature at the mesh points: u0,j=fj=f(xj), j = 1,...,n−1. (14.147) Similarly, the boundary conditions (14.140) require that ui,0=αi=α(ti), ui,n=βi=β(ti), i = 0,1,2,... . (14.148) For consistency, we should assume that the initial and bound ary conditions agree at the corners of the domain: f0=f(0) =u(0,0) =α(0) =α0, fn=f(ℓ) =u(0,ℓ) =β(0) =β0. The three equations (14.145–148) completely prescribe the numerical approximation algo- rithm for solving the initial-boundary value problem (14.1 39–141). Let us rewrite the scheme in a more transparent matrix form. F irst, let u(i)=/parenleftbig ui,1,ui,2,...,ui,n−1/parenrightbigT≈/parenleftbig u(ti,x1),u(ti,x2),...,u(ti,xn−1)/parenrightbigT(14.149) be the vector whose entries are the numerical approximation s to the solution values at timetiat theinterior nodes. We omit the boundary nodes x0= 0, xn=ℓ, since those values are ﬁxed by the boundary conditions (14.140). Then (1 4.145) assumes the compact vectorial form u(i+1)=Au(i)+b(i), (14.150) where A= 1−2µ µ µ1−2µ µ µ1−2µ µ µ...... ......µ µ1−2µ ,b(i)= µαi 0 0 ... 0 µβi .(14.151) 12/11/12 788 c/ci∇cleco√y∇t2012 Peter J. Olver 0.2 0.4 0.6 0.8 1 -1-0.50.51 0.2 0.4 0.6 0.8 1 -1-0.50.51 0.2 0.4 0.6 0.8 1 -1-0.50.51 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 Figure 14.15. Numerical Solutions for the Heat Equation Based on the Explicit Scheme. The coeﬃcient matrix Ais symmetric and tridiagonal. The contributions (14.148) o f the boundary nodes appear in the vector b(i). This numerical method is known as an explicit schemesince each iterate is computed directly without relying on s olving an auxiliary equation — unlike the implicit schemes to be discussed below . Example 14.14. Let us ﬁx the diﬀusivity γ= 1 and the bar length ℓ= 1. For illustrative purposes, we take a spatial step size of h=.1. In Figure 14.15 we compare two (slightly) diﬀerent time step sizes on the same initial d ata as used in (14.22). The ﬁrst sequence uses the time step k=h2=.01 and plots the solution at times t= 0.,.02,.04. The solution is already starting to show signs of instabilit y, and indeed soon thereafter becomes completely wild. The second sequence takes k=.005 and plots the solution at timest= 0.,.025,.05. (Note that the two sequences of plots have diﬀerent verti cal scales.) Even though we are employing a rather coarse mesh, the numeri cal solution is not too far away from the true solution to the initialvalue problem, whi ch can be found in Figure 14.1. In light of this calculation, we need to understand why our sc heme sometimes gives reasonable answers but at other times utterly fails. To this end, let us specialize to homo- geneous boundary conditions u(t,0) = 0 = u(t,ℓ),whereby αi=βi= 0 for all i= 0,1,2,3,... , (14.152) and so (14.150) reduces to a homogeneous, linear iterative s ystem u(i+1)=Au(i). (14.153) According to Proposition 10.11, all solutions will converg e to zero, u(i)→0— as they are supposed to (why?) — if and only if Ais a convergent matrix. But convergence depends upon the step sizes. Example 14.14 is indicating that for mes h sizeh=.1, the time step k=.01 yields a non-convergent matrix, while k=.005 leads to a convergent matrix and a valid numerical scheme. 12/11/12 789 c/ci∇cleco√y∇t2012 Peter J. Olver As we learned in Theorem 10.14, the convergence property of a matrix is ﬁxed by its spectral radius, i.e., its largest eigenvalue in magnit ude. There is, in fact, an explicit formula for the eigenvalues of the particular tridiagonal m atrix in our numerical scheme, which follows from the following general result, which solv es Exercise 8.2.48. It is a direct consequence of Exercise 8.2.47, which contains the explici t formulae for the eigenvectors. Lemma 14.15. The eigenvalues of an (n−1)×(n−1)tridiagonal matrix all of whose diagonal entries are equal to aand all of whose sub- and super-diagonal entries are equal to bare λk=a+2bcosπk n, k = 1,...,n−1. (14.154) In our particular case, a= 1−2µandb=µ, and hence the eigenvalues of the matrix Agiven by (14.151) are λk= 1−2µ+2µcosπk n, k = 1,...,n−1. Since the cosine term ranges between −1 and +1, the eigenvalues satisfy 1−4µ < λk<1. Thus, assuming that 0 < µ≤1 2guarantees that all |λk|<1, and hence Ais a convergent matrix. In this way, we have deduced the basic stability crit erion µ=γk h2≤1 2,or k≤h2 2γ. (14.155) With some additional analytical work, [ 109], it can be shown that this is suﬃcient to conclude that the numerical scheme (14.145–148) converges to the true solution to the initial-boundary value problem for the heat equation. Sincenotallchoicesofspaceandtimestepsleadtoaconverg ent scheme, thenumerical method is called conditionally stable . The convergence criterion (14.155) places a severe restriction on the time step size. For instance, if we have h=.01, andγ= 1, then we can only use a time step size k≤.00005, which is minuscule. It would take an inordinately large number of time steps to compute the value of the solutio n at even a moderate times, e.g.,t= 1. Moreover, owing to the limited accuracy of computers, th e propagation of round-oﬀ errors might then cause a signiﬁcant reduction in t he overall accuracy of the ﬁnal solution values. An unconditionally stable method — one that does not restric t the time step — can be constructed by using the backwards diﬀerence formula ∂u ∂t(ti,xj)≈u(ti,xj)−u(ti−1,xj) k+O(hk) (14 .156) toapproximatethetemporalderivative. Substituting(14. 156)andthesameapproximation (14.143) for uxxinto the heat equation, and then replacing ibyi+1, leads to the iterative system ui+1,j−µ/parenleftbig ui+1,j+1−2ui+1,j+ui+1,j−1/parenrightbig =ui,j,i= 0,1,2,... , j= 1,...,n−1,(14.157) 12/11/12 790 c/ci∇cleco√y∇t2012 Peter J. Olver 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 Figure 14.16. Numerical Solutions for the Heat Equation Based on the Implicit Scheme. where the parameter µ=γk/h2is as above. The initial and boundary conditions also have the same form (14.147,148). The system can be written in the matrix form /hatwideAu(i+1)=u(i)+b(i+1), (14.158) where/hatwideAis obtained from the matrix Ain (14.151) by replacing µby−µ. This deﬁnes animplicit method since we have to solve a tridiagonal linear system at each ste p in order to compute the next iterate u(i+1). However, as we learned in Section 1.7, tridiagonal systems can be solved very rapidly, and so speed does not beco me a signiﬁcant issue in the practical implementation of this implicit scheme. Let us look at the convergence properties of the implicit sch eme. For homogeneous Dirichlet boundary conditions (14.152), the system takes t he form u(i+1)=/hatwideA−1u(i), and the convergence is now governed by the eigenvalues of /hatwideA−1. Lemma 14.15 tells us that the eigenvalues of /hatwideAare λk= 1+2µ−2µcosπk n, k = 1,...,n−1. As a result, its inverse /hatwideA−1has eigenvalues 1 λk=1 1+2µ/parenleftbigg 1−cosπk n/parenrightbigg, k = 1,...,n−1. Sinceµ >0, the latter are alwaysless than 1 in absolute value, and so /hatwideAis always a convergent matrix. The implicit scheme (14.158) is converg ent for any choice of step sizes h,k, and hence unconditionally stable . Example 14.16. Consider the same initial-boundary value problem consider ed in Example14.14. InFigure14.16, weplotthenumerical soluti onsobtainedusingtheimplicit scheme. The initial data is not displayed, but we graph the nu merical solutions at times 12/11/12 791 c/ci∇cleco√y∇t2012 Peter J. Olver 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.2 Figure 14.17. Numerical Solutions for the Heat Equation Based on the Crank–Nicolson Scheme. t=.2,.4,.6 with a mesh size of h=.1. On the top line, we use a time step of k=.01, while on the bottom k=.005. Unlike the explicit scheme, there is very little diﬀere nce between the two — both come much closer to the actual solution than the explicit scheme. Indeed, even signiﬁcantly larger time steps give reasonabl e numerical approximations to the solution. Another popular numerical scheme is the Crank–Nicolson method ui+1,j−ui,j=µ 2/parenleftbig ui+1,j+1−2ui+1,j+ui+1,j−1+ui,j+1−2ui,j+ui,j−1/parenrightbig .(14.159) which can be obtained by averaging the explicit and implicit schemes (14.145,157). We can write the iterative system in matrix form Bu(i+1)=Cu(i)+1 2/parenleftbig b(i)+b(i+1)/parenrightbig , where B= 1+µ−1 2µ −1 2µ1+µ−1 2µ −1 2µ...... ...... , C = 1−µ1 2µ 1 2µ1−µ1 2µ 1 2µ...... ...... .(14.160) Convergenceisgovernedby thegeneralizedeigenvaluesoft hetridiagonalmatrixpair B,C, or, equivalently, the eigenvalues of the product B−1C, cf. Exercise 9.5.33. According to Exercise , these are λk=1−µ/parenleftbigg 1−cosπk n/parenrightbigg 1+µ/parenleftbigg 1−cosπk n/parenrightbigg, k = 1,...,n−1. (14.161) Sinceµ >0, all of the eigenvalues are strictly less than 1 in absolute value, and so the Crank–Nicolson scheme is also unconditionally stable. A de tailed analysis will show that 12/11/12 792 c/ci∇cleco√y∇t2012 Peter J. Olver the errors are of the order of k2andh2, and so it is reasonable to choose the time step to have the same order of magnitude as the space step, k≈h. This gives the Crank–Nicolson scheme one advantage over the previous two methods. However , applying it to the initial value problem considered earlier points out a signiﬁcant we akness. Figure 14.17 shows the result of running the scheme on the initial data (14.22). The top row has space and time step sizes h=k=.1, and does a rather poor job replicating the solution. The se cond row usesh=k=.01, and performs better except near the corners where an anno ying and incorrect local time oscillation persists as the solution d ecays. Indeed, since most of its eigenvalues are near −1, the Crank–Nicolson scheme does not do a good job of damping out the high frequency modes that arise from small scale feat ures, including discontinuities and corners in the initial data. On the other hand, most of the eigenvalues of the fully implicit scheme are near zero, and it tends to handle the high frequency modes better, losing out to Crank–Nicolson when the data is smooth. Thus, a good strategy is to ﬁrst evolve using the implicit scheme until the small scale noise is dissipated away, and then switch to Crank–Nicolson to use a much larger time step for ﬁn al the large scale changes. Numerical Solution Methods for the Wave Equation Let us now look at some numerical solution techniques for the wave equation. Al- though this is in a sense unnecessary, owing to the explicit d ’Alembert solution for- mula (14.124), the experience we gain in designing workable schemes will serve us well in more complicated situations, including inhomogeneous m edia, and higher dimensional problems, when analytic solution formulas are no longer ava ilable. Consider the wave equation ∂2u ∂t2=c2∂2u ∂x2, 0< x < ℓ, t ≥0, (14.162) modeling vibrations of a homogeneous bar of length ℓwith constant wave speed c >0. We impose Dirichlet boundary conditions u(t,0) =α(t), u (t,ℓ) =β(t), t ≥0. (14.163) and initial conditions u(0,x) =f(x),∂u ∂t(0,x) =g(x), 0≤x≤ℓ. (14.164) We adopt the same uniformly spaced mesh ti=ik, xj=jh,where h=ℓ n. In order to discretize the wave equation, we replace the seco nd order derivatives by their standard ﬁnite diﬀerence approximations (14.137), n amely ∂2u ∂t2(ti,xj)≈u(ti+1,xj)−2u(ti,xj)+u(ti−1,xj) k2+ O(k2), ∂2u ∂x2(ti,xj)≈u(ti,xj+1)−2u(ti,xj)+u(ti,xj−1) h2+ O(h2),(14.165) 12/11/12 793 c/ci∇cleco√y∇t2012 Peter J. Olver Since the errors are of orders of k2andh2, we anticipate to be able to choose the space and time step sizes of comparable magnitude: k≈h. Substituting the ﬁnite diﬀerence formulae (14.165) into th e partial diﬀerential equation (14.162), and rearranging terms, we are led to the iterative system ui+1,j=σ2ui,j+1+2(1−σ2)ui,j+σ2ui,j−1−ui−1,j,i= 1,2,... , j= 1,...,n−1,(14.166) for the numerical approximations ui,j≈u(ti,xj) to the solution values at the mesh points. The positive parameter σ=ck h>0 (14 .167) depends upon the wave speed and the ratio of space and time ste p sizes. The boundary conditions (14.163) require that ui,0=αi=α(ti), ui,n=βi=β(ti), i = 0,1,2,... . (14.168) This allows us to rewrite the system in matrix form u(i+1)=Bu(i)−u(i−1)+b(i), (14.169) where B= 2(1−σ2)σ2 σ22(1−σ2)σ2 σ2...... ......σ2 σ22(1−σ2) ,u(i)= ui,1 ui,2 ... ui,n−2 uii,n−1 ,b(i)= σ2αi 0 ... 0 σ2βi . (14.170) The entries of u(i)are, as in (14.149), the numerical approximations to the sol ution values at theinterior nodes. Note that the system (14.169) is a second order iterat ive scheme, since computing the next iterate u(i+1)requires the value of the preceding two, u(i)and u(i−1). The one diﬃculty is getting the method started. We know u(0)sinceu0,j=fj=f(xj) is determined by the initial position. However, we also need to ﬁndu(1)with entries u1,j≈u(k,xj) at time t1=kin order launch the iteration, but the initial velocity ut(0,x) =g(x) prescribes the derivatives ut(0,xj) =gj=g(xj) at time t0= 0 instead. One way to resolve this diﬃcult would be to utilize the ﬁnite d iﬀerence approximation gj=∂u ∂t(0,xj)≈u(k,xj)−u(0,xj) k≈u1,j−fj k(14.171) to compute the required values u1,j=fj+kgj. 12/11/12 794 c/ci∇cleco√y∇t2012 Peter J. Olver However, the approximation (14.171) is only accurate to ord erk, whereas the rest of the scheme has error proportional to k2. Therefore, we would introduce an unacceptably large error at the initial step. To construct an initial approximation to u(1)with error on the order of k2, we need to analyze the local error in more detail. Note that, by Taylo r’s theorem, u(k,xj)−u(0,xj) k≈∂u ∂t(0,xj)+k 2∂2u ∂t2(0,xj) =∂u ∂t(0,xj)+c2k 2∂2u ∂x2(0,xj), where the error is now of order k2, and, in the ﬁnal equality, we have used the fact that u is a solution to the wave equation. Therefore, we ﬁnd u(k,xj)≈u(0,xj)+k∂u ∂t(0,xj)+c2k2 2∂2u ∂x2(0,xj) =f(xj)+kg(xj)+c2k2 2f′′(xj)≈fj+kgj+c2k2 2h2(fj+1−2fj+fj−1), where we can use the ﬁnite diﬀerence approximation (14.137) for the second derivative of f(x) if no explicit formula is known. Therefore, when we initiat e the scheme by setting u1,j=1 2σ2fj+1+(1−σ2)fj+1 2σ2fj−1+kgj, (14.172) or, in matrix form, u(0)=f, u(1)=1 2Bu(0)+kg+1 2b(0), (14.173) we will have maintained the desired order k2(andh2) accuracy. Example 14.17. Consider the particular initial value problem utt=uxx,u(0,x) =e−400(x−.3)2, ut(0,x) = 0, u(t,0) =u(1,0) = 0,0≤x≤1, t≥0, subject to homogeneous Dirichlet boundary conditions on th e interval [0 ,1]. The initial data is a fairly concentrated single hump centered at x=.3, and we expect it to split into two halfsized humps, which then collidewiththe ends. Let us choose a space discretization consisting of 90 equally spaced points, and so h=1 90=.0111.... If we choose a time step ofk=.01, whereby σ=.9, then we get reasonably accurate solution over a fairly lon g time range, as plotted in Figure 14.18 at times t= 0,.1,.2,...,.5. On the other hand, if we double the time step, setting k=.02, soσ= 1.8, then, as plotted in Figure 14.19 at times t= 0,.05,.1,.14,.16,.18, we observe an instability eventually creeping into the picture that eventually overwhelms the numerical solution . Thus, the numerical scheme appears to only be conditionally stable. The stability analysis of this numerical scheme proceeds as follows. We ﬁrst need to recast the second order iterative system (14.169) into a ﬁrs t order system. In analogy with Example 10.6, this is accomplished by introducing the vecto rz(i)=/parenleftbigg u(i) u(i−1)/parenrightbigg ∈R2n−2. Then z(i+1)=Cz(i)+c(i),where C=/parenleftbigg B−I I O/parenrightbigg . (14.174) 12/11/12 795 c/ci∇cleco√y∇t2012 Peter J. Olver 0.2 0.4 0.6 0.8 1 -1-0.75-0.5-0.250.250.50.751 0.2 0.4 0.6 0.8 1 -1-0.75-0.5-0.250.250.50.751 0.2 0.4 0.6 0.8 1 -1-0.75-0.5-0.250.250.50.751 0.2 0.4 0.6 0.8 1 -1-0.75-0.5-0.250.250.50.751 0.2 0.4 0.6 0.8 1 -1-0.75-0.5-0.250.250.50.751 0.2 0.4 0.6 0.8 1 -1-0.75-0.5-0.250.250.50.751 Figure 14.18. Numerically Stable Waves. 0.2 0.4 0.6 0.8 1 -1-0.75-0.5-0.250.250.50.751 0.2 0.4 0.6 0.8 1 -1-0.75-0.5-0.250.250.50.751 0.2 0.4 0.6 0.8 1 -1-0.75-0.5-0.250.250.50.751 0.2 0.4 0.6 0.8 1 -1-0.75-0.5-0.250.250.50.751 0.2 0.4 0.6 0.8 1 -1-0.75-0.5-0.250.250.50.751 0.2 0.4 0.6 0.8 1 -1-0.75-0.5-0.250.250.50.751 Figure 14.19. Numerically Unstable Waves. Therefore, the stability of the method will be determined by the eigenvalues of the coeﬃ- cient matrix C. The eigenvector equation Cz=λz, wherez=/parenleftbigg u v/parenrightbigg , can be written out in its individual components: Bu−v=λu,u=λv. Substituting the second equation into the ﬁrst, we ﬁnd (λB−λ2−1)v=0,or Bv=/parenleftbigg λ+1 λ/parenrightbigg v. The latter equation implies that vis an eigenvector of Bwithλ+λ−1the corresponding eigenvalue. The eigenvalues of the tridiagonal matrix Bare governed by Lemma 14.15, in whicha= 2(1−σ2) andb=σ2, and hence are λ+1 λ= 2/parenleftbigg 1−σ2+σ2cosπk n/parenrightbigg , k = 1,...,n−1. Multiplying both sides by λleads to a quadratic equation for the eigenvalues, λ2−2akλ+1 = 0,where 1 −2σ2< ak= 1−σ2+σ2cosπk n<1.(14.175) 12/11/12 796 c/ci∇cleco√y∇t2012 Peter J. Olver Figure 14.20. The Courant Condition. Each pair of solutions to these n−1 quadratic equations, namely λ± k=ak±/radicalig a2 k−1, (14.176) yields two eigenvalues of the matrix C. Ifak<−1, then one of the two eigenvalues will be both<−1, which means that the linear iterative system has an expone ntially growing mode, and so /ba∇dblu(i)/ba∇dbl → ∞asi→ ∞for almost all choices of initial data. This is clearly incompatible with the wave equation solution that we are try ing to approximate, which is periodic and hence remains bounded. On the other hand, if |ak|<1, then the eigenvalues (14.176) are complex numbers of modulus 1, indicated stabil ity (but not convergence) of the matrix C. Therefore, in view of (14.175), we should require that σ=ck h<1,ork <h c, (14.177) which places a restriction on the relative sizes of the time a nd space steps. We conclude that the numerical scheme is conditionally stable. The stability criterion (14.177) is known as the Courant condition , and can be as- signed a simple geometric interpretation. Recall that the w ave speed cis the slope of the characteristic lines for the wave equation. The Courant con dition requires that the mesh slope, which is deﬁned to be the ratio of the space step size to the ti me step size, namely h/k, must be strictly greater than the characteristic slope c. A signal starting at a mesh point (ti,xj) will reach positions xj±k/cat the next time ti+1=ti+k, which are still between the mesh points xj−1andxj+1. Thus, characteristic lines that start at a mesh point are not allowed to reach beyond the neighboring mesh po ints at the next time step. For instance, in Figure 14.20, the wave speed is c= 1.25. The ﬁrst ﬁgure has equal mesh spacing k=h, and does not satisfy the Courant condition (14.177), where as the second ﬁgure has k=1 2h, which does. Note how the characteristic lines starting at a given mesh point have progressed beyond the neighboring mes h points after one time step in the ﬁrst case, but not in the second. 14.7. A General Framework for Dynamics. According to Section 12.1, the one-dimensional heat and wav e equations are speciﬁc 12/11/12 797 c/ci∇cleco√y∇t2012 Peter J. Olver instances of two broad classes of dynamical systems that inc lude, in a common framework, both discrete dynamics modeled by systems of ordinary diﬀer ential equations and contin- uum systems modeled by (systems of) partial diﬀerential equ ations. In preparation for their multi-dimensional generalizations, it will be usefu l to summarize the general mathe- matical framework, which can be regarded as the dynamical co unterpart to the framework for equilibrium developed in Sections 7.5 and 15.4, which yo u are advised to review before going through this section. Readers more attuned to concret e examples, on the other hand, might prefer to skip this material entirely, referring back as necessary. In all situations, the starting point is a linear function L:U−→V (14.178) that maps a vector space Uto another vector space V. In mechanics, the elements of Urepresent displacements, while the elements of Vrepresent strains (elongations). In electromagnetism and gravitation, elements of Urepresent potentials and elements of V electric or magnetic or gravitational ﬁelds. In thermodyna mics, elements of Urepresent temperature distributions, and elements of Vtemperature gradients. In ﬂuid mechanics, Ucontains potential functions and Vis ﬂuid velocities. And so on. In the discrete, ﬁnite-dimensional setting, when U=RnandV=Rm, the linear func- tionL[u] =Auis represented by multiplication by an m×nmatrixA— the incidence matrix and its generalizations. In the inﬁnite-dimensiona l function space situation pre- sented in Chapter 11 and earlier in this chapter, the linear m apLis a diﬀerential operator; in the case of one-dimensional elastic bars, it is the ﬁrst or der derivative L=Dx, while for beams it is a second order derivative L=D2 x. In the multi-dimensional physics treated in subsequent chapters, the most important example is the grad ient operator, L=∇, that maps scalar potentials to vector ﬁelds. The vector spaces UandVare further each endowed with inner products, that en- capsulate the consitutive assumptions underlying the phys ical problem. Deﬁnition 7.53 explains how to construct the resulting adjoint map L∗:V−→U (14.179) which goess in the reversedirection. In ﬁnite dimensions, when U=RnandV=Rm are both equipped with the Euclidean dot product, the adjoin t corresponds to the sim- ple transpose operation. Modiﬁcations under more general w eighted inner products were discussed in Section 7.5. In inﬁnite-dimensional contexts , the computation of the adjoint presented inSection 11.3reliedonanintegrationby partsa rgument, supplemented by suit- able boundary conditions. In the next chapter, we will adapt this argument to determine adjoints of multi-dimensional diﬀerential operators. The crucial operator underlying the equilibrium and dynami cal equations for a re- markably broad range of physics is the self-adjoint combina tion K=L∗◦L:U−→U. (14.180) According to Theorem 7.60, the operator Kispositive semi-deﬁnite in all situations, and positive deﬁnite if and only if ker L={0}. In the ﬁnite-dimensional context, Kis repre- sented by the symmetric positive (semi-)deﬁnite Gram matri xATAwhen both U=Rn 12/11/12 798 c/ci∇cleco√y∇t2012 Peter J. Olver andV=Rmhave the dot product, by the symmetric combination ATCAwhenVhas a weighted inner product represented by the positive deﬁnite matrixC >0, and the more general self-adjoint form M−1ATCAwhenUalso is endowed with a weighted inner prod- uct represented by M >0. In one-dimensional bars, Kis represented by a self-adjoint second order diﬀerential operator, whose form depends upon the underlying inner prod- ucts., while in beams it becomes a fourth order diﬀerential o perator. The deﬁniteness of K (or lack thereof) depends upon the imposed boundary conditi ons. In higher dimensions, as we discuss in Section 15.4, Kbecomes the Laplacian or a more general elliptic partial diﬀerential operator. With this set-up, the basic equilibrium equation has the form K[u] =f, (14.181) wherefrepresents an external forcing function. In ﬁnite dimensio ns, this isa linear system consisting of nequations in nunknowns with positive (semi-)deﬁnite coeﬃcient matrix. In function space, it becomes a self-adjoint boundary value problem for the unknown function u. IfKis positive deﬁnite, the solution is unique. (But rigorousl y proving the existence of a solution is not trivial, requiring serious an alysis beyond the scope of this text.) Theorem 7.61 says that the solutions are characteriz ed by a minimization principle, as the minimizers of the quadratic function(al) p[u] =1 2/ba∇dblL[u]/ba∇dbl2−/an}b∇acketle{tu;f/an}b∇acket∇i}ht, (14.182) which typically represents the potential energy in the syst em. IfKis only positive semi- deﬁnite, existence of a solution requires that the forcing f unction satisfy the Fredholm condition that it be orthogonal to the the unstable modes, th at is the elements of ker K= kerL With the equilibrium operator Kin hand, there are two basic types of dynamical systems of importance in physical models. Unforced diﬀusion processes are modeled by a dynamical system of the form ut=−K[u],where K=L∗◦L (14.183) is the standard self-adjoint combination of a linear operat orL. In the discrete case, Kis a matrix and so this represents a ﬁrst order system of ordinar y diﬀerential equations, that has the form of a linear gradient ﬂow (9.18), so named because it decreases the energy function q[u] =/ba∇dblL[u]/ba∇dbl2, (14.184) which is (14.182) when f= 0, as rapidly as possible. In the continuous case, Kis a diﬀerential operator, and (14.183) represents a partial di ﬀerential equation for the time- varying function u=u(t,x). The solution to the general diﬀusion equation (14.183) mimi cs earlier our separation of variables method for the one-dimensional heat equation, as well as the original solution technique for linear systems of ordinary diﬀerential equat ions, as developed in Chapter 9. The separable solutions are of exponential form u(t) =e−λtv, (14.185) 12/11/12 799 c/ci∇cleco√y∇t2012 Peter J. Olver wherev∈Uis a ﬁxed element of the domain space — i.e., a vector in the dis crete context, or afunction v=v(x)that onlydepends onthespatialvariablesinthe continuum versions. Since the operator Kis linear and does not involve tdiﬀerentiation, we ﬁnd ∂u ∂t=−λe−λtv, while K[u] =e−λtK[v]. Substituting back into (14.183) and canceling the common ex ponential factors, we are led to the eigenvalue problem K[v] =λv. (14.186) Thus,vmust be an eigenvector/eigenfunction for the linear operat orK, withλthe corre- sponding eigenvalue. Generalizing our earlier observations on the eigenvalues o f positive deﬁnite matrices and theboundary valueproblems associatedwiththe one-dim ensional heat equation, letus establish the positivity of eigenvalues of such general sel f-adjoint, positive (semi-)deﬁnite linear operators. Theorem 14.18. All eigenvalues of the linear operator K=L∗◦Lare real and non- negative:λ≥0. If, moreover, Kis positive deﬁnite, or, equivalently, kerK= kerL= 0, then all eigenvalues are strictly positive :λ >0. Proof: Suppose K[u] =λuwithu/ne}ationslash= 0. Then λ/ba∇dblu/ba∇dbl2=λ/an}b∇acketle{tu;u/an}b∇acket∇i}ht=/an}b∇acketle{tK[u];u/an}b∇acket∇i}ht=/an}b∇acketle{tL∗◦L[u];u/an}b∇acket∇i}ht=/an}b∇acketle{tL[u];L[u]/an}b∇acket∇i}ht=/ba∇dblL[u]/ba∇dbl2≥0, by the deﬁning equation (7.74) of the adjoint operator. Sinc e/ba∇dblu/ba∇dbl2>0, this immediately implies that λ≥0. Furthermore, in the positivedeﬁnite case ker L={0}, and soL[u]/ne}ationslash= 0. Thus,/ba∇dblL[u]/ba∇dbl2>0, proving that λ >0. Q.E.D. We index the eigenvalues in increasing order: 0< λ1≤λ2≤λ3≤ ···. (14.187) where the eigenvalues are repeated according to their multi plicities, and λ0= 0 is an eigenvalue only in the positive semi-deﬁnite case. Each eig envalue induces a separable eigensolution uk(t) =e−λktvk (14.188) to the diﬀusion equation (14.183). Those associated with st rictly positive deﬁnite eigen- values are exponentially decaying, at a rate equal to the eig envalueλk>0, while any null eigenvalue modes correspond to constant solutions u0(t)≡v0wherev0is any null eigenvector/eigenfunction, i.e., element of ker K= kerL. The general solution is built up as a linear combination u(t) =/summationdisplay kckuk(t) =/summationdisplay kcke−λktvk (14.189) of the eigensolutions. Thus, all solutions decay exponenti ally fast as t→ ∞to an equilib- riumsolution, whichis0inthepositivedeﬁnitecase, or, mo regenerally, thenull eigenspace 12/11/12 800 c/ci∇cleco√y∇t2012 Peter J. Olver component c0v0. The overall rate of decay is, generally, prescribed by the s mallest positive eigenvalue λ1>0. In the discrete version, the summation (14.189) has only ﬁni tely many terms, cor- responding to the neigenvalues of the matrix representing K. Moreover, thanks to the Spectral Theorem 8.26, the eigensolutions form a basis for t he solution space to the dif- fusion equation. In inﬁnite-dimensional function space, t here are, in many instances, an inﬁnite number of eigenvalues, with λk→ ∞ask→ ∞. Completeness of the resulting eigensolutions is a more delicate issue. Often, as in the one -dimensional heat equation, the eigenfunctions are complete and the (Fourier) series (1 4.189) converges and represents the general solution. But in situations involving unbounde d domains, like the hydrogen atomtobediscussed inSection18.7, thereareadditionalse parablesolutionscorresponding to the so-called continuous spectrum that are not represented in terms of eigenfunctions, and the analysis is considerably more involved, requiring a nalogs of the Fourier transform. A full discussion of completeness and convergence of eigenf unction expansions must be relegated to an advanced course in analysis, [ 47,154]. Assuming completeness, the eigensolution coeﬃcients ckin (14.189) are prescribed by the initial conditions, which require /summationdisplay kckvk=h, (14.190) whereh=u(0) represents the initial data. To compute the coeﬃcients ckin the eigenfunc- tion expansion (14.190), we appeal, as in the case of ordinar y Fourier series, to orthogo- nality of the eigenvectors/eigenfunctions vk. Orthogonality is proved by a straightforward adaptation of our earlier proof of part ( b) of Theorem 8.20, guaranteeing the orthogonality of the eigenvectors of a symmetric matrix. Theorem 14.19. Two eigenvectors/eigenfunctions u,vof the self-adjoint linear op- eratorK=L∗◦Lthat are associated with distinct eigenvalues λ/ne}ationslash=µare orthogonal. Proof: Self-adjointness of the operator Kmeans that /an}b∇acketle{tK[u];v/an}b∇acket∇i}ht=/an}b∇acketle{tL[u];L[v]/an}b∇acket∇i}ht=/an}b∇acketle{tu;K[v]/an}b∇acket∇i}ht for any for any u,v. In particular, if K[u] =λu, K[v] =µv,are eigenfunctions, then λ/an}b∇acketle{tu;v/an}b∇acket∇i}ht=/an}b∇acketle{tK[u];v/an}b∇acket∇i}ht=/an}b∇acketle{tu;K[v]/an}b∇acket∇i}ht=µ/an}b∇acketle{tu;v/an}b∇acket∇i}ht, which, assuming λ/ne}ationslash=µ, immediately implies orthogonality: /an}b∇acketle{tu;v/an}b∇acket∇i}ht= 0. Q.E.D. If an eigenvalue λadmits more than one independent eigenvector/eigenfuncti on, we can apply the Gram–Schmidt process to produce an orthogonal basis of its eigenspace† Vλ= ker(K−λI). In this way, the entire set of eigenvectors/eigenfuncti ons can be assumed to be mutually orthogonal: /an}b∇acketle{tvj;vk/an}b∇acket∇i}ht= 0, j /ne}ationslash=k. †This assumes that the eigenspace Vλis ﬁnite-dimensional — which is assured whenever the eigenvalues λk→ ∞ask→ ∞. 12/11/12 801 c/ci∇cleco√y∇t2012 Peter J. Olver As a consequence, taking the inner product of both sides of (1 4.190) with the eigenfunction vkleads to the equation ck/ba∇dblvk/ba∇dbl2=/an}b∇acketle{th;vk/an}b∇acket∇i}htand hence ck=/an}b∇acketle{th;vk/an}b∇acket∇i}ht /ba∇dblvk/ba∇dbl2. (14.191) In this manner we recover our standard orthogonality formul a (5.7) for expressing elements of a vector space in terms of an orthogonal basis. The second important class of dynamical systems consists of second order (in time) vibration equations utt=−K[u]. (14.192) Newton’s equations of motion in the absence of non-conserva tive frictional forces, the propagation of waves in ﬂuids and solids, as well as electrom agnetic waves, and many other relatedphysical problemsleadtosuch vibrationalsy stems. Inthiscase, theseparable solutions are of trigonometric form u(t,x) = cos(ωt)v(x) or sin( ωt)v(x). (14.193) Substituting this ansatz back into the vibration equation ( 14.192) results in the same eigenvalue problem (14.186) with eigenvalue λ=ω2equal to the square of the vibrational frequency. We conclude that the normal mode or separable eigensolutions take the form uk(t) = cos(ωkt)vk,/tildewideuk(t) = sin(ωkt)vk,provided λk=ω2 k>0 is a non-zero eigenvalue. In the stable, positive deﬁnite ca se, there are no zero eigenvalues, and so the general solution is built up as a (quasi-)periodic†combination u(t) =/summationdisplay k/bracketleftbig ckuk(t)+dk/tildewideuk(t)/bracketrightbig =/summationdisplay krkcos(ωkt+δk)vk, (14.194) of the eigenmodes. The initial conditions g=u(0) =/summationdisplay kckvk, h =ut(0) =/summationdisplay kdkωkvk,(14.195) are used to specify the coeﬃcients ck,dk, using the same orthogonality formula (14.191): ck=/an}b∇acketle{tf;vk/an}b∇acket∇i}ht /ba∇dblvk/ba∇dbl2, dk=/an}b∇acketle{tf;vk/an}b∇acket∇i}ht ωk/ba∇dblvk/ba∇dbl2. (14.196) In the unstable, positive semi-deﬁnite cases, the null eige nsolutions have the form u0(t) =v0,/tildewideu0(t) =tv0, wherev0∈kerK= kerL, and must be appended to the series solution. The unstable mode/tildewideu0(t) is excited if and only if the initial velocity is not orthogo nal to the kernel element: /an}b∇acketle{th;v0/an}b∇acket∇i}ht /ne}ationslash= 0. †The solution is periodic if and only if the frequencies appearing in the sum are all integer multiples of a common frequency: ωk=nkω⋆fornk∈N. 12/11/12 802 c/ci∇cleco√y∇t2012 Peter J. Olver In classical mechanics, the diﬀusion and vibration equatio ns and their variants (see, for instance, Exercises –for versions with external forcing and Exercise for vibrations with frictional eﬀects) are the most important classes of dy namical systems. In quantum mechanics, the basic dynamical system is known as the Schr¨ odinger equation , ﬁrst written down by the the German physicist Erwin Schr¨ odinger, one of t he founders of quantum mechanics. The abstract form of the Schr¨ odinger equation i s i/planckover2pi1ut=K[u]. (14.197) Herei=√−1, while /planckover2pi1=h 2π≈1.055×10−34Joule seconds (14 .198) isPlanck’s constant , whose value governs the quantization of all physical quant ities. At each time t, the solution u(t,x) to the Schr¨ odinger equation represents the wave function of the quantum system, and so is a complex-valued square inte grable function of constant L2norm:/ba∇dblu/ba∇dbl= 1. (See Sections 12.5 and 13.3 for the basics of quantum mech anics and Hilbert space.) As usual, we interpret the wave fuction as a p robability density on the possible quantum states, and so the Schr¨ odinger equation g overns the dynamical evolution of quantum probabilities. The operator K=L∗◦Lis known as the Hamiltonian for the quantum mechanical system governed by (14.197), and, typic ally represents the quantum energy operator. For physical systems such as atoms and nucl ei, the relevant Hamiltonian operator is constructed from the classical energy through t he rather mysterious process of “quantization”. The interested reader should consult a bas ic text on quantum mechanics, e.g., [124,130], for full details on both the physics and underlying mathem atics. Proposition 14.20. Ifu(t)is a solution to the Schr¨ odinger equation, its Hermitian L2norm/ba∇dblu/ba∇dblis constant. Proof: Since the solution is complex-valued, we use the sesquilin earity of the under- lying Hermitian inner product, as in (3.94), to compute d dt/ba∇dblu/ba∇dbl2=/an}b∇acketle{tut;u/an}b∇acket∇i}ht+/an}b∇acketle{tu;ut/an}b∇acket∇i}ht =/angbracketleftbigg −i /planckover2pi1K[u];u/angbracketrightbigg +/angbracketleftbigg u;−i /planckover2pi1K[u]/angbracketrightbigg =−i /planckover2pi1/an}b∇acketle{tK[u];u/an}b∇acket∇i}ht+i /planckover2pi1/an}b∇acketle{tu;K[u]/an}b∇acket∇i}ht= 0, which vanishes since Kis self-adjoint. Since itsderivativevanishes everywhere , thisimplies that/ba∇dblu/ba∇dbl2is constant. Q.E.D. As a result, if the initial data u(t0) =his a quantum mechanical wave function, meaning that /ba∇dblh/ba∇dbl= 1, then, at each time t, the solution to the Schr¨ odinger equation also has norm 1, and hence remains a wave function for all t. Apart from the extra factor of i /planckover2pi1, the Schr¨ odinger equation looks like a diﬀusion equation (14.183). ( Warning : Despite this superﬁcial similarity, their solutions have radi- cally diﬀerent behavior.) This inspires us to seek separabl e solutions with an exponential ansatz: u(t,x) =eαtv(x). 12/11/12 803 c/ci∇cleco√y∇t2012 Peter J. Olver Substituting this expression into the Schr¨ odinger equati on (14.197) and canceling the com- mon exponential factors reduces us to the usual eigenvalue p roblem K[v] =λv, with eigenvalue λ=−i/planckover2pi1α. Letvk(x) denote the normalized eigenfunction associated with the ktheigenvalue λk. Th corresponding separablesolutionoftheSchr¨ odinger equa tionisthecomplex wavefunctions uk(t,x) =eiλkt//planckover2pi1vk(x). Observe that, in contrast to the exponentially decaying sol utions to the diﬀusion equation, the eigensolutions to the Schr¨ odinger equation are period ic, of frequencies proportional to the eigenvalues: ωk=λk//planckover2pi1. (Along with constant solutions corresponding to the null eigenmodes, if any.) The general solution is a (quasi-)peri odic series in the fundamental eigensolutions. The periodicity of the summands has the add itional implicationthat, again unlike the diﬀusion equation, the Schr¨ odinger equation ca n be run backwards in time. So, we can ﬁgure out both the past and future behavior of a quantum system from its present conﬁguration. Example 14.21. Ina singlespacedimension, thesimplest versionoftheSchr ¨ odinger equation is based on the derivative operator L=Dx, for which, assuming appropriate boundary conditions, the self-adjoint combination K=L∗◦L=−D2 x. Thus, (14.197) reduces to the second order partial diﬀerential equation i/planckover2pi1ut=−uxx. (14.199) Imposing the Dirichlet boundary conditions u(t,0) =u(t,ℓ) = 0, the Schr¨ odinger equation governs the dynamics of a quantum particle conﬁned to the interval 0 < x < ℓ ; the boundary conditions imply that there is zero probabili ty of the particle escaping from the interval. According to Section 14.1, the eigenfunctions of the Dirich let eigenvalue problem vxx+λv= 0, v (0) =v(ℓ) = 0, are vk(x) =/radicalbigg 2 ℓsinkπ ℓxfork= 1,2,... , with eigenvalue λk=k2π2 ℓ2, where the initial factor is ensures that vkhas unit L2norm, and hence is a bona ﬁde wave function. The corresponding separable solutions or eigenm odes are uk(t,x) =/radicalbigg 2 ℓexp/parenleftbigg ik2π2 /planckover2pi1ℓ2t/parenrightbigg sinkπ ℓx. (14.200) The eigenvalues represent the energy levels of the particle , which can be observed from the spectral lines emitted by the system. For instance, when an electron jumps from 12/11/12 804 c/ci∇cleco√y∇t2012 Peter J. Olver one level to another, it conserves energy by emitting a photo n, with energy equal to the diﬀerence in energy between the two quantum levels. These em itted photons deﬁne the observed electromagnetic spectral lines, hence the adopti onofthe physics term “spectrum” to describe the eigenvalues of the Hamiltonian operator K. 12/11/12 805 c/ci∇cleco√y∇t2012 Peter J. Olver