bvz

Chapter 11 BoundaryValueProblemsinOneDimension While its roots are firmly planted in the finite-dimensional w orld of matrices and vectors, the full scope of linear algebra is much broader. It s historical development and, hence, its structures, concepts, and methods, were strongl y influenced by linear analysis — specifically, the need to solve linear differential equation s, linear boundary value problems, linear integral equations, and the like. The time has come fo r us to fully transcend our finite dimensional limitations, and, in prepartion for late r developments, witness linear algebra in action in infinite-dimensional function spaces. In this chapter, we begin to analyze problems arising in cont inuum physics. The equilibrium equation of a one-dimensional continuum — an el astic bar, a bendable beam, and so on — is formulated as a boundary value problem for a scal ar ordinary differential equation. The framework introduced for discrete mechanica l systems in Chapter 6 will carry over, in its essence, to the infinite-dimensional sett ing appropriate to such problems. The underlying Euclidean vector space Rnis replaced by a function space. Vectors become functions, while matrices turn into linear differential ope rators. Physical boundary value problems are based on self-adjoint boundary value problems , founded on a suitable inner product on the function space. As in the discrete context, th e positive definite cases are stable, and the equilibrium solutioncan be characterized b y a minimizationprinciple based on a quadratic energy functional. Finite-dimensional linear algebra not only provides us wit h important insights into the underlying mathematical structure, but also motivates basic analytical and numerical solution schemes. In the function space framework, the gene ral superposition principle is reformulated in terms of the effect of a combination of impu lse forces concentrated at a single point of the continuum. However, constructing a f unction that represents a concentrated impulse turns out to be a highly non-trivial ma thematical issue. Ordinary functions do not suffice, and weareledtodevelop anew calculu sofgeneralized functions or distributions, including theremarkabledeltafunction. T heresponseofthesystemtoaunit impulse force is known as the Green’s function of the boundar y value problem, in honor of the self-taught English mathematician (and miller) Geor ge Green. With the Green’s functioninhand, thegeneralsolutiontotheinhomogeneous systemcanbereconstructedby superimposing the effects of suitably scaled impulses on the entire domain. Understanding thisconstructionwillbecomeincreasinglyimportantaswe progressontopartialdifferential equations, where direct analytical solution techniques ar e far harder to come by.We begin with second order boundary value problems describing the eq uilibria of stretchable bars. We continue on to fourth order boundary value problems that g overn the equilibrium of elastic beams, including piecewise cubic spline interpola nts that play a key role in modern 12/11/12 568 c/circlecopyrt2012 Peter J. Olver computer graphics and numerical analysis, and to more gener al second order boundary value problems of Sturm–Liouville type, which arise in a hos t of physical applications that involve partial differential equations. The simplest boundary value problems can be solved by direct integration. However, more complicated systems do not admit explicit formulae for their solutions, and one must rely on numerical approximations. In the final section, we introduce the powerful finite element method. The key idea is to restrict the infinite -dimensional minimization principle characterizing the exact solution to a suitably c hosen finite-dimensional subspace of the function space. When properly formulated, the soluti on to the resulting finite- dimensional minimization problem approximates the true mi nimizer. As in Chapter 4, the finite-dimensional minimizer is found by solving the induce d linear algebraic system, using either direct or iterative methods. An alternative formulation of the finite element solution, t hat can be applied even in situations where there is no minimum principle available, i s based on the idea of a weak solution to the boundary value problem, where one relaxes th e classical differentiability requirements. 11.1. Elastic Bars. Abaris a mathematical idealization of a one-dimensional linear ly elastic continuum that can be stretched or contracted in the longitudinal dire ction, but is not allowed to bend in a transverse direction. (Materials that can bend are call ed beams, and will be analyzed in Section 11.4.) We will view the bar as the continuum limit o f a one-dimensional chain of masses and springs, a system that we already analyzed in Se ction 6.1. Intuitively, the continuous bar consists of an infinite number of masses co nnected by infinitely short springs. The individual masses can be thought of as the “atom s” in the bar, although one should not try to read too much into the physics behind this in terpretation. We shall derive the basic equilibrium equations for the bar f rom first principles. Recall the three basic steps we already used to establish the corres ponding equilibrium equations for a discrete mechanical system such as a mass–spring chain : (i) First, use geometry to relate the displacement of the masse s to the elongation in the connecting springs. (ii) Second, use the constitutive assumptions such as Hooke’s L aw to relate the strain to the stress or internal force in the system. (iii) Finally, impose a force balance between external and inter nal forces. The remarkable fact, which will, when suitably formulated, carry over to the continuum, is that the force balance law is directly related to the geome trical displacement law by a transpose or, more correctly, adjoint operation. Consider a bar of length ℓhanging from a fixed support, with the bottom end left free, as illustrated in Figure 11.1. We use 0 ≤x≤ℓto refer to the reference or unstressed configuration of the bar, so xmeasures the distance along the bar from the fixed end x= 0 to the free end x=ℓ. Note that we are adopting the convention that the positive xaxis pointsdown. Letu(x) denote the displacement of the bar from its reference configuration. This means that the “atom” that started at position xhas moved to position x+u(x). 12/11/12 569 c/circlecopyrt2012 Peter J. Olver x u(x) Figure 11.1. Bar with One Fixed Support. With our convention, u(x)>0 means that the atom has moved down, while if u(x)<0, the atom has moved up. In particular, u(0) = 0 (11 .1) because we are assuming that the top end is fixed and cannot mov e. Thestrainin the bar measures the relative amount of stretching. Two ne arby atoms, atrespectivepositions xandx+∆x, aremovedtopositions x+u(x)andx+∆x+u(x+∆x). The original, unstressed length of this small section of bar was ∆x, while in the new configuration the same section has length /bracketleftbig x+∆x+u(x+∆x)/bracketrightbig −/bracketleftbig x+u(x)/bracketrightbig = ∆x+/bracketleftbig u(x+∆x)−u(x)/bracketrightbig . Therefore, this segment has been elongated by an amount u(x+ ∆x)−u(x). The di- mensionless strain measures the relative elongation, and s o is obtained by dividing by the reference length: [ u(x+∆x)−u(x)]/∆x. We now take the continuum limit by letting the interatomic spacing ∆ x→0. The result is the strain function v(x) = lim ∆x→0u(x+∆x)−u(x) ∆x=du dx(11.2) that measures the local stretch in the bar at position x. As noted above, we can approximate the bar by a chain of nmasses connected by n springs, letting the bottom mass hang free. The mass–spring chain should also have total lengthℓ, and so the individual springs have reference length ∆x=ℓ n. The bar is to be viewed as the continuum limit , in which the number of masses n→ ∞ and the spring lengths ∆ x→0. Thekthmass starts out at position xk=k∆x=kℓ n, 12/11/12 570 c/circlecopyrt2012 Peter J. Olver and, under forcing, experiences a displacement uk. The relative elongation of the kth spring†is vk=ek ∆x=uk+1−uk ∆x. (11.3) In particular, since the fixed end cannot move, the first value u0= 0 is omitted from the subsequent equations. The relation (11.3) between displacement and strain takes t he familiar matrix form v=Au, v=/parenleftbig v0,v1,...,vn−1/parenrightbigT, u= (u1,u2,...,un)T, where A=1 ∆x 1 −1 1 −1 1 −1 1 ...... −1 1 −→d dx(11.4) is thescaled incidence matrix of the mass–spring chain. As indicated by the arrow, the derivative operator d/dxthat relates displacement to strain in the bar equation (11. 2) can be viewed as its continuum limit, as the number of masses n→ ∞and the spring lengths ∆ x→0. Vice versa, the incidence matrix can be viewed as a discret e, numerical approximation to the derivative operator. Indeed, if we reg ard the discrete displacements and strains as approximations to the sample values of their c ontinuous counterparts, so uk≈u(xk), vk≈v(xk), then (11.3) takes the form v(xk) =u(xk+1)−u(xk) ∆x=u(xk+∆x)−u(xk) ∆x≈du dx(xk). justifying the identification (11.4). The passage back and f orth between the discrete and the continuous is the foundation of continuum mechanics — so lids, fluids, gases, and plas- mas. Discrete models both motivate and provide numerical ap proximations to continuum systems, which, in turn, simplify and give further insight i nto the discrete domain. The next part of the mathematical framework is to use the cons titutive relations to relate the strain to the stress, or internal force experienced by the bar. To keep matters simple, we shall only consider bars that have a linear relati on between stress and strain. For a physical bar, this is a pretty good assumption as long as it is not stretched beyond its elastic limits. Let w(x) denote the stress on the point of the bar that was at referenc e position x. Hooke’s Law implies that w(x) =c(x)v(x), (11.5) †We will find it helpful to label the springs from k= 0 tok=n−1. This will facilitate comparisons with the bar, which, by convention, starts at position x0= 0. 12/11/12 571 c/circlecopyrt2012 Peter J. Olver wherec(x) measures the stiffness of the bar at position x. For a homogeneous bar, made out of a uniform material, c(x)≡cis a constant function. The constitutive function c(x) can be viewed as the continuum limit of the diagonal matrix C= c0 c1 ... cn−1  of individual spring constants ckappearing in the discrete version wk=ckvk,orw=Cv, (11.6) that relates stress to strain (internal force to elongation ) in the individual springs. Indeed, (11.6) can be identified as the sampled version, w(xk) =c(xk)v(xk), of the continuum relation (11.5). Finally, we need to impose a force balance at each point of the bar. Suppose f(x) is an external force at position xon the bar, where f(x)>0 means the force is acting downwards. Physicalexamplesinclude mechanical, gravita tional,or magneticforcesacting solely in the vertical direction. In equilibrium†, the bar will deform so as to balance the external force with its own internal force resulting from st retching. Now, the internal force per unit length on the section of the bar lying between nearby positions xandx+ ∆x is the relative difference in stress at the two ends, namely [ w(x+∆x)−w(x)]/∆x. The force balance law requires that, in the limit, 0 =f(x)+ lim ∆x→0w(x+∆x)−w(x) ∆x=f(x)+dw dx, or f=−dw dx. (11.7) Thiscanbeviewedasthecontinuum limitofthe mass–spring c hainforce balanceequations fk=wk−1−wk ∆x, wn= 0, (11.8) where the final condition ensures the correct formula for the force on the free-hanging bottom mass. (Remember that the springs are numbered from 0 t on−1.) This indicates that we should also impose an analogous boundary condition w(ℓ) = 0 (11 .9) at the bottom end of the bar, which is hanging freely and so is u nable to support any internal stress. The matrix form of the discrete system (11. 8) is f=ATw, †The dynamical processes leading to equilibrium will be discuss ed in Chapter 14. 12/11/12 572 c/circlecopyrt2012 Peter J. Olver where the transposed scaled incidence matrix AT=1 ∆x 1−1 1−1 1−1 1−1 ...... −→ −d dx(11.10) should approximate the differential operator −d/dxthat appears in the continuum force balance law (11.7). Thus, we should somehow interpret −d/dxas the “transpose” or “adjoint”ofthedifferentialoperator d/dx. Thisimportantpointwillbedevelopedproperly in Section 11.3. But before trying to push the theory any furt her, we will pause to analyze the mathematical equations governing some simple configura tions. But first, let us summarize our progress so far. The three basi c equilibrium equations (11.2,5,7) are v(x) =du dx, w (x) =c(x)v(x), f (x) =−dw dx.(11.11) Substituting the first into the second, and then the resultin g formula into the last equation, leads to the equilibrium equation K[u] =−d dx/parenleftbigg c(x)du dx/parenrightbigg =f(x), 0< x < ℓ. (11.12) Thus, the displacement u(x) of the bar is obtained as the solution to a second order ordinary differential equation. As such, it will depend on tw o arbitrary constants, which will be uniquely determined by the boundary conditions†(11.1,9) at the two ends: u(0) = 0, w (ℓ) =c(ℓ)u′(ℓ) = 0. (11.13) Usuallyc(ℓ)>0, in which case it can be omitted from the second boundary con dition, which simply becomes u′(ℓ) = 0. The resulting boundary value problem is to be viewed as the continuum limit of the linear system Ku=ATCAu=f (11.14) modeling a mass-spring chain with one free end, cf. (6.11), i n which A−→d dx, C−→c(x), AT−→ −d dx,u−→u(x),f−→f(x). And, as we will see, most features of the finite-dimensional l inear algebraic system have, when suitably interpreted, direct counterparts in the cont inuous boundary value problem. †We will sometimes use primes, as in u′=du/dx, to denote derivatives with respect to x. 12/11/12 573 c/circlecopyrt2012 Peter J. Olver 0.2 0.4 0.6 0.8 10.20.40.60.81 u(x)0.2 0.4 0.6 0.8 10.20.40.60.81 w(x) Figure 11.2. Displacement and Stress of a Bar with One Fixed End. Example 11.1. Consider the simplest case of a uniform bar of unit length ℓ= 1 subjected to a uniform force, e.g., gravity. The equilibriu m equation (11.12) is −cd2u dx2=f, (11.15) where we are assuming that the force fis constant. This elementary second order ordinary differential equation can be immediately integrated: u(x) =−1 2αx2+ax+b, where α=f c(11.16) is the ratio of the force to the stiffness of the bar. The values of the integration constants aandbare fixed by the boundary conditions (11.13), so u(0) =b= 0, u′(1) =−α+a= 0. Therefore, there is a unique solution to the boundary value p roblem, yielding the displace- ment u(x) =α/parenleftbig x−1 2x2/parenrightbig , (11.17) which is graphed in Figure 11.2 for α= 1. Note the parabolic shape, with zero derivative, indicating no strain, at the free end. The displacement reac hes its maximum, u(1) =1 2α, at the free end of the bar, which is the point which moves downw ards the farthest. The stronger the force or the weaker the bar, the farther the over all displacement. Remark: This example illustrates the simplest way to solve boundar y value problems, which is adapted from the usual solution technique for initi al value problems. First, solve the differential equation by standard methods (if possible) . For a second order equation, the general solution will involve two arbitrary constants. The values of the constants are found by substituting the solution formula into the two boun dary conditions. Unlikeinitial value problems, the existence and/or uniqueness of the solu tion to a general boundary value problem is not guaranteed, and you may encounter situa tions where you are unable to complete the solution; see, for instance, Example 7.42. A more sophisticated method, based on the Green’s function, will be presented in the follo wing section. 12/11/12 574 c/circlecopyrt2012 Peter J. Olver Asinthe discrete situation, thisparticular mechanical co nfigurationis statically deter- minate, meaning that we cansolvedirectly for the stress w(x) intermsof theexternal force f(x) without having to compute the displacement u(x) first. In this particular example, we need to solve the first order boundary value problem −dw dx=f, w (1) = 0, arising from the force balance law (11.7). Since fis constant, w(x) =f(1−x),and v(x) =w(x) c=α(1−x). Note that the boundary condition uniquely determines the in tegration constant. We can then find the displacement u(x) by solving another boundary value problem du dx=v(x) =α(1−x), u (0) = 0, resulting from (11.2), which again leads to (11.17). As befo re, the appearance of one boundary condition implies that we can find a unique solution to the differential equation. Remark: We motivated the boundary value problem for the bar by takin g the contin- uum limit of the mass–spring chain. Let us see to what extent t his limiting procedure can be justified. To compare the solutions, we keep the reference length of the chain fixed at ℓ= 1. So, if we have nidentical masses, each spring has length ∆ x= 1/n. Thekthmass will start out at reference position xk=k/n. Using static determinacy, we can solve the system (11.8), which reads wk−1=wk+f n, wn= 0, directly for the stresses: wk=f/parenleftbigg 1−k n/parenrightbigg =f(1−xk), k = 0,...,n−1. Thus, in this particular case, the continuous bar and the dis crete chain have equal stresses at the sample points: w(xk) =wk. The strains are also in agreement: vk=1 cwk=α/parenleftbigg 1−k n/parenrightbigg =α(1−xk) =v(xk), whereα=f/c, as before. We then obtain the displacements by solving uk+1=uk+vk n=uk+α n/parenleftbigg 1−k n/parenrightbigg . Sinceu0= 0, the solution is uk=α nk−1/summationdisplay i=0/parenleftbigg 1−i n/parenrightbigg =α/parenleftbiggk n−k(k−1) 2n2/parenrightbigg =α/parenleftbig xk−1 2x2 k/parenrightbig +αxk 2n=u(xk)+αxk 2n. (11.18) 12/11/12 575 c/circlecopyrt2012 Peter J. Olver 0.2 0.4 0.6 0.8 10.050.10.150.20.250.3 Figure 11.3. Displacements of a Bar with Two Fixed Ends. The sampled displacement u(xk) is not exactly equal to uk, but their difference tends to zero as the number of masses n→ ∞. In this way, we have completely justified our limiting interpretation. Example 11.2. Consider the same uniform, unit length bar as in the previous example, again subject to a uniform constant force, but now w ith two fixed ends. We impose inhomogeneous boundary conditions u(0) = 0, u (1) =d, so the top end is fixed, while the bottom end is displaced an amo untd. (Note that d >0 means the bar is stretched, while d <0 means it is compressed.) The general solution to the equilibrium equation (11.15) is, as before, given by (11 .16). The values of the arbitrary constants a,bare determined by plugging into the boundary conditions, so u(0) =b= 0, u (1) =−1 2α+d= 0. Thus u(x) =1 2α(x−x2)+dx (11.19) is the unique solution to the boundary value problem. The dis placement is a linear su- perposition of two functions; the first is induced by the exte rnal force f, while the second represents a uniform stretch induced by the boundary condit ion. In Figure 11.3, the dotted curves represent the two constituents, and the solid graph i s their sum, which is the actual displacement. Unlike a bar with a free end, this configuration is statically indeterminate . There is no boundary condition on the force balance equation −dw dx=f, andso theintegrationconstant ainthestress w(x) =a−fxcannotbedetermined without first figuring out the displacement (11.19): w(x) =cdu dx=f/parenleftbig1 2−x/parenrightbig +cd. 12/11/12 576 c/circlecopyrt2012 Peter J. Olver Example 11.3. Finally, consider the case when both ends of the bar are left f ree. The boundary value problem −u′′=f(x), u′(0) = 0, u′(ℓ) = 0, (11.20) represents the continuum limit of a mass–spring chain with t wo free ends and corresponds to a bar floating in outer space, subject to a nonconstant exte rnal force. Based on our finite-dimensional experience, we expect the solution to ma nifest an underlying instability of the physical problem. Solving the differential equation, we find that u(x) =ax+b−/integraldisplayx 0/parenleftbigg/integraldisplayy 0f(z)dz/parenrightbigg dy, where the constants a,bare to be determined by the boundary conditions. Since u′(x) =a−/integraldisplayx 0f(z)dz, the first boundary condition u′(0) = 0 requires a= 0. The second boundary condition requires u′(ℓ) =−/integraldisplayℓ 0f(x)dx= 0, (11.21) which is not automatically valid! The integral represents t he total force per unit length exerted on the bar. As in the case of a mass-spring chain with t wo free ends, if there is a non-zero net force, the bar cannot remain in equilibrium, bu t will move off in space and the equilibrium boundary value problem has no solution. On t he other hand, if the forcing satisfies the constraint (11.21), then the resulting soluti on of the boundary value problem has the form u(x) =b−/integraldisplayx 0/parenleftbigg/integraldisplayy 0f(z)dz/parenrightbigg dy, (11.22) where the constant bis arbitrary. Thus, when it exists, the solution to the bound ary value problem is not unique. The constant bsolves the corresponding homogeneous problem, and represents a rigid translation of the entire bar by a dist anceb. Physically, the free boundary value problem corresponds to an unstable structure: there is a translational instability in which the bar moves o ff rigidly in the longitudinal direction. Only balanced forces of mean zero can maintain eq uilibrium. Furthermore, when it does exist, the equilibrium solution is not unique si nce there is nothing to tie the bar down to any particular spatial position. This dichotomy should remind you of our earlier study of line ar algebraic systems. An inhomogeneous system Ku=fconsisting of nequations in nunknowns either admits a unique solution for all possible right hand sides f, or, when Kis singular, either no solution exists or the solution is not unique. In the latter c ase, the constraints on the right hand side are prescribed by the Fredholm alternative ( 5.80), which requires that f be orthogonal, with respect to the Euclidean inner product, to all elements of coker K. In physical equilibrium systems Kis symmetric, and so coker K= kerK. Thus, the Fredholm alternative requires that the forcing be orthogon al to all the unstable modes. In 12/11/12 577 c/circlecopyrt2012 Peter J. Olver the function space for the bar, the finite-dimensional dot pr oduct is replaced by the L2 inner product /an}bracketle{tf;g/an}bracketri}ht=/integraldisplayℓ 0f(x)g(x)dx. Since the kernel or solution space to the homogeneous bounda ry value problem is spanned by the constant function 1, the Fredholm alternative†requires that the forcing function be orthogonal to it, /an}bracketle{tf;1/an}bracketri}ht=/integraldisplayℓ 0f(x)dx= 0. This is precisely the condition (11.21) required for existence of a (non-unique) equilibrium solution, and s o the analogy between the finite and infinite dimensional categories is complete. Remark: The boundary value problems that govern the mechanical equ ilibria of a simple bar arise in many other physical systems. For example , the equation for the thermal equilibrium of a bar under an external heat source is modeled by the same boundary value problem (11.12); in this case, u(x) represents the temperature of the bar, c(x) represents thediffusivity orthermal conductivity of the material at position x, whilef(x) represents an external heat source. A fixed boundary condition u(ℓ) =acorresponds to an end that is held at a fixed temperature a, while a free boundary condition u′(ℓ) = 0 represents an insulated end that does not allow heat energy to enter or leav e the bar. Details of the physical derivation can be found in Section 14.1. 11.2. Generalized Functions and the Green’s Function. The general superposition principle for inhomogeneous lin ear systems inspires an al- ternative, powerful approach to the solution of boundary va lue problems. This method relies on the solution to a particular type of inhomogeneity , namely a concentrated unit impulse. The resulting solutions are collectively known as the Green’s function for the boundary value problem. Once the Green’s function is known, the response of the sys- tem to any other external forcing can be constructed through a continuous superposition of these fundamental solutions. However, rigorously formu lating a concentrated impulse force turns out to be a serious mathematical challenge. To motivate the construction, let us return briefly to the cas e of a mass–spring chain. Given the equilibrium equations Ku=f, (11.23) let us decompose the external forcing f= (f1,f2,...,fn)T∈Rninto a linear combination f=f1e1+f2e2+···+fnen (11.24) of the standard basis vectors of Rn. Suppose we know how to solve each of the individual systems Kui=ei, i = 1,...,n. (11.25) †In fact, historically, Fredholm first made his discovery while stu dying such infinite-dimen- sional systems. Only later was it realized that the same underlying i deas are equally valid in finite-dimensional linear algebraic systems. 12/11/12 578 c/circlecopyrt2012 Peter J. Olver The vector eirepresents a unit force or, more precisely, a unit impulse , which is applied solely to the ithmass in the chain; the solution uirepresents the response of the chain. Since we can decompose any other force vector as a superposit ion of impulse forces, as in (11.24), the superposition principle tells us that the solu tion to the inhomogeneous system (11.23) is the self-same linear combination of the individu al responses, so u=f1u1+f2u2+···+fnun. (11.26) Remark: The alert reader willrecognize that u1,...,unare the columns of the inverse matrix,K−1, and so weare, infact, reconstructing the solutiontothe li nearsystem (11.23) by inverting the coefficient matrix K. Thus, this observation, while noteworthy, does not lead to an efficient solution technique for discrete systems. In contrast, in the case of continuous boundary value problems, this approach leads to one of the most valuable solution paradigms in both practice and theory. The Delta Function Our aim is to extend this algebraic technique to boundary val ue problems. The key question is how to characterize an impulse force that is conc entrated on a single atom†of the bar. In general, a unit impulse at position x=ywill be described by something called thedelta function , and denoted by δy(x). Since the impulse is supposed to be concentrated solely at x=y, we should have δy(x) = 0 for x/ne}ationslash=y. (11.27) Moreover, since it is a unitimpulse, we want the total amount of force exerted on the bar to be equal to one. Since we are dealing with a continuum, the t otal force is represented by an integral over the length of the bar, and so we also requir e that the delta function satisfy/integraldisplayℓ 0δy(x)dx= 1,provided that 0 < y < ℓ. (11.28) Alas, there is no bona fide function that enjoys both of the required properties! Indee d, according to the basic facts of Riemann (or even Lebesgue) in tegration, two functions which are the same everywhere except at one single point have exactly the same integral, [53,158]. Thus, since δyis zero except at one point, its integral should be 0, not 1. Th e mathematical conclusion is that the two requirements, (11. 27,28) are inconsistent! This unfortunate fact stopped mathematicians dead in their tracks. It took the imagi- nationofa Britishengineer, withthe unlikely nameOliver H eaviside, who was not deterred by the lack of rigorous justification, to start utilizing del ta functions in practical applica- tions — with remarkable effect. Despite his success, Heavisi de was ridiculed by the pure mathematicians of his day, and eventually succumbed to ment al illness. But, some thirty years later, the great theoretical physicist Paul Dirac res urrected the delta function for †As before, “atom” is used in a figurative sense. 12/11/12 579 c/circlecopyrt2012 Peter J. Olver quantum mechanical applications, and this finally made theo reticians sit up and take no- tice. (Indeed, the term “Dirac delta function” is quite comm on.) In 1944, the French mathematician Laurent Schwartz finally established a rigor ous theory of distributions that incorporated such useful, but non-standard objects, [ 116,156]. (Thus, to be more ac- curate, we should really refer to the delta distribution ; however, we will retain the more common, intuitive designation “delta function” throughou t.) It is beyond the scope of this introductory text to develop a fully rigorous theory of dist ributions. Rather, in the spirit of Heaviside, we shall concentrate on learning, through pra ctice with computations and applications, how to tame these wild mathematical beasts. There are two distinct ways to introduce the delta function. Both are important and both worth knowing. Method #1. Limits : The first approach is to regard the delta function δy(x) as a limit of a sequence of ordinary smooth functions†gn(x). These functions will represent more and more concentrated unit forces, which, in the limit, converge to the desired unit impulse concentrated at a single point, x=y. Thus, we require lim n→∞gn(x) = 0, x /ne}ationslash=y, (11.29) while the total amount of force remains fixed at /integraldisplayℓ 0gn(x)dx= 1. (11.30) On a formal level, the limit “function” δy(x) = lim n→∞gn(x) will satisfy the key properties (11.27–28). An explicit example of such a sequence is provided by the rati onal functions gn(x) =n π(1+n2x2). (11.31) These functions satisfy lim n→∞gn(x) =/braceleftbigg0, x/ne}ationslash= 0, ∞, x= 0,(11.32) while‡/integraldisplay∞ −∞gn(x)dx=1 πtan−1nx/vextendsingle/vextendsingle/vextendsingle/vextendsingle∞ x=−∞= 1. (11.33) †To keep the notation compact, we suppress the dependence of the func tionsgnon the point ywhere the limiting delta function is concentrated. ‡For the moment, it will be slightly simpler here to consider the en tire real line — correspond- ing to a bar of infinite length. Exercise discusses how to modify the construction for a finite interval. 12/11/12 580 c/circlecopyrt2012 Peter J. Olver Figure 11.4. Delta Function as Limit. Therefore, formally, we identify the limiting function lim n→∞gn(x) =δ(x) =δ0(x), (11.34) withtheunitimpulsedeltafunctionconcentratedat x= 0. AssketchedinFigure11.4, as n gets larger and larger, each successive function gn(x) forms a more and more concentrated spike, while maintaining a unit total area under its graph. T he limiting delta function can be thought of as an infinitely tall spike of zero width, entire ly concentrated at the origin. Remark: There are many other possible choices for the limiting func tionsgn(x). See Exercise for another useful example. Remark: This construction of the delta function highlights the per ils of interchanging limits and integrals without proper justification. In any st andard theory of integration, the limit of the functions gnwould be indistinguishable from the zero function, so the li mit of their integrals (11.33) would notequal the integral of their limit: 1 = lim n→∞/integraldisplay∞ −∞gn(x)dx/ne}ationslash=/integraldisplay∞ −∞lim n→∞gn(x)dx= 0. The delta function is, in a sense, a means of sidestepping thi s analytic inconvenience. The full ramifications and theoretical constructions underlyi ng such limits must, however, be deferred to a rigorous course in real analysis, [ 53,158]. Once we have found the basic delta function δ(x) =δ0(x), which is concentrated at the origin, we can obtain a delta function concentrated at an y other position yby a simple translation: δy(x) =δ(x−y). (11.35) 12/11/12 581 c/circlecopyrt2012 Peter J. Olver Thus,δy(x) can be realized as the limit of the translated functions /hatwidegn(x) =gn(x−y) =n π/parenleftbig 1+n2(x−y)2/parenrightbig. (11.36) Method #2. Duality : The second approach is a bit more abstract, but much closer to the proper rigorous formulation of the theory of distribu tions like the delta function. The critical observation is that if u(x) is any continuous function, then /integraldisplayℓ 0δy(x)u(x)dx=u(y),for 0 < y < ℓ. (11.37) Indeed, since δy(x) = 0 for x/ne}ationslash=y, the integrand only depends on the value of uat the pointx=y, and so /integraldisplayℓ 0δy(x)u(x)dx=/integraldisplayℓ 0δy(x)u(y)dx=u(y)/integraldisplayℓ 0δy(x)dx=u(y). Equation (11.37) serves to define a linear functional†Ly:C0[0,ℓ]→Rthat maps a con- tinuous function u∈C0[0,ℓ] to its value at the point x=y: Ly[u] =u(y)∈R. In the dual approach to generalized functions, the delta fun ction is, in fact, definedas this particular linear functional. The function u(x) is sometimes referred to as a test function , since it serves to “test” the form of the linear functional L. Remark: If the impulse point ylies outside the integration domain, then /integraldisplayℓ 0δy(x)u(x)dx= 0,when y <0 or y > ℓ, (11.38) because the integrand is identically zero on the entire inte rval. For technical reasons, we will not attempt to define the integral (11.38) if the impulse pointy= 0 ory=ℓlies on the boundary of the interval of integration. The interpretation of the linear functional Lyas representing a kind of function δy(x) is based on the following line of thought. According to Theor em 7.10, every scalar-valued linear function L:V→Ron a finite-dimensional inner product space is given by an inn er product with a fixed element a∈V, so L[u] =/an}bracketle{ta;u/an}bracketri}ht. †Linearity, which requires that Ly[cf+dg] =cLy[f]+dLy[g] for all functions f,gand all scalars (constants) c,d∈R, is easily established; see also Example 7.7. 12/11/12 582 c/circlecopyrt2012 Peter J. Olver In this sense, linear functions on Rnare the “same” as vectors. (But bear in mind that the identification does depend upon the choice of inner product. ) Similarly, on the infinite- dimensional function space C0[0,ℓ], the L2inner product Lg[u] =/an}bracketle{tg;u/an}bracketri}ht=/integraldisplayℓ 0g(x)u(x)dx (11.39) takenwithafixedfunction g∈C0[0,ℓ]definesareal-valuedlinearfunctional Lg:C0[0,ℓ]→ R. However, unlike the finite-dimensional situation, notevery real-valued linear functional has this form! In particular, there is no actual function δy(x) such that the identity /an}bracketle{tδy;u/an}bracketri}ht=/integraldisplayℓ 0δy(x)u(x)dx=u(y) (11 .40) holds for every continuous function u(x). Every (continuous) function defines a linear functional, but not conversely. Or, stated another way, whi le the dual space to a finite- dimensional vector space like Rncan be identified, via an inner product, with the space itself, this is not the case in infinite-dimensional functio n space; the dual is an entirely different creature. Thisdisconcerting facthighlightsyet anotheroftheprofounddifferences between finite- and infinite-dimensional vector spaces! But the dual interpretation of generalized functions acts a s if this were true. Gener- alized functions are real-valued linear functionals on fun ction space, but viewed as a kind of function via the inner product . Although the identification is not to be taken literally, one can, with a little care, manipulate generalized functio ns as if they were actual func- tions, but always keeping in mind that a rigorous justificati on of such computations must ultimately rely on their true characterization as linear fu nctionals. The two approaches — limits and duality — are completely comp atible. Indeed, with a little extra work, one can justify the dual formula (11.37) as the limit u(y) = lim n→∞/integraldisplayℓ 0gn(x)u(x)dx=/integraldisplayℓ 0δy(x)u(x)dx (11.41) of the inner products of the function uwith the approximating concentrated impulse functions gn(x) satisfying (11.29–30). In this manner, the linear functio nalLy[u] =u(y) represented by the delta function is the limit, Ly= lim n→∞Ln, of the approximating linear functionals Ln[u] =/integraldisplayℓ 0gn(x)u(x)dx. Thus, the choice of interpretation of the generalized delta function is, on an operational level, a matter of taste. For the novice, the limit interpret ation of the delta function is perhaps the easier to digest at first. However, the dual, line ar functional interpretation has stronger connections with the rigorous theory and, even in applications, offers some significant advantages. Although on the surface, the delta function might look a litt le bizarre, its utility in modern applied mathematics and mathematical physics more t han justifies including it in 12/11/12 583 c/circlecopyrt2012 Peter J. Olver Figure 11.5. The Step Function. your analytical toolbox. Even though you are probably not ye t comfortable with either definition, you are advised to press on and familiarize yours elf with its basic properties, to be discussed next. With a little care, you usually won’t go fa r wrong by treating it as if it were a genuine function. After you gain more practical exp erience, you can, if desired, return to contemplate just exactly what kind of object the de lta function really is. Remark: If you are familiar with basic measure theory, [ 158], there is yet a third interpretation of the delta function as a point mass or atomi c measure. However, the measure-theoretic approach has definite limitations, and d oes not cover the full gamut of generalized functions. Calculus of Generalized Functions In order to develop a working relationship with the delta fun ction, we need to under- stand how it behaves under the basic operations of linear alg ebra and calculus. First, we can take linear combinations of delta functions. For exampl e, f(x) = 2δ(x)+3δ(x−1) represents a combination of an impulse of magnitude 2 concen trated at x= 0 and one of magnitude 3 concentrated at x= 1. Since δy(x) = 0 for any x/ne}ationslash=y, multiplying the delta function by an ordinary function is the same as multiplying b y a constant: g(x)δy(x) =g(y)δy(x), (11.42) provided that g(x) is continuous at x=y. For example, xδ(x)≡0 is the same as the constant zero function. Warning : It isnotpermissible to multiply delta functions together, or to use more complicated algebraic operations. Expressions like δ(x)2, 1/δ(x),eδ(x), etc., are notwell defined in the theory of generalized functions. This makes th eir application to nonlinear systems much more problematic . 12/11/12 584 c/circlecopyrt2012 Peter J. Olver Figure 11.6. Step Function as Limit. The integral of the delta function is known as a step function . More specifically, the basic formulae (11.37,38) imply that /integraldisplayx aδy(t)dt=σy(x) =σ(x−y) =/braceleftbigg0, a < x < y, 1, x > y > a.(11.43) Figure 11.5 shows the graph of σ(x) =σ0(x). Unlike the delta function, the step function σy(x) is an ordinary function. It is continuous — indeed constant — except at x=y. The value of the step function at the discontinuity x=yis left unspecified, although a popular choice, motivated by Fourier theory, cf. Chapter 12 , is to se tσy(y) =1 2, the average of its left and right hand limits. We note that the integration formula (11.43) is compatible w ith our characteriza- tion of the delta function as the limit of highly concentrate d forces. If we integrate the approximating functions (11.31), we obtain fn(x) =/integraldisplayx −∞gn(t)dt=1 πtan−1nx+1 2. Since lim y→∞tan−1y=1 2π,while lim y→−∞tan−1y=−1 2π, these functions converge to the step function: lim n→∞fn(x) =σ(x) =  1, x > 0, 1 2, x= 0, 0, x < 0.(11.44) A graphical illustration of this limiting process appears i n Figure 11.6. Theintegral ofthediscontinuous stepfunction (11.43)ist hecontinuous ramp function /integraldisplayx aσy(z)dz=ρy(x) =ρ(x−y) =/braceleftbigg0, a < x < y, x−y, x > y > a,(11.45) which is graphed in Figure 11.7. Note that ρ(x−y) has a corner at x=y, and so is not differentiable there; indeed, its derivativedρ dx=σhas a jump discontinuity, and its second 12/11/12 585 c/circlecopyrt2012 Peter J. Olver -1.5 -1 -0.5 0.5 1 1.5 -0.20.20.40.60.811.2 -1.5 -1 -0.5 0.5 1 1.5 -0.20.20.40.60.811.2 Figure 11.7. First and Second Order Ramp Functions. derivatived2ρ dx2=δis no longer an ordinary function. We can continue to integra te; the nthintegral of the delta function is the nthorder ramp function ρn(x−y) =  (x−y)n n!, x > y, 0, x < y.(11.46) What about differentiation? Motivated by the Fundamental Th eorem of Calculus, we shall use formula (11.43) to identify the derivative of th e step function with the delta functiondσ dx=δ. (11.47) This fact is highly significant. In basic calculus, one is not allowed to differentiate a discontinuous function. Here, we discover that the derivat ive can be defined, not as an ordinary function, but rather as a generalized delta functi on. This basicidentity isa particular instanceof a general rul efor differentiating functions with discontinuities. We use f(y−) = lim x→y−f(x), f (y+) = lim x→y+f(x), (11.48) to denote, respectively, the left and right sided limits of a function at a point y. The function f(x) iscontinuous at the point yif and only if its one-sided limits exist and are equal to its value: f(y) =f(y−) =f(y+). If the one-sided limits are the same, but not equal to f(y), then the function is said to have a removable discontinuity , since redefining f(y) =f(y−) =f(y+) serves to make fcontinuous at the point in question. An example is the function f(x) that is equal to 0 for all x/ne}ationslash= 0, but has†f(0) = 1. Removing the discontinuity by setting f(0) = 0 makes f(x)≡0 equal to a continuous constant function. Since removable discontinuities play no role in our theory o r applications, they will always be removed without penalty. Warning : Although δ(0+) = 0 =δ(0−), we will emphatically notcall 0 a removable discontinuityofthedeltafunction. Onlystandardfunctio nshaveremovablediscontinuities. †This function is nota version of the delta function. It is an ordinary function, and its int egral is 0, not 1. 12/11/12 586 c/circlecopyrt2012 Peter J. Olver -1-0.5 0.5 11.5 2 -1-0.50.51 f(x)-1-0.5 0.5 11.5 2 -1-0.50.51 f′(x) Figure 11.8. The Derivative of a Discontinuous Function. Finally, if both the left and right limits exist, but are not e qual, then fis said to have ajump discontinuity at the point y. Themagnitude of the jump is the difference β=f(y+)−f(y−) = lim x→y+f(x)−lim x→y−f(x) (11 .49) between the right and left limits. The magnitude of the jump i s positive if the function jumps up, when moving from left to right, and negative if it ju mps down. For example, the step function σ(x) has a unit, i.e., magnitude 1, jump discontinuity at the ori gin: σ(0+)−σ(0−) = 1−0 = 1, and is continuous everywhere else. Note the value of the func tion at the point, namely f(y), which may not even be defined, plays no role in the specificat ion of the jump. In general, the derivative of a function with jump discontin uities is a generalized function that includes delta functions concentrated at eac h discontinuity. More explicitly, suppose that f(x) is differentiable, in the usual calculus sense, everywhere except at the pointywhere it has a jump discontinuity of magnitude β. We can re-express the function in the convenient form f(x) =g(x)+βσ(x−y), (11.50) whereg(x) is continuous everywhere, with a removable discontinuity atx=y, and differ- entiable except possibly at the jump. Differentiating (11.5 0), we find that f′(x) =g′(x)+βδ(x−y), (11.51) has a delta spike of magnitude βat the discontinuity. Thus, the derivatives of fandg coincide everywhere except at the discontinuity. Example 11.4. Consider the function f(x) =/braceleftbigg−x, x < 1, 1 5x2, x > 1,(11.52) whichwegraphinFigure11.8. Wenotethat fhasasinglejumpdiscontinuityofmagnitude 6 5atx= 1. This means that f(x) =g(x)+6 5σ(x−1),where g(x) =/braceleftbigg−x, x < 1, 1 5x2−6 5, x > 1, 12/11/12 587 c/circlecopyrt2012 Peter J. Olver -1 -0.5 0.5 1 1.5 2 -1-0.50.51 f(x)-1 -0.5 0.5 1 1.5 2 -4-224 f′(x) Figure 11.9. The Derivative of a Discontinuous Function. is continuous everywhere, since its right and left hand limi ts at the original discontinuity are equal: g(1+) =g(1−) =−1. Therefore, f′(x) =g′(x)+6 5δ(x−1),where g′(x) =/braceleftbigg−1, x < 1, 2 5x, x > 1, whileg′(1), and hence f′(1), is not defined. In Figure 11.8, the delta spike in the deri vative offis symbolized by a vertical line — although this pictorial de vice fails to indicate its magnitude of6 5. Sinceg′(x) can be found by directly differentiating the formula for f(x), once we determine the magnitude and location of the jump discontinu ities off(x), we can compute its derivative directly without introducing to the auxilia ry function g(x). Example 11.5. As a second, more streamlined example, consider the functio n f(x) =  −x, x < 0, x2−1,0< x <1, 2e−x, x > 1, which is plotted in Figure 11.9. This function has jump disco ntinuities of magnitude −1 atx= 0, and of magnitude 2 /eatx= 1. Therefore, in light of the preceding remark, f′(x) =−δ(x)+2 eδ(x−1)+  −1, x < 0, 2x, 0< x <1, −2e−x, x > 1, where the final terms are obtained by directly differentiatin gf(x). Example 11.6. The derivative of the absolute value function a(x) =|x|=/braceleftbiggx, x > 0, −x, x < 0, 12/11/12 588 c/circlecopyrt2012 Peter J. Olver Figure 11.10. Derivative of Delta Function as Limit of Doublets. is thesign function s(x) =a′(x) =/braceleftbigg+1, x > 0, −1, x < 0.(11.53) Note that there is no delta function in a′(x) because a(x) is continuous everywhere. Since s(x) has a jump of magnitude 2 at the origin and is otherwise const ant, its derivative s′(x) =a′′(x) = 2δ(x) is twice the delta function. We are even allowed to differentiate the delta function. Its fi rst derivative δ′ y(x) =δ′(x−y) can be interpreted in two ways. First, we may view δ′(x) as the limit of the derivatives of the approximating functions (11.31): dδ dx= lim n→∞dgn dx= lim n→∞−2n3x π(1+n2x2)2. (11.54) The graphs of these rational functions take the form of more a nd more concentrated spiked “doublets”, as illustrated in Figure 11.10. To determine th e effect of the derivative on a test function u(x), we compute the limiting integral /an}bracketle{tδ′;u/an}bracketri}ht=/integraldisplay∞ −∞δ′(x)u(x)dx= lim n→∞/integraldisplay∞ −∞g′ n(x)u(x)dx =−lim n→∞/integraldisplay∞ −∞gn(x)u′(x)dx=−/integraldisplay∞ −∞δ(x)u′(x)dx=−u′(0).(11.55) In the middle step, we used an integration by parts, noting th at the boundary terms at ±∞vanish, provided that u(x) is continuously differentiable and bounded as |x| → ∞. Pay attention to the minus sign in the final answer. 12/11/12 589 c/circlecopyrt2012 Peter J. Olver In the dual interpretation, the generalized function δ′ y(x) corresponds to the linear functional L′ y[u] =−u′(y) =/an}bracketle{tδ′ y;u/an}bracketri}ht=/integraldisplayℓ 0δ′ y(x)u(x)dx,where 0 < y < ℓ, (11.56) that maps a continuously differentiable function u(x) tominusits derivative at the point y. We note that (11.56) is compatible with a formal integratio n by parts /integraldisplayℓ 0δ′(x−y)u(x)dx=δ(x−y)u(x)/vextendsingle/vextendsingle/vextendsingle/vextendsingleℓ x=0−/integraldisplayℓ 0δ(x−y)u′(x)dx=−u′(y). The boundary terms at x= 0 and x=ℓautomatically vanish since δ(x−y) = 0 for x/ne}ationslash=y. Warning : The functions /tildewidegn(x) =gn(x)+g′ n(x) satisfy lim n→∞/tildewidegn(x) = 0 for all x/ne}ationslash=y, while/integraldisplay∞ −∞/tildewidegn(x)dx= 1. However, lim n→∞/tildewidegn= lim n→∞gn+ lim n→∞g′ n=δ+δ′. Thus, our original conditions (11.29–30) are notin fact sufficient to characterize whether a sequence of functions has the delta function as a limit. To be absolute ly sure, one must, in fact, verify the more comprehensive limiting formula (11.41). The Green’s Function To further cement our new-found friendship, we now put the de lta function to work to solve inhomogeneous boundary value problems. Consider a bar of length ℓsubject to a unit impulse force δy(x) =δ(x−y) concentrated at position 0 < y < ℓ. The underlying differential equation (11.12) takes the special form −d dx/parenleftbigg c(x)du dx/parenrightbigg =δ(x−y), 0< x < ℓ, (11.57) which we supplement with homogeneous boundary conditions t hat lead to a unique solu- tion. The solution is known as the Green’s function for the boundary value problem, and will be denoted by Gy(x) =G(x,y). Example 11.7. Let us look at the simple case of a homogeneous bar, of unit len gth ℓ= 1, with constant stiffness c, and fixed at both ends. The boundary value problem for the Green’s function G(x,y) takes the form −cu′′=δ(x−y), u (0) = 0 = u(1), (11.58) where 0< y <1 indicates the point at which we apply the impulse force. The solution to the differential equation is obtained by direct integration . First, by (11.43), u′(x) =−σ(x−y) c+a, whereais a constant of integration. A second integration leads to u(x) =−ρ(x−y) c+ax+b, (11.59) 12/11/12 590 c/circlecopyrt2012 Peter J. Olver 0.2 0.4 0.6 0.8 1 -0.2-0.10.10.20.3 y Figure 11.11. Green’s function for a Bar with Fixed Ends. whereρis the ramp function (11.45). The integration constants a,bare fixed by the boundary conditions; since 0 < y <1, we have u(0) =b= 0, u (1) =−1−y c+a+b= 0,and so a=1−y c. Therefore, the Green’s function for the problem is G(x,y) =−ρ(x−y)+(1−y)x=/braceleftbiggx(1−y)/c, x ≤y, y(1−x)/c, x ≥y,(11.60) Figure 11.11 sketches a graph of G(x,y) whenc= 1. Note that, for each fixed y, it is a continuous and piecewise affine function of x— meaning that its graph consists of connected straight line segments, with a corner where the un it impulse force is being applied. Once we have determined the Green’s function, we are able to s olve the general inho- mogeneous boundary value problem −cu′′=f(x), u (0) = 0 = u(1), (11.61) The solution formula is a consequence of linear superpositi on. We first express the forcing function f(x) as a linear combination of impulses concentrated at points along the bar. Since there is a continuum of possible positions 0 ≤y≤1 at which impulse forces may be applied, we will use an integral to sum them up, thereby writi ng the external force as f(x) =/integraldisplay1 0f(y)δ(x−y)dy. (11.62) In the continuous context, sums are replaced by integrals, a nd we will interpret (11.62) as the (continuous) superposition of an infinite collection of impulses f(y)δ(x−y), of magnitude f(y) and concentrated at position y. The superposition principle states that, for linear system s, linear combinations of inhomogeneities produce linear combinations of solutions . Again, we adapt this principle to the continuum by replacing the sums by integrals. Thus, we claim that the solution to the boundary value problem is the self-same linear superpos ition u(x) =/integraldisplay1 0f(y)G(x,y)dy (11.63) 12/11/12 591 c/circlecopyrt2012 Peter J. Olver of the Green’s function solutions to the individual unit imp ulse problems. For the particular boundary value problem (11.61), we use th e explicit formula (11.60) for the Green’s function. Breaking the integral (11.63) int o two parts, for y < xandy > x, we arrive at the explicit solution formula u(x) =1 c/integraldisplayx 0(1−x)yf(y)dy+1 c/integraldisplay1 xx(1−y)f(y)dy. (11.64) For example, under a constant unit force f, (11.64) reduces to u(x) =f c/integraldisplayx 0(1−x)ydy+f c/integraldisplay1 xx(1−y)dy=f 2c(1−x)x2+f 2cx(1−x)2=f 2c(x−x2), in agreement with our earlier solution (11.19) in the specia l cased= 0. Although this rel- atively simple problem was perhaps easier to solve directly , the Green’s function approach helps crystallize our understanding, and provides a unified framework that covers the full range of linear boundary value problems arising in applicat ions, including those governed by partial differential equations, [ 116,169,187]. Let us, finally, convince ourselves that the superposition f ormula (11.64) does indeed give the correct answer. First, cdu dx= (1−x)xf(x)+/integraldisplayx 0[−yf(y)]dy−x(1−x)f(x)+/integraldisplay1 x(1−y)f(y)dy =−/integraldisplay1 0yf(y)dy+/integraldisplay1 xf(y)dy. Differentiating again, we conclude that −cd2u dx2=f(x), as claimed. Remark: In computing the derivatives of u, we made use of the calculus formula d dx/integraldisplayβ(x) α(x)F(x,y)dy=F(x,β(x))dβ dx−F(x,α(x))dα dx+/integraldisplayβ(x) α(x)∂F ∂x(x,y)dy(11.65) for thederivativeof an integral withvariablelimits, whic h isa straightforward consequence of the Fundamental Theorem of Calculus and the chain rule, [ 9,168]. As with all limiting processes, one must always be careful when interchanging th e order of differentiation and integration. We note the following fundamental properties, that serve to uniquely characterize the Green’s function. First, since the delta forcing vanishes e xcept at the point x=y, the Green’s function satisfies the homogeneous differential equ ation† ∂2G ∂x2(x,y) = 0 for all x/ne}ationslash=y. (11.66) †SinceG(x,y)is afunctionof two variables, weswitchtopartial derivativenotati onto indicate its derivatives. 12/11/12 592 c/circlecopyrt2012 Peter J. Olver Secondly, by construction, it must satisfy the boundary con ditions, G(0,y) = 0 =G(1,y). Thirdly, for each fixed y,G(x,y)is a continuous function of x, but its derivative ∂G/∂xhas a jump discontinuity of magnitude −1/cat the impulse point x=y. The second derivative ∂2G/∂x2has a delta function discontinuity there, and thereby solve s the original impulse boundary value problem (11.58). Finally, we cannot help but notice that the Green’s function is a symmetric function of its two arguments: G(x,y) =G(y,x). Symmetry has the interesting physical consequence that the displacement of the bar at position xdue to an impulse force concentrated at position yis exactly the same as the displacement of the bar at ydue to an impulse of the same magnitude being applied at x. This turns out to be a rather general, although perhaps unanticipated phenomenon. (For the finite-dimensi onal counterpart for mass- spring chains, circuits, and structures see Exercises 6.1. 6, 6.2.18 and 6.3.17.) Symmetry is a consequence of the underlying symmetry or “self-adjoin tness” of the boundary value problem, to be developed properly in the following section. Remark: The Green’s function G(x,y) should be be viewed as the continuum limit of the inverse of the stiffness matrix, G=K−1, appearing in the discrete equilibrium equations Ku=f. Indeed, the entries Gijof the inverse matrix are approximations to the sampled values G(xi,xj). In particular, symmetry of the Green’s function, whereby G(xi,xj) =G(xj,xi), corresponds to symmetry, Gij=Gji, of the inverse of the symmetric stiffnessmatrix. InExercise ,youareaskedtostudythislimitingprocedureinsomedetai l. Let us summarize the fundamental properties that serve to ch aracterize the Green’s function, in a form that applies to general second order boun dary value problems. Basic Properties of the Green’s Function (i) Solves the homogeneous differential equation: −∂ ∂x/parenleftbigg c(x)∂ ∂xG(x,y)/parenrightbigg = 0,for all x/ne}ationslash=y. (11.67) (ii) Satisfies the homogeneous boundary conditions. (iii) Is a continuous function of its arguments. (iv) As a function of x, its derivative∂G ∂xhas a jump discontinuity of magnitude −1 c(y)atx=y. (v) Is a symmetric function of its arguments: G(x,y) =G(y,x). (11.68) (vi) Generates a superposition principle for the solution unde r general forcing functions: u(x) =/integraldisplayℓ 0G(x,y)f(y)dy. (11.69) 12/11/12 593 c/circlecopyrt2012 Peter J. Olver 0.2 0.4 0.6 0.8 1 -0.20.20.40.60.8 y Figure 11.12. Green’s Function for Bar with One Fixed and One Free End. . Example 11.8. Consider a uniform bar of length ℓ= 1 with one fixed and one free end, subject to an external force. The displacement u(x) satisfies the boundary value problem −cu′′=f(x), u (0) = 0, u′(1) = 0, (11.70) wherecis the elastic constant of the bar. To determine the Green’s f unction, we appeal to its characterizing properties, although one could equal ly well use direct integration as in the preceding Example 11.7. First, since, as a function of x, it must satisfy the homogeneous differential equation −cu′′= 0 for all x/ne}ationslash=y, the Green’s function must be of the form G(x,y) =/braceleftbiggpx+q, x≤y, rx+s, x≥y, for certain constants p,q,r,s. Second, the boundary conditions require q=G(0,y) = 0, r =∂G ∂x(1,y) = 0. Continuity of the Green’s function at x=yimposes the further constraint py=G(y−,y) =G(y+,y) =s. Finally, the derivative ∂G/∂xmust have a jump discontinuity of magnitude −1/catx=y, and so −1 c=∂G ∂x(y+,y)−∂G ∂x(y−,y) = 0−p, and so p=s=1 c. We conclude that the Green’s function for this problem is G(x,y) =/braceleftbiggx/c, x ≤y, y/c, x ≥y,(11.71) 12/11/12 594 c/circlecopyrt2012 Peter J. Olver which, for c= 1, is graphed in Figure 11.12. Note that G(x,y) =G(y,x) is indeed sym- metric, which helps check the correctness of our computatio n. Finally, the superposition principle (11.69) implies that the solution to the boundary value problem (11.70) can be written as a single integral, namely u(x) =/integraldisplay1 0G(x,y)f(y)dy=1 c/integraldisplayx 0yf(y)dy+1 c/integraldisplay1 xxf(y)dy. (11.72) The reader may wish to verify this directly, as we did in the pr evious example. 11.3. Adjoints and Minimum Principles. One of the profound messages of this text is that the linear al gebraic structures that were initially†designed for finite-dimensional problems all have direct co unterparts in the infinite-dimensional function spaces. To further develop t his theme, let us now discuss how the boundary value problems for continuous elastic bars fit i nto our general equilibrium framework of positive (semi-)definite linear systems. As we will see, the associated energy minimization principle not only leads to a new mathematical characterization of the equi- librium solution, it also, through the finite element method , underlies the most important class of numerical approximation algorithms for such bound ary value problems. Adjoints of Differential Operators In discrete systems, a key step was the recognition that the m atrix appearing in the force balance law is the transpose of the incidence matrix re lating displacements and elon- gations. In the continuum limit, the discrete incidence mat rix has turned into a differential operator. But how do you take the “transpose” of a differentia l operator? The abstract answer to this quandary can be found in Section 7.5. The trans pose of a matrix is a par- ticular instance of the general notion of the adjointof a linear function, which relies on the specification of inner products on its domain and target spac es. In the case of the matrix transpose, the adjoint is taken with respect to the standard dot product on Euclidean space. Thus, the correct interpretation of the “transpose” of a differential operator is as the adjoint linear operator with respect to suitable inner p roducts on function space. For bars and similar one-dimensional media, the role of the i ncidence matrix is played by the derivative v=D[u] =du/dx, which defines a linear operator D:U→Vfrom the vector space of possible displacements u(x), denoted by U, to the vector space of possible strainsv(x), denoted by V. In order to compute its adjoint, we need to impose inner products on both the displacement space Uand the strain space V. The simplest is to adopt the same standard L2inner product /an}bracketle{tu;/tildewideu/an}bracketri}ht=/integraldisplayℓ 0u(x)/tildewideu(x)dx, /an}bracketle{t/an}bracketle{tv;/tildewidev/an}bracketri}ht/an}bracketri}ht=/integraldisplayℓ 0v(x)/tildewidev(x)dx, (11.73) †Sometimes, the order is reversed, and, at least historically, basic l inear algebra concepts make their first appearance in function space. Examples include the Cauch y–Schwarz inequality, the Fredholm alternative, and the Fourier transform. 12/11/12 595 c/circlecopyrt2012 Peter J. Olver on both vector spaces. These are the continuum analogs of the Euclidean dot product, and, as we shall see, will be appropriate when dealing with ho mogeneous bars. According to the defining equation (7.74), the adjoint D∗of the derivative operator must satisfy the inner product identity /an}bracketle{t/an}bracketle{tD[u];v/an}bracketri}ht/an}bracketri}ht=/an}bracketle{tu;D∗[v]/an}bracketri}htfor all u∈U, v∈V. (11.74) First, we compute the left hand side: /an}bracketle{t/an}bracketle{tD[u];v/an}bracketri}ht/an}bracketri}ht=/angbracketleftbigg /angbracketleftbiggdu dx;v/angbracketrightbigg /angbracketrightbigg =/integraldisplayℓ 0du dxv dx. (11.75) On the other hand, the right hand side should equal /an}bracketle{tu;D∗[v]/an}bracketri}ht=/integraldisplayℓ 0uD∗[v]dx. (11.76) Now, in the latter integral, we see umultiplying the result of applying the linear operator D∗tov. To identify this integrand with that in the previous integr al (11.75), we need to somehow remove the derivative from u. The secret is integration by parts, which allows us to rewrite the first integral in the form /integraldisplayℓ 0du dxvdx=/bracketleftbig u(ℓ)v(ℓ)−u(0)v(0)/bracketrightbig −/integraldisplayℓ 0udv dxdx. (11.77) Ignoring the two boundary terms for a moment, the remaining i ntegral has the form of an inner product −/integraldisplayℓ 0udv dxdx=/integraldisplayℓ 0u/bracketleftbigg −dv dx/bracketrightbigg dx=/angbracketleftbigg u;−dv dx/angbracketrightbigg =/an}bracketle{tu;−D[v]/an}bracketri}ht. (11.78) Equating (11.75) and (11.78), we deduce that /an}bracketle{t/an}bracketle{tD[u];v/an}bracketri}ht/an}bracketri}ht=/angbracketleftbigg /angbracketleftbiggdu dx;v/angbracketrightbigg /angbracketrightbigg =/angbracketleftbigg u;−dv dx/angbracketrightbigg =/an}bracketle{tu;−D[v]/an}bracketri}ht. Thus, to satisfy the adjoint equation (11.74), we must have /an}bracketle{tu;D∗[v]/an}bracketri}ht=/an}bracketle{tu;−D[v]/an}bracketri}htfor all u∈U, v∈V, and so/parenleftbiggd dx/parenrightbigg∗ =D∗=−D=−d dx. (11.79) The final equation confirms our earlier identification of the d erivative operator Das the continuum limit of the incidence matrix A, and its negative −D=D∗as the limit of the transposed (or adjoint) incidence matrix AT=A∗. However, thepreceding argument isvalid onlyiftheboundary termsintheintegration by parts formula (11.77) vanish: u(ℓ)v(ℓ)−u(0)v(0) = 0, (11.80) 12/11/12 596 c/circlecopyrt2012 Peter J. Olver which necessitates imposing suitable boundary conditions on the functions uandv. For example, in the case of a bar with both ends fixed, the boundary conditions u(0) = 0, u (ℓ) = 0, (11.81) will ensure that (11.80) holds, and therefore validate (11. 79). The homogeneous boundary conditions serve to define the vector space U=/braceleftbig u(x)∈C2[0,ℓ]/vextendsingle/vextendsingleu(0) =u(ℓ) = 0/bracerightbig of allowable displacements, consisting of all twice contin uously differentiable functions that vanish at the ends of the bar. The fixed boundary conditions (11.81) are not the only possib ilities that ensure the vanishing of the boundary terms (11.80). An evident alterna tive is to require that the strain vanish at both endpoints, v(0) =v(ℓ) = 0. In this case, the strain space V=/braceleftbig v(x)∈C1[0,ℓ]/vextendsingle/vextendsinglev(0) =v(ℓ) = 0/bracerightbig consists of all functions that vanish at the endpoints. Sinc e the derivative D:U→V must map a displacement u(x) to anallowable strainv(x), the vector space of possible displacements takes the form U=/braceleftbig u(x)∈C2[0,ℓ]/vextendsingle/vextendsingleu′(0) =u′(ℓ) = 0/bracerightbig . Thus, this case corresponds to the free boundary conditions of Example 11.3. Again, restricting D:U→Vto these particular vector spaces ensures that the boundary terms (11.80) vanish, and so (11.79) holds in this situation too. Let us list the most important combinations of boundary cond itions that imply the vanishing of the boundary terms (11.80), and so ensure the va lidity of the adjoint equation (11.79). Self-Adjoint Boundary Conditions for a Bar (a) Both ends fixed: u(0) =u(ℓ) = 0. (b) One free and one fixed end: u(0) =u′(ℓ) = 0 or u′(0) =u(ℓ) = 0. (c) Both ends free: u′(0) =u′(ℓ) = 0. (d) Periodic bar or ring: u(0) =u(ℓ), u′(0) =u′(ℓ). In all cases, the boundary conditions impose restrictions o n the displacement space Uand, in cases ( b–d) when identifying v(x) =u′(x), the strain space Valso. In mathematics, a fixed boundary condition, u(a) = 0, is commonly referred to as aDirichlet boundary condition , to honor the nineteenth century French analyst Lejeune Dirichlet. A free boundary condition, u′(a) = 0, is known as a Neumann boundary con- dition, after his German contemporary Carl Gottfried Neumann. The Dirichlet boundary value problem (a) has both ends fixed, while the Neumann boundary value problem (c) has both ends free. The intermediate case ( b) is known as a mixed boundary value problem . The periodic boundary conditions ( d) represent a bar that has its ends joined together to 12/11/12 597 c/circlecopyrt2012 Peter J. Olver formacircular†elasticring, andrepresents thecontinuum limitoftheperi odicmass–spring chain discussed in Exercise 6.3.11. Summarizing, for a homogeneous bar with unit stiffness c(x)≡1, the displacement, strain, and external force are related by the adjoint formul ae v=D[u] =u′, f =D∗[v] =−v′, provided that we impose a suitable pair of homogeneous boundary condition s. The equi- librium equation has the self-adjoint form K[u] =f, where K=D∗◦D=−D2. (11.82) We note that K∗= (D∗◦D)∗=D∗◦(D∗)∗=D∗◦D=K, (11.83) which proves self-adjointness of the differential operator . In gory detail, /an}bracketle{tK[u];/tildewideu/an}bracketri}ht=/integraldisplayℓ 0/bracketleftbig −u′′(x)/tildewideu(x)/bracketrightbig dx=/integraldisplayℓ 0/bracketleftbig −u(x)/tildewideu′′(x)/bracketrightbig dx=/an}bracketle{tu;K[/tildewideu]/an}bracketri}ht(11.84) for all displacements u,/tildewideu∈U. A direct verification of this formula relies on two inte- gration by parts, employing the selected boundary conditio ns to eliminate the boundary contributions. To deal with nonuniform materials, we must modify the inner p roducts. Let us retain the ordinary L2inner product /an}bracketle{tu;/tildewideu/an}bracketri}ht=/integraldisplayℓ 0u(x)/tildewideu(x)dx, u, /tildewideu∈U, (11.85) on the vector space of possible displacements, but adopt a we ighted inner product /an}bracketle{t/an}bracketle{tv;/tildewidev/an}bracketri}ht/an}bracketri}ht=/integraldisplayℓ 0v(x)/tildewidev(x)c(x)dx, v, /tildewidev∈V, (11.86) on the space of strain functions. The weight function c(x)>0 coincides with the stiffness of the bar; its positivity, which is required for (11.86) to d efine abona fide inner product, is in accordance with the underlying physical assumptions. Let us recompute the adjoint of the derivative operator D:U→V, this time with respect to the inner products (11.85–86). Now we need to comp are /an}bracketle{t/an}bracketle{tD[u];v/an}bracketri}ht/an}bracketri}ht=/integraldisplayℓ 0du dxv(x)c(x)dx,with /an}bracketle{tu;D∗[v]/an}bracketri}ht=/integraldisplayℓ 0u(x)D∗[v]dx. Integrating the first expression by parts, we find /integraldisplayℓ 0du dxcvdx=/bracketleftbig u(ℓ)c(ℓ)v(ℓ)−u(0)c(0)v(0)/bracketrightbig −/integraldisplayℓ 0ud(cv) dxdx=/integraldisplayℓ 0u/bracketleftbigg −d(cv) dx/bracketrightbigg dx, (11.87) †The circle is sufficiently large so that we can safely ignore any curvatur e effects. 12/11/12 598 c/circlecopyrt2012 Peter J. Olver provided that we choose our boundary conditions so that u(ℓ)c(ℓ)v(ℓ)−u(0)c(0)v(0) = 0. (11.88) As you can check, this follows from any of the listed boundary conditions: Dirichlet, Neumann, or mixed, as well as the periodic case, assuming c(0) =c(ℓ). Therefore, in such situations, the weighted adjoint of the derivative operato r is D∗[v] =−d(cv) dx=−cdv dx−c′v. (11.89) The self-adjoint combination K=D∗◦Dis K[u] =−d dx/parenleftbigg c(x)du dx/parenrightbigg , (11.90) and hence we have formulated the original equation (11.12) f or a nonuniform bar in the same abstract self-adjoint form. As an application, let us show how the self-adjoint formulat ion leads directly to the symmetry of the Green’s function G(x,y). As a function of x, the Green’s function satisfies K[G(x,y)] =δ(x−y). Thus, by the definition of the delta function and the self-adj ointness identity (11.84), G(z,y) =/integraldisplayℓ 0G(x,y)δ(x−z)dx=/an}bracketle{tG(x,y);δ(x−z)/an}bracketri}ht=/an}bracketle{tG(x,y);K[G(x,z)]/an}bracketri}ht(11.91) =/an}bracketle{tK[G(x,y)];G(x,z)/an}bracketri}ht=/an}bracketle{tδ(x−y);G(x,z)/an}bracketri}ht=/integraldisplayℓ 0G(x,z)δ(x−y)dx=G(y,z), for any 0 < y,z < ℓ , which validates†the symmetry equation (11.68). Positivity and Minimum Principles We are now able to characterize the solution to a stable self- adjoint boundary value problem by a quadratic minimization principle. Again, the d evelopment shadows the finite-dimensional case presented in Chapter 6. So the first s tep is to understand how a differential operator defining a boundary value problem can b e positive definite. According to the abstract Definition 7.59, a linear operator K:U→Uon an inner product space Uispositive definite , provided that it is (a) self-adjoint, so K∗=K, and (b) satisfies the positivity criterion /an}bracketle{tK[u];u/an}bracketri}ht>0 for all 0 /ne}ationslash=u∈U. Self-adjointness of the product operator K=D∗◦Dwas established in (11.83). Further- more, Theorem 7.62 tells us that Kis positive definite if and only if ker D={0}. Indeed, by the definition of the adjoint, /an}bracketle{tK[u];u/an}bracketri}ht=/an}bracketle{tD∗[D[u]];u/an}bracketri}ht=/an}bracketle{t/an}bracketle{tD[u];D[u]/an}bracketri}ht/an}bracketri}ht=/bardblD[u]/bardbl2≥0, (11.92) †Symmetry at the endpoints follows from continuity. 12/11/12 599 c/circlecopyrt2012 Peter J. Olver soK=D∗◦Dis automatically positive semi-definite. Moreover, /an}bracketle{tK[u];u/an}bracketri}ht= 0 if and only ifD[u] = 0, i.e., u∈kerK. Thus ker D={0}is both necessary and sufficient for the positivity criterion to hold. Now, in the absence of constraints, the kernel of the derivat ive operator Disnot trivial. Indeed, D[u] =u′= 0 if and only if u(x)≡cis constant, and hence ker Dis the one-dimensional subspace of C1[0,ℓ] consisting of all constant functions. However, we are viewing Das a linear operator on the vector space Uof allowable displacements, and so the elements of ker D⊂Umust also be allowable, meaning that they must satisfy the boundary conditions. Thus, positivity reduces, in the pres ent situation, to the question of whether or not there are any nontrivial constant function s that satisfy the prescribed homogeneous boundary conditions. Clearly, the only constant function that satisfies a homogen eous Dirichlet boundary condition is the zero function. Therefore, when restricted to the Dirichlet displacement spaceU={u(0) =u(ℓ) = 0}, the derivative operator has trivial kernel, ker D={0}, soK=D∗◦Ddefines a positive definite linear operator on U. A similar argument applies to the mixed boundary value problem, which is also po sitive definite. On the other hand, any constant function satisfies the homogeneous Neuma nn boundary conditions, and so ker D⊂/tildewideU={u′(0) =u′(ℓ) = 0}is a one-dimensional subspace. Therefore, the Neumann boundary value problem is only positive semi-de finite. A similar argument shows that the periodic problem is also positive semi-defini te. Observe that, just as in the finite-dimensional version, the positive definite cases are stable, and the boundary value problem admitsaunique equilibriumsolutionunder arbitra ryexternalforcing, whereas the semi-definite cases are unstable, and have either no solutio n or infinitely many equilibrium solutions, depending on the nature of the external forcing. In the positive definite, stable cases, we can characterize t he equilibrium solution as the unique function u∈Uthat minimizes the quadratic functional P[u] =1 2/bardblD[u]/bardbl2−/an}bracketle{tu;f/an}bracketri}ht=/integraldisplayℓ 0/bracketleftbig1 2c(x)u′(x)2−f(x)u(x)/bracketrightbig dx. (11.93) A proof of this general fact appears following Theorem 7.61. Pay attention: the norm in (11.93) refers to the strain space V, and so is associated with the weighted inner product (11.86), whereas the inner product term refers to the displa cement space U, which has been given the L2inner product. Physically, the first term measures the inter nal energy due to the stress in the bar, while the second term is the poten tial energy induced by the external forcing. Thus, as always, the equilibrium solutio n seeks to minimize the total energy in the system. Example 11.9. Consider the homogeneous Dirichlet boundary value problem −u′′=f(x), u (0) = 0, u (ℓ) = 0. (11.94) for a uniform bar with two fixed ends. This is a stable case, and so the underlying differential operator K=D∗◦D=−D2, when acting on the space of displacements satisfying the boundary conditions, is positive definite. E xplicitly, positive definiteness 12/11/12 600 c/circlecopyrt2012 Peter J. Olver requires /an}bracketle{tK[u];u/an}bracketri}ht=/integraldisplayℓ 0[−u′′(x)u(x)]dx=/integraldisplayℓ 0[u′(x)]2dx >0 (11 .95) for all nonzero u(x)/ne}ationslash≡0 withu(0) =u(ℓ) = 0. Notice how we used an integration by parts, invoking the boundary conditions to eliminate the boundary contributions, to expose the positivity of the integral. The corresponding energy funct ional is P[u] =1 2/bardblu′/bardbl2−/an}bracketle{tu;f/an}bracketri}ht=/integraldisplayℓ 0/bracketleftbig1 2u′(x)2−f(x)u(x)/bracketrightbig dx. Itsminimumvalue, takenoverallpossibledisplacementfun ctionsthatsatisfytheboundary conditions, occurs precisely when u=u⋆is the solution to the boundary value problem. A direct verification of the latter fact may be instructive. A s in our derivation of the adjoint operator, it relies on an integration by parts. Sinc e−u′′ ⋆=f, P[u] =/integraldisplayℓ 0/bracketleftbig1 2(u′)2+u′′ ⋆u/bracketrightbig dx=u′ ⋆(ℓ)u(ℓ)−u′ ⋆(0)u(0)+/integraldisplayℓ 0/bracketleftbig1 2(u′)2−u′ ⋆u′/bracketrightbig dx =/integraldisplayℓ 01 2(u′−u′ ⋆)2dx−/integraldisplayℓ 01 2(u′ ⋆)2dx, (11.96) where the boundary terms vanish owing to the boundary condit ions onu⋆andu. In the final expression for P[u], the first integral is always ≥0, and is actually equal to 0 if and only ifu′(x) =u′ ⋆(x) for all 0 ≤x≤ℓ. On the other hand, the second integral does not depend upon uat all. Thus, for P[u] to achieve a minimum, u(x) =u⋆(x)+cfor some constant c. But the boundary conditions force c= 0, and hence the energy functional will assume its minimum value if and only if u=u⋆. Inhomogeneous Boundary Conditions So far, we have restricted our attention to homogeneous boun dary value problems. Inhomogeneous boundary conditions are a little trickier, s ince the spaces of allowable displacements and allowablestrains are no longer vector sp aces, and so the abstract theory, as developed in Chapter 7, is not directly applicable. One wayto circumvent thisdifficulty isto slightlymodify the displacement function so asto satisfy homogeneous boundary conditions. Consider, f or example, the inhomogeneous Dirichlet boundary value problem K[u] =−d dx/parenleftbigg c(x)du dx/parenrightbigg =f(x), u (0) =α, u (ℓ) =β. (11.97) We shall choose a function h(x) that satisfies the boundary conditions: h(0) =α, h (ℓ) =β. Note that we are notrequiring hto satisfy the differential equation, and so one, but by no means the only, possible choice is the linear interpolating polynomial h(x) =α+β−α ℓx. (11.98) 12/11/12 601 c/circlecopyrt2012 Peter J. Olver Sinceuandhhave the same boundary values, their difference /tildewideu(x) =u(x)−h(x) (11 .99) satisfies the homogeneous Dirichlet boundary conditions /tildewideu(0) =/tildewideu(ℓ) = 0. (11.100) Moreover, by linearity, /tildewideusatisfies the modified equation K[/tildewideu] =K[u−h] =K[u]−K[h] =f−K[h]≡/tildewidef, or, explicitly, −d dx/parenleftbigg c(x)d/tildewideu dx/parenrightbigg =/tildewidef(x),where /tildewidef(x) =f(x)+d dx/parenleftbigg c(x)dh dx/parenrightbigg .(11.101) For the particular choice (11.98), /tildewidef(x) =f(x)+β−α ℓc′(x). Thus, we have managed to convert the original inhomogeneous problem for uinto a ho- mogeneous boundary value problem for /tildewideu. Once we have solved the latter, the solution to the original is simply reconstructed from the formula u(x) =/tildewideu(x)+h(x). (11.102) We know that the homogeneous Dirichlet boundary value probl em (11.100–101) is positive definite, and so we can characterize its solution by a minimum principle, namely as the minimizer of the quadratic energy functional P[/tildewideu] =1 2/bardbl/tildewideu′/bardbl2−/an}bracketle{t/tildewideu;/tildewidef/an}bracketri}ht=/integraldisplayℓ 0/bracketleftbig1 2c(x)/tildewideu′(x)2−/tildewidef(x)/tildewideu(x)/bracketrightbig dx. (11.103) Let us rewrite the minimization principle in terms of the ori ginal displacement function u(x). Replacing /tildewideuand/tildewidefby their formulae (11.99,101) yields P[/tildewideu] =1 2/bardblu′−h′/bardbl2−/an}bracketle{tu−h;f−K[h]/an}bracketri}ht =/bracketleftbig1 2/bardblu′/bardbl2−/an}bracketle{tu;f/an}bracketri}ht/bracketrightbig −/bracketleftbig /an}bracketle{t/an}bracketle{tu′;h′/an}bracketri}ht/an}bracketri}ht−/an}bracketle{tu;K[h]/an}bracketri}ht/bracketrightbig +/bracketleftbig1 2/bardblh′/bardbl2+/an}bracketle{th;f−K[h]/an}bracketri}ht/bracketrightbig =P[u]−/bracketleftbig /an}bracketle{t/an}bracketle{tu′;h′/an}bracketri}ht/an}bracketri}ht−/an}bracketle{tu;K[h]/an}bracketri}ht/bracketrightbig +C0. (11.104) In the middle expression, the last pair of terms depend only o n the initial choice of h(x), and not on u(x); thus, once hhas been selected, they can be regarded as a fixed constant, here denoted by C0. The first pair of terms reproduces the quadratic energy func tional (11.93) for the actual displacement u(x). The middle terms can be explicitly evaluated: /an}bracketle{t/an}bracketle{tu′;h′/an}bracketri}ht/an}bracketri}ht−/an}bracketle{tu;K[h]/an}bracketri}ht=/integraldisplayℓ 0/bracketleftbig c(x)h′(x)u′(x)+/parenleftbig c(x)h′(x)/parenrightbig′u(x)/bracketrightbig dx =/integraldisplayℓ 0d dx/bracketleftbig c(x)h′(x)u(x)/bracketrightbig dx=c(ℓ)h′(ℓ)u(ℓ)−c(0)h′(0)u(0).(11.105) 12/11/12 602 c/circlecopyrt2012 Peter J. Olver Figure 11.13. Bending of a Beam. In particular, if u(x) satisfies the inhomogeneous Dirichlet boundary condition su(0) =α, u(ℓ) =β, then /an}bracketle{t/an}bracketle{tu′;h′/an}bracketri}ht/an}bracketri}ht−/an}bracketle{tu;K[h]/an}bracketri}ht=c(ℓ)h′(ℓ)β−c(0)h′(0)α=C1 also depends only on the interpolating function hand not on u. Therefore, P[/tildewideu] =P[u]−C1+C0 differ by a constant. We conclude that, if the function /tildewideuminimizes P[/tildewideu], thenu=/tildewideu+h necessarily minimizes P[u]. In this manner, we have characterized the solution to the inhomogeneous Dirichlet boundary value problem by the sameminimization principle. Theorem 11.10. The solution u⋆(x)to the Dirichlet boundary value problem −d dx/parenleftbigg c(x)du dx/parenrightbigg =f(x), u (0) =α, u (ℓ) =β, is the unique C2function that satisfies the indicated boundary conditions a nd minimizes the energy functional P[u] =/integraldisplayℓ 0/bracketleftbig1 2c(x)u′(x)2−f(x)u(x)/bracketrightbig dx. Warning : The inhomogeneous mixed boundary value problem is trickie r, since the extra terms (11.105) willdepend upon the value of u(x). The details are worked out in Exercise . 11.4. Beams and Splines. Unlike a bar, which can only stretch longitudinally, a beamis allowed to bend. To keep the geometry simple, we treat the case in which the beam i s restricted to the xy plane, as sketched in Figure 11.13. Let 0 ≤x≤ℓrepresent the reference position along a horizontal beam of length ℓ. To further simplify the physics, we shall ignore stretchin g, and assume that the “atoms” in the beam can only move in the tra nsverse direction, with y=u(x) representing the vertical displacement of the “atom” that starts out at position x. 12/11/12 603 c/circlecopyrt2012 Peter J. Olver Thestrainin a beam depends on how much it is bent. Mathematically, bend ing is equal to the curvature†of the graph of the displacement function u(x), and is computed by the usual calculus formula κ=u′′ /parenleftbig 1+(u′)2/parenrightbig3/2. (11.106) Thus, for beams, the strain is a nonlinear function of displacement. Since we are still only willing to deal with linear systems, we shall suppress the no nlinearity by assuming that the beam is not bent too far; more specifically, we assume that the derivative u′(x)≪1 is small and so the tangent line is nearly horizontal. Under thi s assumption, the curvature function (11.106) is replaced by its linear approximation κ≈u′′. (11.107) From now on, we will identify v=D2[u] =u′′as thestrainin a bending beam. The second derivative operator L=D2that maps displacement uto strain v=L[u] thereby describes the beam’s intrinsic (linearized) geometry. The next step is to formulate a constitutive relation betwee n stress and strain. Phys- ically, the stressw(x) represents the bending moment of the beam, defined as the pro duct of internal force and angular deflection. Our small bending a ssumption implies an elastic Hooke’s law relation w(x) =c(x)v(x) =c(x)d2u dx2, (11.108) where the proportionality factor c(x)>0 measures the stiffness of the beam at the point x. In particular, a uniform beam has constant stiffness, c(x)≡c. Finally, the differential equation governing the equilibri um configuration of the beam will follow from a balance of the internal and external force s. To compute the internal force, we appeal to our general equilibrium framework, whic h tells us to apply the adjoint of the incidence operator L=D2to the strain, leading to the force balance law L∗[v] =L∗◦L[u] =f. (11.109) Let us compute the adjoint. We use the ordinary L2inner product on the space of displacements u(x), and adopt a weighted inner product, based on the stiffness f unction c(x), between strain functions: /an}bracketle{tu;/tildewideu/an}bracketri}ht=/integraldisplayb au(x)/tildewideu(x)dx, /an}bracketle{t/an}bracketle{tv;/tildewidev/an}bracketri}ht/an}bracketri}ht=/integraldisplayb av(x)/tildewidev(x)c(x)dx. (11.110) According to the general adjoint equation (7.74), we need to equate /integraldisplayℓ 0L[u]vcdx=/an}bracketle{t/an}bracketle{tL[u];v/an}bracketri}ht/an}bracketri}ht=/an}bracketle{tu;L∗[v]/an}bracketri}ht=/integraldisplayℓ 0uL∗[v]dx. (11.111) †By definition, [ 9,168], the curvature of a curve at a point is equal to the reciprocal, κ= 1/r of the radius of the osculating circle; see Exercise A.5.11 for details . 12/11/12 604 c/circlecopyrt2012 Peter J. Olver Simply Supported End Clamped End Free End Sliding End Figure 11.14. Boundary Conditions for a Beam. As before, the computation relies on (in this case two) integ rations by parts: /an}bracketle{t/an}bracketle{tL[u];v/an}bracketri}ht/an}bracketri}ht=/integraldisplayℓ 0d2u dx2cvdx=/bracketleftbiggdu dxcv/bracketrightbigg/vextendsingle/vextendsingle/vextendsingle/vextendsingleℓ x=0−/integraldisplayℓ 0du dxd(cv) dxdx =/bracketleftbiggdu dxcv−ud(cv) dx/bracketrightbigg/vextendsingle/vextendsingle/vextendsingle/vextendsingleℓ x=0+/integraldisplayℓ 0ud2(cv) dx2dx. Comparing with (11.111), we conclude that L∗[v] =D2(cv) provided the boundary terms vanish: /bracketleftbiggdu dxcv−ud(cv) dx/bracketrightbigg/vextendsingle/vextendsingle/vextendsingle/vextendsingleℓ x=0=/bracketleftbiggdu dxw−udw dx/bracketrightbigg/vextendsingle/vextendsingle/vextendsingle/vextendsingleℓ x=0(11.112) =/bracketleftbig u′(ℓ)w(ℓ)−u(ℓ)w′(ℓ)/bracketrightbig −/bracketleftbig u′(0)w(0)−u(0)w′(0)/bracketrightbig = 0. Thus, under suitable boundary conditions, the force balanc e equations are L∗[v] =d2(cv) dx2=f(x). (11.113) A justification of (11.113) based on physical principles can be found in [ 181]. Combining (11.108,113), we conclude that the equilibrium configurati on of the beam is characterized as a solution to the fourth order ordinary differential equat ion d2 dx2/parenleftbigg c(x)d2u dx2/parenrightbigg =f(x). (11.114) As such, the general solution will depend upon 4 arbitrary co nstants, and so we need to impose a total of four boundary conditions — two at each end — i n order to uniquely specify the equilibrium displacement. The (homogeneous) b oundary conditions should be chosen so as to make the boundary terms in our integration by p arts computation vanish, cf. (11.112). There are a variety of ways in which this can be a rranged, and the most important possibilities are the following: 12/11/12 605 c/circlecopyrt2012 Peter J. Olver Self-Adjoint Boundary Conditions for a Beam (a) Simply supported end: u(0) =w(0) = 0, (b) Fixed (clamped) end: u(0) =u′(0) = 0, (c) Free end: w(0) =w′(0) = 0, (d) Sliding end: u′(0) =w′(0) = 0. In these conditions, w(x) =c(x)v(x) =c(x)u′′(x) is the stress resulting from the displace- mentu(x). A second pair of boundary conditions must be imposed at the ot her endx=ℓ. You can mix or match these conditions in any combination — for exa mple, a pair of simply supported ends, or one free end and one fixed end, and so on. Inh omogeneous boundary conditions are also allowed and used to model applied displa cements or applied forces at the ends. Yet another option is to consider a bendable circul ar ring, which is subject to periodic boundary conditions u(0) =u(ℓ), u′(0) =u′(ℓ), w(0) =w(ℓ), w′(0) =w′(ℓ), indicating that the ends of the beam have been welded togethe r. Let us concentrate our efforts on the uniform beam, of unit len gthℓ= 1, choosing units so that its stiffness c(x)≡1. In the absence of external forcing, the differential equation (11.114) reduces to the elementary fourth order or dinary differential equation dxu d4x= 0. (11.115) The general solution is an arbitrary cubic polynomial, u=ax3+bx2+cx+d. (11.116) Let us use this formula to solve a couple of representative bo undary value problems. First, suppose we clamp both ends of the beam, imposing the bo undary conditions u(0) = 0, u′(0) =β, u (1) = 0, u′(1) = 0, (11.117) so that the left end is tilted by a (small) angle tan−1β. We substitute the solution for- mula (11.116) into the boundary conditions (11.117) and sol ve for a=β, b =−2β, c =β, d = 0. The resulting solution u(x) =β(x3−2x2+x) =βx(1−x)2(11.118) is known as a Hermite cubic spline†and is graphed in Figure 11.15. †We first met Charles Hermite in Section 3.6, and the term “spline” wil l be explained shortly. 12/11/12 606 c/circlecopyrt2012 Peter J. Olver 0.2 0.4 0.6 0.8 1 -0.10.10.20.30.40.5 Figure 11.15. Hermite Cubic Spline. As a second example, suppose that we raise the left hand end of the beam without tilting, which corresponds to the boundary conditions u(0) =α, u′(0) = 0, u (1) = 0, u′(1) = 0. (11.119) Substituting (11.116) and solving for a,b,c,d, we find that the solution is u(x) =α(1−x)2(2x+1). (11.120) Observe that if we simultaneously raise and tilt the left end , sou(0) =α,u′(0) =β, then we can simply use superposition to write the solution as the s um of (11.118) and (11.120): u(x) =α(1−x)2(2x+1)+βx(1−x)2. To analyze a forced beam, we can adapt the Green’s function ap proach. As we know, the Green’s function will depend on the choice of (homogeneo us) boundary conditions. Let us treat the case when the beam has two fixed ends, and so u(0) = 0, u′(0) = 0, u (1) = 0, u′(1) = 0. (11.121) To construct the Green’s function, we must solve the forced d ifferential equation dxu d4x=δ(x−y) (11 .122) corresponding to a concentrated unit impulse applied at pos itionyalong the beam. Inte- grating (11.122) four times, using (11.46) with n= 4, we produce the general solution u(x) =ax3+bx2+cx+d+/braceleftbigg1 6(x−y)3, x > y, 0, x < y, to the differential equation (11.122). The boundary conditi ons (11.121) require u(0) =d= 0, u (1) =a+b+1 6(1−y)3= 0, u′(0) =c= 0, u′(1) = 3a+2b+1 2(1−y)2= 0, 12/11/12 607 c/circlecopyrt2012 Peter J. Olver 0.2 0.4 0.6 0.8 1 -0.02-0.010.010.020.03 y Figure 11.16. Green’s Function for a Beam with Two Fixed Ends. and hence a=1 3(1−y)3−1 2(1−y)2, b =−1 2(1−y)3+1 2(1−y)2. Therefore, the Green’s function is G(x,y) =/braceleftigg1 6x2(1−y)2(3y−x−2xy), x < y, 1 6y2(1−x)2(3x−y−2xy), x > y.(11.123) Observe that, as with the second order bar system, the Green’ s function is symmetric, G(x,y) =G(y,x), which is a manifestation of the self-adjointness of the un derlying bound- aryvalueproblem, cf.(11.91). Symmetryimpliesthatthede flectionofthebeamatposition xdue to a concentrated impulse force applied at position yis the same as the deflection atydue to an impulse force of the same magnitude applied at x. As a function of x, the Green’s function G(x,y) satisfies the homogeneous differential equation (11.115) for all x/ne}ationslash=y. Its first and second derivatives ∂G/∂x,∂2G/∂x2are continuous, while ∂3G/∂x3has a unit jump discontinuity at x=y, which then produces the required delta function impulse in ∂4G/∂x4. The Green’s function (11.123) is graphed in Figure 11.16, and appears to be quite smooth. Evidently, t he human eye cannot easily discern discontinuities in third order derivatives! The solution to the forced boundary value problem dxu d4x=f(x), u (0) =u′(0) =u(1) =u′(1) = 0, (11.124) for a beam with fixed ends is then obtained by invoking the supe rposition principle. We view the forcing function as a linear superposition f(x) =/integraldisplayℓ 0f(y)δ(x−y)dx 12/11/12 608 c/circlecopyrt2012 Peter J. Olver 0.2 0.4 0.6 0.8 1 -0.003-0.002-0.0010.0010.0020.003 Figure 11.17. Deflection of a Uniform Beam under Gravity. of impulse delta forces. The solution is the self-same linea r superposition of Green’s func- tion responses: u(x) =/integraldisplay1 0G(x,y)f(y)dy (11.125) =1 6/integraldisplayx 0y2(1−x)2(3x−y−2xy)f(y)dy+1 6/integraldisplay1 xx2(1−y)2(3y−x−2xy)f(y)dy. For example, under a constant unit downwards force f(x)≡1, e.g., gravity, the deflection of the beam is given by u(x) =1 24x4−1 12x3+1 24x2=1 24x2(1−x)2, and graphed in Figure 11.17. Although we could, of course, ob tainu(x) by integrating the original differential equation (11.124) directly, writing the solution formula (11.125) as a single integral has evident advantages. Since the beam operator K=L∗◦Lassumes the standard self-adjoint, positive semi- definite form, the boundary value problem will be positive de finite and hence stable if and only if ker L= kerD2={0}when restricted to the space of allowable displacement functions. Since the second derivative D2annihilates all linear polynomials u(x) =α+βx, positive definiteness requires that no non-zero linear poly nomials satisfy all four homo- geneous boundary conditions. For example, any beam with one fixed end is stable since u(x)≡0 is the only linear polynomial that satisfies u(0) =u′(0) = 0. On the other hand, a beam with two free ends is unstable since every linear polyn omial displacement has zero stressw(x) =u′′(x)≡0, and so satisfies the boundary conditions w(0) =w′(0) =w(ℓ) = w′(ℓ) = 0. Similarly, a beam with a simply supported plus a free end isnot positivedefinite sinceu(x) =βxsatisfies the four boundary conditions u(0) =u′(0) = 0,w(ℓ) =w′(ℓ) = 0. In the stable cases, the equilibrium solution can be charact erized as the unique minimizer 12/11/12 609 c/circlecopyrt2012 Peter J. Olver of the quadratic energy functional† P[u] =1 2/bardblL[u]/bardbl2−/an}bracketle{tu;f/an}bracketri}ht=/integraldisplayb a/bracketleftbig1 2c(x)u′′(x)2−f(x)u(x)/bracketrightbig dx (11.126) among all C4functions satisfying the homogeneous boundary conditions . Inhomogeneous boundary conditions require some extra analysis, since the required integration by parts may introduce additional boundary contributions. Splines In pre–CAD (computer aided design) draftsmanship, a splinewas a long, thin, flexible strip of wood that was used to draw a smooth curve through pres cribed points. The points were marked by small pegs, and the spline rested on the pegs. T he mathematical theory of splines was first developed in the 1940’s by the Romanian ma thematician Isaac Schoen- berg as an attractive alternative to polynomial interpolat ion and approximation. Splines have since become ubiquitous in numerical analysis, in geom etric modeling, in design and manufacturing, in computer graphics and animation, and in m any other applications. We suppose that the spline coincides with the graph of a funct iony=u(x). The pegs are fixed at the prescribed data points ( x0,y0),...,(xn,yn), and this requires u(x) to satisfy the interpolation conditions u(xj) =yj, j = 0,...,n. (11.127) Themesh points x0< x1< x2<···< xnare distinct and labeled in increasing order. The spline is modeled as an elastic beam, and so satisfies the h omogeneous beam equation (11.115). Therefore, u(x) =aj+bj(x−xj)+cj(x−xj)2+dj(x−xj)3,xj≤x≤xj+1, j= 0,...,n−1,(11.128) is a piecewise cubic function — meaning that, between succes sive mesh points, it is a cubic polynomial, but not necessarily the same cubic on each subin terval. The fact that we write the formula (11.128) in terms of x−xjis merely for computational convenience. Our problem is to determine the coefficients aj, bj, cj, dj, j = 0,...,n−1. Since there are nsubintervals, there are a total of 4 ncoefficients, and so we require 4 n equations to uniquely prescribe them. First, we need the spl ine to satisfy the interpolation conditions (11.127). Since it is defined by a different formul a on each side of the mesh point, this results in a total of 2 nconditions: u(x+ j) =aj=yj, u(x− j+1) =aj+bjhj+cjh2 j+djh3 j=yj+1,j= 0,...,n−1,(11.129) †Keep in mind that the norm on the strain functions v=L[u] =u′′is based on the weighted inner product /angbracketleft/angbracketleftv;/tildewidev/angbracketright/angbracketrightin (11.110). 12/11/12 610 c/circlecopyrt2012 Peter J. Olver where we abbreviate the length of the jthsubinterval by hj=xj+1−xj. The next step is to require that the spline be as smooth as poss ible. The interpola- tion conditions (11.129) guarantee that u(x) is continuous. The condition u(x)∈C1be continuously differentiable requires that u′(x) be continuous at the interior mesh points x1,...,xn−1, which imposes the n−1 additional conditions bj+2cjhj+3djh2 j=u′(x− j+1) =u′(x+ j+1) =bj+1, j = 0,...,n−2.(11.130) To make u∈C2, we impose n−1 further conditions 2cj+6djhj=u′′(x− j+1) =u′′(x+ j+1) = 2cj+1, j = 0,...,n−2,(11.131) to ensure that u′′is continuous at the mesh points. We have now imposed a total o f 4n−2 conditions, namely (11.129–131), on the 4 ncoefficients. The two missing constraints will come from boundary conditions at the two endpoints, namely x0andxn. There are three common types: (i)Natural boundary conditions :u′′(x0) =u′′(xn) = 0, whereby c0= 0, cn−1+3dn−1hn−1= 0. (11.132) Physically, this models a simply supported spline that rest s freely on the first and last pegs. (ii)Clamped boundary conditions :u′(x0) =α, u′(xn) =β,whereα,β, which could be 0, are fixed by the user. This requires b0=α, bn−1+2cn−1hn−1+3dn−1h2 n−1=β. (11.133) This corresponds to clamping the spline at prescribed angle s at each end. (iii)Periodic boundary conditions :u′(x0) =u′(xn), u′′(x0) =u′′(xn),so that b0=bn−1+2cn−1hn−1+3dn−1h2 n−1, c0=cn−1+3dn−1hn−1.(11.134) If we also require that the end interpolation values agree, u(x0) =y0=yn=u(xn), (11.135) then the resulting spline will be a periodic C2function, so u(x+p) =u(x) withp=xn−x0 for allx. The periodic case is used to draw smooth closed curves; see b elow. Theorem 11.11. Suppose we are given mesh points a=x0< x1<···< xn=b, and corresponding data values y0,y1,...,yn, along with one of the three kinds of boundary conditions (11.132),(11.133), or(11.134). Then there exists a unique piecewise cubic spline function u(x)∈C2[a,b]that interpolates the data, u(x0) =y0,...,u(xn) =yn, and satisfies the boundary conditions. 12/11/12 611 c/circlecopyrt2012 Peter J. Olver Proof: We first discuss the natural case. The clamped case is left as an exercise for the reader, while the slightly harder periodic case will be t reated at the end of the section. The first set of equations in (11.129) says that aj=yj, j = 0,...,n−1. (11.136) Next, (11.131–132) imply that dj=cj+1−cj 3hj. (11.137) This equation also holds for j=n−1, provided that we make the convention that† cn= 0. We now substitute (11.136–137) into the second set of equati ons in (11.129), and then solve the resulting equation for bj=yj+1−yj hj−(2cj+cj+1)hj 3. (11.138) Substituting this result and (11.137) back into (11.130), a nd simplifying, we find hjcj+2(hj+hj+1)cj+1+hj+1cj+2= 3/bracketleftigg yj+2−yj+1 hj+1−yj+1−yj hj/bracketrightigg =zj+1,(11.139) where we introduce zj+1as a shorthand for the quantity on the right hand side. In the case of natural boundary conditions, we have c0= 0, cn= 0, and so (11.139) constitutes a tridiagonal linear system Ac=z, (11.140) for the unknown coefficients c=/parenleftbig c1,c2,...,cn−1/parenrightbigT, with coefficient matrix A= 2(h0+h1)h1 h12(h1+h2)h2 h22(h2+h3) h3 ......... hn−32(hn−3+hn−2) hn−2 hn−2 2(hn−2+hn−1)  (11.141) and right hand side z=/parenleftbig z1,z2,...,zn−1/parenrightbigT. Once (11.141) has been solved, we will then use (11.136–138) to reconstruct the other spline coefficient saj,bj,dj. †This is merely for convenience; there is no cnused in the formula for the spline. 12/11/12 612 c/circlecopyrt2012 Peter J. Olver 1 2 3 4 -1-0.50.511.52 Figure 11.18. A Cubic Spline. The key observation is that the coefficient matrix Aisstrictly diagonally dominant , cf. Definition 10.36, because all the hj>0, and so 2(hj−1+hj)> hj−1+hj. Theorem 10.37 implies that Ais nonsingular, and hence the tridiagonal linear system has a unique solution c. This suffices to prove the theorem in the case of natural bound ary conditions. Q.E.D. To actually solve the linear system (11.140), we can apply ou r tridiagonal solution algorithm (1.66). Let us specialize to the most important ca se, when the mesh points are equally spaced in the interval [ a,b], so that xj=a+jh,where h=hj=b−a n, j = 0,...,n−1. In this case, the coefficient matrix A=hBis equal to htimes the tridiagonal matrix B= 4 1 1 4 1 1 4 1 1 4 1 1 4 1 .........  that first appeared in Example 1.37. Its LUfactorization takes on an especially simple form, since most of the entries of LandUare essentially the same decimal numbers. This makes theimplementation ofthe Forward andBack Substituti on procedures almost trivial. Figure 11.18 shows a particular example — a natural spline pa ssing through the data points (0 ,0), (1,2), (2,−1), (3,1), (4,0). As with the Green’s function for the beam, the human eye is unable to discern the discontinuities in its thi rd derivatives, and so the graph appears completely smooth, even though it is, in fact, only C2. 12/11/12 613 c/circlecopyrt2012 Peter J. Olver In the periodic case, we set an+k=an, bn+k=bn, cn+k=cn, dn+k=dn, zn+k=zn. With this convention, the basic equations (11.136–139) are the same. In this case, the coefficient matrix for the linear system Ac=z,with c=/parenleftbig c0,c1,...,cn−1/parenrightbigT,z=/parenleftbig z0,z1,...,zn−1/parenrightbigT, is ofcirculant tridiagonal form: A= 2(hn−1+h0)h0 hn−1 h02(h0+h1)h1 h12(h1+h2)h2 ......... hn−32(hn−3+hn−2)hn−2 hn−1 hn−22(hn−2+hn−1) . (11.142) AgainAisstrictlydiagonallydominant, and so thereisa unique sol utionc, from which one reconstructs the spline, proving Theorem 11.11 in the perio dic case. The LUfactorization of tridiagonal circulant matrices was discussed in Exercis e 1.7.14. One immediate application of splines is curve fitting in comp uter aided design and graphics. The basic problem is to draw a smooth parametrized curveu(t) = (u(t),v(t))T that passes through a set of prescribed data points xk= (xk,yk)Tin the plane. We have the freedom to choose the parameter value t=tkwhen the curve passes through the kth point; the simplest and most common choice is to set tk=k. We then construct the functions x=u(t) andy=v(t) as cubic splines interpolating the xandycoordinates of the data points, so u(tk) =xk,v(tk) =yk. For smooth closed curves, we require that both splines be periodic; for curves with ends, either natural or clamped boundary conditions are used. Most computer graphics packages include one or more impleme ntations of parametri- zed spline curves. The same idea also underlies modern font d esign for laser printing and typography (including the fonts used in this book). The grea t advantage of spline fonts over their bitmapped counterparts is that they can be readil y scaled. Some sample let- ter shapes parametrized by periodic splines passing throug h the indicated data points are plotted in Figure 11.19. Better fits can be easily obtained by increasing the number of data points. Various extensions of the basic spline algorithms t o space curves and surfaces are an essential component of modern computer graphics, design , and animation, [ 64,163]. 11.5. Sturm–Liouville Boundary Value Problems. Thesystemsthatgoverntheequilibriumconfigurationsofba rsareparticularinstances of a very general class of second order boundary value proble ms that was first systemat- ically investigated by the nineteenth century French mathe maticians Jacques Sturm and Joseph Liouville. Sturm–Liouville boundary value problem s appear in a very wide range 12/11/12 614 c/circlecopyrt2012 Peter J. Olver Figure 11.19. Three Sample Spline Letters. of applications, particularly in the analysis of partial di fferential equations by the method of separation of variables. A partial list of applications i ncludes (a) heat conduction in non-uniform bars; (b) vibrations of non-uniform bars and strings; (c) quantum mechanics — the one-dimensional Schr¨ odinger equ ation; (d) scattering theory — Hill’s equation; (e) oscillations of circular membranes (vibrations of drums) — Bessel’s equation; (f) oscillations of a sphere — Legendre’s equation; (g) thermodynamics of cylindrical and spherical bodies. Inthissection, wewillshow howtheclassofSturm–Liouvill eboundary valueproblems fits into our general equilibrium framework. However, the mo st interesting cases will be deferred until needed in our analysis of partial differentia l equations in Chapters 17 and 18. The general Sturm–Liouville boundary value problem is based on a second order ordi- nary differential equation of the form −d dx/parenleftbigg p(x)du dx/parenrightbigg +q(x)u=−p(x)d2u dx2−p′(x)du dx+q(x)u=f(x),(11.143) which is supplemented by Dirichlet, Neumann, mixed, or peri odicboundary conditions. To be specific, let us concentrate on the case of homogeneous Dir ichlet boundary conditions u(a) = 0, u (b) = 0. (11.144) To avoid singular points of the differential equation (altho ugh we will later discover that most cases of interest in physics have one or more singul ar points) , we assume thatp(x)>0 for all a≤x≤b. To ensure positive definiteness of the Sturm–Liouville differential operator, we also assume q(x)≥0. These assumptions suffice to guarantee existence and uniqueness of the solution to the boundary val ue problem. A proof of the following theorem can be found in [ 117]. Theorem 11.12. Letp(x)>0andq(x)≥0fora≤x≤b. Then the Sturm– Liouville boundary value problem (11.143–144)admits a unique solution. Most Sturm–Liouville problems cannot be solved in terms of e lementary functions. Indeed, most of the important special functions appearing i n mathematical physics, in- cluding Bessel functions, Legendre functions, hypergeome tric functions, and so on, first arise as solutions to particular Sturm–Liouville equation s, [144]. 12/11/12 615 c/circlecopyrt2012 Peter J. Olver 0.2 0.4 0.6 0.8 1 -0.050.050.10.15 y Figure 11.20. Green’s Function for the Constant Coefficient Sturm–Liouville Problem. Example 11.13. Consider the constant coefficient Sturm–Liouville boundary value problem −u′′+ω2u=f(x), u (0) =u(1) = 0. (11.145) The functions p(x)≡1 andq(x)≡ω2>0 are both constant. We will solve this problem by constructing the Green’s function. Thus, we first conside r the effect of a delta function inhomogeneity −u′′+ω2u=δ(x−y), u (0) =u(1) = 0. (11.146) Rather than try to integrate this differential equation dire ctly, let us appeal to the defining properties of the Green’s function. The general solution to the homogeneous equation is a linear combination of the two basic exponentials eωxande−ωx, or better, the hyperbolic functions coshωx=eωx+e−ωx 2, sinhωx=eωx−e−ωx 2. (11.147) The solutions satisfying the first boundary condition are mu ltiples of sinh ωx, while those satisfying the second boundary condition are multiples of s inhω(1−x). Therefore, the solution to (11.146) has the form G(x,y) =/braceleftbiggasinhωx, x < y, bsinhω(1−x), x > y.(11.148) Continuity of G(x,y) atx=yrequires asinhωy=bsinhω(1−y). (11.149) Atx=y, the derivative ∂G/∂xmust have a jump discontinuity of magnitude −1 in order that the second derivative term in (11.146) match the delta f unction. Since ∂G ∂x(x,y) =/braceleftbiggaωcoshωx, x < y, −bωcoshω(1−x), x > y, the jump condition requires aωcoshωy−1 =−bωcoshω(1−y). (11.150) 12/11/12 616 c/circlecopyrt2012 Peter J. Olver If we multiply (11.149) by ωcoshω(1−y) and (11.150) by sinh ω(1−y) and then add the results together, we find sinhω(1−y) =aω/bracketleftbig sinhωycoshω(1−y)+coshωysinhω(1−y)/bracketrightbig =aωsinhω, where we used the addition formula for the hyperbolic sine: sinh(α+β) = sinhαcoshβ+coshαsinhβ. (11.151) Therefore, a=sinhω(1−y) ωsinhω, b =sinhωy ωsinhω, and the Green’s function is G(x,y) =  sinhωxsinhω(1−y) ωsinhω, x < y, sinhω(1−x) sinhωy ωsinhω, x > y.(11.152) Note that G(x,y) =G(y,x) is symmetric, in accordance with the self-adjoint nature o f the boundary value problem. A graph appears in Figure 11.20; not e that the corner, indicating a discontinuity in the first derivative, appears at the point x=ywhere the impulse force is applied. The general solution to the inhomogeneous boundary value pr oblem (11.145) is given by the basic superposition formula (11.63), which becomes u(x) =/integraldisplay1 0G(x,y)f(y)dy =/integraldisplayx 0sinhω(1−x)sinhωy ωsinhωf(y)dy+/integraldisplay1 xsinhωxsinhω(1−y) ωsinhωf(y)dy. For example, under a constant unit force f(x)≡1, the solution is u(x) =/integraldisplayx 0sinhω(1−x)sinhωy ωsinhωdy+/integraldisplay1 xsinhωxsinhω(1−y) ωsinhωdy =sinhω(1−x)/parenleftbig coshωx−1/parenrightbig ω2sinhω+sinhωx/parenleftbig coshω(1−x)−1/parenrightbig ω2sinhω =1 ω2−sinhωx+sinhω(1−x) ω2sinhω.(11.153) For comparative purposes, the reader may wish to rederive th is particular solution by a direct calculation, without appealing to the Green’s funct ion. To place a Sturm–Liouville boundary value problem in our sel f-adjoint framework, we proceed as follows. (Exercise serves to motivate the construction.) Consider the linear operator L[u] =/parenleftbigg u′ u/parenrightbigg 12/11/12 617 c/circlecopyrt2012 Peter J. Olver that maps u(x) to the vector-valued function whose components are the fun ction and its first derivative. For the homogeneous Dirichlet boundary co nditions (11.144), the domain ofLwill be the vector space U=/braceleftbig u(x)∈C2[a,b]/vextendsingle/vextendsingleu(a) =u(b) = 0/bracerightbig consisting of all twice continuously differentiable functi ons that vanish at the endpoints. Thetargetspaceof L:U→Vconsistsofcontinuouslydifferentiablevector-valuedfun ctions v(x) = (v1(x),v2(x))T; we denote this vector space as V= C1([a,b],R2). Toproceed, wemustcomputetheadjointof L:U→V. TorecovertheSturm–Liouville problem, we use the standard L2inner product (11.85) on U, but adopt a weighted inner product /an}bracketle{t/an}bracketle{tv;w/an}bracketri}ht/an}bracketri}ht=/integraldisplayb a/bracketleftbig p(x)v1(x)w1(x)+q(x)v2(x)w2(x)/bracketrightbig dx,v=/parenleftbigg v1 v2/parenrightbigg ,w=/parenleftbigg w1 w2/parenrightbigg , (11.154) onV. The positivity assumptions on the weight functions p,qensure that this is a bona fideinner product. According to the defining equation (7.74), th e adjoint L∗:V→Uis required to satisfy /an}bracketle{t/an}bracketle{tL[u];v/an}bracketri}ht/an}bracketri}ht=/an}bracketle{tu;L∗[v]/an}bracketri}ht. As usual, the adjoint computation relies on integration by p arts. Here, we only need to manipulate the first summand: /an}bracketle{t/an}bracketle{tL[u];v/an}bracketri}ht/an}bracketri}ht=/integraldisplayb a/bracketleftbig pu′v1+quv2/bracketrightbig dx =p(b)u(b)v1(b)−p(a)u(a)v1(a)+/integraldisplayb au[−(pv1)′+qv2]dx. The Dirichlet conditions (11.144) ensure that the boundary terms vanish, and therefore, /an}bracketle{t/an}bracketle{tL[u];v/an}bracketri}ht/an}bracketri}ht=/integraldisplayb au[−(pv1)′+qv2]dx=/an}bracketle{tu;L∗[v]/an}bracketri}ht. We conclude that the adjoint operator is given by L∗[v] =−d(pv1) dx+qv2. The canonical self-adjoint combination K[u] =L∗◦L[u] =L∗/parenleftbigg u′ u/parenrightbigg =−d dx/parenleftbigg pdu dx/parenrightbigg +qu (11.155) reproduces the Sturm–Liouville differential operator. Mor eover, since ker L={0}is trivial (why?), the boundary value problem is positive definite. The orem 7.62 implies that the solution can be characterized as the unique minimizer of the quadratic functional P[u] =1 2/bardblL[u]/bardbl2−/an}bracketle{tu;f/an}bracketri}ht=/integraldisplayb a/bracketleftbig1 2p(x)u′(x)2+1 2q(x)u(x)2−f(x)u(x)/bracketrightbig dx(11.156) 12/11/12 618 c/circlecopyrt2012 Peter J. Olver among all C2functions satisfying the prescribed boundary conditions. For example, the solution to the constant coefficient Sturm–Liouville proble m (11.145) can be characterized as minimizing the quadratic functional P[u] =/integraldisplay1 0/bracketleftbig1 2u′2+1 2ω2u2−fu/bracketrightbig dx among all C2functions satisfying u(0) =u(1) = 0. 11.6. Finite Elements. Thecharacterizationofthesolutiontoapositivedefiniteb oundaryvalueproblemviaa minimizationprinciple inspires a very powerful and widely used numerical solution scheme, known as the finite element method . In this final section, we give a brief introduction to the finite element method in the context of one-dimensiona l boundary value problems involving ordinary differential equations. Extensions to b oundary value problems in higher dimensions governed by partial differential equations will appear in Section 15.5. The underlying idea isstrikinglysimple. Wearetryingto fin d thesolutiontoa bound- ary value problem by minimizing a quadratic functional P[u] on an infinite-dimensional vector space U. The solution u⋆∈Uto this minimization problem is found by solving a differential equation subject to specified boundary conditi ons. However, as we learned in Chapter 4, minimizing the functional on a finite-dimensional subspace W⊂Uis a prob- lem in linear algebra, and, moreover, one that we already kno w how to solve! Of course, restricting the functional P[u] to the subspace Wwill not, barring luck, lead to the exact minimizer. Nevertheless, if we choose Wto be a sufficiently “large” subspace, the result- ing minimizer w⋆∈Wmay very well provide a reasonable approximation to the actu al solution u⋆∈U. A rigorous justification of this process, under appropriat e hypotheses, requires a full analysis of the finite element method, and we r efer the interested reader to [174,197]. Here we shall concentrate on trying to understand how to ap ply the method in practice. To be a bit more explicit, consider the minimization princip le P[u] =1 2/bardblL[u]/bardbl2−/an}bracketle{tf;u/an}bracketri}ht (11.157) for the linear system K[u] =f,where K=L∗◦L, representing our boundary value problem. The norm in (11.15 7) is typically based on some form of weighted inner product /an}bracketle{t/an}bracketle{tv;/tildewidev/an}bracketri}ht/an}bracketri}hton the space of strains v=L[u]∈V, while the inner product term /an}bracketle{tf;u/an}bracketri}htis typically (although not necessarily) unweighted on the s pace of displacements u∈U. The linear operator takes the self-adjoint form K=L∗◦L, and mustbepositivedefinite—whichrequiresker L={0}. Withoutthepositivityassumption, the boundary value problem has either no solutions, or infini tely many; in either event, the basic finite element method will not apply. Rather than try to minimize P[u] on the entire function space U, we now seek to minimizeitonasuitablychosenfinite-dimensional subspac eW⊂U. Webeginbyselecting 12/11/12 619 c/circlecopyrt2012 Peter J. Olver abasis†ϕ1,...,ϕnofthesubspace W. Thegeneralelementof Wisa(uniquelydetermined) linear combination ϕ(x) =c1ϕ1(x)+···+cnϕn(x) (11 .158) of the basis functions. Our goal, then, is to determine the co efficients c1,...,cnsuch that ϕ(x) minimizes P[ϕ] among all such functions. Substituting (11.158) into (11. 157) and expanding we find P[ϕ] =1 2n/summationdisplay i,j=1mijcicj−n/summationdisplay i=1bici=1 2cTMc−cTb, (11.159) where (a)c= (c1,c2,...,cn)Tis the vector of unknown coefficients in (11.158), (b)M= (mij) is the symmetric n×nmatrix with entries mij=/an}bracketle{t/an}bracketle{tL[ϕi];L[ϕj]/an}bracketri}ht/an}bracketri}ht, i,j = 1,...,n, (11.160) (c)b= (b1,b2,...,bn)Tis the vector with entries bi=/an}bracketle{tf;ϕi/an}bracketri}ht, i = 1,...,n. (11.161) Observe that, once we specify the basis functions ϕi, the coefficients mijandbiare all known quantities. Therefore, we have reduced our original p roblem to a finite-dimensional problem of minimizing the quadratic function (11.159) over all possible vectors c∈Rn. The coefficient matrix Mis, in fact, positive definite, since, by the preceding compu tation, cTMc=n/summationdisplay i,j=1mijcicj=/bardblL[c1ϕ1(x)+···+cnϕn]/bardbl2=/bardblL[ϕ]/bardbl2>0 (11.162) as long as L[ϕ]/ne}ationslash= 0. Moreover, our positivityassumption implies that L[ϕ] = 0 if and only ifϕ≡0, and hence (11.162)is indeed positivefor all c/ne}ationslash=0. We can now invoke the original finite-dimensional minimization Theorem 4.1 to conclude th at the unique minimizer to (11.159) is obtained by solving the associated linear syste m Mc=b. (11.163) Solving (11.163)relies on some form of Gaussian Eliminatio n, or, alternatively, an iterative linear system solver, e.g., Gauss–Seidel or SOR. This constitutes the basic abstract setting for the finite el ement method. The main issue, then, ishow toeffectively choose the finite-dimensio nal subspace W. Two candidates that might spring to mind are the space P(n)of polynomials of degree ≤n, or the space T(n)of trigonometric polynomials of degree ≤n, the focus of Chapter 12. However, for a variety of reasons, neither is well suited to the finite eleme nt method. One criterion is that the functions in Wmust satisfy the relevant boundary conditions — otherwise Wwould †In this case, an orthonormal basis is not of any particular help. 12/11/12 620 c/circlecopyrt2012 Peter J. Olver 0.2 0.4 0.6 0.8 10.20.40.60.8 Figure 11.21. A Continuous Piecewise Affine Function. not be a subspace of U. More importantly, in order to obtain sufficient accuracy, th e linear algebraic system (11.163) will typically be rather large, a nd so the coefficient matrix M should be as sparse as possible, i.e., have lots of zero entri es. Otherwise, computing the solution will be too time-consuming to be of much practical v alue. Such considerations prove to be of absolutely crucial importance when applying t he method to solve boundary value problems for partial differential equations in higher dimensions. The really innovative contribution of the finite element met hod is to first (paradox- ically)enlargethe space Uof allowable functions upon which to minimize the quadratic functional P[u]. The governing differential equation requires its solutio ns to have a cer- tain degree of smoothness, whereas the associated minimiza tionprinciple typicallyrequires only half as many derivatives. Thus, for second order bounda ry value problems, including bars, (11.93), and general Sturm–Liouville problems, (11. 156),P[u] only involves first or- der derivatives. It can be rigorously shown that the functio nal has the sameminimizing solution, even if one allows (reasonable) functions that fa il to have enough derivatives to satisfy the differential equation. Thus, one can try minimiz ing over subspaces contain- ing fairly “rough” functions. Again, the justification of th is method requires some deeper analysis, which lies beyond the scope of this introductory t reatment. For second order boundary value problems, a popular and effec tive choice of the finite- dimensionalsubspaceistousecontinuous, piecewiseaffinef unctions. Recallthatafunction is affine, f(x) =ax+b, if and only if its graph is a straight line. The function is piecewise affineif its graph consists of a finite number of straight line segme nts; a typical example is plotted in Figure 11.21. Continuity requires that the indiv idual line segments be connected together end to end. Given a boundary value problem on a bounded interval [ a,b], let us fix a finite col- lection of mesh points a=x0< x1< x2<···< xn−1< xn=b. The formulas simplify if one uses equally spaced mesh points , but this is not necessary for the method to apply. Let Wdenote the vector space consisting of all continuous, piece - wise affine functions, with corners at the nodes, that satisfy the homogeneous boundary conditions. To be specific, let us treat the case of Dirichlet (fixed) boundary conditions ϕ(a) =ϕ(b) = 0. (11.164) 12/11/12 621 c/circlecopyrt2012 Peter J. Olver 1 2 3 4 5 6 7 -0.20.20.40.60.811.2 Figure 11.22. A Hat Function. Thus, on each subinterval ϕ(x) =cj+bj(x−xj),forxj≤x≤xj+1, j= 0,...,n−1. Continuity of ϕ(x) requires cj=ϕ(x+ j) =ϕ(x− j) =cj−1+bj−1hj−1, j = 1,...,n−1, (11.165) wherehj−1=xj−xj−1denotes the length of the jthsubinterval. The boundary conditions (11.164) require ϕ(a) =c0= 0, ϕ (b) =cn−1+hn−1bn−1= 0. (11.166) The function ϕ(x) involves a total of 2 nunspecified coefficients c0,...,cn−1,b0,...,bn−1. The continuity conditions (11.165) and the second boundary condition (11.166) uniquely determine the bj. The first boundary condition specifies c0, while the remaining n−1 coefficients c1=ϕ(x1),...,cn−1=ϕ(xn−1) are arbitrary. We conclude that the finite element subspace Whas dimension n−1, which is the number of interior mesh points. Remark: Every function ϕ(x) in our subspace has piecewise constant first derivative w′(x). However, the jump discontinuities in ϕ′(x) imply that its second derivative ϕ′′(x) hasadeltafunctionimpulseateachmeshpoint, andistheref orefarfrombeingasolutionto the differential equation. Nevertheless, the finite element minimizer ϕ⋆(x) will, in practice, provide a reasonable approximation to the actual solution u⋆(x). The most convenient basis for Wconsists of the hat functions , which are continuous, piecewise affine functions that interpolate the same basis da ta as the Lagrange polynomials (4.47), namely ϕj(xk) =/braceleftbigg1, j=k, 0, j/ne}ationslash=k,for j= 1,...,n−1, k= 0,...,n. The graph of a typical hat function appears in Figure 11.22. T he explicit formula is easily 12/11/12 622 c/circlecopyrt2012 Peter J. Olver established: ϕj(x) =  x−xj−1 xj−xj−1, xj−1≤x≤xj, xj+1−x xj+1−xj, xj≤x≤xj+1, 0, x ≤xj−1orx≥xj+1,j= 1,...,n−1.(11.167) An advantage of using these basis elements is that the result ing coefficient matrix (11.160) turns out to be tridiagonal. Therefore, the tridiagonal Gau ssian Elimination algorithm in (1.66) will rapidly produce the solutionto the linear syste m (11.163). Since the accuracy of the finite element solution increases with the number of mesh points, this solution scheme allows us to easily compute very accurate numerical approxi mations. Example 11.14. Consider the equilibrium equations K[u] =−d dx/parenleftbigg c(x)du dx/parenrightbigg =f(x), 0< x < ℓ, for a non-uniform bar subject to homogeneous Dirichlet boun dary conditions. In order to formulate a finite element approximationscheme, we begin wi ththe minimizationprinciple (11.93) based on the quadratic functional P[u] =1 2/bardblu′/bardbl2−/an}bracketle{tf;u/an}bracketri}ht=/integraldisplayℓ 0/bracketleftbig1 2c(x)u′(x)2−f(x)u(x)/bracketrightbig dx. Wedividetheinterval[0 ,ℓ]intonequalsubintervals, eachoflength h=ℓ/n. Theresulting uniform mesh has xj=jh=jℓ n, j = 0,...,n. The corresponding finite element basis hat functions are exp licitly given by ϕj(x) =  (x−xj−1)/h, xj−1≤x≤xj, (xj+1−x)/h, xj≤x≤xj+1, 0, otherwise ,j= 1,...,n−1.(11.168) The associated linear system (11.163) has coefficient matrix entries mij=/an}bracketle{t/an}bracketle{tϕ′ i;ϕ′ j/an}bracketri}ht/an}bracketri}ht=/integraldisplayℓ 0ϕ′ i(x)ϕ′ j(x)c(x)dx, i,j = 1,...,n−1. Since the function ϕi(x) vanishes except on the interval xi−1< x < xi+1, whileϕj(x) vanishes outside xj−1< x < xj+1, the integral will vanish unless i=jori=j±1. Moreover, ϕ′ j(x) =  1/h, xj−1≤x≤xj, −1/h, xj≤x≤xj+1, 0, otherwise ,j= 1,...,n−1. 12/11/12 623 c/circlecopyrt2012 Peter J. Olver Therefore, the coefficient matrix has the tridiagonal form M=1 h2 s0+s1−s1 −s1s1+s2−s2 −s2s2+s3−s3 ......... −sn−3sn−3+sn−2−sn−2 −sn−2sn−2+sn−1 ,(11.169) where sj=/integraldisplayxj+1 xjc(x)dx (11.170) is the total stiffness of the jthsubinterval. For example, in the homogeneous case c(x)≡1, the coefficient matrix (11.169) reduces to the very special fo rm M=1 h 2−1 −1 2 −1 −1 2 −1 ......... −1 2 −1 −1 2 . (11.171) The corresponding right hand side has entries bj=/an}bracketle{tf;ϕj/an}bracketri}ht=/integraldisplayℓ 0f(x)ϕj(x)dx =1 h/bracketleftigg/integraldisplayxj xj−1(x−xj−1)f(x)dx+/integraldisplayxj+1 xj(xj+1−x)f(x)dx/bracketrightigg .(11.172) In this manner, we have assembled the basic ingredients for d etermining the finite element approximation to the solution. In practice, we do not have to explicitly evaluate the integr als (11.170,172), but may replace them by a suitably close numerical approximation. W henh≪1 is small, then the integrals are taken over small intervals, and we can use the t rapezoid rule†, [35,168], to approximate them: sj≈h 2/bracketleftbig c(xj)+c(xj+1)/bracketrightbig , bj≈hf(xj). (11.173) †One might be tempted use more accurate numerical integration procedu res, but the im- provement in accuracy of the final answer is not very significant, partic ularly if the step size his small. 12/11/12 624 c/circlecopyrt2012 Peter J. Olver 0.2 0.4 0.6 0.8 10.020.040.060.08 0.2 0.4 0.6 0.8 10.020.040.060.08 0.2 0.4 0.6 0.8 10.020.040.060.08 0.2 0.4 0.6 0.8 10.020.040.060.08 Figure 11.23. Finite Element Solution to (11.175). Remark: Thejthentryoftheresultingfiniteelement system Mc=bis, upondividing byh, given by −cj+1−2cj+cj−1 h2=−u(xj+1)−2u(xj)+u(xj−1) h2=−f(xj).(11.174) The left hand side coincides with the standard finite differen ce approximation to minus the second derivative −u′′(xj) at the mesh point xj. (Details concerning finite differences can be found in Section 14.6.) As a result, for this particular di fferential equation, the finite element and finite difference numerical solution methods hap pen to coincide. Example 11.15. Consider the boundary value problem −d dx(x+1)du dx= 1, u (0) = 0, u (1) = 0. (11.175) The explicit solution is easily found by direct integration : u(x) =−x+log(x+1) log2. (11.176) It minimizes the associated quadratic functional P[u] =/integraldisplayℓ 0/bracketleftbig1 2(x+1)u′(x)2−u(x)/bracketrightbig dx (11.177) over all possible functions u∈C1that satisfy the given boundary conditions. The fi- nite element system (11.163) has coefficient matrix given by ( 11.169) and right hand side (11.172), where sj=/integraldisplayxj+1 xj(1+x)dx=h(1+xj)+1 2h2=h+h2/parenleftbigg j+1 2/parenrightbigg , bj=/integraldisplayxj+1 xj1dx=h. 12/11/12 625 c/circlecopyrt2012 Peter J. Olver The resulting solution is plotted in Figure 11.23. The first t hree graphs contain, respec- tively, 5, 10, 20 points in the mesh, so that h=.2,.1,.05, while the last plots the exact solution (11.176). Even when computed on rather coarse mesh es, the finite element ap- proximation is quite respectable. Example 11.16. Consider the Sturm–Liouville boundary value problem −u′′+(x+1)u=xex, u (0) = 0, u (1) = 0. (11.178) The solution minimizes the quadratic functional (11.156), which in this particular case is P[u] =/integraldisplay1 0/bracketleftbig1 2u′(x)2+1 2(x+1)u(x)2−exu(x)/bracketrightbig dx, (11.179) over all functions u(x) that satisfy the boundary conditions. We lay out a uniform m esh of step size h= 1/n. The corresponding finite element basis hat functions as in ( 11.168). The matrix entries are given by† mij=/integraldisplay1 0/bracketleftbig ϕ′ i(x)ϕ′ j(x)+(x+1)ϕi(x)ϕj(x)/bracketrightbig dx≈  2 h+2h 3(xi+1), i=j, −1 h+h 6(xi+1),|i−j|= 1, 0, otherwise , while bi=/an}bracketle{txex;ϕi/an}bracketri}ht=/integraldisplay1 0xexϕi(x)dx≈xiexih. The resulting solution is plotted in Figure 11.24. As in the p revious figure, the first three graphs contain, respectively, 5, 10, 20 points in the mesh, w hile the last plots the exact solution, which can be expressed in terms of Airy functions, cf. [144]. So far, we have only treated homogeneous boundary condition s. An inhomogeneous boundary value problem does not immediately fit into our fram ework since the set of functions satisfying the boundary conditions does notform a vector space. As discussed at the end of Section 11.3, one way to get around this problem i s to replace u(x) by /tildewideu(x) =u(x)−h(x), where h(x) is any convenient function that satisfies the boundary conditions. For example, for the inhomogeneous Dirichlet c onditions u(a) =α, u (b) =β, we can subtract off the affine function h(x) =(β−α)x+αb−βa b−a. †The integration is made easier by noting that the integrand is zero exce pt on a small subin- terval. Since the function x+1 (but not ϕiorϕj) does not vary significantly on this subinterval, it can be approximated by its value 1+ xiat a mesh point. A similar simplification is used in the ensuing integral for bi. 12/11/12 626 c/circlecopyrt2012 Peter J. Olver 0.2 0.4 0.6 0.8 10.020.040.060.080.1 0.2 0.4 0.6 0.8 10.020.040.060.080.1 0.2 0.4 0.6 0.8 10.020.040.060.080.1 0.2 0.4 0.6 0.8 10.020.040.060.080.1 Figure 11.24. Finite Element Solution to (11.178). Another option is to choose an appropriate combination of el ements at the endpoints: h(x) =αϕ0(x)+βϕn(x). Linearity implies that the difference /tildewideu(x) =u(x)−h(x) satisfies the amended differential equation K[/tildewideu] =/tildewidef,where /tildewidef=f−K[h], now supplemented by homogeneous boundary conditions. The m odified boundary value problem can then be solved by the standard finite element meth od. Further details are left as a project for the motivated student. Finally, one can employ other functions beyond the piecewis e affine hat functions (11.167) to span finite element subspace. Another popular ch oice, which is essential for higher order boundary value problems such as beams, is to use splines. Thus, once we have chosen our mesh points, we can let ϕj(x) be the basis B–splines discussed in Exercise . Sinceϕj(x) = 0 for x≤xj−2orx≥xj+2, the resulting coefficient matrix (11.160) ispentadiagonal , which means mij= 0 whenever |i−j|>2. Pentadiagonal matrices are not quite as pleasant as their tridiagonal cousins, but a re still rather sparse. Positive definiteness of Mimplies that an iterative solution technique, e.g., SOR, ca n effectively and rapidly solve the linear system, and thereby produce the finite element spline approx- imation to the boundary value problem. Weak Solutions There is an alternative way of introducing the finite element solution method, which also applies when there is no convenient minimization princ iple available, based on an important analytical extension of the usual notion of what c onstitutes a solution to a differential equation. One reformulates the differential eq uation as an integral equation. The resulting “weak solutions”, which include non-classic al solutions with singularities and discontinuities, are particularly appropriate in the s tudy of discontinuous and non- smooth physical phenomena, such as shock waves, cracks and d islocations in elastic media, 12/11/12 627 c/circlecopyrt2012 Peter J. Olver singularities in liquid crystals, and so on; see [ 188] and Section 22.1 for details. The weak solution approach has the advantage that it applies eve n to equations that do not possess an associated minimization principle. However, th e convergence of the induced finite element scheme is harder to justify, and, indeed, not a lways valid. The starting point is a trivial observation: the only elemen t of an inner product space which is orthogonal to every other element is zero. More prec isely: Lemma 11.17. IfVis an inner product space, then /an}bracketle{tw;v/an}bracketri}ht= 0for allv∈Vif and only ifw=0. Proof: Choose v=w. The orthogonality condition implies 0 = /an}bracketle{tw;w/an}bracketri}ht=/bardblw/bardbl2, and sow=0. Q.E.D. Note that the result is equally valid in both finite- and infini te-dimensional vector spaces. Suppose we are trying to solve a linear†system K[u] =f, (11.180) whereK:U→Vis a linear operator between inner product spaces. Using the lemma, this can be reformulated as requiring /an}bracketle{t/an}bracketle{tK[u];v/an}bracketri}ht/an}bracketri}ht=/an}bracketle{t/an}bracketle{tf;v/an}bracketri}ht/an}bracketri}htfor all v∈V. According to the definition (7.74), one can replace Kby its adjoint K∗:W→V, and require /an}bracketle{tu;K∗[v]/an}bracketri}ht=/an}bracketle{t/an}bracketle{tf;v/an}bracketri}ht/an}bracketri}htfor all v∈V. (11.181) The latter is called the weak formulation of our original equation. The general philosophy is that one can check whether uis a weak solution to the system by evaluating it on various test elements vusing the weak form (11.181) of the system. In the finite-dimensional situation, when Kis merely multiplication by some matrix, the weak formulation is an unnecessary complication, and no t of use. However, in the infinite-dimensional situation, when Kisadifferentialoperator, thentheoriginalboundary value problem K[u] =frequires that ube sufficiently differentiable, whereas the weak version /an}bracketle{tu;K∗[ϕ]/an}bracketri}ht=/an}bracketle{t/an}bracketle{tf;ϕ/an}bracketri}ht/an}bracketri}htfor all ϕ requires only that the test function ϕ(x) be smooth. As a result, weak solutions are not restricted to be smooth functions possessing the required n umber of derivatives. Example 11.18. Consider the homogeneous Dirichlet boundary value problem K[u] =−d dx/parenleftbigg c(x)du dx/parenrightbigg =f(x),0< x < ℓ, u (0) =u(ℓ) = 0, †The method also straightforwardly extends to nonlinear systems. 12/11/12 628 c/circlecopyrt2012 Peter J. Olver for a nonuniform bar. Its weak version is obtained by integra tion by parts. We initially restrict to test functions which vanish at the boundary ϕ(0) =ϕ(ℓ) = 0. This requirement will eliminate any boundary terms in the integration by part s computation /an}bracketle{tK[u];ϕ/an}bracketri}ht=/integraldisplayℓ 0/bracketleftbigg −d dx/parenleftbigg c(x)du dx/parenrightbigg ϕ(x)/bracketrightbigg dx=−/integraldisplayℓ 0c(x)du dxdϕ dxdx =/integraldisplayℓ 0f(x)ϕ(x)dx=/an}bracketle{tf;ϕ/an}bracketri}ht.(11.182) This“semi-weak”formulationisknowninmechanicsasthe principle of virtual work ,[169]. For example, the Green’s function of the boundary value prob lem does not qualify as a classical solution since it is not twice continuously differ entiable, but can be formulated as a weak solution satisfying the virtual work equation with ri ght hand side defined by the delta forcing function. A second integration by parts produces the weak form (11.181 ) of the differential equation: /an}bracketle{tu;K[ϕ]/an}bracketri}ht=−/integraldisplayℓ 0u(x)d dx/parenleftbigg c(x)dϕ dx/parenrightbigg dx=/integraldisplayℓ 0f(x)ϕ(x)dx=/an}bracketle{tf;ϕ/an}bracketri}ht.(11.183) Now, even discontinuous functions u(x) are allowed as weak solutions. The goal is to findu(x) such that this condition holds for all smooth test function sϕ(x). For example, any function u(x) which satisfies the differential equation except at points o f discontinuity qualifies as a weak solution. Inafinite element orGalerkin approximation to theweak solution, one restrictsatten- tion to a finite-dimensional subspace Wspanned by functions ϕ1,...,ϕn−1, and requires that the approximate solution ϕ(x) =c1ϕ1(x)+···+cn−1ϕn−1(x) (11 .184) satisfy the orthogonality condition (11.181) only for elem entsϕ∈Wof the subspace. As usual, this only needs to be checked on the basis elements. Su bstituting (11.184) into the semi-weak form of the system, (11.182), produces a linear sy stem of equations of the form /an}bracketle{tw;K[ϕi]/an}bracketri}ht=n/summationdisplay i=1mijcj=bi=/an}bracketle{tf;ϕi/an}bracketri}ht, i = 1,...,n. (11.185) The reader will recognize this as exactly the same finite elem ent linear system (11.163) derived through the minimization approach. Therefore, for a self-adjoint boundary value problem, the weak formulation and the minimization princip le, when restricted to the finite-dimensional subspace W, lead to exactly the same equations for the finite element approximation to the solution. In non-self-adjoint scenarios, theweak formulationissti llapplicableeven thoughthere is no underlying minimization principle. On the other hand, there is no guarantee that either the original boundary value problem or its finite elem ent approximation have a solution. Indeed, it is entirely possible that the boundary value problem has a solution, 12/11/12 629 c/circlecopyrt2012 Peter J. Olver but the finite element matrix system does not. Even more worry ing are cases in which the finite element system has a solution, but there is, in fact , no actual solution to the boundary value problem! In such cases, one is usually tipped off by the non-convergence of the approximations as the mesh size goes to zero. Neverthe less, in many situations, the weak solution approach leads to a perfectly acceptable nume rical approximation to the true solution to the system. Further analytical details and applications of weak solutions can be found in [ 76,188]. 12/11/12 630 c/circlecopyrt2012 Peter J. Olver