Biomedical Engineering - Mathematical and Numerical Methods in Engineering

Completed notes of the course - Numerical Methods

Complete course

Numerical Methods Course held by Prof. Christian VergaraYear 2020-2021 Notes by Alessandro M. IppolitiLeonhard Euler 1707-1783. 2 Contents I Finite Dierences 5 1 General Overview 7 2 ODE Approximations 92.1 Approximations of the 1st Derivative . . . . . . . . . . . . . . 9 2.1.1 Forward Euler Method . . . . . . . . . . . . . . . . . . 9 2.1.2 Backward Euler Method . . . . . . . . . . . . . . . . . 10 2.1.3 Central Approximation Formula . . . . . . . . . . . . . 10 2.1.4 Forward and Backward 2nd Order Formulas . . . . . . 12 2.1.5 Approximations of the 2nd Derivative . . . . . . . . . . 12 2.2 Solving the Cauchy Problem . . . . . . . . . . . . . . . . . . . 12 2.3 Model Evaluation Parameters . . . . . . . . . . . . . . . . . . 14 2.3.1 Consistency . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3.2 Zero Stability . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.3 Convergence . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.4 Absolute Stability . . . . . . . . . . . . . . . . . . . . . 17 2.4 The Theta Method . . . . . . . . . . . . . . . . . . . . . . . . 19 2.5 The Advection Problem . . . . . . . . . . . . . . . . . . . . . 19 2.6 Solution through the Forward Euler Centered Method . . . . . 20 2.6.1 Consistency . . . . . . . . . . . . . . . . . . . . . . . . 21 2.6.2 Zero Stability . . . . . . . . . . . . . . . . . . . . . . . 22 2.6.3 Convergence . . . . . . . . . . . . . . . . . . . . . . . . 23 2.6.4 Absolute Stability . . . . . . . . . . . . . . . . . . . . . 23 2.7 The Upwind Method . . . . . . . . . . . . . . . . . . . . . . . 242.7.1 Consistency . . . . . . . . . . . . . . . . . . . . . . . . 24 2.7.2 Absolute Stability . . . . . . . . . . . . . . . . . . . . . 25 2.8 Numerical Viscosity . . . . . . . . . . . . . . . . . . . . . . . . 26 2.9 Backward Euler Central Method . . . . . . . . . . . . . . . . . 292.9.1 Absolute Stability . . . . . . . . . . . . . . . . . . . . . 29 3 4 CONTENTS 3 The Heat Equation 313.1 The Heat Equation . . . . . . . . . . . . . . . . . . . . . . . . 31 3.2 Numerical Approximation . . . . . . . . . . . . . . . . . . . . 323.2.1 Discretization in Space . . . . . . . . . . . . . . . . . . 32 3.2.2 Discretization in Time through FECM . . . . . . . . . 33 3.2.3 BECM . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.3 The Theta Method . . . . . . . . . . . . . . . . . . . . . . . . 343.3.1 Consistency and Zero Stability . . . . . . . . . . . . . . 35 3.3.2 Absolute Stability . . . . . . . . . . . . . . . . . . . . . 35 3.4 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.5 Ghost Node Approximation . . . . . . . . . . . . . . . . . . . 42 3.6 Algebraic Analysis . . . . . . . . . . . . . . . . . . . . . . . . 44 4 Advection with Diusion 474.1 Numerical Discretization . . . . . . . . . . . . . . . . . . . . . 48 4.2 Possible Stabilizations . . . . . . . . . . . . . . . . . . . . . . 514.2.1 Upwind Method . . . . . . . . . . . . . . . . . . . . . . 51 4.2.2 Consistency . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2.3 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.3 Algebraic Problem . . . . . . . . . . . . . . . . . . . . . . . . 52 5 The Wave Equation 555.1 Second-Order Time Discretization . . . . . . . . . . . . . . . . 55 5.2 The Leap Frog Method . . . . . . . . . . . . . . . . . . . . . . 56 5.3 The Newmark Method . . . . . . . . . . . . . . . . . . . . . . 57 5.4 The Wave Equation . . . . . . . . . . . . . . . . . . . . . . . . 585.4.1 Leap-Frog Method applied to the Wave Equation . . . 59 5.4.2 The Newmark Method applied to the Wave Equation . 60 5.4.3 First-Order Methods applied to the Wave Equation . . 61 5.4.4 Elasto-Dynamic Problems . . . . . . . . . . . . . . . . 62 5.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 II Finite Elements 65 6 Elements of Functional Analysis 676.1 Elements of Functional Analysis . . . . . . . . . . . . . . . . . 67 7 The Strong and Weak Formulation 73 CONTENTS 5 8 The Galerkin Approximation 858.1 Cea's Lemma . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 8.2 Solving the Galerkin Approximation . . . . . . . . . . . . . . 88 9 Finite Elements in 1D 939.0.1 Linear Finite Elements . . . . . . . . . . . . . . . . . . 94 9.0.2 Construction of the Linear System . . . . . . . . . . . 95 9.0.3 The Reference Interval . . . . . . . . . . . . . . . . . . 97 9.1 Quadratic Finite Elements . . . . . . . . . . . . . . . . . . . . 1029.1.1 Construction of the Linear System . . . . . . . . . . . 103 10 Finite Elements in 2D 10710.0.1 Linear Finite Elements . . . . . . . . . . . . . . . . . . 108 10.1 Accuracy of Finite Element Methods . . . . . . . . . . . . . . 109 10.2 Increasing r or h? . . . . . . . . . . . . . . . . . . . . . . . . . 110 6 CONTENTS Part I Finite Dierences 7 Chapter 1 General Overview Linear Equations The most general form of aPartial Dierential Equation(PDE) is, @ u(x)@ t r2 u(x) +
ru(x) +u(x) =f(x) In a physical system, each term describes a dierent property regarding the behavior of the functionu(x; t) as it evolves: ˆTheInertial Term@ u@ t describes the temporal evolutionof the sys- tem. ˆTheDiusion/Viscous Termr2 u, with0 describes the evo- lution of a process thatdiuses without a preferential direction of propagation. ˆTheAdvection Term
ru, with= ( 1: 2; :::;  n) models a process that propagates along thedirectiongiven by the vector. ˆTheReaction Termudescribes alocal reaction process. For  >0 the system presents sources, while for 08v6 = 0 and (vjv) = 0()v= 0 72 CHAPTER 6. ELEMENTS OF FUNCTIONAL ANALYSIS Convergence and CompletenessThe ability to dene a norm also allows us to introduce the notion ofconvergenceandcompleteness. A sequence vn2 Vis said to beconvergent in the normkk if, lim n!1k vv nk = 0 for a suitablev2V. A sequencev n2 Vis said to be aCauchy Sequencein the normkk if, lim n;m!inf tyk v n v mk = 0 This second notion is in general weaker than that of convergence, as we could imagine an oscillating sequence whose elements converge around a value but which would not have a denite limit. Finally, if the spaceVis such that all Cauchy sequences in the normkk are convergent with respect to the same norm, then the space is said to becomplete. Hilbert SpacesLooking at the properties of the scalar product we can see that a scalar product always induces a norm. The functional in the spaceV dened as, kvk= (vjv)12 can be shown to satisfy all conditions pertaining to a norm. If a spaceV iscomplete with respect to the norm induced by its scalar product, thenVis said to be anHilbert Space. A particularly useful example of an Hilbert Spaces is the space ofsquare-integrablefunctionsL2 , L2 ( ) =fv(x) :Z j v(x)j2 dx0 which would not be derivable in the ordinary sense. Using the new dis- tributional denition however we have that, Z+1 1@ v@ x jw (x) dx=Z +1 1v (x)@ w@ x d x= Z 0 10 @ w@ x d x+Z +1 0x @ w@ x d x Using integration by parts and noticing that the boundary terms cancel out due tow(x)2Dresults in, Z+1 1@ v@ x jw (x) dx=Z 0 10 w(x) dx+Z +1 01 w(x) dx=( 0x0 1x >0 According to this new denition we are able to assign a piece-wise distri- butional derivative to a piece-wise discontinuous function. It is easy to verify that for continuous functions this new denition of derivative is equivalent to the ordinary one. A useful result that holds for all functional spaces is the so-calledCauchy- Schwartz Inequality, (vjw)  k vk k wk 8 v; w2V Chapter 7 The Strong and Weak Formulation The problem we will be trying to solve is exactly the same as in the previous part of the course except we will always assume homogeneous D.B.C., so that the complete formulation of the problem becomes: Given(x);b(x); (x); f(x) ndu(x) such that, ((x)r2 u(x) +b(x)
ru(x) +(x)u(x) =f(x)x2 u(x) = 0x2@ (7.1) This is called theStrong Formulationof our problem, in which the solutionu(x), called theStrong Solution, satises both relationships for any individual pointx2 . In order for this formulation to make sense we need(x);b(x); (x); f(x) to becontinuous functionsbelonging to the spaceC0 ( ). Moreover, when the strong solutionu(x) exists it must belong to the spaceC2 ( ) oftwice- dierentiable functions. It is however very uncommon to meet both of these conditions in any real application. When studying the ow of blood in the carotid artery for example, the walls of the vessel the blood is traveling in might suddenly change physical composition, as is the case when the ow splits into many smaller vessels made of less elastic material. This sudden (discontinuous) change aects the viscosity of the blood, related to the parameter(x), and its instantaneous velocity, related tob(x), making these quantities discontinuous as well. In order to be able to deal with a larger set of problems in which both the known parameters and the solution might be discontinuous we need to introduce the so-calledWeak Formulationof our problem. 75 76 CHAPTER 7. THE STRONG AND WEAK FORMULATION To do this we rst introduce a subset of the Hilbert spaceH1 ( ), H1 0( ) = fv(x)2H1 ( ) :v(x) = 08x2@ g(7.2) This is the space of functionsv(x) belonging toH1 ( ), meaning that they and their gradient are square-integrable over , whose value on the boundary @ is zero. Let's now take an arbitrary functionv(x)2H1 0( ) and multiply each term of our original problem by it, (x)r2 u(x)v(x) +b(x)
ru(x)v(x) +(x)u(x)v(x) =f(x)v(x) and take the integral of both sides of the equation over the entire space , Z (x)r2 u(x)v(x) dx+Z b (x)
ru(x)v(x) dx+Z  (x)u(x)v(x) dx=Z f (x)v(x) dx Exploiting the denition of the Laplace operator we can then write, Z (x)r
(ru)(x)v(x) dx+Z b (x)
ru(x)v(x) dx+Z  (x)u(x)v(x) dx =Z f (x)v(x) dx Using integration by parts on the rst terms we can write, Z  (x)ru(x)rv(x) dxZ @  (x)v(x)ru(x)
ndx +Z b (x)
ru(x)v(x) dx +Z  (x)u(x)v(x) dx=Z f (x)v(x) dx Sincev(x)2H1 0its values at the boundary is null and the second term in the parenthesis disappears. We can now write theWeak Formulation of our original problem: Findu(x)2H1 0( ) such that, Z  (x)ru(x)rv(x) dx+Z b (x)
ru(x)v(x) dx+Z  (x)u(x)v(x) dx (7.3) 77 =Z f (x)v(x) dx8v(x)2H1 0( ) We ask for the solutionu(x), which we will from now on call theWeak Solution, to belong to the subspaceH1 0( ) because the D.B.C. of the original strong formulation still have to be satised.Some doubts may naturally arise about the relationship of the weak for- mulation to the original dierential problem. The rst question we have to answer is if this formulation even makes sense from a mathematical point of view, and under what conditions. The second, perhaps more important, question regards the relationship between the weak solution and the strong one, if the problem even admits one! Conditions on the Given ParametersStating the rst of our issues in a more mathematically precise way, we would like to know if the func- tions inside the integrals are suciently well-behaved as to avoid any of the integrals in the weak formulation from blowing up to innity. We there- fore have to check that the functionsu(x); v(x) and the given parameters (x);b(x); (x); f(x) are regular enough to make sure that these integrals remain bounded. Looking at the rst integral, Z  (x)ru(x)rv(x) dx we know that, if(x)2L1 , meaning that kk 1= sup x2 j (x)j : (x)r2 u(x) +b(x)
ru(x) +(x)u(x) =f(x)x2 u(x) =g(x)x2 D ru(x)
n= (x)x2 N(7.6) where Dand Nare called, respectively, the Dirichletand theNeu- mann Boundary(In a heat equation the rst condition would describe the temperature of a body at the boundary, while the second would describe the heat- ow across it). In order to nd the weak formulation for this problem we will consider the space, H1 D( ) = fv(x)2H1 ( ) :v(x) = 08x2 Dg (7.7) and the set, Vg= fv(x)2H1 ( ) :v(x) =g(x)8x2 Dg (7.8) Through the same line of reasoning used in the homogeneous case we multiply each term of the dierential equation by an arbitrary functionv(x), this time belonging to the spaceH1 D( ), and we integrate over the the whole space obtaining, 80 CHAPTER 7. THE STRONG AND WEAK FORMULATION Z  (x)ru(x)rv(x)dxZ @  (x)v(x)ru(x)
nd  + Z b (x)ru(x)v(x) dx+Z  (x)u(x)v(x) dx=Z f (x)v(x) dx Because the problem now prescribes the non-trivial behavior of the so- lution at the boundary, the integral over@ can be split into two dierent integrals, one over Dand one over N, where the boundary conditions will have to be satised, Z  (x)ru(x)rv(x)dx Z D (x)v(x)ru(x)
nd D+Z N (x)v(x)ru(x)
nd N + Z b (x)ru(x)v(x) dx+Z  (x)u(x)v(x) dx=Z f (x)v(x) dx The integral over Dis equally zero, as the test functions are null over D, while for the integral over the remaining portion of the boundary we can use the second of our boundary conditions. Bringing all known terms to one side allows us to write the weak formulation of our problem: Findu(x)2V g such that, Z  (x)ru(x)rv(x)dx+Z b (x)ru(x)v(x) dx+Z  (x)u(x)v(x) dx= Z f (x)v(x) dx+Z N (x)v(x) d N8 v(x)2H1 D(7.9) The rst thing to notice is that, unlike before, the solutionu(x)2V g and the test functionv(x)2H1 Ddo not belong to the same space , as their behavior at the boundary is dierent. The solution is equal to the given functiong(x) on D, while its gradient is prescribed over N. The test functions instead vanish over Dwhile their behavior on Nis arbitrary. The second thing to notice is that the imposition of D.B.C. over D in uences the spacesV gand H1 Dto which both solution and test functions belong. The imposition of N.B.C. does not directly in uence the choice of spaces in which to look for the solution, but modies the weak formulation introducing a second forcing term over N. Furthermore, in order for the weak formulation to make sense we must require that (x)2L2 (N). Using the denition of scalar product in a functional space we can rewrite the weak formulation in the presence of non-homogeneous boundary condi- tions in a compact form: Findu(x)2V gsuch that, 81 (ru;rv) L2 ( )+( bru; v) L2 ( )+( u; v) L2 ( )= ( f ; v) L2 ( )+( ; v) L2 ( )8 v2H1 D( ) (7.10) 82 CHAPTER 7. THE STRONG AND WEAK FORMULATION Uniqueness of the Weak SolutionNow that we know the conditions under which the weak formulation makes sense and the relationship between the weak solution and the strong one we can investigate whether or not the problem admits a unique solution. Again we will directly report the result without providing any proof. Let's introduce a theorem which in the literature is actually known as the Lax-Milgram Lemma Theorem 7.1(Lax-Milgram Lemma).Consider a Hilbert spaceVand an applicationa(
;
) :VV!R. Moreover consider a functionalF(
) :V! R. Now suppose that, 1.a(
;
)isbilinear, a(v+; z) =a(v; z) +a(w; z) a(v; w+z) =a(v; w) +a(v; z) 2.a(
;
)iscontinuous, a(v; w)Mkvk Vk wk V 3.a(
;
)iscoercive, a(v; v)kvk2 4.F(v)iscontinuous, F(v)ckvk 5.F(v)islinear, F(v+w) =F(v) +F(w) If these hypotheses are satised, then the problem: Findu(x)2Vsuch thata(u; v) =F(v)8v2Vadmits aunique solution. In order to see how this theorem relates to our dierential problem let's consider the fully homogeneous case with D.B.C.: Findu(x)2H1 0( ) such that, (rujrv) L2 ( )+ ( brujv) L2 ( )+ ( ujv) L2 ( )= ( fjv) L2 ( )8 v(x)2H1 0( ) 83 Intuitively we can identify the functional spaceVneeded in order to apply the Lax-Milgram lemma withH1 0( ). Then we could try choosing, a(v; w) = (rvjrw) L2 ( )+ ( brvjw) L2 ( )+ ( vjw) L2 ( )(7.11) F(v) = (fjv) L2 ( )(7.12) and let's see whether these functionals meet the required hypotheses. 1.a(
;
) should bebilinear, a(v+w; z) = (r(v+w)jrz) L2 ( )+( br(v+w)jz) L2 ( )+( (v+w)jz) L2 ( ) =(rvjrz) L2 ( )+ (rwjrz) L2 ( )+ (brvjz) L2 ( ) +(brwjz) L2 ( )+ (vjz) L2 ( )+ (wjz) L2 ( ) =a(v; z) +a(w; z) and the same holds true for the second argument. 2.a(
;
) should becontinuous. Using the Cauchy-Schwartz inequality we can bounda, a(v; w) = (rvjrw) L2 ( )+ ( brvjw) L2 ( )+ ( vjw) L2 ( )  kk L1 ( )k rvk L2 ( )k rwk L2 ( )+ kbk L1 ( )k rvk L2 ( )k wk L2 ( )+ kk L1 ( )k vk L2 ( )k wk L2 ( ) An important identity that ties the norm inL2 ( ) andH1 ( ) is, kvk H1 ( )= kvk2 L2 ( )+ krvk2 L2 ( ) 12  kvk L2 ( ) so that the right-hand term can be further bounded by,  kk L1 ( )k vk H1 ( )k wk H1 ( )+ kbk L1 ( )k vk H1 ( )k wk H1 ( )+ kk L1 ( )k vk H1 ( )k wk H1 ( ) so that, a(v; w) kk L1 ( )+ kbk L1 ( )+ kk L1 ( ) kvk H1 ( )k wk H1 ( ) 84 CHAPTER 7. THE STRONG AND WEAK FORMULATION 3.a(
;
) must becoercive. For physical reasons we shall consider(x)> 0 and we willsuppose(x)>0. Then, a(v; v) =Z  (x)jrv(x)j2 dx+Z b (x)(rv)v(x) dx+Z  (x)v(x)2 dx min x2Rj (x)jZ j rv(x)j2 dx+Z b (x)(rv)v(x) dx+ min x2Rj (x)jZ v (x)2 dx Callingthe smallest number between min x2Rj (x)jand min x2Rj (x)jand rewriting (rv)v(x) as12 r (v(x))2 we can further bound the last term and integrate by parts, =kvk2 H1 ( )+12 Z b r(v(x))2 dx =kvk2 H1 ( )+12 Z @ v (x)2 b(x)
ndx12 Z v (x)2 r
b(x) dx =kvk2 H1 ( )12 Z v (x)2 r
b(x) dx Supposingthatr
b(x) = 0acan therefore be shown to be coercive. 4.F(v) must becontinuous. Using the Cauchy-Schwartz inequality, F(v) =Z f (x)v(x) dx kfk L2 ( )k vk L2 ( )= Ckvk L2 ( ) 5.F(v) must belinear, F(v+w) =Z f (v(x)+w(x)) dx=Z f (x)v(x) dx+Z f (x)w(x) dx=F(v)+F(w) Using the denitions 7:11,7:12 for the functionalsaandFallow for all of the conditions for the application of the Lax-Milgram lemma to be satised and our weak formulation can therefore be put into the form required for the theorem. 85 SummaryIn summary, our weak formulation of the initial dierential equation makes sense and admits a unique solution if, (x); (x);b(x)2L1 ( ) f(x); (x)2L2 ( )  >0;r
b= 0(7.13) Consider now the non-homogeneous problem: Findu(x) such that, 8 > < > : (x)r2 u(x) +b(x)
ru(x) +(x)u(x) =f(x)8x2 u(x) = 08x2 D (x)(ru)(x)
n= (x)8x2 N It's respective weak formulation would therefore be: Findu(x)2H1 D such that, (rujrv) L2 ( )+( brujv) L2 ( )+( ujv) L2 ( )= ( fjv) L2 ( )+( jv) L2 ( )8 v(x)2H1 D( ) If we were to repeat the demonstration for the coercivity ofa, identied with the left-hand side of the equation, we would get the following inequality, a(v; v)kvk H1 ( )+12 Z Nv (x)2 b(x)
ndx Unlike before the test functionsv(x) do not vanish over N, meaning that in order to be able to apply the Lax-Milgram lemma we would need to meet the further requirement, b
n08x2 N 86 CHAPTER 7. THE STRONG AND WEAK FORMULATION Chapter 8 The Galerkin Approximation The purpose of this course, as was the case with the nite-dierences method, is to modify the relevant equations in order to create a discrete numerical problem whose solution is an approximation of the analytical one. Through the Lax-Milgram lemma we modied the original statement of our problem in order to obtain it's weak formulation. As we will see the purpose of going through these steps is that it allows us to nally introduce the discretizations that allows us to obtain the numerical approximated solution of the problem. Let's consider the general weak problem: Findu(x)2Vsuch that, a(u; v) =F(v)8v(x)2V whereVis an appropriate Hilbert space. We now introduce thenite- dimensional subspaceV h V(remembering thatV, like all Hilbert spaces, is innite-dimensional). Similarly to nite-dierences methods, the parameter his a characteristic quantity (usually themesh size) quantifying the degree to which we are approximating the original spaceVwith our new spaceV h. This means that, ash!0 we must haveV h! V. In other terms we could write, lim h!0inf vh2 V hk v(x)v h( x)k V= 0 8v(x)2V From functional analysis it is possible to demonstrate that a subspace of a Hilbert space is still a Hilbert spacewith the same norm. ThusV his a Hilbert space with normkk V. This fact allows us to write a new discrete version of the weak formulation, called theGalerkin Approximation, of our original problem: Findu h( x)2V hsuch that, a(u h; v h) = F(v h) 8v h( x)2V h(8.1) 87 88 CHAPTER 8. THE GALERKIN APPROXIMATION It is important to note that in this case neither of the functionalsa(
;
) nor F(
) are approximated, choosing to approximateVand its test functionv2 VwithV hand its test functions v h2 V h, exclusively. This presents a trade- o as the smallerV his both the space of to which the test functions belong, meaning that equation 8:1 has to hold for a smaller number of cases, and the space in which to look for the solution, meaning that we are introducing restrictions onu h. Since the existence and uniqueness of the solution to the weak problem depended on the conditions set ona(
;
) andF(
) by the Lax- Milgram lemma, a natural follow-up step is to check whether these conditions still hold in our new subspace. Lax-Milgram ConditionsTaking for example the continuity condition set upona(
;
) in the original Hilbert spaceV, which required that9M >0 such that, a(v; w)Mkvk Vk wk V8 w(x); v(x)2V in the new subspaceV hthis would become, a(v h; w h) ~ Mkv hk Vhk w hk Vh8 w h( x); v h( x)2V h Since the rst equation held true8w(x); v(x)2V, it must also hold true for anyw h( x); v h( x) belonging to a subset ofV, such asV hwith the same value of the constantM. The same line of reasoning can be applied to all other conditions both ona(
;
) andF(
). Thus the Galerking approximation admits a unique solutionu hunder the same conditions of the original weak formulation. Numerical CharacteristicsWe can now ask ourselves what relationships tie the numerical solutionu hto the analytical one u. First we check whether or not this method is consistent,The local truncation error is dened as, h( v h) = F(v h) a(u; v h) but sinceu(x) satised the weak formulation8v(x)2Vwe have that, in particular, a(u; v h) = F(v h) 8v h2 V h meaning that, h( v h) = 0 8v h2 V h 8.1. C  EA'S LEMMA89 The truncation error is zero independently on the value ofh. Now, since the properties ofu(x) hold over any subset ofVwe might be tempted to think ofu h( x) as being exactlyu(x) since the Lax-Milgram theorem guarantees that there is only one function inVthat is a solution to the problem. This is not the case as the Lax-Milgram theorem also requires that the solution be found in thesame spaceas the test functionsv(x) and in the Galerkin approximation we are only considering a a sub-spaceV h and therefore a sub-group of test functionsv h( x)2V h, while the original weak solutionu(x)2Vmight not necessarily belong toV h. What the Lax- Milgram lemma guarantees in this case is thus that, for a givenV hthere will exist a unique functionu h( x) belonging toV hthat will be the numerical solution of our problem. Checking now for the stability of the solution we have that, a(u h; v h) = F(v h) 8v h( x)2V h Sincea(
;
) is coercive andF(
) is bounded settingu h( x) =v h( x) gives us, ku hk2 V a(u h; u h) = F(u h) Cku hk V ku hk VC Since the method is consistent and unconditionally stable it must also be convergent, meaning that, lim h!0k u h uk= 0 This brings us to another import result. 8.1 Cea's Lemma Theorem 8.1(Cea's lemma).Given the weak formulation: Findu(x)2V such that, a(u; v) =F(v)8v(x)2V and the relative Galerkin Approximation: Findu h2 V hsuch that, a(u h; v h) = F(v h) 8v h( x)2V h If both satisfy the Lax-Milgram conditions we have that, 90 CHAPTER 8. THE GALERKIN APPROXIMATION kuu hk VM inf vh2 V hk uv hk V This lemma tells us that the error we make by approximatingu(x) with uh( x) is less than a constant times the minimum distance betweenu(x) and the test functionsv h( x)2V h. As h!0 andV h! Vwe know that the minimum distance of any functionv(x)2Vfromv h( x)2V hmust go to 0. In particular this must hold foru(x)2V, meaning that the approximation error goes to zero we rene the mesh. Proof.Recall that h( v h) = F(v h) a(u; v h) = 0 8v h( x)2V h(8.2) subtracting the statement of the Galerking approximation from the pre- vious equation we get, F(v h) a(u; v h) (a(u h; v h) F(v h)) = 0 and using the bi-linearity ofa(
;
), a(uu h; v h) = 0 8v h2 V h Let's now take the term, a(uu h; u u h) = a(uu h; u u h) a(uu h; v h) = a(uu h; u v h)+ a(uu h; v h u h) The last term is 0 due toGalerkin orthogonality, and, using coercivity along with linearity, leaves us with, kuu hk VM k uv hk V since this hold true for anyv h2 V hthe thesis follows.8.2 Solving the Galerkin Approximation Since we have chosenV hto be a nite-dimensional space there must 9f 1;  2;  2


g 2 RN h such that, vh( x) =N h X j=1 j' j( x)8v h( x)2V h(8.3) 8.2. SOLVING THE GALERKIN APPROXIMATION 91 where the functionsf' j( x)gform abasisofV h. With nite-dierences methods the numerical problem yielded the value of the solution at the nodes and thus the approximated solution was only dened on these points, leaving us without information about the behavior of the solution in-between two nodes. In this casev h( x), andu h( x) in particular, arefunctionscharacterized by a nite amount of information. It is therefore enough to know the value of the nite-dimensional vectorto characterize vh( x). The dependence of the approximated solution onhcomes in only through the coecients that make upand not through the basis functions. We could thus write, uh( x) =N h X j=1u j' j( x) (8.4) so thatu2RN h is the vector of unknowns we are interested in. The Galerkin approximation therefore becomes : Findu h( x)2V hsuch that, a(u h; v h) = F(v h) 8v h( x)2V h Since this is a linear problem it is then enough to verify whether it is satised exclusively for the basis functions, a(u h; ' j) = F(' j) 8j= 1;


; N h writingu h( x) as a linear combination of the basis functions this can then be written as, a(N h X i=1u i' i; ' j) = F(' j) 8j= 1;


; N h Exploiting the linearity ofa(
;
) we can bring the sum outside and the nal statement of the weak formulation using the Galerkin approximation becomes: Findu2RN h such that, Nh X i=1u ia (' i; ' j) |{z} Aj;i= F(' j) |{z} Fj8 j= 1;


; N h(8.5) which is a problem withN hequations in N hunknowns. Using vector notation we can write a linear system corresponding to the Galerkin approx- imation and the statemetn becomes: Find a vectoru2RN h such that, Au=F(8.6) 92 CHAPTER 8. THE GALERKIN APPROXIMATION whereAis the so-calledGalerkin Matrix. We now have to answer the following questions: 1. Is this linear system solvable? 2. What are the characteristics of the Galerkin matrixA? 3. How do we build the matrixA? As for the rst question, take the vectorvrelated to a test functionv h( x). We have that, vT Av=N h X i=1v i( Av) i=N h X i=1v iN h X j=1A ijv j =N h X i=1v iN h X j=1a (' j; ' i) v j=N h X i;j=1v iv ja (' j; ' i) =a Nh X j=1v j' j;N h X i=1v i' i! =a(v h; v h) kv hk2 V and we realize that, vT Av08v6 = 0 vT Av= 0()v= 0 Ais therefore aDenite Positivematrix, making itinvertible. The linear system is solvable. 8.2. SOLVING THE GALERKIN APPROXIMATION 93 SummaryIn summary the process taking us form the original strong for- mulation to a solvable linear system goes as follows:Take the generic dierential problem, (x)r2 u(x) +b
ru(x) +(x)u(x) =f(x) The associated weak formulation is: Findu(x)2H1 0such that, Z  (x)ru(x)gradv(x) +b
ru(x)v(x) +(x)u(x)v(x) dx=Z f (x)v(x) dx Using Galerking this then becomes: Findu2RN h such that Au=F where, Aij=Z  (x)r' j( x)r' i( x) +b
r' j( x)' i( x) +(x)' j( x)' i( x) dx Fj=Z f (x)' j( x) dx and the functions'(x) for a basis of the subspaceV hchosen for the Galerkin approximation. 94 CHAPTER 8. THE GALERKIN APPROXIMATION Chapter 9 Finite Elements in 1D The integrals in the equation relating to the Galerkin approximation problem remain identical whichever problem we might be trying to solve, whether it is studying the diusion of heat on a sheet of metal, the behavior of a vibrating string or a 3D modeling of the blood's velocity in a particular artery. What characterizes these dierent problems is the choice of functional spacesV hwhich the solution and test functions must belong to. The general weak formulation simply requires that the solution and test functions both belong toH1 ( ) in order for the integrals to have nite values and the problem to make sense. Starting with one-dimensional problems we can restrict ourselves to a particular subset ofH1 (; ]), Xr h( ; ) =fv h( x)2C0 ([; ]) :v h( x)2Pr x2I j8 jg(9.1) These are spaces containing continuous functions over [; ] that are ex- pressible as a polynomial of graderover any sub-intervalI j= [ x j; x j+1] within [; ]. Each choice ofrwill generate a dierent functional space: r= 1 will be the space of piece-wise linear functions,r= 2 will also contain quadratic functions,r= 3 will contain cubic functions etc.. It is uncommon for any higher degree to be chosen due to the fact that a higher order also means a more irregular solution within the interval, leading to oscillations that avoid convergence. It is therefore custom to use low-order polynomials for these kinds of problems, otherwise opting for methods such as spectral elements when higher-order solutions are needed.It is possible to prove that,for any value ofr, the following result holds, lim h!0inf vh2 Xr hk v h vk H1 (;)= 0 8v(x)2H1 (; ) (9.2) meaning that an eventual solution is convergent. Therefore, when solving reals problems, the functional spaceV hneeded in the Galerkin approximation 95 96 CHAPTER 9. FINITE ELEMENTS IN 1D can be chosen as the sub-set ofXr hwhose elements also satisfy the restrains imposed by the eventual D.B.C. of the problem, such as having value 0 on a boundary D Setting the ProblemThe general1Dproblem with homogeneous D.B.C, ((x)u00 (x) +b(x)u0 (x) +(x)u(x) =f(x)8x2(; ) u() =u() = 0 generates a weak problem of the type: Findu(x)2H1 0( ; ) such that, Z  (x)u0 (x)v0 (x)+b(x)u0 (x)v(x)+(x)u(x)v(x) dx=Z  f (x)v(x) dx8v(x)2H1 0( ; ) By restricting ourselves to only use test functions inXr h, the Galerkin approximation of the problem becomes: Findu h2 V hsuch that Z  (x)u0 (x)v0 h( x)+b(x)u0 (x)v h( x)+(x)u(x)v h( x) dx=Z  f (x)v h( x) dx8v h( x)2V h where the spaceV his chosen as, Vh= fv h( x)2Xr h( ; ) :v h( ) =v h( ) = 0g And, as in nite-dierences, the domain [; ] is subdivided intoMin- tervals andM+ 1 nodes. Due to the D.B.C. the value of the solution at the rst and last node is known, leaving us withM1 unknowns, which will lead to a (M1)(M1) linear system. 9.0.1 Linear Finite Elements Let's start by settingr= 1, X1 h( ; ) =fv h( x)2C0 ([; ]) :v h( x)2P1 x2I j8 jg(9.3) Our subspace, Vh= X1 h;0= fv h( x)2X1 h: v h( ) =v h( ) = 0g will be the space of piece-wiselinearfunctions whose value at the bound- ary nodesandis null. A possible basis for this space is given by the 97 so-calledLagrangian Basis Functionsof order 1, also known ashat func- tions, 'i( x) =8 > > > > < > > > > :x x i1h i1x 2(x i1; x i] xi+1 xh ix 2(x i; x i+1] 0 otherwise(9.4) so that,Vh3 v h( x) =N h X j=1 j' j( x) We can see that,ˆThese functions are linear over each interval ˆThey are continuous over each interval ˆ' i( x j) =  ij Owing to the last of these properties we have that, vh( x l) =N h X j=1 j' j( x l) =N h X j=1 j lj=  l The value of a random test function at a nodex lis equal to the coecient lrelated to the l-th basis function that makes up that particular function. Therefore, thel-th component of the vector that identies the function is simply the value of the function at the relative node. In particular this hold true foru h, uh( x i) = u i8 j= 1;


; N h This is important as the vectoruthat solves the nal linear system allows us to directly nd the values of the solutionu h. 9.0.2 Construction of the Linear System As we have seen in the previous chapter the Galerkin approximation can be reduced to a linear problem of the type, Au=F 98 CHAPTER 9. FINITE ELEMENTS IN 1D where the elements of this system are given by,Aij= a(' j; ' i) (9.5) Fj= F(' j) (notice the change in order of the indices within the elements of the system matrix). Using our linear Lagrangian basis functions we can split the bilinear form in three dierent matrices, Aij= a(' j; ' i) =Z   (x)'0 j( x)'0 i( x) dx+Z  b (x)'0 j( x)' i( x) dx+Z   (x)' j( x)' i( x) dx =A ij+ Ab ij+ A ij(9.6) called theStiness Matrix, theAdvection Matrixand theMass Matrix. Focusing on the elements of the last matrix, A ij=Z   (x)' j( x)' i( x) dx we can see that, according to the denition of the Lagrangian basis func- tions given by eq. 9:4, the elements on the main diagonal, A ii=Z   (x)' i( x)2 dx=Z xi+1 xi1 (x)' i( x)2 dx6 = 0 (9.7) are obviously not zero. Similarly for the elements belonging to the rst diagonal in both directions, A i;i1=Z   (x)' i1( x)' i( x) dx=Z xi xi1 (x)' i1( x)' i( x)6 = 0 A i;i+1=Z   (x)' i+1( x)' i( x) dx=Z xi+1 xi (x)' i1( x)' i( x)6 = 0 we have an overlap between the two basis functions, resulting in a non- zero integral. In the case of elements further from the main diagonal, A i;i2=Z   (x)' i2( x)' i( x) dx= 0 99 there is no such interval, making all other elements zero. The Mass Matrix A ijis therefore a tridiagonal matrix. By noticing that the rst derivative of the basis functions have the same behavior with respect to the overlap of functions of dierent index we can immediately see that the other two matricesA ijand Ab ijhave the same symmetry, making the overall system matrixAa (sparse) tridiagonal matrixas well. 9.0.3 The Reference Interval Given that we are dealing with piece-wise linear functions it might be useful to consider the individual intervals while calculating the matrix elements. In general the basis functions are not null only for values of the global coordinate xbelonging to a specic interval, called theCurrent Interval, x2[x j; x j+1] = I j(9.8) whose length ish j= x j+1 x j. We now introduce a Reference Interval in the local coordinate is, 2[0;1] =I ref(9.9) whose length is unitary. For a given value ofjwe would like to uniquely associate the points in the coordinatexbelonging to the current interval to the ones in the coordinatebelonging to our reference interval. We can achieve this through the use of alinear map, j( ) : [0;1]![x j; x j+1] 8j= 0;1;


; M1 so that, x() =  j( ) =x j+ h j  2[0;1] (9.10) (x) =  1 (x) =x x jh jx 2[x j; x j+1] (9.11) TheReference Domainwill be generated by a pair of basis functions (^ ' 0( ) = 1 ^ ' 1( ) =(9.12) so that, in the interval [x j; x j+1] the relationship between the global and the local basis functions is, 'j( x) = ^' 0( ) 100 CHAPTER 9. FINITE ELEMENTS IN 1D 'j+1( x) = ^' 1( ) Looking at these relationships and equation 9:11 we can see that the dierent basis functions on each interval are simply the local basis functions transposed and stretched by dierent amounts depending on the value ofj. Assuming that the functions(x) =; b(x) =b; (x) =areconstant across the whole domain we can now rewrite the generic non-zero element Ai;i1=Z  ' 0 i1( x)'0 i( x) +b'0 i1( x)' i( x) +' i1( x)' i( x) dx =Z xi xi1' 0 i1( x)'0 i( x) +b'0 i1( x)' i( x) +' i1( x)' i( x) dx and, exchanging variables, x=x i1+ h i1 dx=h i1d  this becomes, =Z 1 0 '0 i1( x())'0 i( x()) +b'0 i1( x())' i( x()) +' i1( x())' i( x()) hi1d  Before exchanging the current basis functions for the local ones we must calculate their derivatives, '0 i1( x) = ^'0 0( (x))@ @ x = ^ '0 0( (x))1h i1 'i( x)0 = ^'0 1( (x))@ @ x = ^ '0 1( (x))1h i1 which, similarly to their primitives, are simply stretched versions of the reference ones. The generic element thus becomes, Ai;i1=Z 1 0 ( h i1)2^ '0 0( ) ^'0 1( ) +bh i1^ '0 0( ) ^' 1( ) +^ ' 0( ) ^' 1( ) hi1d  =Z 1 0h i1^ '0 0( ) ^'0 1( ) +b^ '0 0( ) ^' 1( ) +h i1^ ' 0( ) ^' 1( ) d 101 We can then split the integral in three and bringing all constant terms outside the integral we obtain, =h i1Z 1 0^ '0 0( ) ^'0 1( ) d+bZ 1 0^ '0 0( ) ^' 1( ) d+h i1Z 1 0^ ' 0( ) ^' 1( ) d Obtaining three equivalent integrals to the ones obtained previously, this time in the reference interval [0;1]. The dependence oni, so the specic interval we are considering, is found only withing the constanth i1. This allows for these integrals to be computed only once at the beginning of the problem and then simply multiplied by a constant when needed, allowing for this method to be particularly robust when it comes to changes in the mesh.We have chosen the elementA i;i1at random between the non-zero ele- ments on the system matrixA, which we know to be tridiagonal. In general, we can write a local version of the previous matrices, such as theLocal Stiness Matrix, ^ A =h i1 ^ A 00^ A 01 ^ A 10^ A 11 =h i10 B B B B @R 1 0( ^ '0 0)2 () dR 1 0^ '0 0( ) ^'0 1( ) d R1 0^ '0 1( ) ^'0 0( ) dR 1 0( ^ '0 1)2 () d1 C C C C A= h i1 11 1 1 (9.13) theLocal Advection Matrix, ^ Ab =b ^ Ab 00^ Ab 01 ^ Ab 10^ Ab 11 =b0 B B B B @R 1 0^ '0 0( ) ^' 0( ) dR 1 0^ '0 0( ) ^' 1( ) d R1 0^ '0 1( ) ^' 0( ) dR 1 0^ '0 1( ) ^' 1( ) d1 C C C C A= b 12 12 12 12  (9.14) and the Local Mass Matrix, ^ A =h i1 ^ A 00^ A 01 ^ A 10^ A 11 =h i10 B B B B @R 1 0^ '2 0( ) dR 1 0^ ' 0( ) ^' 1( ) d R1 0^ ' 1( ) ^' 0( ) dR 1 0^ '2 1( ) d1 C C C C A= h i1 13 16 16 13  (9.15) 102 CHAPTER 9. FINITE ELEMENTS IN 1D These matrices, aside form the constant terms, areindependent on the mesh. 103 ExampleUsing these matrices the elements on the lower and upper diag- onal become, Ai;i1=h i1^ A 10+ b^ Ab 10+ h i1^ A 10= h i1 b2 + h i16 Ai;i+1=h i^ A 01+ b^ Ab 01+ h i^ A 01= h i+ b2 + h i6 For the elements on the main diagonal we get two dierent integrals to compute, Aii=Z   ('0 i( x))2 +b'0 i( x)' i( x) +(' i( x))2 dx =Z xi xi1 ('0 i( x))2 +b'0 i( x)' i( x)+(' i( x))2 dx+Z xi+1 xi ('0 i( x))2 +b'0 i( x)' i( x)+(' i( x))2 dx (9.16) For the rst integral we can identify, 'i( x) = ^' 1( ) while for the second, 'i( x) = ^' 0( ) giving us, Aii=Z 1 0h i1( ^ '0 1)2 +b^ '0 1^ ' 1+ h i1( ^ ' 1)2 d+Z 1 0h i( ^ '0 0)2 +b^ '0 0^ ' 0+ h i( ^ ' 0)2 d which then becomes,=h i1^ A 11+ b^ Ab 11+ h i1^ A 11+h i^ A 00+ b^ Ab 00+ h i^ A 00 = 1h i1+ 1h i +3 ( h i1+ h i) Finite elements poses more importance on the construction of the system matrix rather than the solution of the linear system. At the end of the story this method turn out to be much more complicated than nite dierences but also much more robust method with respect the the geometry. 104 CHAPTER 9. FINITE ELEMENTS IN 1D 9.1 Quadratic Finite Elements As we know, interpolation with higher-order polynomials deals more accurate results when compared to using simply linear functions. We might then be temped to setr= 2, X2 h( ; ) =fv h( x)2C0 ([; ]) :v h( x)2P2 x2I j8 jg(9.17) Our subspace, Vh= X2 h;0= fv h( x)2X2 h: v h( ) =v h( ) = 0g will be the space of quadratic functions whose value atandis zero. Since we are dealing with second-order polynomials for each interval we will need to provide, other than the boundary values, which we will denote asx 2i andx 2i+2, also the mid-point of the interval, which will be x 2i+1, in practice doubling the number of unknowns. While being more accurate, this leads to a bigger quadratic system that is in general harder to solve than the linear one. The Lagrangian basis functions of this new space, whose intervals are given byI i= [ x 2i; x 2i+2], must belong to X2 hand satisfy the property, X2 h3 ' 2i( x j) =( 1j= 2i 0 otherwise The only type of functions that satisfy this condition are parabolic func- tions over the intervalsI i1and I i. These basis functions will now have non-zero values over 5 nodes instead of the 3 of the linear case, aecting the shape of the system matrix. We notice however that if we were to transpose these functions by an interval we would have, '2i1( x j) =( 1j= 2i1 0 otherwise These functions however do not belong to our chosen functional space X2 h, as in the interval I i1= [ x 2i2; x 2i] the function is a piece-wise union of two parabolas. We therefore have to introduce a dierent family of basis functions for odd nodes that satisfy the condition above. 9.1. QUADRATIC FINITE ELEMENTS 105 9.1.1 Construction of the Linear System The process to build the system matrices is essentially the same as the linear case, with the added complexity of having to consider three local basis func- tions. Let's rst consider the elements on aneven row. Taking the Mass matrix as an example, we obviously have that, A2i;2i6 = 0 meaning that the elements on even row of the main diagonal are not null. Due to the fact that the basis functions stretch over 5 nodes we have more overlaps between basis functions with dierent indices and therefore more elements in the matrix to compute. In particular, A2i;2i1=Z 2i1 2i' 0 2i1( x)'0 2i( x) +b'0 2i1( x)' 2i( x) +' 2i1( x)' 2i( x) dx6 = 0 A2i;2i2=Z 2i2 2i' 0 2i2( x)'0 2i( x) +b'0 2i2( x)' 2i( x) +' 2i2( x)' 2i( x) dx6 = 0 Other than the lower and upper diagonal having non zero elements the diagonals at distance 2 from the main one are also not identically zero. On the other diagonals, A2i;2i3=Z 2i3 2i' 0 2i3( x)'0 2i( x) +b'0 2i3( x)' 2i( x) +' 2i3( x)' 2i( x) dx= 0 we have no further overlaps between basis functions. Let's now take the elements onodd rows. As before on the main diagonal we have that, A2i1;2i16 = 0 On the upper and lower diagonals the integrals are also not null A2i;2i1=Z 2i1 2i' 0 2i1( x)'0 2i( x) +b'0 2i1( x)' 2i( x) +' 2i1( x)' 2i( x) dx6 = 0 while for more distant diagonals,A2i1;2i3= 0 106 CHAPTER 9. FINITE ELEMENTS IN 1D On odd rows, only the upper,lower and main diagonal contain non-zero elements. Remembering that all three matrices have the same symmetry we con- clude that the overall system matrixAis penta-diagonal with 5 elements in evenrows and 3 elements inoddrows, A=0 B B B B B B B B B B B B @a 11a 120 0 0 0 0 0 0 a21a 22a 23a 240 0 0 0 0 0a 32a 33a 340 0 0 0 0 0a 42a 43a 44a 45a 460 0 0 0 0 0a 54a 55a 560 0 0 0 0 0a 64a 65a 66a 67a 680 0 0 0 0 0a 76a 77a 780 0 0 0 0 0a 86a 87a 88a 89 0 0 0 0 0 0 0a 98a 991 C C C C C C C C C C C C A As before, we introduce aReference Intervalin the local coordinate is , 2[0;1] =I ref whose length is unitary. We notice that the map j: [0 ;1]!I jfrom the current to the reference interval is the same as in the linear case. What has changed are the basis functions, which in this case will be, 8 > < > :^ ' 0( ) = (1)(12) ^ ' 1( ) = 4(1) ^ ' 2( ) =(21)(9.18) Again, the result of our process is the creation of the Local Stiness, Advection and Mass matrices as before, which will now contain 9 elements each whose values are xed and must only be computed once, ^ A =c0 @A  00A 01A 02 A 10A 11A 12 A 20A 21A 221 A=c0 @73 83 13 83 163 83 13 83 73 1 A(9.19) etc.. Since we have two possibly overlapping integrals for each pair of function the way to compute the elements is to split the initial integral into two parts A 2i;2i=Z x2i+2 x2i2 ( ^'0 2i( ))2 d=Z x2i x2i2 ( ^'0 2i( ))2 d+Z x2i+2 x2i ( 9.1. QUADRATIC FINITE ELEMENTS 107 A 2i;2i+1=Z x2i+2 x2i' 0 2i+2( )'0 2i( ) d=h i^ A 01= 8 3 h i A 2i;2i+2=h i^ A 02=3 h i For the even rows, A 2i1;2i1=Z x2i x2i2 ('0 2i1( ))2 d=h i1^ A 11=16 3 h i1 These matrices will all be sparse and it can be proven that they are badly conditioned, K(A )h 2 K(Ab )h 2 K(A )h0 calling for the use of pre-conditioners. 108 CHAPTER 9. FINITE ELEMENTS IN 1D Chapter 10 Finite Elements in 2D The passage from one to multiple dimensions is, at least conceptually, easy. We have to dene the new domain , create a mesh hand nd an ap- propriate subsetV hof test functions that solve our weak formulation. For two-dimensional problems we usually deal with triangular or, more rarely, quadrilateral meshes, that is, meshes made up of a union of many triangles or quadrilaterals, identied through the use of an indexK, so that[K= h . The subdivision of the domain into regular shapes usually introduces an error. We can imagine that in the case of a linear shape, such as a cylinder or a sphere, the mesh hwill not be able to cover the domain in it's entirety, particularly at the boundary. This error however, called the geometric errorcan be shown to be of the same degree or smaller than the approximation error introduced by the Galerkin method, and is usually not an issue. Alternate meshes can be considered whenever the problem calls for the domain to be covered with greater precision as is the case with the so-called parametric nite elements method, which uses a mesh made up of triangles with curvilinear boundaries. The parameterhquantifying the mesh size can be taken as, h= max Kh K where, hK= max (x;y)2Kk xyk R2 which in a triangular mesh is simply the longest segment making up the side of a triangle in the mesh. Similarly to the one-dimensional case we can take a particular nite element space, subset ofH1 ( ), namely, 109 110 CHAPTER 10. FINITE ELEMENTS IN 2D Xr h= fv k( x; y)2C0 ( K) : v h( x; y)2Pr (x; y)2 K; 8K2 hg where the polynomials are now of the type, Pr =f(x; y) =X i+jra ijxi yj for which it can be shown that,dim(Pr ) =( r+ 1)(r+ 2)2 The space of test functionsV h Xr hwill now be the space of functions that, over each triangleKcan be written as polynomials of at most degree rand that satisfy eventual D.B.C.. 10.0.1 Linear Finite Elements Settingr= 1 we can look for Linear Lagrangian Basis Functions' j( x; y)2V h such that, 'j( x i; y i) =  ij which will generate a sparse system matrixA. To build this matrix we will introduce aReference Trianglein the local coordinates; and it's corresponding Local Basis Functions, 8 > < > :^ ' 1( ; ) = 1 ^ ' 2( ; ) = ^ ' 3( ; ) = Calling ( ^x 1; ^ y 1) ;( ^x 2; ^ y 1) ;( ^x 3; ^ y 3) the coordinates of the vertices of the ref- erence triangle and (x 1K; y 1K) ;(x 2K; y 2K) ;(x 3K; y 3K) those of the current tri- angleK, the map  K( ; ) is given by, K( ; )!( x(; ) =x 1K^ ' 1( ; ) +x 2K^ ' 2( ; ) +x 3K^ ' 3( ; ) y(; ) =y 1K^ ' 1( ; ) +y 2K^ ' 2( ; ) +y 3K^ ' 3( ; ) In order to perform the change in variables within each double integral we will need to compute the Jacobian of this transformation, r(;)' j= jJ Kj r (x;y)' j 10.1. ACCURACY OF FINITE ELEMENT METHODS 111 which is, JK= @ x@  @ y@  @ x@  @ y@  ! = x2 x 1y 2 y 1 x3 x 1y 3 y 1 (10.1) jJ Kj = 2jKj allowing us to exchange, dxdy=12 jKjd d in each integral. The process is then similar to the one-dimensional case and will not be explored further. We can say that the construction of the system matrix is not trivial and requires great attention. 10.1 Accuracy of Finite Element Methods Cea's lemma guarantees that nite elements methods are convergent, kuu hk VM inf vh2 V hk uv hk V It is then possible to prove that, 8v(x)2H1 ( ) : lim h!0inf vh2 Xr hk vv hk H1 = 0 8r so that, kuu hk H1 ( ) Cinfkuv hk H1 ( ) Moreover, if the mesh his regular, that is,8K;9 >0 such that, hk k  where kis the radius of the circle inscribed in the triangle, we have that, kuu hk H1 ( ) ~ C hr Dr +1 u L2 ( )= O(hr ) wherekDr +1 uk L2 ( )is the norm of the derivative of order r+ 1 of the solutionu(x).The Galerkin approximation of orderris a method of orderr. This result however assumes thatDr +1 u(x)2L2 ( ), but the solu- tionu(x) is unknown making it impossible to directly verify this condition. If we were to have for example, 112 CHAPTER 10. FINITE ELEMENTS IN 2D D3 u(x)62L2 ( ) D2 u(x)2L2 ( ) usingr= 2 would not yield better convergence properties than using r= 1. Thankfully, there are a number of important results that allow us to reconstruct the regularity of the solutionu(x)a-priorithrough the regularity of the given data, as well as that of the domain. This allows us to select the appropriate functional space for our numerical problem. Using the example, (r2 u(x) =f(x)x2 u(x) = 0x2@ if is smooth andDs f(x)2L2 ( ) it can be proven that, Ds +2 u(x)2L2 ( ) allowing us to reconstruct the dierential properties of the solution at the beginning of the problem. Intuitively this result is tied to the fact that the Laplace operator acting onu(x), a sort of second derivative, is set by the rst condition of the system equal to a known forcing termf(x). Should this term belong toL2 ( ) we can imagine that the solution has the above mentioned properties. Iff(x)2L2 ( ) butf(x)62H1 ( ) according to the result we would have thatD2 u(x)2L2 ( ) butD3 u(x)62L2 ( ). This result holdonly for smooth . If the domain is not smooth we cannot infer the regularity of the solution fromf(x). 10.2 Increasing r or h? Since, kuu hk H1 ( ) ~ C hr Dr +1 u L2 ( )= O(hr ) it is in principle more advantageous to increase the orderrof the method than using a ner mesh (samllerh), as the error tend to be smaller is the former case while having a system matrix of the same dimensions. The caseBandCfeature the same matrices dimensions. However, should our domain have sharp corners the solutionu(x) might not be regular enough to allow us to suciently increaser, forcing us to rely onhas the only parameter that can increase the accuracy of our method. Moreover we have that, 10.2. INCREASING R OR H? 113 K(A )h 2 r4 which forr= 1 andh0 =h2 becomes, K(A ) h02  2 =4h 2 while forr= 2 andh0 = 2 is, K(A )h0 2 r4 =16h 2 As usual increased accuracy comes with greater computational eort. Diusion-Advection ProblemFormally, Cea's lemma tells us that the error goes to zero ashapproaches zero, which in real scenarios is never the case;his small, but never actually zero. The approximation error then will depend on the factorM . In a Diusion-Advection problem we have that, M=kk L1 + kbk L1 = ~= infkk meaning that, M  k k L1 + kbk L1~  and the amountM , along with the error, ends up being large when the advection dominates over the diusion, similarly to the nite-dierences case. Through the same line of reasoning we can dene a Peclet number P e=bh2 < 1!ok P e=bh2 > 1!oscillation!numerical viscosity 114 CHAPTER 10. FINITE ELEMENTS IN 2D Diusion-Reaction ProblemIn a Diusion-Reaction problem we in- stead have that, M=kk L1 + kk L1 = minfinfkk;infkkg so that the amountM ends up being large if the reaction dominates over the diusion, unlike in the nite-dierences case, where we had no particular condition set on the initial data. We can then dene a dierent kind of Peclet number specic to this kind of problem, P e=h 26  which has the property, P e