An Introduction to the FEM for Differential Equations.pdf -...

This preview shows page 1 out of 307 pages.

Unformatted text preview: An Introduction to the Finite Element Method (FEM) for Differential Equations Mohammad Asadzadeh January 12, 2016 Contents 0 Introduction 0.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.2 Trinities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Partial Differential Equations 1.1 Differential operators, superposition . . . 1.1.1 Exercises . . . . . . . . . . . . . . 1.2 Some equations of mathematical physics 1.2.1 Exercises . . . . . . . . . . . . . . 7 8 9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 19 22 23 33 2 Polynomial Approximation in 1d 2.1 Overture . . . . . . . . . . . . . . . . . . . 2.2 Galerkin finite element method for (2.1.1) 2.3 A Galerkin method for (BVP) . . . . . . 2.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 35 46 50 60 . . . . . . 63 63 71 75 78 82 86 3 Interpolation, Numerical Integration in 1d 3.1 Preliminaries . . . . . . . . . . . . . . . . . . 3.2 Lagrange interpolation . . . . . . . . . . . . . 3.3 Numerical integration, Quadrature rule . . . . 3.3.1 Composite rules for uniform partitions 3.3.2 Gauss quadrature rule . . . . . . . . . 3.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Linear Systems of Equations 91 4.1 Direct methods . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.2 Iterative methods . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.3 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 3 4 5 Two-point boundary value problems 5.1 A Dirichlet problem . . . . . . . . . . . 5.2 A mixed Boundary Value Problem . . . 5.3 The finite element method (FEM) . . . . 5.4 Error estimates in the energy norm . . . 5.5 FEM for convection–diffusion–absorption 5.6 Exercises . . . . . . . . . . . . . . . . . . CONTENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 . 113 . 118 . 121 . 122 . 128 . 137 6 Scalar Initial Value Problems 6.1 Solution formula and stability . . . . . . . . . . . . . 6.2 Galerkin finite element methods for IVP . . . . . . . 6.2.1 The continuous Galerkin method . . . . . . . 6.2.2 The discontinuous Galerkin method . . . . . . 6.3 A posteriori error estimates . . . . . . . . . . . . . . 6.3.1 A posteriori error estimate for cG(1) . . . . . 6.3.2 A posteriori error estimate for dG(0) . . . . . 6.3.3 Adaptivity for dG(0) . . . . . . . . . . . . . . 6.4 A priori error analysis . . . . . . . . . . . . . . . . . 6.4.1 A priori error estimates for the dG(0) method 6.5 The parabolic case (a(t) ≥ 0) . . . . . . . . . . . . . 6.5.1 Some examples of error estimates . . . . . . . 6.6 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177 . 177 . 179 . 185 . 188 . 196 . 198 . 199 . 200 . 203 . 205 . 207 . 209 . 211 . . . . . . . . . . . . . . . . BVPs . . . . . . . . . . . . . . . . 7 Initial Boundary Value Problems in 1d 7.1 Heat equation in 1d . . . . . . . . . . . . . . . . . . . . 7.1.1 Stability estimates . . . . . . . . . . . . . . . . 7.1.2 FEM for the heat equation . . . . . . . . . . . . 7.1.3 Error analysis . . . . . . . . . . . . . . . . . . . 7.1.4 Exercises . . . . . . . . . . . . . . . . . . . . . . 7.2 The wave equation in 1d . . . . . . . . . . . . . . . . . 7.2.1 Wave equation as a system of hyperbolic PDEs 7.2.2 The finite element discretization procedure . . 7.2.3 Exercises . . . . . . . . . . . . . . . . . . . . . . 7.3 Convection - diffusion problems . . . . . . . . . . . . . 7.3.1 Finite Element Method . . . . . . . . . . . . . . 7.3.2 The Streamline-diffusion method (SDM) . . . . 7.3.3 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 146 147 148 151 153 153 160 162 163 163 167 171 174 CONTENTS 5 8 Piecewise polynomials in several dimensions 8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 8.2 Piecewise linear approximation in 2 D . . . . . . . . . . . 8.2.1 Basis functions for the piecewise linears in 2 D . . 8.2.2 Error estimates for piecewise linear interpolation . 8.2.3 The L2 projection . . . . . . . . . . . . . . . . . . 8.3 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 . . . . . . . . . . . . 213 . 213 . 216 . 216 . 223 . 225 . 226 Riesz and Lax-Milgram Theorems 229 9.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . 229 9.2 Riesz and Lax-Milgram Theorems . . . . . . . . . . . . . . . . 234 9.3 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 241 10 The Poisson Equation 10.1 Stability . . . . . . . . . . . . . . . . 10.2 Error Estimates for the CG(1) FEM 10.2.1 Proof of the regularity Lemma 10.3 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 The Initial Boundary Value Problems in RN 11.1 The heat equation in RN . . . . . . . . . . . . . . . . 11.1.1 The fundamental solution . . . . . . . . . . . 11.1.2 Stability . . . . . . . . . . . . . . . . . . . . . 11.1.3 A finite element method for the heat equation 11.1.4 Constructing the discrete equations . . . . . . 11.1.5 An apriori error estimate . . . . . . . . . . . . 11.2 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . 11.3 The wave equation in RN . . . . . . . . . . . . . . . 11.3.1 The weak formulation . . . . . . . . . . . . . 11.3.2 The semi-discrete problem . . . . . . . . . . . 11.3.3 The fully-discrete problem . . . . . . . . . . . 11.3.4 A priori error estimate for the wave equation . 11.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243 . 243 . 245 . 251 . 253 . . . . . . . . . . . . . 257 . 257 . 258 . 259 . 262 . 263 . 264 . 264 . 268 . 269 . 269 . 270 . 271 . 271 12 Algorithms and MATLAB Codes 283 Table of Symbols and Indices 303 6 CONTENTS xs Chapter 0 Introduction This book presents an introduction to the Galerkin finite element method (FEM) as a general tool for numerical solution of differential equations. Our objective is to construct and analyze some simple FEMs for approximate solutions of both ordinary, and partial differential equations (ODEs and PDEs). In its final step, a finite element procedure yields a linear system of equations (LSE) where the unknowns are the approximate values of the solution at certain nodes. Then an approximate solution is constructed by adapting piecewise polynomials of certain degree to these nodal values. The entries of the coefficient matrix and the right hand side of FEM’s final linear system of equations consist of integrals which, e.g. for complex geometries or less smooth data, are not always easily computable. Therefore, numerical integration and quadrature rules are introduced to approximate such integrals. Furthermore, iteration procedures are included in order to efficiently compute the numerical solutions of such obtained matrix equations. Interpolation techniques are presented for both accurate polynomial approximations and also to derive basic a priori and a posteriori error estimates necessary to determine qualitative properties of the approximate solutions. That is to show how the approximate solution, in some adequate measuring environment, e.g. a certain norm, approaches the exact solution as the number of nodes, hence unknowns, increase. Some theoretical aspects as existence, uniqueness, stability and convergence are discussed as well. Mathematically, Galerkin’s method for solving a general differential equation is based on seeking an approximate solution, which is 7 8 CHAPTER 0. INTRODUCTION 1. Easy to differentiate and integrate 2. Spanned by a set of “nearly orthogonal” basis functions in a finitedimensional vector space. 3. Satisfies Galerkin orthogonality relation. Roughly speaking, this means that: the difference between the exact and approximate solutions is orthogonal to the finite dimensional vector space of the approximate solution. 0.1 Preliminaries In this section we give a brief introduction to some key concepts in differential equations. A more rigorgous and thorough introduction will be presented in the following Chapter 1. • A differential equation is a relation between an unknown function u and its derivatives. • If the function u(x) depends on only one variable (x ∈ R), then the equation is called an ordinary differential equation (ODE). Example 0.1. As a simple example of an ODE we mention the population dynamics model du (t) − λu(t) = f (t), t > 0. (0.1.1) dt If f (t) ≡ 0, then the equatios is clalled homogeneous, otherwise it is called inhomogeneous. For a nonnegative λ, the homogeneous equation du −λu(t) = dt 0 has an exponentially growing analytic solution given by u(t) = u(0)eλt , where u(0) is the initial population. • The order of a differential equation is determined by the order of the highest derivative of the function u that appears in the equation. • If the function u depends on more than one variable, and the differential equation posseses derivatives with respect to al least two variables, then the differential equation is called a partial differential equation (PDE), e.g. ut (x, t) − uxx (x, t) = 0, 0.2. TRINITIES 9 is a homogeneous PDE of the second order, whereas for f 6= 0, the equations uxx (x, y) + uyy (x, y) = f (x, y), and ut (x, y, t) − uxx (x, y, t) − uyy (x, y, t) + ux (x, y, t) + uy (x, y, t) = f (x, y, t), are non-homogeneous PDEs of the second order. • A solution to a differential equation is a function; e.g. u(x), u(x, t), u(x, y) or u(x, y, t) in example (0.1), that satisfies in the differential equation. • In general the solution of a differential equation cannot be expressed in terms of elementary functions and numerical methods are the only way to solve the differential equations through constructing approximate solutions. Then, the main questions are: how close is the computed approximate solution to the exact solution? (convergence), and how and in which environment does one measure this closeness? In which extent does the approximate solution preserve the physical properties of the exact solution? How sensitive is the solution to the change of the data (stability) These are some of the questions that we want to deal within this text. • A linear ordinary differential equation of order n has the general form: L(t, D)u = u (n) (t) + n−1 X ak (t)u(k) (t) = b(t), k=0 k where D = d/dt denotes the derivative, and u(k) := Dk u, with Dk := dtd k , 1 ≤ k ≤ n (the k-th order derivative). The corresponding linear differential operator is denoted by L(t, D) = 0.2 dn−1 d dn + a (t) + . . . + a1 (t) + a0 (t). n−1 n n−1 dt dt dt Trinities Below we introduce the so called trinities classifying the main ingredients involved in the process of identifying problemes that are modeled as partial differential equations of the second order, see Karl E. Gustafson [25] for details. 10 CHAPTER 0. INTRODUCTION The usual three operators in differential equations of second order: Laplace operator ∆n := ∂2 ∂2 ∂2 , + + . . . + ∂x21 ∂x22 ∂x2n (0.2.1) ∂ − ∆n , (0.2.2) ∂t ∂2 d’Alembert operator  := 2 − ∆n , (0.2.3) ∂t where we have the space variable x := (x1 , x2 , x3 , . . . , xn ) ∈ Rn , the time variable t ∈ R+ , and ∂ 2 /∂x2i denotes the second partial derivative with respect to xi , 1 ≤ i ≤ n. We also define a first order operator, namely the gradient operator ∇n which is the vector valued operator  ∂ ∂ ∂  . , ,..., ∇n := ∂x1 ∂x2 ∂xn Diffusion operator Classifying general second order PDEs in two dimensions • Second order PDEs in 2D with constant coefficients A general second order PDE in two dimensions and with constant coefficients can be written as Auxx (x, y)+2Buxy (x, y)+Cuyy (x, y)+Dux (x, y)+Euy (x, y)+F u(x, y)+G = 0. Here we introduce the concept of discriminant d = AC − B 2 : The discriminant is a quantity that specifies the role of the coefficients, of the terms with two derivatives, in determining the equation type in the sence that the equation is: Elliptic: if d > 0, Parabolic: if d = 0, Hyperbolic: if d < 0. Example 0.2. Below are the classes of the most common differential equations together with examples of their most simple forms in 2D: Potential equation Heat equation ∂u − ∆u = 0 ∂t ut − uxx = 0 Wave equation ∂ 2u − ∆u = 0 ∂t2 utt − uxx = 0 A = C = 1, B = 0 B = C = 0, A = −1 A = −1, B = 0, C = 1 d = 1 (elliptic) d = 0 (parabolic) d = −1 (hyperbolic). ∆u = 0 uxx + uyy = 0 0.2. TRINITIES 11 • Second order PDEs in 2D with variable coefficients In the variable coefficients case, one can only have a local classification. Example 0.3. Consider the Tricomi equation of gas dynamics: Lu(x, y) = yuxx + uyy . Here the coefficient y is not a constant and we have A = y, B = 0, and C = 1. Hence d = AC − B 2 = y and consequently, the domain of ellipticity is y > 0, and so on, see Figure 1. y elliptic parabolic x hyperbolic Figure 1: Tricomi equation: an example of a variable coefficient classification. The usual three types of problems in differential equations 1. Initial value problems (IVP) The simplest differential equation is u′ (x) = f (x) for a < x ≤ b. But for any such u, also (u(x) + c)′ = f (x) for any constant c. To determine a unique solution a specification of the initial value u(a) = u0 is generally required. For example for f (x) = 2x, 0 < x ≤ 1, we have u′ (x) = 2x and the general solution is u(x) = x2 + c. With an initial value of u(0) = 0 we get u(0) = 02 + c = c = 0. Hence the unique solution to this initial value problem is u(x) = x2 . Likewise for a time dependent differential equation of 12 CHAPTER 0. INTRODUCTION second order (two time derivatives) the initial values for t = 0, i.e. u(x, 0) and ut (x, 0), are generally required. For a PDE such as the heat equation the initial value can be a function of the space variable. Example 0.4. The wave equation, on the real line, augmented with the given initial data: u − u = 0, −∞ < x < ∞, t > 0, tt xx u(x, 0) = f (x), u (x, 0) = g(x), −∞ < x < ∞, t = 0. t Remark 0.1. Note that, here, for the unbounded spatial domain (−∞, ∞) it is required that u(x, t) → 0 (or tou∞ =constant) as |x| → ∞. This corresponds to two boundary conditions. 2. Boundary value problems (BVP) Example 0.5. (a boundary value problems in R). Consider the stationary heat equation:  ′ − a(x)u′ (x) = f (x), for 0 < x < 1. In order to determine a solution u(x) uniquely (see Remark 0.2 below), the differential equation is complemented by boundary conditions imposed at the boundary points x = 0 and x = 1; for example u(0) = u0 and u(1) = u1 , where u0 and u1 are given real numbers. Example 0.6. (a boundary value problems in Rn :) The Laplace equation in Rn , x = (x1 , x2 , . . . , xn ): 2 2 2 ∆ u = ∂ u + ∂ u + . . . + ∂ u = 0, x ∈ Ω ⊂ Rn , n ∂x21 ∂x22 ∂x2n u(x) = f (x), x ∈ ∂Ω (boundary of Ω). Remark 0.2. In general, in order to obtain a unique solution for a (partial) differential equation, one needs to supply physically adequate boundary data. For instance in Example 0.2 for the potential problem uxx + uyy , stated in a rectangular domain (x, y) ∈ Ω := (0, 1) × (0, 1), to determine a unique 0.2. TRINITIES 13 solution we need to give 2 boundary conditions in the x-direction, i.e. the numerical values for u(0, y) and u(1, y), and another 2 in the y-direction: the numerical values of u(x, 0) and u(x, 1). Whereas to determine a unique solution for the wave equation utt − uxx = 0, it is necessary to supply 2 initial conditions in the time variable t, and 2 boundary conditions in the space variable x. We observe that in order to obtain a unique solution we need to supply the same number of boundary (initial) conditions as the order of the differential equation in each spatial (time) variable. The general rule is that one should supply as many conditions as the highest order of the derivative in each variable. See also Remark 0.1, in the case of unbounded spatial domain. 3. Eigenvalue problems (EVP) Let A be a given matrix. The relation Av = λv, v 6= 0, is a linear equation system where λ is an eigenvalue and v is an eigenvector. In the example below we introduce the corresponding terminology for differential equations. Example 0.7. A differential equation describing a steady state vibrating string is given by −u′′ (x) = λu(x), u(0) = u(π) = 0, where λ is an eigenvalue and u(x) is an eigenfunction. u(0) = 0 and u(π) = 0 are boundary values. The differential equation for a time dependent vibrating string with small displacement u(x, t), which is fixed at the end points, is given by utt (x, t) − uxx (x, t) = 0, 0 < x < π, t > 0, u(0, t) = u(π, t) = 0, t > 0, u(x, 0) = f (x), u (x, 0) = g(x). t Using separation of variables, see also Folland [22], this equation is split into two eigenvalue problems: Insert u(x, t) = X(x)T (t) 6= 0 (a nontrivial solution) into the differential equation to get utt (x, t) − uxx (x, t) = X(x)T ′′ (t) − X ′′ (x)T (t) = 0. (0.2.4) Dividing (0.2.4) by X(x)T (t) 6= 0 separates the variables, viz T ′′ (t) X ′′ (x) = =λ T (t) X(x) (a constant independent of x and t). (0.2.5) 14 CHAPTER 0. INTRODUCTION Consequently we get 2 ordinary differential equations (2 eigenvalue problems): X ′′ (x) = λX(x), T ′′ (t) = λT (t). and (0.2.6) The usual three types of boundary conditions 1. Dirichlet boundary condition Here the solution is known at the boundary of the domain, as u(x, t) = f (x), for x = (x1 , x2 , . . . , xn ) ∈ ∂Ω, t > 0. 2. Neumann boundary condition In this case the derivative of the solution in the direction of the outward normal is given, viz ∂u = n · grad(u) = n · ∇u = f (x), x ∈ ∂Ω ∂n n = n(x) is the outward unit normal to ∂Ω at x ∈ ∂Ω, and grad(u) = ∇u =  ∂u ∂x1 , ∂u ∂u  . ,..., ∂x2 ∂xn 3. Robin boundary condition (a combination of 1 and 2), ∂u + k · u(x, t) = f (x), ∂n k > 0, x = (x1 , x2 , . . . , xn ) ∈ ∂Ω. In homogeneous case, i.e. for f (x) ≡ 0, Robin condition means that the heat ∂u flux ∂n through the boundary is proportional to the temperature (in fact the temperture difference between the inside and outside) at the boundary. p Example 0.8. For u = u(x, y) we have n = (nx , ny ), with |n| = n2x + n2y = 1 and n · ∇u = nx ux + ny uy . Example 0.9. Let u(x, y) = x2 + y 2 . We determine the normal derivative of u in (the assumed normal) direction v = (1, 1). The gradient of u is the vector valued function ∇u = 2x · ex + 2y · ey , where ex = (1, 0) and ey = (0, 1) are the unit orthonormal basis in R2 : ex · ex = ey · ey = 1 and 0.2. TRINITIES 15 y n = (nx , ny ) n2 P Ω n1 x Figure 2: The domain Ω and its outward normal n at a point P ∈ ∂Ω. ex · ey = ey · ex = 0. Note that v = ex + ey = (1, 1) is not a unit vector. The ˆ = v/|v|, i.e. normalized v is obtained viz v ˆ= v ex + ey (1, 1) (1, 1) = √ . =√ |ex + ey | 12 + 12 2 Thus with ∇u(x, y) = 2x · ex + 2y · ey , we get ˆ · ∇u(x, y) = v ex + ey · (2x · ex + 2y · ey ). |ex + ey | which gives (1, 1) 2x + 2y (1, 1) ˆ · ∇u(x, y) = √ · [2x(1, 0) + 2y(0, 1)] = √ · (2x, 2y) = √ . v 2 2 2 The usual three issues I. In theory I1. Existence: there exists at least one solution u. I2. Uniqueness, we have either one solution or no solutions at all. I3. Stability, the solution depends continuously on the data. Note. A property that concerns behavior, such as growth or decay, of perturbations of a solution as time increases is generally called a stability property. 16 CHAPTER 0. INTRODUCTION II. In applications II1. Construction, of the solution. II2. Regularity, how smooth is the found solution. II3. Approximation, when an exact construction of the solution is impossible. Three general approaches to analyzing differential equations 1. Separation of Variables Method. The ...
View Full Document

  • Left Quote Icon

    Student Picture

  • Left Quote Icon

    Student Picture

  • Left Quote Icon

    Student Picture