Andreasmihalea@gmail.com 875 Informatii 2d Text

Chapter 12 Elliptic Partial Differential Equations and Difference Schemes ‘This chapter îs he first of three chapters dealing with elliptic partial differential equations and finite difference schemes. We start with a survey of important properties of elliptic ‘equations and the effects caused by boundary conditions. Then we show that finite difference schemes have properties analogous to those of the differential equations. The following two chapters are devoted to methods for obtaining the solution of the finite difference schemes, 12.1 Overview of Elliptic Partial Differential Equations ‘The archetypal elliptic equation in two spatial dimensions is Poisson’s equation Mer tur = FO) [EI] în a domain @ as illustrated in Figure 12.1. The Laplacian operator isthe operator on the left-hand side of (12.1.1), and we will denote it by V?, ie. au = 18 (pda Vo rar ar ‘The homogeneous equation corresponding to (12.1.1) i called Laplace’s equation, vu =0, (02.1.2) ‘The solutions of Laplace’s equation are called harmonic functions and are intimately con- nected with the area of mathematics called complex analysis (see Ahifors (2). ‘To determine completely the solution of (12.1.1) itis necessary to specify a boundary ‘on the solution. Two common boundary conditions for (12.1.1) are the Dirichlet ion, în which the values of the solution are specified on the boundary, i.e, u=b onan, 2.1.3) a 312____ Chapter 12. Elliptic Partial Differential Equations and Difference Schemes Figure 12.1 and the Neumann condition, in which the values of the normal derivative are specified on the boundary, i.e. în = on 38, a2.) where 282 refers tothe boundary of 2. Only one boundary condition can be specified at cach point of the boundary perhaps with (12.1.3) specified on one portion ofthe bound- ary and (12.1.4) specified onthe remaining portion. Elliptic partial differential equations together with boundary condition are called boundary value problems. ‘To gsinaphysicll intuitive understanding of 12.11, we may regarditas describing the steady-state temperature distribution ofan object occupying the domain $2. The solu- tion (x,y) would represent the steady temperature ofthe domain 2 with heat sources and sinks given by f(x, >). The Dirichlet boundary condition (12.13) would represent spec- ified temperatures onthe boundary and the Neumann boundary condition (12.1.8) would represent a specified heat flux. In paicuar, the Neumann boundary condition (12.1.4) with ba equal 1 zero would represent a perfectly insulated boundary An important observation concering (12.1.1) with the Neumann condition (12.18) specified on the boundary is that fora solution o exist the data must satisfy the constraint Î[r- le (215) If the data do not satisfy this constraint, then there is no solution. This relationship is called ‘an integrability condition and is easily proved by the divergence theorem as follows. We have, from equation (12.1.1), the divergence theorem, and (12.1.4), that Ille ‘The vector îi isthe outer unit normal vector tothe boundary 382. Theintegrability condi (12.1.5) has the physical interpretation thatthe heat sources inthe region must balance with the heat flux on the boundary for a steady temperature to exist. Also, the solution to (12.1.1) with (12.1.4) is determined only up to an arbitrary constant. This has the physical 12.1 Overview of Eliptic Equations 313 interpretation that the average temperature of a body cannot be determined from the heat fluxes on the boundary and heat sources and sinks alone. Although Poisson’s equation (12.1.1) îs the most common elliptic equation, there are ‘many other elliptic equations that occur in applications. We now define elliptic equations ‘more generally Definition 12.1.1. The general (quasilinear) second-order elliptic equation in two di- ‘mensions is an equation that may be written as A(X, Vuza + WE, Sute +O, YM yy FAO, Ys My Mey My) = FEY) (12.1.6), where a,c >0 and bi < ac. [Notice thatthe definition requires that the quadratie form atx, DE? + 2b(x, DE + ele, ya be positive for all nonzero values of (£, n) for all values of (x, y) in 2. Other Elliptic Equations We shal be primarily concerned with second-order elliptic equations, but there are elliptic equations of any even order. In addition, elliptic equations in three dimensions are very important in many applications, The biharmonie equation in two space dimensions is Viu = (Vu) = Were + Dune + noa =F (2.4.7) ‘There are also elliptic systems such as the Cauchy-Riemann equations ide (021.8) wy + and the steady Stokes equations Su p= fi Vu py = fr, (12.1.9) DE ‘The steady Stokes equations describe the steady motion of an incompressible highly viscous fluid. The velocity field is given by the velocity components (u,v), and the function p gives the pressure field. The biharmonic equation (12.1.7) is used to describe the vertical displacement of a flexible, thin, nearly horizontal plate, subjected to small vertical displacements and stresses on the boundary. 314____ Chapter 12. Elliptic Partial Differential Equations and Difference Schemes ‘The essential property ofthese equations and systems i thatthe solutions are more 4iferentiable than the data, For example, he solution, u, of (12.11) has two more deriva tives than does f. Similariy, the solution, u, tothe biharmonic equation (12.17) has four ‘more derivatives than does the function J. In particular, the solutions to Laplace's equa- tion (12.1.2) and the Cauchy-Riemann equation (12.8) are iniitely differentiable. This ‘property—that the solution is more differentiable than the data and that his gain in difer- entiabliy ofthe solution is equal tothe order of the differential operator—characteizes an ‘equation orsystem of equations as elliptic, (For ystems such a (12.1.9) some care has tobe taken in appropriately defining the order; see Douglis and Nirenberg [14].) The ellipticity ‘of an equation is often expressed in terms of regularity estimates, as we demonstrate in the next section. {AS willbe shown in the discussion of regularity estimates, the ellipticity ofa single equation depends on the nonvanishing of the symbol of the differential operator. More precisely, if P is differential operator of order 2m, then the operator is elliptic if there isa constant co such thatthe symbol of P, denoted by pts, 6), satisfies Ip. 891 2 cog" (21.10) for values of IE! sucienty large. The symbol of a diferential operator îs defined a în Definition 3.1.4, but for elliptic equations the factor of e! is not required, since elliptic ‘equations do not depend on time. We point out that equations such as Weg yy =O donothave the property atit solutions are more differentiable than the data. This equation is the wave equation, discussed in Chapter 8, and it has discontinuous functions in its class ‘of solutions. It does not satisfy the ellipticity requirement (12.110. Exercises 12.1.1. Verify that if (w, ») îs a solution of the Cauchy-Riemann equations (12.1.8), then wand v also satisfy Laplace's equation (12.1.2). 12.1.2, Show that second-order elliptic equations according to Definition 12.1.1 satisfy the ‘condition (12.1.10). 12.1.3. Show that the biharmonic operator în (12.1.7) satisfies the condition (12.1.10), 12.1.4, Show that the elliptic equation with constant coefficients ates + 2bugy + cun + yu + dy + eu = f can be transformed to an equation of the form tye + mn tv BE. where y is Lor —1, using a linear change of coordinates, ie. (E, ) = (x, y)M for some matrix M, and where (Gm) = Aula, het for some constants A, a, and f. 12.2 Regularity Estimates 315 12.2 Regularity Estimates for Elliptic Equations In this section we prove estimates that show how the smoothness of the solutions of el- liptic equations depends on the data. We prove these estimates only for equations with constant coefficients; similar estimates hold for equations with variable coefficients but the techniques used to prove these estimates are beyond this text, We begin with the constant coefficient equation ugg + Ditsy + yy + dute +dguy teu = f (022. for (x,y) € R?, and we study this equation using the Fourier transform. We have the Fourier transform of the solution 61, LEP atest yen yy ae 2 [eta and he Fourier inversion formula wee n= ze ff een acn en db ae as given in Section 2.1. There i also Parseval’s relation, | [montan os ff ace eo at ae ‘Also note that for nonnegative integers r and s, Ife ye op a Rm dx dy= [fice îti, EF a de (1222) = eterne sor aa ae Applying the transform to (122.1), we obtain (ag? + 2bE Le + €83 — iduti — teak — ei a , it = SG By (1223) OF + Dobbs et — tanti 1083 By the requirements that 6? < ac and a and e be positive, according to Definition 12.1.1, we have abt + bbs tc > cole? +62) (2.24) for some constant co, and hence when JE? = 6? +62 > C3 for some value Co, lat} + 2be 12 + 08} — id\by — idota — el > cil? +) (12.2.5) 316 ___ Chapter 12. Fliptic Partial Differential Equations and Difierence Schemes for some positive constant ci. ‘Therefore, from (12.2.3) there is a constant Cj. such that WG Ci for 67 +3 = CA. Then by Parseval’s relation and (12.2.2), popa i ae ay = [| iehepace.en? ati || atata? as ay= ff ete ace.en? aa ae – mezat en? ae, de PI = cati epitet ace ffl EPălei, EP dă dia = ff. BORAG SOP abs dee a 24 apoi feast +r Îl +8Y LAE: ga)? dă de Ei cae” ff, 1c, Gav? de de +c || ten If we use the norms defined by UPFC, B20? dB de ul XL werapur (see Section 2.1, then the preceding estimate leads to Matia = CAUZE + Il). (226) Estimate (12.2.6) is called a regularity estimare. It states that if a solution to (122.1) ‘exists in L?, ie. if lulo is finite and the function f has all derivatives of order through 5 in L2(R), then the function w has 5 +2 derivatives in 12(R). Notice thatthe relation (12.2.4) is essential to proving the regularity estimate (12.2.6). „A curve on which ag? + 2béi£2 + ct? is constant is an ellipse in the (£1, E2) plane. This îs the historical reason for the name elliptic, although now the name is applied 10 more general equations. (See the discussion of the origin of the names hyperbolic and parabolic in Section 8.1.) ‘The property that characterizes an elliptic equation is that the solution ofthe equation more differentiable than the data and that the increase in differentiablity of the solution is equal to the order of the equation. For a second-order equation the property of ellipticity is expressed by the regularity estimate (12.2.6). The biharmonic equation and other fourth- ‘order eliptic equations satisfy analogous estimates, showing hat the solution has derivatives 12.3 Maximum Principles 317 ‘of order 4 more han the data (ee Exercise 122.2). Epic systems, such asthe Stokes ‘equations, saisfy regularity estimates, but the concept of order must be generalized: see Douglis and Nirenberg [14]. If equation (122.1) holds on a bounded domain 8 in 2, we can easly obtain an ‘meriorreglarity estimate ona subdomain SI, whose boundary is contained in the interior of 2. The interior regularity estimate îs Wl ag, = CAR, 21) (IF .a + ulda) 0229 “This has the same interpretation as (12.2.6), but gives estimates only in the interior of the domain. Norms such as |. are defined asin Section 2.1 for integer values of s, but the integration is only over the domain 2. ‘The estimates (12.2.6) and (12.2.7) also hold if the coefficients of (12.2.1) are func- tions of (x, y), as long asa constant co can be found so that (12.2.4) holds for al (x,y). More sophisticated techniques than those used here must be employed to obtain the cs mates when the coefficients are variable. The theory of pseudodifferential operators has. been developed to extend the techniques used hereto the situation when the coefficients. are not constant; see, for example, Taylor [61]. Exercises 12.21. Prove relation (12.2.4) for equation (12.2.1) from Definition 12.1.1 12.22, For a fourth-order elliptic equation of the form aurar + Denys + curs = Ff with a and e positive and with b > —/@6, prove the regularity estimate lunes = Co (Nfs + lao) 12.23. Prove (12.27) by considering the function gu, where $(x, 9) isa smooth eutofT function such that @ is 1 on 1 and 0 on the boundary of ©. Hint: The function du can be extended (0 all R® by setting it to zero off of 8 and gu satisfies a 0 ina bounded domain ©, then the maximum value of u in 8 is on the boundary of 2. 318____ Chapter 12. fliptic Partial Differential Equations and Diference Schemes is theorem can be regarded as an extension to two dimensions of the following result: Ifa function of one variable has a positive second derivative on a closed interval, then that function must achieve its maximum value at the ends of the interval. Figure 12.2 shows a cartoon illustrating the idea that if a second-order differential elliptic operator is positive, then the surface shape is somewhat upwardly concave and the maxima occur at the boundary. On the other hand, if the operator is negative, then the minima occur at rays Figure 12.2. A cartoon illustrating the maximum principle. ‘Theorem 12.3.2. Ifthe elliptic equation ty + Doty Heyy + ditty + day + eu = 0 holds in a domain ©, with a and c positive and e nonpositive. then the solution w(x, ¥) cannot have a positive local maximum or a negative local minimum in the interior of . We prove both of these theorems under the assumption that w is in C3. We prove ‘Theorem 12.3.1 only in the case when Lu is positive, and we prove Theorem 12.3.2 only in the case when e is negative. The proofs for the general case, when Lu is nonnegative and e is nonpositve, require a more careful analysis; see Garabedian [23], Proof of Theorem 12.3.1. If u is any C? function with a local maximum at (ro, 39), then the gradient of w is zero at (9, Yo), ie. ug 20, 30) = My C0. 0) See the illustration in Figure 12.3. By the Taylor series expansion about (xo, yo), we have that (20+ Ar. yo-+ Av) =u (20,30) + (Arta, + ZA av 0 + 092) + O(max(Ax, Ay))® We use the superscript of O to indicate that the functions are evaluated at (xp, yo). Since s(x + Ax, y0-+ Ay) £ uo, vo) forall sufficiently small values of Ax and Ay, it follows that Atul, + 2Aray 10 + Aud, <0. Since this expression is homogeneous of degree 2in Ax and Ay, we have oud, + 2apu?, + Pub, <0 (23.9 for all real values of e and fi 12.3 Maximum Principles 319 ‘We now prove Theorem 12.3.1 forthe case when Lu > 0. Applying (123.1) twice, first with ce = Va and B = b°/Ja® and then with a =0 and f? = — (60); /a°, wehave Lu = ud, +2602, + cul, = 03) cra (2) aa, + (a 2 (pu, =0. Since this inequality contradicts the assumption that Lu > 0, Theorem 123.1 îs proved. O Proof of Theorem 12.3.2. We prove the theorem only for the case when et. ¥) is strictly negative. The case when er, y) is zero at a maximum requires a more careful analysis, and we will omit it; see Garabedian (23) We first conclude from Theorem 12.3.1 that if has a maximum at (xo, yo), then Lu cannot be positive there, Thus we have Luta, yo) = etro: yo)uCo, Yo) 2 0. Since e(xo. Yo) is negative, it follows that u(x, yo) îs not positive at an interior local maximum, Similarly by considering ~u(x, y), we can show that w is not negative ata local minimum. O ‘The maximum principle applied to Laplace’s equation (12.1.2) on a domain has the physical interpretation that fora steady temperature distribution, both he hottest and coldest ‘temperatures occur at the boundary of the region. This means that harmonic functions, solutions of Laplace’s equation, have their maximum and minimum values on the boundary cof any domain. Figure 12.4 displays a portion ofa surface plot of the harmonic function x? — y? illustrating the locations of the maximum and minimum values. For any domain the highest and lowest points of the surface above the domain must occur on the boundary ‘of the domain. Figure 12.5 is contour plot of the same function shown in Figure 124. The ‘maximum principle can be used to prove the uniqueness of the solution to many elliptic equations. 320____ Chapter 12. fliptic Partial Ditferential Equations and Difierence Schemes Figure 12.4. A surface plot ofa harmonic function, Zz Figure 12.5. A contour plot of a harmonic function. Example 123.1. As an example ofthe application of maximum principles consider the equation Mec tu ou f (232) în a domain © with Dirichlet boundary conditions. Assume that here are wo solutions u and v to (12:32) and assume that u is greater than v somewhere în &. Set w= u — vi then w satisties (123.2) except with f equal to ero, and w is zero on the boundary. Since w is postive somewhere in ȘI andiszeroon 36, w must have a postive interior 12.3 Maximum Principles, 321 local maximum. But this contradicts Theorem 12.3.2, and thus equation (12.3.2) has at ‘most one solution. In fact, (12.3.2) does have a solution if © has a smooth boundary, but we will not prove this. 0 Example 123.2. As an example of an equation with a nonunique solution, consider uns busy +202u = 0 0233) on the unit square with u equal to zero on the boundary. It îs es ly checked that w= A sinxx sin xy isa solution for any value of A. Also, equation (12.3.3) with w equal to | everywhere on the boundary has no solution. O Exercises 123.1. Show that he equation War tay = f ‘on a domain 8 with w equal to zero on the boundary has a unique solution, if a solution exists. Hint: For two functions u and v, the function (e” — e*)/(u — v) îs positive. 123.2. Show that if w satisfies the elliptic equation ite + Dbugy + cuyy = 0 on a domain, where the coefficients a,b, and e are constant, then the quantity 12 +42 takes its maximum on the boundary ofthe domain, 12.33. Prove that if w satisfies the elliptic equation of Exercise 12.3.2 on a domain and P îs a positive definite matrix, then the function Flr, y) = PV DT PVUCe, ¥) takes its maximum on the boundary of the domain, ( Vu is the gradient vector of 4 with components 9u/3x and 3u/3y.) 1234. Show by example that if Viu = 0 in a domain 82, then u can have an interior ‘maximum or minimum. Hint: Consider quadratic functions of (x, y). 12.35. Prove that there are no closed contours in the contour plot ofa harmonic function, 322 ___ Chapter 12. Eliptic Partial Differential Equations and Diference Schemes 12.4 Boundary Conditions for Elliptic Equations We restrict our discussion of boundary conditions to second-order equations and to the Dirichlet condition (12.1.3) the Neumann condition (12.1.4) and the more general Robin condition Faw =o. 24.0) ‘The existence and uniqueness of the solutions ofthe general second-order elliptic equation (12.1.6) given boundary conditions of the form (12.1.3), (12.1.4), and (12.4.1) depend on “global” constraints, such as the integrability condition (12.1.5). For certain equations, especially (12.1.1), on domains with smooth boundaries, the existence and uniqueness questions have been answered. With the Dirichlet boundary condition (12.1.3), Poisson’s equation (12.1.1) has a unique solution, and with the Neumann condition (12.1.4) there is a unique solution, up to the additive constant, if and only if the integrability condition (12.1.5) is satisfied. See Section 13.7 for more on solving boundary value problems with the Neumann boundary condition. General statements can be made about the local behavior of solutions to (12.1.6) given the different types of boundary conditions. If a Dirichlet boundary condition is enforced along a smooth portion of the boundary, then the normal derivative at the boundary will be as well behaved as the derivative of the boundary data function in the direction of the boundary. If the boundary data function is discontinuous, then the normal derivative will be unbounded at discontinuities. Example 12.4.1. As an example of this last statement, consider Laplace’s equation in the upper half-plane, ic., y > 0, with 0 ite > eo-li ite co Asani oo where @ îs the angle in radians measured from the x-axis. The normal derivative of the solution along the houndary is the derivative with respect to y. We have autv.y)_ 1 Ur ay TE OP Atthe boundary we have the normal derivative which is unbounded at the point of discontinuity inthe boundary data function. The normal derivative is unbounded at the point where the tangential derivative ofthe boundary data is, unbounded 12.4 Boundary Conditions 323 Notice also that în the interior of the upper half-plane, the solution and its derivatives are all well-defined and smooth functions. The behavior ofthe solution near a point on the boundary is primarily influenced by the boundary condition and data near that point. The conditions and data at other boundary points have less effect the further they are away. D If either the Neumann or Robin conditions are enforced along a smooth boundary, then the solution will be differentiable up to the boundary and the first derivatives will be as well behaved as the boundary data function ‘Example 12.4.2, As an example, again on the upper half-plane, consider the boundary data 0 itr>0, bi? ifx <0 Laplace's equation has the solution FED) u(x, y) = Ar cos 0 using the polar coordinates of (x,y). The derivatives are given by „ du /ax, has the same qualitative behavior as A serious difficulty occurs at points on the boundary where the boundary conditions Change from Dirichlet to Neumann or Robin type. Example 12.4.3. Consider Laplace's equation in the upper half-plane with the boundary conditions (4,0) =0 forx > 0 and uy(r,0)=0 forx <0. ‘This problem has asa solution us.) = 7"? sin jo, (12.42) Note that w and its frst derivatives are în L2(82) for any bounded domain © in the upper half-plane whose boundary includes a portion of the real axis around zero. This function 1, however, does not have second derivatives in L2(€2) because of their growth near the origin. The first derivatives are also unbounded, but are in L2(0). 324 __ Chapter 12. Elliptic Partial Differential Equations and Difference Schemes Example 1244. Similar difficulties arise at reentrant comers. Consider the domain containing the points (7,8). inpolarcoordinales, with 0 direction, A schematic of the is shown in Figure 12.7. The standard central difference approximations for the second eivatives ead tthe diference formula Sven + 82m = f nese, pi etnia Feta tee + Bent 486m) = Sem 2s.) ‘The difference operator on the left-hand side of (12.5.1) is called the five-point (discrete) Laplacian. We will use the symbol VE to refer tothe discrete five-point Laplacians, i, Vpast 482 Figure 12.7. The five grid points involved in the five-point Laplacian. 326 ___ Chapter 12. Elliptic Partial Differential Equations and Difference Schemes ‘We begin our study of finite difference schemes for elliptic equations by obtaining error estimates forthe solution to (12.5.1). In the next two chapters we consider methods for solving equations such as (12.5.1) The Discrete Maximum Principle ‘We can prove a maximum principle for the discrete five-point Laplacian that is analogous 10 that for the differential equation, ‘Theorem 125.1. Discrete Maximum Principle. If Vu > 0 on a region, then the max- imum value of v on this region is attained on the Boundary. Similarly, if Vw <0, then the minimum value of v is attained on the boundary. Proof. We prove the principle only for the case when Ax = Ay. The condition Viu = 0 is equivalent to Vem £ (vert + Ven bem + Vem i.e, team ithe interior of the region is less than or equal o an average of its four nearest, neighbors. This easly leadsto the conclusion that an interior poin canbe a (local) maximim only if its four neighbors also have his same maximum value and thatthe inequality is actually an equality. This argument then implies that at all grid points, including the boundary points, v must havc the same value, So either there is not a maximum value in the interior or all points have the same value, This proves the principle when Vv = 0 ‘When Viu = 0, by considering Vi(-v) > 0, this ease reduces othe previous case. ‘This completes the poof of the discrete maximum principle. Dl “The maximum norm on a region $2 îs defined by Wwe = Bulle. = max [vem ‘where the maximum is taken over all points in the region. “The chief tool in our error estimates is the following theorem, ‘Theorem 12.5.2. If vim is a discrete function defined on a grid on the unit square with Vem = 0 on the boundary, then 1 Wlso = gIV ul: (252) Proof, Deine the function fe în he interior ofthe uni! square by Vivem = Sem ‘Then, obviously, Wl < Vio < WFlee- (253) 12.5 Schemes for Poisson's Equation 327 ‘To prove the theorem we consider the function wen = [lee 4) + mY] and note that w îs nonnegative and vjw=1 ‘Thus from (12.5.3) we have Vie = Flow) < 0. By the discrete maximum principle and this inequality the function v — J floow has its ‘minimum on the boundary of the square, i.e. =H filo loa < ve = Wflooweom < Yom where Ilea în the maximum value of jw for grid points on the boundary of the Sin, fm (253) we alo cain Vito + feo) > 0, and by the maximum principle, sem £ dem + Mf looweon = Wf ello ‘The value of jelaa. îs 1/8, and so the preceding two inequalities for vem give 1 Wllso = g/l = VE Which proves the theorem. D “Theorem 12.5.2 leads tothe error estimate inthe maximum norm forthe solution of (12.1.1) as approximated by (12.5.1), ‘Theorem 12.5.3. Let u(x, y) be the solution to Vu = f onthe unit square with Dirichlet boundary conditions and et vm be the solution to Vv = f with Vem = uta, Yn) n the boundary. Then ue — Uo < ch2Notulleo» where \iO*ullso is the maximum magnitude of all the fourth derivatives of u over the interior ofthe square. 328 ___ Chapter 12. Eliptic Partial Differential Equations and Difference Schemes Proof. By using the Taylor series for the central difference approximation to the second derivative (see Section 3.3), we have that Viu = f+ 00), where the OC?) terms are hounded by CHA ul for some constant C. Thus AV; — vos = CHP Ulla and u—v is zero on the boundary. Together with Theorem 12.5.2 this estimate proves the theorem, D The Nine-Point Laplacian Another very useful scheme for Poisson's equation (12.1.1) is the fourth-order accurate nine-point Laplacian. To derive this scheme we approximate (12.1.1) by aaa Aa au 4 (+ 8) su (14 ai) ur 018%), which gives Ay a a2 2) zu (14458) stu (1 +See a2) [i ag 4 = (14 Spa) (4 2) FF. - [ + heard + ava] F+0(04 Rearranging this expression we have the fourth-order accurate scheme Vout fh (ax? + ay?) tat = f+ qh (Atat + 22) 7, Inthe case with Ax = Ay = h, this scheme can be written 1 2 10 +3 (Vestn + etn ema + em) = FY (254 R 73 (ferim + Simin + Sem + fom-t + 8fem) 12.5 Schemes fer Poisson's Equation 329 Table 12.5.1 Comparison of second-order and fourth-order schemes. Second-order Fourth-onder a Error Onder Error Order D100 279-5 340-9 0.050, 701-5 1 3.8510 tor 0025, 175-6 2.00 66-11 4.00 ‘This scheme is due to Rosser [54], and it satisfies maximum principles and error estimates similar o the standard five-point Laplacian; sec Exercise 12.5.6, Table 12.5.1 shows the resulis of computations employing both the second-order accurate five-point Laplacian (12.5.1) and the fourth-order accurate nine-point Laplacian (12.5.4) applied to Poisson’s equation on the unit square. (The results were computed using a preconditioned conjugate gradient method, which is discussed in Section 14.5.) ‘The exact solution is given by u = cosxsiny and f = —2cosx sin y. The second and fourth columns give the errors for the two methods, measured in the L? norm due to the difference between the finite difference solution and the solution ofthe differential equation. ‘The third and fifth columns display the order of accuracy for each method, as computed from the approximation that eh) = ch", where e(h) is the error. Thus og(e(h)/e(ha)) Yog(hi/ a) for two successive errors duc to grid spacings fy and hp. The five-point Laplacian (12.5.1) is obviously second-order accurate, and the nine-point formula (12.5.4) i obviously fourth- order accurate. Notice that the error forthe fourth-order scheme with A equal to 1/10 is much smaller than that of the second-order method with equal to 1/40. The gain in accuracy of the nine-point formula is significant compared to the slight increase in work that it requires, ‘Schemes for the general second-order elliptic equation (12.2.1) need not satisfy a ‘maximum principle. For example, ifthe mixed derivative term in (12.2.1) is approximated by Sortorv, then the resulting scheme does not satisfy an obvious maximum principle. If the coefficient x, y) îs positive, then the second-order accurate approximation aac bca SZ Crtdrs +8145) will satisfy a maximum principle if b îs not too large compared with the coefficients. a and e. Schemes that do not satisfy a maximum principle may have solutions and satisfy error estimates such as in Theorem 12.5.3; however, the proofs are not as simple as those just given. We do not need a maximum principle to hold in order to use the scheme. 330____ Chapter 12. Eliptie Partial Diferential Equations and Diference Schemes Regularity Estimates for Schemes We can prove discrete regularity estimates for schemes for elliptic equations, as is done for the differential equation, For example, the scheme atu + Dob oy + e820 + dar + doy + ev =f has the symbol 1 1 sin? $hgy + 2b sin $y sin FhE cos $e, cos $hE, + csin? és a pb sin hey, simhéy et hag +e. +id, ‘The analogue ofthe estimate (12.2.5) for the symbol is, with @ = h&, and @ = hé>, asin? 40 + 2bsin Jo sin Jpcos jocos Jp +csin? $4 sin@ sin |gasin_34 2 is 3900s 3 26 1a, oe ide ae _ > eosh- sin? Yo + sin 40, which holds for some postive constant cy, when Mis small enough and when 07 +63 > ICE for some value Co The interior regularity estimate that follows from this estimate îs Wl s42.0, SC (Ilia + Ilo): (12.5.5) ‘The discrete interior regularity estimate can be used to prove the following result. ‘Theorem 12.54. If the elliptic equation Lu = fis approximated by the scheme Lu = Ryf ona domain 2 such that Wau = Re Lut s-2a = coh us (125.6) and Iu = vityoa < cub fully 0257 and 2 is contained in ©, then NS} u — 8% vlla..a, = coh" lls, where cy depends on the distance between 21 and 39. 12.5 Schemes for Poisson s Equation 331 Proof. The discrete function 1 — v satisfies the scheme Luta —v) = Law — Ruf = Lou — Rau and, by (125.5), Wu VIR aa, £ Ci-a (ULmu = Redo, = Ce-a al, from which the theorem follows. Dl ‘This theorem shows that if the function f is smooth enough, then the divided di ferences of v approximate the divided differences of u to the same order that v itself approximates u. This is demonstrated in Exercises 13.5.3 and 13.5.4 „At frst look, this is a very surprising result. Notice that in genera, if a discrete function vz.» is an approximation to u(x, y) of order 7, then a+ IM vlog) veti — Venton _ dure: Ym) 2h or because the divided difference divides the error by a factor of A. However, Theorem 12.5.4 implies that when vn and ue, ») are solutions of elliptic equations, then the error term can be OCH) rather than O(A7-1) ‘The reason is tat for eliptic equations the error is smooth, and a divided difference ofthe error i again a smooth function. We look more closely at how Theorem 12.5.4 can be used to obtain approximations to derivatives of solutions of eliptic equations that are of the same order of accuracy as the solution itself. Suppose the solutions ofa founth-ordce accurate scheme satisfy estimates (12.5.6) and (12.5.7) with 7 equal to 4. Then l0- Ea) Suatu), = 00% De oa, 2, Fo) aa 4) (see formulas (3.3.6) and (3.3.7). Note, however, that Iu — ull.o.0, = OOF), since 6? is only a second-order accurate approximation to 32. These results also apply to ‘equations with variable coefficients; see Frank [21] or Bube and Strikwerda [7], By comparison, such resulls do not hold for solutions of hyperbolic problems. If vf, îs a solution to a second-order accurate scheme for a hyperbolic equation, such as the one-way wave equation (1.1.1), then, in general, Su, is only a first-order accurate approximation to the ist derivative of u. ‘Also, as was discussed in Section 12.4, the solutions of elliptic equations have, in ‘general, loss of regularity at certain boundary points. The implication for finite difference solutions is thatthe eros willbe largest at these points, such as reentrant corners. Shown in Figure 12.8 îs a contour plot of the error for the solution of Laplace's equation with the boundary data given by the exact solution u(r, 8) = 727% sin($6). As the contour plot shows, the error is concentrated a the reentrant comer. + oh 332 ___ Chapter 12. Elliptic Partial Differential Equations and Difference Schemes Figure 12.8. Contour plot of the error at a reentrant corner Exercises 128.1. Prove the discrete maximum principle on the unit square forthe case with Ax not equal to Ay. 1252. Show that on a domain tha is contained in a square of side d the analogue of the estimate (12.5.2) is eo Melo = Vielen: 1283, Prove the equivalent of Theorems 12.5.1, 12.5.2, and 12.5.3 for the nine-point scheme (12.5.4). 125.4, Consider the domain givenby —I < x < |, —I La wy ar no (2) x 5 Mao ag Since uoj is independent of j—call this value up—we have , 2 w= 5 Dy -10(5) (1263) Using this formula preserves the second-order accuracy of scheme (12.62). ‘For parabolic and hyperbolic equations on a disk, a procedure analogous to the one that gave rise to (12.6.3) can be used to give accurate difference formulas atthe origin. Exercises the discrete maximum principle holds for finite difference scheme (12.6.2) ith formula (12.6.3) used at the origin, 126.2. I we denote the finite difference operator on the lefi-hand side of (12.6.2) by V3, show that the estimate Wise = ANF holds fora disk of radius 1, where the formula (12.6.3) is used atthe origin. 12.7 Coordinate Changes 335 12.7 Coordinate Changes and Differences Frequently partial differential equations must be solved on domains that are not rectangles, disks, or other nice shapes. Sometimes itis possible to change coordinates so that a con: venient coordinate system can be used. To illustrate the techniques and the difficulties, we will work through a relatively simple example. It is not hard to come up with much more difficult examples. oe Pt \ \ \ \ LII Figure 12.10. The trapezoidal region and grid. ‘We consider Poisson’s equation on the trapezoidal domain given by O < x < 1 and 0 e Aem. cet.) 1 1 sh Meme i Amen Ames when these elements are defined. All other elements are 0. Ifthe grid spacing is given by în = 1/N, then the matrices have order K equal to (N ~ 1). We nov, consider the splittings corresponding to the two methods that we are studying in this section, Notice that the B matrix multiplies the unknowns of iteration k + 1 and the C matrix multiplies those of index k. For the Jacobi method we see that B=1 and C=L+U, 0332 ‘The splitting for the Gauss-Scidel method depends on the order of the unknowns. We take the order in which the unknowns are updated to be the same as that used in the vector x. With this proviso, the Gauss-Seidel method has the splitting B -L and C=U. 0333) “The matrix decomposition (133.2) for the Jacobi method isa restatement of (13.1.3), which shows that the updated variables, those multiplied by B, are only the diagonal 346 Chapter 13. Linear iterative Methods elemens. The variables evaluated at step k in formula (13.1.3) are those coresponding to the off-diagonal clements of the matrix. Similarly, the decomposition (13.3.2) forthe Gauss-Seidel method isa restatement of (13.14) in which the variables evaluated at step £-+ 1 are those comesponding tothe elements of the matrix onthe diagonal and below: Notice thatthe matrix B, being a lower triangular matrix, is easy invert. îtisimponan to realize that in the actual implementation ofthese methods in acom- pater program, we do not store the matrices A,B, and C. They ae all quite sparse and itis very inefficient to store them as matrices. The matrices are useful inthe analysis, but the implementation canbe done without explicit reference to them. Thais, a computer im- plementation should not havean (N — 1)? x (N — 1)? array for storage of these matrices. Instead the implementation should use a form such as (13.1.4, n which ony the current values of vf, are stored, There is no reason to store other arrays. Analysis of the Jacobi Method ‘To determine the spectral radius ofthe iteration matrix for each of these methods applied to the five-point Laplacian, we first find the eigenvalues and eigenvectors of the iteration matrix for the Jacobi method (13.1.3). That is, we must find a vector 3 and value yx such that pă = (E+ ODD. IE we represent 5 asa grid function with indices from 0 to N in each direction, with the indices 0 and N corresponding to the boundaries, we have tem = 4 (Venta + Vem at + Vegtm + em) (3.34) forall interior points and vz. equal to zero on the boundaries. ‘As mentioned after equation (13.2.1), it isimportant to make a distinction between the jgenvector 3, which has unknowns corresponding to the (N — 1)? interior grid points, and the grid function v, which has (N + 1)? values corresponding to both the interior and ‘boundary points. Because we specify tha the boundary values of v are zero, we can write the simple formula (13.3.4). The equations for 0 are different than (13.3.4) if (e, m) is ‘next toa boundary, in which case atleast one ofthe terms on the right-hand side of (13.3.4) would not be present, The use of the grid function v in place ofthe eigenvector 3 allows for a simpler way to write the equations. Since L + U isan (N ~ 1)? x (N 1)? matrix, there should be (N — 1)? eigen- values and cigenvectors, ‘Comparing equation (13.3.4) with (13.2.1), we see thatthe eigenvalues of the Jacobi ‘method are related to those of the Laplacian by So the eigenvalues are a oo()] aaa 13.3 The Jacobi and Causs-Seidel Methods 347 for | < a,b < N — 1. By equation (13.2.7) the eigenvectors are given by (13.3.6) “This gives all (N — 1)? cigenvalues and eigenvectors for the Jacobi iteration matrix. See also Exercise 13.3.1 From the formula (1 ), we see that = = cos (BTC) = pL + U) = cos = Since p(L + U) is ess than 1, the Jacobi method will converge; however, since p(L. + U) is very close to Lie. x Es ose 2, ve see that the convergence will be slow. “The relationship 4.%-*-"—» = —y% shows that the nonzero eigenvalues occur in ‘an eigenvalue, then —1 is also an eigenvalue. Notice also that the for a between | and N — 1 are all equal to 0 and these are the only eigenvalues equal to 0. Thus there are N — 1 eigenvalues of L+U that are zero, and ‘consequently there are (N — 1)(N — 2) nonzero cigenvalues. Analysis of the Gauss-Seidel Method We now consider the Gauss-Seidel method. An eigenvector of the iteration matrix (= DAU with eigenvalue A satisfies AU = Ls or, for the grid function vem, we have Avem = 4 (Rtn + Ave m-ai + Yeon + Vei 133.7) forall interior points and ve, equal o zero on the boundaries. Notice thatthe coeficient 2 in(13.3.7) muhiplies only the variables with superscript of k + 1 inthe formula (13.1.4). In the form (13.3.7) the formula is rather intractable; however, there isa substitution that reduces the analysis of this case 1o that of the Jacobi method. If we set Oem (1338) for each nonzero eigenvalue 2, we obtain, after dividing by X(e+m+ 2, Dem = 3 (Went ue mi iti + Meme) (13.3.9) 348 Chapter 13. Linear iterative Methods By comparing (13.3.9) with (13.3.4), we sce that the nonzero eigenvalues A. of the Gauss-Seidel method are related to the eigenvalues pu of the Jacobi method by amr (oarba = i? = 5 (cos = + cos (343.10) aw N In particular, PU — DU = pL + UP, Which shows that the Gauss-Seidel method converges twice as fast as the Jacobi method for the five-point Laplacian. ‘The eigenvalues of the Gauss-Seidel iteration matrix from equation (13.3.10) give ‘only (N — 1)(N ~2)/2 eigenvalues correspondingtothe (N ~ 1)(N ~ 2) nonzerocigen- values of the Jacobi iteration matrix. An cx for the Gauss-Seidel method shows that they are i walues are zero, and they are not semisimple. (See “Appendix A for the definition of a semisimple eigenvalue.) ‘An alternative way to describe the preceding analysis is to consider the equation dat (0 — Dil=0 for the eigenvalues of the Gauss-Seidel iteration matrix. We have the relationship 0 = deta —( — DY1U] = deat — AL — Uy den(d — Ly". ‘The value of det( ~ L)~! is 1, since L is strictly lower triangular. We next transform the matrix AJ —AL —U by a similarity transformation using the diagonal matrix S, where the (£,m)th entry on the diagonal îs ACE. (Recall that the rows and columns of L., U, and 5, are indexed by the ordered pairs of integers corresponding tothe grid i SOI AL USM ALA) corresponding o (1339). Thus dat = AL =U) = det = AL + UI = ANU geti! = (+ UI aie ȚI ate IE) PUII e awh, 2z04sN-1 13.3 The Jacobi and Gauss-Seidel Methods 349 Wherein the last product we used the facts that ye" i zero foreach a and N-%:N-P = =n, This last formula confirms our previous conclusion that there are N(N ~ 1)/2 zero ‘eigenvalues and shows that (133.10) gives the (N — 1)(N — 2)/2 nonzero eigenvalues. ‘An examination of why the substitution (13.3.8) works shows thatthe updating of values in the Gauss-Seidel method can be organized cither în the standard lexicographic ‘order or in the order of increasing values of £+ m. When one updates a value at a grid point with indices (€, m), the computation involves only points of lower value forthe sum. ‘of the indices, the points with “new” values, and points of larger value forthe sum ofthe indices, the points with “old” values ‘The Jacobi method can also be regarded as solving the heat equation Uy = Men + yy Until a steady-state solution is reached using forward-time cental-space differencing and 4)? In general it seems that finding steady-state solutions by solving the correspond- ing time-dependent equations is less efficent than using special methods fo the steady-state ‘equations. The Gauss-Seidel method can be regarded asa finite difference approximation for the time-dependent evolution for the equation My 3 Mex + yy — tura Fum) where Ar = ȘI? and e = AA. The equation shoud be disrctized about (1.9) equal to ((n-4+ Pat, th, mh) oobain (13.14). Methods for Diagonally Dominant Matrices ‘We now state and prove a theorem about the Gauss-Seidel and Jacobi methods for the Class of diagonally dominant matrices. Many schemes for second-order elliptic equations, including the five-point Laplacian, give rise to diagonally dominant matrices. Definition 13.3.1. A matrix is diagonally dominant if Dial < laut 033.1) in for each value of i. A row is strictly diagonally dominant if the inequality in (13.3.11) is a strict inequality and a matrix is strictly diagonally dominant if each row is strictly diagonally dominant. [By a permutation of a matrix A, we mean a simultaneous permutation of the rows ‘and columns of the matrix; ic., aj, is replaced by aat.a(j) for some permutation a. Definition 13.3.2. A matrix is reducible ifthere isa permutation a under which A has the structure Ao " : 133.12) (a 4) map where A, and Az are square matrices. A matrix is irreducible if tis not reducible. 350 Chapter 13. Linear terative Methods For an arbitrary matrix A the Jacobi iterative method for equation (13.1.1) îs sith = DCD — A)! +b) =~ Dat + (033.13) ‘where D is the diagonal matrix with the same diagonal clements as A. If A is written as A=D-L-U, where L and U are strictly lower and upper triangular matrices, respectively, then the Gauss-Seidel method for (13.1.2) is (D= Dx = ui +b, (33.14) Notice that the diagonal dominance of a matrix is preserved if the rows and columns of the matrix are permuted simultancously. The Gauss-Seidel method is dependent on the permutations of the matrix, whereas the Jacobi method is not, and a matrix is reducible if in using the Jacobi method itis possible to have certain components of x be zero for all values of k while x° is not identically zero (see Exercises 13.3.4 and 13.3.5). ‘Theorem 133.1. If A is an irreducible matrix that is diagonally dominant, with at least ‘one row being strictly diagonally dominant, then the Jacobi and Gauss-Seidel methods are converge! Proof. We prove the theorem only for the Gauss-Seidel method; the proof forthe Jacobi method is easier. Our proof is based on that of James [31]. We begin by assuming that there is an eigenvalue of the iteration matrix, 2, that satisfies [Aj > 1. Let x be an eigenvector ofthe iteration matrix with eigenvalue A, and we normalize x so that Ute îs 1 ‘Let x, be a component of x with [xj] equal 10 inequal 1 A Day + DS ayy i eră then we have the series of [Alla tts SI lait + ais bay Zits Ewe aaa e iii) își SAD tal < Pilat ia Since the first and last expressions are the same, each inequality in the preceding sequence must be an equality. This implics that foreach j, either [xj] is 1 or aj, is zero. ‘This conclusion follows for each i with [xj] equal to 1 13.4 Analysis of Point SOR 351 If-we permute the indices of A so that the components with |xj| equal to 1 are placed first and the others, for which ja, is zero, are last then the structure of A is of form (13.3.12). Since A is ireducible, we conclude that Îx,| is 1 for each value of j. By choosing arow thats strictly diagonally dominant, the last inequality of (13.3.15) îs then a strict inequality, which leads to a contradiction. This implies thatthe assumption that A satisfies [Al > 1 is false. Therefore, [AI is less than 1 forthe iteration matrix, and the Gauss-Seidel method is convergent. Dl Exercises 13.3.1. Verity by direct substitution thatthe eigenvalues and eigenvectors forthe Jacobi iteration matrix are given by (133.5) and (13.3.6), respectively 1332. Determine the eigenvalues of the Jacobi iteration matrix when applied to the agonal” five-point Laplacian scheme given by Rm fem (133.16) (04st Vegan t+ Meta + Vent duza) on a uniform grid with Ax = Ay =h. Hint: The eigenvectors for this Jacobi method are the same as for the Jacobi method for the usual five-point Laplacian. The eigenvalues, however, are different. 13.3.3, Verify that zeros not semisimple eigenvaluc ofthe iteration matrix forthe Gauss- ‘Seidel method for the five-point Laplacian on the unit square. 13.3.4, Show thatthe Jacobi method (13.3.13) is not affected by a simultaneous reordering of the rows and columns of a matrix, whereas the Gauss-Seidel method (13.3.14) is affected. Note that such a permutation is equivalent to applying a similarity \ransformation using a permutation matrix P to the matrix A resulting in the matrix PAP™! 133.5. Show that a matrix is reducible ifin using the Jacobi method, itis possible 10 have certain components of x* be zero for all values of k while x° is not identically zero (see Exercise 13.3.4), 133.6. Show that the matrix forthe five-point Laplacian on the unit square is irreducible. 13.4 Convergence Analysis of Point SOR We now analyze the convergence of SOR for the five-point Laplacian as given by (13.1.5). We have to determine the splitting matrices B and C. As before, B multiplies the un- knowns titeration £ + | and C multiplies those a iteration k. We also have the condition that A = B= C = 1 – L~ U. After rearranging the formula (13.1.5) and dividing by c, ‘we obtain thatthe spitting is given by B