Submitted at Institute of Computa-tional Mathematics Supervisor O.Univ.-Prof. Dipl.-Ing. Dr. Ulrich Langer Co-Supervisor Dipl.-Ing. Dr. Martin Neum¨uller November 2017 JOHANNES KEPLER UNIVERSITY LINZ Altenbergerstraße 69 4040 Linz, ¨Osterreich www.jku.at DVR 0093696
Space-Time
Finite Element Methods
for Parabolic
Initial-Boundary Problems
Master Thesis
to obtain the academic degree of
Diplom-Ingenieur
in the Master’s Program
In this thesis, we consider the numerical solution of parabolic initial-boundary value problems with variable in space and time, possibly discontinuous coefficients. Such problems typically arise in the simulation of heat conduction problems, diffusion prob-lems, but also for two-dimensional eddy current problems in electromagnetics. Dis-continuous coefficients allow the treatment of moving interfaces like the rotation of an electrical motor. We recall two different approaches to prove that the continuous problem is well-posed in different settings (spaces) under quite general (physical) as-sumptions. In order to solve parabolic problems numerically, a vertical or horizontal method of lines is traditionally applied. However, in this thesis, an alternative ap-proach is chosen. We treat time just as another variable and derive a conforming space-time finite element method. This introduces some challenges, but enables us to apply results from the existing and well investigated theory on elliptic boundary value problems. We show stability of the method, and, additionally, an a priori error esti-mate is provided. The case of local stabilizations, which are important for adaptivity, is also investigated.
To study the method in practice, we introduced typical model problems in one, two, and three spatial dimensions. The implementation of our space-time finite element method is fully parallelized. The numerical studies were performed on the high per-formance computing cluster RADON1, and the outcomes verify the theoretical results.
In dieser Masterarbeit betrachten wir die numerische L¨osung von parabolischen An-fangsrandwertproblemen mit raum- und zeitabh¨angigen Diffusionskoeffizienten, die m¨oglicherweise unstetig sind. Solche Probleme treten zum Beispiel in der Simulation von W¨armeleitungsproblemen, Diffusionsproblemen, aber auch f¨ur zwei-dimensionale Wirbelstromprobleme in der Elektromagnetik, auf. Hier erlauben unstetige Diffusions-koeffizienten den Fall von sich bewegenden Interfaces, zum Beispiel die Drehung eines Elektromotors. Wir zeigen zwei unterschiedliche Zug¨ange zur Untersuchung der Exis-tenz und Eindeutigkeit einer L¨osung f¨ur das stetige Problem unter allgemeinen (phy-sikalischen) Vorraussetzungen. Zur numerische L¨osung parabolischer Probleme wird ¨
ublicherweise eine vertikale oder horizontale Linienmethode verwendet. Wir w¨ahlen jedoch einen alternativen Zugang. Die Zeit wird einfach als eine zus¨atzliche Variable behandelt und wir leiten eine konforme Raum-Zeit Finite Elemente Methode her. Dies f¨uhrt zu neuen Herrausforderungen, erlaubt uns aber die existierende und gut fun-dierte Theorie elliptischer Randwertrobleme zu verwenden. Wir zeigen die Stabilit¨at der Methode, die es uns erlaubt eine a priori Fehlerabsch¨atzung herzuleiten. Der Fall lokaler Stabilisierungen ist inkludiert. Dies ist wichtig in Hinblick auf Adaptivit¨at. F¨ur die praktische Anwendung definieren wir uns typische Modelprobleme in einer, zwei, und drei Raumdimensionen. Die Implementierung unserer Raum-Zeit Finiten Elemente Methode is komplett parallelisiert. Die numerischen Studien wurden am High-Performance-Computing-Cluster RADON1 durchdef¨uhrt und die Resultate best¨ a-tigen unsere theoretischen Aussagen.
I would like express special thanks to my girlfriend Sarah for always pushing me forward to give my best, and my parents Stefan and Brigitte for enabling me to follow my studies.
I would like to thank my supervisors Prof Ulrich Langer and Dr Martin Neum¨uller for introducing me to the topic and giving me the opportunity to write this thesis. Moreover, I would like to thank them for the valuable discussions, for proofreading, and for organizing financial support.
I also thank the Institute of Computational Mathematics for providing me a place to write this thesis.
This work was supported by a DK-fellowship from the Doctoral Program ”Computa-tional Mathematics” (W1214).
Andreas Schafelner Linz, November 2017
1 Introduction 1
2 Preliminaries 4
2.1 Basic Notation . . . 4 2.2 Function spaces . . . 7 2.3 Additional Theorems . . . 10
3 Space-time variational formulations 11
3.1 Formulation in Sobolev spaces on the space-time cylinder . . . 11 3.2 Alternative approach in Bochner spaces of abstract functions . . . 16
4 Space-time Finite Element Methods 18
4.1 A Space-time FEM based on time upwind stabilization . . . 18 4.2 A Space-time FEM based on stabilizing bubble functions . . . 37
5 Numerical Results 39
5.1 1D Examples . . . 41 5.2 2D Examples . . . 49 5.3 3D Examples . . . 57
6 Conclusions & Outlook 59
Introduction
When we deal with physical problems, for instance, diffusion problems, heat-conduction problems, or simulations of electrical machines, the governing partial differential equa-tions (PDEs) are often of parabolic type. Thus, the development of numerical schemes to solve parabolic equations is of great importance. The standard approach for solving parabolic PDEs is usually some kind of time-stepping method, with semi-discretization in the spatial variables. Another approach would be to first discretize with respect to time and then perform a discretization in the spatial variables. This approach is called Rothe’s method. A more recent and alternative approach consists in a full space-time discretization at once by treating time just as another space variable, i.e., we solve a problem with one dimension more. The basic steps for these methods can be summarized in the following way:
1. Line Variational Formulation and Vertical Method of Lines: • multiply the PDE by an appropriate test-function v(x), • integrate over the spatial computational domain Ω,
• use integration by parts on the highest order spatial derivative,
• discretize first in space by some spatial discretization like finite element method (FEM), and then solve the resulting first-order system of ordinary differential equations in time with an appropriate time-stepping method, e.g., a Runge-Kutta method.
2. Line Variational Formulation and Horizontal Method of Lines (Rothe’s method): • multiply by the PDE an appropriate test-function v(x),
• integrate over the spatial computational domain Ω,
• use integration by parts on the highest order spatial derivative,
• discretize first in time by some time-stepping method like the implicit Euler scheme, and then discretize the resulting sequence of elliptic problems by means of an appropriate discretization method like the FEM.
3. Space-time Variational Formulation:
• multiply the PDE by an appropriate test-function v(x, t), • integrate over the space-time domain (cylinder) Q = Ω × (0, T ),
• use integration by parts, e.g. on the highest order spatial derivative and/or the temporal derivative,
• discretize in space and time simultaneously, e.g., by space-time FEM or Isogeometric Analysis (IgA), and solve the resulting linear system by an efficient solver.
In this master thesis, we will focus on the latter approach. The motivation behind this is that, for elliptic problems, there exist plenty of efficient and, most important, parallel solving methods. If we would be able to derive a stable discrete bilinear form, for which we can prove coercivity (ellipticity) in some mesh-dependent norm in the space-time FE-space, then we can solve the space-time problem fully in parallel. An-other reason for the space-time approach is that we are not restricted to a special structure of the mesh. This means that we can apply adaptive mesh refinement both in space and time simultaneously. Last but not least, we can easily deal with moving interfaces and domains, where the coefficients of the PDE and/or the spatial domain Ωt depend on the time as well. Under certain assumptions imposed on the movement,
we can transform the time dependent spatial domain to a fixed spatial domain via a change of variables (see [14, Chapter III, §1]).
The standard discretization techniques, namely the vertical method of lines and Rothe’s method, and their properties are well investigated, see [31] and [15], respectively. How-ever, their sequential structure complicates the parallel solution of the resulting dis-cretized problem, the development of efficient space-time adaptive methods, as well as the treatment of moving interfaces and spatial domains. The application of a space-time finite element scheme has already a long history, see e.g. [9, 12]. How-ever, the analysis of the equivalent operator equations was done more recently, see, e.g. [26, 18, 33]. Another popular approach are time-parallel multigrid methods [10]. Most of the more recent space-time finite element methods use discontinuous Galerkin methods, at least in time, see, e.g., [19, 20, 21, 30], and the references given therein. But also conforming space-time methods have been developed, e.g., Steinbach intro-duced a stable Petrov-Galerkin method [28], and Toulopoulos used bubble functions to stabilise a Galerkin method [32]. In the context of using Isogeometric Analysis as space-time discretization method, Langer, Moore and Neum¨uller [16] proposed a space-time method for parabolic evolution equations.
The main aims of this thesis are first to generalize the results for a space-time scheme proposed by Langer, Moore and Neum¨uller in [16], where the authors use IgA for the discretization, to the case of moving interfaces, i.e., t-dependent, discontinuous coef-ficients, and, second, to allow local, i.e., element-wise, stabilization. The last issue is important for space-time adaptivity. Instead of IgA, we will use a conforming finite element method (FEM) to discretize the parabolic initial-boundary value problem,
which we specify in the following. Then we will shortly describe a second space-time scheme, introduced by Toulopoulos in [32]. Let Q = QT := Ω × (0, T ) be the
space-time cylinder, with Ω ⊂ Rd, d ∈ {1, 2, 3}, being a sufficiently smooth and bounded
spatial domain, and T > 0 being the final time. Furthermore, let Σ := ∂Ω × (0, T ), Σ0 := Ω × {0} and ΣT := Ω × {T } such that ∂Q = Σ ∪ Σ0∪ ΣT. Then we consider the
following model problem that can formally be written as follows: Given f , g, ν and u0, find u such that (s.t.)
∂u
∂t(x, t) − divx(ν(x, t)∇xu(x, t)) =f (x, t), (x, t) ∈ Q, (1.1) u(x, t) =g(x, t) = 0, (x, t) ∈ Σ, (1.2) u(x, 0) =u0(x) = 0, x ∈ Ω, (1.3)
where the diffusion coefficient (reluctivity in electromagnetics) ν is a given uniformly positive and bounded coefficient. The dependence of ν not only on space but also on time enables us to model moving interfaces. Note that we do not require ν to be smooth. In fact, we will admit discontinuities for ν. For simplicity, we assume homogeneous Dirichlet boundary and initial conditions.
The thesis is structured in the following way: In Chapter 2 we will introduce some basic notation and basic mathematical concepts used throughout the thesis. In Chapter 3, we will consider the existence and uniqueness of a weak solution to the parabolic initial boundary value problem (1.1) – (1.3). In Chapter 4, we will derive a stable discrete variational formulation and the space-time finite element scheme. Moreover, we will prove an a priori error estimate. Additionally, we shortly describe the space-time finite element scheme with stabilizing bubble functions by Toulopoulos. In Chapter 5, we discuss the numerical results. In Chapter 6, we draw some conclusions, and provide an outlook to future work.
Preliminaries
In this chapter, we introduce some definitions and notations, which we will use through-out this thesis. We start with the basic notation.
2.1
Basic Notation
Notation. We denote the Euclidean inner product on Rd by
x · y :=
d
X
i=1
xiyi,
for x, y ∈ Rd. We define the Euclidean standard norm by
|x| = d X i=1 |xi|2 1/2 , for x ∈ Rd.
Notation (Kronecker symbol).
δij =
(
1, for i = j, 0, else.
Notation (Multiindex). Let α = (α1, . . . , αd+1) be a vector. If all components αi are
non-negative integers, we call α a multiindex of order |α| = α1+ · · · + αd+1.
Notation (Partial derivatives). Let u : Rd+1 → R, and α be a multiindex. Then we define Dαu = ∂ |α|u ∂α1 x1 · · · ∂ αd+1 xd+1 . 4
We denote the full gradient of a sufficiently smooth scalar function u : Rd+1→ R by ∇u := (∇xu, ∂tu)>, (2.1) where ∇xu := ∂u ∂x1 , . . . , ∂u ∂xd > , (2.2) ∂tu := ∂u ∂xd+1 . (2.3)
Moreover, we denote the (spatial) divergence of a sufficiently smooth vector field u : Rd+1 → Rd by div u = ∂u1 ∂x1 + · · · + ∂ud+1 ∂xd+1 , (2.4) divxu = ∂u1 ∂x1 + · · · + ∂ud ∂xd . (2.5)
Note that we treat the time t as just another variable xd+1.
Remark 2.1. We denote the derivative of a one-dimensional function u : R → R by
u0(x) := d dxu(x).
Notation. We will use the following spaces of (scalar) continuous functions: C(Ω) := {v : Ω → R : v is continuous on Ω},
Ck(Ω) := {v ∈ C(Ω) : Dαv ∈ C(Ω), for all |α| ≤ k}, C∞(Ω) := {v ∈ C(Ω) : v is infinitely differentiable }, Cc∞(Ω) := {v ∈ C∞(Ω) : v has compact support}, Ck,m(Q) := {v ∈ C(Ω) : D(α1,...,αd,0)v ∈ C(Ω), |(α
1, . . . , αd, 0)| ≤ k
∧ ∂tjv ∈ C(Ω), j ≤ m}
for k, m ∈ N0. We identify C0(Ω) = C(Ω). For vector fields, e.g., u : Rd+1 → Rd, we
write
[C(Ω)]d:= {v : Ω → Rd: vi ∈ C(Ω), i = 1, . . . , d}.
With this notations, we can now formulate smoothness conditions for a solution to the model problem (1.1) – (1.3). If we require the right hand side f to be continuous and the diffusion coefficient ν to be continuously differentiable in Q, then the exact solution
is called a classical solution. However, these are very strong restrictions on both the right hand side and the diffusion coefficient. So it is naturally to ask if we can relax the smoothness conditions on the given functions f and ν, and how does this affect the solution. This motivates the introduction of weak (or generalized) derivatives, which are formally defined as follows.
Definition 2.2 (Weak (or generalized) derivative). Let U ⊂ R open and u, w ∈ L1,loc(U ). If they satisfy the integral identity
Z U w(x) ϕ(x) dx = − Z U u ϕ0(x) dx, for all ϕ ∈ Cc∞(U ), (2.6)
then we call w the weak (or generalized) derivative of u. Here,
L1,loc := {u : U → R : v ∈ L1(V ) for each V ⊂⊂ U }.
One can show that this weak derivative is uniquely defined up to a set of measure zero (see [8]).
Remark 2.3. For functions u : Rd+1 → R, we can define the weak gradient w :
Rd+1 → Rd+1 in a similar manner to Definition 2.2, i.e., if Z U w(x, t) · ϕ(x, t) d(x, t) = − Z U u(x, t) div ϕ(x, t) d(x, t), ∀ϕ ∈ [C∞ c (U )] d+1 (2.7)
then we call w the weak gradient of u. Similarly, if w : Rd+1 → Rd satisfies
Z U w(x, t) · ψ(x, t) d(x, t) = − Z U u(x, t) divxψ(x, t) d(x, t), ∀ψ ∈ [Cc∞(U )] d, (2.8)
then we call w the weak spatial gradient and if w : Rd+1 → R Z U w(x, t) η(x, t) d(x, t) = − Z U u(x, t) ∂tη(x, t) d(x, t), ∀η ∈ Cc∞(U ), (2.9)
then w is called the weak temporal derivative.
Remark 2.4. If w is the weak derivative of u in the sense of Definition 2.2 or Re-mark 2.3, then we use the notation
u0 = w, or ∇u = w, respectively.
2.2
Function spaces
In this section, we will briefly present Banach and Hilbert spaces, which we need throughout the thesis. For further analysis of such spaces and their properties, we refer to e.g. [6], and the references given therein. We start with the Lebesgue spaces of order p, which we introduce in the following definition. In the following and throughout this thesis, Ω ⊂ Rn is always a bounded Lipschitz domain.
Definition 2.5. Let Ω ⊂ Rn. For 1 ≤ p < ∞, the space Lp(Ω) consists of all
measurable functions u, such that Z
Ω
|u|p dx < ∞.
A measurable function u belongs to the space L∞(Ω), if it is essentially bounded, i.e.,
ess sup
x∈Ω
|u(x)| < ∞.
If we equip the Lebesgue spaces Lp(Ω) with the norm
kukLp(Ω) := ( R
Ω|u|
p dx1/p
, for 1 ≤ p < ∞, ess supx∈Ω|u(x)| < ∞, for p = ∞ , one can show the following result.
Theorem 2.6. For 1 ≤ p ≤ ∞, the Lebesgue spaces Lp(Ω) are Banach spaces with the
norm k . kLp(Ω). Moreover, the space L2(Ω) is a Hilbert space equipped with the inner product (u, v) = (u, v)L2(Ω):= Z Ω u(x)¯v(x) dx, ∀u, v ∈ L2(Ω). (2.10) Proof. See [6].
Remark 2.7. Two functions f and g are said to be identical in Lp(Ω), if kf −gkLp(Ω) = 0. This means that if two functions f and g differ on a zero set N ⊂ Ω, we still consider them identical in Lp(Ω), as they coincide almost everywhere (a.e.). Moreover,
we deduce that a function in Lp(Ω) must only be defined a.e., i.e., in general, point
evaluation is not well defined.
We can now combine the weak derivatives with these Lebesgue spaces, leading us to Sobolev spaces of weak derivatives. We again start with a formal definition:
Definition 2.8 (Wpk(Ω)). The Sobolev space Wpk(Ω) consists of all functions u, whose weak derivatives up to index k belong to the space Lp(Ω), i.e.,
Wpk(Ω) := {u ∈ Lp(Ω) : Dαu ∈ Lp(Ω), for all |α| ≤ k},
On Wk
p(Ω), we can define the norm
kukWk p(Ω) := X |α|≤k kDαuk2 Lp(Ω) 1/2 .
Remark 2.9. There is an alternative way to introduce Sobolev spaces by using dense subspaces. Let us denote the closure of C∞(Ω) ∩ Wpk(Ω) wrt the norm k . kWk
p(Ω) by Hk p(Ω), i.e., Hpk(Ω) := C∞(Ω) ∩ Wk p(Ω) k . k W kp (Ω). Meyers and Serrin showed in [17], that
Wpk(Ω) = Hpk(Ω), i.e., the space C∞(Ω) ∩ Wk
p(Ω) is dense in Wpk(Ω).
Theorem 2.10. For each k = 1, . . . and 1 ≤ p ≤ ∞, the Sobolev space Wk
p(Ω) is a
Banach space. Proof. [1, 6, 8].
Remark 2.11. For the special case p = 2, we will use the notation Hk(Ω) := W2k(Ω),
as these spaces are Hilbert spaces [1, 8], equipped with the inner product (u, v)Hk(Ω) :=
X
|α|≤k
(Dαu, Dαv)L2(Ω),
and the norm
kukk := kukHk(Ω) := q
(u, u)Hk(Ω). Moreover, we will use the semi-norm
|u|Hk(Ω) := s X |α|=k (Dαu, Dαu) L2(Ω).
Definition 2.12. We denote the closure of Cc∞(Ω) in Wk
p(Ω) by Wp,0k (Ω), i.e., Wp,0k (Ω) := C∞ c (Ω) k . k W kp (Ω) ,
for 1 ≤ p < ∞ and k ∈ N. Moreover, for p = 2, we use the notation H0k(Ω) := W2,0k (Ω).
This definition may seem rather abstract, however, we will give an equivalent charac-terisation of a function u ∈ Wp,0k (Ω) in the following.
If we want to apply the theory above in the context of boundary value problems, we have consider expressions on the boundary ∂Ω of Ω. But in general, as already mentioned, functions in Sobolev spaces only have to be defined a.e.. Moreover, as the boundary ∂Ω is a zero set, two functions would be considered equal even for different boundary data. However, using the concept of traces, we can again prescribe values on the boundary ∂Ω.
Theorem 2.13 (Trace theorem [8]). Let Ω be bounded with sufficiently smooth ∂Ω. Then there exists a bounded linear operator
T : Wp1(Ω) → Lp(∂Ω)
such that
(i) T u = u|∂Ω, for u ∈ Wp1(Ω) ∩ C(Ω), and
(ii)
kT ukLp(∂Ω) ≤ CkukWp1(Ω), for all u ∈ Wp1(Ω), where C does not depend on u.
Proof. For the existence, see [6, 8], for the norm inequality see [8]. Remark 2.14. For p = 2, we denote
T (H1(Ω)) = H1/2(∂Ω) ⊂ L2(∂Ω).
For details on Sobolev spaces Wk
p(Ω) with non-integer k, we refer to [1, 6, 8].
Definition 2.15. We call T u the trace of u on ∂Ω.
The above theorem gives only information about the traces of functions belonging to W1
p(Ω). However, we might be interested in traces of functions with higher order
weak derivatives. The following theorem gives information about traces of functions belonging to Wpk(Ω).
Theorem 2.16 ([6]). Let Ω be a Lipschitz domain and m ∈ N, 1 ≤ p ≤ ∞. Then Wk
p,0(Ω) consists off all functions u ∈ Wpk(Ω) with T Dαu = 0, for |α| ≤ k − 1.
Proof. See [6].
Remark 2.17. For simplicity, we will not always write the trace operator for expres-sions on the boundary but instead use the same notation as for smooth functions, i.e., u|∂Ω = T u for u ∈ Wp1(Ω).
With these spaces and the definition of traces, we can now formally define special function spaces on the whole space-time cylinder Q (see [14]), which we will need in the analysis of our model problem (1.1) – (1.3).
Definition 2.18. We define the following Sobolev (Hilbert) spaces
H01(QT) = W2,01 (QT) := {u ∈ L2(QT) : ∇u ∈ L2(QT) ∧ u|Σ = 0},
H1,0(QT) = W21,0(QT) := {u ∈ L2(QT) : ∇xu ∈ L2(QT)},
˚
H1,0(QT) = ˚W21,0(QT) := {u ∈ H1,0(QT) : u|Σ= 0},
equipped with the usual inner products and norms, as well as the Banach space V2(QT) := {u ∈ H1,0(QT) : |u|QT < ∞},
with subspaces ˚
V2(QT) := {u ∈ ˚H1,0(QT) : |u|QT < ∞}, V21,0(QT) := {u ∈ V2(QT) : lim
∆t→0ku(·, t + ∆t) − u(·, t)kL2(Ω) = 0, uniformly on [0, T ]},
˚
V21,0(QT) := V21,0(QT) ∩ ˚H1,0(QT),
where the norm | . |QT is defined by |u|Qt := max
0≤τ ≤tku( . , τ )kL2(Ω)+ k∇xukQt. (2.11)
We will make use of the above theory in Chapter 3, where we will introduce a new class of solutions to the model problem (1.1) – (1.3), the so called weak (or generalized) solutions.
2.3
Additional Theorems
When we consider operator equations, e.g., Bu = g,
we want to know whether this equation has a unique solution. The following theorem provides conditions under which there is indeed a unique solution.
Theorem 2.19 ([27, Thm. 3.7]). Let X and Π be Hilbert spaces and let B : X → Π0 be a bounded and linear operator. Further we assume the stability condition
cSkvkX ≤ sup 06=q∈Π
hBv, qi kqkΠ
for all v ∈ (ker B)⊥.
For a given g ∈ ImX(B) there exists a unique solution u ∈ (ker B)⊥ of the operator
equation Bu = g satisfying kukX ≤ 1 cS kgkΠ0. (2.12) Proof. See [27].
Space-time variational formulations
In this chapter, we will focus on deriving a variational formulations for our model problem (1.1) – (1.3). We want to show the existence and uniqueness of a solution, and, moreover, to what class of solutions it belongs. There are different ways to obtain existence and uniqueness results, see e.g. [14, 28]. We will focus on the techniques presented by Ladyˇzhenskaya in [14].
3.1
Formulation in Sobolev spaces on the
space-time cylinder
Now let us reconsider the model problem (1.1) – (1.3) : Find u s.t.
Mu ≡ ∂tu − div(ν∇xu) = f in QT, (3.1)
u = 0 on Σ, u = ϕ on Σ0, (3.2)
with given data
ϕ ∈ L2(Ω) and f ∈ L2,1(QT) := {v :
Z T
0
kf ( . , t)kL2(Ω) dt < ∞}, (3.3)
and a uniformly bounded coefficient
ν ≤ ν(x, t) ≤ ν, for almost all (x, t) ∈ QT, (3.4)
where ν, ν = const. > 0. To show now the existence of a weak solution in an appro-priate function space, we use Galerkin’s method. We formally start with multiply-ing the PDE by the solution u and integrate over the truncated space-time domain Qt= Ω × (0, t), t ∈ (0, T ), i.e., Z Qt Mu · u dxdt = Z Qt f u dxdt. (3.5) 11
Using integration by parts, the homogeneous boundary condition on the lateral bound-ary Σ and nx = 0 on Σ0∪ Σt, where Σt:= Ω × {t}, we obtain
Z Qt Mu · u dxdt = Z Qt ∂tuu − divx(ν(x, t)∇xu) dxdt = Z Qt 1 2∂t(u 2) dxdt + Z Qt ν(x, t)|∇xu|2 dxdt,
for the left hand side of (3.5). Now we use Gauss’ theorem and the fact that nt ≡ 0
on Σ to get rid of the time derivative and we obtain the following identity 1 2ku( . , t)k 2 L2(Ω)+ Z Qt ν(x, t)|∇xu| dxdt = 1 2ku( . , 0)k 2 L2(Ω)+ Z Qt f u dxdt. (3.6)
We call (3.6) energy balance equation. From this equation, we will derive a bound for u in some specific norm | . |Qt. First, we estimate the left hand side (lhs) of (3.6) from below by 1 2ku( . , t)k 2 L2(Ω)+ Z Qt ν(x, t)|∇xu|2 dxdt ≥ 1 2ku( . , t)k 2 L2(Ω)+ ν Z Qt |∇xu|2 dxdt,
and the right hand side (rhs) of (3.6) from above by 1 2ku( . , 0)k 2 L2(Ω)+ Z Qt f u dxdt ≤ 1 2ku( . , 0)k 2 L2(Ω)+ Z t 0 Z Ω f (x, τ )u(x, τ ) dx dτ ≤ 1 2ku( . , 0)k 2 L2(Ω)+ Z t 0 kf (·, τ )kL2(Ω)ku(·, τ )kL2(Ω) dτ ≤ 1 2ku( . , 0)k 2 L2(Ω)+ Z t 0 kf (·, τ )kL2(Ω) max σ∈[0,t]ku(·, σ)kL2(Ω)dτ ≤ 1 2ku( . , 0)k 2 L2(Ω)+ kf k2,1,Qt max τ ∈[0,t]ku(·, τ )kL2(Ω).
Combining these two estimates gives us 1 2ku( . , t)k 2 L2(Ω)+ νk∇xuk 2 Qt ≤ ku( . , 0)k 2 L2(Ω)+ max0≤τ ≤tku(·, τ )kL2(Ω)kf k2,1,Qt. (3.7) Denoting max0≤τ ≤tku(·, τ )kL2(Ω) by y(t) and multiplying (3.7) by 2, we obtain
ku( . , t)k2
L2(Ω)+ 2νk∇xuk
2
Qt ≤ y(t)ku( . , 0)kL2(Ω)+ 2y(t)kf k2,1,Qt = j(t), where we used the estimate ku( . , 0)k2L
2(Ω) ≤ max0≤τ ≤tku(·, τ )kL2(Ω)ku( . , 0)kL2(Ω). From this, we deduce two inequalities, i.e,
y(t)2 ≤ j(t) and k∇xuk2Qt ≤ (2ν)
−1
The second estimate can easily be verified, whereas the first one is obtained by esti-mating the lhs from below by ku( . , t)k2L
2(Ω). This expression holds for any τ ∈ [0, t], hence it holds also for the maximum. However, the only terms in j(t) depending on t are y(t), where we already take a maximum over [0, t]. Thus, the first expression of (3.8) follows. We take the square-root of both expressions and add them up to obtain
|u|Qt = y(t) + k∇xukQt ≤ (1 + 1 2ν) −1/2|u|1/2 Qt(ku( . , 0)kL2(Ω)+ 2kf k2,1,Qt) 1/2 . (3.9)
We bring similar terms on the same side and take the square on each side of the inequality. Thus we have obtained an upper bound for |u|Qt in the form
|u|Qt ≤ (1 + 1 2ν)
−1
(ku( . , 0)kL2(Ω)+ 2kf k2,1,Qt) =: cF (t), (3.10)
which holds for any t ∈ [0, T ]. However, this bound requires a solution where point evaluation with respect to (wrt) time is well defined. Before we proof that our problem (3.1) has such a weak solution, we have to introduce a suitable definition of the weak solution.
Definition 3.1. A function u ∈ ˚H1,0(Q
T) is called a generalized (weak) solution in
H1,0(Q
T) of the parabolic initial-boundary value problem (3.1) – (3.2) if it satisfies the
identity M(u, v) ≡ Z QT −u∂tv + ν(x, t)∇xu∇xv dxdt = Z Ω ϕv( . , 0) dx + Z QT f v dxdt, (3.11) for all v ∈ ˆH1 0(QT) := {v ∈ H01(QT) : v = 0 on ΣT}.
To proof solvability of (3.1) in this class, i.e. solvability of (3.11), we will use Galerkin’s method. Let {ϕj} be a L2-orthonormal fundamental system in ˚W21(Ω). In (3.1), we
substitute u with an appropriate test function uN, multiply the obtained equation by each ϕj for j = 1, . . . , N and integrate wrt x over Ω. We use integration by parts in
the principle term, and obtain a system of N equations
(∂tuN, ϕj) + (ν( . , t)∇xu, ∇xϕj) = (f, ϕj), (3.12)
where ( . , . ) = ( . , . )L2(Ω) is the standard L2(Ω) inner product. In (3.12), we express uN with the fundamental system {ϕj}, i.e., uN(x, t) := PNj=1cNj (t)ϕj(x). We can
rewrite (3.12) wrt the coefficient functions ci(t) = cNi (t), N X j=1 d dtcj(t) (ϕ| j{z, ϕi}) δi j + N X j=1 cj(t)(ν( . , t)∇xϕj, ∇xϕi) = (f ( . , t), ϕi), for i = 1, ..., N, (3.13)
with the initial condition
N
X
j=1
cj(0)(ϕj, ϕi) = (ϕ, ϕi), (3.14)
that is nothing but the L2-projection to span{ϕ1, . . . , ϕN}. This system is a system
of N linear ordinary differential equations with principal terms dtdci(t) and bounded
coefficient functions in front of the zero-order terms ci(t). This system has a unique
solution of absolutely continuous functions cN
1 (t), . . . , cNN(t) (see [11, IX, 60.2]), which
uniquely define the approximate solution uN. Now, we want to derive a bound for the
series uN. To do so, we first multiply each equation in (3.12) with the corresponding solution coefficient cN
l (t), l = 1, . . . , N and sum these equations up from 1 to N . We
proceed by integrating this sum over (0, t) and obtain an equation of the form (3.6) for uN, i.e., 1 2ku N( . , t)k2 L2(Ω)+ Z Qt ν(x, t)|∇xuN| dxdt = 1 2ku N( . , 0)k2 L2(Ω)+ Z Qt f uN dxdt. (3.15)
We can derive a bound for |uN|
QT from (3.15) in the same manner as we did for |u|QT from (3.6), i.e.,
|uN|Qt ≤ c ku
N
( . , 0)kL2(Ω)+ 2kf k2,1,Qt. Furthermore, we know the upper bound kuN( . , 0)k
L2(Ω) ≤ kϕkL2(Ω). Therefore, we obtain the bound
|uN|QT ≤ ˜c, (3.16)
where ˜c is a constant independent of N . Hence the sequence {uN} is a bounded se-quence in the Hilbert space L2(QT). This can easily be deduced from (3.16) and
Def-inition 2.18. Hilbert spaces are reflexive spaces, see e.g. [6]. Thus, {uN} has a weakly
convergent subsequence {uNk}. The same holds true for its derivatives {∇
xuN} with
{∇xuNk}. Therefore, {uNk} and {∇xuNk} converge weakly to some unique element
u ∈ ˚H1,0(Q
T). Is this u the desired generalized (weak) solution of our model problem
(3.1) - (3.2)? Let us again multiply (3.12) by some arbitrary absolutely continuous functions dl(t) with dtddl ∈ L2(0, T ), with dl(T ) = 0. We sum the obtained equations
up from 1 to N , integrate over the interval (0, T ) and perform integration by parts with respect to time. The resulting equation is
Z QT −uN∂ tΦ + ν(x, t)∇xuN∇xΦ dxdt = Z Ω uNΦ|t=0dx + Z QT f Φ dxdt, (3.17) for Φ(x, t) = PN
k=1dk(t)ϕk(x). The set of all such functions Φ with the desired
properties of dl is denoted by MN. The superset
S∞
p=1Mp is dense in ˆH01(QT) (see
[14]). We fix a Φ ∈ Mp and take the limit of (3.17) for Nk≥ p, i.e.,
Z QT − uNk∂ tΦ | {z } →u ∂tΦ +ν(x, t) ∇xuN∇xΦ | {z } →∇xu∇xΦ dxdt = Z Ω uNΦ|t=0 | {z } →ϕΦ|t=0 dx + Z QT f Φ dxdt. (3.18)
We obtain exactly the definition of a generalized (weak) solution (3.11), with v = Φ ∈ M. As these union of all such spaces is dense in ˆH01(QT), the equation (3.18) holds for
any v ∈ ˆH1
0(QT). Thus u is indeed a generalized solution of our model problem (3.1)
– (3.2). We gather these results in the following theorem.
Theorem 3.2 ([14, Chapter III, Thm. 3.1]). Under the conditions (3.3) and (3.4), the problem (3.1) – (3.2) has at least one generalized (weak) solution in ˚H1,0(Q
T), as
defined in Definition 3.1.
Proof. Follows from the derivation above.
We know now that at least one solution u exists, but is this solution unique? To prove this, we will make again use of the results presented by Ladyˇzhenskaya in [14, Chapter III,§2]. First, we consider our generalized solution u as a generalized solution in L2(QT) of the problem
∂tu − ∆u = ˜f + divx(F ) in QT, (3.19)
u|t=0 = φ(x) for x ∈ Ω, u|Σ = 0, (3.20)
with ˜f ≡ f and Fi = ν(x, t)∇xu − ∇xu. Hence, by [14, Chapter III, Thm. 2.2 & Thm.
2.3], it follows that u(x, t) is a generalized solution of (3.19) – (3.20) in ˚V21,0(QT). By
this, we can define a new class of generalized solutions.
Definition 3.3 ([14, Chapter III]). A generalized solution u ∈ ˚H1,0 of (3.1) - (3.2)
is a called a generalized solution of (3.1) - (3.2) in ˚V21,0(QT), if u ∈ ˚V 1,0
2 (QT) and it
fulfils the energy-balance equation (3.6) and the identity Z Ω u(x, t)v(x, t) dx − Z Ω ϕv(x, 0) dx + Z Qt −u∂tv + ν∇xu∇xv dxdt = Z Qt f v dxdt, (3.21)
for all v ∈ H01(QT) and any t ∈ (0, T ).
We will show uniqueness of the problem (3.1) – (3.2) in H1,0(Q
T) as usual by
contra-diction. Let u1 6= u2 ∈ ˚H1,0(QT) be two generalized solutions of (3.1) – (3.2), then
the difference u := u1 − u2 is also a generalized solution of (3.1) - (3.2), but with
homogeneous initial data and zero right hand side. Moreover, by what we have shown above, it is also a generalized solution in ˚V21,0(QT), so it satisfies (3.6) with zero right
hand side. If it satisfies (3.6), we have shown that its norm |u|QT is subject to the bound (3.10), but also with zero right hand side. We obtain u = u1− u2 ≡ 0, which is
a contradiction to our assumption u1 6= u2. Moreover, the operator B, which assigns
each tuple (f, ϕ) its generalized solution in ˚V21,0(QT) is linear and the energy balance
equation (3.6) can be obtained from the identity (3.21) (see [14]). We can summarise the results in the following theorem.
Theorem 3.4 ([14, Chapter III, Thm. 3.2]). If the assumptions (3.3) and (3.4) are fulfilled, then any generalized solution of (3.1) – (3.2) in ˚H1,0(QT) is the generalized
solution in ˚V21,0(QT) and it is unique in ˚H1,0(QT).
Corollary 3.5. If the assumptions (3.3) and (3.4) hold, then there exists a unique generalized solution u ∈ ˚H1,0(QT) ∩ ˚V
1,0
2 (QT) to the problem (3.1) – (3.2).
3.2
Alternative approach in Bochner spaces of
ab-stract functions
In the previous section we ensure the existence and uniqueness of a solution of (1.1) – (1.3) using Sobolev spaces defined on the whole space-time cylinder Q. An alternative approach makes use of Bochner spaces of abstract functions, see, e.g. [18, 26, 28]. In this section, we will briefly summarise the results of Steinbach in [28]. In order to do that, we first need to recall the definition of Bochner spaces.
Definition 3.6 (Bochner space). Let X be a Banach space with norm k . k. Then we define the spaces
C(0, T ; X) := {u : [0, T ] → X : u is continuous and kukC(0,T ;X) < ∞}, (3.22)
Lp(0, T ; X) := {u : (0, T ) → X : u is measurable and kukLp(0,T ;X) < ∞}, (3.23) of abstract functions mapping [0, T ] into some Banach space X, with the respective norms kukC(0,T ;X) := max 0≤t≤Tku(t)k, (3.24) kukLp(0,T ;X) := Z T 0 ku(τ )kp dt 1/p . (3.25)
For further details on Bochner spaces, see e.g. [8, 34].
Let us consider once more the model problem (1.1) – (1.3), i.e. find u s.t. ∂tu − divx(ν∇xu) = f in Q,
u = 0 on Σ, u = u0 on Σ0
for given data f and u0. We derive the weak formulation in the usual manner and
obtain the variational problem: find u ∈ L2(0, T ; H01(Ω)) ∩ H1(0, T ; H
−1(Ω)), with u( . , 0) = u0(·), s.t. a(u, v) := Z T 0 Z Ω ∂tu v + ν∇xu · ∇xv dx dt = Z T 0 Z Ω f v dx dt =: l(v), (3.26)
for all v ∈ L2(0, T ; H01(Ω)), where ν is as before, f ∈ L2(0, T ; H−1(Ω)) and u( . , 0) ∈
H01(Ω). We can not show any form of ellipticity for this Petrov-Galerkin formula-tion. Hence, Steinbach used a stability or inf-sup-condition to guarantee existence and uniqueness of an weak solution. The stability condition uses the corresponding energy norm kuk2 L2(0,T ;H01(Ω))∩H1(0,T ;H−1(Ω)) = k∂tuk 2 L2(0,T ;H−1(Ω))+ kuk 2 L2(0,T ;H01(Ω)). Theorem 3.7 ([28, Thm. 2.1]). For all u ∈ L2(0, T ; H01(Ω)) ∩ H1(0, T ; H
−1(Ω)), with
u( . , 0) = 0 in Ω, there holds the stability condition 1 2√2kukL2(0,T ;H 1 0(Ω))∩H1(0,T ;H−1(Ω)) ≤ sup 06=v∈L2(0,T ;H10(Ω)) a(u, v) kvkL2(0,T ;H1 0(Ω)) . (3.27)
Proof. The proof makes use of the so called quasi-static elliptic Dirichlet boundary value problem, for further details see [28].
The initial condition can be seen as an Dirichlet datum for the space-time domain Q. Hence,vim we can introduce the splitting
u( . , . ) = ¯u( . , . ) + ¯u0( . , . ) in Q, (3.28)
where ¯u0 ∈ L2(0, T ; H01(Ω)) ∩ H1(0, T ; H−1(Ω)) is an extension of the initial datum
u0 ∈ H01(Ω).
Now we homogenise the variational problem (3.26), i.e. find u ∈ X s.t.
a(u, v) = hf, viQ− a(u0, v) for all v ∈ Y, (3.29)
where u0 is a given extension of the initial datum u ∈ H01(Ω) and f ∈ L2(0, T ; H−1(Ω))
is a given right hand side. We use the spaces X := {u ∈ L2(0, T ; H01(Ω)) ∩ H
1
(0, T ; H−1(Ω)) : u( . , 0) = 0 in Ω}, Y := L2(0, T ; H01(Ω)).
The unique solvability of (3.29) follows now from Theorem 2.19. Corollary 3.8 ([28, Col. 2.2]). Let ¯u0 ∈ L2(0, T ; H01(Ω)) ∩ H1(0, T ; H
−1(Ω)) be some
given extension of the initial datum u0 ∈ H01(Ω) and assume f ∈ L2(0, T ; H−1(Ω)).
Then the bilinear form a( . , . ) as given in (3.26) is bounded and satisfies the stability condition (3.27). Hence there exists a unique solution ¯u ∈ X of the variational problem (3.29) satisfying the a priori estimate
k¯ukL2(0,T ;H1 0(Ω))∩H1(0,T ;H−1(Ω))≤ 2 √ 2kf kL2(0,T ;H−1(Ω)) +√2k¯u0kL2(0,T ;H01(Ω))∩H1(0,T ;H−1(Ω)). (3.30)
Space-time Finite Element
Methods
4.1
A Space-time Finite Element Method based on
time upwind stabilization
From the previous chapter, we know that there exists a unique generalized solution of the initial-boundary value problem (1.1) – (1.3) in ˚H1,0(Q) ∩ ˚V21,0(Q). The goal of this chapter is to derive a stable space-time finite element scheme with a coercive (elliptic) discrete bilinear form, and, therefore, to ensure existence and uniqueness of a finite element solution. Similar to Langer et.al. in [16], we use special time-upwind test functions that are locally scaled in our case. First, we need a regular triangulation Th of our space-time domain Q (for details, see e.g. [2, 4]). We formally define this
triangulation as
Th := {E : E ⊂ Q and E open} (4.1)
with the properties
Q = [
E∈Th
E and E ∩ E0 = ∅ for all E0 ∈ Th with E 6= E0. (4.2)
On each of these elements E, we define individual time upwind test functions
vh,t(x, t) := vh(x, t) + θEhE∂tvh(x, t), for all (x, t) ∈ E, (4.3)
where θE is a positive parameter that will be defined later, and hE := diam(E).
Here, vh is some test function from a standard conform finite element space V0h, e.g.,
V0h = {v ∈ C(Q) : v|E ∈ Pp ⊂ Qp} that is considered in this thesis. From now
on, unless specified otherwise, all functions depend on both space and time variables. For simplicity, we can omit the arguments. In this chapter, we will make use of the
following spaces: V0 = H0,01,1(Q) := {u ∈ L2(Q) : ∇xu ∈ L2(Q), ∂tu ∈ L2(Q) and u|Σ∪Σ0 = 0}, (4.4) H0,02,1(Th) := {v ∈ H 1,1 0,0(Q) : v|E ∈ H2,1(E), ∀E ∈ Th}, (4.5) W∞1(Th) := {v ∈ L∞(Q) : v|E ∈ W∞1(E), ∀E ∈ Th}. (4.6) We assume that ν ∈ W1
∞(Th) and that the PDE has a sufficiently smooth solution
u, e.g., u ∈ H0,02,1(Th). Then we proceed in the usual manner, i.e., we first multiply
the PDE (1.1) by our space-time test function vh,t, and then integrate over a single
element E, obtaining Z E (∂tu − divx(ν∇xu))vh,t d(x, t) = Z E f vh,t d(x, t).
Summing up over all elements and applying integration by parts on the principle term, we obtain X E∈Th Z E ∂tu vh,t+ ν∇xu∇xvh,t d(x, t) − Z ∂E ν∇xu · nxvh,t ds(x,t) = X E∈Th Z E
∂tuvh+ θEhE∂tu∂tvh+ ν∇xu∇xvh+ θEhEν∇xu∇x(∂tvh) d(x, t)
− Z
∂E
ν∇xu · nxvh+ θEhEν∇xu · nx∂tvh ds(x,t),
for the left hand side, while the right hand side remains unchanged. For the exact solution u of (1.1) – (1.3), we know that the fluxes have to be continuous, i.e., let E and E0 be two adjacent elements, then
(ν∇xu · nx)|E = (ν∇xu · nx)|E0. (4.7) From this, we know that one part of the boundary terms vanishes from all inner edges, i.e. we obtain
X
E∈Th Z
E
∂tu vh+ θEhE∂tu∂tvh+ ν∇xu∇xvh+ θEhEν∇xu∇x(∂tvh) d(x, t)
−X E∈Th ∂E∩∂Q6=∅ Z ∂E ν∇xu · nxvh ds(x,t) − X E∈Th Z ∂E νθEhE∇xu · nx∂tvh ds(x,t) = X E∈Th Z E f vh,t d(x, t).
We require vh to be zero on Σ, and know that nx vanishes on Σ0 and ΣT. Therefore,
the first boundary term completely disappears from our equation, and we obtain X
E∈Th Z
E
[∂tuvh+ θEhE∂tu∂tvh+ ν∇xu∇xvh+ θEhEν∇xu∇x(∂tvh)] d(x, t)
− X E∈Th Z ∂E νθEhE∇xu · nx∂tvh ds(x,t) = X E∈Th Z E f (vh+ θEhE∂tvh) d(x, t).
We now arrived at the consistency identity for (1.1)
ah(u, vh) = lh(vh), ∀vh ∈ V0h, (4.8)
that holds for a sufficiently smooth solution u, e.g., u ∈ H0,02,1(Th), where
ah(u, vh) := X E∈Th Z E ∂tu vh + θEhE∂tu ∂tvh d(x, t) + Z E ν∇xu · ∇xvh+ θEhEν∇xu · ∇x(∂tvh) d(x, t) (4.9) − Z ∂E θEhEν∇xu · nx∂tvh ds(x,t), lh(vh) := X E∈Th Z E f (vh+ θEhE∂tvh) d(x, t), (4.10) with given ν ∈ W1 ∞(Th) and f ∈ L2(Q).
Remark 4.1. We can derive an equivalent scheme to (4.9). In particular, we perform the same steps as above, but instead of applying integration by parts on both principal terms, we only apply it to the first principal term and keep the second. Hence we obtain another consistency identity for (1.1)
˜
ah(u, vh) = lh(vh), ∀vh ∈ V0h,
that holds for a sufficiently smooth solution u, e.g., u ∈ H0,02,1(Th), where
˜ ah(u, vh) := X E∈Th Z E ∂tu vh+ θEhE∂tu ∂tvh d(x, t) + Z E ν∇xu · ∇xvh+ θEhE divx(ν∇xu)∂tvh d(x, t) with given ν ∈ W1 ∞(Th) and f ∈ L2(Q), and lh as in (4.10).
Remark 4.2. If the test functions vh ∈ V0h are continuous and piecewise linear (p =
1), then the term in (4.9) containing ∇x(∂tvh) vanishes in all elements E ∈ Th, since
it only contains mixed second order derivatives.
Now we look for a Galerkin approximation uh ∈ V0h to the generalized solution u of
our initial boundary value problem (1.1) – (1.3) using the variational identity (4.8), i.e., find uh ∈ V0h such that
ah(uh, vh) = lh(vh), ∀vh ∈ V0h, (4.11)
with ah and lh as defined above by (4.9) and (4.10), respectively. In Chapter 3, we
value problem (1.1) – (1.3). However, our discrete variational problem (4.11) is of a different form. Thus, we have to investigate the stability of the space-time finite element scheme. More precisely, we will even show ellipticity of the bilinear form ah( . , . ) : V0h× V0h→ R wrt the mesh-dependent norm
kvhk2h := X E∈Th kν1/2∇ xvhk2L2(E)+ θEhEk∂tvhk 2 L2(E) + 1 2kvhk 2 L2(ΣT). (4.12)
For the following derivations, we assume that our triangulation Th of Q is shape regular
such that the local approximation error estimates are available, [2, 4]. The triangula-tion Th of Q is called quasi-uniform, if there exists a constant cu such that
hE ≤ h ≤ cuhE, for all E ∈ Th, (4.13)
where h = maxE∈ThhE. Moreover, we introduce localised bounds for our coefficient function ν, i.e.,
νE ≤ ν(x, t) ≤ νE, for almost all (x, t) ∈ E and for all E ∈ Th, (4.14)
where νE and νE = const. > 0. In the following, we need some inverse inequalities for
functions from finite element spaces.
Lemma 4.3. There exist generic positive constants cI,1 and cI,2, such that
kvhkL2(∂E) ≤ cI,1h −1/2 E kvhkL2(E), (4.15) k∇vhkL2(E) ≤ cI,2h −1 E kvhkL2(E) (4.16)
for all vh ∈ V0h and for all E ∈ Th.
Proof. For (4.15), see e.g. [22, 5], and for (4.16) see e.g. [2, 4, 5]. From ∇ = (∇x, ∂t)> and (4.16), we can immediately deduce
k∂tvhkL2(E) ≤ cI,2h
−1
E kvhkL2(E). (4.17)
The above inequalities hold for the standard norms. However, we will also need such a result in some scaled norm.
Lemma 4.4. Let ν ∈ W∞1(Th) be a given uniformly positive function. Then
kvk2 Lν 2(E) = Z E ν(x, t) |v(x, t)|2 d(x, t)
is a norm and there holds the inverse estimate
k∂tvhkLν2(E) ≤ k∇vhkLν2(E)≤ cI,νh
−1
E kvhkLν2(E), (4.18) for all vh ∈ V0h and for all E ∈ Th.
Proof. If ν = νE = const > 0 on E, then (4.18) is nothing else than the classical
inverse inequality (4.16). In general, we can at least assume that (4.14) holds. Using (4.14) and (4.16), we obtain k∇vhkLν 2(E) ≤ √ νEk∇vhkL2(E)≤ √ νEcI,2h−1E kvhkL2(E) ≤ νE νE 1/2 cI,2 | {z } =:cI,ν h−1E kvhkLν 2(E)
It is clear that in 1 ≤ νE/νE is close to 1 in practical applications.
Below, we will need the estimate
k∂t∂xivhkLν2(E) ≤ cI,νh
−1
E k∂xivhkLν2(E), (4.19) which obviously holds for all vh ∈ V0h and for all E ∈ Th. Moreover, we need the
following inverse inequality. Lemma 4.5. Let ν ∈ W1
∞(Th) be a given uniformly positive function. Let Wh|E :=
{wh : wh = ∇xvh, vh ∈ V0h|E}. Then there holds the inverse estimate
k divx(νwh)kL2(E) ≤ cI,3h
−1
E kνwhkL2(E), ∀wh ∈ Wh|E, (4.20) where cI,3 is a positive constant, independent of hE.
Proof. First, we know that V0h|E is a finite space spanned by the local shape functions
{p(i)}
i∈ωE. Hence the space Wh|E is also finite and spanned by the generating system {∇xp(i)}i∈ωE. Moreover, for a fixed ν, each product zh := ν wh can be represented by means of a non-necessary unique linear combination {ν ∇xp(i)}i∈ωE on E. We denote this space by Zh(E) := spani∈ωE{ν ∇xp
(i)}. Using Cauchy’s inequality, we obtain
k divxzhk2L2(E) = Z E | divxzh|2 d(x, t) = Z E | d X i=1 ∂xizh,i| 2 d(x, t) ≤ d Z E d X i=1 |∂xizh,i| 2 d(x, t) = d d X i=1 k∂xizh,ik 2 L2(E),
for all zh ∈ Zh(E). Now, by a simple scaling argument, we can estimate each element
in the sum and obtain
d d X i=1 k∂xizh,ik 2 L2(E) ≤ d d X i=1 C2h−2E kzh,ik2L2(E) = dC2h−2E kzhk2L2(E).
Indeed, transforming to the reference triangle, using the norm equivalence on finite dimensional spaces, and transforming back to E, we obtain
k∂xizh,ik 2 L2(E)≤ k∇zh,ik 2 L2(E)= Z E |∇zh,i|2 d(x, t) ≤ chd+1E Z ∆ |∇ˆzh,i|2d(ξ, τ ) ≤ chd+1E h −2 E Z ∆ | ˆ∇ˆzh,i|2d(ξ, τ ) ≤ Ch−2E Z E
|zh,i|2 d(x, t) = Ch−2E kzh,ikL2(E).
Taking the square root and setting cI,3 := C
√
d closes the proof.
Lemma 4.5 gives information how the two norms involved scale wrt the mesh-size hE.
However, the estimate (4.20) is not sharp wrt the constant. Lemma 4.6. Let the assumptions of Lemma 4.5 hold. Then
k divx(νwh)kL2(E)≤ coptkνwhkL2(E), ∀wh ∈ Wh|E (4.21)
with c2
opt = sup06=zh∈Zh(E)
k divx(zh)k2L2(E)
kzhk2L2(E) .
Proof. From Lemma 4.5 we know that there must be a constant c such that k divx(zh)kL2(E)≤ ckzhkL2(E) ∀zh ∈ Zh(E).
With the assumption zh 6= 0 we can rewrite the inequality above as
k divx(zh)k2L2(E)
kzhk2L2(E)
≤ c2.
Now we immediately see that the optimal value for c is nothing else than the supremum of the expression on left hand side, i.e.,
c2opt := sup
06=zh∈Zh(E)
k divx(zh)k2L2(E)
kzhk2L2(E)
.
What remains is to ensure that this supremum is finite. We start by identifying the kernel of kν∇x. kL2(E). Using the notation of the proof of Lemma 4.5, we know
0 =kzhkL2(E)= k X i∈ωE ziq(i)kL2(E) =kX i∈ωE ziν∇xp(i)kL2(E) = kν∇x X i∈ωE zip(i) | {z } =ϕ kL2(E).
This identity holds if and only if ∇xϕ ≡ 0, i.e., if ϕ = ϕ(t). Now let ϕ ∈ ker kν∇x. kL2(E). Then we immediately deduce that
k divx(zh)kL2(E)= k divx(ν∇x X
i∈ωE
zip(i))kL2(E) = k divx(ν ∇xϕ | {z }
=0
kL2(E) = 0,
i.e., ker kν∇x· kL2(E) ⊂ ker k divx(ν∇x. )kL2(E).
Remark 4.7. Note that the constant copt in Lemma 4.6 is not only optimal but also
computable Let zh ∈ Zh(E), then by definition we have
zh(x, t) =
X
j∈˜ωE ˜ zjq˜(j).
Here we assume that the {˜q(j)}
j∈˜ωE form a basis of Zh(E). Moreover, we know kzhk2L2(E) = (zh, zh)L2(E) and k divxzhk
2
L2(E)= (divxzh, divxzh)L2(E)
| {z }
=:b(zh,zh)
.
As our space Zh(E) is finite, we can further rewrite these inner products. Let yh, zh ∈
Zh(E), then (yh, zh)L2(E) = X i∈˜ωE (yh, ˜q(i))L2(E)z˜i = X i,j∈˜ωE ˜ yj(˜q(j), ˜q(i))L2(E)z˜i.
This can be interpreted as
(yh, zh)L2(E)= (Mhy, z)`2, with (Mh)ij = (˜q
(j), ˜q(i)) L2(E),
where y and z are the vector of coefficients wrt the basis. By the same argument we obtain
b(yh, zh) = (Bhy, z)`2, with (Bh)ij = b(˜q
(j), ˜q(i)) L2(E). Combining the above identities, we get with NE = |ωE|
c2opt = sup 06=zh∈Zh(E) k divx(zh)k2L2(E) kzhk2L2(QT)(E) = sup z∈RNE (Bhz, z)`2 (Mhz, z)`2 . (4.22)
Hence, c2opt is the largest eigenvalue of the generalized eigenvalue problem Bhz = λMhz.
Now, we are able to proof the following lemma. Lemma 4.8. There exits a constant µa such that
ah(vh, vh) ≥ µakvhk2h, ∀vh ∈ V0h, (4.23)
with µa = minE∈Th1 − cI,3 q νEθE 4hE ≥ 1 2 for θE ≤ hE c2 I,3νE , i.e., µa= 12 for θE = c2hE I,3νE .
Proof. We first do integration by parts at the last term, obtaining ah(vh, vh) = X E∈Th Z E 1 2∂t(v 2 h) + θEhE(∂tvh)2+ ν|∇xvh|2 d(x, t) + Z E θEhEν∇xvh∇x∂tvh d(x, t) − Z ∂E θEhEν∇xvhnx∂tvh ds(x,t) = X E∈Th Z E 1 2∂t(v 2 h) d(x, t) + θEhEk∂tvhk2L2(E)+ Z E ν|∇xvh|2 d(x, t) − Z E θEhE divx(ν∇xvh)∂tvh d(x, t)
Now using Gauss’ theorem and the facts that vh is continuous across the element
boundary and that nt = 0 on Σ, we obtain
ah(vh, vh) = X E∈Th Z ∂E 1 2v 2 hntds(x,t)+ θEhEk∂tvhk2L2(E) + Z E ν|∇xvh|2− θEhE divx(ν∇xvh)∂tvh d(x, t) = 1 2 kvhk 2 L2(ΣT)− kvhk 2 L2(Σ0) + X E∈Th θEhEk∂tvhk2L2(E) + Z E ν|∇xvh|2− θEhE divx(ν∇xvh)∂tvh d(x, t)
The first, second and third term already appear in the definition of our discrete norm (4.12). What remains is to estimate the last term. Using the Cauchy-Schwarz inequal-ity, Lemma 4.5 and a scaled Young’s inequalinequal-ity, we arrive at the estimates
|θEhE
Z
E
divx(ν∇xvh)∂tvh d(x, t)| ≤ θEhEk divx(ν∇xvh)kL2(E)k∂tvhkL2(E) ≤ θEhEcI,3h−1E kν∇xvhkL2(E)h −1/2 E h 1/2 E k∂tvhkL2(E) ≤ cI,3 ενEθE 2hE k∇xvhk2Lν 2(E)+ 1 2εθEhEk∂tvhk 2 L2(E). Using this estimate in the equality above and the fact that vh = 0 on Σ0, we get
ah(vh, vh) ≥ 1 2kvhk 2 L2(ΣT)+ X E∈Th 1 −cI,3 2εθEhEk∂tvhk 2 L2(E) + 1 − εcI,3νEθE 2hE k∇xvhk2Lν 2(E).
Now we choose ε =phE/(θEνE) and obtain ah(vh, vh) ≥ min E∈Th 1 − cI,3 r θEνE 4hE × X E∈Th k∇xvhk2Lν 2(E)+ θEhEk∂tvhk 2 L2(E) + 1 2kvhk 2 L2(ΣT) ≥ µakvhk2h,
which concludes the first part of the proof. The second assertion can be shown by a simple calculation, i.e.,
1 − cI,3 r θEνE 4hE ≥ 1 2 ⇔ cI,3 r θEνE 4hE ≤ 1 2 ⇔ c2 I,3 θEνE hE ≤ 1 ⇔ θE ≤ hE νEc2I,3 .
Remark 4.9. The above proof does hold for any polynomial degree p ≥ 1 of vh and
any fixed, uniformly positive ν ∈ L∞(Q). However, for the special case p = 1 and
ν|E = const, the above proof is trivial, since
∂t(∇xvh) ≡ 0 and ν|E∆xvh ≡ 0.
Hence, there holds the identity
ah(vh, vh) = X E∈Th Z E ∂tvhvh+ θEhE(∂tvh)2+ ν|∇xvh|2 d(x, t) = X E∈Th 1 2 Z ∂E vh2ntds(x,t)+ θEhEk∂tvhk2L2(E)k∇xvhk 2 Lν 2(E) =kvhk2h,
i.e., µa= 1. Moreover, we immediately deduce that for this special case, the choice of
θE has no influence on the ellipticity of the space-time finite element method.
Remark 4.10. An alternative approach to the proof of Lemma 4.8 consists of not applying integration by parts on the last two terms of (4.9), but instead estimate
θEhE Z E ν∇xvh∇x(∂tvh) d(x, t) and θEhE Z ∂E ν∇xvhnx∂tvh ds(x,t) separately.
Lemma 4.8 already ensures uniqueness of the finite element solution uh ∈ V0h.
Fur-thermore, since we use the same trial- and test-space V0h, and this space is finite
dimensional, uniqueness implies existence of finite element solution uh ∈ V0h of (4.8).
For the special case of uniform meshes and uniform θ, i.e., hE = h and θE = θ for
all E ∈ Th, and and ν ≡ 1, a proof for ellipticity with a mesh-independent constant
was done by Langer et.al. [16]. For a second special case, where θE vanishes, i.e.,
θE = θ = 0 for all E ∈ Th, Steinbach in [28] has shown existence and uniqueness of
discrete version of (4.8). In addition, both papers include also a priori error estimates, where Steinbach’s estimate is based on a discrete inf-sup condition.
To show an a priori error estimate wrt the mesh dependent norm (4.12), we need to show that our bilinear form ah( . , . ) is uniformly bounded on V0h,∗× V0h, where
V0h,∗= H01,0(Q) ∩ H2(Th) + V0h with the norm
kvk2 h,∗ =kvk 2 h+ X E∈Th (θEhE)−1kvk2L2(E)+ θEhE|v| 2 H2(E) =1 2kvk 2 L2(ΣT)+ X E∈Th θEhEk∂tvk2L2(E)+ k∇xvkLν2(E) + (θEhE)−1kvk2L2(E)+ θEhE|v| 2 H2(E) (4.24)
Moreover, we will make use of the following scaled trace inequality. Lemma 4.11. There exists a positive constants cT r > 0 such that
kvk2 L2(∂E) ≤ 2c 2 T rh −1 E kvk 2 L2(E)+ h 2 Ek∇vk 2 L2(E) (4.25) for all v ∈ H1(E), ∀E ∈ Th.
Proof. See e.g. [22].
Lemma 4.12. The discrete bilinear form ah( . , . ) is uniformly bounded on V0h,∗× V0h,
i.e.,
|ah(u, vh)| ≤ µbkukh,∗kvhkh, (4.26)
where µb = maxE∈Th2(1 + θEh
−1 E c2T r ν2 E νE), 2c 2 T rν 2 E, 2 + c2I,1, 1 + (cI,νθE)2 1/2 that is bounded provided that θE = O(hE).
Proof. We will estimate the bilinear form (4.9) term by term. For the first term, since V0h⊂ H
1,1
0,0(Q), we can apply integration by parts and the Cauchy-Schwarz inequality,
and obtain X E∈Th Z E ∂tuvh d(x, t) = X E∈Th − Z E u∂tvh d(x, t) + Z ∂E untvh ds(x,t) ≤ X E∈Th (θEhE)−1kuk2L2(E) 1/2 (θEhEk∂tvhk2L2(E) 1/2 + kuk2L2(ΣT)1/2 kvhk2L2(ΣT) 1/2 .
For the second and third term, applying the Cauchy-Schwarz inequality for each term of the sum yields
θEhE Z E ∂tu∂tvh d(x, t) ≤ θEhEk∂tuk2L2(E) 1/2 θEhEk∂tvhk2L2(E) 1/2 , Z E ν∇xu∇xvh d(x, t) ≤ k∇xuk2Lν 2(E) 1/2 k∇xvhk2Lν 2(E) 1/2 ,
respectively. For the fourth term, we use again Cauchy-Schwarz’ inequality, the inverse estimate (4.19), and obtain
θEhE Z E ν∇xu∇x(∂tvh) d(x, t) ≤ k∇xuk2Lν 2(E) 1/2 (θEhE)2k∂t∇xvhk2Lν 2(E) 1/2 = k∇xuk2Lν 2(E) 1/2 (θEhE)2 d X i=1 k∂t(∂xivh)k 2 Lν 2(E) 1/2 ≤ k∇xuk2Lν 2(E) 1/2 (θEhE)2 d X i=1 c2I,νh−2E k∂xivhk 2 Lν 2(E) 1/2 = k∇xuk2Lν 2(E) 1/2 (cI,νθE)2k∇xvhk2Lν 2(E) 1/2 .
For the last term, we apply Cauchy-Schwarz and the trace inequalities (4.17) and (4.25), and get θEhE Z ∂E ν∇xunx∂tvh ds(x,t) ≤ θEν2Ek∇xuk2L2(∂E) 1/2 θEh2Ek∂tvhk2L2(∂E) 1/2 ≤ 2θEν2Ec 2 T rh −1 E k∇xuk2L2(E)+ h 2 E d X i=1 k∇∂xiuk 2 L2(E) 1/2 × θEhEc2I,1k∂tvhk2L2(E) 1/2 ≤ 2θEc2T r ν2E νEh −1 E k∇xuk2Lν 2(E)+ 2c 2 T rν2EθEhE|u|2H2(E) 1/2 × c2I,1θEhEk∂tvhk2L2(E) 1/2 .
items, i.e., |ah(u, vh)| ≤ kuk2L2(ΣT) 1/2 kvhk2L2(ΣT) 1/2 + X E∈Th (θEhE)−1kuk2L2(E) 1/2 θEhEk∂tvhk2L2(E) 1/2 + θEhEk∂tuk2L2(E) 1/2 θEhEk∂tvhk2L2(E) 1/2 + k∇xuk2Lν 2(E) 1/2 k∇xvhk2Lν 2(E) 1/2 + k∇xuk2Lν 2(E) 1/2 (cI,νθE)2k∇xvhk2Lν 2(E) 1/2 + 2θEc2T r νE νEh −1 E k∇xuk 2 Lν 2(E)+ 2c 2 T rθEνEhE|u|2H2(E) 1/2 × c2 I,1θEhEk∂tvhk2L2(E) 1/2 ≤ kuk2L2(ΣT)+ X E∈Th θEhEk∂tuk2L2(E)+ 2(1 + θEc 2 T r ν2 E νEh −1 E )k∇xuk 2 Lν 2(E) + (θEhE)−1kuk2L2(E)+ 2c 2 T rν 2 EθEhE|u|2H2(E) 1/2 × kvhk2L2(ΣT)+ X E∈Th (2 + c2 I,1)θEhEk∂tvhk2L2(E) + 1 + (cI,1θE)2k∇xvhk2Lν 2(E) 1/2 ≤ max E∈Th 2(1 + θEh−1E c 2 T r ν2 E νE), 2c 2 T rν 2 E, 2 + c 2 I,1, 1 + (cI,νθE)2 1/2 | {z } =:µb kukh,∗kvhkh.
Choosing now θE = O(hE) ensures the boundedness of the constant µb.
Remark 4.13. Choosing θE as in Lemma 4.8, i.e., θE = hE/(c2I,3νE), we obtain
µa = 1/2 and µb = maxE∈Th2(1 +
νEc2T r νEc2 I,3 ), 2c2T rν2E, 2 + c2I,1, 1 + (cI,νhE c2 I,3νE )2 1/2.
Remark 4.14. As in Remark 4.9, we can provide a simplified estimate for the special case p = 1 and ν|E = νE = const. The first three terms can be estimated as in the
above proof. The fourth term completely vanishes, since ∇x(∂tvh) = 0. For the fifth
inequality, obtaining θEhE Z ∂E νE∇xu · nx∂tvh ds(x,t) = θEhEνE∂tvh Z ∂E ∇xu · nx ds(x,t) = θEhEνE∂tvh Z E divx(∇xu) d(x, t) = θEhEνE Z E ∆xu∂tvh d(x, t) ≤ θEhEνE2k∆xuk2L2(E) 1/2 θEhEk∂tvhk2L2(E) 1/2 . Gathering the terms from the proof and the above estimate, we obtain
|ah(u, vh)| ≤ kuk2L2(ΣT) 1/2 kvhk2L2(ΣT) 1/2 + X E∈Th (θEhE)−1kuk2L2(E) 1/2 θEhEk∂tvhk2L2(E) 1/2 + θEhEk∂tuk2L2(E) 1/2 θEhEk∂tvhk2L2(E) 1/2 + k∇xuk2Lν 2(E) 1/2 k∇xvhk2Lν 2(E) 1/2 + θEνE2k∆xuk2L2(E) 1/2 θEhEk∂tvhk2L2(E) 1/2 ≤ kuk2 L2(ΣT)+ X E∈Th θEhEk∂tuk2L2(E)+ k∇xuk 2 Lν 2(E) + (θEhE)−1kuk2L2(E)+ ν 2 EθEhEk∆xukL2(E) 1/2 × kvhk2L2(ΣT) + X E∈Th 3θEhEk∂tvhk2L2(E)+ k∇xvhk 2 Lν 2(E) 1/2 ≤ max E∈Th {3, ν2 E} 1/2kuk h,∗∗kvhkh.
We immediately deduce that this new constant ˜µb = maxE∈Th{3, ν
2
E}1/2 is also
inde-pendent of hE.
To obtain a priori error estimate wrt to the mesh dependent norm (4.12), we need interpolation error estimates for finite elements wrt (4.24), which we summarise in the next Lemmata. Moreover, we need the broken Sobolev space
Hs(Th) := {v ∈ L2(Q) : v|E ∈ Hs(E)}, (4.27)
equipped with the broken Sobolev (semi-)norm |v|2 Hs(T h) := X E∈Th |v|2 Hs(E) and kvk2Hs(T h):= X E∈Th kvk2 Hs(E), (4.28) where s is some positive integer. For further details on such spaces, we refer to [5, 22].
Lemma 4.15. Let s and k be positive integers with s ∈ [2, p + 1] and k > (d + 1)/2, respectively. Let v ∈ V0∩ Hk(Q) ∩ Hs(Th). Then there exists an interpolation operator
Πh, mapping from V0∩ Hk(Q) to V0h, such that
kv − ΠhvkL2(E)≤ Ch s+1 E |v|Hs(E), (4.29) k∇(v − Πhv)kL2(E)≤ Ch s E|v|Hs(E), (4.30) |v − Πhv|H2(E)≤ Chs−1E |v|Hs(E), (4.31) where C is some generic constant independent of v. Here, p denotes the polynomial degree of the finite element basis functions.
Proof. See e.g. [3, Theorem 4.4.4] or [4, Theorem 3.1.6].
Lemma 4.16. Let the assumptions of Lemma 4.15 hold. Then the following interpo-lation error estimates hold:
kv − ΠhvkL2(ΣT) ≤ c1 X E∈Th ∂E∩ΣT6=∅ h2s−1E |v|2 Hs(E) 1/2 , (4.32) kv − Πhvkh ≤ c2 X E∈Th h2(s−1)E |v|Hs(E) 1/2 , (4.33) kv − Πhvkh,∗ ≤ c3 X E∈Th h2(s−1)E |v|Hs(E) 1/2 . (4.34)
The constants c1, c2, c3 do not depend on hE or v, provided that θE = O(hE).
Proof. We start with the first estimate (4.32). We use the scaled trace inequality (4.25), and the interpolation error estimates (4.29) and (4.30), obtaining
kv − Πhvk2L2(ΣT) = X E∈Th ∂E∩ΣT6=∅ kv − Πhvk2L2(∂E∩ΣT) ≤ X E∈Th ∂E∩ΣT6=∅ kv − Πhvk2L2(∂E) ≤X E∈Th ∂E∩ΣT6=∅ 2c2 T rh −1 E (kv − Πhvk2L2(E)+ h 2 Ek∇(v − Πhv)k2L2(E)) ≤c2 T r X E∈Th ∂E∩ΣT6=∅ C2h2s−1 E |v| 2 Hs(E)+ C2h2s−1E |v|2Hs(E) ≤c2T rC2X E∈Th ∂E∩ΣT6=∅ [h2s−1E |v|Hs(E)].
(4.30), and the above estimate (4.32), and obtain kv − Πhvk2h = X E∈Th θEhEk∂t(v − Πhv)k2L2(E)+ k∇x(v − Πhv)k 2 Lν 2(E) + 1 2kv − Πhvk 2 L2(ΣT) ≤ X E∈Th θEC2h 2(s−1) E |v| 2 Hs(E)+ νEC2h 2(s−1) E |v| 2 Hs(E) + 1 2c 2 1 X E∈Th h2s−1E |v|2 Hs(E) ≤ X E∈Th (C2θEhE + νEC2+ c21hE)h 2(s−1) E |v| 2 Hs(E) .
For the last estimate (4.34), we use definition (4.24), the above estimate (4.33), and the interpolation error estimate (4.31), obtaining
kv − Πhvk2h,∗ =kv − Πhvk2h+ X E∈Th (θEhE)−1kv − ΠhvkL22(E)+ θEhE|v − Πhv| 2 H2(E) ≤ X E∈Th c2 2h 2(s−1) E |v| 2 Hs(E)+ C2θ −1 E h 2s−1 E |v| 2 Hs(E)+ C2θEhEh 2(s−2) E |v| 2 Hs(E) ≤ X E∈Th c22+ hEθE−1C 2+ θ Eh−1E C 2h2(s−1) E |v| 2 Hs(E).
The special choice θE = O(hE) ensures that the constant c3 is independent of hE.
Remark 4.17. The strong assumption v ∈ Hk(Q) with k > (d + 1)/2 is needed for
the interpolation error estimates for the Lagrange interpolation operator. However, in practical application this requirement is most likely never met. However, in such a practical application, the space-time cylinder Q =SM
i=1Qi can be split into subdomains
Qi, which correspond e.g. to different materials. On each such subdomain Qi, we can
assume some regularity for v ∈ Hs(T (Q)) := {v ∈ L
2(Q) : v|Qi ∈ H
s(Q
i), for all i =
1, . . . , M } and the diffusion coefficient ν ∈ W1
∞(T (Q)) := {v ∈ L∞(Q) : v|Qi ∈ W∞1(Qi), for all i = 1, . . . , M }. For a similar case, Duan et.al. [7] have shown an
interpolation error estimate of the form
k∇(v − Ihv)kL2(Q) ≤ Ch s−1 M X i=1 kvkHs(Q i). Now we can formulate the following a priori estimate for the error.
Theorem 4.18. Let s and k be positive integers with s ∈ [2, p + 1] and k > (d + 1)/2. Furthermore, let u ∈ V0 ∩ Hk(Q) ∩ Hs(Th) be the exact solution, and uh ∈ V0h
the solution of the finite element scheme (4.11). Then there holds the a priori error estimate ku − uhkh ≤ c X E∈Th h2(s−1)E |u|2 Hs(E) 1/2 . (4.35)
Proof. First, we know from the consistency identity (4.8) that ah(u, vh) = lh(vh), and,
since uh is the approximate solution of (4.11), that ah(uh, vh) = lh(vh). Hence we have
Galerkin orthogonality for our bilinear form ah( . , . ), i.e.
ah(u − uh, vh) = 0, ∀vh ∈ V0h. (4.36)
We start with the triangle inequality for the discretization error, i.e., ku − uhkh ≤ ku − Πhukh+ kΠhu − uhkh.
We continue by estimating the second term. Using the ellipticity proved in Lemma 4.8, the Galerkin orthogonality and the generalized boundedness from Lemma 4.12, we obtain
µakΠhu − uhk2h ≤ ah(Πhu − uh, Πhu − uh) = ah(Πhu − u, Πhu − uh)
≤ µbkΠhu − ukh,∗kΠhu − uhkh.
We insert this estimate in the triangle inequality above, use the interpolation error estimates (4.33) and (4.34), and obtain
ku − uhkh ≤ ku − Πhukh+ µb µa kΠhu − ukh,∗ ≤ c2 X E∈Th h2(s−1)E |u|2Hs(E) 1/2 + c3 µb µa X E∈Th h2(s−1)E |u|2Hs(E) 1/2 ≤ c2+ c3 µb µa X E∈Th h2(s−1)E |u|2 Hs(E) 1/2 ,
which proves the estimate (4.35) with c = c2+ c3(µb/µa).
The Finite Element Method
Now we proceed with solving the discrete variational problem (4.11) that is nothing but a huge system of linear algebraic equations. Indeed, let {p(i) : i ∈ Ih} be some basis
of V0h, where Ih is some index set, which we will specify later. Then we can express
the approximate solution uh in terms of this basis, i.e. uh(x, t) =
P
i∈Ihuip
(i)(x, t).
Furthermore, each basis function is a valid test function. Thus, we obtain Nh equations
from (4.11),
ah(uh, p(i)) = lh(p(i)), for all i ∈ Ih, (4.37)
where Nh = |Ih| is the dimension of V0h. Now we replace uh by its basis representation,
which yields X j∈Ih ui ah(p(j), p(i)) | {z } =:(Kij) = lh(p(i)) | {z } =:(fi) , for all i ∈ Ih. (4.38)
We can rewrite this system in terms of a system of linear algebraic equations
where Kh = (Kij), uh = (ui) and fh = (fi). The system matrix is non-symmetric,
but positive definite due to Lemma 4.8. Indeed,
(Khvh, vh) = ah(vh, vh) ≥ µakvhk2h > 0 (4.40)
for all V0h 3 vh ↔ vh ∈ RNh : vh 6= 0. The linear system (4.39) can be solved
efficiently and most important in parallel by either a sparse direct solver (e.g sparse LU-factorisation) or an iterative solver (e.g. preconditioned GMRES).
But how to construct such a basis {p(i)} of V0h? We need again the regular
triangu-lation Th of our space-time domain Q, which we already introduced in (4.1). We now
define shape functions and the corresponding function set F (E). This function set is either a subspace or equal to the following function spaces
Pk := { X |α|≤k cαxα : cα ∈ R}, Qk := { X αi≤k cαxα: i = 1, . . . , d + 1, cα ∈ R}.
Furthermore, we need some degrees of freedom, denoted by l(E,α), which are functionals in the dual space F (E)∗ of F (E). If these functionals span the whole dual space, i.e. they are a basis, or, equivalently, have the interpolation property
l(E,α)(v) = cα, for v ∈ F (E),
then we can uniquely determine all coefficients of v ∈ F (E) on an element E. A particular example for a degree of freedom is the point evaluation, i.e., let v ∈ F (E), then
l(E,α)(v) = v(x(E,α)), α ∈ AE, (4.41)
where x(E,α) ∈ E is called node and A
E is a set of local indices. Combining all of the
above, we can now introduce a formal definition for a finite element.
Definition 4.19 (Finite Element [2]). Let Th be an admissible triangulation, i.e., for
E, E0 ∈ Th, E ∩ E0 = ∅, common vertex, common edge, common face (3D). (4.42)
A triple (E, F (E), {l(E,α)}) is called a finite element with the following properties:
(i) E is a polyhedron in Rd+1.
(ii) F (E) is a subspace of C(E) with finite dimension p.
(iii) {l(E,α)} is a set of p linear independent functionals over F (E) which uniquely
Hence, we have local nodal basis of shape functions, i.e.,
{p(E,α): p(E,α)∈ F (E), α ∈ AE}, (4.43)
with the property l(E,α)(p(E,β)) = δ αβ.
Now we define the global set of nodes {x(i) : i ∈ I
h}, and if x(i) ∈ E, then x(i) = x(E,α)
for some α. Furthermore, we need a global set of degrees of freedom {l(i) : i ∈ Ih},
where l(i) = l(E,α) on E, and a global nodal basis {p(i) : i ∈ I
h}, with l(i)(p(j)) = δij
and p(i) = p(E,α) on E. Then this global nodal basis spans our discrete function space
V0h= span{p(i) : i ∈ Ih}.
For the rest of this thesis, we will restrict ourselves to the polynomial space Pp on
d + 1-simplices, e.g., triangles in 2D and tetrahedrons in 3D. The degree of freedom is the point evaluation. This class of finite elements belongs to the class of Lagrange finite elements.
(a) 2D
(b) 3D
Figure 4.1: The linear and quadratic Lagrange finite element for different dimensions. The black dot denotes the point evaluation in the vertex.
To efficiently compute the entries in Kh and fh, we observe that the nodal basis
we can write each entry as (Kh)ij = ( 0, if Bij = Bi∩ Bj = ∅, P E∈Bijah,E(p (j), p(i)), else , (fh)i = X E∈Bi lh,E(p(i)),
where Bi = {E ∈ Th : x(i) ∈ E} is the neighbourhood of a node x(i) and
ah,E(uh, vh) := Z E ∂tuhvh+ θEhE∂tuh∂tvh d(x, t) + Z E ν∇xuh· ∇xvh+ θEhEν∇xuh· ∇x(∂tvh) d(x, t) (4.44) − Z ∂E θEhEν∇xuh· nxvh ds(x,t), lh,E(vh) := Z E f (vh+ θEhE∂tvh) d(x, t), (4.45)
are the integrals over one element.
In order to compute the entries of Kh, we will assemble the stiffness matrix Kh and
the load vector fh element-wise, i.e., on each element E, we have to identify i ↔ α and j ↔ β, so we obtain the local element matrix K(E)h and the element load vector f(E)h , with
(K(E)h )αβ = ah,e(p(E,β), p(E,α)), and (f (E)
h )α = lh,e(p(E,α)).
In order to avoid the computation of the coefficients of p(E,α) on each element, we
will instead transform the arbitrary element E to a canonical element, the so called reference element ∆, with
∆ := {(ξ1, . . . , ξd+1) : d+1
X
i=1
ξi ≤ 1 ∧ ξi ≥ 0, for i = 1, . . . , d + 1}.
For our finite elements, it is sufficient to do this transformation via an affine mapping XE, which is defined as
ξ 7→ XE(ξ) := x(E,1)+
x(E,2)− x(E,1) ... . . . ... x(E,d+1)− x(E,1)
| {z }
JE
ξ,
where x(E,α)are the vertices of the element E ∈ T
h. For instance, for d = 1, the shape
functions p(α) can be easily computed on the reference element ∆, i.e., p(1)(ξ) = 1 − ξ1− ξ2, p(2) = ξ1, and p(3)= ξ2.
