### Create a StudySoup account

#### Be part of our community, it's free to join!

Already have a StudySoup account? Login here

# 641 Class Note for MATH M0050 at Purdue

### View Full Document

## 17

## 0

## Popular in Course

## Popular in Department

This 562 page Class Notes was uploaded by an elite notetaker on Friday February 6, 2015. The Class Notes belongs to a course at Purdue University taught by a professor in Fall. Since its upload, it has received 17 views.

## Similar to Course at Purdue

## Reviews for 641 Class Note for MATH M0050 at Purdue

### What is Karma?

#### Karma is the currency of StudySoup.

#### You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!

Date Created: 02/06/15

The Numerical Approximation of Solutions of Partial Differential Equations of Elliptic and Parabolic Types Jim Douglas Jr SER Visiting Professor of Mathematics University of Wyoming Laramie WY USA University of Wyoming Spring Semester 2009 Lecture 1 DJ Two Point Boundary Problems Throughout most of the lectures on two point boundary problems we shall consider the problem given by Lu auxxcu f XE01 u g XEB01 where for some oz gt 0 aX 2 oz and CX Z 0 At times we shall assume that CX 2 oz A Finite Difference Method h1N17 xih xih Ii Xiilyxih 1109 TJi The object is to approximate u i 0 7 N by VI7 i 0 7 N where W is determined as the solution of a linear algebraic system of the form ik ik LhWi ZBIJMj Z Yij 177N7 jiik jiik W0 5 07 WN gm Let m 11007 Ogn grnW npO39 From Taylor series expansions 1 ww E 77 710Wkoohk 1 k2 or 3 1 1 wXI EWHH 7 11739 wmjhz 011114ooh37 1 1 dim1 117397 217 mjh200 14mh31 Set 61e671e 8iie X i h ul1 LII 7 X i h LII ul71 7 X i ul1 171 7 1 6Xa6yu p aiui1 7 Lli 7 ai7lLI39 7 LI 1 Apply these relations to 11 aux and 11 u 7auxxi 6xa6fui 5h M S Miiaii3ooiiuii4ooh2 Define Lh by LhW39 76Xa67Wi i cW7 1 Then Lhuf i c7 177N717 U0 g07 LIN gN7 The finite difference equation LhW 17 W0 g07 WN gN Existence and uniqueness of a solution will follow from either of the convergence arguments to come A Maximum Principle Convergence Proof Assume here that CX 2 oz gt 0 and let z u W39 Then LhZ39 E 177N 17 20 ZN0 Either z 0 i 0 7 N in which case the solution of the difference equation coincides with that of the differential equation on the nodes X or there must exist an integer k 1 S k S N 1 such that miaxz zk If zk gt 0 then the piecewise linear interpolant Wh of Wi is concave at in 6Xa6ywk 07 since zk 2 zkil So the error equation 6xa6fzk ka 6k implies that Mk S ka S 6k S Miiaii3ooiiuii4ooh2 If zk lt 0 then Wh is convex at Xk and azk 2 Miiaiisooiiuii4ooh2 Thus if CX 2 Oz and ax gt O on I 1 miaxiu 7 Wii MHaiisooiiuii4ooh2 Theorem Assume that ax gt O and CX 2 04 gt O forx E I and that Ha 3 00 is finite IfHu 4 00 is finite then Huiwh 2 400h 7 1 000 S WHE 3mm 04 Where u is the solution of the two point boundary problem and Wh is the piecewise linear interpolant of the solution of the finite difference equation A Discrete Energy Estimate Proof Let L2 f fgfx2dxltoo H10 f1 Hin Hng llfxllg lt 00 H30 f 6 H10 1 0 1 0 with associated inner products and norms given by f7g0 fgdx7 f7g1f7g0 flt7gX07 and We W02 x fx2dxgt27 ml m2 yam 7 l Hle Hme HfH3 W A The Poincar inequality states that HfHo Owl C0117 f E H6I Denote spaces corresponding to L2I and H1I by 2 and hl with inner products and norms respectively being given by N71 1 1 f7goh E figih 5f0go ngNh7 iifii0h f7 08w i1 N fig1h fgohZ6yf6ygih i1 l figoh5y 6ygoh7 iifii1hf7fih39 Denote the hl seminorm on hi1 f E h1 f0 fN O by mm maxim Lemma Let f E htl Then its 2 and I norms are bounded by its h1 seminorm Wow S Wm Hmem S Wm Moreover if g 6 htl the discrete Green39s formula N 6Xa6yfg Egg 6715 67g h i1 holds Proof Express f in the form Then a 20 h x MM j1 and the first two inequalities follow The third conclusion can be shown by summation by parts but perhaps a more enlightening approach is to recognize the relation between 6Xa67fg and a corresponding L2 inner product Extend 3171 to the piecewiseconstant function AX by setting 2 AX 37 X E X1X7 I391N7 2 I and let F and G be the piecewiselinear functions interpolating and g at the nodes Then 6Xa6yf g AFXXG dx 7 AFXGX dx I I N 7 Z ai76yfi67gi hi i1 as was to be shown An immediate corollary is that llfl m f 6 hi 1h 2lf so that the seminorm on hi1 is equivalent to the full norm there The 2 convergence argument is easy Test the error equation against 2 and sum over i 1 N7 1 and apply Green39s formula 3 6y22 h l 62 2 6 2 7 N 1 E 1 I Here we need assume only that aiilzagt0 and L720 Vi 2 7 Then by Cauchy Schwarz alzlih S llEllohllleoh S llEllohlZl1h7 so that lZl1h S a lllEllow S a 1MallUll4ooh2 Theorem Assume that Ha gm and HuHMX are finite Then Hui WHOM Hui WH0h u7 Wm S MH3H3ooHUH4ooh2 The analysis above for error has been presented in essentially the simplest possible way consequently the theorem is not close to best possible Since the primary aim of these lectures is to treat a variety of finite element procedures for partial differential equations we shall not pursue deeper results here However it is valuable to consider briefly some extensions and possible improvements related to the finite difference procedure Inclusion of a first order term Consider the two point boundary problem Lu7auxxbuxcu f XEIO17 u g7 X613I017 where Redefine Lh as LhW 76Xa6yw biQXW C W 7 177N 1 where 1 8xWi Wi1 W121 The finite difference equation Then with 6 O u kmhk l k 2 or 37 Lh ifi5i7 177N17 U0 g07 N gN7 Again with z u 7 W LZ39lt 397 l39l77Nl7 20 ZN 0 Rewrite Lhz l l LhZi plhg Ehbilwiu 34 ai7Wi 1 ai EhbiWi71l CiWi 6i If Mh lt a then the maximum principle holds and z 9h2 as before We will defer the 2 argument as it is essentially the same as for the Galerkin method Variable Grid Spacing In many practical problems the use of a regular grid frequently is inappropriate Consider the irregular grid 1 O X0 lt X1 lt X2 lt lt XN 17 hi X7X17 xiii XltlX17 and with 1544 fX 7 B 07 i extend the definition of the second order difference operator to be 7 2 Wi1 Wi Wi W121 6Xla6ywli 7 hi hi1 ai hi1 i ali hi gt 39 The finite difference method is defined as before LhW39 76Xa6yw l cW 1 W0 g07 WN gN A simple but tedious Taylor series expansion argument shows that LhUi i7 l il S Mlllal 2Oollul 3m maxhj7 i1N71 J The maximum principle convergence analysis is essentially unchanged except for the loss of accuracy caused by the poorer estimate of 8 Only relatively minor modifications in the definitions of the terms entering the 2 argument are required in order to carry that proof through as well Hence it can be shown that llui Wllooh llui Wl 0h l lUW1h Mllal zoolllul 300 mjaxhj Since the point of using a variable grid is to localize the difficulties in approximating u the global arguments and global error estimates above fail to indicate the true value of varying the grid however we shall defer discussing arguments that localize the regularity requirements to the treatment of finite element methods Discontinuous coefficients Consider the simple example when 31 O Xlt 7 ax N l 32 XltX 1 For the original differential equation written in divergence form consistency at the point X Y requires that lim ux lim ux and lim vx lim vx7 xi XL xi XL where VX 3XUXX represents the flux associated with the differential operator If X1 lt x lt X then introduce a Lagrange multiplier to replace the nonexistent value Wm this multiplier serves to enforce the continuity of the approximate solution at x Then approximate the left and right fluxes V and Vi at x by W V7 731N77 VJr 732 i X 7 Xiil X 7 X ie instead of thinking of the approximate solution as being linear between X1 and Xi consider it to be a broken line maintaining the continuity of the flux There are two ways to impose this concept One would be by adding the point i to the grid and require that with V and VJr as defined above VJ V For the two point boundary problem this is feasible though not altogether convenient however for several space variables it could cause a considerable modification of a computer code based on an orderly grid Alternatively we can eliminate algebraically by finding an effective value aeff iii of the coefficient a such that 2 Set i 7 X1 39yh and solve for by equating V and Vi W 1 W 31 w a217 Wh implies that The effective conductivity coefficient is given by a lt31 32 gt71 31 32 3132 1 i i 7 em v 17v v 17v 17v31v32 which is more easily recognized as the length weighted harmonic average of 31 and 32 h h 1 h L aeffi7 31 32 This can be generalized to treat variable ax discontinuous or not on I by taking The definition above for an effective conductivity coefficient extends to multiple dimensions This evaluation will be seen later to correspond to the use of a mixed finite element or finite volume method Now having found aeffIil replace aiil in the definition of 2 2 6Xa6yw to obtain a definition of the second difference compatible with discontinuous a ie let Wi1 Wi Wi W121 7 2 6xa6YWi 7 Beam aeffJi h I 7 hi hi1 hi1 and define the finite difference procedure accordingly Tihonov and Samarskii published a series of papers in the 6039s showing that this choice of the conductivity coefficient leads to improved accuracy the method then coincides with the piecewiselinear Galerkin method Tridiagonal Linear Algebraic Equations Tridiagonal systems of linear algebraic equations will arise in the discussion of a very large number of numerical methods to be treated in these lectures including the finite difference method Their solution is of course a simple exercise in Gaussian elimination Let ao cV 0 and write the linear system as ailV121 biWi CiWi1fi397 i 0 Forward elimination 1 Normalize the first equation 50 b6 17 V0 50 607 70 Bo fo 2quot For1N 5i bi ai Yi71717 39Yi ich WI 7 WWIAW Then WiViWi1 W7 077N17 Backward substitution 3 ForN71O Wi 4Pi YiWi1 Let us count the operations In the forward elimination there are slightly less than N 1 divisions l 4N l 1 multiplies l 2N l 1 additions In the back solution there are N multiplies l N additions Thus there are essentially nine arithmetic operations per unknown If the same tridiagonal matrix arises with several right hand sides then the factorization need not be repeated and only the right hand side operations are required in this case only three multiplies and two additions per unknown are needed to solve the equations Lecture 2 Finite Differences for a Mildly Nonlinear Problem Consider the mildly nonlinear two point boundary problem given by 7auxx l X7 u 0 u O XEI7 X68 Here a ax as before and fx has been absorbed in the the function X7 u Assume that The finite difference approximation of u is the solution of 76Xa6yw l cx7 W 07 177N717 W0WN 07 which generates a nonlinear algebraic problem It is necessary to establish existence and uniqueness separately We shall demonstrate uniqueness first and then determine existence as a consequence of the convergence of an iteration for its solution Lemma There exists at most one solution of the finite difference equation Proof Let W0 and W2 be two solutions and set 2 W17 W9 Then 76Xa6yz l cx7 W1 7 cx7 W2 07 i 17 7 N 71 Since by the mean value theorem cx7 W1 7 cx7 W2 oil2 it follows that 2 satisfies the linear equation 76Xa6yz l oilz 07 and since cf 2 0 then uniqueness is a consequence of previous results To establish the existence of a solution of the difference equation let us define and prove convergent an iterative procedure for it We should like to choose 3 such that B gt maXcuXC X 6 I7 C 6 IR however this would require us to assume that cLl is bounded on IR Let us derive a rough bound for W Assume for the moment that W is a solution of the difference equation and let v O for i 07 N Then v satisfies the equation 76X36yVi l cx7 vi cx7 07 i 17 N 71 So W W 7 v satisfies an equation of the form 76Xa6yWcW 7CX39707 i1N71 Now test the last equation against W 357W757WD 63W7 W CX7 0 W7 Thus 1 W1lt7C70 007 so that any solution of the difference equation is bounded in hi1 and by the Poincare inequality uniformly on I by p Hence it is possible to modify the definition of X7 v when lvl gt p to X7 v X7 signvp without affecting any eventual solution of either the differential or difference problem Then it is feasible to require 3 to satisfy 3 gt maXcuxC X 6 I7 lCl p Initialize the iteration arbitrarily but such that iwfmi p i1 N71 W30 WIS 0 7 7 Then with W k wk O for all k let 76Xa6ywk i wim wlik lh 7 cx7 W k71 zk W00 W0917 then 76X36yzki 32 lt5 7 of 2W1 i 1 N 71 Test this against 20 wasz lmizm 6 a c zltk1gt7zltkgt l k k k l k H H 2lt iz 26 72 7 so that a yz Mk g 6 c3552 2 l 7 k k4 k4 2 3 culz 7z Apply the inequality 367272 2 39yizig 2 39szHg along with the non negativity of cu 5 iizmiig miizw 1iig Thus if 7392 B i 2y then 0 lt 739 lt 1 and iizmiio S Tiizk 1ii0 S S Tkiziizmiio imply that k Wm Wm 20 a W j1 Moreover 76Xa6ywk i 3Wk 7 Wk71 i CX7 Wk71 76Xa6yWi CXh Wi7 so that W is the solution of the difference equation Essentially the same iterative argument produces a solution of the differential problem Note also that lkl1 a O in hl analogously the iteration converges in H6 in the differential problem Then it will follow from a lemma which appears early in the treatment of finite element methods to follow that u E H2I A better iteration can be based on the expansion 86 k 7 RA 7 RA 0 7 k4 cxW 7 cxW 8uXW W W gt 826 2 k 7 k4 0 lt auz W W gt This leads to the Newton iteration given by 5X357Wk CUWk1Wk CWk1 CUWk1Wk17 and quadratic convergence ie Hzkii0 Ciizk 1iig can be expected if the initial guess is sufficiently close to the solution W A safe strategy is to take a few of the iterations described earlier and then to shift to the Newton iteration The analysis of the error e u 7 W is trivial since e satisfies 6Xa679i l cjiei 8i7 I39lN717 808N 07 where lsil Ma7 ch2 for uniform subintervals Thus testing against e leads immediately to the error estimate llelll cmiamihz This finishes the treatment of the two point boundary problem by finite differences in these lectures The Piecewise Linear Finite Element Method We continue treating the usual two point boundary problem Lu auxxcu f XE01 u g XEB01 now by the finite element method that is the simplest analogue of the finite difference procedures considered above Assume for the time being that the coefficients aX and CX are smooth and that f 6 L2 Also assume that aX 2 oz gt 0 and CX Z 0 with norm k i z zlxzdx ilk gt and k seminorm lzllt llzkll0 As above let HampI z E H1I 20 21 0 The Poincare inequality in the lemma below is slightly more precise than that given earlier Lemma lfz E HampI then 1 llzllo S ilzll thusll H1 and l l1 are equivalent on HampI Proof The inequality is an immediate consequence of the string X l X 2 l 1 zXtdt gxz zXt dt glezll O X E 0 0 1 and the corresponding Inequality related to the Integral from to 1 lZXl lt suffices to treat vanishing boundary values If not let E u 7 g017 X l gNX7 so that 30 31 O and E f Lg01 X gNX f gN a goiax a Clg01 x W For smooth 3 and c E L2 if f 6 L2 so that E satisfies a boundary value problem of the same form as u and vanishes on 8 In what follows the assumptions on the coefficients will not be repeated unless a change in the assumptions is desired Lemma If f E L2I then the solution u of the two point boundary problem with g 0 satisfies u e H2I m H30 Proof This is a formal proof to convert it to a rigorous proof approximate f E L2I by a sequence of functions fm C C OI converging to f in L2I obtain the a priori estimates below for the solution um corresponding to n and then let m a 00 Test the differential equation against H and integrate the second order term by parts aux7 ux i cu7 u f u7 so that by the Cauchy Schwarz and Poincare inequalities aiuii CLAN S iifiioiiuiio S Ciifiioiuil Thus Huiil Ciifiio Now for smooth coefficients auxx f 7 cu 7 axux7 and iiUxxiio S C iifiio iiuii1 C iifiim as was to be shown Lecture 3 Next let us consider the approximation of a function u E H2I by continuous piecewise linear functions Let 0X0ltX1ltltXN17 I39X391X397 hiX397X391 Define the piecewiselinear interpoant Ju of u over the partition Xi to be Jux u1 i Wwi X17 X E I7 i 17 7 N i Lemma If u E H2I and Ju is its piecewiselinear interpoant over Xi then N uiJu C u2dxh1 cu hz H H0 L xx i i2 N cltZugxdxh2gt Ciuigh i1 miH miH iuiJuil i A where h max hi I Proof This lemma is a very special case of a general approximation lemma called the BrambleHilbert lemma Set Ju v and consider the approximation on the subinterval I By the mean value theorem there exists a point 5 E I such that LI 7 LI uX VXX7 X E I7 so that vX u1 l LIXEX7 X17 X E I Then the Taylor expansion with integral remainder gives 00 iil l uxi71X Xi71 l X t7X1uXXtdt VX Um71 7 uXEx7x1 t7X1uXXt dt g X vX7X7X1X uxxtdtX t7X1uXXtdt Then by the trivial inequality 3 l b2 2a2 l b2 and Ca uchy Schwarz X X H 7 v2x 2 t 7 XH2 dt ugh dt X71 xiii E 2X 7 X12E 7 X1 Mix df X71 2 EX7X13hX7X12 LI X df L Integrate this inequality over the left half lilof I 3 7 2 lt74 2 IZu v XdX7 64h IZ uxxdx A similar bound holds for the right half of I so that adding over the subintervals leads to h max hi 3 V 3 7 2lt7 h 2d lt7 2h llu Vll0764 IUxx X764lul2 7 Hence 2 N 3 Hui viio lt2 hm u dxgt Ciuighz i1 Next let us estimate iu 7 vil Differentiating the expansion for ux gives 5 u 7 vxx iX um dt x 7 x1uxxx so that uxivx2x 2 gagg u xdtX7X12u xx Integrate over I Iuxivx2 dxg2I hilugdthi2uixx dx 4hi2Iu x dx Adding over subintervals and taking the square root gives l 2 N u7v1 2 hIZu2 dx 2u2h Ixx H Now let us turn to the formulation and analysis ofa Galerkin finite element method based on the piecewise linear approximation space Mh v E C0I v is piecewise linear over X A basis for Mh is given by the hat functions 11 EMh i0N where 11 Xj6ij7 i7j077N Let Mhov Mh v0vN0 dith N l 1 and dithp N71 Galerkin finite element methods are associated with weak forms of the differential problems If our equation is multiplied by a function v and the second order term is integrated by parts and the homogeneous boundary condition applied the result is the equation Au7 v aux7 VX l cu7 v f7 V In order that the equation make sense u and v should belong to H1I For the boundary condition to be satisfied by u we need to have u E HampI Since the boundary condition is strongly imposed on u it suffices to require that v E HampI So the weak form of the two point boundary problem related to the bilinear form A7 is given by seeking u E HampI such that Auv 1 7v7 VV 6 HampI Other weak forms of the differential problem will be introduced later and associated with other types of finite element methods It is a standard topic in the study of differential equations ordinary or partial to determine the relationship between the solutions of weak forms of boundary problems and of their strong forms The hypotheses that we have imposed on the coefficients and the data function f are adequate to imply that the unique solution of the weak problem is such that u E H2I HampI and consequently satisfies the original strong problem The basic idea in the Galerkin method is to take a subspace of HampI such that functions in H U or in H2I H3I can be approximated closely We have seen that a function in H2 O H6 can be approximated by its linear interpolant which is in Mhp to 9h in H1 and 9h2 in L2 Thus it is reasonable to seek an approximation of the solution of the weak form of the boundary problem in Mm though we shall temporarily delay consideration of the rate of convergence of the approximate solution as the mesh spacing h a 0 Repeating the Galerkin finite element approximation Wh 6 M of the solution u is the solution of AW17 v awhx7 VX l CW17 v f7 v7 VV 6 Mm The problem is linear and finite dimensional LemmaThere eXists a unique solution of the Galerkin problem on Mhp Moreover the Galerkin solution is bounded 1 lWhl1 S llfllo 04 Proof Test against Wh alWhlg S aWhX7 Wxh CWIM Wh f7 W S llfllOlWhlh and all three conclusions follow Before proving convergence as h H 0 let us consider the algebraic problem generated by this Galerkin method The matricial equivalent of the inner product form of the Galerkin problem can be derived as follows The solution has the representation N71 Wh Z W UlX j1 this leads to the equations N71 N71 ameJMLX chwjwi fw i1N71 j1 11 Set O iJ a jX7 ixC j71Ji7 WI leiii Let A be the matrix with entries 04 and let W and 4p be vectors corresponding to W and 4 respectively Then the inner product form is equivalent to the linear algebraic equations AW 4p Since the supports of 15 and 1 do not intersect unless li ijl 1 the only nonzero elements in the ith row of A are 0431217 04313 afld aii1 Thus the tridiagonal matrix algorithm applies to the matricial formulation It will be convenient to consider either of the two equivalent formulations of the Galerkin method interchangeably in the discussion to follow Note first that both formulations are in terms of positivedefinite symmetric operators symmetry is obvious positive definiteness follows from the fact that alulg Au7 Comparison of finite difference and Galerkin methods Let us construct the linear equations explicitly in the case of uniform subintervals Note that the integrals in matricial formulation all are over intervals of length h or 2h and consequently contain a factor h and we need to compensate for this factor to give a closer comparison between the finite difference equations and the Galerkin equations Since h 1X 7 X171 X14 lt X lt X17 lox a x 11 XL X lt X lt XJ17 then 1 awiilx phx p c x ix 1 3 i1x7 ix pI axdx7 1 aipix lzlhx P UI 3XdX ai iemiiJLX 3 i1X739 lix7 and l ammo p I cxx ixxx 7mm 1 2 1 2 c1111 p I CXX7X1 dxp I CXX3917X dx7 1 1 cm17 p I cxx 7XiX391 7XdX 1 Let 1 3 Ea i71X7wix7 17m7N7 and 7 1 CiJEC LZJj7 LZJi7 jl71ll1l177N71 Consider the special case of the two point boundary problem given by quXu f xel u 07 X68 Then the Galerkin equation reduces to 1 pWi1 2W W121 l 2 l 1 5Wi1 gWi EWH EUMML I 177N17 whereas the finite difference equation reduces to using W to denote its solution Wi12Wi M471Wf7 1N71 Thus the second order term is the same in the two equations the zero order term in the finite difference equation is replaced by a Simpson quadrature and if the function f is interpreted to be its piecewiselinear interpoant the right hand side is also replaced by a Simpson quadrature The effect of shifting from finite differences to Galerkin finite elements is that the lower order terms in the equation get smeared a bit While this appears to be a minor perturbation we shall see that the solution of the Galerkin method maintains second order accuracy 9h2 for variable spacing in the grid Xi which is not in general true for the finite difference method It also requires less smoothness of the data in order to obtain this accuracy H1 and L2 Analyses of the Error in the Galerkin Approximation Since MILO C H U it follows that the error 2 u7 Wh E H6I satisfies the equations Au7 v 7 AW17 v Az7 v azx7 VX l cz7 v 07 V 6 MM We cannot assume that z 6 MM so that we cannot decide immediately that z 0 However we can write 2 in the form Z CCWh7 EAhp Then AWh7C7 v Au7C7 v au 7 OX VXcu 7 C v7 V 6 MM Since Wh 7 C 6 MM it can be employed as a test function in the equation above Q Wh CE S AWhC7WhC A C7Wh 0 MH CH1HWh CH17 where we have assumed that maxax7 cx M lt 00 Thus there exists a constant K such that HWhiC h S KHU CHh VC Mh0 So UC1 HWh7CH1 1 KHUCH17 VC Mh07 HLHWhHl M M Since the piecewiselinear interpolant Ju of u is in Mhp N llui Whlll S 1 KllUJUll1 S C Zlu i1 ciuizh where h max h and 1 dku 2 g kBltBltWgt dxgt BCI Note that the error estimate naturally localizes by involving estimates subinterval by subinterval l u The Hl inequality is the energy norm estimate of the error and frequently is the error estimate of greatest interest in engineering or scientific computations Note that the error is bounded in terms of the local smoothness of the solution and is optimal in the sense that it can be proved that if there eXists a sequence Mhmo of piecewiselinear spaces on I for which 1 hmaXhHO and inf ui 7H0 n I Wm H CH1 h 7 then u must itself be piecewise linear and the points at which the slope of u is discontinuous must belong to Xmi for n sufficiently large Let us look for a bound on llu 7 While This can be accomplished by means of a duality argument as follows Our current differential operator L is self adjoint so that its adjoint L is equal to L Let 4p 6 H2 O H6 be the solution of E0 W Xk l C4 L Wh Then ll pllg Cllu7 while So ll Whllg Uquot Wm WM Uquot W awxlx CW Au 7 Wh74p Au 7 Wh74p 7 v7 V 6 MM Thus ll Whllg S Mllui Whlll ll 7 Vll17 V E Mme Since 391 7 ltC hltC 7 h39 vewwlw viil lwlz iiu wile Hence N iiu7 While ciiu7 Whilth c Du i1 5 any h Clulghz Again this is an optimal order estimate for the error in L2 since llu 7 While 0h2 implies that u is piecewise linear We are unable to localize the h lifting of the H1 error estimate to an L2 estimate by this argument The analysis was given under the simplification that u O on 8 if instead u g7 X E 8 all that needs be done is to impose the boundary values on Wh and modify the Galerkin method to require finding Wh E Mh such that AWhv v7 VGAhp Wh g7 X68 Note two things First the test space remains Mm since the two additional degrees of freedom in Mh over MILO are specified by satisfying the boundary values Second a glance through the error analysis shows that the same error estimates will follow Let us collect the error bounds into a theorem Theorem Let u be the solution of the two point boundary problem and Wh that of the Galerkin approximation Assume that the coefficients a and c are smooth and that ax 2 04 gt 07 CX 2 07 maxax7 cx M lt 00 Then there exist constants C such that out N lluiwhlll S CZlu i1 llui Whllo S Cllui Whlllh 3117 S Clulzhi 1 2 gm h Clulghz l A Lecture 4 L analysis of the Galerkin error Continue studying the piecewise linear Galerkin approximation to the solution of the boundary value problem Lu7auxxcu f XEI7 u 07 X687 and continue assuming that a and c are smooth and that ax 2 04 gt 07 CX 2 Oand mlaxax7 cx M Recall the Green39s function CE7 X l X where I E H6I satisfies the equation L39X LU 5amp7 where E 6 01 and 6g is the Dirac measure at the point 5 More conventionally I can be defined as the solution of the system 7aI39XX i ci 07 X 6 05 U 517 7 717 0 It is easy to see that 39 6 H6 H207 H2 717 though it is not in H2I as a consequence of thejump in aI39X at E Let J H1I a M denote the operator interpolating functions in H1I into Mh for convenience we shall let JV denote the piecewiselinear interpolant of the grid function vi as well Let E X a node in our partition of I Then it follows from our approximation lemmas that now with GX7 miH hZ lll Jl llo hlll Jl ll1 S C ll 640 t ll lgx1gt where as before h max h Now note that for any v E H6I vx V6X v7 L39 Av7 I39 Thus if Wh is the solution of the Galerkin equations then Hi Whli A Whi r A Whi 39 i 7 C E Mao and i A u 7 W C U7 W inf I 7 h H hiilCEMWii CH1 S Ciui2h iii Jrih CiuiZirig0mirigm1gt h miH If in each partition XJn with h a 0 there exists a node n Xi we have proved that the error in the Galerkin X approximation is 9h 2 at that fixed point what is needed is a uniform bound on the term in the parenthesis in the inequality above Lemma There is a constant K such that mam mam K 0 lt 5 lt1 Proof Since 390 O and I is not identically zero aI39XO cannot vanish if it did then I would vanish on 05 and aI39X 71 Then there eXists an interval 57 7 gt E on which I lt 0 Then aI XXCl ltO7 E X 7I7 which implies that I39X decreases on 57 which in turn implies that I decreases on 57 so that it follows that 7 1 and that 391 lt 0 which is a contradiction A similar argument shows that aI39XO gt 0 Working back from X 1 shows that aI39X1 lt 0 Hence I increases from zero at X O to a maximum at X E and then decreases to zero at X 1 Thus it follows from the jump condition on aI39X that 39XgtO7 Olta39XXlt17 0ltXlt 7 I39xgt0 71lta39xxlt0 ltxlt1 Thus I39X and I are bounded on I independently of E 6 01 Then 1 I39XX 73axl39x 7 cl39 is bounded for all 5 6 01 and all X 6 05 U 51 which is stronger than we set out to prove This lemma leads to cun1HN7L Mn 7 Whi 7 where the constant C is independent of the partition With JWh indicating the piecewise linear interpolant of WM then Hwimmmgqm since the extrema of piecewiselinear functions occur at nodes So HuiJWh HuiJu HuiJu 000 iiJu 7 th 000 Ciui2h27 000 000 S S and we need to bound ui Ju in L I Let djv Koo 7 co 00 397 W I7 V6 L I TX 6 L 01717k with associated norm k k 39 de iiViikoo iVijoo TX 0 10 10 The restriction of these norms to a set B C I will be indicated by iiViikpoJB Lemma If u E W2gt I then HuiJuHOpO mlax uiJuHOpOJ C miaxu2OOJ h2 Qu gmhz Proof Let X E I Then g X UiJLIXX 7X7Xi uxxdt fX391LIXde7 1 X 1 x where 5 lies in I Thus 1 U JUXXN S X XFIHU ZWJWI39 E u ZpoJhiz S C u ZpoJhizv as was to be shown The immediate consequence of the above lemmas is the following theorem Theorem Let the solution u of the two point boundary problem be in W2gt I Then llu 7 JWhllopo lt Clulgmhz This completes the discussion of the piecewiselinear Galerkin method except for an observation concerning the very simple differential equation Lu iuxx f The Green39s function for this operator is piecewise linear with 39 G E 517 5 Obviously I39XX O for X 7 5 Thus for this equation Wh39LI397 Higher Order Galerkin Methods We shall treat Galerkin methods for the two point boundary problem based on piecewise polynomial subspaces of Hampl again associated with a partition I X1X7 i 17 7 N Let Mln v Ckl leEPnUi 1N n2k1 where is the set of polynomials of degree at mostj on the set B The most interesting cases are given by k O C0 piecewise 7 polynomial spaces7 k n 7 1 Hermite piecewise 7 polynomial spaces7 n 7 1 Spline spaces CO Piecewise Polynomial Galerkin Approximation Consider the space Mao v 6 M9 v0 v1 0 over the partition with nodes 0 X0 lt X1 lt lt XN 1 Let us define a local basis for Mao which we shall do first on the reference interval 711 Let 1 190 51X7 190 1X7 19 17x1xx1117x2x11 j1n71 Ml Note that 1940 and 190 represent the right and left halves of the piecewise linear basis functions on 711 respectively Lemma Pn711 Span190190191 19n11 Proof lt suffices to consider Xk k 07 n Now 1 1940 190 and x 190 7 1940 If k 2m n is even then Xquot 1X2m 71 1717X2X2 quot2 l sz 4 l 1 Similarly if k 2m 1 n is odd Xk X l XX2m 71 X 7 17 X2X2 quot1 l sz 3 l X7 which by induction establishes the lemma Bases for M9 and Mao can be described as follows 139 For i 07 7 N let who be the hat function at x W909 65 j0N 2 Set M Spanwm i0N Mao Span psi 1N71 3 For1N let wiJX19jltXX1ilgt7 jl7n71 4 Then MSMSPan U7 1N7 j17m7n17 BaSi5Mgo BaSiSMg1Joo7 o In order to estimate the error in the Galerkin approximation to be defined below it is necessary to understand the best approximation properties of the space M20 These can be summarized in the following lemma the proof of which will be given ater Lemma Let v E Hkl H U for some k such that 1g k g n i 1 Let the partition Xi be such that max h g h Then there exists a constant C independent of the partition Xi and v such that ian iv He hiiv 7 CH1 S Civith EAmo If v E Wkgt l HampI where 1 lt k n1 then there eXists a constant C again independent of the partition and v such that 39f 7mh7mltc oohk Ewgy iiv He iiv CH1 ivik The Galerkin approximation Wh E Mao to the solution of the self adjoint problem Lu7auxxcu f XEI7 07 u X687 is defined in like manner to that in the piecewise linear case find Wh 6 M20 such that AW17 v awhx7 VX l cw17 v f7 v7 V E Mao The H1 and L2 error analyses given in the piecewise linear case can be repeated very easily Note that again Au7 Wh7 v au 7 Whx7 VX l CLI 7 Wh7 v 07 V E Mao so that it will follow that lluiwhlllgC inf lluiglllgClulkhk l 1 k n17 6M0 no if u 6 Hr n Hg The duality argument can also be repeated in an analogous manner to that for the piecewise linear case so that llui While Clulkhk 1lt k n 1 Let us generalize the problem somewhat by adding first order derivative terms to the equation Lu7auxxib1uxb2uxcu f XEI7 u 07 X68 The bilinear form on H6 x H6 associated with it is A 7V auX7VX b1u7VX bQUX7V CU7V7 which is no longer necessarily symmetric or positive definite Assume however that the differential problem possesses a unique solution for f E L2I and that the regularity result llullz Cllfllo holds The Galerkin approximation Wh remains being found as the solution of AW17 v f7 v7 v 6 M20 Since it cannot be assumed that Av7 v 2 alvlg it is necessary to prove that the Galerkin problem has a solution or equivalently that uniqueness holds for its solutions So assume that z 6 M20 satisfies Az7 v 07 v E Mao Then 1 0 14272 2 will Clzllllzllo llleg Z allzllg Cllzllgi or lllel S Cllzllo Now invoke a duality argument Let 4p 6 H2 O H6 be the solution of PM aka wa 7 bwk cw Z By the regularity assumption above ll pllg Cllzllo and MS 27 PM AZ7 w AZ7 7 XL X E Milo Thus 2 z lt C 2 inf 7 lt Ch 2 z ll lloi ll1 6Mioll Xllli ll0 Thus llleo S Chllzlll S Chllzllo Thus for h sufficiently small llzllo 0 and uniqueness is valid and we have proved the following lemma Lemma There exists a unique solution to the Galerkin problem for the two point boundary problem including first order terms for h sufficiently small Lecture 5 A new argument along the lines of the uniqueness argument just presented must be given to bound the error this argument will develop the bounds in L2 and H1 simultaneously and begins with a duality argument Let 4p 6 H6 be the solution of L4p u 7 Wh Then ll pllg Cllui while and llui Whllg U 7 Why PM AU WW AluiwmwiCL EA1207 so that lluiwhllg S Clluiwhlll inf llwi lll S Clluiwhlll ChlluiwhllO EAlie Hence Hu 7 WhHO CHu7 Whth C u 7 Whhh Now the fundamental error equation Au 7 Wh7 v 07 v E Mao implies that AW 7 c v 7 Au 7 c v v e M2107 c 6 M29 So with v Wh 7C AWh C7Wh CA C7Wh O S ClluiClllllWhiClll S CluiclllWh Cll For v 6 H1 a simple calculation shows that AV7 V Z M CllVllglVl llVllS alVlF CllvllS Relations of this type one of which also occurred in the uniqueness argument above are known as Garding inequalities in the study of partial differential equations another of their typical uses will follow 1 galWrCl CluiClllWrCl llWrCllS CluiClllWrCl llui ll3llu7Whll3 ialWrCl CUluiCll lluiwhll sothat 1 ZalWh 7 cl C llu cl hzlu Wm Then lWrCl lui lg C llui ll hzlui W l Whl S S Thus for sufficiently small h lui Whll S ClluiClll Hence by the approximation lemma if uermHg 13kgn1 then lui Whll C inf lluiClll S Clulkhk l EMU mo We again have an optimal order energy estimate for the error in the Galerkin solution Moreover an optimal order estimate in L2 is an immediate consequence of the estimates above We can summarize these bounds in a theorem Theorem There exists h gt 0 such that for O lt h lt h there exists a unique Galerkin approximation Wh E Mao If u E Hkl H6I is the solution of the boundary value problem for some k such that 1 lt k g n 1 then there exists a constant C such that for O lt h lt h llui While hllui Whlll Clulkhk 1lt k n 1 The Matricial Problem for the CO Galerkin Method Recall that MSOM0Span19J i1Nj1n71 The support of 197 j 1n 71 overlaps exactly those of 1914 190 and 19A 1n 7 The matrix generated by the Galerkin problem is block tridiagonal and the resulting linear algebraic equations can be solved by an obvious block generalization of the Gaussian elimination scheme discussed earlier A more efficient solution technique which is known in the engineering literature on finite element methods as static condensation can be devised quite easily For i 17 N consider the set of linear equations corresponding to AWh719ixjf719ixj7 j1n71 These equations can be solved for the coefficients of 19 j 17 n 71 in terms of the coefficients of the hat functions at X1 and X plus right hand side quantities These can then be substituted into the equations generated by the hat functions Then the remaining equations for the coefficients of the hat functions is just an ordinary tridiagonal system Ie change the order of elimination in the Gaussian procedure The Hermite Piecewise Polynomial Spaces One important property of the CO piecewise polynomial spaces M9 0 n 2 1 is that the basis elements have support on either one or tIO subintervals ln constructing a code with one of these spaces it is easy to associate n basis elements with each node X1 namely the elements 19F and 19 j 1n 71 The Hermite piecewisepolynomial spaces have the property that the support of each basis function consists of two adjacent subintervals The degrees of freedom are associated with nodes Fewer degrees of freedom need be associated with each subinterval for the Hermite case than in the C0 case The standard Hermite spaces are given by HkM n n2k1 k12 The most commonly used Hermite spaces are the Hermite cubics H1 and the Hermite quintics H2 and most of the discussion below will be limited to these two examples Note that the local degrees of the Hermite spaces are the minimum possible in order to have Ck continuity at the nodes since k l 1 constraints are generated at each node by the the continuity requirements Thus there are 2k l 2 constraints imposed on the polynomial on each subinterval thereby requiring a polynomial of degree at least 2k l 1 on each subinterval Let us construct local bases for H1 and H2 Consider the cubic case first Let 7 01 be the reference element Then we need to represent a cubic function on T by the degrees of freedom consisting of its values at X O and 1 along with the values of its first derivative at the same points It is easy to see that the two basis cubics associated with the left hand end of are given by VAX 17 X21 2X and SAX X17 X27 the value and slope basis elements at X 0 respectively ie VAU 1 and Slf0 1 and the remaining three degrees of freedom for each vanish The basis elements associated with X 1 are given by VX X23 7 2X and X X2X 71 Now let us scale these basis functions to be applicable on I U I1 X1X1 where we do not assume that the subintervals have equal lengths The value basis function at X is straight forward V XELF WA X72 V X 6 1 hi1 The slope basis function 5 must be scaled so as to have derivative one from both sides of X thus XXi71 his XE 5139X I X 7 X hi151 lt h 7 X E Ii1 i1 Next consider the basis for H2 on the reference element Here we need three basis functions at each end of the reference element corresponding to the value the slope and the second derivative of the quintic polynomial at an end point The left value basis function can be assumed to have the form V1X17 X31 ax l bxz A short calculation shows that to have Vl 0 Vl 0 0 it is necessary that V1X17X31 3X12X2 and VX Vl17 X7 where VX is obtained by symmetry The slope basis function SI takes the form SAX X17 X31 ax Again a short calculation and symmetry show that SAX X17 X31 3X7 X 717 X Let Kg and Kr denote the second derivative or curvature basis functions Then 71 KAX 2 X217X37 KX Kl17x Now consider scaling these functions to the two subintervals abutting at X Then set X7X3971 Vr I 7 X 6 iv hi V2JX X 7 X V gt 7 X E li1 hi1 his X 3 x e I 52iX X L XI hi151 lt h 7 X E li1 i1 W X 7 X x e I 7 hi H2iX X X hIZJrlKl 7 X E I 1 Given these basis functions it is easy to write down the H1 and H2 interpoants of a function u E C2I ll NH quot 0 k JkuXx Z uivkX LIISkiX u kix7 k i0 Optimal order approximation properties for these interpoants will be proved later but will be used in the convergence analyses below The bases for 77 and Hg are obtained by deleting v and kaV from the bases for 775quot Note that the slope and second derivative basis functions at the endpoints of I are retained in the basis for H3 The Hermite Galerkin finite element method is defined by seeking a solution Wh 6 Hquot of AWh72 aWhX7 Zx b1Wh72x bZWhX7 Z CWhiz 1 727 z 6 H5 The analysis of the error in L2 and H1 was reduced to approximation estimates in the treatments of the Mao Galerkin methods both for n O and for n 2 1 and nothing changes in those arguments when repeated for the Hermite spaces except for the approximation results The following approximation result holds inf lluiCllo hluiCl1 Clulzh EH k 0 ifueH forany such thatk nl2k2 Consequently the following theorem is immediate Theorem Assume that the differential problem is HZ regular ie if f 6 L2 then u 6 H2 Then there exists h gt 0 such that the Hermite Galerkin problem has a unique solution for O lt h lt h Moreover there exists a constant C independent of the partition and u such that ll Whllo hl Whll S Clulzhli if0lthlthueH and k nl2k2 The L00 error estimate requires a different approach from that presented for the self adjoint problem as approximated by the piecewiselinear Galerkin process and will be delayed The linear algebraic equations generated by the Hermite Galerkin method are block tridiagonal since each basis function has support in two intervals adjacent to a node The blocks are 2 x 2 for H1 and 3 x 3 for H2 or more generally k 1 x k 1 for M20 n 2k 1 gt Some Extensions Let us address the question of discontinuous ax again For the concepts it suffice to consider ax to be piecewise constant and consistent with the partition ie let 3X3i7 xeli i1N Let us look at the Hermite cubic space 77 There are two possible remedies for requiring the continuity of the first derivative across a jump in the coefficient a The obvious one is to relaX the smoothness at a node across which a is discontinuous and just to require continuity of the solution at these nodes It is trivial to modify the definition of the interpolant in the intervals abutting this relaxed constraint with the result that the convergence theorem can be restated with the error bounds being treated locally subinterval by subinterval as obtained in the approximation arguments earlier for the Hermite case with smooth ax However a more elegant and computationally simpler remedy is to modify the space 77 by replacing the slope basis function 5X by lhisr 7 X 61h a h X 7 X ihi1sllt 7 X E Ii1 3i1 hi1 Making this modification in the basis is simply imposing the correct flux consistency condition on the finite element space There is no difficulty in extending the convergence theorem to this modified Hermite cubic Galerkin approximation again the approximation results need only to be applied locally Galerkin Approximation by Cubic Splines Given the usual partition 0X0ltX1lt ltXN1 Let us define the Sn spline space over Xi to be M l for n 2 3 Call 3 3 M the space of cubic splines We shall confine our treatment of splines to the cubic case Let us construct a basis for 33 Note that the constraint that v E 3 3 lie in C2I requires three constraints on v at each endpoint of a subinterval I thus it is not possible to give these values for a piecewisecubic function at a point X and then have enough free parameters in the two adjacent subintervals I and I1 to localize a basis function on I U Ii1 Consequently the support of a basis function to be associated with X must spread over XijXj wherej is an integer to be determined by the Cz constraint The number of parameters needed to define the cubics in a basis function bX needs to equal the number of constraints related to the nodes between Xiij and Xj gt Total parameters 4 per intervalx 2 intervals 8139 gt Constraints i At X Xiij biXiij blXiij biXiij 07 ii Atxi Z1111j71 bllklxf bfker k 012 12 11 1 71 iii At xi bib biX 17 14 b ker k 12 iv Total constraints 6j 4 Thuswe need 8 6 i 4 orj 2 so that a basis function to be associated with X should be constructible on the set X2X2 The basis function bX on lisl X2X1 must take the form biX 041X Xi7237 X E 1217 since blikX2 O for k 012 Then on I biX 04lX X1223 BAX Xi7137 X E Ii since this combination agrees through the second derivative at X1 Similarly biX 04rXi1 X37 X 61427 b X OzX2 7 X3 BX 1 7 X37 X E hurl There are four parameters remaining in the expressions above and there remain the four conditions at X X A short calculation shows that the constraints lead to the following four equations IA721 hi3 zh 1 04rhi1 hi23 rm31 1 Old721 hi2 5N7 04rhi1 hi22 BrhiZJrl 07 04lhi71 hi Blhi 04rhi1 hi2 rhi1 0 For the special case of a uniform partition by symmetry bX b 7X and bX 0 so that the system above reduces to the pair of equations 800 l Br hisv 40 l Br 07 which has the solution Thus bX is given by 1 lxixl 3 lxixl 3 727 39gt 7lt1ii39gt 0ltlX7Xillth7 b 7 4 h h I00 1 3 Zlt2ilxihxllgt hltlxixillt2h Thus the basis function b is a smooth bell shaped function and it should be possible to interpolate a smooth function in terms of cubic splines Next consider a different special case In many if not most practical problems described by differential equations the solution has regions in which it varies slowly and regions in which it varies rapidly and can have rapidly varying derivatives For these problems it is natural to employ a variable grid spacing and it is of interest to discover the behavior of b where the spacing changes Let hi4 hi 5 lt17 hi1 hi2 1 Then the equations for the four parameters reduce to 804 Bl 5 37 800 Br 17 462W 62m 4a B 07 26041 6B i 200 Br So 726 7 126 047712637 7 7363 7 7126 7 26 6 126 3639 Consequently as 6 a 0 b X391 X 1126 whereas for a uniform partition ie for 6 1 bX1 25 As a result it is most likely that an attempt to represent the approximate solution of a boundary value problem in the neighborhood of a strong change in spacing in the partition will lead to oscillatory behavior which is of little practical use Thus we are motivated to restrict our consideration to the case of a uniform partition and relatively smooth solutions of the underlying boundary value problem Lecture 6 The transparencies for Lectures 1 5 are available at wwwmathpurdueedudouglas Click on math598C and then on Lectures1 5pdf It is also difficult to take a discontinuous ax into account Thus it seems not too likely that a cubic spline Galerkin method is the ideal way to approximate the solution of a realistic physical problem having a solution with significant character However for the sake of historical completeness we shall analyze the procedure Note that there are many applications of splines in one or more dimensions that are very appropriate though most do not involve differential equations Before going to the differential equations problem let us consider interpolation by cubic splines Assume a uniform grid X ih Note that the basis functions corresponding to X71 and XN1 two nodes external to I are nontrivial in I There eXist a large set of suggested remedies to account for having these basis elements in the interpolation process most of these suggestions involve making some assumption as to the behavior at the end points of I of the function being interpolated Seek an interpolant that is given a function fx on I let us look for a cubic spline N1 Bx Z emuX i71 such that Bxi f i O N 7 Clearly these relations do not uniquely determine B since we have only N l 1 equations in N l 3 unknowns If f is reasonably smooth then we can augment the equations by requiring that the slopes be matched at the end points then the interpolant would be determined by the equations Bf7 I39ON7 BX fX7 X 01 Order these equations by taking the B O f 0 equation first then the B f equations in order and finally the B 1 f 1 equation last It is easy to see that this system of linear equations is fivediagonal and is nonsingular A simple extension of the tridiagonal Gauss solver handles these equations efficiently It is however nontrivial to show that this interpolant has optimal order approximation and a different type of approximant will be constructed to demonstrate the optimal order approximation properties associated with cubic splines Let us develop a basis for 3 8 v E B3 v0 v1 0 Now the basis elements b1b0b1 can be replaced by two entries b3 and bi that vanish at X 0 b3x boX 2b71X 17100 1 X 3 3 x x 3 l X 3 32 Zgt3 lt1 Egt2ltEgt lt1Zgt OltXlth I lt3iggt hltxlt2h Similarly replace bN1bNbN1 by bivi1 and where ivalx bN71X bN1X7 b7iX lex 2bN71XbN1X Then a basis for the N 1 dimensional space 3 8 is given by 3 8 Spanb87 bf b2 7 qu7 by17 biv7 and the cubic spline Galerkin approximation can be defined as the solution Wh 6 3 8 of the equations AW17 v aWthX l b1Wh7VX bgwhx7 v l cw17 v f7 v7 v E 38 The analysis of the error ui Wh can as usual be reduced to approximation and the bounds to be derived later can be applied to show that ll Whllo hllui Whlll S Clul4h47 when h is sufficiently small as to assure the existence of a solution of the Galerkin equations input part3 Elliptic and Parabolic Finite Difference Methods A Model Elliptic Problem Consider the Dirichlet problem iAu f xeo012 u 07 X E 89 Let 1 h N Xi 177 0 J h Xij thj39 quot039 thj39 Denote the discrete Laplacian by Ah where 1 AhZL j 5yx5yyzij p lieu zi1J ZAH ziJ1 42039 We shall consider the simplest finite difference approximation of the model problem iAhWij 1 X0quot 6 Qh X0 6 Q W0quot 07 X0quot 6 89h X0 6 Convergence analysis 2 estimates As in the treatment of the two point boundary value problem we shall not attempt to give the sharpest possible error estimates for the finite difference approximation Now Ahuij fij 807 X E 917 where 6039 0 UL 49072 Then if e0 My 7 W0 the error equations are expressed by iAheU 8039 X0 6 Qh e0 07 X0 6 89h Set Vh ivh 5757 5W 5757570 N71 N71 N N 2 67555777th Z Z 57506777039h239 i1 11 i1 j1 The Poincare inequality holds again N71 Mile 5 SClthEllm sens iJ1 where hi1 E E 2 EU 07 xij 6 89h In a single space variable lthEllo dominates llEllom however in more than one variable it does not Now test the error equation against e Aheie Vheivhe 876 S llEllollello By the Poincare inequality lthello S We S Clul4ooh27 and llello S Clul4ooh2 These bounds can be collected in the following theorem Theorem Assume that lu 490 lt 00 where u is the solution of model Dirichlet problem Thenif W is the solution of the finite difference approximation llui Wllo WW 7 Mile S Cl 40th Lecture 7 Convergence analysis an 600 estimate The 6 analysis is based on the following maximum principle Lemma Let 2 be defined on Qh U 89h and assume that Ahzij gt O on Qh Then maX Zquot maxzij mam am Proof If 2 were to have a maximum at a point xij E Qh then Ahzij 0 Now consider the auxiliary function Z i0 W for which Ahzij 1 Note that iAheiji S Ciuizumh2 1 Let C 39y i 72 i e where 7 gt 0 then AhC vnAheZ vnvgt 07 so that 1 23th 38 7 7 r332 EWJF n Thus 1 GUS vnzueij 5vn7 Vngt 0 Analogously 79 y n2 and 1 Heiiom 5v ciuiimhz Theorem If the solution u of the model Dirichlet problem is such that lu 4 00 lt 00 and W is the solution of its finite difference approximation then llui W 4mh2 000 Cl Actually we proved a bit more Replace 39y by v ng xlAheo39ly which for less regular u can fail to be 9h2 Then we have shown that llu 7 W 000 S Y The Finite Difference Algebraic Equations Order a vector corresponding to WU as follows 7 tr W 7 W117 W127 7 W1N717 W217 7 W2N717 7 WN71N71 Then the difference equations for the model problem can be written in the form Awf where each row in A has at most five nonzero entries however two of the diagonals are removed by N 71 places from the principal diagonal with the other two being adjacent to the principal diagonal Thus if the equations were to be solved by standard Gaussian elimination the factorization of A would require 9N2 arithmetic operations for nearly all of the N7 12 equations The back solution requires 9N operations per equation Consequently the total work required for Gaussian elimination would be 9N4 operations This leads to investigating iterative methods to see if cheaper techniques can be devised 1 There are many possible iterative methods 2 Each method needs to be effective for the algebraic equations associated with more general problems than the simple equations given by the model problem 3 It can be useful to consider the model problem since both the formulation and the analysis of an iterative method are easier to understand for this problem 4 More general and more efficient iterative methods for treating algebraic equations of this nature will be discussed later El Jacobi Iteration By considering the unknowns in the order specified above the ijth equation can be written in the form Wixjil Wiiu 4W0quot Wi1J WiJ1 h2 j7 and the Jacobi iteration can be defined by the following recursion gt Let w0 be chosen arbitrarily gt For n 2 0 let 1 1 Win 1ltW n39 u i1J u1 n n n 1 U W wig Wltgt Wltgt gtZh2fif Let us prove that this iteration is convergent though not rapidly Let Cn Wn1 Wn indicate the difference between successive iterates if 2 llCWllo lt 007 n0 then convergence of the iteration is assured Now n 1 n n n n n n 3 H Z Ch11 LL cltld Cid11 c3 l jAth x0 6 m If A were to be replaced by the differential Laplacian we would be led to using an eigenfunction expansion with eigenfunctions sin 7rpxsin 7rqy7 p7 q 127 Fortunately the same technique works in the finite difference case as a consequence of the following lemma Lemma The eigenfunctions of the discrete Laplace operator Ah subject to vanishing boundary values on 012 are given by sinIrpxsinn39qy7 p7 q 1N71 moreover 4 h h iAh sin 7rpxsin 7rqy p sin2 sin2 sin 7rpxsin 7rqy7 sin 7rplxsin 7rqu7 sin 7rng7 sin wqu H6p1p26q1q27 where H gt O is independent h Proof The equation 6Xsin 7rpx 7 sin2 WTph sin 7rpx results from standard trigonometric identities and the evaluation of the eigenvalue follows The orthogonality claimed in the lemma is implied by the inequality of the eigenvalues for p1 7 p2 or ql 7 qg That H is independent of h is another consequence of standard trigonometric inequalities Set Appq sin 7rpxsin 7rqy and expand W in terms of Appq pq1N71 N71 Cn Z ath quxU pq1 Then 04 17 sin2 WTph i sin2 04 ppqagg Note that 71ltpN71JV71ltltppqltltp11lt1 n so that 3qu a O for all p and q thus the Jacobi iteration converges Let us estimate the rate at which the iteration converges Now 712 27rh maxippqi pm pN71JV71 17 2sn 17 Let us determine how many iterations m are required to assure that m 04 1 max 7 pq aw 2 Pa By taking logarithms and expanding the logarithm of 17 713922N2 as N a 00 we see that it is necessary to have to assure that desired reduction is satisfied The number of iterations required to reduce the error by a factor 7 is 9N2 logn 1 which implies that the total number of arithmetic operations required for this reduction is 9N4 logn 1 which is not competitive with some other also simple but more sophisticated iterative methods However this very simple iteration can be modified slightly to be quite effective in reducing the error in the span of the more highly oscillatory eigenfunctions As a result the modified Jacobi iteration is very commonly used in the most efficient iterative method known for the algebraic equations associated with numerical methods for approximating elliptic equations this method known as the multigrid method will be discussed later Let us consider the modification of the Jacobi iteration known as Richardson extrapolation The recursion can be changed to 2 n1 n h n W 7 WU l VT AhWIj 7 f where the parameter 39y gt O is introduced to accelerate the convergence rate for the portion of the error associated with one set of eigenvalues while possibly slowing the rate for another portion of the error We shall concentrate on the very oscillatory modes given by 0 wpq i maXP7q Z N27 which includes almost exactly threefourths of the eigenfunctions Then 2 7rph 2 7rqh 1 27139 sin Egsin Tsin T 27 AppqEO lf0lty 2 then h h pqu 17yltsin2sin2gt E 17371 C 7171 qu O7 while 1 9pm 6 1 2 1 5v wpq 6 0 Since lppq yl remains uniformly less than or equal to one for Appq 9 for O lt 39y lt 2 we can choose 39y so as to minimize maqueo lppq yl as a function of 39y This can be accomplished by finding V such that 1V12 or PE for which 4 3 maX i We pp 5 5 Thus for the more oscillatory three quarters of the eigenfunctions this optimal choice of the parameter 39y causes the reduction of the coefficients of these eigenfunctions in the expansion of the error at least by a factor of three fifths in one iteration instead of taking 9N2 iterations to effect that reduction It will be this characteristic that will allow the use of the Richardson iteration to produce a very effective multigrid iteration A Model Parabolic Problem The model parabolic problem is given by the initial boundary value problem for the heat equation 7Au fxt7 XEQ7 t 0TJ7 u 07 x E 89 t 6 J7 ux0 gx7 x E Q Let us apply the techniques discussed above for the Dirichlet problem to the heat equation First denote the value of a function fx7 t at the point xij7 t where At gt O and t nAt by f0 The Forward or Explicit Difference Method The simplest finite difference procedure for the heat equation is given by approximating the differential equation at xij7 t by a forward difference in time for the time derivative and by Ah for the Laplacian thus the approximate solution Wig must satisfy the system for 0 g t lt T m 7 W n n TiAhWIj 1 X0quot 69h wgll 07 X0 6 89h W3 gij x0 6 Qh Since we can solve explicitly for Wig1 as W1531 Wigl l AtAhWIJr7 x0 6 Qh this procedure is called the explicit finite difference analogue of the heat equation clearly the evaluation of the solution at the new time level is a trivial calculation It is easy to see that u agrees with W initially and on the boundary discretely and that un1 7 un At where the subscript U will be suppressed where not needed for clarification and iAhLlnfn n7 XEQh7 8339 0lluttllooAt llLl 49072 Thus if e u 7 W denotes the error in the approximate solution then en1 7 6In 77A e 6 XEQ At h 7 h We can employ the same eigenfunction analysis procedure as in the treatment of the Jacobi iteration to obtain an Abound on llenllo Expand e and 6 in terms of Appq N71 N71 n 7 n n 7 n e i Z 0 ququ 5 7 Z 6pq qu pgtq1 Pq1 Then 1qu lt 7 4 sinszph i sin2 a2qAt6gq Now 4At 27rph 27rqh At 1 7 7 7 7 lt 7 7 nlgx 1 hz sin 2 sin 2 71 Ifand only If hz 7 4 The time step restriction is called the Courant Friedrichs Lewy CFL stability condition its consequence is that each error mode evolves boundedly the condition can be replaced by At h 2 g 1 CAt4 where C is independent of I77 h and At It will be used to establish the convergence of the approximate solution to the exact solution Set 4At 239Ii39ph 239Ii39qh ppq17ltsin Tsin 7 Let the CFL condition be satisfied Since e0 0 agq O and then n71 an Zqu 163q H0 By the orthogonality of the eigenfunctions N71 n71 2 quot1 2 ieni Z Zippqt qur Zazqm pa 0 pq1 n0 ism in pq n0 n0 n0 pq n71 n71 t Hg iigAt g t C1u2At2 C2u2h4 At H0 H0 S tn2 C1U2Af2 C2200 i A Thus l A lle llo tquotC1LIAf C2Uh2 l A Ct At h2 gg unl g lUWllumyl gt Theorem Let Wh be the solution of the explicit finite difference 1 equation and assume that the CFL limit At ihz has been imposed Let u be the solution of the heat equation and assume that um um and um are bounded on Q x 07 T Then the error e u 7 Wh satisfies the bound above It is also quite easy to prove an ZW convergence theorem Rewrite the error equation in the form n n n n n n n 1 17 Y 4 me endwish sisii e J11At637 and note that all five coefficients of the eff terms are nonnegative if the CFL condition is satisfied moreover they sum to one Thus 1 is i mgxieih Armkgxism so that n iiGWDiiom S iiequotii0oo iilt quotii0ooAt E Z ii miioooAf7 m0 and under the hypotheses of the last theorem MMquot 7 WWW om hz The 6 argument can be extended to the parabolic equation with variable coefficients with little difficulty Consider the initial boundary problem cgivlaVu f XEQOlt1 ltT7 81 u 07 X689 Olttlt T7 u L107 XES1 O7 where c CX 2 6min gt O and a ax 2 amin gt 0 The corresponding forward difference equation is given by 1 Wig39n 7 m Ci where V thW 6ya6xw l 6ya6yw Let Then the error equation can be written in the form n1 At n gt 1 icvh ah N n n n n lt3iJei1 39 aii xjeiilj aij eij1 aiji eijilgt v J and it follows that the error satisfies the same bound as for constant coefficients provided that the CFL condition 1 lt 3039 7 mln h 7 53 hods Remarks 1 If 5a is adjusted to be the sum of the off diagonal values of the coefficient a then the resultjust obtained is valid in n dimensions for n 12 2 The CFL condition imposes a severe constraint on the time step in the forward difference procedure it arises from the fact that the solution ux7 t At of the heat equation depends on all of the values of u at time level 1 If 739 gt O is fixed then the CFL condition requires that the approximate solution Wx7 t 739 to depend on the values WC7 t for K 7 xi 4Th as h a 0 A Disappointing Attempt at Higher Order Accuracy Since the time step allowable for the forward difference equation must be 9h2 it would be of significant interest to obtain higher order accuracy with respect to the time discretization A first try would be to approximate the time derivative to second order Restrict attention to the heat equation in a single space variable and apply a centered first difference to the time derivative Wn1 Woel 2At 6yXWn Consider the evolution of the eigenfunction sin 7rpx Wn 7 n 7 p sm 7rpx7 where Let 4A up 7t sin27rph2 Then p2 Mpp 7 1 0 lf p1 and p2 are the roots of this quadratic equation their product is 71 and neither has absolute value one If p1 gt 0 then 1 1 12 E P2 lt Mpltl plgt 0lt7r2Atwlu1ltIu2ltltluN1li7r2At Now 1 So when up 2 i 1 pup 71 V lt 7125 and this explicit finite difference equation is unconditionally unstable Lecture 8 The Backward Difference Equation The object here is to eliminate the CFL condition thus we need wx7 t l 739 to depend on values of WC7 t for C well away from x One way is to replace the forward difference in time by a backward in time difference as the approximation of the time derivative Thus we are led to considering the backward or backward Euler difference equation Wn 7 Wnil At 7 AhWn f X E Qh w 07 x 6 89h w0 g7 x E Qh It is easy to see that Mr 7 unil n n n T 7 Ahu f l 6 7 63 0 lLIttlAt lUxml lUyyyyllhzl If e u 7 W and 6 are expanded in terms of the same eigenfunctions Appq sin 7rpx sin 7rqy so that N71 N71 n 7 n n 7 n e Z O quPqi 6 Z 6pq Pqi pgtq1 Pq1 then 1 n n71 6n At 1Apuqmtltam W gt7 where 4 h p E sin2 Now 1 0 1 ltqu 1pqAtlt for any At gt 0 thus stability holds without the imposition of a CFL limitation The practicality of this method will be discussed after we show that its solution is convergent An analogous convergence argument as that given for the forward difference procedure can be given though without a CFL constraint Again e0 0 Thus n 2 n He ii3 ZltZppq 16qmgt tquotZZ6gq2At pq n1 n1 pa 3 t H6 H3At S t 2 C1LI2At2 C2002 n1 and as before n lt n 2 39 iie Ho Cr gjvjiuuiAH gg iuWi iuymiih Theorem The same error bound holds for the solution of the backward difference equation under the same regularity assumptions on the solution of the heat equation as for the forward difference equation without the imposition of a CFL condition Consider the evaluation of the solution of the backward difference equation At each time step we must solve the equations 1 1 1 EWyiAhWUBU JUEW5 xijEQh WU 07 X0quot 6 89h an algebraic system closely resembling that generated by the finite difference approximation of the Dirichlet problem Clearly this is likely to cost more arithmetic operations than the evaluation of the solution of the forward difference equation lterative methods need to be applied to it and there are two features that distinguish the difficulty of the algebraic solves here from those coming from the elliptic problem First in an independent elliptic problem there is no easily identifiable good initial guess in general here however we can find initial guesses that are very good For example if W0 denotes the initial guess for W W then taking n71 0 7 W0 7 W0 is very reasonable In fact un 7 Wnil unil 7 Wn71un 7 unil 0At 7 72 ie Wrquot1 is already to the same order in At and h as accurate an approximation of u as W thus in some sense the iteration serves only to preserve the stability of the backward difference procedure For n 2 2 an even better initial guess can be obtained by extrapolating the approximate solutions for the last two time levels 07 n71 n72 W0 72W 7W0 In either case it can be shown that a fixed number of iterations of a simply preconditioned conjugate gradient procedure this iterative method will be discussed later suffice to maintain the asymptotic rate of convergence of the approximate solution 9At l hz to the solution of the differential problem Second the time difference term induces a small diagonal dominance in the equations which speeds up any of the iterative procedures One Hint at a More Efficient Difference Method What the backward difference equation indicates is that moving the spatial part of the difference operator to the advanced time level eliminates the CFL constraint What we saw in studying the two point boundary problem is that the banded matrices generated in their approximate solutions led to easily solved systems Consider splitting the application of the discrete Laplacian into its X and y parts Wn1 7 Wn T 7 ywn 7 ywn fquot x 6 oh W 1 07 xea h The algebraic equations can be ordered so as to give N 71 sets of tridiagonal equations in the X direction for each fixed yj so that the evaluation of their solution is essentially no more costly than that for the forward difference equation The corresponding error equation is of the form 6ln1 7 eN At 7 gas 1 7 f yen c with 6 9At l hz If the separation of variables argument is repeated then 1 ApAtagq1 17 Athagq gq and 7 likth ppquot 71 ApAt39 If p is small and q large it is clear that a CFL condition must be invoked to ensure stability the admissible time step is twice that for the forward difference equation Now the hint is to think about what happens if at the next time step the roles of X and y are interchanged This will be taken up after a second hint as to better accuracy is presented Better Accuracy the Crank Nicolson Method All of the difference methods above have led to first order accuracy in the time step All have replaced the differential equation at the point xij7 t by onesided differences in the time direction while using centered differences in the spatial directions Consider the Crank Nicolson difference method which centers the difference expansions at the point xij7 t n1 7 n 7 AMwnll l W f x E Qh W 1 07 xea h lts error equation is given by en1 7 6In l n1 n 7 n At 72Ahe e76 Since 6 0UmAf2 0AuttAf2 0Uxxxx Uyyyyh27 U second order accuracy in the time step has been achieved locally The usual separation of variables argument leads to the relation p 1 WP mm Pq 7 1 axp qAt so that 71lt ppq lt1 VAt gt 0 Consequently the Crank Nicolson equation is stable without a CFL constraint Another repetition of the convergence analysis can be made to establish the global error estimate llenllo S CUf At2 hz Additional regularity of the solution of the differential problem must be required in comparison with that required for the explicit and backward difference equations if the optimal order convergence rate is to be obtained Convergence can be established at reduced rates with reduced regularity The algebraic problem generated each time step is essentially equivalent to that occurring for the backward difference equation and requires an iterative solution Extension to Variable Coefficients Let us again consider the initial boundary problem 6amp7V3VU f XEQOlt1 ltT7 8139 u 07 X689 Olttlt T7 u uo7 XES1 O7 where c CX E cmin7 cmax C 000 and a ax 2 amin gt 0 Its Crank Nicolson approximation is given by the difference equation mmquot 7 vz ltth WW1 We n3 where Wn1 7 WW n 6 At If e again represents the error u 7 W it follows that n 1 n n n l c6te 7 EVE th lte 1 l e l 8 2 Since em 0 Lea v ltth 90 8 Test this relation against them and sum the second term by parts 1 theUWhe n gltgt6telt0gt c teml7 6te0 l m By Cauchy Schwarz catemk me iiavheakvhem clewii Now test the error equation against 6te and again sum the second term by parts 1 c6te 6te mavhen17vhenl l im avhewhvhem sltquotgt76teltquotgt ic temtew gciidn ii Combine the first terms on each side of this equation multiply by 2At and add for n0to nm71 mil mil Dene 6teltquotgtAttheltmgtView C E Hawt hi m n0 n0 where the initial inequality has been applied to bound the07 Vhe0 Under the assumption that n1 50 2 O At2 hz it follows that 2 TAt mnax thn 7 WM l0 lt Z ll6tu 7 Wquotgt5At C At2 hz n0 By the Poincare inequality this implies the additional error bound mraxHu7 W 0 g C At2 h The 2 error bound can be derived directly by testing the error equation against e 1 Then 5ten7 en1 avhen17 Vhen1 6quot7 en1 Since 1 n n1 7 n1 n1 7 n n c6te 7e 22Atce 7e ce 7e i thequot17Vhequot1 Z aminithequot1ii37 sltquotgteltquot1gt celtquot1gteltquot1gt 2LHEW9H37 Cmin Mi then 17 Atce 1e 1ice equot 2aminiivhen1iigAt 1 iiequotgtiigm min Multiplying by 2At and summing on n from one to m 71 leads to the inequality ceW em 2amin Z livhequotiigm n1 g cew eltngtm i Z Hawaiigm n1 C min quot1 The following lemma which is a variant of the familiar Gronwall lemma that is frequently used in bounding solutions of ordinary differential equations will permit us to bound the error em Lemma Let f 2 Oand g 2 Ofor n 1 and 0lt MAtlt1 If n71 f g MkaAt k1 then 1 n71 k nlt n 39 f g 17MAtgg At Proof Note that f1 g1 f2 g g2 i f1MAt and f3 g g3 i g1MAt i MAG i g2 MAt An induction argument shows that n71 f g g ng ltMA1 MA1 Rgt7 k1 from which the conclusion follows Applying the Gronwall lemma implies that 1 n 2 iieHiisAr l 1 Cmin 7 n71 R71 1 l m E E lie At quot k1 1 ce en Thus for sufficiently smooth solutions u the 2 error bound has again been demonstrated In the finite element analyses to follow the regularity required of the solution of the differential equation will be traced more precisely and it will be seen that less regularity could have been required of u in the 2 error derivation given immediately above in comparison to the minimum necessary in the h1 derivation that preceded it it is feasible to carry out these more precise arguments in the finite difference case but it will not be done here An Alternating Direction Difference Procedure The object here is to produce a finite difference method that is both second order correct in time and space and is comparable with the forward difference equation in ease of evaluation This will be our first example of an operator splitting method We shall define the difference procedure for a full time step through a two sweep process first an X sweep followed by a y sweep Let 1 Wn Wn 1 1At 7 fawnli 7 67Wn Fl 7 x E Qh 5 1 W 1 7 Wn 2 n 1 7 2 n1 7 n lAt 73wa 2 76WW 7 F27 x E Qh 2 1 Wn W 1 07 x 6 89h W0 g7 x E Qh where F1 and F2 will be determined momentarily Each half step leads to a family of two point boundary problems indexed by the y coordinate for the X sweep and by X for the y sweep hence the algebraic problem consists of solving tridiagonal matrix problems which is easy and efficient 1 Difference the first two equations and solve for W n 1 n n At n n At n n W2 W1W T6VWI WTF1 F2 1 Then eliminate W from either of the first two equations to obtain the equation n1 7 n AIKWn1 Wquot 6 Y67Wn1 Wquot 1 n n n n F1 F2 T6 YF1 F2 The left hand side is centered on xij7 tnl and is a second order correct in time and space approximation of the heat operator So the right hand side should be chosen so as to be equal to fxj7 tnl Clearly the simplest way to do this is to set am F2xu my M W 2 Then n1 7 n 1 At W At W 7 EAhW 1 l W l Tz yz wnll 7 W f and the perturbing term is OAt2 7 The error equation for this alternating direction method e 1 7 e 1 At T EAME HH 9quot T en en EH7 where 63 ammo 0uW uymw 0UWAt2 Again we can apply the eigenfunction expansion technique Then it follows that 1 n1 7 n n 7 6 apquot ppqapquot 1 pAt1 th pquot ppq 7 1 th 1 pAt 7 1 pAt 1 th 1 ApAt 1 th 7 1 ApAt 1 th39 Each of the fractions lies in 711 for any At gt 0 so that the alternating direction method is stable without a CFL restriction Yet another repetition of the convergence analysis shows that lle llo t Cu my 12 Lecture 9 This particular alternating direction method would seem to have an extension to three space variables by taking three time steps of length At3 and moving one part of the Laplacian ahead on each but the resulting equation is stable only under a CFL restriction However there are other ways to obtain both algebraic simplicity and unconditional stability as we shall see neXt A General Approach to Alternating Direction Methods Let us consider the parabolic initial boundary value problem 8H 5 Au 7 f x E Q u 07 x E 89 70 g x 6 97 where d A ZAi i1 with 82 8 A iaxi i bxi i cx 8X2 8X being a differential operator in the single variable X The leading terms are often found in divergence form so that 8 8H Aiu 787m gm87m i biugt i cu The development below is applicable to either divergence or nondivergence form We start from an underlying satisfactory discretization of the differential problem such as the Crank Nicolson method Let z 26 be the solution of zn1 7 zn 1 At t 5M2 z f where d Ah ZAhm H1 with Aim being an appropriate discrete approximation to An The general alternating direction procedure is given by the d sweep procedure where for n 1d7 Wn1n 7 Wn H d At zAhJWn1J Wn Z AhJWn fn 7 j1 jn1 where d 2 AhJW E O and Wnll d Wnll jd1 Computationally the most convenient way to carry out this algorithm is to employ the equation above for n 1 and then to shift to the equation n1n 7 n1 71 1 which is obtained by subtracting the n 71 equation from the nth equation The equation for a 27 d can also be used to eliminate the intermediary solutions theoretically starting from a d and working down to obtain the perturbation of the Crank Nicolson equation satisfied by the alternating direction solution For d 2 and 3 we obtain the equations W 1 7 W 1 At T l AWn1 W TAh1Ah2Wn1 7 Wquot r d 2 and W 1 7 W 1 7 n1 n At 2 4W W At T 4h1 4h2 Ah2Ah3 Ah1Ah3Wn1 7 Wquot At 2 AMAMAMWW 7 Wquot Hi d 3 8 For the special case of A iAh the general procedure coincides with the alternating direction equation with Fj properly chosen for the two space variable problem so that the method for d 3 gives the correct extension of the alternating direction concept to higher dimensions There eXists a quite general convergence analysis for this alternating direction method when d 2 It is not necessary that Q be a rectangle and the convergence argument extends not only to variable coefficients but also to semilinear parabolic equations For higher dimensions there are fragmentary convergence results If as in the case of the heat equation on a cube a separation of variables procedure can be applied convergence can again be established This requires that the the operators AH commute which then requires that Q be rectangular Alternating direction versions of other fundamental differencing procedures such as the backward difference method can be constructed in a similar fashion Indeed for many practical problems having not very smooth solutions the corresponding alternating direction backward method gives more satisfactory computational results than the one above based on Crank Nicolson since the backward difference method attenuates the more highly oscillatory components of the error quicker than the Crank Nicolson based procedures do Improved Accuracy for Alternating Direction Methods We have shown that the error in the alternating direction method has the same asymptotic order of accuracy as the Crank Nicolson equation namely 9At2 l hz However numerical experiments will show that the actual error can be considerably larger for any fixed choice of At and h if the derivative uxxyyt of the solution of the differential problem is large as it can well be for practical problems Fortunately there is an inexpensive modification that changes the perturbation term we shall consider only the d 2 case here At TAh1Ah2Wn1 7 w into a term that is formally OAt3 Essentially we should like to modify the right hand side of the standard alternating direction method to be At f TAMAle 7 Wquot if we could do this we would recover the Crank Nicolson equation Since we do not know Wnll we cannot but we can approximate the modification by using the most recently obtained timedifference of the approximate solution thus let us replace the right hand side fnl by f lt f Ah1Ah2Wn WWI This can be accomplished trivially just replace the right hand side in the n 1 equation and continue using the simplified equations for n 27 7 d This requires a few extra arithmetic operations per grid point to evaluate the right hand side but no extra operations to solve the single space dimensional algebraic systems If this modification is made then elimination of the intermediate Wr39l f39s leads to W 1 7 W 1 quot1 n At n1 n n71 nl AW W 7 4h1Ah2W 2W W f 2 At 2 4 and the new perturbation term involves the second difference in time of the approximate solution multiplied by At Thus the local truncation error now takes the form 63 0umAt2 0uW umw 0uWAt3 For the model heat equation a simple modification in the convergence analysis given above proves that the error associated with the modified alternating direction method is again 9At2 l hz but experiments with the heat equation and more general parabolic equations have shown that actual errors for it are almost identical to those for the Crank Nicolson equation and are significantly smaller than those associated with the original method with an unmodified right hand side The idea of adding a term to the right hand side that will nearly cancel the perturbation term or more precisely to transform it into a term of higher order in At can be applied to other alternating direction procedures and to a collection of so called locally one dimensional methods due to several Russian authors Alternating Direction Methods as Iterators 1 The Jacobi iteration can be interpreted as an application of the explicit finite difference method for the heat equation to the elliptic difference equation to find the solution of the elliptic problem as the steady state solution of the parabolic problem 2 Any of the difference methods for parabolic problems could be substituted for the forward difference scheme however nothing would be gained of practical value if either the backward or Crank Nicolson procedures were applied since the algebraic equations generated at each time step for each of these are effectively of the same difficulty as those for the elliptic problem 3 This is not the case for alternating direction methods since only singlespace dimension algebraic equations arise for them 4 For the Laplace difference equation the Jacobi method corresponds to taking a time step equal to the critical CFL time step that is the virtual time step is given by At h24 With this step the coefficients of the low frequency error components are reduced quite slowly leading to the requiring of OM iterations to reduce the error by a factor of two Now the lack of a CFL constraint for alternating direction methods gives us the opportunity to take larger time steps thereby approaching the steady state somewhat quicker 5 We shall consider two techniques one employing a fixed virtual time step and the second a cycle of different virtual steps The fixed step iterator is provany applicable to a wide variety of elliptic problems though here we shall limit ourselves to the two dimensional model problem The technique employing a cycle of parameters has been analyzed completely only for commuting Aim though it has proved effective in much less restrictive situations An Optimal Fixed Virtual Time Step Consider the discrete Dirichlet problem iAhWij fij XiJ39EQh7 WU 07 X68917 where Q 012 Let v3 be an arbitrary initial guess subject to satisfying the boundary condition Then choose a virtual time step At and iterate using the original alternating direction algorithm Let n W 7 Vn denote the difference between the solution of the difference equation and the nth iterate 6ln1 7 6In At 1 At 7 EAh enll l e l Tz yz y enll 7 e 0 We shall again make use of an eigenfunction expansion to analyze the behavior of e Let N h 1 and set 1 4pm dimicy sm PX E 539 7 4 7rph 17 At A i 2 7 7p P hZSm 2 pp 1pAt If e0 has the expansion N71 0 Z 0 e 7 quLme pq1 e has the the expansion N71 n n e 7 E quappq pq1 where pq1N71 7 n 7 n71 7 n 0 qu pppquq 7 pp051 qu39 Now pp is a decreasing function of p with l AlAfi 717N71 1 if gt gt 777 gt11At 4quot le 1N1 gt 71 The optimum fixed At is that for which p1 pN71 thus a trivial calculation shows that 1N1Afopt2 17 from which it follows that h At 7 opt 2 Then liwh2 xixli h 3901 17rh2 7quot sothat maxlpppql p7 q 17 N 71 g 17 27rh l 9h2 if Atopt is employed If so then ignoring the 9h2 term llenllllo S 1 27rhllenllo S 1 27rhquotlleollO Thus I 2 1 0g 0 ngti Im liesthat e lt7 e m p H ll072ll iioi in contrast with the requirement that n 0h 2 in order that the Jacobi iteration assuredly reduces the difference between the solution of the Laplace difference equation and the nth iterate by a factor of two thereby leading to an 9N3 work estimate The alternating direction iteration has been analyzed here only for the Dirichlet problem on a rectangle however a more sophisticated argument can be given to show that the 17 Ch error reduction factor holds for variable coeficients and on a nonrectangular domain for iterating with a fixed virtual time step Iteration Based on a Cycle of Virtual Time Steps Both the Jacobi iteration and the iteration described just above approximate the solution of the initial value problem for the heat equation with initial values given by the initial guess for the steady state solution but if only the steady state solution is of interest then a different strategy can prove to be more efficient Let us consider taking virtual time steps Af17Af27 and set 17 ApAt ppAt m 6 711 VAt If the iteration is carried out employing Atk in the kth iteration then qu PpAtnPqAtnEZ1 rl 1ppAtk rl 1pqAtk qu39 Choose Atp such that PplAtp 0 that is let 1 h2 Ati77 1N71 p p 4sin27r7 27h p Then 625 07 n 2 minp7 q7 and we have obtained except for rounding error the exact solution of the Laplace difference equation on 012 in N 71 iterations each being stable Theorem Alternating direction iteration can produce the exact solution with exact arithmetic of the Laplace difference equation on the unit square in not more than N 71 iterations where N h l Since the number of arithmetic operations needed to carry out a single iteration is 9N2 this means that an exact solution can be obtained in 9N3 operations whereas Gaussian elimination would require 9N4 and the optimal fixed parameter alternating direction iteration would require 9N3 to reduce the initial error by a factor of two In the argument above the effect of an iteration with parameter Atp on modes 4pm with r or 5 close to p has been ignored Let us take up consideration of these modes Set liAAt p7Af 71 Ami and let 6 6 01 A simple calculation will show that 176 16 7 7 lt lt 7 pAt 6 66 If and only If 1 6 AAt 1 6 Recall that 1 z 71392 and ANA x 4h so that we need a cycle 73917 7392 7Tm such that gt If E 7r274h 2 there eXistsj7 1 j m such that lP7739jl S 5 Let us define 7 as follows Set 1 71392 and determine 7391 so that pp1T1 6 01 1M1T1 17 1 7 6 7 or 739 7 pm 1 Hmn 1 1 m 01 Then find p2 so that JOHN1 76 7 1M2T1 5 0012771 1 MT17 which implies that 716 if 16 2 quot2 176 n 176 quot139 Proceed recursively by requiring that pm7n5 and pw17n5 Then Thus 16 2 16 2quot 2 l f6 H m 7 76 2 176 2H 2 T1 1 l Tj 1lt16gt T 39 The length m of the cycle is determined by choosing the firstj such that m 2 4h 2 consequently choose m such that 2m72 176 l 4 wzhz39 Thus let 16 1 1 4 gt 7 7 m71ltlogli6gt lt2loghlog7r2gt With this choice of m HemHo S 52H80H07 since qu S n21pPTlpQ7 l quv and both pp7j 6 and pq739k 6 for somej and k not greater than m Next let us consider the choice of 6 We know that the the difference between the solution of the differential problem and the discrete one is 9h2 so that reducing llenllo below say hzlleollo is not producing a better approximate solution of the differential problem So a first thought as to picking 6 might be 6 h so thatjust one cycle of parameters is used and the iteration terminated For 6 h 16 176262h log and a number that exceeds N71 h 1 71 Clearly this is a poor choice for 6 since it leads to more arithmetic operations than required to find the exact solution using the parameters Atp p1N71 Instead let us try 6 so that Hemiio eoiio Then mXCIog7 a number far lower than the 0h 1 required by the optimal fixed parameter alternating direction iterative method and even more so than the OM iterations for the Jacobi method In order that iieniio S hziieoiim it suffices to repeat the cycle 7177739m M times where 2 l ighz or 2 Ogh 2 og2 With these choices for 6 and M then the total number of alternating direction iterations required to reduce the error by 9h2 is m M O ltltIoggt2gt O ltIogN2gt7 and the total number of arithmetic operations generated by this iteration is 0N2Iog N2 The results concerning the original alternating direction method both for the heat equation and as an iterator for elliptic equations are due to Peaceman Rachford and the lecturer and were introduced in the summer of 1954 by the lecturer at a meeting of the American Mathematical Society at the University of Wyoming III 3 Lecture 10 Iterative Methods Basic Iterative Methods Consider the system of linear equations AX b where A aban is a specified real matrix and b bnX1 a given vector respectively Given the initial guess XOld to get a better approximation to the exact solution X of the linear equations it is desirable to update old 7 x Xnew Xold X ixoldgt X5 XOld by adding the error e X ie Unfortunately this is equivalent to solving the equations On the other hand it is inexpensive to compute the residual r A ltX 7 XOldgt b 7 Axold It is natural to update the current approximation by B lr Xnew Xold Bilr Xold Bil b 7 AXoldgt7 where B is an nonsingular matrix to be chosen Repeating this procedure defines an iterative method xltk1 x 341 7 M k 1 2 for a given initial guess X0 Now the issue is how to choose B so that B ly can be computed inexpensively for any y E R and that B 1 approximates A 1 as closely as possible Usually these two criteria are inconsistent since the better approximation of A 1 that B 1 is the more expensive is the solving of the equations By r We begin with the decomposition A L D U7 where L D and U are the strict lower triangular diagonal and strict upper triangular parts of A respectively That is 0 0 0 a O O L 21 Ddiaga an1 ann71 0 0 312 a1n U i 0 0 an71n Taking B D and 0 gives the Jacobi and Richardson methods respectively which were discussed earlier I One step Jacobi iteration Xnew Xold 071b7AXold ll One step Richardson iteration Xnew Xold wDilb7 AXold The componentwise form of the Jacobi method is as follows 1 n Xinew Xiold E bi 7 Zaijxjoldgt Ii 11 for I 17 n When we compute the new approximation to the i th component the current approximations to the first if 1 components are Xnew forj 17 ii 1 If we make use of those we then have the GaussSeidel GS method iil n GS new old 1 new old XI XI 1 f bl anxj Eauxj H j1 jl for I 17 n Based on the GS method we have approximations XiOId and XGS for the i th component By taking an appropriate weighted average of them it leads to the Successive OverReIaxation SOR method GS 50R new 17 wxod wxi Xi Xi I iil n LU 17 wXOd i 7 b 7 E ainneW 7 E aijxOld aii J J j1 ji1 The matricial forms for the GS and SOR are as follows III One step Gauss Seidel iteration Xnew Xold 0 L71b7AXod39 IV One step SOR iteration Xnew Xold 440 wL71b7 Axold7 where w is the relaxation factor to be chosen When w 1 SOR is the same as GS The SOR iteration is called underrelaxation for w lt 1 and overrelaxation for w gt 1 The Jacobi iteration does not depend on the ordering of unknowns but GS and SOR do There is a natural ordering for the unknowns in elliptic difference equations ie in a two dimensional domain order the unknowns from left to right and from bottom to top However the ordering can have an influence on the rate of convergence of some iterative methods For the five point difference scheme for the Poisson equation another natural choice is red black ordering ie the chess board like coloring of the grid of unknowns Then its GS iteration for red black ordering is given by gt For the red points 1 xgew I Xiojiij Xicirliij Xicjjigl X gl hZ j gt For the black points 1 new new new new new 2 Xi 7 Z xi71 Xi1J XIJ1 XiJ1 h fu Some Convergence Results Consider the general iterative method 7 7 NH xkB 1b7Axk I7B 1AxquotB 1b k 1 2 If X denotes the exact solution of AX b then it satisfies X X i B 1b7 AX I 7 B IAX i B lb Set em X 7 XU then the difference of the two equations above gives the error equation NH I 7 B lAek I 7 B lAke0 Denote the iterative matrix by M I 7 B lA and its spectral radius by pM maxlgign lMl where AM are the eigenvalues of M Theorem The general iterative method converges for any X0 and any b if and only if pM lt 1 Proof Since X0 is arbitrary by the error equation convergence is equivalent to Iim Mk 0 kaoo which is equivalent to pM lt 1 Assume that M has a complete orthonormal eigenbasis gtf391 ie M gti i i7 177n7 m lt15quot 5039 where 7 is the 62 inner product Expand the initial error with respect to the above eigenbasis 90 X X0 Zea1 i1 Then the error for the kth iterate is 5 00 ElliWi i1 thus the amplitude of the ith component in the initial error is reduced by a factor in k iterations By the orthogonality of the 1 llemll S PUlllkllewllh where is the 2 norm Now pM is called the convergence factor and the convergence rate is the natural logarithm of the inverse of the convergence factor RUVI a lnpUVI ie the minimum number of iterations required to assure reducing the initial error by a factor of le The next lemma gives a necessary condition for the convergence of the SOR method Lemma If the 50R method converges then to E 07 2 Proof Denote the SOR iterative matrix by MSOR I 7wD wL 1A Then MSOR DWLrliilwwa AdrianoWU where Z D lL and U D lU Since detl wl 1 then the characteristic polynomial of MSOR is p det lwlI7MSOR detxw71lwZwl1 Hence i xiimsom p0 det w 1Iw0 w 1quot i1 which implies that pUlISOR 2 lw 7 ll hence the lemma follows from the previous theorem Theorem If A is symmetric and positive definite then PMSOR lt1 fora0ltwlt2 Proof Since A is symmetric then U Lt Then D wLek1 1 WAD ithek Set 500 ek ek139 Straightforward calculations give 7me uAeW and uAeWI 17wD6k7th6k Taking the inner products of As and AeU H with em and ek1 respectively differencing the two inner products and simplifying leads to Auk 9kAek17 emu aunt 5m 2 0 w This indicates that the sequence Aea 7 ek is monotone decreasing and bounded below by zero Hence the limit of AeU 7 em exists which together with the last inequality implies that 060 7 600 has limit zero Since D is positive definite limkH00 60 0 Then the positive definiteness of A implies that em a 0 and convergence is assured An Optimal Relaxation Factor for SOR Let us determine the optimal relaxation factor of the SOR method for the finite difference equations for the model Dirichlet problem If A ithh then its ij th equation is 2 Awlij 4Wiv39 WHJ Wi1J W071 WiJ1 7 fiv The matricial form of the SOR method is the same as that described previously ie WW1 MSORWW wD wL 1h2f where the iterative matrix is given by MSOR D wL 117wD in To estimate the convergence factor we need to estimate eigenvalues of MSOR which satisfy the fundamental equation for the eigenvalue and the eigenvalue v MSORV v Multiplying the basic iterative equation by D l wL gives 17 wD 7 wUv D wLv7 whose component wise form is 41 MW wVi1J V041 4ViJ wlvieu mill Considering the eigenvector of the form v Vij AzJ leads to LU LU 1 WZiJ1Zi1j l zij1 Zij 7ZZi71xj 21921 1 l Zij 1Zi1j ZiJ1 Ziilj Zixjil This means that u Ali is the eigenvalue of the corresponding 5 I I I uJA Jacobi Iterative matrix Lemma The eigenvalues of MSOR are determined by Amimi 1710 where u is an eigenvalue of the corresponding Jacobi iterative matrix Proof We have just seen that the eigenvalues of MSOR satisfy the required equation of the Lemma Aw71 1 wA Z is the eigenvalue of the Jacobi iterative matrix with its associated eigenvector z 2 which in turn implies that A is the eigenvalue of MSOR with the associated eigenvector VViJltZiJgt Conversely if A satisfies this equation and A 7 0 then u It is easy to see that O is an eigenvalue of MSOR when w 1 Otherwise we get i Adlai xwzluz 7 4w 71 7 2 It is then easy to show that ii lt 1 if and only if M lt 1 and w E 07 2 Theorem If the Jacobi method converges then the SOR method converges for w E 07 2 Now we are ready to determine the optimal relaxation factor Assume that M lt 1 and w E 07 2 When 4w 71 wzpz ie 272i17p2 2 w eo2 2 1 17p2 there is only a real double root to the characteristic polynomial 1 7 7 2 5 i ilwwiw1 t 2 1 17p2 When 4w 71 gt wzpz there are a pair of conjugate complex roots with moduli u i f 4 4 W m an increasing function ofw When 4w 71 lt wzpz there are two distinct real roots and the largest one has absolute value 1 1 1 134ij wawn which is a decreasing function of w since lk naxl lt 0 1 Hence the largest root as a function of w lAEnaXl reaches its minimum at the w given by the first case above Thus lei172 2 2 when w HuiL 1 17 minimaxl LL It is obvious that min lAmaXl is a increasing function of hence it LL reaches maximum when M takes maximum value 1 vlil brznax M max min POP SOR FL w i maxi 1 m with the optimal relaxation factor 2 we t p 1vlil br2nax gt1 We computed the spectral radius of the Jacobi iterative matrix to be h th umax cos7rh 17 2sin2 7 z 17 WT 17 0h2 Hence the spectral radius of the SOR iterative matrix with the optimal relaxation factor is 17 17 cos27rh M 7 Poptl SOR 1 licos27rh lei1717L2 2 1 17 lt17quot2h2gt2 22 2 172m 17 Oh 22 This means that an SOR iteration should require 0N 1log1n iterations to reduce the residual by a factor of 7 thereby requiring 9N3 log1n arithmetic operations Lecture 1 1 Conjugate Gradient CG and Preconditioned CG PCG Iteration We will treat the conjugate gradient CG method and the preconditioned CG for a symmetric positive definite system of linear equations Conjugate Gradient Iteration Assume that the matrix A is symmetric and positive definite It is known that the system of linear equation Axb is equivalent to the minimization problem find X E R such that x min y with fy Ay7yeb7y yER where 7 is the 62 inner product Given an initial approximation X0 E R an iterative method finds successive approximations X00 E R according to the recursion Xk1 Xk akdk7 k 1 2 7 where 10 E R is the search direction and 04k gt 0 determines the step length Different methods differ in the choice of the search direction and the step length For example the search directions for the GaussSeidel and gradient methods are the coordinate direction and the direction opposite to that of the gradient respectively The equipotential surfaces of fx form an ellipsoid and it is almost always a long thin American football in shape Consequently if one tries to minimize fx by taking the steepest descent search directions the approximate solution just bounces back and forth reducing f trivially on each iteration This clearly can lead to accepting a solution that is far from the actual minimum So some control must be imposed on the selection of the search directions and the conjugate gradient method provides one such control The method will first be defined in the form leading to the most efficient algorithm for its application then it will be derived in a second fashion that more clearly indicates its properties Once the search direction is chosen the optimal 04k satisfies fxk l akdm mig fxk l 041k 04gt This is a one dimensional optimization problem whose solution is No7 dk 07 7 k Adk7dk r 7b AX 04k The search directions for the CG are chosen as follows 1 10 r0 2 For k 07 17 the k 1 th search direction of the form dk1 rk1 Bkdm is to be A orthogonal to dm ie dk17 dkA Adk17 dk 07 so that Ark17 dk Bk Adm dk 39 This can be summarized in the following algorithm The Conjugate Gradient Algorithm Given the initial approximation X0 E R let 10 r0 b AX07 For k 0717 let M k7 dkAdk7 dk7 Xk1 X00 0441007 rk1 b 7 AXk17 5k Ak17 dkAdk7 dk7 dk1 rH1 BMW Theorem The conjugate gradient algorithm has the following properties AMT7 11 O for j 0717quot397 k Adk17 11 0 for j 017 k rk17 0 0 for j 01 k The proof of these orthogonalities will be given shortly but first let us consider some practical aspects of the CG algorithm The algorithm seems to require three matrix vector multiplications per iteration but by using the orthogonalities stated in the theorem the number of multiplications can be reduced to one Since 10 NO Bkildwiny ak k7 kAdk7 100 Since A is symmetric and rk1 b 7 AXk1 rm 7 akAdU then Ak17 dk rk17 Adk rk17 k k1ak k17 k1ak which implies that Bk ra il7 rk1rl 7 rm Hence we can obtain the following version of the CG algorithm Given the initial approximation X0 E R set 10 rm b7 AX0 and then for k 07 17 compute ak No7 kAdk7 dk Xk1 Xk WM rk1 rm WAN 5k k17 k1rk7 k dk1 rk1 gkdm VVVVV Note that the only matrix vector multiplication is given by Adm There are two inner products three multiplications of a vector by a constant and three vector additions per iteration Now let us investigate the properties of the CG iteration The fundamental concept in this iteration is the A orthogonality of the search directions So let us begin by choosing H ado k k M d ngdoxdondj where rm b7 AXk as usual Then we have the desired A orthogonality Adkgtd0gt 0 j 0k71 Since 10 rm A17 d0 o o Ad07d0 d 7quot 0 Adar 10 A17 10 7 Then by induction and the symmetry of A for Z lt k Adk7 10 Ak7 10 Ak7 11 7 as was to be shown if Ak7 dj j z Adj7djAd 7quot l ll 39 0 1 6MArk 10 0 ll Now let UU Spand07 7 dm clearly UU Spand07 7 104 rm Then dimUU 1 unless rm 0 in which case XU X the solution of the linear equations Since the linear space has dimension n the CG iteration with exact arithmetic produces the solution of the equations in not more than n iterations Since Xk1 X00 l akda the optimization problem forces the choice 04 rm7 dl Adl 7 dm thereby coinciding with our first definition of the CG algorithm It follows trivially that rk1 rm 7 akAdU Next the A orthogonality of dm immediately implies that dk1 rk1 l Bkda with Bk evaluated as earlier Thus we have shown the equivalence of the two ways we have used to define the CG procedure Left to be shown is that rm7 rm 07139 07 7 k 71 for which it suffices to show that rm L U07 j lt k First rm rm rm rm 7 emuAdm d 0 by the definition of ak1 Next k7 11 k17 djak71Adk17 11 0 j 07 H k2 by induction and A orthogonality Forj k7 1 k7 dk1 k17 dk1 ak71Adk17 dk1 k4 k17 dk1 k1 k17 Zyjdm 07 10 and we have shown that rm L Uk 1 so that rm7 rm 07139 lt k This completes the demonstration of the properties of the CG iteration listed in the theorem stated earlier The Convergence Factor for CG Iteration Note first that ltAxk 7 W 10 7 Aw 7 W 10 j 7 0 k 71 Thus X00 7 X0 is the projection of X 7 X0 onto Spand07 dk 1with respect to the A inner product That is X00 7 X0 is the best approximation of X 7 X0 in Spand07 dk 1with respect to the A norm iix7xltkgtiiA 7 W 7xlt gt7 W 7xlt gtiiA W 7xlt gt7yiiA for any y E Spand07 104 It is easy to show that Spand07 dk 1 Spanr07 rk 1 Spanr07 Arm7 Ak 1r0 These relations follow from the representations k I k I I k ZakJAjr07 dk Z mi4amp07 j0 10 which can be shown by induction k Let y Z ajAjX 7 X0 then 11 k x e M e y I e Z aw x e M MAW e X 7 j1 where pkt is a polynomial of degree k satisfying pk0 1 Denote the eigenvalues of A by AJ and order them from the smallest to largest 1 S S n Let 7 denote the collection of polynomials of degree k and let 7 be its subset taking the value one at t 0 For any pk 6 75k llPkAX X llA 7 0 lgagnlpkh llx X llA llX XmllA lt t 7 0 A1r t22ltAnlivkllllt X llA Optimize over pkt 6 75k 1 I f W i W i A result from approximation theory states that k e 1 t 2 Wk AEEZ Aankl l lt E1 where a is the spectral condition number of A defined by r n1 Theorem Let X and X00 be the exact solution and the k th CG iterate respectively Then k 71 llx XkHA g 2 ltg 1 th X0HA Preconditioned Conjugate Gradient Iteration The convergence factor of CG iteration depends on the condition number of A The condition number of the five point discrete Laplacian is Oh 2 This implies that its CG convergence factor is 17 Oh which is close to one when the mesh size is small and is approximately the same as for SOR Here we shall study preconditioning techniques which lead to new equivalent linear systems with smaller condition numbers Again consider the system of linear equations AX b7 A symmetricpositive definite The preconditioning technique is based on replacing AX b by B le B lb7 where the preconditioner B is an approximation to A Assume that B is also symmetric and positive definite To be of practical value it must be inexpensive to solve the system By d and B lA must be better conditioned than A Note that B lA is not symmetric with respect to the 62 inner product but it is self adjoint with respect to the B inner product 7 B B lAX7 yB AX7 y X BB 1AyxB 1AyB Hence the CG method for the preconditioned system can be easily described by replacing the 2 inner product in the CG algorithm by the B inner product Denote by rm b7 AXUO and 2m Bilrm the residuals for the original and the preconditioned systems respectively Then ak 2007zkBB1Adk7100 W zltkgtAdltkgt7 m 5k 2k17 zk1Bzk7 200 k17 zk1rk7 ZOO The algorithm for the preconditioned CG iteration PCG can be formulated as follows The Preconditioned Conjugate Gradient Algorithm Given the initial approximation X0 E R set rm b7 AX0 and 10 20 B lrm Then for k 07 17 compute ak No7 zkAdk7 dk Xk1 Xk WM rH1 rk WAN zk1 371M 5k k17 zk1k7 2 dk1 rH1 gkdm VVVVVV The only extra calculation over the CG iteration is the solution of the 82 r equations which we assume to be not too difficult Since 1B 1AX7 XB AX7 X llxllA the convergence factor for the preconditioned CG is then R W71 llX ixmllA S 2 W W 7 X llA where 5 is the condition number of B lA How to choose a good preconditioner is still the topic of current research For the system of linear equations arising from the discretization of elliptic partial differential equations it is known that multigrid method generates a good preconditioner as well as being a superior iterative method on its own Multigrid Methods This topic will be the subject of a guest lecture by Professor Craig Douglas next week Lecture 12 Galerkin Methods for Elliptic Problems Some Typical Planar Finite Element Spaces Let us indicate a few possible finite elements that can be applied in the approximation of second order elliptic problems in the plane An element will be called conforming if it can used to generate a function in H1Q where Q is the domain on which the elliptic problem is defined If not it will be called nonconforming it is most common to treat conforming methods though nonconforming methods have serious applications as well The examples will be presented graphically In the figures the degrees of freedom will be indicated as follows 0 A black dot will indicate that the value of the polynomial is to specified at that point 0 A small circle about a point indicates that the gradient is to be given there 0 Two small circles about a point means that both the gradient and the three second derivatives are to be specified 0 A short segment orthogonal to an edge of the figure indicates that the normal derivative must be specified there The Conforming Piecewise Linear Triangle The simplest two dimensional finite element space is given by Mg0 v 6 com v he P1T T 6 Th where 711 T j 17 N is a partition of Q into nonoverlapping triangles of radius not greater than h Figure Conforming PiecewiseLinear Triangle A family 71 hk a O of partitions is called quasiregular if diam Tm lt lt max max 7 00 k TmTneThk diam Tn 0 and the minimum vertex angle in any triangle in U k is greater than some 0min gt 0 The family of partitions is called shape quasiregular if the minimum angle condition is met Usually we shall not take the trouble to embed a partition in a family but the existence of the family should be understood The piecewise linear element is defined by v l7 v E P1T Span1xy vqj 093139 172737 where qj7 j 127 3 is the set of vertices of T Clearly v reduces to a linear function in the arclength on each edge so that its edge values are uniquely determined by the corresponding vertex values Thus if Tm and Tn share an edge em the piecewise linear functions determined by the vertex degrees of freedom coincide on em and continuity across the edge follows so that the global extension of this piecewise linear interpolation of the degrees of freedom is in H1Q Also continuity is forced in the neighborhood of a vertex of T Thus we have constructed a conforming finite element space Mao El A Nonconforming Piecewise Linear Space Shift the degrees of freedom from the vertices of T to the midpoints of its edges Figure Nonconforming PiecewiseLinear Triangle Clearly v E P1T is again uniquely determined on T by these degrees of freedom however only the average values of the interpolations in adjacent triangles will necessarily agree on their common edge so that the approximation space generated by these parameters will not be in H1Q Thus the resulting finite element space is nonconforming so that a standard Galerkin procedure cannot be based on this space Later we shall see that a useful method can be based on it The Conforming CO Piecewise Quadratic Triangle Let v be P2T Span1xyx2xyy2 so that there must be siX degrees of freedom associated with T and v le the restriction of v to the edge e is a quadratic function of the arclength on e In order to require continuity in the neighborhoods of the vertices it is necessary to impose value degrees of freedom at the vertices Now to obtain continuity across e we must impose another value degree of freedom on e which we naturally place at its midpoint Figure Conforming PiecewiseQuadratic Triangle Lemma The siX degrees of freedom listed above uniquely determine a quadratic function on T Proof Denote the three edges of T by e17 e2 e3 where ej is opposite to the vertex qj and let denote the linear function giving the distance to the line passing along ej Now it suffices to establish uniqueness of the quadratic function So let v E P2T vanish at the siX loci of degrees of freedom Then v vanishes on el and eg hence vx cd1xd2 Let r3 be the midpoint of eg neither d1 nor d2 vanishes at r3 Since vr3 0 c O and v vanishes identically as was to be shown A local basis for the piecewisequadratic triangle can be constructed easily Consider a basis element that can be associated with the point r3 d d br3x 1xl 20 d1r3d2r3 b3 vanishes at the other five support points for the degrees of freedom and equals one at r3 as desired Next let d12x be the distance of x from the line passing through r1 and r2 Then the interpolatory basis function for 73 is given by d3Xd12X bq3lx 7 d3q3d12q339 The remaining four basis elements can be constructed similarly A Conforming CO Piecewise Cubic Triangle The dimension of P3T is ten and it takes four value degrees of freedom on each edge located at the two vertices and at two interior points on the edge to assure that the global interpolation based on value degrees of freedom leads to a continuous function The one remaining degree of freedom can be taken as the value at the center of the triangle Figure Conforming PiecewiseCubic Triangle Lemma The ten degrees of freedom indicated above determine a unique cubic Proof Again let us establish uniqueness So if v is a cubic vanishing at the points associated with the degrees of freedom v vanishes on all three edges so that vx c3931dx Then the vanishing of v at the center of T implies that v is identically zero The interpolatory basis on T can be constructed normalizing the product of distances to lines as in the case of the quadratic element In each case either the distance to two edges and one distance to a line associated with the center and two boundary nodes are employed or for the center all three edge distances are used Another Conforming CO Piecewise Cubic Triangle Figure A Second Conforming PiecewiseCubic Triangle Knowing the components of the gradient at a vertex implies that the derivative tangential to any edge associated with the vertex is evaluable Thus these degrees of freedom determine the restriction of the cubic function on T to the edges Thus the uniqueness argument is unaltered consequently these degrees of freedom also suffice to define a cubic element for which the global interpolation operator is continuous so that this element generates a conforming finite element space The basis element corresponding to the center is again the cubic bubble function cd1xd2xd3x The construction of the remaining basis elements can be made by treating the isocees triangle Tmf with vertices 07 07 01 and 17 0 and then mapping this reference triangle affiner onto the desired T and transferring the basis to T as follows Let 1 be the affine map such that Tref Ti if 86 is the nodal basis on Tmf then the basis on T is given by 66T gt 1oww 6 Eu It is actually possible to employ both of these conforming cubic elements in constructing t he global finite element space by using transition elements to bridge the change from one to the other These elements take the forms indicated below It is easy to see that the same uniqueness argument shows that both transition elements are well defined Figure Transition Elements for Conforming Cubic Triangles Some Conforming Higher Order Triangles Let us consider several quintic finite elements for polynomials of higher degree there can be many possible choices of degrees of freedom Note that the dimension of P5T is twenty one and that it takes siX degrees of freedom on an edge to insure continuity across the edge Figure Conforming Quintic Triangles These two finite elements are well defined Note that vanishing boundary degrees of freedom lead to P5X d1Xd2Xd3XP2X7 and the interior degrees of freedom kill a quadratic polynomial Figure More Conforming Quintic Triangles These two elements are again well defined as the boundary nodes imply the vanishing of the element on the boundary for vanishing data on the boundary then the interior data again defines a linear function Figure Other Conforming Quintic Triangles The vertex degrees of freedom in either of this pair give uniqueness on the edges The three internal degrees of freedom suffice to determine a linear function which is insufficient to force uniqueness of the quintic Thus a different argument must be given to see that we obtain a well defined quintic The simplest if least attractive argument is to write out a quintic polynomial in the form VX Z aaxalya27 04 0417042 lal 041 0427 lalS5 nnA nnll I39Virl inn Annw m AF meA m A ii The degrees of freedom at the origin force aw O for iai 2 The degrees of freedom at 10 kill the terms X37 X4 and X5 those at 01 eliminate the corresponding y terms This leaves nine monomias and nine constraints a straightforward argument will show that the matrix generated by these nine degrees of freedom is nonsinguar which proves the uniqueness The details of this argument will be omitted Some Conforming Cl Piecewise Polynomial Elements There are no conforming C1 triangular elements of low order Perhaps the simplest such element is the quintic triangle defined in the following figure If the indicated degrees of freedom a vanish then the quintic polynomial v vanishes on the boundary along with its normal derivative Thus VX C d1Xd2Xd3X27 so that v is either a sixth order polynomial or it vanishes identically Since v is a quintic it must vanish This shows that the element can be used to construct a conforming C1 finite element space Figure A Conforming C1 Quintic Triangle It is possible to construct a conforming C1 element based on cubics but only by building a so called composite element Consider a reference triangle and divide it into three triangles by introducing an interior vertex and connecting it to the three original vertices The object of the CIough Tocher triangle is put together the thirty available degrees of freedom from a cubic in each of the three subtriangles in a way that connects the three cubics across the internal edges in a continuously differentiable fashion and that allows the specification of boundary degrees of freedom that define the function and its normal derivative on the boundary Figure The CIoughTocher Composite Triangle Note that there are six degrees of freedom specified at each external vertex since the three at each apply to two triangles Thus there are twenty one boundary degrees of freedom leaving nine to impose the internal consistencies across the internal edges Equality of the values and gradients at the interior vertex imposes six conditions not nine and the final three are required to impose equality of the normal derivatives at the midpoints of the internal edges After taking into account the interior structural requirements for the element there are twelve degrees of freedom that enter the calculations in a Galerkin approximation Conforming Elements on Rectangles The standard construction of conforming elements on rectangles is by means of taking the tensor product of the single space variables finite elements studied in relation to the two point boundary problem The bilinear biquadratic and bicubic conforming C0 elements are pictured below Their bases are conveniently given by taking the tensor product of the one dimensional bases which also demonstrates that the illustrated nodal degrees of freedom are independent Figure The C J Bilinear Biquadratic and Bicubic Elements It is equally easy to construct conforming C1 and C2 rectangular elements through tensor products of the Hermite cubic and quintic functions in a single variable The Hermite bicubic element is defined by specifying four degrees of freedom at each vertex xj VXj7 VXXj7 Vyxj7 nylle These degrees of freedom determine the Hermite bicubic function on the reference element R 711 To see this let the sixteen degrees of freedom vanish for v E H3 H3R On the edge X 1 v and vy vanish at its endpoints so that v E O on the edge Moreover VX and VXy vanish at the endpoints so that VX also vanishes on the edge and v contains the factor dxy2 where d is the distance function to the line X 71 0 Thus nil1diX2 l VX7 so that v would have to be a polynomial of degree at least eight to be nontrivial Hence v is identically zero as we set out to prove The obvious basis for the Hermite cubic space on the reference element R is given by the tensor product of two onedimensional bases for H3 thus if 1117 7114 VASL Vr75r7 then BH3 H3R w xi 1X k 113 k 1774 is this local basis An analogous oca basis can be constructed for the Hermite biquintic space Some Very Simple Three Dimensional Finite Elements Piecewise Linear Tetrahedral Elements The conforming piecewiselinear tetrahedral element has its four degrees of freedom as vertex values The nonconforming element corresponding to the planar nonconforming element has the facial midpoint values as degrees of freedom A PiecewiseQuadratic Tetrahedral Element The dimension of P2IR3 is ten Let the degrees of freedom for the piecewisequadratic tetrahedral element be the vertex values and the values at the midpoints of the edges These degrees of freedom define a conforming element Note that vanishing degrees of freedom imply that the quadratic vanish on all four faces as the function reduces to a quadratic in two variables on each face and the degrees of freedom associated with a face are exactly those of the planar quadratic element Hence not to vanish identically would require the polynomial to be at least quartic Tensor Product Elements on Cubic Elements Tensor product elements in three or n dimensions can be constructed on the unit cube from any of the elements defined on an interval in the same fashion as for the plane They can be extended using linear transformations to elements having flat faces Lecture 13 An Approximation Lemma One version of the Bramble Hilbert Lemma will be stated below The proof will be delayed The version to be given here is by no means the most general form of the lemma as we shall relate the statement of the lemma to the conforming finite element spaces described previously A partition Th of a two dimensional domain Q will be called shape regular if there eXists a positive number V such that the following conditions hold gt For each triangle T 6 Th with vertex angles 04 17273 min minagt Ten quotV gt For each rectangle R 6 Th with side lengths hX and hy hX hy if gt Sl39 lhy hXV In dimension three interpret a vertex angle as the solid angle at the vertex for a simplex and take the ratios of all three sides for a rectangular element The partition 7 Tk is said to be quasiregular if it is shape regular and there exists p gt 0 such that mkin diam TR 2 meaXdiam Tk ph Quasiregularity can be defined equivalently by requiring that if k is the diameter of the circle sphere inscribed in Tk Ck mEXd39amTk h and 1mm Zpgt0 Let Mh C C0Q be a finite element space over the shape regular partition 7 based on one of the triangular or rectangular elements previously described Let dof degree of freedom alalv Wk is a dof for Mh in T E n maxlal In order that a degree of freedom of order n be evaluable for a function v it suffices to require that v E H illlg for some 6 gt 0 when Q is planar and v E H ilg when Q is threedimensional Let If denote a 1 6 or n as appropriate Let dimQ 2 or 3 and assume that PmTClhiT7 TETh Let v E HkQ where k 2 if then define the Mh interpolant th of v by requiring that dOthTJV dOfVlmTV7 T 6771 where the meaning of this constraint is that the Mh degrees of freedom of JV and v coincide on all elements T of the partition Th Then two conclusions of the BrambleHilbert Lemma are that for T 6 Th ijhz if g k m1 llV JVl aT l thV JV T g ClV and that llviJvllol hlviJvll Clvlkhk7 If g kg m 17 where the constant C depends on the parameter 39y related to the partition and on 6 in the planar case The approximation estimates above are valid when elements having a curved face are employed along the boundary 89 when these faces are smooth and these elements remain shape regular and have diameters bounded by h A Trace Theorem Two quite technical results concerning Sobolev spaces are that o If v E H1Q and if 39yv denotes its trace on 89 that is its restriction to 89 then weH wn and iwila gciiviim 2 where C depends on Q o If g E H a then g 39yG where G E H1Q Moreover lgl 30 can be defined by the relation 2 lgl Q infllGll1 G E H1Q and g39yG To be valid these two results do impose conditions on 89 which will not be specified here For the analyses in these lectures the following simpler trace theorem frequently is adequate Lemma Let Q 01 7 n 2 or 3 and let v E H1Q Then w E L28Q and there eXists a constant C such that l l llWllom S Cllvll gllvlliq Proof It will suffice to prove the lemma for v E C OQ since the constant C determined below will be independent of smoothness of v and C OQ is dense in H1Q For n 2 express v0y2 in the form my vim a X2v5yvX5ycs 0 Thus 1 1 1 v07y2dy VO7y2dxdy 0 0 0 1 1 S vxy2dxdy 0 0 1 1 1 1 8v 2 2 7 7 ll A 1 MW 8 835 ddedy l 3 1 llvll3n Ellelng ll 6 llVVllOQllVllOQ then 1 A V07y2 dy S llVll3 2 2llVllo 2llVVllo 27 and the lemma follows with C x 2 For n 3 express v0yz2 in an analogous fashion and proceed as above in this case C m 3 It is clear from the proof that the lemma is valid for Q a compact subset of IR having a piecewisesmooth boundary and we shall apply it under those conditions Galerkin Methods for the Neumann Problem Let us begin our study of Galerkin methods for elliptic problems by considering the Neumann problem 7VaVucu f XEQ7 8H 35 where 881 denotes derivation in the direction of the outer normal to 89 g7 x689 Assume that Oltat ax a7 OltQt cx c7 and that 04 minat7 cc and B maxa7 6 Also assume that the Neumann problem has a unique solution for every g 6 89 x Ham and that Hullz c iifiio llgll ml The HZ regularity does impose a constraint on the shape of Q but it holds for many domains that arise in practical computations lt suffices that 89 be in C2 or that Q be convex however this regularity does not hold for an L shaped domain It also imposes some smoothness constraints on a and c For convenience in our presentation we shall take the coefficients and the boundary to be however smooth necessary for the solution to have the regularity used in the analysis Thus the results to be stated and proved below represent the most rapid rates of convergence that can be associated with the problem and the finite element procedure under consideration We have chosen the Neumann problem to be our first problem since the Neumann boundary condition is natural when the weak form of the differential problem is based on the Dirichlet integral and can be imposed weakly whereas a Dirichlet condition is essential and most be imposed strongly and thereby leads to some technical difficulties that we can postpone treating The weak form of the Neumann problem is given by seeking u E H1Q such that AU7 V aVLhVV Gui V 5 V ltg7 Vgt7 V 6 419 where aVu7Vv should be interpreted as aVuVv1 and gv ds 30 Note that the bilinear form A is both coercive and continuous on H1Q ie for uv E H1Q AW 2 allulli lAU7vl S llulllllvlll ltg7 Vgt Let 7 be drawn from a quasiregular family of partitions of Q with h maxdiam TeTh Assume that Mh is a conforming ie Mh C H1Q finite element space such that inf HuixiiOJr hiiuixiil Ciuishs 1 s r 1 XEMh This corresponds to requiring that Mh D PT for all T 6 7L and to a somewhat more precise version of the BrambleHilbert lemma Then the Galerkin approximation Wh 6 7 to the solution to the Neumann problem is the solution of the finite dimensional problem AW17 v f7 v l g v v E Mh Theorem Under the assumptions laid out above there eXists a unique solution of the Galerkin approximation and ll Whllo hll Whlll S Clulshs for any 5 E 17 r l 1 for which u E H5Q Proof The argument is essentially unchanged from that given in the case of the two point boundary problem Let e u 7 Wh7 so that Ae7 v 07 v EMh Let X E Mh and take v X7 Wh O Au7 Wh7X7 Wh Auixxi WhAX7 Wh7X7 Wh7 so that allX Whll S Bllx 7 lllllX 7 Whll17 or 5 llX Whlll S lluixlll 04 Hence llui Whlll S llui Xlll llX Whlll S Mllui Xllh and it follows that Hui Whll M inf lluixlll Clulshs l xth Next let us apply the duality argument that we introduced for the two point boundary problem Note first that e E H1Q Thus there exists 4p 6 H2Q such that Av7w v7 8 v 6 H167 llwllz S Cllello Then for 11 E Mh llellg New New 11 S Mllelllllw fi lll so that 2 lt M 39 f 7 lt C h5 llello llelllwlethllw ll1 7 luls llelloy and the conclusion follows Robin Boundary Conditions The Galerkin method described above and its analysis can be extended easily to apply to a stable Robin boundary condition Let 7V aVu cu f x E Q 3 du g7 x E 89 81 where in addition to the hypotheses imposed on a and 6 above we require that 0 dx lt 3 Modify the definition of the bilinear form A to be Au7 v aVu7Vv cu7 v du7 v Then we see that the weak form of the Robin problem is again the seeking of u E H1Q such that Au7 v f7 v l g v v E H1Q Note that the error equation is the same as for the Neumann problem Thus if both ax and Cx are bounded below by 04 gt O and if the HZ regularity result holds the entire convergence analysis above is applicable without change A more interesting case arises if the bounds on the coefficients are changed to Then Av7 v can fail to dominate Instead the argument leads to the inequality HVWMMHMM Cwuixhluix C39f 7 w u mu 030 S 030 S and the same optimum order error estimates follow again provided that the HZ regularity holds for the differential problem Lecture 14 Nonsymmetric Differential Operators Let us extend the above analysis in a different direction to cover the not necessarily symmetric problem given by Lu7VaVub1ub2Vucu f XEQ7 aVu blu1 g7 x E 89 for which the corresponding bilinear form on H1Q x H1Q is given by Au7 v aVu b1u7Vv l b2 Vu cu7 v Assume that O lt a ax 3 as before and that lbjxl and assume that there exists a unique solution u E H1Q of Au7 v f7 v l g v v E H1Q7 and that llullz c me Mum Note that we are assuming existence uniqueness and regularity for the differential problem this allows us to drop the imposition of a positive lower bound on c The Galerkin approximation remains being defined as the solution Wh E Mh of AWh7 V f7 V Q Vgt7 v 6 NM Theorem For h sufficiently small there exists a unique solution Wh of the Galerkin equations and Hui Whll0 hllu a Will Clulshs for s E 27 r l 1 when u E H5Q Proof We shall prove the error estimate first the solvability of the equations will be shown afterwards Let us note that lAV7 W S CllVll1 2llWll1 27 V7 W 6 H167 While A is not necessarily coercive on H1Q or Mh it satisfies a Garding inequality on H1Q Av7 v aVv l b1v7 Vv l b2 Vv l cv7 v Z allVVllg CllVVllollVllo 7 CllVllg 1 Z gallVVllg MllVllg 1 Z gallvllfimllvllg this inequality will suffice in place of coercivity Assume for the moment that the Galerkin problem is solvable and let e u 7 Wh then Aev 07 v EMh Since the adjoint problem possesses the same regularity as the original problem there exists 4p 6 H2Q such that Av4p v7e7 v E H1Q Then llellgg New New 7 X S Cllellmllw 7 Xllm so that llello Cllelllh The equality AX 7 Whix 7 Wh 7AU7 xx 7 Wh7 X 6 Mm implies that 1 gallX7 Will 7 MllX7 Whllg S Cl 7 XlllllX7 Whllli so that llX 7 Whllg S Clllu 7 xii Czllx 7 Whllg Now we follow the argument used to establish the convergence of the Galerkin approximation to the solution of a nonsymmetric two point boundary problem Then 2llu 7 xii 2llX 7 Whlli S Call 7 Xlli C4llgtlt 7 Whllg Callu7gtltll 264 Hells llu7xll3 Csllu 7 Xlli Cahzllellg Hell lA lA lA Thus for properly chosen X E Mh 17 Ceh2llell Clul h2 1 so that for h sufficiently small llu 7 Whlll Clulshs 1 and llu 7 Whllo Clulshs7 as was to be shown The existence and uniqueness of a solution Wh can be demonstrated by a very similar argument based first on duality and then on the Garding inequality To show uniqueness assume the data functions vanish giving AWh7 v 07 v E Mh Again there eXists 4p 6 H2Q such that Av4p v7 Wh7 v E H1Q This leads immediately to the inequality llWhllo ChllWhlll Then the Garding inequality shows that 1 1 0 Aiwm mi 2 gallWth a MiiWhii 2 5a a Mm iiWhiii which forces Wh to vanish for sufficiently small h Galerkin Methods for the Dirichlet Problem Consider the Dirichlet problem 7V aVu f X E Q u g7 X E 89 and assume that the problem has a unique solution and that Hullz c Hfllo iigiigm or more generally Hulls cs Hfllsiz iigii ml s 2 0 The regularity assumption imposes conditions on the coefficient a and on the boundary 89 that depend on the choice of 5 To avoid technical difficulties we shall assume these conditions to be met for whatever choices of s that are made below In particular if Q is a bounded convex domain and a E C1Q then the HZ regularity holds Let HglQ v E H1Q v lag g Then the weak form of the Dirichlet problem is to find u 6 H36 such that AuvaVuVv1 v7 v 6 H662 Note that the Dirichlet condition which is an essential boundary condition is imposed strongly on the trial space Hail whereas the Neumann condition which is a natural boundary condition was imposed weakly through a boundary integral without affecting either the trial space or the test space Note also that the test space H662 differs from the trial space for the Dirichlet problem We shall consider several finite element methods for approximating the solution of the Dirichlet problem some based on having the trial and test spaces equal and some where they differ There are advantages to both choices The Nitsche Method The first technique due to Joachim Nitsche 35 is based on having the trial and test spaces equal Let 7 be drawn from a quasiregular family of partitions by triangles an analogous procedure can be defined for partitions by rectangles with the elements being bilinear functions of Q and let h maxdiam T T E Boundary triangles can have one curved edge Let Mh v E C0Q v ire PkT7 TeTh hus 39 f 7 h 7 lt C hkl l XIEFM llu Xllo ll Xlll 7 lulk1 2 7 for any function u E Hk1Q Note that no boundary conditions are imposed on Mh so that it is being considered as a subspace of H1Q not of Him Now multiply the differential equation by a test function v E H1Q and integrate over Q aVuVv 7 Gig vgt f v v e H1Q There are difficulties with this relation First the data function g does not appear second the bilinear form on the left hand side is not symmetric though the problem corresponds to a self adjoint problem Let us address the second difficulty first by symmetrizing the form by subtracting equal terms from both sides to obtain 3WWltagvgtltuaggt fvltgaggt mm Now the bilinear form is symmetric but it is coercive over neither H1Q nor Mh So let us add a second pair of equal terms to get the Nitsche form 8U 8v 0 Au7 v aVu7Vv 7 lt35 vgt 7 ltu7 35gt l E u v 8v va ltg7aggt ltgvgt v E H1Q The form Au7 v is defined over H2Q x H1Q in the sense of distributions however it is not coercive there Fortunately it is coercive over Mh X M for properly chosen a as we shall show Moreover both data functions appear in the right hand side of Nitsche39s method which is defined for g E H a and v E H1Q Now let us show that again for properly chosen a there eXists a unique Nitsche approximation uh the solution in Mh of 8v 0 Am v f v a ltg7agt g ltg vgt7 v 6 M and that U converges to u the solution of the Dirichlet problem Lemma Let the family Mh be quasiregular ie there exist positive constants 39y and 6 such that if T E M and p7 is the diameter of its inscribed ball and hT its diameter then pr 7 gt 6h lt h lt h mhianinhTi39y7 7 77 Then there exists a constant C such that for V E Mh 8v llvlll Ch lllvllo and Cir iivillo 030 Proof Let Tmf be a reference triangle or rectangle for the partition 7 so that quasiregularity implies that T fTre 7 where a Mm B det2 e V142 0 lt V1 lt w lt 00 for all T E 7 Now Mmf MTref C PkTref 0r QkTref for some integer k Since both Hoffef and HLTref are norms on Mmf the finite dimensionality of Mmf implies that there exists a constant C such that iiwii m S Ciiwiioim Since for T 6 7L MhTvf 1o Mref7 it follows that l hZ71VQLl va JacViap This in turn implies that ini2dx iWiZiJacidigi New T Tm V 7 Tm 1 C 2 A C 2 lt 7 dxlt v dx7 W7 Trefiwi Vlhy Ti i and Hviil Ch liiviio for V E Mh Repeating the argument shows that C iVih2 2 S Tivim v E Mm where 1 5 iVihk 2ltE Mir T Recall the proof of the simple trace result Only the elements having a nontrivial edge in the boundary actually are required in the bound on its right hand side so that we can combine the last inequality and the trace result to see that 8v 1 1 5 329 g Ch QlV l lt ClVleV 107 030 7 which demonstrates the second inequality in the lemma Lemma Under the hypotheses on Th and Mh made above there exists a constant 00 such that A is coercive on Mh with respect to the norm 2 l 1 2 2 l 7 v hll 039 lllvlll iv 030 for a 2 00 Moreover A0 is continuous on H x H with respect to the ill norm where 8v Mia H z e H1Q g e Ham 1 Consequently the linear equations of the Nitsche method are uniquely solvable for each a 2 00 Proof If v E Mh the previous lemma implies and the trace inequality for the normal derivative that 8v Aaiviv 2 alVli025H 039 llVll030 EllVllg30 030 2 575 M2 030 h 5 0 30 04 a 7 Cl 3 ElVli0 TllVllg30 8v 2 alvli 7 ch 77 C 2 Z TllVll030 04 8V Zlvlig 6h 2 030 for properly chosen 6 and 6 Then for sufficiently large a the desired dominance follows The continuity with respect to this norm is trivial The existence and uniqueness of a solution uh to the Nitsche equations follows immediately from the coercivity and continuity Let a 2 00 and consider the error e u 7 uh Clearly A0097 V 07 V E Mh As usual let e U7 x l X 7 uh where X E Mh Then AdxiumxiumAAX7wxium Cwui hhwiuMh so that lllxi Uhlll S C llluixlllv X 6 Mh Thus we need to obtain an optimal bound for 8U i X 81 2 HwewrweubhH 1 2 Ell Xlloa 2 030 By the trace theorem and the remarks above MM 81 2 S QU Xhmu Xh2 2 030 and H 7 XHgm S CH 7 XHomU Xha Thus mu 7 W cu e X im Mu e xmue mm hilHuigtltHo 2uix1Q A C U Min hzhl X m h ZHU XH3 27 so that 39 f 7 lt hk X34th U HLQ The Nitsche energy norm error estimate follows immediately k lllui Uhlll S Inf llluixlll lllX Uhlll S Clulkmh XEMh This estimate is in a not quite standard energy norm An L2Q error estimate can be derived by duality Let 7V aVAp e7 x E Q 4p 07 x E 89 the assumed HZ regularity for the Dirichlet problem implies that llwllm S Cllellm Now 8 2 7 7 47 llellog aVein lte7afaygt AaeiwAaeiw7X7 XEMh Thus 2 lt C 39 f 7 lt C h iieiio i Xiiii iieiiOQ v and ii Uhiio 2 S Ciuik1 2hk1 Theorem Let the Dirichlet problem be H2Q regular Let a 2 00 Then there exists a unique Nitsche approximation uh to u and ii Uhiizm hiiiui Uhiii S Ciuir hh provided that u E H Q for some r such that 2 g r k 1 Lecture 15 Interpolation of Boundary Values Let us again consider the Dirichlet problem except that we shall restrict the domain Q to be a convex polygon The convexity of Q implies that so 6 mm ligiisg g s 172 W Assume that 7 is a quasiregular triangulation of Q consistent with the boundary vertices that is if x E 8Q is a vertex at least two triangles in T have x as a vertex Let knvecommTeP nTe Let gh be the piecewiselinear interpolant of g over the boundary vertices of 7 Then 2 2p3 2h7 1S P S 00 opa 2 S Cpllg lwia Let uh be the solution of the Dirichlet problem 7VaVuh 1 X697 uh gh X689 By the maximum principle for elliptic equations llU Uhl 2003072 0oo 2 S Cllgl It follows from the L Q bound that H 7 Uh Logan72 s s 01 s52 S Cllgl Actually more technical results on elliptic equations would lead to requiring less regularity on g to obtain the rates above Set HglQ v6 H1Q vl39g Mgv Mh v laggh The Galerkin approximation to the solution of the Dirichlet problem based on interpolating the boundary values is given by seeking Wh 6 Mi such that AW17 v aVWva f7 v7 V 6 M2 Note that the trial and test spaces now differ since the degrees of freedom for Wh at boundary nodes have been determined by the Dirichlet data It should noted also that Wh is the Galerkin approximation to uh as well as to u and that Wh coincides with uh on 89 Since Av7 v 2 aiv 0 uniqueness and therefore existence is trivial for the approximate problem Then note that Aui uh7 v 07 v 6 H662 Au 7 Wh7 v Auh 7 Wh7 v 07 V 6 M2 Then for X 6 Mg Wh 7X 6M2 uhi Wh 6 H662 and Oliuh WhiQ S AL h Wmuh WhAuh X7uh Wh AU7X7uhiwh iuiXlQuh Whl y so that iuh Whl gciU leb xeM h The piecewise linear interpolants Ju 6 Mi and Juh E M differ by 9h2 in L Q and consequently by 9h in W1gt Q Thus iUh Whim Ciui Juhim S C U 7 Juim iJLI Uhi1 2 S 3 Ciuim iigii2oo3 2h7 an optimal order estimate with respect to rate but not quite with respect to regularity By the Lm bound on u 7 uh iui Whim S C iuizg iigii2oo3 2 h S C iifiio iigii2oo3 2 h Optimality with respect to regularity requires Hgii 30 in place of 2 iigiizpom Let us invoke a duality argument to bound the error in L2Q Let 4p 6 H66 satisfy 7V aVap uh 7 Wh x E Q Then Hapiigg CHuh 7 WhHOQ and iiuh 7 Whiigu AUh 7 Why A h Whiw XL X 6 M2 Thus Huh 7 whiisg BM 7 Wm W W mi gt XEMg so that Huh 7 Whii0 2 Chiuh 7 Wm Ciuim Again Hui WhiiO S C iuim iigiizpom h2 Theorem Let T be a quasiregular partition of the convex polygonal domain Q and let M be the piecewiselinear space over 7 Then 10 S C lLI 29030 hzi if u is the solution of the Dirichlet problem and Wh 6 Mi is the solution of interpolated boundary values Galerkin method 20 llgl 052 hlui Wh llU7Whl Some Extensions 0 The method involving the interpolation of boundary values could be extended to higher order piecewise polynomial spaces but only HZ regularity can be expected of the solution of the Dirichlet problem on a convex polygonal domain 0 If the polygonal domain is not convex the regularity of the solution of the differential problem will reduce to HP where 39y depends on the greatest interior angle at a vertex The proof will go through with the error estimate reduced to llui Whllo hlllui Whlll S Cllulllwhh o If 89 is smooth and convex there is a more interesting extension of the method above using an inscribed polygon There is an analysis for piecewise linear interpolation of the boundary values when a piecewiselinear approximation space is employed and also another more complicated analysis when higher order polynomials are employed along with a boundary condition that transfers the values on 89 to the the polygonal boundary using both the values and the normal derivative of the solution on 89 0 Clearly the proof above holds without essential modification for three dimensional convex polygonal domains These topics will be omitted The Effect of Quadratures It is not always possible to compute the inner products in Galerkin methods exactly or in other cases because they would be computationally overly costly Thus it is common to replace one or more of the inner products in a finite element procedure by quadratures We shall consider some simple examples to indicate the important aspects of employing quadrature rules to evaluate an inner product approximately Let Q be a convex polygon in R2 and consider the homogeneous Dirichlet problem 7VaVu 1 X697 u 07 X689 so that its weak form is to find u 6 H66 such that AuvaVuVv1 v7 v 6 H662 Let 7 be drawn from a quasiregular family of simplicial partitions of Q such that any vertex of Q is a vertex of at least two triangles in 711 Then let Mh C HMS be such that XQLUiV XHo hiivixiil S Civik1hk17 v 6 Hk19 HampQ The Galerkin approximation to u is given by Wh M2V Mh I Vi390 such that AW17 v f7 v7 v 6 M2 Now given a triangle T define a quadrature rule QT on T by selecting a set 9 i 17 7 m of points in T and a collection of weights w gt 07 i1m and setting QTF Z FCiwilTl7 i1 where lTl is the area of T The rule QT is said to exact on PnT the set of polynomials of total degree n on T if QTlPTPXdX7 P PnT If 6 711 it is an affine map of T let U E correspond to C E T by this same map Then define an approximate bilinear form AQ by AQWV Z QTBVWVV TEMh The Q Galerkin approximation ie the Galerkin approximation subject to the quadrature rule QT is given by zh 6 M2 such that AQz17 v f7 v7 v 6 M2 It is also common to apply a quadrature rule to evaluate the right hand side term and the error analysis to be presented below can easily be extended to cover this perturbation of the Galerkin procedure The quadrature rule for the right hand side or to treat lower order terms in the bilinear form normally would be different from the one used for the second order term Let us discuss some examples 1 The Midpoint Rule Let T be the midpoint of the triangle T and set QTF FCTiTi7 AQW7 V Z 3CTVWCT VVCTiTi TQM1 It is easy to see that the midpoint rule is exact on P1T since the midpoint is defined to be that point 7 E T such that Tx 7 CTdx 0 If Mh v E C0Q v 76 P1T7 T 6 711 then AQWV Z aCTVWVvTiTi7 W7V 6M0 T6771 since VW Vv is constant on each T If we maintain the assumption that O lt 04 ax B lt 00 then it is clear that for W and v in Mh AQW7 W Z 64ng and lAQW7 Vl S lW 1 2lV 10 Thus it follows immediately that there eXists a unique solution Wh 6 M2 to these Q Galerkin equations The error analysis for this approximation will be deferred momentarily The proof will be applicable to the piecewise linear procedure on tetrahedra as well as the two dimensional case Consider however the application of this quadrature rule when the trial space is changed to piecewisequadratic functions Mhv C Q vlTE P2T TeTh Here AQW7 W Z aCTlVWCTl2lTl TeTh Coercivity over M2 has been lost that is having VWCT 0 VT 6 7 does not force W to vanish identically An example can be constructed as follows Let Q be the diamond with vertices i10 and 0i1 and divide Q into four triangles by inserting the diagonals Consider the triangle in the first quadrant with vertices at 00 10 and 10 Assume W E PQT x P2T to vanish on the edge X y 1 this leaves free the three parameters associated with the points x1 00 x2 0 and X3 0 Unnormalized basis functions associated with these three nodes are given by 1 whey Xy7 5 Xy17 whey XXy17 whey yXy1 Set W 0449109 5 2X7Y 7309 If 204 7 B 0 then VW7 0 so that the gradient of a nontrivial W vanishes at the center of this triangle It is easy to extend W from this triangle to the other three triangles by symmetry to obtain a nontrivial W 6 M2 for which AQW7 W 0 Thus the Q Galerkin equations are singular and we have established the nonsolvability of the Q Galerkin equations for this choice of the quadrature rule Consequently this quadrature rule fails to be valuable except for the CO piecewise linear case 2 The Midedge Rule Let the midpoints of the edges of the triangle T be denoted by 7 17273 and let QTF wii 3 ZFCTiiTi39 i1 Lemma If F E P2T then QTF T Fxdx Proof t suffices to demonstrate the lemma on the unit right triangle that is to show that T pxdx g ltp07 0 pig p07 i for p E 1X7 yxzxy7 yz Just evaluate the integrals and sums Now consider Mh V E C0Q v 39 O and v ire F 2T7 T 6 AQWV Z QTaVWVv TeTh If again 0 lt 04 ax B lt 00 then it is trivial to see that aiwig WEJVIh7 iwilivil W7 v E Mh AQW7 W 2 iAQW7 M S Thus with this pair MmAQ there is a unique solution Wh E Mh to the Q Galerkin equations AQW17 v f7 v7 v E Mh In practice the inner product on the right hand side would usually be replaced by a quadrature which would usually be a different rule than that applied to A7 We shall not take this into account here Lecture 16 Error Analysis with Quadratu res The analysis of the error associated with quadratures is based on Strang39s Lemma which is useful in the analysis of nonconforming Galerkin methods as well The perturbed bilinear form AQ is restricted only by the hypotheses in the lemma Lemma Assume coercivity and continuity for the bilinear AQ form Let Wh E Mh be the solution of the Q Galerkin equations AQWh7 V f v17 V E Mh39 Let u 6 H66 be the solution of Au7 v f7v7 V 6 H662 Then l Xl1 2 sup W Mb W152 sup lf7X f7Xhl39 XEMq le1 2 7 lt C 39f l Whl1 2 Jerkh Remark Note that quadrature has been allowed in the evaluation of the right hand side We are implicitly assuming that Mh C H662 in the application of this lemma to nonconforming methods the l LQ seminorm usually is replaced by a broken norm of the form 2 1T7 352 Z lX TeTh 2 llX 1h llX and the coercivity and continuity are interpreted with respect to this norm The proof is modified somewhat in the case that Mh g H362 300w rmw X m gt 25 3 w vs S xv H 23 w vs S xv 3 w vs S xv gt25 vs S w x xx X4 E Xv gt E Xv L E Xv L5 vs S xv 3 w vs S xv Mkcwx rw wFwavgtcswav 3 3 xv gt0 5 xv H k w vs S xv 09 3 xv gt00 5 w x E S xv C S x xi Since alwh 7 xli AQWh 7 x Wh 7 x A 7W7 7A 7W7 Q Whxll g mu7leW lWh Xlljl lf7Wh X f7Wh Xhl lWh Xlljl A 7A f 7 f muix m r sup l X7111 QX7 l sup l 711 MMl 1pth W152 1pth W152 from which the lemma easily follows Let us apply the Strang Lemma to our examples The midpoint rule and piecewise linear approximation Let QT denote the midpoint rule and let Mh V E C0Q v 39 O and VT 6 F 1T7 T 6 Th Thus for X71J E Mh AQX7 236 VX Vi TNTL T Now by the BrambleHilbert Lemma and the constancy of VX and Viz on T i3VX7V T 7 QT3VX WM S ChziaVXViMzJ S ChZHaHNHVwaHoJ 2ther W S Ch2a2ooTVXTV T Ti Thus lAlX J 7 AQlX Jll l A Chzllallzoo z lVXT W1TllTl T l A Chzllall2oo 2lxl1 2l l1 2 An analogous argument shows that lAlX J 7 AQX7 l S Chllall1oo 2le1 2l1Jl1 2 Thus A 7A inf uix m sup l xiiJ 00mm Sc u mh39 Xth 1pth W152 If the right hand side is evaluated exactly then it follows from the triangle inequality that l Whl1 2 S Clul2 2h It is also possible to apply the midpoint quadrature rule to the right hand side In that case for V E Mh iifi VT 7 mm T fve amva dx life fltrllorllvllor llfCTllorllv a vltrllor It can be shown by ordinary calculus that llfi fCTll3T S Chzlflio Then llfCTll3T S 2llfll3r2llffCTll3T S C llfllgj hzlflirl Thus by the Poincare inequality lf7 V 5 Mi S Chllfll1 2lvl1 2 Consequently if the midpoint quadrature rule is applied to the right hand side Strang39s Lemma plus another use of the triangle inequality implies that 252 llf lLli Wh 1 9 h 10 S CU Note that f must have point values in order to apply any quadrature rule to the f7 v inner product it suffices to assume that f E H1SQ for some 6 gt O or f E H Q in three space The two error estimates derived above are optimal order estimates in H1Q If the right hand side is evaluated exactly then only minimal regularity of the solution is required However if the midpoint quadrature rule is applied to it then additional regularity is required of the data function f so that f will have point values El It is of interest to compare the exact Galerkin approximation and the Q Galerkin solutions Let Wh be the Q Galerkin solution and let zh be the exact Galerkin solution Az17 v f7 v7 v 6 M2 Let Eh Wh 7 zh so that if the right hand side is evaluated exactly AQ5h7 V A AQXZM V Z aVzh7 VVT 7 QTaVzh Vv7 V 6 M2 T Thus Aoishish 2 vsm Tax a aCTdx T Now the integral of ax 7 3CT vanishes for a E P1T hence S Cllall2ooTh2lTl ax a aCrdx T Consequently aim ChZZ lVEhVZhTllTl Ci hilizhilhz T so that lEhlm S Cthl1 2l72 S Cllflloth If the midpoint rule is applied to the right hand side then an additional term must be added AQ h7 h 2Vzh39vghT 3x aCTldX hh f7Eh T Now W7 hh 125m Tfx 7 mnmxwx 7 and Tfx 7 fCT5hxdx T VfCT x7 T 90x ZT Zsz Ehxdx TVfCT x 7 Ch hxdx Oa szthOth TWltT x7e mm m r xicmdx Of2TH hH0Th2 Twch x7 mm W xichdx oufwuushuwhz ouwmmm T WW oufwzjughumrhz By writing VHCT Vfx 7 sz x 7 T where the last term is to be interpreted as an integral remainder in a Taylor series expansion it follows that T ZlVfCTl2l7lgt 2 S Clflm lflzgh where the norms can be the broken norms related to 7 Then Km 7 mm clifiimishimhz again with a broken norm on 1 Lemma Let 2 and Wh denote the solutions of the exact Galerkin and Q Galerkin equations respectively Then izh a Whim c iifiim minim hz where 9 0 if the right hand sides are computed exactly and t9 1 if the midpoint rule is applied to the right hand side in the Q Galerkin equations An immediate consequence is that llui Wh 1252 72 09 S Cllf 00 0llf An 9h2 estimate in L2Q for u 7 Wh can also be derived by a modified duality argument an example of such arguments will be given for the midedge rule and local quadratic approximation The midedge rule and piecewise quadratic approximation Consider again the choices M2 ve com v lre P2T TeTh v lag 0 3 QTF 2Fyriln i1 where yr is the midpoint of the ith edge of T Again let 2 and Wh denote the exact Galerkin and Q Galerkin solutions respectively where we apply the midedge rule to both sides of the Q Galerkin equations and set 5 Wh 7 zh Also as above Otl hlg S AQEh7 Eh Z BVZm VENT QT3VZh V50 f75hh fill T By the Bramble Hilbert Lemma and the Pg exactness of QT l3VZh7 VENT QT3VZh VEhl S ClaVZh VEhl3Th3 S Cllzhll2Tll5hll2Th37 since zh and Eh are quadratic polynomials on T By the assumed quasiregularity of Th and the Poincare inequality the inverse property llEhllzj S ClEhlLTh l holds Thus lf75hh Ehll39 lEhlm S Clthllh2 2l72 l5hl1 2 Recall that we know that llu 7 Zhllog 7l7 hllu7 zhlllg Cluls hs7 1lt s g 3 For X 6 M2 properly chosen write zh as Zh l X l zh X7 so that 1 1 5 lt2 lthll Jgt S llull2nlt llu7gtltll 1gt T T l 2 ltE lth7gtltll Jgt T Cllullm Ch zlth 7 Xll0 2 Cllullm Ch 2 llu 7 Zhllon llu 7 Xll0 2 Cllullm where quasiregularity has been used several times lA lA lA Next let Jf be the piecewise linear interpolant of f determined by the values of f at yr7 139 127 3 Then QTUEh QTUf Eh and W7 EMT J hhl S lel2Tll hlloTh2 Since by BrambleHilbert and quasiregularity lQTJf39Eh Jf75hTl S ClJffhlsjh3 S CllJfll1Tll hll2Th3 S Cllfll1Tl hl1Th27 it follows that W7 Ehlh f7 Ehll l A C lllfll1Tl hl1T l1 l2Tll hlloTll72 T l A Cllfllh2 2l5hl1 2h2 Combining several of the inequalities above leads to the estimate h2 2h2i 19 Cllu 2 2t9llf th Wh 10 lEh here 9 1 if the right hand side is evaluated using the midedge rule and t9 0 if it is evaluated exactly This inequality implies an optimal order energy estimate lU Wh 27 19 S Cllu 30 0llf h2 2 h for the Q Galerkin approximation it requires minimal regularity if the right hand side is computed exactly However it is not adequate to imply an optimal order estimate in L2Q and an additional argument must be sought Let us extend the duality argument to cover the introduction of quadratures Let 7VaVltp uiw17 XEQ7 4p 07 X689 Then for X 6 M2 iiu7 Whiigg Au 7 WW Au7 w 7 AoWh7 w 7 A 7 AoWh7 w AUW 7 X 7 AoWh7 w 7 X X 7 f7 Xh 7 A 7 AoWw7 AU 7 WW7 7 X 7 A 7 AQWh7X 1 X 7 f7 Xh Observing the derivation of the bound on gh 7 f75hh given above allows us to see that lf7gtlt xhl S Cllfllh2 2llgtltllh2 2h3 Also lA AQWh7Xl S ClWthh3 2h3 S CllWhllh2 2llelh2 2h3 Now there exists X 6 M2 such that ll P Xll1 2 hllelhm S Cllu Whllo 2h Thus as llwhllhgg Cllullzg lluiwhllgn S C llu Whllmh t9llquotllh2 2 73 llull3 h3 lluiwhllonv and ll Whllosl S C llull3 2 t9ll1 llh2 2 h37 which is an optimal order L2Q error estimate A trivial observation is that if ax E 1 and the right hand side is evaluated exactly both components measuring the difference between zh and Wh vanish and Wh zh as should happen since all quadratures are then exact Indeed if f is piecewise constant with respect to 711 the quadrature error due to the right hand side disappears Lecture 17 A Negative Norm Error Estimate Negative norm estimates are estimates in Sobolev spaces dual to spaces of positive index These results will be referenced later in attempting to reduce regularity requirements in the analysis of the errors in Galerkin methods for parabolic equations They can be used in post processing Galerkin approximations to obtain additional accuracy Let 7 E L2Q Then define the H SQ norm of 7 by l77711l max weHsm llw 59 llnllisn E In many applications in the theory of partial differential equations it is advantageous to impose boundary conditions on 7 or 11 or both however the definition above is the most useful one in studying finite element methods Given the Dirichlet problem Lu7VaVuv1uv2Vucu 1 X697 u g7 X689 the solution u is said to satisfy the j shift property or to be j l 2 regular if Hulle c llfllm lab339 for all pairs f g e Him x Warm Other boundary conditions lead to analogous definitions Each such definition requires appropriate smoothness conditions on 89 and the coefficients in the differential equation In particular 89 being in C and the coefficients in C Q plus solvability of the Dirichlet problem imply j l 2 regularity for all j Let us develop negative norm error estimates for the Nitsche method for the Dirichlet problem here we shall assume that all integrals in the procedure are evaluated exactly though it is possible to obtain analogous if more complicated estimates when quadrature rules are taken into account We shall also set v1 v2 O and assume that c 2 0 in addition to assuming whatever smoothness requirements necessary to imply the j l 2 regularity needed in the analysis Let the partition 7 be taken from a quasiregular family and let Mh C H1Q have the approximation property 3 live Xll0 2 hllvixllm civimhh 2 g k r1 h for all v E HkQ Then let a be sufficiently large that the Nitsche procedure determined by 039 8 Aaiwhivaveltgia87gthltg7vgt7 mm is uniquely solvable for Wh E Mh We have already proved that H Whll0 2 hllu Whllm S Clulk hk7 2 S k S r17 if the differential problem has the 0 shift property Now assume that it has thej shift property for somej gt 0 If eh U Wm then Aaeh7 X 07 X E Mh We shall extend the duality argument used to derive the L2Q error estimate in the earlier analysis of the Nitsche method Let 1 E HjQ by the j shift property there exists 4p 6 H1426 H66 such that V 3W A x 6 97 llWlle S Cll llj By application of integration by parts and the vanishing of 4p on 89 em0 em 7V BWVD Aaeh7 w Aaeh7 7 7 X for any X E Mh It is necessary to bound each of the four terms in Aaeh4p 7 x separately The following lemma is needed to treat the last three of these terms its proof can be found in 31 Lemma For sufficiently smooth ax and 89 the following inequalities hold for v and z in H1Q iltvzgti llviiagiizllgaw llVlla 2 S Cllvllmi 8v a S CllVll1 27 V 7530 lMLgau S Cllvllzm Consider A00917 4p i X3 aVemvwix S CHehHl W XHLQv 8 gt 8wix emai S CHehHla 37 7 81 2 8V 7539 89h 8eh 7 7 lt C 7 lt38V7Lp Xgt 7 3131 laQHLp XH w 2 g Ewehwix g Ch 1HehH gH p XH7 Q39 By the lemma above Helm1N CHeHm w m h lexHom Now use the j shift property and the approximation properties of Mh to see that New0i Ciieh mhj iw for O gj r71 where r71 represents the maximum rate of approximability for 4p Thus Hem0i S CiUs 2hH iii1 and 10 Lg 2 s r170 j r717 Hehii Q Cius hs 397 2 sg r17 O gjg r71 The maximum rate of convergence is given by H Whiierm S Ciuirmh2r Under certain circumstances these more rapid rates of convergence in Sobolev spaces of negative indices can be employed in a post processing calculation to obtain a better approximation to the solution of the differential equation than that given by the Galerkin solution Unfortunately these constructions are quite technical in nature and will not be treated here Lecture 18 A Galerkin Method for a Nonlinear Dirichlet Problem This lecture will present an analysis Douglas Dupont 16 of the method of Nitsche 35 as applied to the nonlinear Dirichlet problem 7V ax7 uVu fx7 X E Q ux gm x e am where Q is a bounded domain in IR 7 n g 3 with smooth Le C boundary Also assume that ax7 u is twice continuously differentiable with bounded derivatives through second order and that aXLIZ XlR gt040Oq7 0lta0 a1ltoo Moreover assume that f E C for some 04 6 01 and that g can be extended to Q as a function in C2D It follows from known regularity results Morrey 32 and a uniqueness result Douglas Dupont Serrin 17 that there exists a unique weak solution of the Dirichlet problem and that u E C2D A critical observation used in proving the convergence of the Galerkin approximation is that the formal Fr chet derivative of the nonlinear equation LW 7V ax7 uVW l waLx7 uVu is a linear elliptic operator for which the Dirichlet problem has a unique solution 17 Let WkgtPG denote the Sobolev space of functions with LPG integrable derivatives through order k with norm 1p Hwiimc i Dawipci iDtiSk Let H v e W1gt6 2 8v81 e L289 and for almost all x e 89 there is an open ball BX about X such that v E H2Q For0lthlt1 normey WW lellli Hle h llwlz Maw81427 where lwl llwllom Consider a family Mh0lthlt1 of finitedimensional subspaces of H such that i There exists a positive integer r and a constant 1 such that if 2 g s r 1 and v E H5Q then Xierk hiiiv Xiii hZiiV Xiiw1w6 2 S Clhsiiviis ii For v E C8061 inf iiiiv 7 Xiil 1X 6 Mm suppx CC 9 01 as h a 0 iii For all X E Mh 3 7 7 ah lZiixiil Let us extend the definition of the Nitsche bilinear form to cover the nonlinear differential operator suppressing writing the dependence on X 32 v7 w siziw W a ltag wgt The weak solution u E H of the Dirichlet problem satisfies Bu u W f W 7 ltg agE 7 yh 1wgt w e H As in the linear case the Nitsche form is coercive over Mh with constants depending on Q n7 cl 0407 and 041 for 39y 2 yo PlllVlll2 BZV7V7 VGA1m 2639 moreover it is continuous on H lBzv7Wl czlllvllllllwllly viweni 2689 Fix 39y 2 yo and define an approximate solution uh E Mh of the Dirichlet problem by the equation 8 Biuhuhwfiweltgage h 1wgt wth To see that there exists a solution uh of the Galerkin equations consider the map 5 Mh HMh defined by 8W 3y 5y W 5 W ltgiagg vh 1Wgt7 W 6 Mh The coercivity of B over Mh implies that the linear equations for 5y are solvable Now By y75y S C glll5ylll implies that 5y is contained in a ball in H The Brouwer fixed point theorem implies that there exists at least one solution of the Galerkin equations Turn to the convergence of the approximate solution uh Assume that for some 04 6 01 f 6 Cam and g can be extended to a function in C2D Also assume that u E HSQ for some 5 E 27 r l 1 Finite dimensionality implies that there exists X E Mh such that ll 7 Xll hlllui Xlll wrgmhllllui i lll hllu filillwmm Then Pllluh Xlll2 S 3Uh Uh X7 Uh X 3 Bl U X7 h i X 300 30in WW 7 X S Czlllu Xlll llalu 3uhllL3 2llelW15 2lll h Xlll Since llxllwmm is bounded for our choice of x HW MHW XW q 4l 7 WMmml HW WW S S We will apply the following interpolation relation Adams 1 12 12 0 1 lMWQSCWHHW Then S 1 2 1 2 s lllLHIhlll Clh 1llu7LIhllo llLHIhlll i Clh 1llLHIhlloi Duality will be used to bound llui uhll First note that any Dirichlet problem for the linear elliptic operator L defined above as the Fr chet derivative of the nonlinear elliptic operator has a unique weak solution standard elliptic regularity plus Theorem 2 of 17 Let C ui uh and LW 7V ax7 uVW l auVu VW Let L4ponQ and 4pOon 89 Then llwllz S CllCllo Set 1 sum aux7uxerltx dr7 1 auux 17fauuX7LlX7fCXdf 0 Then llg ll2 Bu u74p 7 Buh uh74p EUCVC75uuC2VU7VL7I llCll2 7 15 3UWC7W7 auUC7Vu W7 as all boundary terms vanish or cancel Continuing HCH2 3u 7 w auUC7 Vu W 3U W 7 3Uh NW 7 30 7 8UhVUh7W7 7auuC7Vu W7 3U W 7 3Uh NW 7 uUCVUh7W7 7auuC7Vu W7 3U W 7 3Uh NW 7 5uUCVC7W7 au i EUXVU th 3U 7 7 3Uh NW 5uuCVC7W7 EuuC2VL 7 V4 3u u w 7 X 7 BM um 7 7 X 5uuCVC7 W7 EWCZVU7 W llCllZ 3U Uwi X BUh um X 5uUCVC7W7 l EUUCZVU7 V4 3Uh 7 w i X EuCVLh V0 7 X EuuCZVLu W7 ClllClll lllwiXlll CllCllW16 2lllClll llwllwmm ClllClll lh llCllL3 2lll ll27 lA lA ll as EuuC39CVLhV P S CllCllL3 2llCll llwllwlm S llCllL3 2llCll lllwlll Hence llCll ClllClll lh llCllL3 2l S CllllClllh lllClllwllClllzl CllllClllhlllClll3l lA lA For h sufficiently small llCll hlllClll S Clhs llCll3l The last inequality would prove the convergence of the Nitsche method for our nonlinear Dirichlet problem if we knew that gt O as h a 0 We will establish this by means of a compactness argument We have shown earlier that lll Uhlll S Clhs 1 llau 3UhlllL3 2l since ax7 v was assumed to be bounded for all X7 v lll Uhlll S C for some C This implies that lth iglloao S Chlz and lluhlll Cllull1 Thus from any sequence hj a 0 we can extract a subsequence hk hk a 0 such that there exists W E H1Q for which uhk a W in L2Q and uhk a W weakly in H1Q Thus Wg7 X689 Our object is to show that W u The first step is to show that W is a weak solution of the Dirichlet problem Let v 6 C306 and let vhk E th be such that each vhk has compact support and live vhklll 01 as ka oo Then 3WVW7VV 5 V 3WVW7VV Vh aWVW 7 auhkVuhk7 Vvhk 7 f7 v 7 vhk ch 7 mm awww 7 awkwthan Now Hv 7 vhk l 7 O and aWVW 7 auhkVuhk7 Vv aWVW 7 uhkVv aW 7 auhkVuhkVv 7 07 by the strong convergence in L2Q of uhk to W and its weak convergence in H1Q along with the boundedness of Vuhk Thus W is a weak solution in H1Q of the Dirichlet problem Let us assume temporarily that the weak solution is unique so that W u Consequently uh a u in L2Q as h a 0 Thus a O as h a O and ll Uhll hlllui hill S Chg under the uniqueness assumption To show that weak solutions of the nonlinear Dirichlet problem are unique i weak solutions are smooth solutions ii smooth solutions under the assumptions on f and g are in C2D Q iii such solutions are unique i Since W is an H1Q solution of the Dirichlet problem the famous results of di Giorgi and Nash see Morrey 32 Theorem 537 show that W is Holder continuous on 5 ii By Morrey Theorem 554 W E W1gtPQ for 1 lt p lt 00 Then note that W is a weak solution of l 3W w g X 6 897 AW am ValtpVLp 1 X697 which can be shown by testing the equation against 1 aWv for v E C8062 This implies that AW 6 L2Q from which it follows that W E H2Q Finally Schauder estimates Morrey Theorem 563 imply that W E CHOKE iii Theorem 1 of DouglasDupont Serrin 17 shows that the C2 solution of the Dirichlet problem is unique which completes the proof the error estimate derived above Under slightly more restrictive conditions on the subspaces Mh one can show that a Newton iteration for uh converges for sufficiently small h see 17 Lecture 19 Galerkin Methods for Parabolic Problems A Model Heat Conduction Problem We shall begin with the simplest problem associated with parabolic equations namely the initial boundary value problem for the heat equation with variable coefficients cgivlaVu f XEQt OT7 u 07 x689 t 0T7 u L107 x69t07 where 0ltoz cx ltoo7 0ltoz ax ltoo7 ffxt7 u0 u0x E L2Q There are several reasonable weak forms of the heat equation the first choice is between a global weak form obtainable by integration against a test function over Q x 07 T or one that is local in time resulting from testing either at each 1 E 07 T or over 9 x t 7 At 1 Each choice is useful but under different circumstances The simplest of these is obtained by testing for each value of the time and that will be what we consider Let us define the weak solution of the differential equation to be the continuous map u 0T H66 such that 6 v l aVu7Vv f7 v7 v 6 H662 t E OT7 u u0 t0 This weak form can be approximated either continuously in time or discretely at times t nAt where At gt 0 We shall consider both possibilities Continuous in time Galerkin approximation Let 7 be taken from a quasiregular family of partitions of Q and let Mh C H66 have the approximation property for 1 lt s g k l 1 that o 2hllVX 39f 7 ltc XIEth lv X 152k lquot 5375 w e Mum439 i e Mhv com v lre PkT TeTh vlag 0 over simplices or the tensor product equivalent over rectangles or cubes We are implicitly assuming Q to be polygonal if it is not then it would be necessary to employ the Nitsche bilinear form or boundary interpolation methods to account for the Dirichlet boundary condition The extension of the arguments below to their use is straight forward The continuous in time Galerkin approximation of the model parabolic problem is given by the continuous map Wh 07 T aMh such that 8 lt Wh v l aVWva f7 v7 v E Mh t E OT7 Ci 8t Wh W2 6 Mh t 07 where W approximates uo There are many choices for W2 but usually in practical problems u0 is reasonably simple and can be approximated easily by interpolation or projection We shall employ an elliptic projection 3VW27VV aVu07Vv7 v 6 NM This choice would seem to impose more regularity on the initial condition than belonging to L2Q but this is rarely a practical constraint Before taking up the convergence of the Galerkin approximation to u let us discuss some computational aspects of the continous in time procedure Clearly the method generates a system of ordinary differential equations first order in time for the degrees of freedom of Mh Consider the following example Let Q 01 Cx ax1 h N 1 and Xj jh MhV C0Q v lTE F 1Tj7 73g1g v0vN0 With the usual hat function basis Ag7 j 17 7 N7 1 for Mh the Galerkin procedure in this case is equivalent to the system d 1 2 1 W391 7 2W l W391 1 E EMJ21 EWj 6Wj1gt j 1 2 1 x 56710 36t 56w J 17N7 1 if Simpson39s rule is applied to the right hand side or equivalently f is piecewiselinearly interpolated Note that the matrix multiplying the derivative terms is not in diagonal form In a single space variable this poses no serious difficulty however in two or more space variables it is likely that any ODE package will be forced repeatedly to solve what amounts to an elliptic difference equation in order to advance in time Thus the continuousin time Galerkin procedure is of more theoretical than practical interest since the ODE solver is unlikely to have efficient elliptic difference equation solvers built in The convergence analysis for the Galerkin procedure contains almost all of the important concepts of the analyses for the discretetime procedures we shall introduce below and is technically simpler hence it is advantageous to present it Assume that as stated above Mh is quasiregular and that the Dirichlet problem 7VaV1J g7 X E Q 1 07 X E 89 is HZ regular ie llwl 20 S Cllgl 00 The error analysis is indirect and is based on using the information already obtained concerning the convergence of Galerkin methods for elliptic equations Wheeler 45 Introduce the elliptic projection of u7 t as follows zh 0T AM where aVzh7Vv aVu7Vv7 v E Mh t E 07 T Assume that 8 u 7 0 T L2 Q vat 7 1e Then if elliptic projection equation is differenced in time and At sent to zero it will follow that there eXists 82h E OT th Also 82h 7 8H BVE7VVgt 7 lt3V 7vvgtv V E M t E 07 Moreover from the known error estimates for elliptic Galerkin methods iiuizhiio i hiiuizhiil S Ciuis hsi 2 SSS k17 au 82h 8U 82h 8H 77 777 C 7 1 0 8139 i A h 2 r k1 r 2 8t 81 0Q i8t 8t Nowset ehWh Zh Mh7 772h U Since aVrLVv O for v E Mh it follows that 88h CE v aVeh7Vv H A G J lt v lt m 3 Choose v eh Note that 68 ili e e 8t7h 2dt h7h7 so that integrating with respect to time gives gosh eht 0taVehVeh7d7 70tlt57ehgt7d1 By the bounds on the coefficients iiehiio 2dT 0 2 t 2 t Eiiehtii00a iehi1 2dT 0 0 0 09 t 8 2 1 20 iiehiim at 2 dT 00 87 8t i A 2 dT 00 In particular t an iiehtii30 c Meni n HE 0 The continuous version of the discrete Gronwall Lemma is needed to bound llehllg Lemma Let F and G be nonnegative monotone increasing differentiable functions such that t Ft g Ct M FTd7 0 for some M gt 0 Then if Ft Ct eMtFm M eMt TGTd7 0 Proof It is clear that F is maximized when equality holds in the inequality thus assume that Ft Ct MOt FTd7 Differentiate this equation and multiply the result by e Mt to see that Hwy G The solution of this differential equation can be written in the form Ft Ct eMtFO 7 eMtGO MOt eMtTgtGTdT which implies the lemma Now set t an FrHehH3Q7 6106 a 2 d7 and M C 00 Then by the lemma above 15 HehHSQm c t c 0 t C 0 t 2CeCt 0 2 t 739 dT CZ eqt 00 0 0 2 2 t 2 t CtiT d7 C d5 e d7 0 0 at 00 6 Q 0 8t 8t 87 2 E 005d5dT Td7 Q 8777 2 8t Td7 00 Apply the bound for llQnatllo 15 052E Cl 0 This error estimate can give an optimal order convergence rate but the regularity required of the solution is not minimal 8U llehf a 2 Td7 gt hr7 2 r kl r 2 By utilizing the negative norm error estimate 877 8t 7 we can reduce the required regularity This bound is valid for k gt 1 if Q and the coefficients are such that H3 regularity holds that is the elliptic projection zh satisfies the shift result hr1 r 2 2 r k17 lg 10 8t 30 S Cllg lth 10 For H3 regularity to hold 89 cannot be polygonal consequently the Galerkin method under consideration should be based on say the Nitsche bilinear form with the Dirichlet boundary condition imposed through integrals rather than strongly The negative norm inequality is invalid for piecewise linear approximation and as stated earlier is associated with piecewisequadratic or higher order approximation Since optimal order error estimates were proved for elliptic equations with lower order terms these terms could have been added to the spatial operator in the parabolic equation without changing the convergence analysis except to require that h be sufficiently small which is irrelevant in an asymptotic argument Lecture 20 Two Level Time Discretization We shall treat the model problem again for simplicity since the analysis will be based on the elliptic projection lower order spatial terms can be included without essential change in the analysis for small h However practically it is necessary to distinguish between problems for which the solution is dominated by diffusion the second order elliptic term and those for which transport or convection the first order spatial terms dominate The treatment just below is aimed at diffusive domination transport dominated problem will be treated later Also we will simplify by taking vanishing boundary values which can be treated by the use of a Nitsche method as usual The timediscretization problem was treated in considerable detail in Douglas Dupont 14 the analysis preceded the introduction of the elliptic projection approach of Wheeler and could be simplified with newer techniques The Differential Problem X7 t 7 V ax7 tVu fx7 X E Q u 07 X689 The Discrete Problem Again let 7 be taken from a quasiregular family of partitions of Q and let M0 C H66 have the approximation property for 1 lt s g k l 1 that inf NV 7 Xllo th 7 xllm Clvls hs w e HSQmH3Q XEME Assume Q to be polygonal if it is not then it would be necessary to employ the Nitsche bilinear form or boundary interpolation methods to account for the Dirichlet boundary condition The two level or 0 Galerkin approximation of the model problem is given by W tn AM Whquot WWquot 1 such that Wn 17Wn WWW aquot9VWhquot9Vv Wu veM I Wh 7 w eMO t 07 where gnle 0g 1 l 17 t9g 7 12 9 g 1 and where W Hg and uh is the elliptic projection of u given by aVWh7Vv aVu7Vv7 v E 107 t E 07 T If 0 g 9 lt12 a CFL constraint would be required for stabilty There are other equally reasonable ways to discretize the coefficients a and c i EMU 0aquot1x 1 t9quot007 ii a 9x aixi Wt iii aVWhquot9 oaquot1xVWhquot1 17 0aquotxVWhquot and i Cquot9X 0Cquot1X 1 t9CquotX7 ii cn9x Cx We 7 The analysis to follow covers all cases related to these choices with minor modifications since n9un1 7 un n9 n9 n9 n 0 c Tm a Vu 7Vvf l 6v7 VGAt7 for all choices where 6 x o 0 7 gm Arm Thus with 7 Hg 7 u and choices for the interpolations n1 7 n lt6fn97 v 051 ng 17 manwg Vv n17 n Ctquot9 7 W9 6quot v v 6 M2 At Set 9 Hg 7 W Since 77 1 7 n At 9h5 if u E L O7 T HSQ n1 7 n lt6fn97 v 0a 1Ve 1 17 t9a Ve 7 Vv lt0 hs 0 7 gm At v v 6 M0 Let v e 1 and note that n1 7 6In tn9 Q n1 d 7At 7e 1 gt7 2At Ctn9en17 6ln1 7 Ctn9en7 en 7 oaquot1Vequot1 17 0a Ve Ve 1 2 530 7 1aquot1 7 39yAtVe 1Ve 1gt 7 7 t9a Ve 7 We Thus it follows after multiplication by 2At that Ctn9en17 6ln1 7 Ctn9en7 en 7 30 71 1 7 2yAtVe 1Ve 1 7 17 93 Ve Ve At c hs a 7 gm 7 my Henliim 2 g C h5 0 7 At At Ar 39yllenllllet Add on time after some manipulation it follows that a 7 ClAflle 1ll2 7 62 llekllzAt k1 30 71a 7 2yAtVe 1Ve 1 At 7 ma 7 239yAtVek Wk At k1 2 g C h5 0 7 At At At k1 l l quot 1 hs a 7 gm my 2 t Since 9 2 there exists h such that for h lt h n n71 2 1 llenll2k 71llVekll2At c1 El lleklletJng hs a 7 Em m At for t g T An application of the discrete Gronwall lemma implies that l 2 TAt 1 n k 2 s 2 Inglis ii k llVe ii At clth 07 AtAt ln equivalent notation lth WhllzwoTL2 2 lth Whllz2oTH1 2 g c h5 a 7 gm my First the estimate for U 7 Wh in H1Q is a superconvergence result as we can imply only that ll WhlllltgtltgtoTL2 2 hll Whlll0TH1Q C hs 0 7 gm my from the estimate above for U 7 Wh and that for u 7 uh The error bounds for u 7 U are of optimal order but a more detailed derivation would show that they are not optimal with respect to regularity For a smooth boundary and with piecewisepolynomials of degree at least two optimal regularity can be obtained using the Nitsche method for the elliptic part of the differential operator Let us derive another error estimate that gives a uniform error bound for the flux aVu Denote the timedifference of a function z by dtz zr391 7 z At Then with e uh 7 Wh and 7 uh 7 u as before choose the test function dte and note that CHMquot 6 aquot on a 7 gm my The auxiliary error equation becomes cr39iedter397 dte 0a 1Ve 1 i 17 93 Ve dte an7 dte Now cr39Jredter397 dte 2 oz dter39HZ Q E oan1Ven1 17 0a Ve oan1Ven1 17 0a 1Vequot Ve At gquot 01 Thus Q We Q7 dtVe amp0a 1Ve 1Ve 1 17 20aquot1Vequot1 V9 717 0a 1VenVequot 0 dtVe At Since M17 20a 1Ven1Veni 0 7 a 1Ve 1Ve 1a 1Ve Ve 2oziidtenii2At a 1Ve 1Ve 1ia Ve Ve 2 CiiVen1ii2 HVeniiz At c hs a 7 gm my At Adding on time the Gronwall lemma and the known estimate for HeH give the estimate iidteii22onL2 2 iieiizwamnwm C h5 07 gm At2 Again the bound for Heniil is superconvergent The full estimate above implies that ii Whiil0TL 2 hii Whiil 0TH1 2 g C h 07 At At2 which is optimal with respect to rate Time Stepping Galerkin Methods for Nonlinear Parabolic Problems Consider the nonlinear parabolic problem given by X7 Mg 7 V ax7 uVu fx7 t7 u7 X f gX7 f ux7 0 u x7 X697 X689 X697 for 0 g t g T The simplest discrete Galerkin method for this problem is based on backward differencing the timederivative term and using the Nitsche bilinear form Bz v7 W azVvVW 7 ltag8V7 for the elliptic part of the operator Wgt Thus the discrete problem requires the solution of the nonlinear equations CWEHMWL V BWE1WE17V 1 Wh 17v7ltg 17ag 1 7 8V 39yh 1vgt v E Mh at each time step ie a full nonlinear elliptic Galerkin problem at each step All that can be expected from this time discretization is that the error is 9At plus that due to spatial discretization The object of this discussion is to avoid having to solve nonlinear equations to advance in time So what can we do to linearize the method above First suggestion Replace W ll by W in the evaluation of the functions Cu7 au and fu so that the algebraic problem reduces to the linear problem CW dtW 7 V BWWIZ3917V awn v a ltgn17agquotl a w v e M The boundary values g could also be evaluated at time t The linearized procedure can be proved to converge in optimal order Second Suggestion For some physical problems it is highly desirable to have the coefficients evaluated at the same time level as the approximate solution Clearly this should be accomplishable by iterating using updated evaluations of the coefficients but this is expensive Consider instead a predictor corrector approach This amounts to updating the coefficients just once so that the elliptic problem at time t 1 is solved twice The difference between the predictor corrector solution and that of the fully nonlinear problem can be shown to be OAt2 Third Suggestion At times t 7 n 2 2 extrapolate linearly in the evaluation of the coefficients n3 7 n n71 Wh 7 2Wh 7 Wh for constant At This leads to the linear problem CW dtW 7 V Blwhn i Willllv V Wm v a ltgn17agquot1 7 8V 39yh lvgt v E Mh Again one can expect that the difference between the solution of these equations and that of the fully nonlinear problem to be OAt2 Linearization of the Crank Nicolson Method The fully nonlinear Crank Nicolson method is given by the equations 1 l l CWhnlg tw v BM WE v l awful v a ltgnialtgn ig a wlvgt7 v 6 M In order to retain the second order accuracy in At a predictor corrector corrector procedure can be used Alternatively a four level extrapolation can be applied in the evaluation of the functions a7 c7 and f This linearization is analyzed in Douglas Dupont Ewing 15 Lecture 2 1 Mixed Finite Element Methods A general reference for mixed methods is the book of Brezzi Fortin 8 The object here is to analyze the mixed method of Raviart Thomas 37 or Raviart Thomas N delec 34 for second order elliptic equations on bounded domains in R2 Several other methods will be discussed briefly Let Lu 7V aXVu bxu l cxu fx7 X E Q u gx7 X689 where we assume the problem to be solvable uniquely for g 6 89 x H32ao with M2 S Qlllfllo lglszl Denote the flux by v 7aVu l bu7 and let Rewrite the elliptic equation as a first order system ozv l Vu l uO7 Vvcuf Note that many if not most mathematical models of physical processes have such a first order system as their natural descriptions with the second order models being derived by combining the equations of the first order models Spaces v Hdiv 9 v 6 39 V v 6 39 W 39 Weak form of the Dirichlet problem Find v7 u E V x W such that avyCliluvv ClerMCltg7C Vgt7 C V7 Vv7 W l cu7 W f7 W7 W E W These equations do not represent a minimization problem instead they lead to a saddle point problem Mixed Finite Element Spaces We need finite element spaces V x W approximating V x W They cannot be chosen independently as then the resulting linear equations will almost always be singular A thorough discussion of the requirements for a pair of function spaces to be applicable in the formulation of saddle point problems both differential and finite element is given in Brezzi Fortin The primary requirement is an inf sup condition to be satisfied between the spaces The treatment in these lectures will be based on a commuting diagram property which implies the inf sup condition The RTN spaces will be shown to have the commuting diagram property The Commuting Diagram Property Given a pair Vh x Wh of finitedimensional spaces such that Vh C V and Wh C W Vh x Wh is said to have the commuting diagram property if there exist projections I39ll1 V HVh and Ph W a Wh with Ph being the L2Q projection such that div I39I1 Phdiv V a onto Wh Defining the I39Ih projection for a prospective mixed finite element space is essentially the determining factor in validating a MFEM procedure The Raviart Thomas N delec Spaces On triangles or tetrahedra Let PkT PAT 7 n 2 or 37 and x X177Xn Set VAT SpanPkT7xF kT7 WhT PkT7 and Vk v E Hdiv 9 v 76 VkT7 VT 6 7L Wh W 6 29 2 W iTE WhT7 VT 6 Boundary triangles triangles or tetrahedra are allowed to have one curvilinear face the spaces VkT and WkT are unaltered The consistency condition for v E Vk across element faces is that the normal component of the vector be continuous across the face The most popular of the RTN spaces is RTNO for which the vector space has dimension n l 1 on each T The useful degrees of freedom on T 6 7 are the constant normal components of the vector vh on T and the constant value of the scalar Wh on T Thus the projection Pk can be defined by requiring VI7nhVId507 1 j n17 5 where is the jth face of T Thus the consistency condition on T1 T2 is that VT1 V1 VT2 I2 07 v EVk Lemma For k 2 0 dimVkTk1k37 n27 dimVkT k1k 2k 4 n 3 Proof For n 2 The dimension of P2 is k 12 and since only the part of kaT involving homogeneous polynomials of degree k add to the basis of VkT 2k l 1 new basis functions complete the basis of VkT Lemma If v E VAT then V v E PAT7 VI FE RAF7 VF C 8T Moreover V VAT a PAT is surjective Proof Clearly V v E PAT and V xPAT PAT Now write v p17 p2 p3x where pj E PAT Then VVp17p2vp3XI7 xeF Since x I is constant on the face F v 11 E PAF as was to be shown Convenient degrees of freedom for VkT are determined by the following uniqueness lemma Lemma For k 2 O and v E VkT v 0 if ltVV7pgt 07 p E PkI 7 VF C 8T V7QT 07 q E Pk71Tquot n Proof For k 0 the proof of the lemma is trivial since the dimension of VkT is equal to the number of faces edges and P1T is empty For k gt 0 v 11 O on 8T Now VPkT Pk1T so that the second constraint in the lemma is equivalent to 0V7VqVv7qltVI7qgt7 qe PAT HenceVv0andv0 The I39Ih projection can be defined locally as the map I39ll1 VT a VkT such that v7 I39Ihv 1 pgtF 07 p E PkI 7 VF C 8T v7 I39Ihv7qT 07 q E Pk1T For k 2 0 dimVkTk1k37 n27 dimVkT k1k 2k 4 n 3 For n 2 the dimension of P2 is k 12 and since only the part of kaT involving homogeneous polynomials of degree k add to the basis of VkT 2k l 1 new basis functions complete the basis of VkT Since on T1 T2 I39Ihv 11 l I39Ihv 12 0 I39I1 is eXtendible globally A modification of its local definition is necessary on the curved portion of a boundary elem ent one keeps the constraints on the internal flat faces of the element and requires that v 7 I39Ihv be orthogonal to the complement of VkT with respect to those vectors having vanishing normal components on the internal faces Lemma div I39I1 Ph div maps V onto Wh Proof It follows from the proof of the uniqueness of the degrees of freedom that V v 7 I39Ihv L PkT and this implies the lemma Since v I39Ihv for v 6 Pk the Bramble Hilbert lemma implies that live I39Ihvllo Cllvllrh 7 1g r k 1 It also follows from duality with respect to LZ projection that W v7 Miles CHV viirh39n o k1 llwiPhwllssg Cllwllh 7 0 r75 k1 Redundantly V V nhv7W07 WE Whi VV7W7PhW07 VEVh Lecture 22 Let V Vh7 UnhV Vh7 7u7u17 739F huiuh7 pu7Phu Then a6 v i V 39 V777 5717 v 07 v 6 Vh V 57WCU7W07 WeWm or equivalently O 7v7v v7T T7v 7 pvv7 VEVIM V 7 Wc7397 W 7co7 W7 WE Wh The remainder of the convergence analysis follows Douglas Roberts 21 Duality Lemmas These lemmas will allow us to devise an argument analogous to that of Schatz 41 for noncoercive bilinear forms in Galerkin methods Define negative norms by duality with respect to HSQ llwllisllwllis2n sup W 6H5 2 lli lls As earlier we say that Q is s l 2 regular if the Dirichlet problem L 7V3V bV c 6L2Qxen 4p 07 X689 is uniquely solvable with llwlls2 S Cllwllsv 1 6 H56 If V is the dual space of V and f E V fv f0v i f1V v where f0 6 L2Q2 and f1 6 L2Q Lemma Let k 2 land 0 g s g k7 1 Assume that Q is 5 2 reguar Let C E V f f0 f1 6 V and g E W L2Q Let z 6 Wk satisfy 04Cv7Vvz zvfv VEVk V C W i 62 W g W W 6 Wk Then for sufficiently small h iiziiis S Qihstliiciio hSHiiV Ciio iifii7571 hsiliifiio iif1iiis hsii iio iigiiisiz hsiziigiioi Proof Let 11 E HSQ and let 4p 6 HS2Q be the solution of the dual Dirichlet problem Since AV aVapD z7 39haV4p 27 27 7V 3W7 3 3W7 W f39hanaC739han 0827 an han 527 Since v7 W 7 Phw O for W E Wh 274 62 V 39 Cw 7 PM 7 PM i V 39 C74 and V39Ciw CW7QC78W7 a9 FlhaW7 6 ani FlhaW7 Substituting 2711 fquot39haW7 g PW 04C Z EVA7 7 39haV V C 627 7 lamp Now f haW7 S WE VW WEWF haW7 S HfoHisAHaVWHH H LsHVfanms HfoHoHaV MEWVHO Hf1HoHV 39 EWF MBWVDHO S QWoHisil hSHHfoHo WM th Ho W s Similarly gy PM E QUMLH h52HgH0 W5 Also Ma zvawi haW l W C c2477 PW Q thHUlCllo llzllo WWW He llzllol llwlls Consequently lzlis ths llCllo WW He llfoll7571 ll llis h1llfollo hsll llo llgllisiz Wiigiio hs iiziioi For sufficiently small h the hsllllzllo term on the right hand side can be absorbed in the left hand side and the conclusion of the lemma holds The cases k O and s k when k 2 1 are not covered by the last lemma However reading the proof and noting the limit on approximability for 4p 7 Phtp leads to the following lemma Lemma Let k 2 O and let Q be k 2 regular Let C E V f f07 O E V and g E L2Q If 2 E Wh satisfies the equation of the previous lemma then for h sufficiently small llzllik S thkHUlCllo llVCllo llfollo llgllo llfollik71 llgllikizl Since the projections I39I1 and Ph were defined locally the bounds for 2 derived above can be localized Error Estimates Recall the error equations O 7v7v v7T T7v 7 pvv7 VEVIM Va WCT7 W CP7 W7 W6 Wm where v7v17 03939hv7v17 7u7u17 739F huiuh7 pu7Phu For small h and Q 2 regular the last two lemmas imply that We S thll llo h2 5 llV Elle llpll71 hllpllol LZ projection p u 7 Phu HpH71 thHo QHuHWE 0 r k1 So 739 Phu 7 uh HTHo S QWEHO hZ MHV 5H0 HUHrh39HJ Then 7 u 7 uh W0 S W0 M0 S QWEHO hHmHV EH0 MW Witho 39hv7vh and v7vh by Vv739hvw07 WE Wh VUWV 7W7 WE Wh By the second error equation and the choice W V 039 W 0H0 S QHUHm and W Ho S QUMHO HV VthqL 0 S q S k1 Recall again the first error equation a iv V VJ 577V PNL V E Vhi Take the test function to be v 039 After a bit of manipulation Mi 0 V OJ aV hVL 0 7 5777 0 so that llUllo Qlllnllo liviitht 1 t ML and ll llo llvi Vhllo S Qlllnllo llVlltht 1S f S k1 Combining several of these estimates gives llnllo S thllnllo llullrh39 llVlltht1 llV 39 Vllqhq2 6 7 forOgrgkll7 O q kl and 1 t kl For small h and r t 1 q 2 7 6K0 since llVllrel llV 39 Vllrez S Qllullri lelt WW kq Qllullh 7 iszland2 r k1 Orfork0and0 q 2 Mbnenmmmmwz2eqmux It follows that ll llo llviwllo Qllullr1th 17 q r k17 llV llo llV vivmlo Qiiuiirmz o g r k 1 This completes the analysis of the error in L2Q for the RTN mixed method for the Dirichlet problem The analysis above was made under the assumption that there existed a unique solution to the RTN mixed method It suffices to establish uniqueness Let f O and g 0 Take W V vk in the second equation defining the mixed finite element method V vhV vh cumv vh 07 so that llV39Vhllo S QlluhllO The duality lemmas imply that lthllo S ththllo for small h Then choose vh as the test function in the first equation in the mixed method thNh 7 V Vhi Uh uhivh 07 so that ithiio S Qiiuhiio S Qhiivhiim and vh O and then uh O for small h Thus the RTN mixed finite element method possesses a unique solution Lecture 23 Error Bounds in H SQ The following error bounds are proved in DouglasRoberts 21 under the assumption that Q is s i 2 regular HuiLIhHis Qiiuiih O s k7172 r k17 Qiiuiir1hrk7 5k7 1 r k17 Qiiuii2h k1 sk10 r kl Hviwiiis lt Qiiuii1h 5 0 s k 1 r kl 7 Qiiuiir2hrk17 5 k17 0 S r S k 1 iiVVVkiis Qiiuii2h 0 s k170 r k1 Rectangular RTN Mixed Elements Let x 6 IR and let Pk1k2mknx Pk1X1 Pk2X2 PknXn Spanxflxggxfn 0 S j S kjv j 17n Then on a rectangle R in IR set Pk1kX Pkk1 2 RTN R gt gt 7 7 kll Pk1kk X Pkk1k X Pkkk17 n 7 l 9quot If T R is a partition of a domain Q by rectangles let V E I V lRE E T Qm W E 129 i W lRE Pk1kgmknR7 kj k7 VR 6 T It is clear that dimRT Vk1R 3k1 Hg and that V RTNkR QlkR If e is a face of R and v E RTNkR then e W 27 It is convenient to define WAR Pk71kR X Pkk71R7 397 27 Pk71kkR X Pkk71kR X Pkkk71R7 n 3 Lemma If v E RTNkR then 7 Vap E Pke7 n 27 ltv VWIgt 7 07 l W7 E Qkei7 397 37 V7 R 07 W E WAR imply that v 0 The proof of this lemma is straight forward since the boundary integrals for example reduce the first component to a polynomial in Pk1kk and the integrals over the interior kill what is left The lemma above indicates convenient degrees of freedom on R as well as leading to an obvious definition of the desired I39Ihm projection The analysis of the convergence of the mixed finite element method based on RTNM x QM is essentially unchanged from that for the triangle based RTN method Let us consider some other mixed finite element spaces The BDM Spaces 5 BDM Spaces Over Triangles n2 BDMlt lies between RTNk1 and RTNk has smaller dimension than RTNk and provides the same order of asymptotic convergence for the vector variable as does RTNk If a hybridization procedure to be defined later is applied to the BDMlt method a superconvergence phenomenon related to the scalar variable can be used to produce an approximation to the scalar variable having the same asymptotic rate of convergence as that coming from the RTNlt method Vector Space VkT BDMkT PRUk Scalar Space WAT Pk71T V BDMkT WkT and dimBDMkT k 1k 2 Let RkT 6 L28T gt ee PkeVe 6 8T Then v1 E RAT7 Ve 6 8T Lemma For k 2 land v E BDMkT Vp e 07 p E Pke7 Ve 6 8T V7VqT 07 q E Pk71T7 v TO7 eBDMkT vVgtVOvz3TO7 imply that v 0 Proof Note that V vVvT 7VVVv7 l ltv1Vvgt3T7 so that the first two sets of relations imply that V v 0 and the third set of relations implies that v 0 As before the lemma leads to a definition of a I39ll1 projection on BDMlt v E Hdiv Q v lTE BDMkT7 VT 6 T The commuting diagram property holds with this projection so that an error analysis can be made along the same lines as presented for the RTN mixed method The following theorem is proved in Brezzi Douglas Marini The argument uses the first duality lemma derived in the analysis of the RTN MFEM and the rest of the proof is a slightly simpler version of the RTN proof Theorem Let V17 uh E BDMlt x Wk be the solution of the MFEM problem 04Vh74P V 4P7Uh ltg7 vgt7 veBDMm Vvhw w7 WEWk Then Hui uhiiis Oliiflriz igiHzih ti 2 s r k2 0 g s ke2 llV Vhllis Qlllfllril lglr12lh i 1 r ML 0 g s k 17 llVv7vhllS Qllfllh s Ogrgk 0 s k Hybridization of the BDM MFEM Solving the linear equations generated by a mixed method or any other saddle point problem is usually a more expensive computation whether directly or iteratively than for a minimization problem Fraeijs de Veubeke 28 29 proposed a procedure that leads to a minimization problem by relaxing the requirement that the vectors have continuous normal components across element faces and meet this requirement by introducing Lagrange multipliers for the scalar variable on these faces Considerably later Arnold and Brezzi 3 showed that these multipliers can be used in many cases to obtain a more accurate approximation to the scalar variable than that which comes directly from the mixed method This hybridization will be developed for the BDM MFEM Let e denote the set of all edges of triangles T E T Then let MkmmleePkeifeCQandmle0ifeC89 Modify the definition of BDMlt to eliminate the constraints on the normal components across internal edges Vk v v lTE BDMkT7 VT 6 T7 and keep unchanged the choice of Wk Wk W I W lTE 3 i1T7 VTET Lemma If v E Vk then v E Vk if and only if ZltVV77mgt37 07 m6 Mk T Consider the simple Dirichlet problem in mixed form avVu 07 XEQ7 Vv 1 X697 u g7 X689 Now test the first equation above over T E T against 4p O VvhpT 7 V LPvuT ltV VT7UgtT 039 The approximate solution uh is discontinuous across element edges Fraeijs de Veubeke replaced uh in the boundary integrals by Lagrange multipliers and added over elements Thus the hybridized BDM MFEM is defined as finding Vh7 Uh m 6 Vk x Wk x Mk such that 0mm72V39W7UhTZltLP VT7mhgtaT lt4P39V7ggt7 4P Vk7 T T ZVvhwr w7 WEWk7 T EmlynOT 07 PEAIo T where the third equation forces vh to belong to Vk as desired Remark 1 The solution vh7 uh considered as functions of the hybridized problem satisfies the original BDM MFEM simply restrict the test function LP to Vk Thus the portion of the solution of the hybridized problem representing the vector and scalar variables is unique The Lagrange multipliers are also unique as will be shown below Remark 2 While the vector solutions as functions coincide their parametric representations are different dimVk gt dimVk Consider the algebraic problem generated by the hybridized method its matricial form is given by Av Dtu Etm g7 DV 1 Ev The matrix A is now block diagonal by T E T and readily invertible thus it is practical to eliminate v v A lg 7 A lDtu 7 A lEtm Substituting this relation into the Dv f results in the equation DA lDtu DA lEtm 1quot which can be solved for u since DA lDt is again block diagonal this leads to the equation 5m1 EA lDtDA lDtXEA lDtYm fquot where 5 is symmetric positive definite and sparse Thus m can be evaluated at more or less the same cost as for finding the solution of a Galerkin finite element method Then u and v can be evaluated by back solving Following Arnold Brezzi let us make use of the Lagrange multipliers mh to define a function uz that converges to the solution u of the differential problem more rapidly then does uh Lemma Let k 2 2 be an even integer Then there eXists u be Pk1T determined by the relations lt I mh7Pgte 07 P E Pkei7 172737 LIZ Um qT 07 q E Pk72T This extension is not known for odd integers except for k 1 and k 3 For k 1 let u l7 be determined by ltuzimh1gte 0 123 Uziumph 07 p P1T For k 3 let u 76 P4 T satisfy ltle 7 mh7qgte 07 q E P20577 i1237 LIZ Uh7PT 07 P E P2T7 where mh g if e C 89 Then 3 Uuolt QHUW k17 h QHuHmhw k2346 There is an alternate way to make use of the Lagrange multipliers to obtain a better approximation of u define a Nitsche approximation locally on each T E T Let T E T and define H E Pk1T as the solution of 8u AaT 7P Vu VP 7 ltaailfmgt 3T 8p lt 7 35gt3T Uh71ltmh7PgtaT 3 7 fpeltmha8igt ahT1ltmhpgtaT pe Plum V 3T An argument slightly varying from that for the convergence of the Nltsche method shows that For k 2 2 Hui uhio Q llfllk igiMz W2 Rectangular BDM Elements Let k 1 Let T where R is a rectangle in R2 such that diam R hR h and the ratio of the lengths of its sides is bounded by a constant independent of R and h Set WkR lquotz1 7 Wk W 6 29 2 W lRE WAR7 R E Using VkR PkR fails to provide a stable pair of spaces for the mixed method however augmenting this choice by exactly two polynomials of degree k l 1 suffices Let VkR PkR 63Spancurlxkllycurlxyk17 Vk v E Hdiv Q v lRE VkR7 R E T As usual define the Ph projection to be the L2Q projection into Wk For a normal rectangle define the 39hR projection into VkR by the relations p E Pke 12 34 v E Pk2R ltVinhvIe7pgte 07 Vinhvile 07 If R has a curved edge say e4 modify the relations for 39hR to ltV nhV Vs Pgte 07 P E Pkei7 17273 V Vi hvli WR 07 W 6 Pk71R vil39lhvvR 0 v y VkRVy0 andyIROon e1Ue2Ue3 The same error estimates hold for the BDM MFEM based on rectangular elements as for one based on triangles Indeed they are valid as well for a MFEM based on mixing triangular and rectangular elements with triangular elements around the boundary of Q and rectangular elements in the interior It is also advantageous to hybridize the method with rectangular or a mixture of elements The Lagrange multiplier space on an edge e is Pke The multipliers can be used to produce an improved approximation of u either by the local Nitsche approach or through a rather complicated analogue of the first method treated for triangular elements see The BDDF Elements for Elliptic Problems in R3 Reference Brezzi Douglas Dura n Fortin 6 These mixed finite elements for second order elliptic equations on domains in R3 lead to the same asymptotic rates of convergence as do the corresponding RTN elements and are defined by considerably fewer degrees of freedom than those required by the RTN elements Simplicial Elements VkT PAT7 Vh V E HdV Q I V lTE VkT7 VT7 WkT 3 z1T7 Wh W 6 29 2 W lTE WkT7 VT Cubic Elements Let ZkR Span curl007inlycurl0xzk107 curlyk12OOcurl00xyi12k 17 curl0x l1ykilzOcurlx lt 1yzi1007 i1k VAR PAR 692W Vh v E Hdiv Q v lRE VkR7 VR7 Wh W 6 39 W he Pk1R VR The analysis of the BDDF MFEM is essentially unchanged from that of the BDM MFEM since the critical projection I39Ih is constructed in It is perhaps even more advantageous to hybridize this MFEM than the BDM MFEM since three dimensional problems are frequently expensive to evaluate In addition there is an alternating direction iteration discussed in The BDFM Rectangular Mixed Elements Reference Brezzi Douglas Fortin Marini 7 The bases for the BDM and BDDF mixed finite element spaces while of lower dimension than the bases for the RTN spaces do not have the intuitive feel of rectangular relationship The BDFM spaces in both two and three dimensions offer bases with more intuitively rectangular character than those for the BDM and BDDF spaces while remaining more efficient than the RTN spaces Planar Elements BDFMkR v E BDMkR VVR E F k1e7 i12347 BDFMh v E HdivQ v lRE BDFMkR7 VR7 Wh W 6 29 2 W lRE lquoti1 7 VR Cubic Elements The cubic BDFM elements are formed by restrictions of the BDDF elements Let X X17 X1X17 X3 and denote by F khom7 i the set of homogeneous polynomials of degree k in x BDFMkR Pk Pkme7 1 gtlt Pk Pkme7 2 gtlt Pk Pkme7 37 BDFMh v e HdivQ v lRe BDFMkR VR Wh W 6 39 W lRe Pk1R VR The analysis of the BDFM MFEM is very little different from that of the previously defined MFEM39s Again hybridization is valuable and another alternating direction iteration procedure is discussed in 7 Various Other Properties and Extensions Related to BDFM Elements Superconvergence Let Q be a rectangle in R2 Then it can be shown that convergence at a rate hk in L occurs for the vector variable on horizontal and vertical lines through certain Gauss quadrature points in the rectangular elements partitioning Q whereas convergence in LPQ7 1 g p g 07 cannot be faster than hk1 See 24 Also see 33 44 for related results for RTN rectangular spaces Finite Volume Methods The simplest finite volume method or cell centered finite difference method is derivable from the lowest index RTN space by the choice of a quadrature rule for the vector vector inner product The BDFM mixed finite element spaces of index 2 in either two or three dimensions can also be used as a basis for deriving finite volume methods again all that is necessary is to choose quadrature rules that diagonalize the matrix corresponding to the vector vector inner product In the planar case the following quadrature rule h2 7 4 Rfxy dxdy z E lt8fxR Z fxRgt7 i1 where xR is the center point of R and xR is the ith vertex of R For the cubic BDFM element the appropriate quadrature rule is unusual in that it involves not only point values of the integrand but also the two dimensional Laplacian at the midpoints of the faces of the element See 10 where the BDFM based finite volume procedures are treated along with several others based on mixed elements In particular the asymptotic accuracy of the BDFM based methods for both the scalar and the vector variables is preserved In all of the finite volume methods derived in 10 the volume element coincides with the mixed element R E T Other finite volume methods frequently involve two partitions of Q often with one being vertex centered with respect to the other Prismatic Mixed Finite Elements Reference Chen Douglas 11 Partition Q into prisms of the form K T x 272 hz7 where T is a triangle in the X7y 7 plane Three families of mixed elements are developed in 11 with the RTN BDDF and BDFM mixed spaces serving as sources for the three families All are analyzed in the citation with the proofs being variants of those in 21 An Iterative Technique for Solving the MFEM Equations References Brezzi Fortin Fortin Glowinski 27 and Saad 40 Write the MFEM in matricial form AvDtBu f DvCug7 where A relates to the avhv inner product B to the v7 u inner product C to the cuh7 W inner product and D to the V vh W inner product Then note that A is invertible Thus v A lf 7A 1Dt A lBu7 and Mu DA lDt c DA lBu BA lf 7g DA lDt is positive semi definite DA lDth A lth m2 2 vllthllZ so that DA 1 l C is positive definite However it is not necessarily the case that DA lDt l C l DA lB is invertible Consider the special case in which B 0 Theoretically the equation for u could be solved directly but A 1 is almost inevitably a full matrix rather than being a sparse matrix so that it is usually impractical to proceed directly If a well chosen quadrature rule is applied to the avhv inner product for some MFEM39s to be defined later the matrix A will reduce to a diagonal matrix and then the matrix M remains sparse and the direct method which is then a finite volume method 10 is effective and efficient The Uzawa Iterative Algorithm Limit the discussion to the case B O and C 0 The Uzawa Algorithm 8 40 1 Choose Ugo E Wh 2 Forj 2 0 solve for vg17u j1 Aw 7 Drug Mom um wDvhj1 g 3 Iterate until convergence is obtained for v597 HEP If DA lDt is symmetric positive definite then the Uzawa algorithm converges if and only if 2 0ltw lt maXDAilDt The Uzawa algorithm is more or less a time discretization of the algebro differential system AvDtu 7 1 8H iip 7 7 8t quot g If the first equation is tested against v and the second against u then 1d MW Aviv in v a g u If one asks to bound the error between the iterates and the exact solution of the MFEM equations it suffices to take f O and g 0 with initial errors v0 and um Clearly the continuous in time problem has a decaying solution The constraint on w reflects both the explicit treatment of u and the difference in time levels for v and u in the v equation The Arrow Hurwicz Iterative Algorithm 40 1 Choose v07 um 2 Forj 2 0 compute Vj17 uj1 V041 W g f 7 AVG Drum 7 Lj1 HO w g Dvmn 3 Iterate until convergence is obtained for v597 HEP Intuitively the Arrow Hurwicz algorithm corresponds to timestepping the parabolic system gAvDtu 1 8H E FDV g7 for which Lecture 24 V VV V V A Convection Dominated Parabolic Equation Let Q C foo7 00 7 n 2 or 3 Let fs 2 O and set Fs f HE d5 Assume that j gtX and u ux7 t are specified and that Vuq07 XEQ7 X E 89 where q 0 corresponds to an incompressible fluid in many applications u 11 0 corresponds to no transport across the boundary uIO7 The model equation in nondivergence form a gt8 fsuV57Vst 0 In divergence form V Fsu 7 st Fsq 0 gt Initial condition 5X7 0 50X gt O gt To obtain a no flow boundary condition for s Fsu 7 st 11 7st 11 0 An Evolution of Eulerian Lagrangian Methods The Modified Method of Characteristics References 22 38 39 gt The MMOC is based on the nondivergence form of the equation gt Rewrite the transport part of the equation as a directional derivative 8 8 s s i f i at 5quot V5 wow where w w snunzri T a muquot gt Thus gt It suffices to consider discretization in the time variable to present the main concepts of the MMOC gt The differential equation 85 11 7 V st 7 0 gt Let At gt 0 t nAt gt g gt Consider the following fractional step method 1 Set S J so 50X 2 For n 2 1 solve for transport fltuVltO ertquot 1lttgtquot ltx t Squot 1x x e 0 u1O7 X6190 and set gquot 091 7 X E Q 3 Solve for diffusive effects 7Vdvn0 x69 t 1ltt t 71X7 f quotX7 X 6 Q7 an 11 07 X E 89 and set X 77X7 t 7 X E Q 4quot If t lt T loop back to step 2 Discretization of the fractional step procedure gt To discretize the transport part of the time step we have two choices 1 Approximate 64787 0 2 Approximate 1564781 and f u V separately in the equation 6lt6t fltu Vlt 0 2C 2 gt In trans ort dominated roblems 7 is small and 7 and p p M at 82C w are usually both large somewhere We discretize along characteristics gt The characteristic through X7 1 is given by y y t fltu gt7 yX t X Approximate this characteristic on tn 17tn by its approximate tangent at X7 1 and define its foot at time tnil by Y Y X X 7 f 1xu xAt gt gt Complete the transport step by approximating gquot by quotX s WOO gt Next discretize the diffusive step in time by backward differencing n 9397 At 7VdV O7 XEQ7 dV I 07 X689 gt This completes a time step gt Advantages to the MMOC 1 The transport step is explicit and is very fast 2 Any reasonable discretization of the diffusive step leads to a symmetric matrix in the induced linear equations instead of a nonsymmetric matrix which would have resulted from choosing to discretize the time and spatial derivatives separately gt The Disadvantage to the MMOC 3 It can fail to preserve as algebraic identities important conservation principles associated with many differential problems The Modified Method of Characteristics with Adjusted Advection Reference 18 19 gt The conservation error in the MMOC is entirely in the transport step which being explicit requires only a minor part of the total computation gt To adjust the advection execute the following algorithm 1 Let Q 1 Q 15974 clx 2 Compute Y and Equot as in the MMOC 7 fS 1Xu XAt 7 Y X X 7 6 X X x squot1ynxi 3 Compute 6 6 dx Q and compare an and Qn l If they are equal within a very small tolerance accept gquot as g and proceed to the diffusive step 6 below If not it never isl go to 4 4quot Let 6quotX v6 XAf7 v31 gtix roo ng maxSn 1 gt 1 gt 7 ifanltQ 1 min Sn 1 gt 1 gt if 6 gt QH an n dx 2 5 Find 9 such that Qnil enan 7 t9n n7 and set A i N X 0 nx i 17 0 X g defines the adjusted advection 6 Do the diffusive step 5quot 7 Equot At gt 7V dV dV u X697 0 07 X E 89 gt Clearly Q gt dxQ gt dxQ gtS 1dx and conservation now holds for the approximate solution gt A practical note all of the above was presented with the time steps for the transport and for the diffusion equal in practice we normally microstep the transport using 10 15 transport microsteps per diffusive step gt In a waterflood problem the difference in computing time between the MMOC and the MMOCAA was less than 3 thus the second objective of not significantly increasing the computing time while accomplishing global conservation was achieved A remarkable increase in accuracy was also obtained As a consequence it appears that equivalent accuracy using the MMOC would require at least 10 times as much computing time as the MMOCAA for the immiscible displacement problem when the physical parameters are realistically modelled as fractals A Locally Conservative Eulerian Lagrangian Method References 2 20 23 gt The LCELM is based on the divergence form of the differential equation rewritten as gts Vt 7VXstFsq XEQ7 tgt07 Fsu stIuIO7 X613Qt gt07 sx7 0 50X7 X E Q 1 gt Assume that u ux7 t is known and is such that vuqqx7r x697 uIO7 X689 gt The differential fractional step procedure for the LCELM Vt lt PM gt FCq x 6 Q t 1 lt t g t 0quot uI7 X68971 rquot1ltt t 7 CX7 f 1 5X7 f 17 X 6 Q7 X CX7 t 7 X E Q 80 gt57VdV00 x69 t 1ltt t dve 11 0 x e 89 t 1 lt t g t 0X7tn 1 E X7 X E Q 5X7 t 9X7 t 7 X E Q gt The associated Conservation Principle 1 Let n Ufil ClosureK with K n Kj ab 2 Define W C Q as follows Let X 6 8K and let X gt X 2X tn l where f 7 Fltu t gt t 2 t H dti lt 2X t X X 2X tn l E 8 3 Let Di be the tube domain determined by having K as top W the interior of 6 as bottom and the surface Ei generated by the integral curves 2X t as lateral surface tnjuI n K ypw K X Figure The SpaceTime Domain D 4 Integrate the transport equation over Di Since 154 VD i Fu 0 Einv 7V FZngxdt K CndxiV7 CnildxEn Fiat VDndg gt ndxi gtsn71dx FCqudt K K D Thus gt ndx gtsn71dx FCqudt K K D gt The LCELM is based on using this conservation relation to determine the approximation E of the transport gt Fundamental questions include the choice of the sets K and the algorithm for finding the corresponding Kn sets gt After experiments with several choices of K0 we settled upon gt Let us detail a Mixed Finite Element Method on the simple domain Q 01 partitioned into squares of side length h M0 X1X x ykhyj X ih7 yj jh for the differential system in mixed form v i st 07 v Fsu v Fisiq by applying the RTO mixed finite element space gt To define K n gt u uX t is assumed to be specified gt The approximation 5 of s is discontinuous on 6K0 Use the bilinear interpolar S 1 of Sn 1 related to cellcenter values 53771 in defining Kg gt Let the vertices of K be Xijm H 1 i i i 4 and let F Enil N n KER Xijm Xumu le At H 1 i i i 4 Xij 1xu gt Let V be the quadrilateral having the four vertices Y0 gt Given K0 and consider the discrete transport fractional step Theoretically we should have dx gt 1dx FCqu dt KU KU Dn Practically 1 En a constant on Kg 2 The integral of 1594 over K can be evaluated exactly 3 Approximate the integral over D by the trapezoidal rule D n Fltquf KU F qndx73 FSn71qn1dx 4quot Thus the transport fractional step is given by h2 0 g A gt 1dx KU At in 7 aspK U qndx I 1q 1dx7 75 OI 7 At 7 2 n n n h g07 7FijKU q dx gt 1dxg F 1q 1dx V 2 V 5 The right hand side above can be evaluated exactly as can fKU qndx Hence the equation is a simple nonlinear equation for solvable element by element This completes the transport step Usually for most Kij q O and g3 is found directly gt Consider now the discretization of the diffusive step given by v dV0 XEQ7 n7 n71 Vv 0 XEQ7 v IO7 X689 0n 1x E X7 X E Q gt Apply the RTo space over the partition M MU ie seek Vr397 5 E RT0M VM x such that d 1V 7yV39y75 07 y 6 WM st gtT7W VV 7W 07 W6 WM7 V 11 07 X E 89 gt This completes the definition of the diffusive fractional step and consequently that of the full time step gt It is possible and in fact desirable to employ the approximate flux to evaluate the gradient of the scalar approximation and to post process the approximate scalar to be linear on each element using this gradient before initiating the next time step This adds almost nothing to the computational requirements and does lead to a more accurate approximate solution gt A convergence proof for the LCELM when it includes the post processing has been found by Douglas Spagnuolo and Yi 23 gt For a linear parabolic equation Arbogast and Wheeler defined a special case of the LCELM and both proved convergence for it and applied it to a miscible displacement model for flow in porous media An Application to a Waterflood Problem 300 A00 500 Figure Comparison of MMOC MMOCAA and LCELM Simulations top to bottom at 200 days using a 256 X 64 grid for a small CV 49 500 Figure Comparison of MMOC MMOCAA and LCELM Simulations top to bottom at 200 days using a 256 X 64 grid for a large CV References R A ADAMS Sobolev Spaces Academic Press New York 1975 C ALMEIDA J DOUGLAS F PEREIRA AND L C ROMAN Algorithmic aspects ofa locally conservative eulerian lagrangian method for transport dominated diffusive systems In quotFluid Flow and Transport in Porous Media Mathematical and Numerical Treatment CONM Book Series Z Chen and R Ewing eds American Mathematical Society 2002 v 295 p 37 48 D N ARNOLD AND F BREZZI Mixed and nonconforming finite element methods implementation postprocessing and error estimates RAIRO Mod lisation Math matique et Analyse Num rique 19 1985 7 32 S C BRENNER AND L R SCOTT The Mathematical Theory of Finite Element Methods 3rd edition Springer Verlag New York 2008 F BREZZI J DOUGLAS JR AND L D MARINI Two families of mixed finite elements for second order elliptic problems Numer Math 47 1985 217 235 F BREZZI J DOUGLAS JR R DURAN AND M FORTIN Mixed finite elements for second order elliptic problems in three variables Numer Math 51 1987 237 250 F BREZZI J DOUGLAS JR M FORTIN AND L D MARINI Efficient rectangular mixed finite elements in two and three space variables RAIRO Mod lisation Math matique et Analyse Num rique 21 1987 581 604 F BREZZI AND M FORTIN Mixed and Hybrid Finite Element Methods Springer Verlag New York 1991 W L BRIGGS V E HENSON AND S F MCCORMICK A Multigrid Tutorial 2nd edition SIAM Philadelphia 2000 Z CA1 J DOUGLAS AND M PARK Development and analysis of higher order finite volume methods over rectangles NEW for elliptic equations Advances in Computational Mathematics 19 2003 3 33 Z CHEN AND J DOUGLAS JR Prismatic mixed finite elements for second order elliptic problems Calcolo 26 1989 135148 P G CIARLET The Finite Element Method for Elliptic Problems North Holland New York 1978 J W DEMMEL Applied Numerical Linear Algebra SIAM Philadelphia 1997 J DOUGLAS JR AND T DUPONT Galerkin methods for parabolic problems SIAM J Numer Anal 7 1970 575 626 J DOUGLAS JR T DUPONT AND R E EWING Incomplete iteration for time stepping a Galerkin method for a quasilinear parabolic problem Slam J Numer Anal 16 1979 503522 J DOUGLAS JR AND T DUPONT A Galerkin method for a nonlinear Dirichlet problem Math Comp 29 1975 689 696 J DOUGLAS JR T DUPONT AND J SERRIN Uniqueness and comparison theorems for nonlinear elliptic equations in divergence form Arch Rational Mech Anal 42 1971 157168 J DOUGLAS JR F FURTADO AND F PEREIRA On the numerical simulation of waterflooding of heterogeneous petroleum reservoirs Computational Geosciences 1 1997 155 190 J DOUGLAS JR C S HUANG AND F PEREIRA The modified method of characteristics With adjusted advection Numerische Mathematik 83 1999 353 369 J DOUGLAS F PEREIRA AND L M YEH A locally conservative Eulerian Lagrangian method and its application to nonlinear transport in porous media Computational Geosciences 4 2000 1 40 J DOUGLAS JR AND J E ROBERTS Global estimates for mixed methods for second order elliptic equations Math Comp 44 1985 39 52 J DOUGLAS JR AND T F RUSSELL Numerical methods for convection dominated diffusion problems based on combining the method of characteristics With finite element or finite difference procedures SIAM J Numer Anal 19 1982 871 885 J DOUGLAS JR A M SPAGNUOLO AND S Y Y1 The Convergence ofa Multidimensional Locally Conservative Eulerian Lagrangian Finite Element Method for a Semiinear Parabolic Equation To appear J DOUGLAS7 JR7 AND J WANG Superconvergence of mixed finite element methods on rectangular domains Calcolo 26 1989 121133 R FALK AND J OSBORN Error estimates for mixed methods RARO Anal Numer 4 1980 249 277 M FORTIN An analysis of the convergence of mixed finite element methods RAIRO Anal Numer 11 1977 341354 M FORTIN AND R GLOWINSKI Augmented Lagrangian Methods North Holland Amsterdam 1983 B X FRAEIJS DE VEUBEKE Displacement and equilibrium models in the finite element method Stress Analysis 0 C Zienkiewicz and G Hollister eds Wiley New York 1965 B X FRAEIJS DE VEUBEKE Stress function approach International Congress on the Finite Element Method in Structural Mechanics Bournemouth 1975 C JOHNSON Numerical Solution of Partial Differential Equations by the Finite Element Method Cambridge University Press Cambridge 1987 J L LIONS AND E MAGENES Non Homogeneoous Boundary Value Problems and Applications l Springer Verlag New York 1972 C B MORREY7 JR Multiple Integrals in the Calculus of Variations Die Grundlehren der math Vissenschaften Band 130 Springer Verlag New York 1966 M NAKATA A WEISER AND M F WHEELER Some superconvergence results for mixed finite element methods for elliptic problems on rectangular domains In The Mathematics of Finite Elements and Applications V MAFELAP 1985 J Whiteman ed 1985 J C NEDELEC Mixed finite elements in R3 Numer Math 35 1980 315 341 J A NITSCHE Uber ein Variationsprinzip zur Losung von Dirichlet Problemen von Teilraumen Randbedingungen unterworfen sind Abh Math Sem Univ Hamburg 36 197071 915 A QUARTERONI AND A VALLI Numerical Approximation of Partial Differential Equations Springer Verlag Berlin 1994 P A RAVIART AND J M THOMAS A mixed finite element method for second order elliptic problems In Mathematical Aspects of the Finite Element Method Lecture Notes in Mathematics volume 606 pages 292 315 Springer Verlag Berlin New York 1977 l Galligani and E Magenes eds T F RUSSELL An incompletely iterated characteristic finite element method for a miscible displacement problem PhD thesis University of Chicago Chicago 1980 T F RUSSELL Time stepping along characteristics With incomplete iteration for a Galerkin approximation of miscible displacement in porous media SIAM J Numer Anal 22 1985 970 1013 Y SAAD Iterative Methods for Sparse Linear Systems 2nd edition SIAM Philadelphia 2003 A H SCHATZ An observation concerning Ritz Galerkin methods With indefinite bilinear forms Math Comp 28 1974 959962 J C STRIKWERDA Finite Difference Schemes and Partial Differential Equations Wadsworth amp Brooks Pacific Grove California 1989 R V VARGA Matrix Iterative Analysis 2nd edition Springer Verlag Berlin 2000 J WANG Asymptotic expansions and L error estimates for mixed finite element methods for second order elliptic problems Thesis University of Chicago 1988 M F WHEELER A priori L2 error estimates for Galerkin approximation to parabolic partial differential equations 1973 723 759SIAM J Numer Anal 10

### BOOM! Enjoy Your Free Notes!

We've added these Notes to your profile, click here to view them now.

### You're already Subscribed!

Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'

## Why people love StudySoup

#### "There's no way I would have passed my Organic Chemistry class this semester without the notes and study guides I got from StudySoup."

#### "I made $350 in just two days after posting my first study guide."

#### "Knowing I can count on the Elite Notetaker in my class allows me to focus on what the professor is saying instead of just scribbling notes the whole time and falling behind."

#### "Their 'Elite Notetakers' are making over $1,200/month in sales by creating high quality content that helps their classmates in a time of need."

### Refund Policy

#### STUDYSOUP CANCELLATION POLICY

All subscriptions to StudySoup are paid in full at the time of subscribing. To change your credit card information or to cancel your subscription, go to "Edit Settings". All credit card information will be available there. If you should decide to cancel your subscription, it will continue to be valid until the next payment period, as all payments for the current period were made in advance. For special circumstances, please email support@studysoup.com

#### STUDYSOUP REFUND POLICY

StudySoup has more than 1 million course-specific study resources to help students study smarter. If you’re having trouble finding what you’re looking for, our customer support team can help you find what you need! Feel free to contact them here: support@studysoup.com

Recurring Subscriptions: If you have canceled your recurring subscription on the day of renewal and have not downloaded any documents, you may request a refund by submitting an email to support@studysoup.com

Satisfaction Guarantee: If you’re not satisfied with your subscription, you can contact us for further help. Contact must be made within 3 business days of your subscription purchase and your refund request will be subject for review.

Please Note: Refunds can never be provided more than 30 days after the initial purchase date regardless of your activity on the site.