**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