leq

Chapter 15 ThePlanarLaplaceEquation The fundamental partial differential equations that govern the equilibrium mechanics of multi-dimensional media are the Laplace equation and its inhomogeneous counterpart, the Poisson equation. The Laplace equation is arguably the m ost important differential equation in all of applied mathematics. It arises in an aston ishing variety of mathematical and physical systems, ranging through fluid mechanics, elec tromagnetism, potential the- ory, solid mechanics, heat conduction, geometry, probabil ity, number theory, and on and on. The solutions to the Laplace equation are known as “harmo nic functions”, and the discovery of their many remarkable properties forms one of t he most significant chapters in the history of mathematics. In this chapter, we concentrate on the Laplace and Poisson eq uations in a two-dim- ensional (planar) domain. Their status as equilibrium equa tions implies that the solutions are determined by their values on the boundary of the domain. As in the one-dimensional equilibrium boundary value problems, the principal cases a re Dirichlet or fixed, Neumann or free, and mixed boundary conditions arise. In the introdu ctory section, we shall briefly survey the basic boundary value problems associated with th e Laplace and Poisson equa- tions. We also take the opportunity to summarize the crucial ly important tripartite clas- sification of planar second order partial differential equat ions:elliptic, such as the Laplace equation; parabolic , such as the heat equation; and hyperbolic , such as the wave equation. Each species has quite distinct properties, both analytica l and numerical, and each forms an essentially distinct discipline. Thus, by the conclusio n of this chapter, you will have encountered all three of the most important genres of partia l differential equations. The most important general purpose method for constructing explicit solutions of linear partial differential equations is the method of separ ation of variables. The method will be applied to the Laplace and Poisson equations in the tw o most important coordinate systems — rectangular and polar. Linearity implies that we m ay combine the separable solutions, and the resulting infinite series expressions wi ll play a similar role as for the heat and wave equations. In the polar coordinate case, we can , in fact, sum the infinite series in closed form, leading to the explicit Poisson integ ral formula for the solution. More sophisticated techniques, relying on complex analysis, bu t (unfortunately) only applicable to the two-dimensional case, will be deferred until Chapter 16. Green’s formula allows us to properly formulate the Laplace and Poisson equations in self-adjoint, positive definite form, and thereby characte rize the solutions via a minimiza- tion principle, first proposed by the nineteenth century mat hematician Lejeune Dirichlet, who also played a crucial role in putting Fourier analysis on a rigorous foundation. Mini- mization forms the basis of the most important numerical sol ution technique — the finite 12/11/12 806 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.1. Planar Domain. element method that we first encountered in Chapter 11. In the final section, we discuss numerical solution techniques based on finite element analy sis for the Laplace and Poisson equations and their elliptic cousins, including the Helmho ltz equation and more general positive definite boundary value problems. 15.1. The Planar Laplace Equation. The two-dimensional Laplace equation is the second order linear partial differential equation ∂2u ∂x2+∂2u ∂y2= 0. (15.1) It is named in honor of the outstanding eighteenth century Fr ench mathematician Pierre– Simon Laplace. Along with the heat and wave equations, it com pletes the trinity of truly fundamental partial differential equations. A real-valued solutionu(x,y) to the Laplace equation is known as a harmonic function . The space of harmonic functions can thus be identified as the kernel of the second order linear partial di fferential operator ∆ =∂2 ∂x2+∂2 ∂y2, (15.2) known as the Laplace operator , orLaplacian for short. The inhomogeneous or forced version, namely −∆[u] =−∂2u ∂x2−∂2u ∂y2=f(x,y) (15 .3) is known as Poisson’s equation , named for Sim´ eon–Denis Poisson, who was taught by Laplace. Poisson’s equation can be viewed as the higher dime nsional analogue of the basic equilibrium equation (11.12) for a bar. The Laplace and Poisson equations arise as the basic equilib rium equations in a re- markable variety of physical systems. For example, we may in terpretu(x,y) as the dis- placement of a membrane , e.g., a drum skin; the inhomogeneity f(x,y) in the Poisson equation represents an external forcing. Another example i s in the thermal equilibrium of flat plates; here u(x,y) represents the temperature and f(x,y) an external heat source. In fluid mechanics, u(x,y) represents the potential function whose gradient v=∇uis the velocity vector of a steady planar fluid flow. Similar cons iderations apply to two- dimensional electrostatic and gravitational potentials. The dynamical counterparts to the 12/11/12 807 c/ci∇clecopy∇t2012 Peter J. Olver ∂Ωh(x,y) Figure 15.2. Dirichlet Boundary Conditions. Laplace equation are the higher dimensional versions of the heat and wave equations, to be analyzed in Chapter 17. SinceboththeLaplaceandPoissonequationsdescribeequil ibriumconfigurations, they arise in applications in the context of boundary value probl ems. We seek a solution u(x,y) to the partial differential equation defined on a fixed bounded , open domain†(x,y)∈ Ω⊂R2. The solution is required to satisfy suitable conditions on the boundary of the domain, denoted ∂Ω, which will consist of one or more simple, closed curves, as illustrated in Figure 15.1. As in one-dimensional equilibria, there are three especially important types of boundary conditions. The first are the fixedorDirichlet boundary conditions , which specify the value of the functionuon the boundary: u(x,y) =h(x,y) for ( x,y)∈∂Ω. (15.4) The Dirichlet conditions (15.4) serve to uniquely specify t he solution u(x,y) to the Laplace or the Poisson equation. Physically, in the case of a free or f orced membrane, the Dirichlet boundary conditions correspond to gluing the edge of the mem brane to a wire at height h(x,y) over each boundary point ( x,y)∈∂Ω, as illustrated in Figure 15.2. Uniqueness meansthattheshapeoftheboundarywirewillunambiguously specifytheverticaldisplace- ment of the membrane in equilibrium. Similarly, in the model ing of thermal equilibrium, a Dirichlet boundary condition represents the imposition o f a prescribed temperature dis- tribution, represented by the function h, along the boundary of the plate. The second important class are the Neumann boundary conditions ∂u ∂n=∇u·n=k(x,y) on ∂Ω, (15.5) in which the normal derivative of the solution uon the boundary is prescribed. For exam- ple, in thermomechanics, a Neumann boundary condition spec ifies the heat flux into the †See Appendix A for the precise definitions of the terms “domain”, “bound ed”, “boundary”, etc. 12/11/12 808 c/ci∇clecopy∇t2012 Peter J. Olver platethroughitsboundary. The“no-flux”orhomogeneous Neu mannboundary conditions, wherek(x,y)≡0, correspond to a fully insulated boundary. In the case of a m embrane, homogeneous Neumann boundary conditions correspond to an u nattached edge of the drum. In fluid mechanics, the no-flux conditions imply that th e normal component of the velocity vector v=∇uvanishes on the boundary, and so no fluid is allowed to flow acro ss the solid boundary. Finally, one can mix the previous two sorts of boundary condi tions, imposing Dirichlet conditionsonpart oftheboundary, andNeumannonthecomple mentary part. Thegeneral mixed boundary value problem has the form −∆u=fin Ω, u =honD,∂u ∂n=konN, (15.6) with the boundary ∂Ω =D∪Nbeing the disjoint union of a “Dirichlet part”, denoted by D, and a “Neumann part” N. For example, if urepresents the equilibrium temperature in a plate, then the Dirichlet part of the boundary is where th e temperature is fixed, while the Neumann part is insulated, or, more generally, has presc ribed heat flux. Similarly, when modeling the displacement of a membrane, the Dirichlet part is where the edge of the drum is attached to a support, while the homogeneous Neum ann part is where it is left hanging free. Classification of Linear Partial Differential Equations in T wo Variables We have, at last, encountered all three of the fundamental li near, second order, partial differential equations for functions of two variables. The h omogeneous versions of the trinity are (a) The wave equation: utt−c2uxx= 0,hyperbolic, (b) The heat equation: ut−γuxx= 0,parabolic, (c) Laplace’s equation: uxx+uyy= 0,elliptic. The last column specifies the equations’ type, in accordance withthe standard taxonomy of partial differential equations. An explanation of the termi nology will appear momentarily. The wave, heat and Laplace equations arethe prototypical re presentatives of the three fundamental genres of partial differential equations, each with its own intrinsic features and physical manifestations. Equations governing vibrati ons, such as the wave equation, are typically hyperbolic. Equations governing diffusion, s uch as the heat equation, are parabolic. Hyperbolic and parabolic equations both govern dynamical processes, and one of the variables is identified with the time. On the other h and, equations model- ing equilibrium phenomena, including the Laplace and Poiss on equations, are typically elliptic, and only involve spatial variables. Elliptic par tial differential equations are asso- ciated with boundary value problems, whereas parabolic and hyperbolic equations require initial-boundary value problems, with, respectively, one or two required initial conditions. Furthermore, each type requires a fundamentally different k ind of numerical solution al- gorithm. While the initial tripartite classification first appears in partial differential equations in two variables, the terminology, underlying properties, and assocaited physical models 12/11/12 809 c/ci∇clecopy∇t2012 Peter J. Olver carry over to equations in higher dimensions. Most of the imp ortant partial differential equations arising in applications are of one of these three t ypes, and it is fair to say that the field of partial differential equations breaks into three major, disjoint subfields. Or, rather four subfields, the last being all the equations, incl uding higher order equations, that do not fit into this preliminary categorization. The classification of linear, second order partial different ial equations for a scalar- valued function u(x,y) of two variables†proceeds as follows. The most general such equa- tion has the form L[u] =Auxx+2Buxy+Cuyy+Dux+Euy+Fu=f, (15.7) where the coefficients A,B,C,D,E,F are all allowed to be functions of ( x,y), as is the inhomogeneity or forcing function f=f(x,y). The equation is homogeneous if and only iff≡0. We assume that at least one of the leading coefficients A,B,Cis nonzero, as otherwise the equation degenerates to a first order equation . The key quantity that determines the typeof such a partial differential equation is its discriminant ∆ =B2−AC. (15.8) This should (and for good reason) remind the reader of the dis criminant of the quadratic equation Q(ξ,η) =Aξ2+2Bξη+Cη2+Dξ+Eη+F= 0. (15.9) The solutions ( ξ,η) describes a plane curve — namely, a conic section. In the non degen- erate cases, the discriminant determines its geometrical t ype; it is •a hyperbola when ∆ >0, •a parabola when ∆ = 0, or •an ellipse when ∆ <0. This tripartite classification provides the underlying mot ivation for the terminology used to classify second order partial differential equation s. Definition 15.1. At a point ( x,y), the linear, second order partial differential equa- tion (15.7) is called (a) hyperbolic (b) parabolic (c) elliptic (d) degenerateif and only if∆(x,y)>0, ∆(x,y) = 0,butA2+B2+C2/\e}atio\slash= 0, ∆(x,y)<0, A=B=C= 0. In particular: •The wave equation uxx−uyy= 0 has discriminant ∆ = 1, and is hyperbolic. •The heat equation uxx−uy= 0 has discriminant ∆ = 0, and is parabolic. •The Poisson equation uxx+uyy=−fhas discriminant ∆ = −1, and is elliptic. †For dynamical equations, we will identify yas the time variable t. 12/11/12 810 c/ci∇clecopy∇t2012 Peter J. Olver Example 15.2. Since the coefficients in the partial differential equation ar e allowed to vary over the domain, the type of an equation may vary from p oint to point. Equations that change type are much less common, as well as being much ha rder to handle. One example arising in the theory of supersonic aerodynamics is theTricomi equation yuxx−uyy= 0. (15.10) Comparing with (15.7), we find that A=y, C =−1,andB=D=E=F=f= 0. The discriminant in this particular case is ∆ = y, and hence the equation is hyperbolic wheny>0, ellipticwhen y<0, and parabolic on the transition line y= 0. The hyperbolic region corresponds to subsonic fluid flow, while the superson ic regions are of elliptic type. The transitional parabolic boundary represents the shock l ine between the sub- and super- sonic regions. Characteristics In Section 14.5, we learned the importance of the characteri stic lines in understanding the behavior of solutions to the wave equation. Characteris tic curves play a similarly fun- damental role in the study of more general linear hyperbolic partial differential equations. Indeed, characteristics are another means of distinguishi ng between the three classes of second order partial differential equations. Definition 15.3. A smooth curve x(t)⊂R2is called a characteristic curve for the second order partial differential equation (15.7) if its tan gent vector/squaresmallsolidx= (/squaresmallsolidx,/squaresmallsolidy)Tsatisfies the quadratic characteristic equation A(x,y)/squaresmallsolidy2−2B(x,y)/squaresmallsolidx/squaresmallsolidy+C(x,y)/squaresmallsolidx2= 0. (15.11) Pay careful attention to the form of the characteristic equa tion — the positions of/squaresmallsolidx and/squaresmallsolidyare the opposite of what you might expect, while a minus sign a ppears in front of B. Furthermore, only the highest order terms in the original p artial differential equation play a role; the first and zerothorder terms are irrelevant as far as its characteristics go. For example, consider the hyperbolic wave equation† −c2uxx+uyy= 0. In this case, A=−c2, B= 0,C= 1, and so (15.11) takes the form −c2/squaresmallsolidy2+/squaresmallsolidx2= 0,which implies that/squaresmallsolidx=±c/squaresmallsolidy. All solutions to the latter ordinary differential equations are straight lines x=±cy+k, (15.12) †Warning : Here, we regard yas the “time” variable in the differential equation, rather than t, which assumes the role of the curve parameter. 12/11/12 811 c/ci∇clecopy∇t2012 Peter J. Olver wherekis an integration constant. Therefore, the wave equation ha s two characteristic curves passing through each point ( a,b), namely the straight lines (15.12) of slope ±c, in accordance with our earlier definition of characteristic s. In general, a linear partial differential equation is hyperbolic at a point ( x,y) if and only if there are two characteristic curves passing through it. Moreover, as with the wave equati on, disturbances that are concentrated near the point will tend to propagate along the characteristic curves. This fact lies at the foundation of geometric optics. Light rays m ove along characteristic curves, and are thereby subject to the optical phenomena of refracti on and focusing. On the other hand, the elliptic Laplace equation uxx+uyy= 0 has no (real) characteristic curves since the characterist ic equation (15.11) reduces to /squaresmallsolidy2+/squaresmallsolidx2= 0. Elliptic equations have no characteristics, and as a conseq uence, do not admit propagating signals; the effect of a localized disturbance, say on a membr ane, is immediately felt everywhere. Finally, for the parabolic heat equation uxx−uy= 0, the characteristic equation is simply /squaresmallsolidy2= 0, and so there is only one characteristic curve through each po int (a,b), namely the hori- zontal line y=b. Indeed, our observation that the effect of an initial concen trated heat source is immediately felt all along the bar is in accordance with propagation of localized disturbances along the characteristics. In this manner, elliptic, parabolic, and hyperbolic partia l differential equations are distinguished by the number of (real) characteristic curve s passing through a point — namely, zero, one and two, respectively. Further discussio n of characteristics and their applications to solving both linear and nonlinear partial d ifferential equations can be found in Section 22.1. 15.2. Separation of Variables. Oneoftheoldest—andstilloneofthemostwidelyused—techn iquesforconstructing explicit analytical solutions to partial differential equa tions is the method of separation of variables . We have, in fact, already used separation of variables to co nstruct particular solutions to the heat and wave equations. In each case, we sou ght a solution in the form of a product, u(t,x) =h(t)v(x), of scalar functions of each individual variable. For the heat and similar parabolic equations, h(t) was an exponential, while the wave equation chose a trigonometric function. In more general situations , we might not know in advance which function h(t) is appropriate. When the method succeeds (which is not guar anteed in advance), both factors are found as solutions to certain o rdinary differential equations. 12/11/12 812 c/ci∇clecopy∇t2012 Peter J. Olver Turning to the Laplace equation, the solution depends on xandy, and so the multi- plicative separation of variables ansatz has the form u(x,y) =v(x)w(y). (15.13) Let us see whether such a function can solve the Laplace equat ion by direct substitution. First of all, ∂2u ∂x2=v′′(x)w(y),∂2u ∂y2=v(x)w′′(y), where the primes indicate ordinary derivatives, and so ∆u=∂2u ∂x2+∂2u ∂y2=v′′(x)w(y)+v(x)w′′(y) = 0. The method will succeed if we are able to separate the variabl es by placing all of the terms involvingxon one side of the equation and all the terms involving yon the other. Here, we first write the preceding equation in the form v′′(x)w(y) =−v(x)w′′(y). Dividing both sides by v(x)w(y) (which we assume is not identically zero as otherwise the solution would be trivial) yields v′′(x) v(x)=−w′′(y) w(y), (15.14) which effectively “separates” the xandyvariables on each side of the equation. Now, how could a function of xalone be equal to a function of yalone? A moment’s reflection should convince the reader that this can happen if and only if the two functions are constant†, so v′′(x) v(x)=−w′′(y) w(y)=λ, where we use λto indicate the common separation constant . Thus, the individual factors v(x) andw(y) satisfy ordinary differential equations v′′−λv= 0, w′′+λw= 0, as promised. We already know how to solve both of these ordinary differenti al equations by elemen- tary techniques. There are three different cases, depending on the sign of the separation constantλ, each leading to four different solutions to the Laplace equa tion. We collect the entire family of separable harmonic functions together in t he following table. †Technical detail: one should assume that the underlying domain be conn ected for this to be valid; however, in practical analysis, this technicality is irrel evant. 12/11/12 813 c/ci∇clecopy∇t2012 Peter J. Olver Separable Solutions to Laplace’s Equation λ v(x) w(y) u(x,y) =v(x)w(y) λ=−ω2<0cosωx,sinωx e−ωy, eωy,eωycosωx, e−ωycosωx,eωysinωx, e−ωysinωx λ= 0 1, x 1, y 1, x, y, xy λ=ω2>0e−ωx, eωxcosωy,sinωyeωxcosωy, e−ωxcosωy,eωxsinωy, e−ωxsinωy Since Laplace’s equation is a homogeneous linear system, an y linear combination of solutions is also a solution. Thus, we can try to build genera l solutions as finite linear com- binations, or, provided we pay proper attention to converge nce issues, infinite series in the separable solutions. To solve boundary value problems, one must ensure that the resulting combination satisfies the boundary conditions. This is not e asy, unless the underlying domain has a rather specific geometry. In fact, the only domains for which we can explicitly solve bo undary value problems using the separable solutions constructed above are rectan gles. In this manner, we are led to consider boundary value problems for Laplace’s equation ∆u= 0 on a rectangle R={0<x<a, 0<y<b}. (15.15) To be completely specific, we will focus on the following Diri chlet boundary conditions: u(x,0) =f(x), u(x,b) = 0, u(0,y) = 0, u(a,y) = 0. (15.16) It will be important to only allow a nonzero boundary conditi on on one of the four sides of the rectangle. Once we know how to solve this type of proble m, we can employ linear superposition to solve the general Dirichlet boundary valu e problem on a rectangle; see Exercise for details. Other boundary conditions can be treated in a si milar fashion — with the proviso that the condition on each side of the rectan gle is either entirely Dirichlet or entirely Neumann. We will ensure that the series solution we construct satisfie s the three homogeneous boundary conditions by only using separable solutions that satisfy them. The remaining nonzero boundary condition will then specify the coefficient s of the individual summands. The function u(x,y) =v(x)w(y) will vanish on the top, right and left sides of the rectangle provided v(0) =v(a) = 0,andw(b) = 0. Referring to the preceding table, the first condition v(0) = 0 requires v(x) =  sinωx, λ =ω2>0, x, λ = 0, sinhωx, λ =−ω2<0, 12/11/12 814 c/ci∇clecopy∇t2012 Peter J. Olver where sinh z=1 2(ez−e−z) is the usual hyperbolic sine function. However, the second and third cases cannot satisfy the second boundary conditio nv(a) = 0, and so we discard them. The first case leads to the condition v(a) = sinωa= 0,and hence ωa=π,2π,3π,... is an integral multiple of π. Therefore, the separation constant λ=ω2=n2π2 a2,wheren= 1,2,3,... , (15.17) and the corresponding functions are v(x) = sinnπx a, n = 1,2,3,... . (15.18) Note: We have merely recomputed the known eigenvalues and eigenf unctions of the familiar boundary value problem v′′+λv= 0, v(0) =v(a) = 0. Sinceλ=ω2>0, the third boundary condition w(b) = 0 requires that, up to constant multiple, w(y) = sinhω(b−y) = sinhnπ(b−y) a. (15.19) Therefore, each of the separable solutions un(x,y) = sinnπx asinhnπ(b−y) a, n = 1,2,3,... , (15.20) satisfies the three homogeneous boundary conditions. It rem ains to analyze the inhomo- geneous boundary condition along the bottom edge of the rect angle. To this end, let us try a linear superposition of the separable solutions in the form of an infinite series u(x,y) =∞/summationdisplay n=1cnun(x,y) =∞/summationdisplay n=1cnsinnπx asinhnπ(b−y) a, whose coefficients c1,c2,...are to be prescribed by the remaining boundary condition. At the bottom edge, y= 0, we find u(x,0) =∞/summationdisplay n=1cnsinhnπb asinnπx a=f(x),0≤x≤a, (15.21) which takes the form of a Fourier sine series for the function f(x). According to (12.84), the coefficients bnof the Fourier sine series f(x) =∞/summationdisplay n=1bnsinnπx aare given by bn=2 a/integraldisplaya 0f(x)sinnπx adx.(15.22) Comparing (15.21,22), we discover that cnsinhnπb a=bnorcn=bn sinhnπb a=2 asinhnπb a/integraldisplaya 0f(x)sinnπx adx. 12/11/12 815 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.3. Square Membrane on a Wire. Therefore, the solution to the boundary value problem takes the form of an infinite series u(x,y) =∞/summationdisplay n=1bnsinnπx asinhnπ(b−y) a sinhnπb a, (15.23) wherebnare the Fourier sine coefficients (15.22) of f(x). Does this series actually converge to the solution to the bou ndary value problem? Fourier analysis says that, under very mild conditions on th e boundary function f(x), the answer is “yes”. Suppose that its Fourier coefficients are uni formly bounded, |bn| ≤Mfor alln≥1, (15.24) which, according to (14.23) is true whenever f(x) is piecewise continuous or, more gen- erally, integrable:/integraldisplaya 0|f(x)|dx <∞. Boundedness is also satisfied by many generalized functions, such as the delta function. In this case, as you ar e asked to prove in Exercise , the coefficients of the Fourier sine series (15.23) Bn=sinhnπ(b−y) a sinhnπb abn−→0 as n−→ ∞ (15.25) exponentially fast for all 0 <y≤b. Thus, according to Section 12.3, the solution u(x,y) is an infinitely differentiable function of xat each point in the rectangle, and can be well approximated by partial summation. The solution is also infi nitely differentiable with respect toy; see Exercise . In fact, as we shall see, the solutions to the Laplace equati on arealways analytic functions insidetheir domain ofdefinition —even when their boundary values are rather rough. Example 15.4. A membrane is stretched over a wire in the shape of a unit squar e 12/11/12 816 c/ci∇clecopy∇t2012 Peter J. Olver with one side bent in half, as graphed in Figure 15.3. The prec ise boundary conditions are u(x,y) =  x, 0≤x≤1 2, y= 0, 1−x,1 2≤x≤1, y= 0, 0, 0≤x≤1, y= 1, 0, x = 0, 0≤y≤1, 0, x = 1, 0≤y≤1. The Fourier sine series of the inhomogeneous boundary funct ion is readily computed: f(x) =/braceleftigg x, 0≤x≤1 2, 1−x,1 2≤x≤1, =4 π2/parenleftbigg sinπx−sin3πx 9+sin5πx 25−···/parenrightbigg =4 π2∞/summationdisplay m=0(−1)msin(2m+1)πx (2m+1)2. Specializing (15.23) when a=b= 1, we conclude that the solution to the boundary value problem is given by the Fourier series u(x,y) =4 π2∞/summationdisplay m=0(−1)msin(2m+1)πxsinh(2m+1)π(1−y) (2m+1)2sinh(2m+1)π. In Figure 15.3 we plot the sum of the first 10 terms in the series . This gives a reason- ably good approximation to the actual solution, except when we are very close to the raised corner of the boundary wire — which is the point of maxi mal displacement of the membrane. Polar Coordinates The method of separation of variables can be successfully ex ploited in certain other very special geometries. One particularly important case i s a circular disk. To be specific, let us take the disk to have radius 1 and centered at the origin . Consider the Dirichlet boundary value problem ∆u= 0, x2+y2<1,andu=h, x2+y2= 1,(15.26) so that the function u(x,y) satisfies the Laplace equation on the unit disk and satisfies the specified Dirichlet boundary conditions on the unit circ le. For example, u(x,y) might represent the displacement of a circular drum that is attach ed to a wire of height h(x,y) =h(cosθ,sinθ)≡h(θ),0≤θ≤2π, (15.27) above each point ( x,y) = (cosθ,sinθ) on the unit circle. The rectangular separable solutions are not particularly h elpful in this situation. The fact that we are dealing with a circular geometry inspires us to adopt polar coordinates x=rcosθ, y=rsinθ,orr=/radicalbig x2+y2, θ= tan−1y x, and write the solution u(r,θ) as a function thereof. 12/11/12 817 c/ci∇clecopy∇t2012 Peter J. Olver Warning : We will retain the same symbol, e.g., u, when rewriting a function in a different coordinate system. This is the convention of tenso r analysis and differential geometry, [ 2], that treats the function or tensor as an intrinsic object, which is concretely realized through its formula in any chosen coordinate syste m. For instance, if u(x,y) = x2+2yin rectangular coordinates, then u(r,θ) =r2cos2θ+2rsinθ— andnotr2+2θ — is its expression in polar coordinates. This convention av oids introducing new symbols when changing coordinates. We also need to relate derivatives with respect to xandyto those with respect to r andθ. Performing a standard chain rule computation, we find ∂ ∂r= cosθ∂ ∂x+sinθ∂ ∂y, ∂ ∂θ=−rsinθ∂ ∂x+rcosθ∂ ∂y,so∂ ∂x= cosθ∂ ∂r−sinθ r∂ ∂θ, ∂ ∂y= sinθ∂ ∂r+cosθ r∂ ∂θ.(15.28) These formulae allow us to rewrite the Laplace equation in po lar coordinates; after some calculation in which many of the terms cancel, we find ∆u=∂2u ∂x2+∂2u ∂y2=∂2u ∂r2+1 r∂u ∂r+1 r2∂2u ∂θ2= 0. (15.29) The boundary conditions are imposed on the unit circle r= 1, and so, by (15.27), take the form u(1,θ) =h(θ). (15.30) Keep in mind that, in order to be single-valued functions of x,y, the solution u(r,θ) and its boundary values h(θ) must both be 2 πperiodic functions of the angular coordinate: u(r,θ+2π) =u(r,θ), h (θ+2π) =h(θ). (15.31) Polar separation of variables is based on the ansatz u(r,θ) =v(r)w(θ) (15 .32) that assumes that the solution is a product of functions of th e individual polar variables. Substituting (15.32) into the polar form (15.29) of Laplace ’s equation, we find v′′(r)w(θ)+1 rv′(r)w(θ)+1 r2v(r)w′′(θ) = 0. Wenow separatevariablesby movingallthe termsinvolving ronto oneside oftheequation and all the terms involving θonto the other. This is accomplished by first multiplying the equation by r2/v(r)w(θ), and then moving the last term to the right hand side: r2v′′(r)+rv′(r) v(r)=−w′′(θ) w(θ)=λ. 12/11/12 818 c/ci∇clecopy∇t2012 Peter J. Olver As in the rectangular case, a function of rcan equal a function of θif and only if both are equal to a common separation constant, which we call λ. The partial differential equation thus splits into a pair of ordinary differential equations r2v′′+rv′−λv= 0, w′′+λw= 0, (15.33) that will prescribe the separable solution (15.32). Observ e that both have the form of eigenfunction equations in which the separation constant λplays the role of the eigenvalue, and we are only interested in nonzero solutions or eigenfunc tions. We have already solved the eigenvalue problem for w(θ). According to (15.31), w(θ+2π) =w(θ) must be a 2 πperiodic function. Therefore, according the discussion in Section 12.1, this periodic boundary value problem has the n onzero eigenfunctions 1, sinnθ, cosnθ, for n= 1,2,... . (15.34) corresponding to the eigenvalues (separation constants) λ=n2, wheren= 0,1,2,.... Fixing the value of λ, the remaining ordinary differential equation r2v′′+rv′−n2v= 0 (15.35) has the form of a second order Euler equation for the radial co mponentv(r). As discussed in Example 7.35, its solutions are obtained by substituting the power ansatz v(r) =rk. We discover that this is a solution if and only if k2−n2= 0,and hence k=±n. Therefore, for n/\e}atio\slash= 0, we find two linearly independent solutions, v1(r) =rn, v2(r) =r−n, n = 1,2,.... (15.36) Ifn= 0, there is an additional logarithmic solution v1(r) = 1, v2(r) = logr, n = 0. (15.37) Combining(15.34)and(15.36–37),weproduceacompletelis tofseparablepolarcoordinate solutions to the Laplace equation: 1, rncosnθ, rnsinnθ, logr, r−ncosnθ, r−nsinnθ,n= 1,2,3,.... (15.38) Now, the solutions in the top row of (15.38) are continuous (i n fact analytic) at the origin, whereas the solutions in the bottom row have singularities a sr→0. The latter are not relevant since we require the solution uto remain bounded and smooth — even at the center of the disk. Thus, we should only use the former to conc oct a candidate series solution u(r,θ) =a0 2+∞/summationdisplay n=1/parenleftbig anrncosnθ+bnrnsinnθ/parenrightbig (15.39) 12/11/12 819 c/ci∇clecopy∇t2012 Peter J. Olver to the Dirichlet boundary value problem. The coefficients an,bnwill be prescribed by the boundary conditions (15.30). Substituting r= 1, we find u(1,θ) =a0 2+∞/summationdisplay n=1/parenleftbig ancosnθ+bnsinnθ/parenrightbig =h(θ). We recognize this as a standard Fourier series for the 2 πperiodic function h(θ). Therefore, an=1 π/integraldisplayπ −πh(θ)cosnθdθ, bn=1 π/integraldisplayπ −πh(θ)sinnθdθ, (15.40) are precisely its Fourier coefficients, cf. (12.28). Remark: Introducing the complex variable z=reiθ=x+ iyallows us to write zn=rneinθ=rncosnθ+ irnsinnθ. (15.41) Therefore, the non-singular separable solutions are nothi ng but the harmonic polynomials we first found in Example 7.52, namely rncosnθ= Rezn, rnsinnθ= Imzn. (15.42) Exploitation of the remarkable connections between the sol utions to the Laplace equation and complex functions will form the focus of Chapter 16. In view of (15.42), the nthorder term in the series solution (15.39), anrncosnθ+bnrnsinnθ=anRezn+bnImzn= Re/bracketleftbig (an−ibn)zn/bracketrightbig , is, in fact, a homogeneous polynomial in ( x,y) of degreen. This means that, when written in rectangular coordinates xandy, (15.39) is, in fact, a power series for the function u(x,y). Proposition C.4 implies that the power series is, in fact, theTaylor series for u(x,y) based at the origin, and so its coefficients are multiples of t he derivatives of u atx=y= 0. Details are worked out in Exercise . Thus, the fact that u(x,y) has a convergent Taylor series implies that it is an analytic func tion at the origin. Indeed, as we will see, analyticity holds at any point of the domain of defin intion of a harmonic function . Example 15.5. Consider the Dirichlet boundary value problem on the unit di sk with u(1,θ) =θfor−π<θ<π. (15.43) The boundary data can be interpreted as a wire in the shape of a single turn of a spiral helix sitting over the unit circle, with a jump discontinuit y, of magnitude 2 π, at (−1,0). The required Fourier series h(θ) =θ∼2/parenleftbigg sinθ−sin2θ 2+sin3θ 3−sin4θ 4+···/parenrightbigg 12/11/12 820 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.4. Membrane Attached to a Helical Wire. ψ(x,y) Figure 15.5. Geometrical Construction of the Solution. was computed in Example 12.2. Therefore, invoking our solut ion formula (15.39–40), u(r,θ) = 2/parenleftbigg rsinθ−r2sin2θ 2+r3sin3θ 3−r4sin4θ 4+···/parenrightbigg (15.44) is the desired solution, and is plotted in Figure 15.4. In fac t, this series can be explicitly summed. In view of (15.42), u= 2 Im/parenleftbigg z−z2 2+z3 3−z4 4+···/parenrightbigg = 2 Im log(1+ z) = 2 ph(1+ z) = 2ψ,(15.45) where ψ= tan−1y 1+x(15.46) is the angle that the line passing through the two points ( x,y) and (−1,0) makes with the x-axis, as sketched in Figure 15.5. You should try to convince yourself that, on the unit 12/11/12 821 c/ci∇clecopy∇t2012 Peter J. Olver circle, 2ψ=θhas the correct boundary values. Obserfve that, even though the boundary values are discontinuous, the solution is an analytic funct ion inside the disk. Unlike the rectangular series solution (15.23), the polar s eries solution (15.39) can, in fact, be summed in closed form! If we substitute the explicit Fourier formulae (15.40) into (15.39) — remembering to change the integration variable to , say,φto avoid a notational conflict — we find u(r,θ) =a0 2+∞/summationdisplay n=1/parenleftbig anrncosnθ+bnrnsinnθ/parenrightbig =1 2π/integraldisplayπ −πh(φ)dφ (15.47) +∞/summationdisplay n=1/bracketleftbiggrncosnθ π/integraldisplayπ −πh(φ)cosnφdφ+rnsinnθ π/integraldisplayπ −πh(φ)sinnφdφ/bracketrightbigg =1 π/integraldisplayπ −πh(φ)/bracketleftigg 1 2+∞/summationdisplay n=1rn/parenleftbig cosnθcosnφ+sinnθsinnφ/parenrightbig/bracketrightigg dφ =1 π/integraldisplayπ −πh(φ)/bracketleftigg 1 2+∞/summationdisplay n=1rncosn(θ−φ)/bracketrightigg dφ. We next show how to sum the final series. Using (15.41), we can w rite it as the real part of a geometric series: 1 2+∞/summationdisplay n=1rncosnθ= Re/parenleftigg 1 2+∞/summationdisplay n=1zn/parenrightigg = Re/parenleftbigg1 2+z 1−z/parenrightbigg = Re/parenleftbigg1+z 2(1−z)/parenrightbigg = Re/parenleftbigg(1+z)(1−z) 2|1−z|2/parenrightbigg =Re(1+z−z−|z|2) 2|1−z|2=1−|z|2 2|1−z|2=1−r2 2(1+r2−2rcosθ). Substituting back into (15.47) leads to the important Poisson Integral Formula for the solution to the boundary value problem. Theorem 15.6. The solution to the Laplace equation in the unit disk subject to Dirichlet boundary conditions u(1,θ) =h(θ)is u(r,θ) =1 2π/integraldisplayπ −πh(φ)1−r2 1+r2−2rcos(θ−φ)dφ. (15.48) Example 15.7. A particularly important case is when the boundary value h(θ) =δ(θ−φ) is a delta function concentrated at the point (cos φ,sinφ),−π<φ≤π, on the unit circle. The solution to the resulting boundary value problem is the Poisson integral kernel u(r,θ) =1−r2 2π/bracketleftbig 1+r2−2rcos(θ−φ)/bracketrightbig=1−|z|2 2π|1−ze−iφ|2. (15.49) 12/11/12 822 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.6. The Poisson Kernel. The reader may enjoy verifying that this function does indee d, solve the Laplace equation and has the correct boundary values in the limit as r→1. Physically, if u(r,θ) represents the equilibrium temperature of the disk, then the delta func tion boundary data corre- spond to a concentrated unit heat source applied to a single p oint on the boundary. The resulting solution is sketched in Figure 15.6. Thus, the Poi sson kernel plays the role of the fundamental solution for the boundary value problem. Indee d, Poisson integral formula (15.48) follows from our general superposition principle, writing the boundary data as a superposition of delta functions: h(θ) =/integraldisplayπ −πh(φ)δ(φ−θ)dφ, Averaging and the Maximum Principle If we setr= 0 in the Poisson formula (15.48), then we obtain u(0,θ) =1 2π/integraldisplayπ −πh(φ)dφ. (15.50) The left hand side is the value of uat the origin — the center of the disk; the right hand side is the average of its boundary values around the unit cir cle. This is a particular instance of an important general fact. Theorem 15.8. Letu(x,y)be harmonic inside a disk of radius acentered at a point (x0,y0)with piecewise continuous (or, more generally, integrable )boundary values on the circleC={(x−x0)2+(y−y0)2=a2}. Then its value at the center of the disk is equal to the average of its values on the boundary circle : u(x0,y0) =1 2πa/contintegraldisplay Cuds=1 2π/integraldisplay2π 0u(x0+acosθ,y0+asinθ)dθ. (15.51) Proof: We use the scaling and translation symmetries of the Laplac e equation to map the disk of radius rcentered at ( x0,y0) to the unit disk centered at the origin. Specifically, we set U(x,y) =u(x0+ax,y0+ay). (15.52) 12/11/12 823 c/ci∇clecopy∇t2012 Peter J. Olver An easy chain rule computation proves that U(x,y) is harmonic on the unit disk, with boundary values h(θ) =U(cosθ,sinθ) =u(x0+acosθ,y0+asinθ). Therefore, by (15.50) , U(0,0) =1 2π/integraldisplayπ −πh(θ)dθ=1 2π/integraldisplayπ −πU(cosθ,sinθ)dθ. ReplacingUby its formula (15.52) produces the desired result. Q.E.D. An important consequence of the integral formula (15.51) is theMaximum Principle for harmonic functions. Theorem 15.9. Ifuis a nonconstant harmonic function defined on a domain Ω, thenudoes not have a local maximum or local minimum at any interior point ofΩ. Proof: The average of a continuous real function lies strictly bet ween its maximum and minimum values — except in the trivial case when the funct ion is constant. Since uis harmonic, it is continuous inside Ω. So Theorem 15.8 impli es that the value of uat (x,y) lies strictly between its maximal and minimal values on any small circle centered at (x,y). This clearly excludes the possibility of uhaving a local maximum or minimum at (x,y). Q.E.D. Thus, ona bounded domain, aharmonic function achieves itsm aximum and minimum valuesonlyatboundarypoints. Anyinteriorcriticalpoint , where∇u=0, mustbeasaddle point. Physically, if we interpret u(x,y) as the vertical displacement of a membrane, then Theorem 15.9 says that, in the absence of external forcing, t he membrane cannot have any internal bumps — its highest and lowest points are necess arily on the boundary of the domain. Thisreconfirmsour physical intuition: therestori ngforceexertedby thestretched membrane will serve to flatten any bump, and hence a membrane w ith a local maximum or minimum cannot be in equilibrium. A similar interpretati on holds for heat conduction. A body in thermal equilibrium can achieve its maximum and min imum temperature only on the boundary of the domain. Again, physically, heat energ y would flow away from any internal maximum, or towards any local minimum, and so if the body contained a local maximum or minimum on its interior, it could not be in thermal equilibrium. This concludes our discussion of separation of variables fo r the planar Laplace equa- tion. The method works in a few other special coordinate syst ems. See Exercise for one example, and [ 131,134,136] for a complete account, including connections with the underlying symmetries of the equation. 15.3. The Green’s Function. Now we turn to the Poisson equation (15.3), which is the inhom ogeneous form of the Laplace equation. In Section11.2, welearned how to solveon e-dimensional inhomogeneous boundary value problems by contructing the associated Gree n’s function. This important technique can be adapted to solve inhomogeneous boundary va lue problems for elliptic 12/11/12 824 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.7. Gaussian Distributions Converging to the Delta Function. partial differential equationsin higher dimensions, inclu ding Poisson’s equation. Asbefore, the Green’s function is characterized as the solution to the homogeneous boundary value problem in which the inhomogeneity is a concentrated unit im pulse — a delta function. The solution to the general forced boundary value problem is then obtained via linear superposition, that is, as a convolution integral with the G reen’s function. The first order of business is to establish the proper form for a unit impulse in our two-dimensional situation. We denote the delta function concentrated at position ξ= (ξ,η)∈R2by δξ(x) =δ(ξ,η)(x,y) =δ(x−ξ). (15.53) The delta function δ0(x) =δ(x,y) at the origin can be viewed as the limit, as n→ ∞, of a sequence of more and more highly concentrated functions gn(x,y), with lim n→∞gn(x,y) = 0,for (x,y)/\e}atio\slash= (0,0),while/integraldisplay/integraldisplay Ωgn(x,y)dxdy= 1. A good example of a suitable sequence is provided by the radial Gaussian distributions gn(x,y) =n πe−n(x2+y2), (15.54) which relies on the fact that /integraldisplay/integraldisplay R2e−n(x2+y2)dxdy=π n, established in Exercise A.6.6. As plotted in Figure 15.7, as n→ ∞, the Gaussian pro- files become more and more concentrated at the origin, while m aintaining a unit volume underneath their graphs. Alternatively, one can assign the delta function a dual inte rpretation as a linear func- tional on the vector space of continuous scalar-valued func tions. We formally prescribe the delta function by the integral formula /a\}b∇acketle{tδ(ξ,η);f/a\}b∇acket∇i}ht=/integraldisplay/integraldisplay Ωδ(ξ,η)(x,y)f(x,y)dxdy=/braceleftbiggf(ξ,η),(ξ,η)∈Ω, 0, (ξ,η)/\e}atio\slash∈Ω,(15.55) which holds for any continuous function f(x,y) and any domain Ω ⊂R2. As in the one-dimensional situation, we will avoid defining the integ ral when the delta function is concentrated at a boundary point, ( ξ,η)∈∂Ω, of the integration domain. 12/11/12 825 c/ci∇clecopy∇t2012 Peter J. Olver Since double integrals can be evaluated as repeated one-dim ensional integrals, we can conveniently view δ(ξ,η)(x,y) =δξ(x)δη(y) =δ(x−ξ)δ(y−η) (15 .56) as the product of a pair of one-dimensional delta functions. Indeed, if (ξ,η)∈R=/braceleftbig a<x<b, c<y<d/bracerightbig ⊂Ω is contained in a rectangle inside the domain Ω, then /integraldisplay/integraldisplay Ωδ(ξ,η)(x,y)f(x,y)dxdy=/integraldisplay/integraldisplay Rδ(ξ,η)(x,y)f(x,y)dxdy =/integraldisplayb a/integraldisplayb aδ(x−ξ)δ(y−η)f(x,y)dydx=/integraldisplayb aδ(x−ξ)f(x,η)dx=f(ξ,η). To find the Green’s function, we must solve the equilibrium eq uation subject to a concentrated unit delta force at a prescribed point ξ= (ξ,η)∈Ω inside the domain. In the case of Poisson’s equation, the partial differential equ ation takes the form −∆u=δξ,or−∂2u ∂x2−∂2u ∂y2=δ(x−ξ)δ(y−η),(x,y)∈Ω,(15.57) and the solution is subject to homogeneous boundary conditi ons, either Dirichlet or mixed. (The nonuniqueness of solutions to the pure Neumann boundar y value problem precludes the existence of a Green’s function.) The resulting solutio n to the Poisson boundary value problem is denoted as G(x;ξ) =G(x,y;ξ,η), (15.58) and called the Green’s function . Thus, the Green’s function (15.58) measures the effect, at position x= (x,y), of a concentrated force applied at position ξ= (ξ,η). Once we know the Green’s function, the solution to the genera l Poisson boundary value problem −∆u=fin Ω, u = 0 on ∂Ω (15 .59) is reconstructed through a superposition principle. We reg ard the forcing function f(x,y) =/integraldisplay/integraldisplay Ωξηδ(x−ξ)δ(y−η)f(ξ,η) as a superposition of delta impulses, whose strength at each point equals the value of f there. Linearity implies that the solution to the boundary v alue problem is the correspond- ing superposition of Green’s function responses to each of t he constituent impulses. The net result is the fundamental superposition formula u(x,y) =/integraldisplay/integraldisplay ΩξηG(x,y;ξ,η)f(ξ,η) (15 .60) 12/11/12 826 c/ci∇clecopy∇t2012 Peter J. Olver for the solution. This can be verified by direct evalution: −∆u(x,y) =/integraldisplay/integraldisplay Ωξη[−∆G(x,y;ξ,η)]f(ξ,η) =/integraldisplay/integraldisplay Ωξηδ(x−ξ,y−η)f(ξ,η) =f(x,y), as claimed. As in the one-dimensional situation, self-adjointness of t he boundary value problem is manifested in the symmetry of the Green’s function under i nterchange of its arguments: G(ξ,η;x,y)=G(x,y;ξ,η). (15.61) The general proof of symmetry follows as in the one-dimensio nal version (11.91); see Ex- ercise. Symmetry has the following intriguing physical interpret ation: Let x,ξ∈Ω be any pair of points in the domain. We apply a unit impulse to the membrane at the first point, and measure its deflection at the second; the result is exactly the same as if we apply the impulse at the second point, and measure the deflection at the first! (On the other hand, the deflections at other points in the domain will typic ally bear very little connec- tion with each other.) Similarly, in electrostatics, the so lutionu(x,y) is interpreted as the electrostatic potential for a system in equilibrium. A delt a function corresponds to a point charge, e.g., an electron. The symmetry property says that t he electrostatic potential at xdue to a point charge placed at position ξis exactly the same as the potential at ξdue to a point charge at x. The reader may wish to meditate on the physical plausibilit y of these remarkable facts. Unfortunately, most Green’s functions — with a few notable e xceptions — cannot be written down in closed form. However, their intrinsic form c an be based on the following construction. Let us begin by considering the solution to th e required Poisson equation −∆u=δ(x−ξ,y−η) (15 .62) where (ξ,η)∈Ω is the point that the unit impulse force is being applied. As usual, the general solution to an inhomogeneous linear equation is a su m u(x,y) =u⋆(x,y)+z(x,y) (15 .63) of a particular solution u⋆combined with the general solution zto the corresponding homogeneous equation, namely −∆z= 0. That is,z(x,y) is an arbitrary harmonic function. We shall assume that the particular solutionu⋆(x,y) is due to the effect of the unit impulse, irrespective of any i mposed boundary conditions. Once we have determined u⋆, we shall use the freedom inherent in the harmonic constituent z(x,y) to ensure that the sum (15.63) satisfies the required boundary conditions. One way to find a particular solution u⋆is to appeal to physical intuition. First, since the delta function is concentrated at the point ξ, the solution u⋆must solve the homogeneous Laplace equation ∆ u⋆= 0 except at the point x=ξ, where we expect 12/11/12 827 c/ci∇clecopy∇t2012 Peter J. Olver it to have some sort of discontinuity. Second, since the Pois son equation is modeling a homogeneous, uniform medium (membrane, plate, gravitatio nal potential in empty space, etc.), intheabsence ofboundary conditions, theeffect ofau nitimpulseshouldonlydepend upon on the distance away from the source of the impulse. Ther efore, we expect that the desired particular solution will depend only on the radial v ariable: u⋆=u⋆(r),wherer=/ba∇dblx−ξ/ba∇dbl=/radicalbig (x−ξ)2+(y−η)2. According to (15.37), the only radially symmetric solution s to the Laplace equation are u(r) =a+blogr, (15.64) whereaandbare constants. The constant term aissmoothand harmonic everywhere, and so cannot contribute to a delta function singularity. There fore, our only chance to produce a solution with such a singularity at the point ξis to take a multiple of the logarithmic potential: u⋆=blogr. We claim that, modulo the determination of b, this gives the correct formula, so −∆u⋆=−b∆(logr) =δ(x−ξ), r =/ba∇dblx−ξ/ba∇dbl. (15.65) is the delta function for an appropriate constant b. To justify this claim, and so determine the proper value of b, we first note that, by construction, log rsolves the Laplace equation everywhere except at r= 0, i.e., at x=ξ: ∆logr= 0, r /\e}atio\slash= 0. (15.66) Secondly, if Da=/braceleftbig 0≤r≤a/bracerightbig =/braceleftbig /ba∇dblx−ξ/ba∇dbl ≤a/bracerightbig is any disk centered at ξ, then, by the divergence form (A.60) of Green’s Theorem, /integraldisplay/integraldisplay Da∆(logr)dxdy=/integraldisplay/integraldisplay Da∇·∇(logr)dxdy =/contintegraldisplay Ca∂(logr) ∂nds=/contintegraldisplay Ca∂(logr) ∂rds=/contintegraldisplay Ca1 rds=/integraldisplayπ −πdθ= 2π, whereCa=∂Da=/braceleftbig /ba∇dblx−ξ/ba∇dbl=a/bracerightbig is the boundary of the disk, i.e., the circle of radius acentered at ξ. (The identification ∂/∂n=∂/∂ron a circle can be found in Exercise A.7.7.) Thus, if Ω is any domain, then /integraldisplay/integraldisplay Ω∆(logr)dxdy=/braceleftbigg2π,ξ∈Ω, 0,ξ/\e}atio\slash∈Ω.(15.67) In the first case, when ξ∈Ω, (15.66) allows us replace the integral over Ω by an integra l over a small disk centered at ξ, and then apply the preceding identity; in the second case, (15.68) implies that the integrand vanishes on all of the dom ain, and so the integral is 0. Equations (15.66–67) are the defining properties for 2 πtimes the delta function, so ∆(logr) = 2πδ(x−ξ). (15.68) 12/11/12 828 c/ci∇clecopy∇t2012 Peter J. Olver Comparing (15.68) with (15.65), we conclude that u⋆(x,y) =−1 2πlogr=−1 2πlog/ba∇dblx−ξ/ba∇dbl=−1 4πlog/bracketleftbig (x−ξ)2+(y−η)2/bracketrightbig (15.69) is a particular solution to the Poisson equation (15.62) wit h a unit impulse force. Thelogarithmic potential (15.69) represents the gravitational potential in empty tw o- dimensional space due to a unit point mass at position ξ, or, equivalently, the two- dimensional electrostatic potential due to a point charge a tξ. The corresponding gravita- tional (electrostatic) force field is obtained by taking its gradient: F=∇/parenleftbigg −1 2πlog/ba∇dblx−ξ/ba∇dbl/parenrightbigg =−x−ξ 2π/ba∇dblx−ξ/ba∇dbl2. Note that /ba∇dblF/ba∇dbl= 1/(2π/ba∇dblx−ξ/ba∇dbl) is proportional to the inverse distance, which is the two-dimensional form of Newton’s (Coulomb’s) three-dimen sional inverse square law. The gravitational potential due to a mass, e.g., a plate, in the s hape of a domain Ω ⊂R2can be obtained by superimposing delta function sources with st rengths equalo to the density of the material at each point. The result is the potential fun ction u(x,y) =−1 4π/integraldisplay/integraldisplay Ωξηρ(ξ,η) log/bracketleftbig (x−ξ)2+(y−η)2/bracketrightbig dξdη, (15.70) in whichρ(ξ,η) denotes the density of the body at position ( ξ,η). For example, the gravitational potential due to the unit disk D={x2+y2≤1}with unit density ρ≡1 is u(x,y) =−1 4π/integraldisplay/integraldisplay Dξηlog/bracketleftbig (x−ξ)2+(y−η)2/bracketrightbig . Returning to our boundary value problem, the general soluti on to the Poisson equa- tion (15.62) can, therefore, be written in the form u(x,y) =−1 2πlog/ba∇dblx−ξ/ba∇dbl+z(x,y), (15.71) wherez(x,y) is an arbitrary harmonic function. To construct the Green’ s function for a prescribed domain, we need to choose the harmonicfunction z(x,y)so that (15.71)satisfies the relevant homogeneous boundary conditions. Let us state this result for the Dirichlet problem. Proposition 15.10. The Green’s function for the Dirichlet boundary value probl em −∆u=fonΩ, u = 0on∂Ω, has the form G(x,y;ξ,η)=−1 4πlog/bracketleftbig (x−ξ)2+(y−η)2/bracketrightbig +z(x,y) (15 .72) wherez(x,y)istheharmonicfunctionthathasthesameboundaryvaluesas thelogarithmic potential function : ∆z= 0onΩ, z (x,y) =1 4πlog/bracketleftbig (x−ξ)2+(y−η)2/bracketrightbig for(x,y)∈∂Ω. 12/11/12 829 c/ci∇clecopy∇t2012 Peter J. Olver Let us conclude this subsection by summarizing the key prope rties of the Green’s functionG(x,ξ) for the two-dimensional Poisson equation., which (a) Solves Laplace’s equation, ∆ G= 0, for all x/\e}atio\slash=ξ. (b) Has a logarithmic singularity†atx=ξ. (c) Satisfies the relevant homogeneous boundary conditions. (d) Is symmetric: G(ξ,x) =G(x,ξ). (e) Establishes the superposition formula (15.60) for a gener al forcing function. The Method of Images The preceding analysis exposes the underlying form of the Gr een’s function, but we are still left with the determination of the harmonic compon entz(x,y) required to match the logarithmic potential boundary values. There are three principal analytical techniques employed to produce explicit formulas. The first is an adapta tion of the method of separa- tionof variables, andleadsto infiniteseries expressions, similarto thoseofthe fundamental solution for the heat equation derived in Chapter 14. We will not dwell on this approach here, although a couple of the exercises ask the reader to fill in the details. The second is themethod of images and will be developed in this section. The most powerful is ba sed on the theory of conformal mappings, but must be deferred until we have learned the basics of complex analysis; the details can be found in Section 16.3 . While the first two methods only apply to a fairly limited class of domains, they do adapt straightforwardly to higher dimensionalproblems, aswellascertainothertypesofelli pticpartialdifferentialequations, whereas the method of conformal mapping is, unfortunately, restricted to two-dimensional problems involving the Laplace and Poisson equations. WealreadyknowthatthesingularpartoftheGreen’sfunctio nforthetwo-dimensional Poisson equation is provided by a logarithmicpotential. Th e problem, then, is to construct the harmonic part, called z(x,y) in (15.72), so that the sum has the correct homogeneous boundary values, or, equivalently, that z(x,y) has the same boundary values as the log- arithmic potential. In certain cases, z(x,y) can be thought of as the potential induced by one or more hypothetical electric charges (or, equivalen tly, gravitational point masses) that are located outsidethe domain Ω, arranged in such a manner that their combined electrostatic potential happens to coincide with the logar ithmic potential on the boundary of the domain. The goal, then, is to place the image charges of suitable strength in the proper positions. Here, we will only consider the case of a single image charge, located at a position η/\e}atio\slash∈Ω. We scale the logarithmic potential (15.69) by the charge s trength, and, for added flexibility, include an additional constant — the charge’s p otential baseline: z(x,y) =alog/ba∇dblx−η/ba∇dbl+b, η ∈R2\Ω. This function is harmonic inside Ω since the logarithmic pot ential is harmonic everywhere except at the singularity η, which is assumed to lies outside the domain. For the Dirichl et †Note that this is in contrast to the one-dimensional situation, where t he Green’s function is continuous at the impulse point. 12/11/12 830 c/ci∇clecopy∇t2012 Peter J. Olver ξηx Figure 15.8. Method of Images for the Unit Disk. boundary value problem, then, for each point ξ∈Ω, we must find a corresponding image pointη∈R2\Ω and constants a,b∈R, such that‡ log/ba∇dblx−ξ/ba∇dbl=alog/ba∇dblx−η/ba∇dbl+bfor all x∈∂Ω, or, equivalently, /ba∇dblx−ξ/ba∇dbl=λ/ba∇dblx−η/ba∇dblafor all x∈∂Ω, (15.73) whereλ= logb. For each fixed ξ,η,λ,a, the equation in (15.73) will, typically, implicitly prescribe a plane curve, but it is not clear that one can alway s arrange that these curves all coincide with the boundary of our domain. In order to make further progress, we appeal to a geometrical construction based upon similar triangles. We select η=cξto be a point lying on the ray through ξ. Its location is fixed so that the triangle with vertices 0,x,ηis similar to the triangle with vertices 0,ξ,x, noting that they have the same angle at the common vertex 0— see Figure 15.8. Similarity requires that the triangles’ corresponding sid es have a common ratio, and so /ba∇dblξ/ba∇dbl /ba∇dblx/ba∇dbl=/ba∇dblx/ba∇dbl /ba∇dblη/ba∇dbl=/ba∇dblx−ξ/ba∇dbl /ba∇dblx−η/ba∇dbl=λ. (15.74) The last equality implies that (15.73) holds with a= 1. Consequently, if we choose /ba∇dblη/ba∇dbl=1 /ba∇dblξ/ba∇dbl,so that η=ξ /ba∇dblξ/ba∇dbl2, (15.75) then /ba∇dblx/ba∇dbl2=/ba∇dblξ/ba∇dbl /ba∇dblη/ba∇dbl= 1. Thusxlies on the unit circle, and, as a result, λ=/ba∇dblξ/ba∇dbl. The map taking a point ξinside the disk to its image point ηdefined by (15.75) is known as inversion with respect to the unit circle. ‡To simplify the formulas, we have omitted the 1 /(2π) factor, which can easily be reinstated at the end of the analysis. 12/11/12 831 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.9. Green’s Function for the Unit Disk. We have now demonstrated that the functions 1 2πlog/ba∇dblx−ξ/ba∇dbl=1 2πlog/parenleftbig /ba∇dblξ/ba∇dbl /ba∇dblx−η/ba∇dbl/parenrightbig =1 2πlog/ba∇dbl/ba∇dblξ/ba∇dbl2x−ξ/ba∇dbl /ba∇dblξ/ba∇dblwhen/ba∇dblx/ba∇dbl= 1, (15.76) has the same boundary values on the unit circle. Consequentl y, their difference G(x;ξ) =−1 2πlog/ba∇dblx−ξ/ba∇dbl+1 2πlog/ba∇dbl/ba∇dblξ/ba∇dbl2x−ξ/ba∇dbl /ba∇dblξ/ba∇dbl=1 2πlog/ba∇dbl/ba∇dblξ/ba∇dbl2x−ξ/ba∇dbl /ba∇dblξ/ba∇dbl /ba∇dblx−ξ/ba∇dbl(15.77) has the required properties for the Green’s function for the Dirichlet problem on the unit disk. In terms of polar coordinates x= (rcosθ,rsinθ),ξ= (ρcosϕ,ρsinϕ), applying the Law of Cosines to the triangles in Figure 15.8 le ads to the explicit formula G(r,θ;ρ,ϕ) =1 4πlog/parenleftbigg1+r2ρ2−2rρcos(θ−ϕ) r2+ρ2−2rρcos(θ−ϕ)/parenrightbigg . (15.78) InFigure15.9wesketchtheGreen’s functioncorresponding toaunitimpulsebeing applied at a point half way between the center and the edge of the disk. Remark: Unlikeone-dimensionalboundaryvalueproblems,theGree n’sfunction(15.78) has a singularity and is not continuous at the impulse point x=ξ. Applyingthegeneral superposition rule(15.60),wearrive atasolutiontotheDirichlet boundary value problem for the Poisson equation in the unit d isk. Theorem 15.11. The solution to the homogeneous Dirichlet boundary value pr ob- lem −∆u=f,forr=/ba∇dblx/ba∇dbl<1, u = 0,forr= 1, is, when expressed in polar coordinates, u(r,θ) =1 4π/integraldisplay2π 0/integraldisplay1 0f(ρ,ϕ) log/parenleftbigg1+r2ρ2−2rρcos(θ−ϕ) r2+ρ2−2rρcos(θ−ϕ)/parenrightbigg ρdρdϕ. (15.79) 12/11/12 832 c/ci∇clecopy∇t2012 Peter J. Olver The Green’s function was originally designed for the homoge neous boundary value problem. Interestingly, it can also be used to handle inhomo geneous boundary conditions. Theorem 15.12. LetG(x;ξ)denote the Green’s function for the homogeneous Dirichlet boundary value problem for the Poisson equation o n a domain Ω⊂R2. Then the solution to the inhomogeneous Dirichlet problem −∆u=f, x∈Ω, u =h, x∈∂Ω, (15.80) is given by u(x) =/integraldisplay/integraldisplay ΩξηG(x;ξ)f(ξ)−/contintegraldisplay ∂Ω∂G(x;ξ) ∂nh(ξ)ds. (15.81) For example, applying (15.81) to the Green’s function (15.7 8) for the unit disk with f≡0 recovers the Poisson integral formula (15.48). Proof: Letψ(x) be any function such that ψ=hforx∈∂Ω. Setv=u−ψ, so thatvsatisfies the homogeneous boundary value problem −∆v=f+∆ψin Ω, v = 0 on ∂Ω. We can therefore express v(x) =/integraldisplay/integraldisplay ΩξηG(x;ξ)/bracketleftbig f(ξ)+∆ψ(ξ)/bracketrightbig . (15.82) Integrationby parts, based onthe second formula inExercis eA.7.6, canbe used tosimplify the integral: /integraldisplay/integraldisplay ΩG(x;ξ)∆ψ(ξ)dξdη=/integraldisplay/integraldisplay Ωξη∆G(x;ξ)ψ(ξ)+ +/contintegraldisplay ∂Ω/parenleftbigg G(x;ξ)∂ψ(ξ) ∂n−∂G(x;ξ) ∂nψ(ξ)/parenrightbigg ds. Since the Green’s function solves −∆G=δξ, the first term reproduces −ψ(x). Moreover, G= 0 andψ=hon∂Ω, and so the right hand side of (15.82) reduces to the desired formula (15.81). Q.E.D. 15.4. Adjoints and Minimum Principles. In thissection, weexplainhow theLaplace and Poissonequat ionsfit into our universal self-adjoint equilibrium framework. The most important ou tcome will be to establish a very famous minimization principle characterizing the equ ilibrium solution, that we will exploit in the design of the finite element numerical solutio n method. The one-dimensional version of the Poisson equation, −d2u dx2=f, 12/11/12 833 c/ci∇clecopy∇t2012 Peter J. Olver is the equilibrium equation for a uniform elastic bar. In Sec tion 11.3, we wrote the under- lying boundary value problems in self-adjoint form K[u] =D∗◦D[u] =f based on the product of the derivative operator Du=u′and its adjoint D∗=−Dwith respect to the standard L2inner product. For the two-dimensional Poisson equation −∆[u] =−∂2u ∂x2−∂2u ∂y2=f(x,y) the role of the one-dimensional derivative Dwill be played by the gradient operator ∇u= gradu=/parenleftbigg ux uy/parenrightbigg . The gradient ∇defines a linear map that takes a scalar-valued function u(x,y) to the vector-valued function consisting of its two first order par tial derivatives. Thus, its domain is the vector space U= C1(Ω,R) consisting of all continuously differentiable functions u(x,y) defined for ( x,y)∈Ω. The target space V= C0(Ω,R2) consists of all continuous vector-valued functions v(x,y) = (v1(x,y),v2(x,y))T, also known as vector fields . (By way of analogy, scalar-valued functions are sometimes refe rred to as scalar fields .) Indeed, ifu1,u2∈Uare any two scalar functions and c1,c2∈Rany constants, then ∇(c1u1+c2u2) =c1∇u1+c2∇u2, which is the requirement for linearity as stated in Definitio n 7.1. In accordance with the general Definition 7.53, the adjoint o f the gradient must go in the reverse direction, ∇∗:V−→U, mapping vector fields v(x,y) to scalar functions z(x,y) =∇∗v. The defining equation for the adjoint /a\}b∇acketle{t/a\}b∇acketle{t∇u;v/a\}b∇acket∇i}ht/a\}b∇acket∇i}ht=/a\}b∇acketle{tu;∇∗v/a\}b∇acket∇i}ht (15.83) depends on the choice of inner products on the two vector spac es. The simplest inner product between real-valued scalar functions u(x,y),/tildewideu(x,y) defined on a domain Ω ⊂R2 is given by the double integral /a\}b∇acketle{tu;/tildewideu/a\}b∇acket∇i}ht=/integraldisplay/integraldisplay Ωu(x,y)/tildewideu(x,y)dxdy. (15.84) As in the one-dimensional case (3.12), this is often referre d to as the L2inner product between scalar fields, with associated norm /ba∇dblu/ba∇dbl=/radicalbig /a\}b∇acketle{tu;u/a\}b∇acket∇i}ht=/radicaligg/integraldisplay/integraldisplay Ωu(x,y)2dxdy. 12/11/12 834 c/ci∇clecopy∇t2012 Peter J. Olver Similarly, the L2inner product between vector-valued functions (vector fiel ds) defined on Ω is obtained by integrating their usual dot product: /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)dxdy=/integraldisplay/integraldisplay Ω/bracketleftbig v1(x,y)/tildewidev1(x,y)+v2(x,y)/tildewidev2(x,y)/bracketrightbig dxdy. (15.85) These form the two most basic inner products on the spaces of s calar and vector fields, and are the ones required to place the Laplace and Poisson equ ations in self-adjoint form. The adjoint identity (15.83) is supposed to hold for all appr opriate scalar fields uand vector fields v. For the L2inner products (15.84,85), the two sides of the identity rea d /a\}b∇acketle{t/a\}b∇acketle{t∇u;v/a\}b∇acket∇i}ht/a\}b∇acket∇i}ht=/integraldisplay/integraldisplay Ω∇u·vdxdy=/integraldisplay/integraldisplay Ω/parenleftbigg∂u ∂xv1+∂u ∂yv2/parenrightbigg dxdy, /a\}b∇acketle{tu;∇∗v/a\}b∇acket∇i}ht=/integraldisplay/integraldisplay Ωu∇∗vdxdy. Thus, to equate these two double integrals, we must soomehow remove the derivatives from the scalar field u. As in the one-dimensional computation (11.74), the secret is integration by parts. For single integrals, the integration by parts formula is fo und by applying the Fun- damental Theorem of Calculus to Leibniz’s rule for the deriv ative of the product of two functions. According to Appendix A, Green’s Theorem A.26 pl ays the role of the Funda- mental Theorem when dealing with double integrals. We will fi nd the divergence form /integraldisplay/integraldisplay Ω∇·vdxdy=/contintegraldisplay ∂Ωv·nds, (15.86) as in (A.60), the more convenient for the present purposes. P roceeding in analogy with the one-dimensional argument, we replace the vector field vby the product uvof a scalar fielduand a vector field v. An elementary computation proves that ∇·(uv) =u∇·v+∇u·v. (15.87) As a result, we deduce what is usually known as Green’s formula /integraldisplay/integraldisplay Ω/bracketleftbig u∇·v+∇u·v/bracketrightbig dxdy=/contintegraldisplay ∂Ωu(v·n)ds, (15.88) which is valid for arbitrary bounded domains Ω, and arbitrar y scalar and vector fields defined thereon. Rearranging the terms in this integral iden tity produces the required integration by parts formula for double integrals: /integraldisplay/integraldisplay Ω∇u·vdxdy=/contintegraldisplay ∂Ωu(v·n)ds−/integraldisplay/integraldisplay Ωu∇·vdxdy. (15.89) The terms in this identity have direct counterparts in our on e-dimensional integration by parts formula (11.77). The first term on the right hand side of this identity is a boundary term, just likethe first terms onthe right hand side of the one -dimensional formula (11.77). Moreover, the derivative operation has moved from a gradien t on the scalar field in the 12/11/12 835 c/ci∇clecopy∇t2012 Peter J. Olver doubl.e integral on the left to a divergence on the vector fiel d in the double integral on the right — even the minus sign is in place! Now, the left hand side in the integration by parts formula (1 5.89) is the same as the left hand side of (15.83). If the boundary integral vanishes , /contintegraldisplay ∂Ωuv·nds= 0, (15.90) then the right hand side of formula (15.89) also reduces to an L2inner product −/integraldisplay/integraldisplay Ωu∇·vdxdy=/integraldisplay/integraldisplay Ωu(−∇·v)dxdy=/a\}b∇acketle{tu;−∇·v/a\}b∇acket∇i}ht between the scalar field uand minus thedivergence of thevector field v. Therefore, subject to the boundary constraint (15.90), the integration by part s formula reduces to the inner product identity /a\}b∇acketle{t/a\}b∇acketle{t∇u;v/a\}b∇acket∇i}ht/a\}b∇acket∇i}ht=/a\}b∇acketle{tu;−∇·v/a\}b∇acket∇i}ht. (15.91) Comparing (15.83) with (15.91), we conclude that ∇∗v=−∇·v, (15.92) and hence, when subject to the proper boundary conditions, t he adjoint of the gradient operator is minus the divergence: ∇∗=−∇·. In this manner, we are able to write the two-dimensional Poisson equation in the standard self-adj oint form −∆u=∇∗◦∇u=−∇·(∇u) =f (15.93) subject to an appropriate system of boundary conditions tha t justify (15.91). The vanishing of the boundary integral (15.90) will be ensur ed by the imposition of suitable homogeneous boundary conditions on the scalar fiel duand/or the vector field v. Clearly the line integral will vanish if either u= 0 orv·n= 0 at each point on the boundary. These lead immediately to the three principle typ es of boundary conditions. The first are the fixed or Dirichlet boundary conditions , which require u= 0 on ∂Ω. (15.94) Alternatively, we can require v·n= 0 on ∂Ω, (15.95) which requires that vbe tangent to ∂Ω at each point, and so there is no net flux across the (solid) boundary. If we identify v=∇u, then the no flux boundary condition (15.95) translates into the Neumann boundary conditions ∂u ∂n=∇u·n= 0 on ∂Ω. (15.96) One can evidently also mix the boundary conditions, imposin g Dirichlet conditions on part of the boundary, and Neumann on the complementary part: u= 0 onD,∂u ∂n= 0 onN, where ∂Ω =D∪N(15.97) 12/11/12 836 c/ci∇clecopy∇t2012 Peter J. Olver is the disjoint union of the Dirichlet and Neumann parts. More generally, when modeling inhomogeneous membranes, he at flow through inho- mogeneous media, and similar physical equilibria, we repla ce the L2inner product between vector fields (15.85) by a weighted version /a\}b∇acketle{t/a\}b∇acketle{tv;/tildewidev/a\}b∇acket∇i}ht/a\}b∇acket∇i}ht=/integraldisplay/integraldisplay Ω/bracketleftbig p(x,y)v1(x,y)/tildewidev1(x,y)+q(x,y)v2(x,y)/tildewidev2(x,y)/bracketrightbig dxdy, (15.98) in whichp(x,y),q(x,y)>0 are strictly positive functions on the domain ( x,y)∈Ω. These functions are determined by the underlying physical p roperties of the medium being modeled. RetainingtheusualL2innerproduct(15.84)betweenscalarfields, letuscompute the weighted adjoint of the gradient operator, as defined by ( 15.83). As before, we use the basic integration by parts formula (15.89) to remove the der ivatives from the scalar field u, and so /a\}b∇acketle{t/a\}b∇acketle{t∇u;v/a\}b∇acket∇i}ht/a\}b∇acket∇i}ht=/integraldisplay/integraldisplay Ω/parenleftbigg pv1∂u ∂x+qv2∂u ∂y/parenrightbigg dxdy =/contintegraldisplay ∂Ω/parenleftbig −uqv2dx+upv1dy/parenrightbig −/integraldisplay/integraldisplay Ωu/parenleftbigg∂(pv1) ∂x+∂(qv2) ∂y/parenrightbigg dxdy.(15.99) Equating the left hand side to /a\}b∇acketle{tu;∇∗v/a\}b∇acket∇i}ht, we deduce that, provided the boundary integral vanishes, the weighted adjoint of the gradient operator wit h respect to the inner products (15.84), (15.98) is given by ∇∗v=−∂(pv1) ∂x−∂(qv2) ∂y=−p∂v1 ∂x−q∂v2 ∂y−v1∂p ∂x−v2∂q ∂y. (15.100) This holds provided the scalar and vector fields satisfy suit able boundary conditions; for example, requiring either u(x,y) = 0 or v(x,y) = 0 at every boundary point ( x,y)∈∂Ω will cause the boundary integral in (15.99) to vanish, and he nce justify (15.100). As a result, all of the usual homogeneous boundary conditions — Dirichlet, Neumann or mixed — retain their validity in this more general context. T he corresponding self-adjoint boundary value problem takes the form ∇∗◦∇u=−∂ ∂x/parenleftbigg p(x,y)∂u ∂x/parenrightbigg −∂ ∂x/parenleftbigg q(x,y)∂u ∂x/parenrightbigg =f(x,y),(x,y)∈Ω,(15.101) along with the chosen boundary conditions. The partial differential equation (15.101) arises in many ot her contexts. For example, consider a steady-state fluid flow described by a vector field vmoving in a domain Ω ⊂R2. The flow is called irrotational if has zero curl, ∇×v=0, and hence, assuming Ω is simply connected, is a gradient v=∇u. The function u(x,y) is known as the fluid velocity potential . The constitutive assumptions connect the fluid velocity wi th its rate of flow w=ρv, whereρ(x,y)>0 is the scalar density of the fluid. Conservation of mass prov ides the final equation, namely ∇·w+f= 0, where f(x,y) represents fluid sources f >0 or sinksf <0. Therefore, the basic equilibrium equations take the form −∇·(ρ∇u) =f,or−∂ ∂x/parenleftbigg ρ(x,y)∂u ∂x/parenrightbigg −∂ ∂y/parenleftbigg ρ(x,y)∂u ∂y/parenrightbigg =f(x,y),(15.102) 12/11/12 837 c/ci∇clecopy∇t2012 Peter J. Olver which is (15.101) with p=q=ρ. The most common case of a homogeneous (constant density) fluid thus reduces to the Poisson equation (15.3), w ithfreplaced by f/ρ. In electrostatics, the gradient equation v=∇urelates the voltage drop to the elec- trostatic potential u, and is the continuous analog of the circuit formula (6.18) r elating potentials to voltages. The continuous version of Kirchhoff ’s Voltage Law (6.20) that the net voltage drop around any loop is zero is the fact that any gr adient vector has zero curl, ∇×v=0, i.e., the flow is irrotational. Ohm’s law (6.23) has the form y=Cvwhere the vector field yrepresents the current, while C= diag(p(x,y),q(x,y))represents th conduc- tance of the medium; in the case of Laplace’s equation, we are assuming a uniform unit conductance. Finally, the equation f=∇·y=∇∗vrelating current and external current sources forms the continuous analog of Kirchhoff’s Current L aw (6.26) — the transpose of the discrete incidence matrix translates into the adjoin t of the gradient operator is the divergence. Thus, our discrete electro-mechanical analog y carries over, in the continu- ous realm, to a tripartite electro-mechanical-fluid analog y, with all three physical systems leading to the exact same general mathematical structure. Positive Definiteness and the Dirichlet Principle In conclusion, as a result of the integration by parts calcul ation, we have formulated the Poisson and Laplace equations (as well as their weighted counterparts) in positive (semi-)definite, self-adjoint form −∆u=∇∗◦∇u=f, when subject to the appropriate homogeneous boundary condi tions: Dirichlet, Neumann, or mixed. A key benefit is, in the positive definite cases, the c haracterization of the solutions by a minimization principle. According to Theorem 7.60, the self-adjoint operator ∇∗◦∇is positive definite if and only if the kernel of the underlying gradient operator — rest ricted to the appropriate space of scalar fields — is trivial: ker ∇={0}. The determination of the kernel of the gradient operator relies on the following elementary fact. Lemma 15.13. Ifu(x,y)is aC1function defined on a connected domain Ω, then ∇u≡0if and only if u≡cis a constant. This result can be viewed as the multi-variable counterpart of the result that the only function with zero derivative is a constant. It is a simple co nsequence of Theorem A.20; see Exercise . Therefore, the only functions which could show up in ker ∇, and thus prevent positive definiteness, are the constants. The boundary cond itions will tell us whether or not this occurs. The only constant function that satisfies ei ther homogeneous Dirichlet or homogeneous mixed boundary conditions is the zero function , and thus, just as in the one- dimensional case, the boundary value problem for the Poisso n equation with Dirichlet or mixed boundary conditions is positive definite. On the other hand, any constant function satisfies the homogeneous Neumann boundary conditions ∂u/∂n= 0, and hence such boundary value problems are only positive semi-definite. In the positive definite cases, the equilibrium solution is c haracterized by our basic minimization principle (7.81). For the Poisson equation, t he result is the justly famous Dirichlet minimization principle . 12/11/12 838 c/ci∇clecopy∇t2012 Peter J. Olver Theorem 15.14. The function u(x,y)that minimizes the Dirichlet integral P[u] =1 2/ba∇dbl∇u/ba∇dbl2−/a\}b∇acketle{tu;f/a\}b∇acket∇i}ht=/integraldisplay/integraldisplay Ω/parenleftbig1 2u2 x+1 2u2 y−fu/parenrightbig dxdy (15.103) among all C1functions that satisfy the prescribed homogeneous Dirichl et or mixed bound- ary conditions is the solution to the corresponding boundar y value problem for the Poisson equation −∆u=f. In physical applications, the Dirichlet integral (15.103) represents the energy in the system. As always, Nature chooses the equilibrium configura tion so as to minimize the energy. AkeyapplicationoftheDirichletminimumprincipl eisthefiniteelementnumerical solution scheme, to be described in detail in Section 15.5. Remark: The fact that a minimizer to the Dirichlet integral (15.103 ) satisfies the Poisson equation is an immediate consequence of our general Minimization Theorem 7.62. However, unlike the finite-dimensional situation, proving theexistence of a minimizing function is a non-trivial issue. This was not immediatiely r ecognized: Dirichlet originally thought this to be self-evident, but it then took about 50 yea rs until Hilbert supplied the first rigorous existence proof. In this introductory tre atment, we adopt a pragmatic approach, concentrating on the computation of the solution — reassured, if necessary, by the theoreticians’ efforts in establishing its existence. The Dirichlet minimization principle (15.103) was derived under the assumption that the boundary conditions are homogeneous — either pure Diric hlet or mixed. As it turns out, the principle, as stated, also applies to inhomogeneou s Dirichlet boundary conditions. However, if we have a mixed boundary value problem with inhom ogeneous Neumann con- ditions on part of the boundary, then we must include an addit ional boundary term in the minimizing functional. The general result can be stated as f ollows: Theorem 15.15. The solution u(x,y)to the boundary value problem −∆u=fin Ω, u =honD,∂u ∂n=konN, with∂Ω =D∪N, andD/\e}atio\slash=∅, is characterized as the unique function that minimizes the modified Dirichlet integral /hatwideP[u] =/integraldisplay/integraldisplay Ω/parenleftbig1 2/ba∇dbl∇u/ba∇dbl2−fu/parenrightbig dxdy+/integraldisplay Nukds (15.104) among all C1functions that satisfy the prescribed boundary conditions . The inhomogeneous Dirichlet problem has N=∅andD=∂Ω, in which case the boundary integral does not appear. Exercise asks you to prove this result. As we know, positive definiteness is directly related to the s tability of the physical system. The Dirichlet and mixed boundary value problems are stable, and can support any imposed force. On the other hand, the pure Neumann bounda ry value problem is unstable, owingtotheexistenceofanontrivialkernel—the constantfunctions. Physically, 12/11/12 839 c/ci∇clecopy∇t2012 Peter J. Olver the unstable mode represents a rigid translation of the enti re membrane in the vertical direction. Indeed, the Neumann problem leaves the entire bo undary of the membrane unattached to any support, and so the unforced membrane is fr ee to move up or down without affecting its equilibrium status. Furthermore, asinfinite-dimensionallinearsystems, non- uniquenessandnon-existence of solutions go hand in hand. As we learned in Section 11.3, th e existence of a solution to a Neumann boundary value problem is subject to the Fredholm alternative , suitably adapted to this multi-dimensional situation. A necessary c ondition for the existence of a solution is that the forcing function be orthogonal to the e lements of the kernel of the underlying self-adjoint linear operator, which, in the pre sent situation requires that fbe orthogonal to the subspace consisting of all constant funct ions. In practical terms, we only need to check orthogonality with respect to a basis for the su bspace, which in this situation consists of the constant function 1. Theorem 15.16. The Neumann boundary value problem −∆u=f,inΩ,∂u ∂n= 0,on∂Ω, (15.105) admits a solution u(x,y)if and only if /a\}b∇acketle{t1;f/a\}b∇acket∇i}ht=/integraldisplay/integraldisplay Ωf(x,y)dxdy= 0. (15.106) Moreover, whenitexists,thesolutionisnotuniquesincean yfunctionoftheform u(x,y)+c, wherec∈Ris an arbitrary constant, is also a solution. Forcing functions f(x,y) which do not satisfy the orthogonality constraint (15.106 ) will excite the translational instability, and no equilibr ium configuration is possible. For example, if we force a free membrane, (15.106) requires that the net force in the verti- cal direction be zero; otherwise, the membrane will start mo ving and cannot be in an equilibrium. 15.5. Finite Elements. As the reader has no doubt already surmised, explicit soluti ons to boundary value problems for the Laplace and Poisson equations are few and fa r between. In most cases, exact solution formulae are not available, or are so complic ated as to be of scant utility. To proceed further, one is forced to design suitable numeric al approximation schemes that can accurately evaluate the desired solution. Anespeciallypowerful classofnumerical algorithmsfor so lvingellipticboundary value problems are the finite element methods. We have already lear ned, in Section 11.6, the key underlying idea. One begins with a minimization princip le, prescribed by a quadratic functional defined on a suitable vector space of functions Uthat serves to incorporate the (homogeneous) boundary conditions. The desired soluti on is characterized as the unique minimizer u⋆∈U. One then restricts the functional to a suitably chosen finit e- dimensional subspace W⊂U, and seeks a minimizer w⋆∈W. Finite-dimensionality of 12/11/12 840 c/ci∇clecopy∇t2012 Peter J. Olver Whas the effect of reducing the infinite-dimensional minimiza tion problem to a finite- dimensional problem, which can then be solved by numerical l inear algebra. The resulting minimizerw⋆will — provided the subspace Whas been cleverly chosen — provide a good approximation to the true minimizer u⋆on the entire domain. Here we concentrate on the practical design of the finite element procedure, and refer t he reader to a more advanced text, e.g., [ 174], for the analytical details and proofs of convergence. Mos t of the multi- dimensional complications are not in the underlying theory , but rather in the realms of data management and organizational details. In this section, we first concentrate on applying these ideas to the two-dimensional Poisson equation. For specificity, we concentrate on the hom ogeneous Dirichlet boundary value problem −∆u=fin Ω u= 0 on∂Ω. (15.107) According to Theorem 15.14, the solution u=u⋆is characterized as the unique minimizing function for the Dirichlet functional (15.103) among all sm ooth functions u(x,y) that satisfytheprescribedboundary conditions. Inthefiniteel ement approximation, werestrict the Dirichlet functional to a suitably chosen finite-dimens ional subspace. As in the one- dimensionalsituation,themostconvenientfinite-dimensi onalsubspacesconsistoffunctions that may lack the requisite degree of smoothness that qualifi es them as possible solutions to the partial differential equation. Nevertheless, they do provide good approximations to the actual solution. An important practical considerati on, impacting the speed of the calculation, is to employ functions with small support, as i n Definition 13.5. The resulting finite element matrix will then be sparse and the solution to t he linear system can be relatively rapidly calculate, usually by application of an iterative numerical scheme such as the Gauss–Seidel or SOR methods discussed in Chapter 10. Finite Elements and Triangulation For one-dimensional boundary value problems, the finite ele ment construction rests on the introduction of a mesh a=x0<x1<···<xn=bon the interval of definition. The mesh nodes xkbreak theinterval into a collectionof small subintervals. In two-dimensional problems, a meshconsists of a finite number of points xk= (xk,yk),k= 1,...,m, known asnodes, usually lying inside the domain Ω ⊂R2. As such, there is considerable freedom in the choice of mesh nodes, and completely uniform spacing i s often not possible. We regard the nodes as forming the vertices of a triangulation of the domain Ω, consisting of a finite number of small triangles, which we denote by T1,...,TN. The nodes are split into two categories — interior nodes andboundary nodes , the latter lying on or close to the boundary of the domain. A curved boundary is approximate d by the polygon through the boundary nodes formed by the sides of the triangles lying on the edge of the domain; see Figure 15.10 for a typical example. Thus, in computer imp lementations of the finite element method, thefirst moduleisaroutinethat willautoma ticallytriangulateaspecified domain in some reasonable manner; see below for details on wh at “reasonable” entails. As in our one-dimensional finite element construction, the f unctionsw(x,y) in the finite-dimensional subspace Wwill be continuous and piecewise affine . “Piecewise affine” 12/11/12 841 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.10. Triangulation of a Planar Domain. means that, on each triangle, the graph of wis flat, and so has the formula† w(x,y) =αν+βνx+γνy,for (x,y)∈Tν. (15.108) Continuity of wrequires that its values on a common edge between two triangl es must agree, and this will impose certain compatibility conditio ns on the coefficients αµ,βµ,γµ andαν,βν,γνassociated with adjacent pairs of triangles Tµ,Tν. The graph of z=w(x,y) forms a connected polyhedral surface whose triangular face s lie above the triangles in the domain; see Figure 15.10 for an illustration. The next step is to choose a basis of the subspace of piecewise affine functions for the given triangulation. As in the one-dimensional version, th e most convenient basis consists ofpyramid functions ϕk(x,y) which assume the value 1 at a single node xk, and are zero at all the other nodes; thus ϕk(xi,yi) =/braceleftbigg1, i=k, 0, i/\e}atio\slash=k.(15.109) Note thatϕkwill be nonzero only on those triangles which have the node xkas one of their vertices, and hence the graph of ϕklooks like a pyramid of unit height sitting on a flat plane, as illustrated in Figure 15.12. The pyramid functions ϕk(x,y) corresponding to the interior nodes xkautomatically satisfy the homogeneous Dirichlet boundary conditions on t he boundary of the domain — or, more correctly, on the polygonal boundary of the triang ulated domain, which is †Here and subsequently, the index νis a superscript, not a power! 12/11/12 842 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.11. Piecewise Affine Function. Figure 15.12. Finite Element Pyramid Function. supposed to be a good approximation to the curved boundary of the original domain Ω. Thus, the finite-dimensional finite element subspace Wis the span of the interior node pyramid functions, and so general element w∈Wis a linear combination thereof: w(x,y) =n/summationdisplay k=1ckϕk(x,y), (15.110) where the sum ranges over the ninterior nodes of the triangulation. Owing to the original specification (15.109) of the pyramid functions, the coeffici ents ck=w(xk,yk)≈u(xk,yk), k = 1,...,n, (15.111) are thesameas the values of the finite element approximation w(x,y) at the interior 12/11/12 843 c/ci∇clecopy∇t2012 Peter J. Olver nodes. This immediately implies linear independence of the pyramid functions, since the only linear combination that vanishes at all nodes is the tri vial onec1=···=cn= 0. Thus, theinteriornodepyramidfunctions ϕ1,...ϕnformabasisforfiniteelementsubspace W, which therefore has dimension equal to n, the number of interior nodes. Determining the explicit formulae for the finite element bas is functions is not difficult. On one of the triangles Tνthat has xkas a vertex, ϕk(x,y) will be the unique affine function (15.108) that takes the value 1 at the vertex xkand 0 at its other two vertices xl andxm. Thus, we are in need of a formula for an affine function or element ων k(x,y) =αν k+βν kx+γν ky, (x,y)∈Tν, (15.112) that takes the prescribed values ων k(xi,yi) =ων k(xj,yj) = 0, ων k(xk,yk) = 1, at three distinct points. These three conditions lead to the linear system ων k(xi,yi) =αν k+βν kxi+γν kyi= 0, ων k(xj,yj) =αν k+βν kxj+γν kyj= 0, ων k(xk,yk) =αν k+βν kxk+γν kyk= 1.(15.113) The solution†produces the explicit formulae αν k=xiyj−xjyi ∆ν, βν k=yi−yj ∆ν, γν k=xj−xi ∆ν, (15.114) for the coefficients; the denominator ∆ν= det 1xiyi 1xjyj 1xkyk =±2areaTν (15.115) is, up to sign, twice the area of the triangle Tν; see Exercise . Example 15.17. Consider an isoceles right triangle Twith vertices x1= (0,0),x2= (1,0),x3= (0,1). Using (15.114–115) (or solving the linear systems (15.113) directly), we immediately pro- duce the three affine elements ω1(x,y) = 1−x−y, ω2(x,y) =x, ω3(x,y) =y. (15.116) As required, each ωkequals 1 at the vertex xkand is zero at the other two vertices. †Cramer’s Rule (1.88) comes in handy here. 12/11/12 844 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.13. Vertex Polygons. Thefiniteelementpyramidfunctionisthenobtainedbypieci ngtogethertheindividual affine elements, whence ϕk(x,y) =/braceleftbiggων k(x,y),if (x,y)∈Tνwhich has xkas a vertex, 0, otherwise.(15.117) Continuity of ϕk(x,y) is assured since the constituent affine elements have the sam e values at common vertices. The support of the pyramid function (15. 117) is the polygon suppϕk=Pk=/uniondisplay νTν (15.118) consisting of all the triangles Tνthat have the node xkas a vertex. In other words, ϕk(x,y) = 0 whenever ( x,y)/\e}atio\slash∈Pk. We will call Pkthekthvertex polygon . The node xk lies on the interior of its vertex polygon Pk, while the vertices of Pkare all those that are connected to xkby a single edge of the triangulation. In Figure 15.13 the sha ded regions indicate two of the vertex polygons for the triangulation in Figure 15.10. Example 15.18. The simplest, and most common triangulations are based on regular meshes. Suppose that the nodes lie on a square grid, a nd so are of the form xi,j= (ih+a,jh+b) whereh>0 is the inter-node spacing, and ( a,b) represents an over- all offset. If we choose the triangles to all have the same orie ntation, as in the first picture in Figure 15.14, then the vertex polygons all have the same sh ape, consisting of 6 triangles of total area 3 h2— the shaded region. On the other hand, if we choose an alterna ting, perhaps more æsthetically pleasing triangulation as in the second picture, then there are two types of vertex polygons. The first, consisting of four tr iangles, has area 2 h2, while the second, containing 8 triangles, has twice the area, 4 h2. In practice, there are good reasons to prefer the former triangulation. In general, in order to ensure convergence of the finite eleme nt solution to the true minimizer, one should choose a triangulation with the follo wing properties: (a) The triangles are not too long and skinny. In other words, th eir sides should have comparable lengths. In particular, obtuse triangles shoul d be avoided. 12/11/12 845 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.14. Square Mesh Triangulations. (b) The areas of nearby triangles Tνshould not vary too much. (c) The areas of nearby vertex polygons Pkshould also not vary too much. For adaptive or variable meshes, one might very well have wid e variations in area over the entire grid, with small triangles in regions of rapid change in the solution, and large ones in less interesting regions. But, overall, the sizes of the tri angles and vertex polygons should not dramatically vary as one moves across the domain. The Finite Element Equations WenowseektoapproximatethesolutiontothehomogeneousDi richletboundaryvalue problem by restricting the Dirichlet functional to the sele cted finite element subspace W. Substituting the formula (15.110) for a general element of Winto the quadratic Dirichlet functional (15.103) and expanding, we find P[w] =P/bracketleftiggn/summationdisplay i=1ciϕi/bracketrightigg =/integraldisplay/integraldisplay Ω /parenleftiggn/summationdisplay i=1ci∇ϕi/parenrightigg2 −f(x,y)/parenleftiggn/summationdisplay i=1ciϕi/parenrightigg dxdy =1 2n/summationdisplay i,j=1kijcicj−n/summationdisplay i=1bici=1 2cTKc−bTc. Here,K= (kij) is the symmetric n×nmatrix, while b= (b1,b2,...,bn)Tis the vector that have the respective entries kij=/a\}b∇acketle{t∇ϕi;∇ϕj/a\}b∇acket∇i}ht=/integraldisplay/integraldisplay Ω∇ϕi·∇ϕjdxdy, bi=/a\}b∇acketle{tf;ϕi/a\}b∇acket∇i}ht=/integraldisplay/integraldisplay Ωfϕidxdy.(15.119) 12/11/12 846 c/ci∇clecopy∇t2012 Peter J. Olver Thus, to determine the finite element approximation, we need to minimize the quadratic function P(c) =1 2cTKc−bTc (15.120) over all possible choices of coefficients c= (c1,c2,...,cn)T∈Rn, i.e., over all possible functionvaluesattheinterior nodes. Restricting tothefin iteelement subspace hasreduced us to a standard finite-dimensional quadratic minimization problem. First, the coefficient matrixK >0 is positive definite due to the positive definiteness of the o riginal functional; the proof in Section 11.6 is easily adapted to the present sit uation. Theorem 4.1 tells us that the minimizer is obtained by solving the associated lin ear system Kc=b. (15.121) The solution to (15.121) can be effected by either Gaussian el imination or an iterative technique. To find explicit formulae for the matrix coefficients kijin (15.119), we begin by noting that the gradient of the affine element (15.112) is equal to ∇ων k(x,y) =aν k=/parenleftbigg βν k γν k/parenrightbigg =1 ∆ν/parenleftbigg yi−yj xj−xi/parenrightbigg , (x,y)∈Tν, (15.122) which is a constant vector inside the triangle Tν, while outside ∇ων k=0. Therefore, ∇ϕk(x,y) =/braceleftigg ∇ων k=aν k,if (x,y)∈Tνwhich has xkas a vertex, 0, otherwise,(15.123) reduces to a piecewise constant function on the triangulati on. Actually, (15.123) is not quite correct since if ( x,y) lies on the boundary of a triangle Tν, then the gradient does not exist. However, this technicality will not cause any diffi culty in evaluating the ensuing integral. We will approximate integrals over the domain Ω by integrals over the triangles, which relies on our assumption that the polygonal boundary of the t riangulation is a reasonably close approximation to the true boundary ∂Ω. In particular, kij≈/summationdisplay ν/integraldisplay/integraldisplay Tν∇ϕi·∇ϕjdxdy≡/summationdisplay νkν ij. (15.124) Now, according to (15.123), one or the other gradient in the i ntegrand will vanish on the entiretriangle Tνunlessbothxiandxjarevertices. Therefore, theonlytermscontributing tothesumarethosetriangles Tνthathaveboth xiandxjasvertices. If i/\e}atio\slash=jthereareonly two such triangles, while if i=jevery triangle in the ithvertex polygon Picontributes. The individual summands are easily evaluated, since the gra dients are constant on the triangles, and so, by (15.123), kν ij=/integraldisplay/integraldisplay Tνaν i·aν jdxdy=aν i·aν jareaTν=1 2aν i·aν j|∆ν|. 12/11/12 847 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.15. Right and Equilateral Triangles. LetTνhave vertices xi,xj,xk. Then, by (15.122,123,115), kν ij=1 2(yj−yk)(yk−yi)+(xk−xj)(xi−xk) (∆ν)2|∆ν|=−(xi−xk)·(xj−xk) 2|∆ν|, i/\e}atio\slash=j, kν ii=1 2(yj−yk)2+(xk−xj)2 (∆ν)2|∆ν|=/ba∇dblxj−xk/ba∇dbl2 2|∆ν|. (15.125) In this manner, each triangle Tνspecifies a collection of 6 different coefficients, kν ij=kν ji, indexed by its vertices, and known as the elemental stiffnesses ofTν. Interestingly, the elemental stiffnesses depend only on the three vertex anglesin the triangle and not on its size. Thus, similar triangles have the sameelemental stiffnesses. Indeed, if θν i,θν j,θν k denote the angles in Tνat the respective vertices xi,xj,xk, then, according to Exercise , kν ii=1 2/parenleftbig cotθν k+cotθν j/parenrightbig ,whilekν ij=kν ji=−1 2cotθν k, i/\e}atio\slash=j.(15.126) Example 15.19. The right triangle with vertices x1= (0,0),x2= (1,0),x3= (0,1) has elemental stiffnesses k11= 1, k22=k33=1 2, k12=k21=k13=k31=−1 2, k23=k32= 0. (15.127) The same holds for any other isoceles right triangle, as long as we chose the first vertex to be at the right angle. Similarly, an equilateral triangle has all 60◦angles, and so its elemental stiffnesses are k11=k22=k33=1√ 3≈.577350, k12=k21=k13=k31=k23=k32=−1 2√ 3≈ −.288675.(15.128) Assembling the Elements The elemental stiffnesses of each triangle will contribute, through the summation (15.124), to the finite element coefficient matrix K. We begin by constructing a larger matrixK∗, which we call the full finite element matrix , of sizem×mwheremis the total number of nodes in our triangulation, including both interi or and boundary nodes. The rows and columns of K∗are labeled by the nodes xi. LetKν= (kν ij) be the corresponding m×mmatrixcontaining theelemental stiffnesses kν ijofTνintherowsandcolumnsindexed 12/11/12 848 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.16. The Oval Plate. by its vertices, and all other entries equal to 0. Thus, Kνwill have (at most) 9 nonzero entries. The resulting m×mmatrices are all summed together over all the triangles, K∗=N/summationdisplay ν=1Kν, (15.129) to produce the full finite element matrix, in accordance with (15.124). The full finite element matrix K∗is too large, since its rows and columns include all the nodes, whereas the finite element matrix Kappearing in (15.121) only refers to the n interior nodes. The reducedn×nfinite element matrix Kis simply obtained from K∗by deleting all rows and columns indexed by boundary nodes, ret aining only the elements kij when both xiandxjareinterior nodes. (Thismayremind thereader ofour constr uction of thereduced incidence matrixfor a structure inChapter 6.) F or thehomogeneous boundary value problem, this is all we require. As we shall see, inhomo geneous boundary conditions are most easily handled by retaining (part of) the full matri xK∗. The easiest way to digest the construction is by working thro ugh a particular example. Example 15.20. A metal plate has the shape of an oval running track, consisti ng of a rectangle, with side lengths 1m by 2m, and two semicircul ar disks glued onto its shorter ends, as sketched in Figure 15.16. The plate is subje ct to a heat source while its edges are held at a fixed temperature. The problem is to find the equilibrium temperature distribution within the plate. Mathematically, we must sol ve the Poisson equation with Dirichlet boundary conditions, for the equilibrium temper atureu(x,y). Let us describe how to set up the finite element approximation to such a boundary value problem. We begin with a very coarse triangulation of t he plate, which will not give particularly accurate results, but does serve to illustrat e how to go about assembling the finite element matrix. We divide the rectangular part of the p late into 8 right triangles, while each semicircular end will be approximated by three eq uilateral triangles. The tri- angles are numbered from 1 to 14 as indicated in Figure 15.17. There are 13 nodes in all, numbered as in the second figure. Only nodes 1 ,2,3 are interior, while the boundary nodes are labeled 4 through 13, going counterclockwise around the boundary starting at the top. 12/11/12 849 c/ci∇clecopy∇t2012 Peter J. Olver 1 2 34 567 8910 1112 13 14 Triangles1 2 34 5 6 7 8 9 10111213 Nodes Figure 15.17. A Coarse Triangulation of the Oval Plate. The full finite element matrix K∗will have size 13 ×13, its rows and columns labeled by all the nodes, while the reduced matrix Kappearing in the finite element equations (15.121) consists of the upper left 3 ×3 submatrix of K∗corresponding to the three interior nodes. Each triangle Tνwill contribute the summand Kνwhose values are its elemental stiffnesses, as indexed by its vertices. For example, the firs t triangleT1is equilateral, and so has elemental stiffnesses (15.128). Its vertices are labe led 1, 5, and 6, and therefore we place the stiffnesses (15.128) in the rows and columns numb ered 1,5,6 to form the summand K1= .577350 0 0 0 −.288675 −.288675 0 0 ... 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 ... −.288675 0 0 0 .577350 −.288675 0 0 ... −.288675 0 0 0 −.288675.577350 0 0 ... 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 ... ........................... , where all the undisplayed entries in the full 13 ×13 matrix are 0. The next triangle T2 has the same equilateral elemental stiffness matrix (15.128 ), but now its vertices are 1 ,6,7, and so it will contribute K2= .577350 0 0 0 0 −.288675 −.288675 0 ... 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 ... −.288675 0 0 0 0 .577350 −.2886750 0 ... −.288675 0 0 0 0 −.288675.5773500 0 ... 0 0 0 0 0 0 0 0 ... ........................... . 12/11/12 850 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.18. A Square Mesh for the Oval Plate. Similarly for K3, with vertices 1 ,7,8. On the other hand, triangle T4is an isoceles right triangle, and so has elemental stiffnesses (15.127). Its ver tices are labeled 1, 4, and 5, with vertex 5 at the right angle. Therefore, its contribution is K4= .5 0 0 0 −.5 0 0 0 ... 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 ... 0 0 0 .5−.5 0 0 0 ... −.5 0 0 −.5 1.0 0 0 0 ... 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 ... ........................... . Continuing in this manner, we assemble 14 contributions K1,...,K14, each with (at most) 9 nonzero entries. The full finite element matrix is the sum K∗=K1+K2+···+K14 = 3.732−1 0 0 −.7887−.5774−.5774 −1 4 −1−1 0 0 0 0−1 3.732 0 0 0 0 0−1 0 2 −.5 0 0 −.7887 0 0 −.5 1.577−.2887 0 −.5774 0 0 0 −.2887 1 .155−.2887 −.5774 0 0 0 0 −.2887 1 .155 −.7887 0 0 0 0 0 −.2887 0−1 0 0 0 0 0 0 0 −.7887 0 0 0 0 0 0 −.5774 0 0 0 0 0 0 −.5774 0 0 0 0 0 0 −.7887−.5 0 0 0(15.130) 12/11/12 851 c/circlecopyrt2012 Peter J. Olver −.7887 0 0 0 0 0 0−1 0 0 0 0 0 0 −.7887−.5774−.5774−.7887 0 0 0 0 0 −.5 0 0 0 0 0 0 0 0 0 0 0 0 −.2887 0 0 0 0 0 1.577−.5 0 0 0 0 −.5 2 −.5 0 0 0 0−.5 1.577−.2887 0 0 0 0 −.2887 1 .155−.2887 0 0 0 0 −.2887 1 .155−.2887 0 0 0 0 −.2887 1 .577 . Since only nodes 1 ,2,3 are interior nodes, the reduced finite element matrix only u ses the upper left 3 ×3 block ofK∗, so K= 3.732−1 0 −1 4 −1 0−1 3.732 . (15.131) It is not difficult to directly construct K, bypassing K∗entirely. For a finer triangulation, the construction is similar, but t he matrices become much larger. The procedure can, of course, be automated. Fortuna tely, if we choose a very regular triangulation, then we do not need to be nearly as met iculous in assembling the stiffness matrices, since many of the entries are the same. Th e simplest case is when we use a uniform square mesh, and so triangulate the domain into isoceles right triangles. This is accomplished by laying out a relatively dense square grid over the domain Ω ⊂R2. The interior nodes are the grid points that fall inside the ov al domain, while the boundary nodes are all those grid points lying adjacent to one or more o f the interior nodes, and are near but not necessarily precisely on the boundary ∂Ω. Figure 15.18 shows the nodes in a square grid with intermesh spacing h=.2. While a bit crude in its approximation of the boundary of the domain, this procedure does have the ad vantage of making the construction of the associated finite element matrix relati vely painless. For such a mesh, all the triangles are isoceles right triangl es, with elemental stiffnesses (15.127). Summing the corresponding matrices Kνover all the triangles, as in (15.129), the rows and columns of K∗corresponding to the interior nodes are seen to all have the same form. Namely, if ilabels an interior node, then the corresponding diagonal en try is kii= 4, while the off-diagonal entries kij=kji,i/\e}atio\slash=j, are equal to either −1 when node i is adjacent to node jon the grid, and is equal to 0 in all other cases. Node jis allowed to be a boundary node. (Interestingly, the result does not depe nd on how one orients the pair of triangles making up each square of the grid, which only pla ys a role in the computation of the right hand side of the finite element equation.) Observ e that the same computation applies even to our coarse triangulation. The interior node 2 belongs to all right isoceles triangles, and the corresponding entries in (15.130) are k22= 4, andk2j=−1 for the four adjacent nodes j= 1,3,4,9. Remark: Interestingly, the coefficient matrix arising from the finit e element method on a square (or even rectangular) grid is the same as the coeffic ient matrix arising from a 12/11/12 852 c/ci∇clecopy∇t2012 Peter J. Olver finitedifferencesolutiontotheLaplaceorPoissonequation , asdescribed inExercise . The finiteelement approach hastheadvantageofapplying tomuch moregeneral triangulations. In general, while the finite element matrix Kfor a two-dimensional boundary value problem is not as nice as the tridiagonal matrices we obtaine d in our one-dimensional problems, it is still very sparse and, on regular grids, high ly structured. This makes solution of the resulting linear system particularly amena ble to an iterative matrix solver such as Gauss–Seidel, Jacobi, or, for even faster convergen ce, successive over-relaxation (SOR). The Coefficient Vector and the Boundary Conditions So far, we have been concentrating on assembling the finite el ement coefficient matrix K. Wealsoneed tocomputetheforcingvector b= (b1,b2,...,bn)Tappearingontheright hand side of the fundamental linear equation (15.121). Acco rding to (15.119), the entries biare found by integrating the product of the forcing function and the finite element basis function. As before, we will approximate the integral over t he domain Ω by an integral over the triangles, and so bi=/integraldisplay/integraldisplay Ωfϕidxdy≈/summationdisplay ν/integraldisplay/integraldisplay Tνfων idxdy≡/summationdisplay νbν i. (15.132) Typically, the exact computation of the various triangular integrals is not convenient, and so we resort to a numerical approximation. Since we are as suming that the individual trianglesaresmall, wecanadopt averycrude numerical inte grationscheme. Ifthefunction f(x,y) does not vary much over the triangle Tν— which will certainly be the case if Tνis sufficiently small — we may approximate f(x,y)≈cν ifor (x,y)∈Tνby a constant. The integral (15.132) is then approximated by bν i=/integraldisplay/integraldisplay Tνfων idxdy≈cν i/integraldisplay/integraldisplay Tνων i(x,y)dxdy=1 3cν iareaTν=1 6cν i|∆ν|.(15.133) The formula for the integral of the affine element ων i(x,y) follows from solid geometry. Indeed, it equals the volume under its graph, a tetrahedron o f height 1 and base Tν, as illustrated in Figure 15.19. How to choose the constant cν i? In practice, the simplest choice is to let cν i=f(xi,yi) bethevalueofthefunctionatthe ithvertex. Withthischoice, thesum in(15.132)becomes bi≈/summationdisplay ν1 3f(xi,yi) areaTν=1 3f(xi,yi)areaPi, (15.134) wherePiis the vertex polygon (15.118) corresponding to the node xi. In particular, for the square mesh with the uniform choice of triangles, as in Ex ample 15.18, areaPi= 3h2for alli, and so bi≈f(xi,yi)h2(15.135) is well approximated by just h2times the value of the forcing function at the node. This is the underlying reason to choose the uniform triangulatio n for the square mesh; the alternating version would give unequal values for the biover adjacent nodes, and this would introduce unnecessary errors into the final approxima tion. 12/11/12 853 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.19. Finite Element Tetrahedron. Example 15.21. For the coarsely triangulated oval plate, the reduced stiffn ess ma- trix is (15.131). The Poisson equation −∆u= 4 models a constant external heat source of magnitude 4◦over the entire plate. If we keep theedges oftheplatefixed at 0◦, thenwe need to solvethefinite element equation Kc=b, whereKis the coefficient matrix (15.131), while b=4 3/parenleftig 2+3√ 3 4,2,2+3√ 3 4/parenrightigT = (4.39872,2.66667,4.39872)T. The entries of bare, by (15.134), equal to 4 = f(xi,yi) times one third the area of the corresponding vertex polygon, which for node 2 is the square consisting of 4 right triangles, each of area1 2, whereas for nodes 1 and 3 it consists of 4 right triangles of a rea1 2plus three equilateral triangles, each of area√ 3 4; see Figure 15.17. The solution to the final linear system is easily found: c= (1.56724,1.45028,1.56724)T. Its entries are the values of the finite element approximatio n at the three interior nodes. The finite element solution is plotted in the first illustrati on in Figure 15.20. A more accurate solution, based on a square grid triangulation of s izeh=.1 is plotted in the second figure. Inhomogeneous Boundary Conditions So far, we have restricted our attention to problems with hom ogeneous Dirichlet boundary conditions. According to Theorem 15.15, the solut ion to the inhomogeneous Dirichlet problem −∆u=fin Ω, u =hon∂Ω, is also obtained by minimizing the Dirichlet functional (15 .103). However, now the min- imization takes place over the affine subspace consisting of a ll functions that satisfy the 12/11/12 854 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.20. Finite Element Solutions to Poisson’s Equation for an Oval P late. inhomogeneous boundary conditions. It is not difficult to fit t his problem into the finite element scheme. Theelementscorresponding totheinteriornodesofourtria ngulationremainasbefore, but now we need to include additional elements to ensure that our approximation satisfies the boundary conditions. Note that if xkis a boundary node, then the corresponding boundary element ϕk(x,y) satisfies the interpolation condition (15.109), and so has the same piecewise affine form (15.117). The corresponding finite element approximation w(x,y) =m/summationdisplay i=1ciϕi(x,y), (15.136) has the same form as before, (15.110), but now the sum is over a ll nodes, both interior and boundary. As before, the coefficients ci=w(xi,yi)≈u(xi,yi) are the values of the finite element approximation at the nodes. Therefore, in ord er to satisfy the boundary conditions, we require cj=h(xj,yj) whenever xj= (xj,yj) is a boundary node .(15.137) Remark: If the boundary node xjdoes not lie precisely on the boundary ∂Ω, we need to approximate the value h(xj,yj) appropriately, e.g., by using the value of h(x,y) at the nearest boundary point ( x,y)∈∂Ω. The derivation of the finite element equations proceeds as be fore, but now there are additional terms arising from the nonzero boundary values. Leaving the intervening details to the reader, the final outcome can be written as follows. Let K∗denote the full m×m finite element matrix constructed as above. The reduced coeffi cient matrix Kis obtained by retaining the rows and columns corresponding to only inte rior nodes, and so will have sizen×n, wherenis the number of interior nodes. The boundary coefficient matrix /tildewideK is then×(m−n) matrix consisting of the entries of the interior rows that d o not appear inK, i.e., those lying in the columns indexed by the boundary nod es. For instance, in the the coarse triangulation of the oval plate, the full finite el ement matrix is given in (15.130), and the upper 3 ×3 subblock is the reduced matrix (15.131). The remaining ent ries of the 12/11/12 855 c/ci∇clecopy∇t2012 Peter J. Olver Figure 15.21. Solution to the Dirichlet Problem for the Oval Plate. first three rows form the boundary coefficient matrix /tildewideK= 0−.7887−.5774−.5774−.7887 0 0 0 0 0 −1 0 0 0 0 −1 0 0 0 0 0 0 0 0 0 0 −.7887−.5774−.5774−.7887 . (15.138) We similarly split the coefficients ciof the finite element function (15.136) into two groups. We letc∈Rndenote the as yet unknown coefficients cicorresponding to the values of the approximation at the interior nodes xi, whileh∈Rm−nwill be the vector of boundary values (15.137). The solution to the finite element approxim ation (15.136) is obtained by solving the associated linear system Kc+/tildewideKh=b,orKc=f=b−/tildewideKh. (15.139) Example 15.22. For the oval plate discussed in Example 15.20, suppose the ri ght handsemicircularedgeisheldat10◦, thelefthandsemicircularedgeat −10◦, whilethetwo straight edges have a linearly varying temperature distrib ution ranging from −10◦at the leftto10◦attheright, asillustratedinFigure15.21. Our taskistoco mputeitsequilibrium temperature, assuming no internal heat source. Thus, for th e coarse triangulation we have the boundary nodes values h= (h4,...,h13)T= (0,−1,−1,−1,−1,0,1,1,1,1,0)T. Using the previously computed formulae (15.131,138) for th e interior coefficient matrix K and boundary coefficient matrix /tildewideK, we approximate the solution to the Laplace equation by solving (15.139). We are assuming that there is no externa l forcing function, f(x,y)≡ 0, and so the right hand side is b=0, and so we must solve Kc=f=−/tildewideKh= (2.18564,3.6,7.64974)T. The finite element function corresponding to the solution c= (1.06795,1.8,2.53205)Tis plotted in the first illustration in Figure 15.21. Even on s uch a coarse mesh, the approximation is not too bad, as evidenced by the second illustration, which plots the finite element solution for a square mesh with spacingh=.2 between nodes. 12/11/12 856 c/ci∇clecopy∇t2012 Peter J. Olver