hwy

Chapter 17 DynamicsofPlanarMedia In this chapter, we continue our ascent of the dimensional la dder for linear systems. In Chapter 6, we embarked on our journey with equilibrium con figurations of discrete systems — mass–spring chains, circuits, and structures — wh ich are governed by certain linear algebraic systems. In Chapter 9, the dynamical behav ior of such discrete systems was modeled by systems of linear ordinary differential equat ions. Chapter 11 began our treatment of continuous media with the boundary value probl ems that describe the equi- libria of one-dimensional bars, strings and beams. Their dy namical motions formed the topic of Chapter 14, in the simplest case leading to two funda mental partial differential equations: the heat equation describing thermal diffusion, and the wave equation model- ing vibrations. In Chapters 15 and 16, we focussed our attent ion on the boundary value problems describing equilibria of planar bodies — plates an d membranes — with primary emphasis on solving the ubiquitous Laplace equation, both a nalyticallyor numerically. We now turn to the analysis of their dynamics, as governed by the two-dimensional†forms of the heat and wave equations. The heat equation describes diff usion of, say, heat energy, or population, or pollutants in a homogeneous two-dimensiona l domain. The wave equation models small vibrations of two-dimensional membranes such as a drum. Although the increase in dimension does challenge our analy tical prowess, we have, in fact, already mastered the key techniques: separation of va riables and fundamental solu- tions. (Disappointingly, conformal mappings are not parti cularly helpful in the dynamical universe.) When applied to partial differential equation in higher dimensions, separation of variables often leads to new ordinary differential equation s, whose solutions are no longer elementary functions. These so-called special functions , which include the Bessel functions appearing in the present chapter, and the Legendre function s, spherical harmonics, and spherical Bessel functions in three-dimensional problems , play a ubiquitous role in more advanced applications in physics, engineering and mathema tics. Basic series solution tech- niques for ordinary differential equations, and key propert ies of the most important classes of special functions, can be found in Appendix C In Appendix C, we collect together the required results abou t the most important classes of special functions, including a short presentati on of the series approach for solving non-elementary ordinary differential equations. †Throughout, “dimension” refers to the number of space variables. In Ne wtonian dynamics, the time “dimension” is accorded a separate status, which distinguis hes dynamics from equilib- rium. Of course, in the more complicated relativistic universe, t ime and space must be regarded on an equal footing, and the dimension count modified accordinagly. 12/11/12 936 c/ci∇cleco√y∇t2012 Peter J. Olver Numerical methods for solving boundary value and initial va lue problems are, of course, essential in all but the simplest situations. The tw o basic methods — finite element and finite difference — have already appeared, and the only new aspect is the (substantial) complicationofworkinginhigher dimensions. Thus, inthei nterestsofbrevity,wedefer the discussion of the numerical aspects of multi-dimensional p artial differential equations to more advanced texts, e.g., [ 109], and student projects outlined in the exercises. However, the student should be assured that, without knowledge of the qualitative features based on direct analysis and explicit solutions, the design, implem entation, and testing of numerical solution techniques would be severely hampered. 17.1. Diffusion in Planar Media. As we learned in Chapter 15, the equilibrium temperature u(x,y) of a homogeneous plate is governed by the two-dimensional Laplace equation ∆u=uxx+uyy= 0. In conformity with our general framework, the dynamical diff usion of such a plate will be modeled by the two-dimensional heat equation ut=γ∆u=γ/parenleftbig uxx+uyy/parenrightbig , (17.1) where the diffusivity coefficient γ >0 measures the relative speed of diffusion of heat energy throughout the plate; its positivity is required on p hysical grounds, and also avoids ill-posedness of the dynamical system. In this simplest mod el of two-dimensional diffusion, weareassumingthattherearenolossofheatorexternalheat sourcesintheplate’sinterior, which can be arranged by covering it with insulation. The solution u(t,x) =u(t,x,y) to (17.1) measures the temperature at time tat each pointx= (x,y)∈Ω in the domain Ω ⊂R2occupied by the plate. To uniquely specify u(t,x,y) at each point ( x,y)∈Ω and each positive t >0, we must impose both initial and boundary conditions. The most important are: (a)Dirichlet boundary conditions : u=h on ∂Ω, (17.2) which fix the temperature on the boundary of the plate. (b)Neumann boundary conditions : ∂u ∂n=∂u ∂n=k on ∂Ω, (17.3) that prescribe the heat flux along the boundary; the case k= 0 corresponds to an insulated boundary. (c)Mixed boundary conditions : we impose Dirichlet conditions on part of the boundary D/subsetnoteql∂Ω and Neumann conditions on the remainder N=∂Ω\D. For instance, the homogeneous mixed boundary conditions u= 0 on D,∂u ∂n= 0 on N, (17.4) 12/11/12 937 c/ci∇cleco√y∇t2012 Peter J. Olver correspond to freezing part of the boundary and insulating t he remainder. The initial conditions specify the temperature of the plate u(0,x,y) =f(x,y), (x,y)∈Ω, (17.5) at an initial time, which for simplicity, we take as t0= 0. Under reasonable assumptions, e.g., the domain Ω is bounded, its boundary is piecewise smoo th, and the boundary data is piecewise continuous, say, a general theorem, [ 47], guarantees the existence of a unique solution u(t,x,y) to any of these initial-boundary value problems for all sub sequent times t >0. Our practical goal is to both compute and understand the be havior of the solution in specific situations. Derivation of the Diffusion Equation This section is for those interested in understanding how th e heat equation arises as a model for heat flow. The physical derivation of the two-di mensional (and three- dimensional)heat equationreliesuponthesametwo basicth ermodynamicallawsthatwere used, in Section 14.1, to establish the one-dimensional ver sion. The first principle is that heat energy flows from hot to cold as rapidly as possible. Acco rding to the multivariable calculus Theorem 19.38, the negative temperature gradient −∇upoints in the direction of the steepest decrease in the function uat a point, and so heat energy will flow in that direction. Therefore, the heat flux vector w, which measures the magnitude and direction of the flow of heat energy, should be proportional to the tempe rature gradient: w(t,x,y) =−κ(x,y)∇u. (17.6) Thescalar quantity κ(x,y)>0 measures the thermal conductivity ofthematerial, so(17.6) is the multi-dimensional form of Fourier’s Law of Cooling (14.3). We are assuming that the thermal conductivity depends only on the position ( x,y)∈Ω, which means that the material in the plate (a) is not changing in time; (b) isisotropic , meaning that its conductivity is the same in all directions , and, further, (c) does not depend upon temperature. Dropping either ( b) or (c) would result in a much more complicated nonlinear diffusion equation. The second thermodynamical principle is that, in the absenc e of external heat sources, heat can only enter any subregion D⊂Ω through its boundary ∂D. (Keep in mind that the plate is insulated above and below.) In other words, the r ate of change of the heat energy in Dis prescribed by the heat flux across its boundary ∂D. Letε(t,x,y) denote the heat energy density at each time and point in the domain, s o that/integraldisplay/integraldisplay Dε(t,x,y)dxdy represents thetotalheatcontained withintheregion Dattimet. Theamount ofadditional heat energy entering Dat a boundary point x∈∂Dis the inward normal component of the heat flux vector: −w·n, wherendenotes the outward unit normal to ∂D. Thus, the total heat flux entering the region Dis given by the flux line integral −/contintegraldisplay ∂Dw·nds, 12/11/12 938 c/ci∇cleco√y∇t2012 Peter J. Olver cf. (A.42). Equating the rate of change of heat energy to the h eat flux yields ∂ ∂t/integraldisplay/integraldisplay Dε(t,x,y)dxdy=−/contintegraldisplay ∂Dw·nds=−/integraldisplay/integraldisplay D∇·wdxdy, where we applied the divergence form of Green’s Theorem, (A. 60), to convert the flux line integral into a double integral. We bring the time derivativ e inside the first integral and collect the terms, whence /integraldisplay/integraldisplay D/parenleftbigg∂ε ∂t+∇·w/parenrightbigg dxdy= 0. (17.7) Keep in mind that this integral formula must hold for anysubdomain D⊂Ω. Now, the only way in which an integral of a continuous function can van ish for all subdomains is if the integrand is identically zero, cf. Exercise , and so ∂ε ∂t+∇·w= 0. (17.8) In this way, we derive the basic conservation law relating heat energy εand heat flux w. As in our one-dimensional model, (14.2), the heat energy ε(t,x,y) at each time and point in the domain is proportional to the temperature, so ε(t,x,y) =σ(x,y)u(t,x,y),where σ(x,y) =ρ(x,y)χ(x,y) (17 .9) is the product of the densityand theheat capacity of the material. Combining this with the Fourier Law (17.6) and the energy balance equation (17.9 ) leads to the general two- dimensional diffusion equation σ∂u ∂t=∇·/parenleftbig κ∇u/parenrightbig (17.10) governing the temperature dynamics of an isotropic medium i n the absence of external heat sources. In full detail, this second order partial diffe rential equation is σ(x,y)∂u ∂t=∂ ∂x/parenleftbigg κ(x,y)∂u ∂x/parenrightbigg +∂ ∂y/parenleftbigg κ(x,y)∂u ∂y/parenrightbigg . (17.11) In particular, if the body is homogeneous, then both σandκare constant, and so general diffusion equation (17.10) reduces to the heat equation (17. 1) with thermal diffusivity γ=κ σ=κ ρχ. (17.12) The heat and diffusion equations are also used to model moveme nts of populations, e.g., bacteria in a petri dish or wolves in the Canadian Rocki es, [biol]. The solution u(t,x,y) represents the number of individuals near position ( x,y) at time t, and diffuses over the domain due to random motions of the individuals. Sim ilar diffusion processes model the mixing of chemical reagents in solutions, with the diffusion induced by the randomBrownianmotionfrommolecularcollisions. Convect ionduetofluidmotionand/or changes due to chemical reactions lead to the more general cl ass ofreaction–diffusion and convective–diffusion equations , [chem]. 12/11/12 939 c/ci∇cleco√y∇t2012 Peter J. Olver Self-Adjoint Formulation Thegeneraldiffusionequation(17.11)canbereadilyfitinto theself-adjoint framework of Section 14.7, taking the form ut=−K[u] =−∇∗◦∇u. (17.13) The gradient operator ∇maps scalar fields uto vector fields v=∇u; its adjoint ∇∗, which goes in the reverse direction, is taken with respect to the weighted inner products /a\}b∇acketle{tu;/tildewideu/a\}b∇acket∇i}ht=/integraldisplay/integraldisplay Ωu(x,y)/tildewideu(x,y)σ(x,y)dxdy,/a\}b∇acketle{t/a\}b∇acketle{tv;/tildewidev/a\}b∇acket∇i}ht/a\}b∇acket∇i}ht=/integraldisplay/integraldisplay Ωv(x,y)·/tildewidev(x,y)κ(x,y)dxdy, (17.14) between, respectively, scalar and vector fields. A straight forward integration by parts argument, done in detail in Section 15.4, tells us that ∇∗v=−1 σ∇·(κv) =−1 σ/bracketleftbigg∂(κv1) ∂x+∂(κv2) ∂y/bracketrightbigg ,when v=/parenleftbigg v1 v2/parenrightbigg .(17.15) Therefore, the right hand side of (17.13) is equal to −K[u] =−∇∗◦∇u=1 σ∇·(κ∇u), (17.16) which recovers (17.10). As always, we need to impose suitabl e homogeneous boundary conditions — Dirichlet, Neumann or mixed — to ensure the vali dity of the integration by parts argument used to establish the adjoint formula (17.15 ). In particular, to obtain the heat equation, we take σandκto be constant, and so (17.14) reduce, up to a constant factor, to the usual L2inner products between scalar and vector fields. In this case, the adjoint of the gradient is , up to a scale factor, minus the divergence: ∇∗=−γ∇·, whereγ=κ/σ. In this scenario, (17.13) reduces to the two-dimensional heat equation (17.1). The self-adjoint operator K=∇∗◦∇is always positive semi-definite, and is positive definite if and only if ker ∇={0}. As we saw in Section 15.4, positive definiteness holds in the Dirichlet and mixed cases, whereas the Neumann bounda ry conditions lead to only a positive semi-definite operator. Indeed, assuming Ω is con nected, the only functions in the kernel of the gradient operator are the constants. Moreo ver, only the zero constant function satisfies the Dirichlet or mixed boundary conditio ns, and hence the gradient’s kernel is trivial, ensuring positive definiteness. The heat and diffusion equations are examples of parabolic partial differential equa- tions, the terminology being adapted from Definition 15.1 to apply to partial differential equations in more than two variables. As we will see, all of th e basic qualitative features we learned when studying the one-dimensional heat equation carry over to solutions to parabolic partial differential equations in higher dimensi ons. Separation of Variables In Section 14.1, we applied the method of separation of varia bles to express the so- lution to the one-dimensional heat equation as a series invo lving the eigenfunctions of an 12/11/12 940 c/ci∇cleco√y∇t2012 Peter J. Olver associated boundary value problem. Section 14.7 argued tha t the eigenfunction series solu- tionmethod carries over, essentially unchanged, to thegen eral class of self-adjoint diffusion equations. In particular, the solutions to the two-dimensi onal heat and diffusion equations are expressed in series form based on the eigenfunctions of a n associated two-dimensional boundary value problem. As we know, the separable solutions to any diffusion equation (17.13) are of exponen- tial form u(t,x,y) =e−λtv(x,y). (17.17) Sincethelinearoperator Konlyinvolvesdifferentiationwithrespecttothespatialva riables x,y, we find ∂u ∂t=−λe−λtv(x,y),while K[u] =e−λtK[v]. Substituting back into the diffusion equation (17.13) and ca nceling the exponentials, we conclude that K[v] =λv, (17.18) Thus,v(x,y) must be an eigenfunction for the linear operator K, subject to the relevant homogeneous boundary conditions. In the case of the heat equ ation (17.1), K[u] =−γ∆u, and hence the eigenvalue equation (17.18) takes the form γ∆v+λv= 0,or, in detail, γ/parenleftbigg∂2v ∂x2+∂2v ∂y2/parenrightbigg +λv= 0.(17.19) This generalization of the Laplace equation is known as the Helmholtz equation , and was briefly discussed in Example . According to Theorem 14.18, the eigenvalues of the self-adj oint operator K=∇∗◦∇ are all real and non-negative: λ≥0. When Kispositive definite, they are strictly positive: λ >0. Let us index the eigenvalues in increasing order: 0< λ1≤λ2≤λ3≤ ···. (17.20) Under reasonable conditions on the underlying linear bound ary value problem, it can be proved, [ 47;Chapter V], that there are, in fact, an infinite number of eig envalues, and, moreover, their size increases without bound, so λk→ ∞ask→ ∞. Each eigenvalue is repeated according to the number (which is necessarily fin ite) of linearly independent eigenfunctions it admits. The problem has a zero eigenvalue ,λ0= 0, corresponding to the constant null eigenfunction v(x,y)≡1, if and only if Kis not positive definite, i.e., only in the case of pure Neumann boundary conditions. Each eigenvalue and eigenfunction pair will produce a separ able solution uk(t,x,y) =e−λktvk(x,y) to the diffusion equation (17.13). The solutions correspond ing to positive eigenvalues are exponentially decaying in time, while a zero eigenvalue, wh ich only occurs in the Neumann 12/11/12 941 c/ci∇cleco√y∇t2012 Peter J. Olver problem, produces a constant solution. The general solutio nto the homogeneous boundary value problem can then be built up as a linear combination of t hese basic solutions, in the form of an eigenfunction series u(t,x,y) =∞/summationdisplay k=1ckuk(t,x,y) =∞/summationdisplay k=1cke−λktvk(x,y), (17.21) which is a form of generalized Fourier series. The eigenfunc tion coefficients ckare pre- scribed by the initial conditions, which require ∞/summationdisplay k=1ckvk(x,y) =f(x,y). (17.22) Theorem 14.19 guarantees orthogonality of the eigenfuncti ons, and hence the coefficients are given by our standard orthogonality formula ck=/a\}b∇acketle{tf;vk/a\}b∇acket∇i}ht /ba∇dblvk/ba∇dbl2=/integraldisplay/integraldisplay Ωf(x,y)vk(x,y)σ(x,y)dxdy /integraldisplay/integraldisplay Ωvk(x,y)2σ(x,y)dxdy, (17.23) where the weighting function σ(x,y)wasdefined in (17.9). Inthe caseof theheat equation, σis constant and so can be canceled from both numerator and den ominator, leaving the simpler formula ck=/integraldisplay/integraldisplay Ωf(x,y)vk(x,y)dxdy /integraldisplay/integraldisplay Ωvk(x,y)2dxdy. (17.24) Under fairly general hypotheses, it can be shown that the eig enfunctions form a com- plete system , which means that the eigenfunction series (17.22) will con verge (at least in norm) to the function f(x,y), provided it is not too bizarre. Moreover, the resulting se ries (17.21) converges to the solution to the initial-boundary v alue problem for the diffusion equation. See [ 47;p. 369] for a precise statement and proof of the general theo rem. Qualitative Properties Before tackling simple examples where we are able to constru ct explicit formulae for the eigenfunctions and eigenvalues, let us see what the eige nfunction series solution (17.21) can tell us about general diffusion processes. Based on our ex perience with the case of a one-dimensional bar, the final conclusions will not be espe cially surprising. Indeed, they also apply, word for word, to diffusion processes in thre e-dimensional solid bodies. A reader who prefers to see explicit solution formulae may wis h to skip ahead to the following section, returning here later. Keep in mind that we are still dealing with the solution to the homogeneous boundary value problem. The first observation is that all terms in the s eries solution (17.21), with the possible exception of a null eigenfunction term that app ears in the semi-definite case, 12/11/12 942 c/ci∇cleco√y∇t2012 Peter J. Olver are tending to zero exponentially fast. Since most eigenval ues are large, all the higher order terms in the series become almost instantaneously neg ligible, and hence the solution can be accurately approximated by a finite sum over the first fe w eigenfunction modes. As time goes on, more and more of the modes can be neglected, an d the solution decays to thermal equilibrium at an exponentially fast rate. The ra te of convergence to thermal equilibrium is, for most initial data, governed by the small est positive eigenvalue λ1>0 for the Helmholtz boundary value problem on the domain. In the positive definite cases of homogeneous Dirichlet or mi xed boundary conditions, thermal equilibrium is u(t,x,y)→u⋆(x,y)≡0. Thus, in these cases, the equilibrium temperature is equal to the boundary temperature — even if th is temperature is only fixed on a small part of the boundary. The initial heat is eventuall y dissipated away through the non-insulated part of the boundary. In the semi-definite Neumann case, corresponding to a completely insulated plate, In this case, the general so lution has the form u(t,x,y) =c0+∞/summationdisplay k=1cke−λktvk(x,y), (17.25) where the sum is over the positive eigenmodes, λk>0. Since all the summands are exponentially decaying, the final equilibrium temperature u⋆=c0is the same as the con- stant term in the eigenfunction expansion. We evaluate this term using the orthogonality formula (17.23), and so, as t→ ∞, u(t,x,y)−→c0=/a\}b∇acketle{tf;1/a\}b∇acket∇i}ht /ba∇dbl1/ba∇dbl2=/integraldisplay/integraldisplay Ωf(x,y)σ(x,y)dxdy /integraldisplay/integraldisplay Ωσ(x,y)dxdy, representing a suitably weighted average of the initial tem perature over the domain. In particular, in the case of the heat equation, the weighting f unctionσis constant, and so the equilibrium temperature u(t,x,y)−→c0=1 area Ω/integraldisplay/integraldisplay Ωf(x,y)dxdy (17.26) equals the average initial temperature distribution. In th is case, the heat energy is not allowed to escape through the boundary, and thus redistribu tes itself in a uniform manner over the domain. Diffusion has a smoothing effect on the initial temperature di stribution f(x,y). As- sume that the eigenfunction coefficients are uniformly bound ed, so|ck| ≤Mfor some constant M. This will certainly be the case if f(x,y) is piecewise continuous, but even holds for quite rough initial data, including delta functio ns. Then, at any time t >0 after the initial instant, the coefficients cke−λktin the eigenfunction series solution (17.21) are exponentially small as k→ ∞, which is enough to ensure smoothness†of the solution u(t,x,y) for each t >0. Therefore, the diffusion process serves to immediately sm ooth out †Forageneraldiffusionequation, thisrequiresthatthefunctions σ(x,y)andκ(x,y)besmooth. 12/11/12 943 c/ci∇cleco√y∇t2012 Peter J. Olver Figure 17.1. Smoothing a Grey Scale Image. jumps, corners and other discontinuities in the initial dat a. As time progresses, the local variations in the solution become less and less pronounced, as it asymptotically reaches a constant equilibrium state. As a result, diffusion processes can be effectively applied to clean and denoise planar images. The initial data f(x,y) represents the grey-scale value of the image at position (x,y), so that 0 ≤f(x,y)≤1 with 0 representing black, and 1 representing white. As time progresses, the solution u(t,x,y) represents a more and more smoothed version of the image. Although this has the effect of removing unwanted nois e from the image, there is also a gradual blurring of the actual features. Thus, the “ti me” or “multiscale” parameter tneeds to be chosen to optimally balance between the two effect s — the larger tis the more noise is removed, but the more noticeable the blurring. A representative illustration appears in Figure 17.1. To further suppress undesirable blu rring effects, recent image processing filters are based on anisotropic (and thus nonlinear ) diffusion equations. See Sapiro, [162], for a survey of recent progress in this active field. Since the forward heat equation effectively blurs the featur es in an image, we might be tempted to try going backwards in time to sharpen the image . However, the argu- ment presented in Section 14.1 tells us that the backwards he at equation is ill-posed, and hence cannot be used directly for this purpose. Various “reg ularization” strategies have been devised to circumvent this mathematical barrier, and t hereby design effective image enhancement algorithms, [ reg]. Inhomogeneous Boundary Conditions and Forcing Finally, let us briefly mention how to incorporate inhomogen eous boundary conditions and external heat sources into the problem. Consider, as a sp ecific example, the forced heat equation ut=γ∆u+F(x,y),for ( x,y)∈Ω, (17.27) whereF(x,y) represents an unvarying external heat source, subject to i nhomogeneous Dirichlet boundary conditions u=hfor ( x,y)∈∂Ω, (17.28) that fixes the temperature of the plate on its boundary. When t he external forcing is fixed for allt, we expect the solution to eventually settle down to an equil ibrium configuration: u(t,x,y)→u⋆(x,y) ast→ ∞, which will be justified below. 12/11/12 944 c/ci∇cleco√y∇t2012 Peter J. Olver The time-independent equilibrium temperature u=u⋆(x,y) satisfies the equation obtained by setting ut= 0 in the differential equation (17.27), namely Poisson equa tion −γ∆u⋆=F, for ( x,y)∈Ω, (17.29) and subject to the same inhomogeneous Dirichlet boundary co nditions (17.28). Positive definiteness of the Dirichlet boundary value problem implie s that there is a unique equi- librium solution. Once we have determined the equilibrium solution — usually t hrough a numerical approximation — we set v(t,x,y) =u(t,x,y)−u⋆(x,y), so thatvmeasures the deviation of the solution ufrom its eventual equilibrium. By linearity v(t,x,y) satisfies the unforced heat equation subject to homogeneou s boundary conditions: vt=γ∆v,(x,y)∈Ω, u = 0,(x,y)∈∂Ω. (17.30) Therefore, vcan be expanded in an eigenfunction series (17.21), and will decay to zero, v(t,x,y)→0, at a exponentially fast rate prescribed by the smallest ei genvalue λ1of the associated homogeneous Helmholtz boundary value problem. Consequently, the solution to the forced, inhomogeneous problem u(t,x,y) =v(t,x,y)+u⋆(x,y)−→u⋆(x,y) will approach thermal equilibrium, namely u⋆(x,y), at exactly the same exponential rate as its homogeneous counterpart. 17.2. Explicit Solutions for the Heat Equation. Thus, solving the two-dimensional heat equation in series f orm requires knowing the eigenfunctions for the associated Helmholtz boundary valu e problem. Unfortunately, as with any partial differential equation, explicit solution f ormulae are few and far between. In thissection, wediscuss two specific cases where the requi red eigenfunctions can be found in closed form. The calculations rely on separation of varia bles, which, as we know, works in only a very limited class of domains. Nevertheless, inter esting solution features can be gleaned from these particular cases. Heating of a Rectangle A homogeneous rectangular plate R=/braceleftbig 0< x < a, 0< y < b/bracerightbig is heated to a prescribed initial temperature u(0,x,y) =f(x,y),for ( x,y)∈R, (17.31) and then insulated. The sides of the plate are held at zero tem perature. Our task is to determine how fast the plate returns to thermal equilibrium . 12/11/12 945 c/ci∇cleco√y∇t2012 Peter J. Olver The temperature u(t,x,y) evolves according to the two-dimensional heat equation ut=γ(uxx+uyy),for ( x,y)∈R, t > 0, (17.32) subject to homogeneous Dirichlet conditions u(0,y) =u(a,y) = 0 =u(x,0) =u(x,b),0< x < a, 0< y < b, (17.33) along the boundary of the rectangle. As in (17.17), the eigen solutions to the heat equation are obtained from the usual exponential ansatz u(t,x,y) =e−λtv(x,y). Substituting this expression into the heat equation, we conclude that the function v(x,y) solves the Helmholtz eigenvalue problem γ(vxx+vyy)+λv= 0,(x,y)∈R, (17.34) subject to the same homogeneous Dirichlet boundary conditi ons v(x,y) = 0,(x,y)∈∂R. (17.35) To solve the rectangular Helmholtz eigenvalue problem (17. 34–35), we shall, as in (15.13), introduce a further separation of variables, writing v(x,y) =p(x)q(y) as the product of functions depending upon the individual Ca rtesian coordinates. Substi- tuting this ansatz into the Helmholtz equation (17.34), we fi nd γp′′(x)q(y)+γp(x)q′′(y)+λp(x)q(y) = 0. To effect the variable separation, we collect all terms invol vingxon one side and all terms involving yon the other side of the equation. This is accomplished by div iding by v=pq and rearranging the terms; the result is γp′′(x) p(x)=−γq′′(y) q(y)−λ≡ −µ. The left hand side of this equation only depends on x, whereas the right hand side only dependson y. AsarguedinSection15.2,theonlywaythiscanoccurisifth etwosidesequal a common separation constant , denoted by −µ. (The minus sign is for later convenience.) In this manner, we reduce our partial differential equation t o a pair of one-dimensional eigenvalue problems γp′′+µp= 0, γq′′+(λ−µ)q= 0, each of which is subject to homogeneous Dirichlet boundary c onditions p(0) =p(a) = 0, q (0) =q(b) = 0, stemming from the boundary conditions (17.35). To obtain a n ontrivial solution to the Helmholtz equation, we seek nonzero solutions to these two s upplementary eigenvalue problems. The fact that we are dealing with a rectangular dom ain is critical to the success of this procedure. 12/11/12 946 c/ci∇cleco√y∇t2012 Peter J. Olver We have already solved these particular two boundary value p roblems many times; see, for instance, (14.17). The eigenfunctions are, respec tively, pm(x) = sinmπx a, m= 1,2,3,..., qn(y) = sinnπy b, n= 1,2,3,..., with µ=m2π2γ a2, λ−µ=n2π2γ b2,so that λ=m2π2γ a2+n2π2γ b2. Therefore, the separable eigenfunction solutions to the He lmholtz boundary value problem (17.33–34) have the doubly trigonometric form vm,n(x,y) = sinmπx asinnπy b, (17.36) with associated eigenvalues λm,n=m2π2γ a2+n2π2γ b2=/parenleftbiggm2 a2+n2 b2/parenrightbigg π2γ. (17.37) Each of these corresponds to an exponentially decaying, sep arable solution um,n(t,x,y) =e−λm,ntvm,n(x,y) = exp/bracketleftbigg −/parenleftbiggm2 a2+n2 b2/parenrightbigg π2γt/bracketrightbigg sinmπx asinnπy b (17.38) to the original rectangular boundary value problem for the h eat equation. Using the fact that the univariate sine functions form a comp lete system, it is not hard to prove, [ 187], that the separable eigenfunction solutions (17.38) are c omplete, which implies that there are no non-separable eigenfunctions. As a consequence, the general solution to the initial-boundary value problem can be expre ssed as a linear combination u(t,x,y) =∞/summationdisplay m,n=1cm,num,n(t,x,y) =∞/summationdisplay m,n=1cm,ne−λm,ntvm,n(x,y) (17 .39) of our eigenfunction modes. The coefficients cm,nare prescribed by the initial conditions, which take the form of a double Fourier sine series f(x,y) =u(0,x,y) =∞/summationdisplay m,n=1cm,nvm,n(x,y) =∞/summationdisplay m,n=1cm,nsinmπx asinnπy b. Self-adjointness of the Laplacian coupled with the boundar y conditions implies that the eigenfunctions vm,n(x,y) are orthogonal with respect to the L2inner product on the rectangle: /a\}b∇acketle{tvk,l;vm,n/a\}b∇acket∇i}ht=/integraldisplayb 0/integraldisplaya 0vk,l(x,y)vm,n(x,y)dxdy= 0 unless k=mandl=n. 12/11/12 947 c/ci∇cleco√y∇t2012 Peter J. Olver Figure 17.2. Heat Diffusion in a Rectangle. (The skeptical reader can verify the orthogonality relatio ns directly from the eigenfunction formulae (17.36).) Thus, we can appeal to our usual orthogon ality formula (17.24) to evaluate the coefficients cm,n=/a\}b∇acketle{tf;vm,n/a\}b∇acket∇i}ht /ba∇dblvm,n/ba∇dbl2=4 ab/integraldisplayb 0/integraldisplaya 0f(x,y) sinmπx asinnπy bdxdy, (17.40) where the formula for the norms of the eigenfunctions /ba∇dblvm,n/ba∇dbl2=/integraldisplayb 0/integraldisplaya 0vm,n(x,y)2dxdy=/integraldisplayb 0/integraldisplaya 0sin2mπx asin2nπy bdxdy=1 4ab.(17.41) follows from a direct evaluation of the double integral. (Un fortunately, while the orthogo- nality is automatic, the computation of the norm must inevit ably be done “by hand”.) The rectangle approaches thermal equilibrium at the rate eq ual to the smallest eigen- value: λ1,1=/parenleftbigg1 a2+1 b2/parenrightbigg π2γ, (17.42) i.e., the sum of the reciprocals of the squared lengths of its sides multiplied by the diffusion coefficient. The larger the rectangle, or the smaller the diffu sion coefficient, the smaller λ1,1, and hence slower the return to thermal equilibrium. The exp onentially fast decay rate of the Fourier series implies that the solution immedia tely smooths out any initial discontinuites in the initial temperature profile. Indeed, the higher modes, with mandn large, decay to zero almost instantaneously, and so the solu tion immediately behaves like a finite sum over a few low order modes. Assuming that c1,1/\e}atio\slash= 0, the slowest decaying mode in the Fourier series (17.39) is c1,1u1,1(t,x,y) =c1,1exp/bracketleftbigg −/parenleftbigg1 a2+1 b2/parenrightbigg π2γt/bracketrightbigg sinπx asinπy b. (17.43) Thus, in the long run, the temperature is of one sign — either p ositive or negative de- pending upon the sign of c1,1— throughout the rectangle. This observation is, in fact, 12/11/12 948 c/ci∇cleco√y∇t2012 Peter J. Olver indicative of the general phenomenon that the eigenfunctio n associated with the smallest positive eigenvalue of a self-adjoint elliptic operator is of one sign throughout the domain. A typical solution is plotted in Figure 17.2. Heating of a Disk Let us perform a similar analysis of the thermodynamics of a c ircular disk. For simplicity (or by choice of suitable physical units), we wil l assume that the disk D={x2+y2≤1} ⊂R2 has unit radius and unit diffusion coefficient γ= 1. We shall solve the heat equation on Dsubject to homogeneous Dirichlet boundary values of zero te mperature at the circular edge ∂D=C={x2+y2= 1}. Thus, the full initial-boundary value problem is ∂u ∂t=∂2u ∂x2+∂2u ∂y2, x2+y2<1, t > 0, u(t,x,y) = 0, x2+y2= 1, u(0,x,y)=f(x,y), x2+y2≤1.(17.44) We remark that a simple rescaling of space and time, as outlin ed in Exercise , can be used to recover the solution for an arbitrary diffusion coeffic ient and a disk of arbitrary radius from this particular case. Since we are working in a circular domain, we instinctively p ass to polar coordinates (r,θ). In view of the polar coordinate formula (15.29) for the Lap lace operator, the heat equation and boundary conditions assume the form ∂u ∂t=∂2u ∂r2+1 r∂u ∂r+1 r2∂2u ∂θ2, 0≤r <1, t > 0, u(t,1,θ) = 0, u (0,r,θ) =f(r,θ), r ≤1,(17.45) where the solution u(t,r,θ) is defined for all 0 ≤r≤1 andt≥0. To ensure that the solution represents a single-valued function on the entire disk, it is required to be a 2 π periodic function of the angular variable: u(t,r,θ+2π) =u(t,r,θ) To obtain the separable solutions u(t,r,θ) =e−λtv(r,θ), (17.46) we need to solve the polar coordinate form of the Helmholtz eq uation ∂2v ∂r2+1 r∂v ∂r+1 r2∂2v ∂θ2+λv= 0,0≤r <1, 0≤θ≤2π,(17.47) 12/11/12 949 c/ci∇cleco√y∇t2012 Peter J. Olver subject to the boundary conditions v(1,θ) = 0, v (r,θ+2π) =v(r,θ). (17.48) We invoke a further separation of variables by writing v(r,θ) =p(r)q(θ). (17.49) Substituting this ansatz into the polar Helmholtz equation (17.47), and then collecting together all terms involving rand all terms involving θ, we are led to the pair of ordinary differential equations r2p′′+rp′+(λr2−µ)p= 0, q′′+µq= 0, (17.50) whereλis the Helmholtz eigenvalue, and µthe separation constant. Let us start with the equation for q(θ). The periodicity condition (17.48) requires that q(θ) be 2πperiodic. Therefore, the required solutions are the elemen tary trigonometric functions q(θ) = cosmθor sin mθ, where µ=m2, (17.51) withm= 0,1,2,...a non-negative integer. Substitutingtheformula, µ=m2,fortheseparationconstant, thedifferentialequation forp(r) takes the form r2d2p dr2+rdp dr+(λr2−m2)p= 0,0≤r≤1. (17.52) Ordinarily, one imposes two boundary conditions in order to pin down a solution to such a second order ordinary differential equation. But our Dirich let condition, namely p(1) = 0, only specifies its value at one of the endpoints. The other end point is a singular point for the ordinary differential equation, because the coefficient o f the highest order derivative, namelyr2, vanishes at r= 0. This situation might remind you of our solution to the Eul er differential equation (15.35) in the context of separable so lutions to the Laplace equation on the disk. As there, we only require the solution to be bound ed atr= 0: |p(0)|<∞, p (1) = 0. (17.53) As we now show, this pair of boundary conditions suffices to dis tinguish the relevant eigenfunction solutions to (17.52). Althoughtheordinarydifferentialequation(17.52)appear sinavarietyofapplications, this may be the first time that you have encountered it. It is no t an elementary equation; indeed, most solutions cannot be written in terms of the elem entary functions you see in first year calculus. Nevertheless, owing to their significan ce in a wide range of physical applications, the solutions have been extensively studied and tabulated, and so are, in a sense, well-known. After some preparatory manipulations , we shall summarize their relevant properties, relegating details and proofs to Appe ndix C. To simplify the analysis, we make a preliminary rescaling of the independent variable, replacing rby z=√ λr. 12/11/12 950 c/ci∇cleco√y∇t2012 Peter J. Olver 5 10 15 20 -0.4-0.20.20.40.60.81 J0(z)5 10 15 20 -0.4-0.20.20.40.60.81 J1(z)5 10 15 20 -0.4-0.20.20.40.60.81 J2(z) Figure 17.3. Bessel Functions. Note that, by the chain rule, dp dr=√ λdp dz,d2p dr2=λd2p dz2, and hence rdp dr=zdp dz, r2d2p dr2=z2d2p dz2. The net effect is to eliminate the eigenvalue parameter λ(or, rather, hide it in the change of variables), so that (17.52) assumes the slightly simpler form z2d2p dz2+zdp dz+(z2−m2)p= 0. (17.54) The ordinary differential equation (17.54) is known as Bessel’s equation , named after the early 19thcentury astronomer Wilhelm Bessel, who first used its soluti ons to analyze planetary orbits. The solutions to Bessel’s equation have b ecome an indispensable tool in applied mathematics, physics, and engineering. To begin, the one thing we know for sure is that, as with any sec ond order ordinary differential equation, there are two linearly independent s olutions. However, it turns out that, up to a constant multiple, only one solution remains bo unded as z→0. This solution is known as the Bessel function oforderm, and is denoted by Jm(z). Applying the general systematic method for finding power series solution s to linear ordinary differential equations presented in Appendix C, it can be shown that the Be ssel function of order m has the Taylor expansion Jm(z) =∞/summationdisplay k=0(−1)kzm+2k 2m+2kk!(m+k)!(17.55) =zm 2mm!/bracketleftbigg 1−z2 4(m+1)+z4 32(m+1)(m+2)−z6 384(m+1)(m+2)(m+3)+···/bracketrightbigg . A simple application of the ratio test tells us that the power series converges for all (com- plex) values of z, and hence Jm(z) is everywhere analytic. Indeed, the convergence is quite rapid when zis of moderate size, and so summing the series is a reasonably effective method for computing the Bessel function Jm(z) — although in serious applications one adopts more sophisticated numerical techniques based on as ymptotic expansions and in- tegral formulae, [ 3,144]. Moreover, verification that it indeed solves the Bessel eq uation (17.54) is a straightforward computation. Figure 17.3 disp lays graphs of the first three 12/11/12 951 c/ci∇cleco√y∇t2012 Peter J. Olver Bessel functions for z >0. Most software packages, both symbolic and numerical, inc lude routines for accurately evaluating and graphing Bessel fun ctions. Reverting back to our original radial coordinate r=z/√ λ, we conclude that every solution to the radial equation (17.52) which is bounded at r= 0 is a constant multiple p(r) =Jm/parenleftbig√ λr/parenrightbig (17.56) of the rescaled Bessel function of order m. So far, we have only dealt with the boundary condition at the singular point r= 0. The Dirichlet condition at the other end requires p(1) =Jm/parenleftbig√ λ/parenrightbig = 0. Therefore, in order that λbe a legitimate eigenvalue,√ λmust be a rootof themthorder Bessel function Jm. Remark: We already know, thanks to the positive definiteness of the D irichlet bound- ary value problem, that the Helmholtz eigenvalues λ >0 must be positive, and so there is no problem taking the square root. Indeed, it can be proved , [186], that the Bessel functions do not have any negative roots! The graphs of Jm(z) strongly indicate, and, indeed, it can be rigorously prove d, that each Bessel function oscillates between positive and negat ive values as zincreases above 0, with slowly decreasing amplitude. In fact, it can be prove d that, asymptotically, Jm(z)∼/radicalbigg 2 πzcos/parenleftbig z−/parenleftbig1 2m+1 4/parenrightbig π/parenrightbig as z−→ ∞, (17.57) and so the oscillations become essentially the same as a (pha se-shifted) cosine whose am- plitude decreases, relatively slowly, like z−1/2. As a consequence, there exists an infinite sequence of Bessel roots , which we number in the order in which they appear: Jm(ζm,n) = 0,where 0< ζm,1< ζm,2< ζm,3<···withζm,n−→ ∞ asn−→ ∞.(17.58) ItisworthnotingthattheBessel functions are notperiodic, andtheir rootsarenotinitially evenly spaced. Owing to their physical importance in a wide range of problem s, the Bessel roots have been extensively tabulated in the literature, cf. [ 3,145]. The accompanying table displays all Bessel roots that are <12 in magnitude. The columns of the table are indexed by m, the order of the Bessel function, and the rows by n, the root number. 12/11/12 952 c/ci∇cleco√y∇t2012 Peter J. Olver Table of Bessel Roots ζm,n n/backslashbigg m 0 1 2 3 4 5 6 7 ... 1 2.4048 3 .8317 5 .1356 6.3802 7.5883 8 .7715 9.9361 11 .0860... 2 5.5201 7 .0156 8 .4172 9 .761 11 .0650......... 3 8.6537 10 .1730 11 .6200...... 4 11.7920...... ...... Remark: According to (17.55), Jm(0) = 0 for m >0,while J0(0) = 1. However, we do not count 0 as a bona fide Bessel root, since it does not lead to a valid eigenfunction for the Helmholtz boundary value problem. Summarizing our progress, the eigenvalues λm,n=ζ2 m,n, n = 1,2,3,..., m = 0,1,2,..., (17.59) of the Bessel boundary value problem (17.52–53) are the squa res of the roots of the Bessel function of order m. The corresponding eigenfunctions are wm,n(r) =Jm(ζm,nr), n = 1,2,3,..., m = 0,1,2,..., (17.60) defined for 0 ≤r≤1. Combining (17.60) with the formula (17.51) for the angula r compo- nents, we conclude that the separable solutions (17.49) to t he polar Helmholtz boundary value problem (17.47) are v0,n(r,θ) =J0(ζ0,nr), vm,n(r,θ) =Jm(ζm,nr) cosmθ, /hatwidevm,n(r,θ) =Jm(ζm,nr) sinmθ,wheren= 1,2,3,..., m= 0,1,2,....(17.61) These solutions define the so-called normal modes for the unit disk, and Figure 17.4 plots the first few of them. The eigenvalues λ0,nare simple, and contribute radially symmetric eigenfunctions, whereas the eigenvalues λm,nform >0 are double, and produce two lin- early independent separable eigenfunctions, with trigono metric dependence on the angular variable. We have at last produced the basic separable solutions u0,n(t,r) =e−ζ2 0,ntJ0(ζ0,nr), um,n(t,r,θ) =e−ζ2 m,ntJm(ζm,nr) cosmθ, /hatwideum,n(t,r,θ) =e−ζ2 m,ntJm(ζm,nr) sinmθ,n= 1,2,3,..., m= 1,2,....(17.62) 12/11/12 953 c/ci∇cleco√y∇t2012 Peter J. Olver v0,1(r,θ) v0,2(r,θ) v0,3(r,θ) v1,1(r,θ) v1,2(r,θ) v1,3(r,θ) v2,1(r,θ) v2,2(r,θ) v2,3(r,θ) v3,1(r,θ) v3,2(r,θ) v3,3(r,θ) Figure 17.4. Normal Modes for a Disk. to the homogeneous Dirichlet boundary value problem for the heat equation on the unit disk (17.45). The general solution is a linear superpositio n, in the form of an infinite series u(t,r,θ) =1 2∞/summationdisplay n=1a0,nu0,n(t,r)+∞/summationdisplay m,n=1/bracketleftbig am,num,n(t,r,θ)+bm,n/hatwideum,n(t,r,θ)/bracketrightbig ,(17.63) where the initial factor of1 2is included, as with ordinary Fourier series, for later con- venience. As usual, the coefficients am,n,bm,nare determined by the initial condition, 12/11/12 954 c/ci∇cleco√y∇t2012 Peter J. Olver Norms of the Fourier–Bessel Eigenfunctions /ba∇dblvm,n/ba∇dbl=/ba∇dbl/hatwidevm,n/ba∇dbl n/backslashbigg m0 1 2 3 4 5 6 7 1 .9202.5048.4257.3738.3363.3076.2847.2658 2 .6031.3761.3401.3126.2906.2725.2572.2441 3 .4811.3130.2913.2736.2586.2458.2347.2249 4 .4120.2737.2589.2462.2352.2255.2169.2092 5 .3661.2462.2353.2257.2171.2095.2025.1962 so u(0,r,θ) =1 2∞/summationdisplay n=1a0,nv0,n(r)+∞/summationdisplay m,n=1/bracketleftbig am,nvm,n(r,θ)+bm,n/hatwidevm,n(r,θ)/bracketrightbig =f(r,θ). (17.64) Thus, we must expand the initial data into a Fourier–Bessel series in the eigenfunctions. As in the rectangular case, it is possible to prove, [ 47], that the separable eigenfunctions arecomplete — there are no other eigenfunctions — and, moreover, every (r easonable) function defined on the unit disk can be written as a convergen t series in the Bessel eigenfunctions. Theorem 14.19 gurantees that the eigenfunctions are orthog onal†with respect to the standard L2inner product /a\}b∇acketle{tu;v/a\}b∇acket∇i}ht=/integraldisplay/integraldisplay Du(x,y)v(x,y)dxdy=/integraldisplay1 0/integraldisplay2π 0u(r,θ)v(r,θ)rdθdr on the unit disk. (Note the extra factor of rcoming from the polar coordinate form (A.51) of the infinitesimal element of area dxdy=rdrdθ.) The L2norms of the Fourier–Bessel functions are given by the interesting formulae /ba∇dblv0,n/ba∇dbl=√π/vextendsingle/vextendsingleJ1(ζ0,n)/vextendsingle/vextendsingle,/ba∇dblvm,n/ba∇dbl=/ba∇dbl/hatwidevm,n/ba∇dbl=/radicalbiggπ 2/vextendsingle/vextendsingleJm+1(ζm,n)/vextendsingle/vextendsingle,(17.65) that involves the value of the Bessel function of the next hig her order at the appropriate Bessel root; numerical values appear in the accompanying ta ble. A proof of this formula can be found in C.3.65–66. Orthogonality of the eigenfunctions implies that the coeffic ients in the Fouier–Bessel series †For the two eigenfunctions corresponding to one of the double eigenvalu es, orthogonality must be verified by hand. 12/11/12 955 c/ci∇cleco√y∇t2012 Peter J. Olver Figure 17.5. Heat Diffusion in a Disk. (17.64) are given by a0,n= 2/a\}b∇acketle{tf;v0,n/a\}b∇acket∇i}ht /ba∇dblv0,n/ba∇dbl2=2 π J1(ζ0,n)2/integraldisplay1 0/integraldisplayπ −πf(r,θ)J0(ζ0,nr)rdθdr, am,n=/a\}b∇acketle{tf;vm,n/a\}b∇acket∇i}ht /ba∇dblvm,n/ba∇dbl2=2 π Jm+1(ζm,n)2/integraldisplay1 0/integraldisplayπ −πf(r,θ)Jm(ζm,nr)rcosmθdθdr, bm,n=/a\}b∇acketle{tf;/hatwidevm,n/a\}b∇acket∇i}ht /ba∇dbl/hatwidevm,n/ba∇dbl2=2 π Jm+1(ζm,n)2/integraldisplay1 0/integraldisplayπ −πf(r,θ)Jm(ζm,nr)rsinmθdθdr.(17.66) In accordance with the general theory, each individual sepa rable solution (17.62) to the heat equation decays exponentially fast, at a rate λm,n=ζ2 m,nprescribed by the square of the corresponding Bessel root. In particular, the domina nt mode, meaning the one that persists the longest, is u0,1(t,r,θ) =e−ζ2 0,1tJ0(ζ0,1r). (17.67) Its decay rate ζ2 0,1≈5.783 (17 .68) is the square of the smallest non-zero root of the Bessel func tionJ0(z). The dominant eigenfunction v0,1(r,θ) =J0(ζ0,1r)>0 is strictly positive within the entire disk and ra- dially symmetric. Consequently, for most initial conditio ns (specifically those for which a0,1/\e}atio\slash= 0), the disk’s temperature distribution eventually becom es entirely of one sign and radially symmetric, while decaying exponentially fast to zero at the rate given by (17.68). See Figure 17.5 for a plot of a typical solution, dis played as successive times t= 0,.04,.08,.12,.16,.2. Note how, in accordance with the theory, the solution almo st immediately acquires a radial symmetry, followed by a fairl y rapid decay to thermal equi- librium. 12/11/12 956 c/ci∇cleco√y∇t2012 Peter J. Olver 17.3. The Fundamental Solution. As we learned in Section 14.1, the fundamental solution to the heat equation measures the temperature distribution resulting from a concentrate d initial heat source, e.g., a hot soldering iron applied instantaneously at one pointe. The p hysical problem is modeled mathematically by imposing a delta function as the initial c ondition for the heat equation, along with the chosen homogeneous boundary conditions. Kno wledge of the fundamental solution enables one to use linear superposition to recover the solution for any other initial data. As in our one-dimensional analysis, we shall concentrate on the most tractable case, when the domain is the entire plane: Ω = R2. Our first goal will be to solve the initial value problem ut=γ∆u, u (0,x,y) =δ(x−ξ,y−η) =δ(x−ξ)δ(y−η), (17.69) fort >0 and (x,y)∈R2. The initial data is a delta function representing a concent rated unit heat source placed at position( ξ,η). The resulting solution u=F(t,x,y;ξ,η)is called thefundamental solution for the heat equation on R2. The quickest route to the desired solution relies on the foll owing simple lemma that combines solutions of the one-dimensional heat equation to produce solutions of the two- dimensional version. Lemma 17.1. Ifv(t,x)andw(t,x)are any two solutions to the one-dimensional heat equation ut=γuxx, then the product u(t,x,y) =v(t,x)w(t,y) (17 .70) is a solution to the two-dimensional heat equation ut=γ(uxx+uyy). Proof: Our assumptions imply that that vt=γvxx, whilewt=γwyywhen we write w(t,y) as a function of tandy. Therefore, differentiating (17.70), we find ∂u ∂t=∂v ∂tw+v∂w ∂t=γ∂2v ∂x2w+γv∂2w ∂y2=γ/parenleftbigg∂2u ∂x2+∂2u ∂y2/parenrightbigg , and hence u(t,x,y) solves the heat equation. Q.E.D. For example, if v(t,x) =e−γω2tsinωx, w (t,y) =e−γν2tsinνy, are separable solutions of the one-dimensional heat equati on, then u(t,x,y) =e−γ(ω2+ν2)tsinωxsinνy are the separable solutions we used to solve the heat equatio n on a rectangle. A more interesting case is to chose v(t,x) =1 2√πγte−(x−ξ)2/(4γt), w (t,y) =1 2√πγte−(y−η)2/(4γt),(17.71) to be the fundamental solutions (14.59) to the one-dimensio nal heat equation at respec- tive locations x=ξandy=η. Multiplying these two solutions together produces the fundamental solution for the two-dimensional problem. 12/11/12 957 c/ci∇cleco√y∇t2012 Peter J. Olver Figure 17.6. Fundamental Solution to the Planar Heat Equation. Proposition 17.2. The fundamental solution to the heat equation ut=γ∆ucorre- sponding to a unit delta function placed at position (ξ,η)∈R2at the initial time t0= 0 is F(t,x,y;ξ,η)=1 4πγte−[(x−ξ)2+(y−η)2]/(4γt). (17.72) Proof: Since we already know that both function (17.71) are soluti ons to the one- dimensionalheatequation,Lemma17.1guaranteesthatthei rproduct u(t,x,y) =v(t,x)w(t,y), which equals (17.72), solves the two-dimensional heat equa tion fort >0. Moreover, at the initial time u(0,x,y) =v(0,x)w(0,y)=δ(x−ξ)δ(y−η) is a product of delta functions, and hence the result follows . Indeed, the total heat /integraldisplay/integraldisplay u(t,x,y)dxdy=/integraldisplay∞ −∞v(t,x)dx/integraldisplay∞ −∞w(t,y)dy= 1, t ≥0, remains constant, while lim t→0+u(t,x,y) =/braceleftbigg∞,(x,y) = (ξ,η), 0,otherwise . has the standard delta function limit at the initial time ins tant. Q.E.D. Figure 17.6 depicts the evolution of the fundamental soluti on when γ= 1 at successive timest=.01,.02,.05,.1. Observe that the initially concentrated heat spreads out in a radially symmetric manner. The total amount of heat remains constant. At any individual point (x,y)/\e}atio\slash= (0,0), the initially zero temperature rises slightly at first, b ut then decays monotonically back to zero at a rate proportional to 1 /t. 12/11/12 958 c/ci∇cleco√y∇t2012 Peter J. Olver Both the one- and two-dimensional fundamental solutions ha ve a bell-shaped profile known as a Gaussian function. The most important difference is the initial facto r. In a one-dimensional medium, the fundamental solution decays i n proportion to 1 /√ t, whereas in the plane the decay is more rapid, being proportional to 1 /t. The physical explanation is that the energy is able to spread out in two independent dir ections, and hence diffuses away from its initial source more rapidly. As we shall see, th e decay in three-dimensional space is more rapid still, being proportional to t−3/2for similar reasons; see (18.103). The principal purpose of the fundamental solution is to solv e the general initial value problem. We express the initial temperature distribution a s a superposition of delta func- tion sources, u(0,x,y) =f(x,y) =/integraldisplay/integraldisplay f(ξ,η)δ(x−ξ,y−η)dξdη, where, at the point ( ξ,η)∈R2, the source has magnitude f(ξ,η). Linearity implies that the solution is then given by the same superposition of funda mental solutions. Theorem 17.3. The solution to the initial value problem ut=γ∆u, u (t,x,y) =f(x,y), (x,y)∈R2, for the planar heat equation is given by the linear superposi tion formula u(t,x,y) =1 4πγt/integraldisplay/integraldisplay f(ξ,η)e−[(x−ξ)2+(y−η)2]/(4γt)dξdη. (17.73) We can interpret the solution formula (17.73) as a two-dimen sionalconvolution u(t,x,y) =F(t,x,y)∗f(x,y) (17 .74) of the initial data with a one-parameter family of progressi vely wider and shorter Gaussian filters F(t,x,y) =F(t,x,y;0,0)=1 4πγte−[x2+y2]/(4γt). (17.75) As in (13.128), such a convolution can be interpreted as a wei ghted averaging of the function, which has the effect of smoothing out the initial si gnalf(x,y). Example 17.4. If our initial temperature distribution is constant on a cir cular region, say u(0,x,y) =/braceleftbigg1x2+y2<1, 0,otherwise , Then the solution can be evaluated using (17.73), as follows : u(t,x,y) =1 4πt/integraldisplay/integraldisplay De−[(x−ξ)2+(y−η)2]/(4t)dξdη, where the integral is over the unit disk D={ξ2+η2≤1}. Unfortunately, the integral cannot be expressed in terms of elementary functions. On the other hand, numerical evaluation of the integral is straightforward. A plot of the resulting radially symmetric solution appears in Figure h2disk . 12/11/12 959 c/ci∇cleco√y∇t2012 Peter J. Olver For more general configurations, when analytical formulas a re no longer available, one turns to numerical approximation methods. The most popu lar are based on a two- dimensional variant of the Crank–Nicholson scheme (14.159 ), relying on either finite dif- ferences or finite elements to discretize the space coordina tes. We do not have space to develop the details, and refer the interested reader to [ 35,153,109]. 17.4. The Planar Wave Equation. The second important class of dynamical equations are those governing vibrational motions. The simplest planar system of this type is the two-d imensional wave equation ∂2u ∂t2=c2∆u=c2/parenleftbigg∂2u ∂x2+∂2u ∂y2/parenrightbigg , (17.76) which models the free (unforced) vibrations of a uniform two -dimensional membrane, e.g., a drum. Here u(t,x,y)represents the displacement of the membrane at time tand position (x,y)∈Ω, where the domain Ω ⊂R2represents the undeformed shape of the membrane. The constant c2>0 encapsulates the membrane’s physical properties — densit y, tension, stiffness, etc.; its square root cis called, as in the one-dimensional case, the wave speed , since it is the speed at which localized signals propagate th rough the membrane. Remark: In this simplified model, we are only allowing small, transv erse (vertical) displacements of the membrane. Large elastic vibrations le ad to the nonlinear partial differential equations of elastodynamics, [ 8]. The bending vibrations of a flexible plate, which can be viewed as the two-dimensional version of a beam, are governed by a more complicated fourth order partial differential equation; se e Exercise . The solution u(t,x,y) to the wave equation will be uniquely specified once we impos e suitable boundary and initial conditions. The Dirichlet co nditions u=hon∂Ω, (17.77) correspond to gluing our membrane to a fixed boundary — a rim. O n the other hand, the homogeneous Neumann conditions ∂u ∂n= 0 on ∂Ω, (17.78) represent a free boundary where the membrane is not attached to any support. Mixed boundary conditions attach part of the boundary and leave th e remaining portion free to vibrate: u=honD,∂u ∂n= 0 on N, (17.79) where∂Ω =D∪NwithDandNnon-overlapping. Since the wave equation is second order in time, we also need to impose two initial conditions: u(0,x,y) =f(x,y),∂u ∂t(0,x,y) =g(x,y), (x,y)∈Ω.(17.80) 12/11/12 960 c/ci∇cleco√y∇t2012 Peter J. Olver The first specifies the initial displacement of the membrane, while the second prescribes its initial velocity. The wave equation is the simplest example of a general second order system of New- tonian form ∂2u ∂t2=−K[u] =−∇∗◦∇u, (17.81) as presented in Section 14.7. As in (17.15), adopting genera l weighted inner products /a\}b∇acketle{tu;/tildewideu/a\}b∇acket∇i}ht=/integraldisplay/integraldisplay Ωu(x,y)/tildewideu(x,y)ρ(x,y)dxdy,/a\}b∇acketle{t/a\}b∇acketle{tv;/tildewidev/a\}b∇acket∇i}ht/a\}b∇acket∇i}ht=/integraldisplay/integraldisplay Ωv(x,y)·/tildewidev(x,y)κ(x,y)dxdy, (17.82) on, respectively, the spaces of scalar and vector fields, the adjoint to the gradient is a rescaled version of the divergence operator ∇∗v=−1 ρ∇·(κv). Therefore, the system (17.81) assumes the self-adjoint for m utt=−K[u] =1 ρ∇·(κv), or, in full detail, ρ(x,y)∂2u ∂t2=∂ ∂x/parenleftbigg κ(x,y)∂u ∂x/parenrightbigg +∂ ∂y/parenleftbigg κ(x,y)∂u ∂y/parenrightbigg . (17.83) Theresultinghyperbolicpartialdifferentialequationmod elsthesmalltransversevibrations of a nonuniform membrane, in which ρ(x,y)>0 represents the density of the membrane at the point ( x,y)∈Ω, while κ(x,y)>0 represents its stiffness — in direct analogy with the one-dimensional version (14.192). In particular, if the ma terial is homogeneous, then both ρandκare constant, and (17.83) reduces to the two-dimensional wa ve equation (17.76) with wave speed c=/radicalbigg κ ρ. (17.84) Separation of Variables AccordingtothegeneralframeworkestablishedinSection1 4.7,theseparablesolutions to the vibration equation (17.81) have the trigonometric fo rm uk(t,x,y) = cosωkt vk(x,y) and /tildewideuk(t,x,y) = cosωkt vk(x,y),(17.85) in which vk(x,y) is an eigenfunction of the assocaited boundary value probl em K[v] =ω2 kv=λkv. (17.86) The eigenvalue λk=ω2 kequals the square of the vibrational frequency. As in (17.20 ), there are an infinite number of such normal modes , having progressively faster and faster vibrational frequencies: ωk→ ∞ask→ ∞. In addition, in the positive semi-definite case 12/11/12 961 c/ci∇cleco√y∇t2012 Peter J. Olver — which occurs under homogeneous Neumann boundary conditio ns — there is a single constant null eigenfunction, leading to the additional sep arable solutions u0(t,x,y) = 1 and /tildewideu0(t,x,y) =t. (17.87) The first represents a stationary membrane that has been disp laced by a fixed amount in the vertical direction, while the second represents a membr ane that is moving off in the vertical direction with constant unit speed. (Think of the m embrane moving in outer space unaffected by any external gravitational force.) Specializ ing to the wave equation (17.76), the required eigenvalue problem (17.86) reduces to the Helm holtz eigenvalue problem c2∆v+λv=c2(vxx+vyy)+λv= 0 (17 .88) that we analyzed earlier in this chapter. Inthestablecases(Drichlet ormixed), thegeneral solutio ntotheinitialvalueproblem can be built up as a quasi-periodic eigenfunction series u(t,x,y) =∞/summationdisplay k=1akuk(t,x,y)+bk/tildewideuk(t,x,y) =∞/summationdisplay k=1/parenleftbig akcosωkt+bksinωkt/parenrightbig vk(x,y) (17.89) in the fundamental vibrational modes. The coefficients ak,bkare prescribed by the initial conditions: ∞/summationdisplay k=1akvk(x,y) =f(x,y),∞/summationdisplay k=1ωkbkvk(x,y) =g(x,y), (17.90) whence, by orthogonality of the eigenfunctions, ak=/a\}b∇acketle{tf;vk/a\}b∇acket∇i}ht /ba∇dblvk/ba∇dbl2=/integraldisplay/integraldisplay Ωf vkρ dxdy /integraldisplay/integraldisplay Ωv2 kρ dxdy, bk=1 ωk/a\}b∇acketle{tg;vk/a\}b∇acket∇i}ht /ba∇dblvk/ba∇dbl2=/integraldisplay/integraldisplay Ωgvkρ dxdy ωk/integraldisplay/integraldisplay Ωv2 kρ dxdy.(17.91) In the case of the wave equation, the density ρis constant, and hence can be canceled from the numerator and denominator of the orthogonality formula e (17.90). For the Neumann boundary value problem, the eigenfunction s eries solution takes an amended form u(t,x,y) =a0+b0t+∞/summationdisplay k=1/parenleftbig akcosωkt+bksinωkt/parenrightbig vk(x,y).(17.92) The coefficients ak,bkfork >0 are given by the same orthogonality formulae (17.91). The only unstable, non-periodic mode is the linearly growing te rmb0t; its coefficient b0=/a\}b∇acketle{tg;1/a\}b∇acket∇i}ht /ba∇dbl1/ba∇dbl2=/integraldisplay/integraldisplay Ωgρ dxdy /integraldisplay/integraldisplay Ωρ dxdy, 12/11/12 962 c/ci∇cleco√y∇t2012 Peter J. Olver is a weighted average of the initial velocity g(x,y) =ut(0,x,y) over the domain. In the case of the wave equation, the density ρis constant, and hence b0=1 area Ω/integraldisplay/integraldisplay Ωg(x,y)dxdy equals the average initial velocity. If the (weighted) aver age initial velocity b0/\e}atio\slash= 0 is nonzero, thenthemembranewillmoveoffatanaveragevertica lspeedb0—whilequasiperi- odically vibrating in any of the normal modes that have been e xcited by the initial dis- placement and/or initial velocity. Again, this is a two-dim ensional formulation of our observations of a free, vibrating bar — which in turn was the c ontinuum version of an unsupported mass–spring chain. Remark: An interesting question is whether two different drums can h ave identical vibrational frequencies. Or, more descriptively, can one h ear the shape of a drum? The answer turns out to be “no”, but for quite subtle reasons. See [drum] for a discussion. 17.5. Analytical Solutions of the Wave Equation. The previous section summarized the general, qualitative f eatures exhibited by solu- tions to two-dimensional vibration and wave equations. Exa ct analytical formulas are, of course, harder to come by. In this section, we analyze the two most important special cases — rectangular and circular membranes. Vibration of a Rectangular Drum Let us first consider the vibrations of a membrane in the shape of a rectangle R=/braceleftbig 0< x < a, 0< y < b/bracerightbig with side lengths aandb, whose sides are fixed to the ( x,y)–plane. Thus, we seek to solve the wave equation utt=c2∆u=c2(uxx+uyy),0< x < a, 0< y < b, (17.93) subject to the initial and boundary conditions u(t,0,y)=v(t,a,y)= 0 =v(t,x,0) =v(t,x,b), u(0,x,y)=f(x,y), ut(0,x,y) =g(x,y),0< x < a, 0< y < b.(17.94) As we saw in Section 17.2 the eigenvalues and eigenfunctions for the associated Helmholtz equation c2(vxx+vyy)+λv= 0, (x,y)∈R, (17.95) on a rectangle, subject to the homogeneous Dirichlet bounda ry conditions v(0,y) =v(a,y) = 0 =v(x,0) =v(x,b), 0< x < a, 0< y < b, (17.96) are vm,n(x,y) = sinmπx asinnπy b,where λm,n=π2c2/parenleftbiggm2 a2+n2 b2/parenrightbigg ,(17.97) 12/11/12 963 c/ci∇cleco√y∇t2012 Peter J. Olver withm,n= 1,2,.... The fundamental frequencies of vibration are the square ro ots of the eigenvalues, so ωm,n=/radicalBig λm,n=πc/radicalbigg m2 a2+n2 b2. (17.98) The frequencies will depend upon the underlying geometry — m eaning the side lengths — of the rectangle, as well as the wave speed c, which is turn is a function of the membrane’s density and stiffness, (17.84). The higher the wave speed, or the smaller the rectangle, the faster the vibrations. In layman’s terms, (17.98) quantifie s the observation that smaller, stiffer drums made of less dense material vibrate faster. According to (17.85), the normal modes of vibration of our re ctangle are um,n(t,x,y) = cosπc/radicalbigg m2 a2+n2 b2tsinmπx asinnπy b, /tildewideum,n(t,x,y) = sinπc/radicalbigg m2 a2+n2 b2tsinmπx asinnπy b.(17.99) The general solution can then be written as a double Fourier s eries u(t,x,y) =∞/summationdisplay m,n=1/bracketleftbig am,num,n(t,x,y)+bm,n/tildewideum,n(t,x,y)/bracketrightbig . in the normal modes. The coefficients am,n,bm,nare fixed by the initial displacement u(0,x,y) =f(x,y) and the initial velocity ut(0,x,y) =g(x,y), as in (17.90). The usual orthogonality relations among the eigenfunctions imply am,n=/a\}b∇acketle{tvm,n;f/a\}b∇acket∇i}ht /ba∇dblvm,n/ba∇dbl2=4 ab/integraldisplayb 0/integraldisplaya 0f(x,y) sinmπx asinnπy bdxdy, (17.100) bm,n=/a\}b∇acketle{tvm,n;g/a\}b∇acket∇i}ht ωm,n/ba∇dblvm,n/ba∇dbl2=4 πc√ m2b2+n2a2/integraldisplayb 0/integraldisplaya 0g(x,y) sinmπx asinnπy bdxdy. Since the fundamental frequencies are not rational multipl es of each other, the general solution is a genuinely quasi-periodic superposition of th e various normal modes. In Figure 17.7, we plot the solution resulting from the initi ally concentrated displace- ment† u(0,x,y) =f(x,y) =e−100[(x−.5)2+(y−.5)2] at the center of a unit square, so a=b= 1 and the wave speed c= 1. The plots are at successive times 0 ,.02,.04,...,1.6. Note that, unlike a one-dimensional string where a †The alert reader may object that the initial displacement f(x,y) does not exactly satisfy the Dirichlet boundary conditions on the edges of the rectangle. But this d oes not prevent the exis- tence of a well-defined solution to the initial value problem, whose i nitial boundary discontinuities will subsequently propagate inside the rectagle. However, they are so tiny as to be unnoticeable in the solution graphs. 12/11/12 964 c/ci∇cleco√y∇t2012 Peter J. Olver Figure 17.7. Vibrations of a Square. concentrated displacement remains concentrated at all sub sequent times and periodically repeats, the initial displacement spreads out in a radially symmetric manner and propa- gatesto the edges of therectangle, where itreflects and then interactswith itself. However, owing to the quasi-periodicity of the solution, the displac ement of the drum never exactly repeats itself, and the initial concentrated signal never q uite reforms in the rectangle’s center. Vibration of a Circular Drum Let us next analyze the vibrations of a circular membrane. As always, we build up the solution as a quasi-periodic linear combination of the n ormal modes, which, by (17.85), are specified by the eigenfunctions for the associated Helmh oltz boundary value problem. As we saw in Section 17.2, the eigenfunctions of the Helmholt z equation on a disk of radius 1, say, subject to homogeneous Dirichlet boundary conditions, are products of trigonometric and Bessel functions: vm,n(r,θ) =Jm(ζm,nr)cosmθ, /tildewidevm,n(r,θ) =Jm(ζm,nr)sinmθ,m= 0,1,2,... , n= 1,2,3,... .(17.101) Herer,θare the usual polar coordinates, while ζm,n>0 denotes the nth(positive) root of themthorder Bessel function Jm(z), cf. (17.58). The corresponding eigenvalue is its square,λm,n=ζ2 m,n, and hence the natural frequencies of vibration are equal to the Bessel roots, scaled by the wave speed: ωm,n=c/radicalbig λm,n=cζm,n. (17.102) 12/11/12 965 c/ci∇cleco√y∇t2012 Peter J. Olver Atableoftheirvalues(forthecase c= 1)canbefoundinthepreceding section. TheBessel roots do not follow any easily discernible pattern, and are c ertainly not rational multiples of each other. Thus, the vibrations of a circular drum are als o truly quasi-periodic. The frequencies ω0,n=cζ0,ncorrespond to simple eigenvalues, with a single radially symmetriceigenfunction J0(ζ0,nr), whilethe“angularmodes” ωm,n, form >0, aredouble, each possessing two linearly independent eigenfunctions ( 17.101). According to the general formula (17.85), each eigenfunction engenders two indepen dent normal modes of vibration, having the explicit forms coscζm,ntcosmθ Jm(ζm,nr),coscζm,ntsinmθ Jm(ζm,nr), sincζm,ntcosmθ Jm(ζm,nr),sincζm,ntsinmθ Jm(ζm,nr).(17.103) The general solution is written as a Fourier–Bessel series: u(t,r,θ) =1 2∞/summationdisplay n=1/bracketleftbig a0,ncoscζ0,nt+c0,nsincζ0,nt/bracketrightbig Jm(ζm,nr) +∞/summationdisplay m,n=1/bracketleftbig/parenleftbig am,ncoscζm,nt+cm,nsincζm,nt/parenrightbig cosmθ +/parenleftbig bm,ncoscζm,nt+dm,nsincζm,nt/parenrightbig sinmθ/bracketrightbig Jm(ζm,nr),(17.104) whose coefficients am,n,bm,n,cm,n,dm,nare determined, as usual, by the initial displace- ment and velocity of the membrane. In Figure vdisk , the vibrations due to an initially concentrated displacement are displayed. Again, the motio n is only quasi-periodic and never quite returns to the original configuration. Remark: As we learned in Section 14.4, the natural frequencies of vi bration a (homo- geneous) one-dimensional medium, e.g., a violin string or a column of air in a flute, are integer multiples of each other. This has the important cons equence that any resulting vibration is periodic in time. Musically, the overtones in s uch a one-dimensional instru- ment are integer multiples of each other, and so the music sou nds harmonic to our ear. On the other hand, the natural frequencies of circular and rect angular drums are irrationally related, and the vibrations are only quasi-periodic. As a re sult, we hear a percussive sound! Thus, for some reason, our appreciation of music is ps ychologically attuned to the differences between rationally related/periodic and ir rationally related/quasi-periodic vibrations. Scaling and Symmetry Symmetry methods can be effectively employed in the analysis of the wave equation. Let us consider the simultaneous rescaling t/ma√sto−→αt, x /ma√sto−→βx, y /ma√sto−→βy, (17.105) of time and space, whose effect is to change the function u(t,x,y) into a rescaled version U(t,x,y) =u(αt,βx,βy ). (17.106) 12/11/12 966 c/ci∇cleco√y∇t2012 Peter J. Olver The chain rule is employed to relate their derivatives: ∂2U ∂t2=α2∂2u ∂t2,∂2U ∂x2=β2∂2u ∂x2,∂2U ∂y2=β2∂2u ∂y2. Therefore, if usatisfies the wave equation utt=c2∆u, thenUsatisfies the rescaled wave equation Utt=α2c2 β2∆U=/tildewidec2∆U,where the rescaled wave speed is /tildewidec=αc β.(17.107) In particular, rescaling only time by setting α= 1/c, β= 1, results in a unit wave speed /tildewidec= 1. In other words, we are free to choose our unit of time measu rement so as to fix the wave speed equal to 1. If we set α=β, scaling space and time in the same proportion, then the wave speed does not change, /tildewidec=c, and so t/ma√sto−→βt, x /ma√sto−→βx, y /ma√sto−→βy, (17.108) defines a symmetry transformation for the wave equation: If u(t,x,y) is any solution to the wave equation, then so is its rescaled version U(t,x,y) =u(βt,βx,βy ) (17 .109) for any choice of scale parameter β/\e}atio\slash= 0. Observe that if u(t,x,y) is defined on a domain Ω, then the rescaled solution U(t,x,y) will be defined on the rescaled domain /tildewideΩ =1 βΩ =/braceleftbigg/parenleftbiggx β,y β/parenrightbigg/vextendsingle/vextendsingle/vextendsingle/vextendsingle(x,y)∈Ω/bracerightbigg ={(x,y)|(βx, βy)∈Ω}.(17.110) For instance, the scaling parameter β= 2 has the effect of halving the size of the domain. The normal modes for the rescaled domain have the form Un(t,x,y) =un(βt,βx,βy ) = cos(βωnt)vn(βx,βy), /tildewideUn(t,x,y) =/tildewideun(βt,βx,βy ) = sin(βωnt)vn(βx,βy), and hence the vibrational frequencies /tildewideωn=βωnare scaled by the same overall factor. Thus, when β <1, the rescaled membrane is larger by a factor 1 /β, and its vibrations are slowed down by the same factor β. For instance, a drum that is twice as large will vibrate twice as slowly, and hence have an octave lower overall tone. Musically, this means that all drums of a similar shape have the same pattern of overtones, d iffering only in their overall pitch, which is a function of their size, tautness and densit y. In particular, choosing β= 1/Rwill rescale the unit disk into a disk of radius R. The fundamental frequencies of the rescaled disk are /tildewideωm,n=β ωm,n=c Rζm,n, (17.111) 12/11/12 967 c/ci∇cleco√y∇t2012 Peter J. Olver wherecis the wave speed and ζm,nare the Bessel roots, defined in (17.58). Observe that the ratios ωm,n/ωm′,n′between vibrational frequencies remain the same, independ ent of the size of the disk Rand the wave speed c. We define the relative vibrational frequencies ρm,n=ωm,n ω0,1=ζm,n ζ0,1,in proportion to ω0,1=cζ0,1 R≈2.4c R,(17.112) which is the drum’s dominant or lowest vibrational frequenc y. The relative frequencies ρm,nare independent of the size, stiffness or composition of the d rum membrane. In the following table, we display a list ofall relativevibration alfrequencies (17.112)that are <6. Once the lowest frequency ω0,1has been determined — either theoretically, numerically or experimentally — all the higher overtones ωm,n=ρm,nω0,1are obtained by rescaling. Relative Vibrational Frequencies of a Circular Disk n/backslashbigg m0 1 2 3 4 5 6 7 8 9 ... 1 1.000 1.593 2.136 2.653 3.155 3.647 4.132 4.610 5.084 5.553... 2 2.295 2.917 3.500 4.059 4.601 5.131 5.651......... 3 3.598 4.230 4.832 5.412 5.977...... 4 4.903 5.540......... ......... 17.6. Nodal Curves. When a membrane vibrates, the individual points move up and d own in a quasi- periodic manner. As such, correlations between the motions of different points are not immediately evident. However, if the membrane is set to vibr ate in a pure eigenmode, say un(t,x,y) = cos(ωnt)vn(x,y), then all points move up and down at a common frequency ωn=/radicalbig λn, which is the square root of the eigenvalue corresponding to the eigenfunction vn(x,y). The exceptions are the points where the eigenfunction vanishes: vn(x,y) = 0. (17.113) Such points remain stationary. The set of all points ( x,y)∈Ω that satisfy (17.113) is known as the nthnodal set of the domain Ω. Scattering small particles (e.g., fine sand) over the membrane performing such a pure vibration, will ena ble us to see the nodal set, because the particles will, though random movement over the oscillating regions of the membrane, tend to accumulate along the stationary nodal cur ves. 12/11/12 968 c/ci∇cleco√y∇t2012 Peter J. Olver 1.000 1.593 2.136 2.295 2.653 2.917 3.155 3.500 3.598 Figure 17.8. Nodal Curves and Relative Vibrational Frequencies of a Circular Membrane. It can be shown that, in general, each nodal set consists of a fi nite system of nodal curves. The nodal curves intersect at critical points of the eigenf unction, where ∇vn=0, and thereby partition the membrane into nodal regions . Points lying in a common nodal region all vibrate in tandem, so that all points in a common no dal region are either up or down, except, momentarily, when the entiremembrane has zero displacement. Adjacent nodalregions, lyingontheoppositesidesofa nodalcurve, v ibrateinopposing directions— when one side is up, the other is down, and then, as the membran e becomes momentarily flat, simultaneously switch directions. Example 17.5. Circular Drums . Since the eigenfunctions (17.101) for a disk are products of trigonometric functions in the angular variabl e and Bessel functions of the radius, the nodal curves for the normal modes of vibrations o f a circular membrane are rays emanating from and circles centered at the origin. Thus , the nodal regions are annu- 12/11/12 969 c/ci∇cleco√y∇t2012 Peter J. Olver lar sectors. Pictures of the nodal curves for the first nine no rmal modes indexed by their relative frequencies are plotted in Figure 17.8. Represent ative displacements of the mem- brane in each of the first twelve modes can be found in Figure 17 .4. The dominant (lowest frequency) mode is the only one that has no nodal curves; it ha s the form of a radially symmetric bump where the entire membrane flexes up and down. E very other mode has at least one nodal curve. For instance, the next lowest modes vibrate proportionally faster at a relative frequency ρ1,1≈1.593. The most general solution with this vibrational fre- quency is a linear combination αu1,1+β/tildewideu1,1of the two eigensolutions. Each combination has a single diameter as a nodal curve, whose slope depends up on the coefficients α,β. Here, the two semicircular halves of the drum vibrate in oppo sing directions — when the top half is up, the bottom half is down and vice versa. The next set of modes have two perpendicular diameters as nodal curves; the four quadrant s of the drum vibrate in tan- dem, with opposite quadrants having the same displacements . Next in increasing order of vibrational frequency is a single mode, with a circular noda l curve whose (relative) radius ζ0,2/ζ0,1≈.43565 equals the ratio of the first two roots of the order zero B essel function; see Exercise for a justification. In this case, the inner disk and the outer annulus vibrate in opposing directions. Example 17.6. Rectangular Drums . For a general rectangular drum, the nodal curves are relatively uninteresting. Since the normal mode s (17.99) are separable products of trigonometric functions in the coordinate variables x,y, the nodal curves are regularly equi-spaced straight lines parallel to the sides of the rect angle. The internodal regions are small rectangles, all of the same size and shape, with adj acent rectangles vibrating in opposite directions. A more interesting collection of nodal curves occurs when th e rectangle admits mul- tiple eigenvalues — so-called accidental degeneracies . Two of the eigenvalues (17.97) co- incide,λm,n=λk,l, if and only if m2 a2+n2 b2=k2 a2+l2 b2(17.114) where (m,n)/\e}atio\slash= (k,l) are distinct pairs of positive integers. In such situation s, the distinct eigenmodes happen to vibrate with a common frequency ω=ωm,n=ωk,l. Consequently, any linear combination of the eigenmodes, e.g., cosωt/parenleftbigg αsinmπx asinnπy b+βsinkπx asinlπy b/parenrightbigg , α,β ∈R, is a pure vibration, and hence also qualifies as a normal mode. The associated nodal curves αsinmπx asinnπy b+βsinkπx asinlπy b= 0 (17 .115) have a more intriguing geometry, which can change dramatica lly asα,βvary. For example, on a unit square R=/braceleftbig 0< x,y < 1/bracerightbig , an accidental degeneracy occurs whenever m2+n2=k2+l2(17.116) 12/11/12 970 c/ci∇cleco√y∇t2012 Peter J. Olver Figure 17.9. Some Nodal Curves for a Square Membrane. for distinct pairs of positive integers ( m,n)/\e}atio\slash= (k,l). The simplest possibility arises when- everm/\e}atio\slash=n, in which case we can merely reverse the order, setting k=n, l=m. In Figure 17.9 we plot three sample nodal curves sin4πxsinπy+βsinπxsin4πy= 0,with β=.2,.5,1, corresponding to the three different linear combinations of the eigenfunctions with m= l= 4,n=k= 1. The associated vibrational frequency is ω4,1=πc√ 17. Remark: Classifying such accidental degeneracies takes us into th e realm of number theory, [ 10,40]. The simplest (square) case (17.116) asks one to determine all integer points that lie on a common circle. Remark: An interesting question is whether a circular disk has acci dental degenera- cies, which would occur if two different Bessel roots were to c oincide. However, it is known, [ 186; p. 129], that ζm,n/\e}atio\slash=ζk,lwhenever ( m,n)/\e}atio\slash= (k,l). Thus, a disk has no such degeneracies, and all nodal curves are circles and rays arou nd the origin. 12/11/12 971 c/ci∇cleco√y∇t2012 Peter J. Olver