Numerical Methods

for General Types of PDEs

Lecture Notes

Mario Bebendorf
Contents

1 Non-Conforming and Mixed Methods 1

1.1 Non-Conforming Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 The Crouzeix-Raviart element . . . . . . . . . . . . . . . . . . . . . . 4
1.1.2 Polygonal Approximation of Curvilinear Boundaries . . . . . . . . . . 8
1.2 Mixed Finite Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.2.1 The inf-sup Condition . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3 Saddle Point Problems and Constrained Variational Problems . . . . . . . . 15
1.3.1 Conforming Approximation of Mixed Variational Problems . . . . . . 20
1.4 Solvability of the Mixed Formulation of the Poisson Problem . . . . . . . . . 22
1.4.1 The Raviart-Thomas Element . . . . . . . . . . . . . . . . . . . . . . 24
1.5 Computing the Solution of the Discrete Problem . . . . . . . . . . . . . . . . 27
1.5.1 The Uzawa Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
1.5.2 The Bramble-Pasciak CG Method . . . . . . . . . . . . . . . . . . . . 29

2 Stokes’ Equations 34
2.1 Variational Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.2 The Discrete LBB Condition for Stokes’ Equations . . . . . . . . . . . . . . 36
2.3 Solution of the Discrete Problems . . . . . . . . . . . . . . . . . . . . . . . . 43

3 Finite Element Methods for Electro-Magnetic Problems 45

3.1 The Space H(curl; Ω) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.2 Variational Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.3 Approximation of the Variational Problem . . . . . . . . . . . . . . . . . . . 53
3.4 Conforming Finite Elements for H(curl; Ω) . . . . . . . . . . . . . . . . . . . 55
3.5 Assembling Finite Element Stiffness Matrices . . . . . . . . . . . . . . . . . . 58

4 Generalized Finite Elements 60

4.1 Partition of Unity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.2 GFEM and Classical FEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.3 Choice of Local Approximation Spaces . . . . . . . . . . . . . . . . . . . . . 71
4.4 Implementational Issues in the GFEM . . . . . . . . . . . . . . . . . . . . . 73

i
Preface

These notes summarize a lecture series given as a master’s course at the University of
Bayreuth in summer 2025. Any comments and corrections are appreciated.

Bayreuth, July 24, 2025

Literature:

• D. Braess: Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics.
Springer

• O. Steinbach: Numerical Approximation Methods for Elliptic Boundary Value Prob-

lems: Finite and Boundary Elements, Springer

• C. Großmann and H.-G. Roos: Numerical Treatment of Partial Differential Equations.

Teubner

• V. Girault and P.-A. Raviart: Finite Element Methods for Navier-Stokes Equations.
Springer, 1986

• D. Boffi, F. Brezzi, and M. Fortin: Mixed Finite Element Methods and Applications.
Springer, 2013

• P. Monk: Finite Element Methods for Maxwell’s Equations. Clarendon Press, 2003

ii
1 Non-Conforming and Mixed
Methods

1.1 Non-Conforming Methods

Conforming finite elements approximate a continuous problem

u∈V : a(u, v) = ℓ(v) for all v ∈ V (1.1)

by choosing a finite-dimensional subspace Vh ⊂ V with a discrete problem

uh ∈ Vh : a(uh , vh ) = ℓ(vh ) for all vh ∈ Vh . (1.2)

This has turned out to be useful in particular for the so-called Galerkin orthogonality. How-
ever, in the following situations, for instance, conforming discretizations turn out to be
inappropriate:

• the construction of conforming function spaces Vh ⊂ V in the case of equations of

higher order can be complicated;

• if quadrature formulas are used then the bilinear form a and the functional ℓ can only
be evaluated approximately;

• curvilinear boundaries of the computational domain do not allow to represent boundary

conditions exactly.

If the bilinear form a or the linear functional ℓ is perturbed, for instance, by numerical
integration, then (1.1) has to be replaced with the finite-dimensional problem

uh ∈ Vh : ah (uh , vh ) = ℓh (vh ) for all vh ∈ Vh . (1.3)

Here, we assume that ℓh ∈ Vh′ and ah : Vh × Vh → R are continuous on the Hilbert space Vh
and that ah is equi-coercive, i.e. there is a constant cK > 0 independent of h such that

ah (vh , vh ) ≥ cK ∥vh ∥2Vh for all vh ∈ Vh .

Due to the Lax-Milgram theorem, (1.3) is uniquely solvable. In addition to the FE approx-
imation error resulting from (1.2), (1.3) introduces a consistency error. The impact of these
perturbations on the solution is described by the following two Strang lemmas, which can
be viewed as generalizations of Céa’s lemma. As a first step, we assume that Vh ⊂ V and
∥ · ∥Vh = ∥ · ∥V .

1
1 Non-Conforming and Mixed Methods

Lemma 1.1 (first Strang lemma). Let u ∈ V be a solution of (1.1). Let the bilinear form
ah be equi-coercive and let a : V × V → R be continuous with continuity constant cS > 0.
Then the solution uh ∈ Vh ⊂ V of (1.3) satisfies
  
cS −1 |a(vh , wh ) − ah (vh , wh )
∥u − uh ∥V ≤ inf 1+ ∥u − vh ∥V + cK sup
vh ∈Vh cK 0̸=wh ∈Vh ∥wh ∥V
|ℓ(wh ) − ℓh (wh )|
+ c−1K sup .
0̸=wh ∈Vh ∥wh ∥V

Proof. From the equi-coercivity of ah we have for vh ∈ Vh and wh := uh − vh

cK ∥uh − vh ∥2V ≤ ah (uh − vh , uh − vh ) = ah (uh , wh ) − ah (vh , wh )

= ℓh (wh ) + a(u, wh ) − ℓ(wh ) − ah (vh , wh )
= a(u − vh , wh ) + {a(vh , wh ) − ah (vh , wh )} + {ℓh (wh ) − ℓ(wh )}.

Dividing by ∥wh ∥V yields

 
|a(u − vh , wh )| |a(vh , wh ) − ah (vh , wh )| |ℓ(wh ) − ℓh (wh )|
cK ∥uh − vh ∥V ≤ sup + +
0̸=wh ∈Vh ∥wh ∥V ∥wh ∥V ∥wh ∥V
 
|a(vh , wh ) − ah (vh , wh )| |ℓ(wh ) − ℓh (wh )|
≤ cS ∥u − vh ∥V + sup + .
0̸=wh ∈Vh ∥wh ∥V ∥wh ∥V

The assertion follows from the triangle inequality ∥u − uh ∥V ≤ ∥u − vh ∥V + ∥uh − vh ∥V .

Example 1.2 (Mass Lumping). Let Vh denote the space of continuous, piecewise linear
functions S 1,0 (Th ). We choose the quadrature formula
n
X
Q(f ) = |Dk | f (pk ),
k=1

where pk denotes the corner points of the triangles τ ∈ Th and Dk is a neighborhood of pk

(e.g. a cell in a Voronoi decomposition) such that {Dk , k = 1, . . . , n} is a partition of Ω ⊂ Rd ,
the so-called dual mesh.

Replacing ℓ(v) = Ω f v dx with the quadrature formula ℓh (v) := Q(f v), we obtain
R

X
|ℓ(vh ) − ℓh (vh )| ≤ chd+1 ∥∇vh ∥L2 (τ )
τ ∈Th

provided that f is Lipschitz-continuous. This special substitution of functionals ℓ with ap-

proximations ℓh is referred to as mass lumping. It is particularly useful for the discretization
of boundary value problems

−∆u + c(x)u = f in Ω,
u = 0 on ∂Ω.

2
 1.1 Non-Conforming Methods

The corresponding bilinear form

Z Z
a(u, v) = ∇u · ∇v dx + c uv dx
Ω Ω

is replaced with
Z n
X
ah (uh , vh ) = ∇uh · ∇vh dx + c(pk )|Dk |uh (pk )vh (pk ).
Ω k=1

If ah is applied to the Lagrange functions φi , φj ∈ Vh then we obtain

Z
ah (φj , φi ) = ∇φj · ∇φi dx + c(pi )|Di |δij .
Ω

Hence, using this technique has the particular advantage that the mass matrix M ∈ Rn×n
having the entries mij = (φj , φi )L2 becomes diagonal.

Finite element methods which do not satisfy Vh ⊂ V are referred to as non-conforming.

In contrast to conforming methods it cannot be assumed that ∥ · ∥V is defined on Vh . Hence,
grid-dependent norms ∥·∥Vh are employed which are defined on Vh +V . For instance, instead
of ∥ · ∥H 1 (Ω) the grid-dependent norm
!1/2
X
2
∥vh ∥Vh = ∥vh ∥H 1 (τ ) , vh ∈ Vh ,
τ ∈Th

can be used if Vh restricted to each τ ∈ Th is in H 1 (τ ). Notice that ∥v∥H 1 (Ω) = ∥v∥Vh for all
v ∈ H 1 (Ω).

Additionally we assume that the bilinear form ah is well-defined for functions in Vh + V .

The difference u − uh of the solution u of (1.1) and the solution uh ∈ Vh of (1.3) is estimated
in the following lemma with respect to the norm ∥ · ∥Vh of Vh .

Lemma 1.3 (second Strang lemma). Let ah be equi-coercive and continuous w.r.t. ∥·∥Vh .
Then
 
cS |ah (u, wh ) − ℓh (wh )|
∥u − uh ∥Vh ≤ 1 + inf ∥u − vh ∥Vh + c−1
K sup .
cK vh ∈Vh 0̸=wh ∈Vh ∥wh ∥Vh

Proof. Let vh ∈ Vh . From the equi-coercivity we obtain

cK ∥uh − vh ∥2Vh ≤ ah (uh − vh , uh − vh )
= ah (u − vh , uh − vh ) + {ℓh (uh − vh ) − ah (u, uh − vh )}.
Using wh := uh − vh , we obtain after division by ∥wh ∥Vh
|ah (u, wh ) − ℓh (wh )|
cK ∥uh − vh ∥Vh ≤ cS ∥u − vh ∥Vh + .
∥wh ∥Vh
Similar to the first Strang lemma, the assertion follows from the triangle inequality.

3
1 Non-Conforming and Mixed Methods

For non-conforming discretezations the Aubin-Nitsche lemma (cf. Theorem 2.44 in Numerical
Methods for Partial Differential Equations) contains two additional terms.

Lemma 1.4 (Aubin-Nitsche). Let W be a Hilbert space in which V is continuously

embedded. Furthermore, let Vh ⊂ W and let the V + Vh -continuous bilinear form ah coincide
with a on V × V . Then
1 n
∥u − uh ∥W ≤ sup cS ∥u − uh ∥Vh ∥uφ − uφ,h ∥Vh + |ah (u − uh , uφ ) − (u − uh , φ)W |
0̸=φ∈W ∥φ∥W
o
+ |ah (u, uφ − uφ,h ) − ℓ(uφ ) + ℓh (uφ,h )| .

Here, uφ ∈ V and uφ,h ∈ Vh for each φ ∈ W denote the solutions of the dual problems
a(w, uφ ) = (w, φ)W , w ∈ V , and ah (w, uφ,h ) = (w, φ)W , w ∈ Vh , respectively.

Proof. Due to the embedding of V in W and Vh ⊂ W we have u − uh ∈ W . The norm of

u − uh can be expressed in terms of the dual norm

(u − uh , φ)W
∥u − uh ∥W = sup .
0̸=φ∈W ∥φ∥W

According to the definition of uh , uφ , and uφ,h , for φ ∈ W it holds

(u − uh , φ)W = ah (u, uφ ) − ah (uh , uφ,h )

= ah (u − uh , uφ − uφ,h ) + ah (uh , uφ − uφ,h ) + ah (u − uh , uφ,h ).

Using

ah (uh , uφ − uφ,h ) = ah (uh − u, uφ ) + ah (u, uφ ) − ah (uh , uφ,h )

= ah (uh − u, uφ ) + (u − uh , φ)W

and

ah (u − uh , uφ,h ) = ah (u, uφ,h − uφ ) + ah (u, uφ ) − ah (uh , uφ,h )

= ah (u, uφ,h − uφ ) + ℓ(uφ ) − ℓh (uφ,h ),

we obtain

(u − uh , φ)W = ah (u − uh , uφ − uφ,h ) − {ah (u − uh , uφ ) − (u − uh , φ)W }

− {ah (u, uφ − uφ,h ) − ℓ(uφ ) + ℓh (uφ,h )}.

The continuity of ah yields the assertion.

1.1.1 The Crouzeix-Raviart element

We are going to apply the previous results to the easiest non-conforming element, the follow-
ing Crouzeix-Raviart element, in the context of the solution of Poisson’s equation with

4
 1.1 Non-Conforming Methods

homogeneous boundary conditions. We assume a polygonal domain Ω ⊂ R2 . As the ansatz

space we use
Vh = {v ∈ L2 (Ω) : v|τ is linear for each τ ∈ Th ,
v is continuous at the edge centers}.
In order to satisfy the vanishing boundary conditions, we additionally require that v = 0 at
the edge centers on ∂Ω. A function v ∈ Vh is uniquely determined by its values at the edge

Figure 1.1: Crouzeix-Raviart element

centers m0 , m1 , m2 of τ ∈ Th . To see this, let z0 , z1 , z2 denote the vertices of τ ∈ Th and

φ0 , φ1 , φ2 the associated nodal basis functions. Then it is easily seen that
2
X
v|τ = v(mi ) ψi ,
i=0

where ψi := φi+1 + φi+2 − φi . Here, we have assumed that the i-th edge lies opposite the
i-th vertex and that all vertex indices are understood modulo 3. Notice that ψi (mj ) = δij .

We have Vh ⊂ L2 (Ω) and for v ∈ Vh

Z Z
1 1
v|τ ds = v(me ) = v|τ ds, e = τ 1 ∩ τ2 . (1.4)
|e| e 1 |e| e 2
However, the functions in Vh are not continuous and hence Vh ̸⊂ H 1 (Ω). Notice that the
piecewise linear conforming elements are contained in Vh because they are in particular
continuous at the edge centers.

Hence, we define
XZ Z
ah (u, v) = ∇u · ∇v dx, ℓh (v) = f v dx for u, v ∈ H01 (Ω) + Vh
τ ∈Th τ Ω

and !1/2
X
∥v∥Vh := |v|2H 1 (τ ) .
τ ∈Th

In view of the second Strang lemma we prove the following lemma.

Lemma 1.5. Let u ∈ H01 (Ω) ∩ H 2 (Ω) be the solution of (1.1) and Rh [u](w) := ah (u, w) −
ℓh (w). If Th is non-degenerate, then for all w ∈ Vh it holds

|Rh [u](w)| ≤ c h |u|H 2 (Ω) ∥w∥Vh .

5
1 Non-Conforming and Mixed Methods

Proof. Let w ∈ Vh . Then

XZ Z
Rh [u](w) = ah (u, w) − ℓh (w) = ∇u · ∇w dx − f w dx
τ ∈Th τ Ω
X Z Z  Z
= w ∂ν u ds − w∆u dx − f w dx
τ ∈Th ∂τ τ Ω
XZ Z XZ
= w ∂ν u ds − (∆u + f )w dx = w ∂ν u ds.
τ ∈Th ∂τ τ τ ∈Th ∂τ

Due to (1.4), we define for e = τ1 ∩ τ2

Z Z
1 1
we := w|τ1 ds = w|τ2 ds.
|e| e |e| e

Since ν u only changes

∂P the sign depending from which side an inner edge is seen, we have
e∈∂τ we e ∂ν u ds = 0 and thus
P R
τ ∈Th

X XZ
Rh [u](w) = (w − we ) ∂ν u ds.
τ ∈Th e∈∂τ e

Since the nodal interpolant satisfies Ih u ∈ Vh , we have that ∂ν Ih u is constant on each edge
and thus Z
(w − we )∂ν Ih u ds = 0
e
from which by the Cauchy-Schwarz inequality it follows
X XZ
|Rh [u](w)| = | (w − we ) ∂ν (u − Ih u) ds|
τ ∈Th e∈∂τ e
XX
≤ ∥w − we ∥L2 (e) ∥∇(u − Ih u)∥L2 (e) .
τ ∈Th e∈∂τ

For û ∈ H 2 (τ̂ ) we obtain from the trace theorem and Corollary 2.62 in Numerical Methods
for Partial Differential Equations on the reference element τ̂

∥∇(û − Îh u)∥L2 (∂ τ̂ ) ≤ ∥∇(û − Îh u)∥H 1/2 (∂ τ̂ ) ≤ c∥û − Îh u∥H 2 (τ̂ ) ≤ c′ |û|H 2 (τ̂ ) . (1.5)

Using the transformation formula, this implies for τ ∈ Th

∥∇(u − Ih u)∥L2 (∂τ ) ≤ ch1/2 |u|H 2 (τ ) .

In addition, the Bramble-Hilbert lemma (see Lemma 2.61 in Numerical Methods for Partial
Differential Equations) implies for ê ∈ ∂ τ̂ and w ∈ Π1

∥ŵ − ŵê ∥L2 (ê) ≤ c|ŵ|H 1 (τ̂ ) ,

because the left-hand side vanishes for constant functions. Again, the transformation formula
yields for e ∈ ∂τ and w ∈ Vh

∥w − we ∥L2 (e) ≤ ch1/2 |w|H 1 (τ ) .

6
 1.1 Non-Conforming Methods

From this we obtain using the Cauchy-Schwarz inequality for Euclidean scalar products
!1/2 !1/2
X X X
|Rh [u](w)| ≤ ch|u|H 2 (τ ) |w|H 1 (τ ) ≤ ch |u|2H 2 (τ ) |w|2H 1 (τ )
τ ∈Th τ ∈Th τ ∈Th

≤ ch|u|H 2 (Ω) ∥w∥Vh .

Remark. Inspecting the proof of the previous lemma, we see that the lemma also holds for
w ∈ Vh + H01 (Ω).

Keeping in mind that the piecewise linear conforming elements are contained in Vh , The-
orem 2.64 of Numerical Methods for Partial Differential Equations and the second Strang
lemma together with the previous lemma yield

∥u − uh ∥Vh ≤ c h |u|H 2 (Ω) ≤ c h ∥f ∥L2 (Ω) . (1.6)

Notice that the last estimate follows from the continuity of the inverse Laplacian ∆−1 :
L2 (Ω) → H01 (Ω) ∩ H 2 (Ω).

In order to measure the error w.r.t. the L2 -norm, we apply the Aubin-Nitsche lemma.

Theorem 1.6. Let Th be a family of non-degenerate triangulations of Ω ⊂ R2 . Furthermore,

let u ∈ H 2 (Ω) and f ∈ L2 (Ω). Then the finite-element error for the Poisson problem using
Crouzeix-Raviart elements satisfies

∥u − uh ∥L2 (Ω) + h∥u − uh ∥Vh ≤ c h2 |u|H 2 (Ω) .

Proof. In order to apply the Aubin-Nitsche lemma, we let V = H01 (Ω) and W = L2 (Ω). We
estimate the difference uφ − uφ,h of the solutions of the dual problems using (1.6):

∥uφ − uφ,h ∥Vh ≤ ch∥φ∥L2 (Ω) .

Furthermore, for the two additional terms in the Aubin-Nitsche lemma it follows from
Lemma 1.5

|ah (u − uh , uφ ) − (u − uh , φ)W | = |Rh [uφ ](u − uh )| ≤ ch|uφ |H 2 (Ω) ∥u − uh ∥Vh

≤ c′ h∥φ∥L2 (Ω) ∥u − uh ∥Vh

and

|ah (u, uφ − uφ,h ) − (f, uφ − uφ,h )W | = |Rh [u](uφ − uφ,h )| ≤ ch|u|H 2 (Ω) ∥uφ − uφ,h ∥Vh .

Hence we obtain
∥u − uh ∥L2 (Ω) ≤ ch∥u − uh ∥Vh + c′ h2 |u|H 2 (Ω)
and with (1.6) the assertion.

7
1 Non-Conforming and Mixed Methods

1.1.2 Polygonal Approximation of Curvilinear Boundaries

We have always assumed that Ω ⊂ R2 is polygonal. Domains Ω that are not polygonal
require curvilinear elements in the decomposition Th of Ω.

In the following discussion we confine ourselves to linear triangular elements

Vh := {v ∈ C(Ω) : v|τ is linear for every τ ∈ Th ,
v(x) = 0 for every vertex x ∈ ∂Ω of the decomposition}
for domains Ω ⊂ R2 . Then Vh ̸⊂ H01 (Ω). Due to Vh ⊂ H 1 (Ω) we do not need grid-dependent
norms. We can set ah = a and ℓh = ℓ.

Lemma 1.7. Let Ω be a C 2 domain and let Th be a family of non-degenerate decompositions.

Then
∥vh ∥L2 (∂Ω) ≤ c h3/2 |vh |H 1 (Ω) for all vh ∈ Vh .

Proof. Let τ ∈ Th be an element with curvilinear boundary Γτ := τ ∩ ∂Ω. We show

Z Z
vh ds ≤ chτ |∇vh |2 dx.
2 3
(1.7)
Γτ τ

The assertion follows from summation over all elements τ ∈ Th .

We may assume that the two vertices in Γτ are located on the x1 -axis. Let its coordinates be
(ξ, 0) and (ξ ′ , 0) and let the boundary Γτ be described by x2 = ϕ(x1 ). Due to ϕ(ξ) = ϕ(ξ ′ ) = 0
and |ξ − ξ ′ | ≤ hτ , we obtain1 for ξ ≤ x1 ≤ ξ ′
3
|ϕ(x1 )| ≤ c h2τ , c := max |ϕ′′ (x1 )|. (1.8)
2 ξ≤x1 ≤ξ′
Since vh ∈ Vh is linear on τ and thus vanishes on the x1 -axis, we have vh (x1 , x2 ) = b x2 and
hence |∇vh | = |b|. From |ϕ′ | ≤ c it follows that
Z Z ξ′ Z ξ′
vh ds = |b ϕ(x1 )| 1 + |ϕ (x1 )| dx1 ≤ c b hτ 1 dx1 = c′ b2 h5τ .
2 2
p ′ 2 4
′ 2
Γτ ξ ξ

The area of τ can be estimated from below by the area of its incircle. Hence,
h2τ 2
Z
b ≤ |∇vh |2 dx.
c2a τ

Comparing the previous two estimates yields (1.7).

Let u ∈ H01 (Ω) be the solution of the Poisson problem and let uh ∈ Vh denote the associated
weak solution, i.e.
a(uh , vh ) = (f, vh )L2 (Ω) for all vh ∈ Vh .
We estimate the finite-element error w.r.t. the H 1 -norm.
1
this follows from Taylor expansions for ϕ and ϕ′ taking into account the zeros of ϕ and ϕ′ .

8
 1.1 Non-Conforming Methods

Theorem 1.8. Let Ω be a C 2 domain and f ∈ L2 (Ω). Then for a non-degenerate decom-
position Th
∥u − uh ∥H 1 (Ω) ≤ c h ∥u∥H 2 (Ω) ≤ c′ h ∥f ∥L2 (Ω) .

Proof. The smoothness of the boundary shows u ∈ H 2 (Ω) ∩ H01 (Ω). Integration by parts
yields for vh ∈ Vh ⊂ H 1 (Ω)
Z
(f, vh )L2 (Ω) = (−∆u, vh )L2 (Ω) = a(u, vh ) − vh ∂ν u ds.
∂Ω

Using Cauchy-Schwarz inequality, the trace theorem, and Lemma 1.7, shows

|a(u, vh ) − (f, vh )L2 (Ω) | ≤ c∥∇u∥L2 (∂Ω) ∥vh ∥L2 (∂Ω) ≤ c h3/2 ∥u∥H 2 (Ω) ∥vh ∥H 1 (Ω) .

The assertion is a consequence of the second Strang lemma. Notice that the higher conver-
gence rate of the previous estimate is dominated by the approximation error

inf ∥u − vh ∥H 1 (Ω) ∼ h ∥u∥H 2 (Ω) .

vh ∈Vh

If we replace curvilinear edges in Th with secants then a polygonal approximation Ωh of the

original domain Ω is obtained. Correspondingly, we set Vh := S01,0 (Ωh ). Notice that due
to (1.8) an area τ ′ := τ ∩ (Ω \ Ωh ) or τ ′ := τ ∩ (Ωh \ Ω), resp., with

|τ ′ | ≤ c h |τ | (1.9)

is cut off. Therefore, the estimates of the previous theorem remain valid if a is replaced with
Z
ah (u, v) := ∇u · ∇v dx
Ωh

and ℓ with ℓh (v) = f v dx. This follows from d (1.9) with some hand waving as
R
Ωh
Z
|ah (u, vh ) − a(u, vh )| = |∇u · ∇vh | dx ≤ ch ∥u∥H 1 (Ω) ∥vh ∥H 1 (Ω)
Ω\Ωh

and Z
|ℓh (vh ) − ℓ(vh )| = | ∆u vh dx| ≤ ch ∥u∥H 2 (Ω) ∥vh ∥L2 (Ω)
Ω\Ωh

which leads to an error estimate w.r.t. the H 1 -norm as in Theorem 1.8.

Since Lemma 1.7 does not hold for vh ∈ Vh + H01 (Ω), we cannot apply the Aubin-Nitsche
lemma in order to obtain an L2 -estimate. Hence, we have to present a separate proof for
this situation.

9
1 Non-Conforming and Mixed Methods

Theorem 1.9. With the assumptions of Theorem 1.8 it holds that ∥u − uh ∥L2 (Ω) ≤
c h3/2 ∥u∥H 2 (Ω) .

Proof. Setting w := u − uh , let φ be the solution of

−∆φ = w in Ω, φ = 0 on ∂Ω.

Since Ω is assumed to be smooth, we have φ ∈ H 2 (Ω) ∩ H01 (Ω) and ∥φ∥H 2 (Ω) ≤ c∥w∥L2 (Ω) .
In contrast to the case of conforming elements, the Green formula shows boundary terms
due to w ̸∈ H01 (Ω)
Z
2
∥w∥L2 (Ω) = (w, −∆φ)L2 (Ω) = a(w, φ) − w ∂ν φ ds.
∂Ω

Since for arbitrary vh ∈ Vh from φ ∈ H01 (Ω)

Z
a(w, vh ) = a(u, vh ) − a(uh , vh ) = (−∆u, vh )L2 (Ω) + vh ∂ν u ds − (f, vh )L2 (Ω)
Z Z ∂Ω

= vh ∂ν u ds = (vh − φ) ∂ν u ds,
∂Ω ∂Ω

it follows that
Z Z
∥w∥2L2 (Ω) = a(w, φ − vh ) + (vh − φ) ∂ν u ds − w ∂ν φ ds. (1.10)
∂Ω ∂Ω

According to Theorem 1.8, the choice vh = Ih φ leads to

a(w, φ − vh ) ≤ cS ∥w∥H 1 (Ω) ∥φ − vh ∥H 1 (Ω) ≤ ch∥w∥H 1 (Ω) ∥φ∥H 2 (Ω)

≤ ch∥w∥H 1 (Ω) ∥w∥L2 (Ω) ≤ ch2 ∥u∥H 2 (Ω) ∥w∥L2 (Ω) .

In order to estimate the second term in (1.10), we need the approximation result ∥φ −
Ih φ∥L2 (∂Ω) ≤ c h3/2 ∥φ∥H 2 (Ω) , which can be proved using the trace theorem and the transfor-
mation formula; see (1.5). Applying the trace theorem to ∇u gives

|(∂ν u, φ − vh )L2 (∂Ω) | ≤ ∥∇u∥L2 (∂Ω) ∥φ − vh ∥L2 (∂Ω) ≤ c h3/2 ∥u∥H 2 (Ω) ∥φ∥H 2 (Ω)
≤ c′ h3/2 ∥u∥H 2 (Ω) ∥w∥L2 (Ω) .

The last term in (1.10) is estimated using Lemma 1.7, Theorem 1.8, and the trace theorem:

|(w, ∂ν φ)L2 (∂Ω) | ≤ ∥u − uh ∥L2 (∂Ω) ∥∇φ∥L2 (∂Ω) = ∥uh ∥L2 (∂Ω) ∥∇φ∥L2 (∂Ω)
≤ c h3/2 ∥uh ∥H 1 (Ω) ∥φ∥H 2 (Ω) ≤ c h3/2 ∥u∥H 1 (Ω) + ∥u − uh ∥H 1 (Ω) ∥w∥L2 (Ω)

≤ c′ h3/2 ∥u∥H 2 (Ω) ∥w∥L2 (Ω) .

In summary, we have

∥u − uh ∥2L2 (Ω) = ∥w∥2L2 (Ω) ≤ c h3/2 ∥u∥H 2 (Ω) ∥w∥L2 (Ω) ,

which leads to the assertion.

10
 1.2 Mixed Finite Elements

1.2 Mixed Finite Elements

As an example we consider the Poisson equation with homogeneous Dirichlet boundary

conditions
−∆u = f in Ω, u = 0 on ∂Ω.
This boundary value problem can be interpreted as a simple model for the displacement of
a clamped membrane under exterior load f . The numerical methods we have presented lead
to approximations uh of the displacement u. For technical applications the stress ∇u is at
least as important. An approximation to the stress can be obtain by differentiation of the
approximation uh . This will, however, lead to a loss of accuracy. In the following mixed
finite element method we introduce σ := ∇u ∈ R3 as a new variable, i.e. we consider the
first order system for [σ, u]

−div σ = f in Ω,
−∇u + σ = 0 in Ω,
u = 0 on ∂Ω

with the divergence operator div v := ∂1 v1 + ∂2 v2 + ∂3 v3 . The direct approximation

guarantees a higher accuracy than the indirect approach via differentiation of uh . In order
to derive a variational formulation, we integrate the first estimate

(σ, ∇v)L2 (Ω) = (f, v)L2 (Ω) for all v ∈ H01 (Ω), (1.11a)
−(∇u, τ )L2 (Ω) + (σ, τ )L2 (Ω) = 0 for all τ ∈ [L2 (Ω)]3 . (1.11b)

The solution is searched in the product space X := [L2 (Ω)]3 ×H01 (Ω), which can be equipped
with the norm
∥[σ, u]∥X := ∥σ∥L2 (Ω) + |u|H 1 (Ω) .
An equivalent formulation of (1.11) is the variational problem

[σ, u] ∈ X : b([σ, u], [τ, v]) = ℓ([τ, v]), [τ, v] ∈ X,

with the bilinear form b : X × X → R

b([σ, u], [τ, v]) := (σ, ∇v)L2 (Ω) − (∇u, τ )L2 (Ω) + (σ, τ )L2 (Ω)

and the linear form ℓ : X → R defined as ℓ([τ, v]) = (f, v)L2 (Ω) . From

|b([σ, u], [τ, v])| ≤ ∥σ∥L2 (Ω) |v|H 1 (Ω) + |u|H 1 (Ω) ∥τ ∥L2 (Ω) + ∥σ∥L2 (Ω) ∥τ ∥L2 (Ω)
≤ {∥σ∥L2 (Ω) + |u|H 1 (Ω) }{∥τ ∥L2 (Ω) + |v|H 1 (Ω) }
= ∥[σ, u]∥X ∥[τ, v]∥X

we see the b is continuous. However, the bilinear form is not coercive as

b([τ, v], [τ, v]) = (τ, ∇v)L2 (Ω) − (∇v, τ )L2 (Ω) + ∥τ ∥2L2 (Ω) = ∥τ ∥2L2 (Ω)

does not depend on v. Hence, the previous techniques cannot be applied to this kind of
problem.

11
1 Non-Conforming and Mixed Methods

1.2.1 The inf-sup Condition

In the following section, we present a theory for the solvability of general variational problems
u∈V : a(u, w) = ℓ(w), w ∈ W, (1.12)
with a continuous bilinear form a : V ×W → R and Hilbert spaces V, W . To this end, we will
use the following closed range theorem known from functional analysis. Let T : X → Y be a
linear operator between Banach spaces X and Y . Then the adjoint operator T ′ : Y ′ → X ′
of T is defined as
(T ′ φ)(x) = φ(T x), x ∈ X and φ ∈ Y ′ .
Additionally, we define the annihilator of a closed subspace U ⊂ Y ′
U 0 := {y ∈ Y : φ(y) = 0 for all φ ∈ U }.
The annihiliator should not be confounded with the orthogonal complement
U ⊥ = {v ∈ V : (v, u) = 0 for all u ∈ U }
of U ⊂ V .

Theorem 1.10 (Closed range theorem). Let T : X → Y be linear and continuous

between two Banach spaces X and Y . Then the following statements are equivalent:

(i) The image Im T is closed in Y .

(ii) It holds that Im T = (Ker T ′ )0 .

Using the previous theorem, the solvability of (1.12) can be investigated. To this end, define
the linear operator A : V → W ′ as
(Av)(w) = a(v, w), v ∈ V, w ∈ W.
The adjoint operator A′ of A satisfies (A′ w)(v) = a(v, w), v ∈ V , w ∈ W ′′ .

Theorem 1.11. The problem (1.12) has a unique solution u ∈ V for each ℓ ∈ W ′ if the
following conditions are satisfied:

(i) The bilinear form a is continuous, i.e. |a(v, w)| ≤ cS ∥v∥V ∥w∥W for all v ∈ V , w ∈ W ;

(ii) The inf-sup condition

a(v, w)
inf sup ≥ cE > 0 (1.13)
0̸=v∈V 0̸=w∈W ∥v∥V ∥w∥W

holds;

(iii) for every 0 ̸= w ∈ W there is v ∈ V such that a(v, w) ̸= 0.

In this case, we have Im A = W ′ and the solution u satisfies ∥u∥V ≤ c−1

E ∥ℓ∥W ′ .

12
 1.2 Mixed Finite Elements

Proof. The continuity of a carries over to A. Additionally, A is injective as Av1 = Av2

implies sup0̸=w∈W a(v1 − v2 , w) = 0 and (1.13) shows v1 − v2 = 0. Hence, for every φ ∈ Im A
there is a unique inverse A−1 φ. According to (1.13), we have ∥Av∥W ′ ≥ cE ∥v∥V for all v ∈ V
and thus the continuity of A−1 on Im A:
1
∥A−1 φ∥V ≤ ∥φ∥W ′ for all φ ∈ Im A. (1.14)
cE

Due to the continuity of A−1 , the preimage Im A of the closed set V is closed in W ′ . Theo-
rem 1.10 shows Im A = (Ker A′ )0 and with (iii) the reflexivity2 of Hilbert spaces shows

Ker A′ = {w ∈ W : a(v, w) = 0 for all v ∈ V } = {0}.

Hence, Im A = {0}0 = W ′ , which proves the surjectivity of A. For a given but arbitrary
ℓ ∈ W ′ we thus obtain the solution u := A−1 ℓ, which is unique due to the injectivity of A.
The stability estimate follows from (1.14).

The inf-sup condition is thus equivalent with the continuity of A−1 .

Corollary 1.12. If a is continuous then the following statements are equivalent.

(i) The inf-sup condition (1.13) holds;

(ii) The operator A : V → (Ker A′ )0 is an isomorphism with ∥Av∥W ′ ≥ cE ∥v∥V , v ∈ V ;

(iii) The operator A′ : (Ker A′ )⊥ → V ′ is an isomorphism with ∥A′ w∥V ′ ≥ cE ∥w∥W , w ∈

(Ker A′ )⊥ .

Proof. From the proof of Theorem 1.11 it is visible that (i) implies (ii).

If (ii) holds then ∥Av∥W ′ ≥ cE ∥v∥V for all v ∈ V implies (1.13). Additionally, for every q ∈
(Ker A′ )⊥ we define a functional φ : w 7→ (w, q)W . Then φ ∈ (Ker A′ )0 and ∥φ∥W ′ = ∥q∥W .
Since A is an isomorphism, there is p ∈ V with a(p, w) = (w, q)W for all w ∈ W . Hence,
∥q∥W = ∥φ∥W ′ = ∥Ap∥W ′ ≥ cE ∥p∥V and with a(p, q) = ∥q∥2W we obtain

a(v, q) a(p, q) ∥q∥2W

∥A′ q∥V ′ = sup ≥ = ≥ cE ∥q∥W .
0̸=v∈V ∥v∥V ∥p∥V ∥p∥V

Since ∥Av∥W ′ > 0 for all v ̸= 0, A′ : (Ker A′ )⊥ → V ′ satisfies the three assumptions of
Theorem 1.11. Conversely, if A′ : (Ker A′ )⊥ → V ′ is an isomorphism then for given v ∈ V

φ(v) (A′ w)(v)

∥v∥V = sup = sup ′
0̸=φ∈V ′ ∥φ∥V ′ 0̸=w∈(Ker A′ )⊥ ∥A w∥V ′

a(v, w) a(v, w)
= sup ′
≤ c−1
E sup
0̸=w∈(Ker A′ )⊥ ∥A w∥V ′ 0̸=w∈W ∥w∥W

and we obtain (i).

2
A Banach space V is called reflexive if V ′′ = V .

13
1 Non-Conforming and Mixed Methods

Remark. Theorem 1.11 is a generalization of the Lax-Milgram theorem as in the case

V = W and a coercive bilinear form we obtain the inf-sup condition
a(v, w) a(v, v)
inf sup ≥ inf ≥ cK .
0̸=v∈V 0̸=w∈V ∥v∥V ∥w∥V 0̸=v∈V ∥v∥2 V

Condition (iii) in Theorem 1.11 follows from choosing v = w and applying the coercivity.

In the case of positive-semidefinite bilinear forms the inf-sup condition and coercivity are
equivalent.

Lemma 1.13. Let a : V × V → R be a symmetric and continuous bilinear form satisfying

a(v, v) ≥ 0 for all v ∈ V . If a satisfies the inf-sup condition

a(v, w)
sup ≥ cE ∥v∥V , v∈V
0̸=w∈V ∥w∥V

then a is coercive with cK = c2E /cS .

Proof. Due to the assumptions a is a semi-scalar product. For v ∈ V we obtain from the
Cauchy-Schwarz inequality

a2 (v, w) a(v, v) a(w, w) a(v, v) cS ∥w∥2V

c2E ∥v∥2V ≤ sup ≤ sup ≤ sup = cS a(v, v).
0̸=w∈V ∥w∥2V w∈V ∥w∥2V w∈V ∥w∥2V

Remark. Let a : V × W → R be a continuous bilinear form and T : V → W be the linear

operator defined by (T v, w)W = a(v, w) for all v ∈ V and w ∈ W . If a satisfies the inf-sup
condition then
∥T v∥ (T v, w) a(v, w)
inf = inf sup = inf sup ≤ cE
0̸=v∈V ∥v∥ 0̸=v∈V 0̸=w∈W ∥v∥∥w∥ 0̸=v∈V 0̸=w∈W ∥v∥∥w∥

and
∥T v∥ (T v, w) a(v, w)
sup = sup sup = sup sup ≤ cS .
0̸=v∈V ∥v∥ 0̸=v∈V 0̸=w∈W ∥v∥∥w∥ 0̸=v∈V 0̸=w∈W ∥v∥∥w∥

Hence, for every singular value σ of T , i.e. T ∗ T v = σ 2 v with some vector v ̸= 0 we have

(T v, T v) (T ∗ T v, v) (T ∗ T v, v) ∥T v∥2
c2E ≤ inf = inf ≤ σ 2 ≤ sup = sup 2
≤ c2S .
0̸=v∈V (v, v) 0̸=v∈V (v, v) 0̸=v∈V (v, v) 0̸=v∈V ∥v∥

Thus cE ≤ σ ≤ cS .

If V = W and a is coercive and continuous then for every eigenvalue λ of T , i.e. T v = λv

with some vector v ̸= 0 we observe from
a(v, v) (T v, v)
cK ≤ = =λ
∥v∥2 ∥v∥2

14
 1.3 Saddle Point Problems and Constrained Variational Problems

and
a(v, v)
λ= ≤ cS
∥v∥2
that cK ≤ λ ≤ cS .

For the variational problem (1.12) we can prove a Céa-type convergence result for the solution
of the discrete problem

uh ∈ Vh : a(uh , wh ) = ℓ(wh ), wh ∈ Wh , (1.15)

with conforming discretizations Vh ⊂ V and Wh ⊂ W .

Theorem 1.14. Let the assumptions of Theorem 1.11 be valid for both the Hilbert
spaces V, W and the conforming finite element spaces Vh ⊂ V and Wh ⊂ W . Then the
unique solution uh ∈ Vh satisfies
 
cS
∥u − uh ∥V ≤ 1 + inf ∥u − vh ∥V .
cE vh ∈Vh

Proof. The linear functional φ : Wh → R defined by φ(w) := a(u − vh , w) satisfies ∥φ∥W ′ ≤

cS ∥u − vh ∥V . Furthermore, for Ah : Vh → Wh′ defined by (Ah vh )(wh ) := a(vh , wh ) for vh ∈ Vh
and wh ∈ Wh we obtain that ∥A−1 h ∥ ≤ 1/cE . From the Galerkin orthogonality it follows that

φ(wh ) = a(u − vh , wh ) = a(uh − vh , wh ) = (Ah (uh − vh ))(wh ), wh ∈ Wh .

Thus we have φ = Ah (uh − vh ) or equivalently uh − vh = A−1

h φ. This shows

1 cS
∥uh − vh ∥V ≤ ∥φ∥W ′ ≤ ∥u − vh ∥V
cE cE
and using the triangle inequality the assertion.

1.3 Saddle Point Problems and Constrained Variational

Problems

Let a(v, w) = (v, w)L2 (Ω) , b(v, q) = −(v, ∇q)L2 (Ω) and V = [L2 (Ω)]3 , W = H01 (Ω) then (1.11)
takes the form: find (u, p) ∈ V × W such that

a(u, v) + b(v, p) = f (v) for all v ∈ V, (1.16a)

b(u, q) = g(q) for all q ∈ W. (1.16b)

Here, f ∈ V ′ , g ∈ W ′ and a : V × V → R, b : V × W → R are continuous bilinear forms

with continuity constants cS and c′S , resp. Problems of this kind are called saddle point
problems or mixed variational problems.

15
1 Non-Conforming and Mixed Methods

For the solvabilty of saddle point problems we investigate the constrained variational
problem
u ∈ G : a(u, z) = f (z) for all z ∈ Z, (1.17)
where the set
G := {v ∈ V : b(v, q) = g(q) for all q ∈ W }
incorporates the constraint and Z := {v ∈ V : b(v, q) = 0 for all q ∈ W }.

Theorem 1.15. Let a be coercive on Z. If G ̸= ∅ then (1.17) has a unique solution u ∈ G

and u satisfies  
1 cS
∥u∥V ≤ ∥f ∥V ′ + + 1 sup ∥v∥V .
cK cK v∈G

Proof. Let v ∈ G. Due to the bilinearity of b we have v + z ∈ G for all z ∈ Z. On the other
hand, u ∈ G can be represented as
u = v + z̃ (1.18)
with some z̃ ∈ Z. The problem (1.17) is thus equivalent with the problem: find z̃ ∈ Z such
that
a(z̃, z) = f (z) − a(v, z) for all z ∈ Z. (1.19)
Since a is assumed to be coercive, we can apply the Lax-Milgram theorem. Hence, (1.19)
and therefore also (1.17) are uniquely solvable. Finally, (1.19) implies for z = z̃
1
∥z̃∥V ≤ (∥f ∥V ′ + cS ∥v∥V ). (1.20)
cK

Using (1.18) and the triangle inequality leads to the assertion.

Example 1.16. We investigate the Poisson boundary value problem with inhomogeneous
boundary condition
−∆u = f in Ω, u = g on ∂Ω
with f ∈ H −1 (Ω) and g ∈ H 1/2 (∂Ω). We set V = H 1 (Ω) and W = H −1/2 (∂Ω) and choose
Z
a(v, w) = ∇v · ∇w dx, f (v) = ⟨v, f ⟩Ω
Ω
b(v, q) = ⟨v, q⟩∂Ω , g(q) = ⟨g, q⟩∂Ω

with the duality products ⟨·, ·⟩Ω and ⟨·, ·⟩∂Ω . With this we obtain

G = {v ∈ H 1 (Ω) : ⟨v, q⟩∂Ω = ⟨g, q⟩∂Ω for all q ∈ H −1/2 (∂Ω)}

and Z = H01 (Ω). According to Theorem 1.15, the boundary value problem with inhomoge-
neous boundary conditions is uniquely solvable.

The next step is to investigate the connection of the solvability of (1.16) and (1.17).

16
 1.3 Saddle Point Problems and Constrained Variational Problems

Lemma 1.17. If (u, p) ∈ V × W is a solution of the saddle point problem (1.16) then u
solves the constrained variational problem (1.17).

Proof. The second part of (1.16) is equivalent with u ∈ G. For z ∈ Z ⊂ V we have b(z, p) = 0
due to p ∈ W . Hence, the first part of (1.16) yields
a(u, z) = f (z) for all z ∈ Z.
Therefore, u is a solution of (1.17).

In the following discussion, we try to revert the previous lemma, i.e. it will be shown that
for every solution u of (1.17) there is p ∈ W such that (u, p) solves the saddle point prob-
lem (1.16). The orthogonal complement Z ⊥ of Z is a closed subspace of V . Hence, the
space V can be decomposed into a direct sum
V = Z ⊕ Z ⊥.
Due to (1.17), (u, p) ∈ V × W solves (1.16) iff
b(v, p) = f (v) − a(u, v) for all v ∈ Z ⊥ . (1.21)
For the solvability of the previous equation w.r.t. p we assume the following LBB condition
for b, which is named after Ladyshenskaja, Babuška, and Brezzi.

Lemma 1.18. Let the bilinear form b satisfy the LBB condition

b(v, q)
sup ≥ cLBB ∥q∥W for all q ∈ Y ⊥ , (1.22)
0̸=v∈V ∥v∥V

where Y := {q ∈ W : b(v, q) = 0 for all v ∈ V }. Then (1.21) has a solution p ∈ W .

Proof. Due the bilinearity and the continuity of b, by (Bv)(q) := b(v, q) we define a contin-
uous linear operator B : V → W ′ . Corollary 1.12 (iii) applied to A = B ′ and V = Y ⊥ shows
due to Ker A′ = Z and the LBB condition that
∥Bv∥W ′ ≥ cLBB ∥v∥V for all v ∈ Z ⊥ . (1.23)
Denote by R : W ′ → W the Riesz represention operator and define a symmetric and
continuous bilinear form d : V × V → R via
d(w, v) := (Bv)(RBw) = (RBv, RBw)W for all v, w ∈ V.
From (1.23) it follows that
d(v, v) = ∥RBv∥2W = ∥Bv∥2W ′ ≥ c2LBB ∥v∥2V for all v ∈ Z ⊥ .
Hence, d is Z ⊥ -coercive and the Lax-Milgram theorem shows the existence of y ∈ Z ⊥ satis-
fying
b(v, RBy) = d(y, v) = f (v) − a(u, v) for all v ∈ Z ⊥ .
The element p := RBy ∈ W satisfies (1.21).

17
1 Non-Conforming and Mixed Methods

We summarize the previous results.

Theorem 1.19. Let a be coercive on Z. Let b satisfy the LBB condition (1.22). If G ̸= ∅
then the saddle point problem (1.16) has a solution (u, p) ∈ V × W . Its first component
u ∈ V is unique and we have the estimates
 
1 1 cS
∥u∥V ≤ ∥f ∥V ′ + + 1 ∥g∥W ′
cK cLBB cK

and   
1 cS cS
inf ∥p + y∥W ≤ +1 ∥f ∥V ′ + ∥g∥W ′ .
y∈Y cLBB cK cLBB

Proof. According to Theorem 1.15, the constrained variational problem (1.17) has a unique
solution u ∈ G. Lemma 1.18 shows the existence of a component p ∈ W such that (u, p) ∈
V × W solves the saddle point problem (1.16). Lemma 1.17 states that the V -component of
every solution of (1.16) also solves (1.17). Hence, the component u is uniquely determined.

The element u ∈ V can be uniquely represented as u = v + z̃ with z̃ ∈ Z, v ∈ Z ⊥ . The

bilinearity of b, (1.16), and (1.23) show that

∥g∥W ′ = ∥Bu∥W ′ = ∥Bv∥W ′ ≥ cLBB ∥v∥V .

Hence, ∥v∥V ≤ c−1

LBB ∥g∥W ′ . Since with u ∈ G it also holds that v ∈ G, from (1.20) we obtain
the estimate
 
1 1 cS
∥u∥V ≤ ∥v∥V + ∥z̃∥V ≤ ∥f ∥V ′ + + 1 ∥g∥W ′ .
cK cLBB cK

We investigate the second component p ∈ W . From (1.21) we obtain

|b(v, p)| ≤ |f (v)| + |a(u, v)| ≤ (∥f ∥V ′ + cS ∥u∥V )∥v∥V for all v ∈ Z ⊥

and thus
b(v, p)
sup ≤ ∥f ∥V ′ + cS ∥u∥V .
0̸=v∈V ∥v∥V
Since y ∈ Y exists such that p + y ∈ Y ⊥ , we obtain from (1.22)

b(v, p + y) b(v, p)
cLBB ∥p + y∥W ≤ sup = sup ≤ ∥f ∥V ′ + cS ∥u∥V .
0̸=v∈V ∥v∥V 0̸=v∈V ∥v∥V

Taking into account the first estimate, we obtain the assertion.

Remark. If Y = {0} holds then from the proof of Lemma 1.18 it is visible that the W -
component and thus the complete solution (u, p) ∈ V × W of (1.16) is unique. Then the
two estimates in the previous theorem can be viewed as stabilty estimates on the influence
of perturbations of f and g on the solution (u, p).

18
 1.3 Saddle Point Problems and Constrained Variational Problems

Saddle points and optimality criteria

If the bilinear form a is symmetric and a(v, v) > 0 for all 0 ̸= v ∈ Z then u ∈ G solves the
constrained variational problem (1.17) iff it solves the minimization problem
1
minimize J(v) := a(v, v) − f (v) over v ∈ G; (1.24)
2
see the exercises. We define the Lagrange functional associated with this constrained mini-
mization problem

L(v, q) := J(v) + b(v, q) − g(q) for all v ∈ V, q ∈ W.

Using the Lagrange functional L, a sufficient optimality criterion for (1.24) in saddle point
form can be presented. Here, (u, p) ∈ V × W is called saddle point of L if

L(u, q) ≤ L(u, p) ≤ L(v, p) for all v ∈ V, q ∈ W. (1.25)

Lemma 1.20. Let (u, p) ∈ V × W be a saddle point of L. Then the component u ∈ V

solves the minimization problem (1.24).

Proof. First we show that u ∈ G is valid. From the left part of (1.25) it follows that

b(u, q) − g(q) ≤ b(u, p) − g(p) for all q ∈ W. (1.26)

Since W is a linear space, we have q := p + y ∈ W for all y ∈ W . Taking into account the
linearity of b(u, ·) and g, (1.26) yields the estimate

b(u, y) − g(y) ≤ 0 for all y ∈ W.

With y ∈ W we also have −y ∈ W . Hence, b(u, y) − g(y) = 0 for all y ∈ W and thus u ∈ G.
The right part of (1.25) shows for all v ∈ G

J(u) = J(u) + b(u, p) − g(p) = L(u, p)

≤ L(v, p) = J(v) + b(v, p) − g(p) = J(v).

Hence, u solves the variational problem (1.24).

In the lecture Numerical Methods for Partial Differential Equations we have shown for
symmetric bilinear forms that the minimum of the energy functional J coincides with the
solution of the variational problem. The following theorem states that the solution of (1.16)
is the saddle point (the minimum for v = u and the maximum for q = p).

Theorem 1.21. Let a(v, v) ≥ 0 for v ∈ V . (u, p) ∈ V × W is a solution of (1.16) iff (u, p)
is a saddle point.

Proof. see the exercises.

19
1 Non-Conforming and Mixed Methods

1.3.1 Conforming Approximation of Mixed Variational Problems

We investigate a conforming finite element discretization of the mixed variational prob-

lem (1.16). Accordingly, let Vh ⊂ V and Wh ⊂ W . We search for a solution (uh , ph ) ∈ Vh ×Wh
of the discrete mixed variational problem

a(uh , vh ) + b(vh , ph ) = f (vh ) for all vh ∈ Vh , (1.27a)

b(uh , qh ) = g(qh ) for all qh ∈ Wh . (1.27b)

This finite element discretization of (1.16) is called method of mixed finite elements.
The solvability and stability of the solutions of (1.27) will be investigated analogous to the
continuous case. Let

Gh := {vh ∈ Vh : b(vh , qh ) = g(qh ) for all qh ∈ Wh }

and Zh := {vh ∈ Vh : b(vh , qh ) = 0 for all qh ∈ Wh }. Note that despite the assumption
Vh ⊂ V it holds that
Gh ̸⊂ G and Zh ⊂ ̸ Z
in general. Hence, the Z-coercivity of a does not automatically imply the Zh -coercivity. We
assume that cK > 0 exists such that

a(vh , vh ) ≥ cK ∥vh ∥2V for all vh ∈ Zh . (1.28)

Furthermore, we assume that there is a constant cLBB > 0 such that

b(vh , qh )
sup ≥ cLBB ∥qh ∥W for all qh ∈ Yh⊥ (1.29)
0̸=vh ∈Vh ∥vh ∥V

with Yh := {qh ∈ Wh : b(vh , qh ) = 0 for all vh ∈ Vh }. With these assumptions Theorem 1.19
can be carried over to the discrete problem (1.27). We obtain

Theorem 1.22. Let Gh ̸= ∅. Furthermore, let a and b satisfy the conditions (1.28) and
(1.29). Then (1.27) has a solution (uh , ph ) ∈ Vh ×Wh . Its first component uh ∈ Vh is uniquely
determined and the following estimates hold
 
1 1 cS
∥uh ∥V ≤ ∥f ∥V ′ + + 1 ∥g∥W ′
cK cLBB cK

and   
1 cS cS
inf ∥ph + yh ∥W ≤ +1 ∥f ∥V ′ + ∥g∥W ′ .
yh ∈Yh cLBB cK cLBB

For the sake of simplicity we investigate the convergence of the method of mixed finite
elements in the case Yh = {0}. In this case, the convergence analysis of Galerkin methods
can be applied.

20
 1.3 Saddle Point Problems and Constrained Variational Problems

Theorem 1.23. Let the bilinear forms a and b satisfy the assumptions of Theorems 1.19
and 1.22. Furthermore, let the discrete LBB condition (1.29) be valid with Yh = {0}, i.e.

b(vh , qh )
sup ≥ c̃LBB ∥qh ∥W for all qh ∈ Wh . (1.30)
0̸=vh ∈Vh ∥vh ∥V

Then the difference of the solution (uh , ph ) of (1.27) and the solution (u, p) of (1.16) satisfies
 
∥u − uh ∥V + ∥p − ph ∥W ≤ c inf ∥u − vh ∥V + inf ∥p − qh ∥W
vh ∈Vh qh ∈Wh

with a constant c > 0.

Proof. For arbitrary ṽh ∈ Vh , q̃h ∈ Wh we obtain

a(uh − ṽh , vh ) + b(vh , ph − q̃h ) = a(u − ṽh , vh ) + b(vh , p − q̃h ) =: f˜(vh ) for all vh ∈ Vh ,
b(uh − ṽh , qh ) = b(u − ṽh , qh ) =: g̃(qh ) for all qh ∈ Wh .
Taking into account Yh = {0}, Theorem 1.22 applied to the previous saddle point problem
for the errors uh − ṽh and ph − q̃h with right-hand sides f˜ and g̃ implies
c′S c′
  
cS cS
∥uh − ṽh ∥ ≤ + +1 ∥u − ṽh ∥ + S ∥p − q̃h ∥
cK c̃LBB cK cK
and
c′S
    
1 cS ′
∥ph − q̃h ∥ ≤ +1 cS ∥p − q̃h ∥ + cS 1 + ∥u − ṽh ∥ .
c̃LBB cK cLBB
The triangle inequalities ∥u − uh ∥ ≤ ∥u − ṽh ∥ + ∥uh − ṽh ∥ and ∥p − ph ∥ ≤ ∥p − q̃h ∥ + ∥ph − q̃h ∥
lead to the assertion.

Checking the validity of the discrete LBB condition (1.30) directly is usually technically
difficult. The following criterion gives an easy sufficient condition.

Lemma 1.24 (Fortin criterion). Let the bilinear form b satisfy the LBB condition (1.22)
with Y = {0}. Let Vh × Wh ⊂ V × W be a family of subspaces with projectors Πh : V → Vh
satisfying
b(v − Πh v, qh ) = 0 for all v ∈ V and qh ∈ Wh
and ∥Πh v∥V ≤ cΠ ∥v∥V for all v ∈ V with an h-independent constant cΠ > 0. Then b satisfies
also the discrete LBB condition (1.30).

Proof. Since Πh v ∈ Vh from the assumptions we have for all qh ∈ Wh ⊂ W

b(v, qh ) b(Πh v, qh ) b(Πh v, qh )
cLBB ∥qh ∥W ≤ sup = sup ≤ cΠ sup
0̸=v∈V ∥v∥V 0̸=v∈V ∥v∥V 0̸=v∈V ∥Πh v∥V
b(vh , qh )
≤ cΠ sup .
0̸=vh ∈Vh ∥vh ∥V

This shows the assertion with c̃LBB := cLBB /cΠ .

21
1 Non-Conforming and Mixed Methods

1.4 Solvability of the Mixed Formulation of the Poisson

Problem

We have started the chapter with a mixed formulation (1.11)

(u, v)L2 (Ω) − (v, ∇p)L2 (Ω) = 0 for all v ∈ [L2 (Ω)]3 , (1.31a)
(u, ∇q)L2 (Ω) = (f, q)L2 (Ω) for all q ∈ H01 (Ω). (1.31b)

of the Poisson problem. This is also called primal mixed formulation. After the general
theory on existence and uniqueness of mixed variational problems we will check whether (1.31)
and a second formulation satisfy the assumptions of this theory.

Theorem 1.25. The primal mixed formulation (1.31) has a unique solution (u, p) ∈
[L2 (Ω)]3 × H01 (Ω) for every f ∈ H −1 (Ω).

Proof. We set V = [L2 (Ω)]3 , W = H01 (Ω) (equipped with the norm | · |H 1 (Ω) ) and

a(v, w) := (v, w)L2 (Ω) , b(v, q) := −(v, ∇q)L2 (Ω) .

Then a is coercive on the whole space V . In order to check the LBB condition, we choose
v := −∇q ∈ V with q ∈ W . For q ̸= 0 we have v ̸= 0. To see this, suppose that v = 0, Then
q would be constant, which due to q ∈ H01 (Ω) implies that q = 0. Hence, we obtain

b(v, q) b(−∇q, q) (∇q, ∇q)L2 (Ω)

inf sup ≥ inf = inf = 1.
0̸=q∈W 0̸=v∈V ∥v∥V ∥q∥W 0̸=q∈W ∥∇q∥V ∥q∥W 0̸=q∈W |q|2H 1 (Ω)

As a consequence, the assumptions of Theorem 1.19 are satisfies and the assertion follows.

In order to be able to introduce the second formulation of the Poisson problem, we introduce
the Hilbert space
H(div; Ω) = {v ∈ [L2 (Ω)]3 : div v ∈ L2 (Ω)}
with the scalar product (v, w)H(div;Ω) := (v, w)L2 (Ω) + (div v, div w)L2 (Ω) . Here, the divergence
div v of v is defined by

(div v, q)L2 (Ω) = −(v, ∇q)L2 (Ω) for all q ∈ C0∞ (Ω).

The following theorem guarantees the existence of a trace operator for H(div; Ω) functions.

Theorem 1.26. Let Ω ⊂ R3 be a bounded Lipschitz domain with unit outer normal ν. Then
the mapping γν : [C ∞ (Ω)]3 → C ∞ (∂Ω) with γν v = ν · v|∂Ω can be extended to a continuous
linear operator γν : H(div; Ω) → H −1/2 (∂Ω) and the Green’s identity for v ∈ H(div; Ω) and
w ∈ H 1 (Ω) reads Z Z
⟨γν v, γ0 w⟩∂Ω = w div v dx + v · ∇w dx.
Ω Ω

22
 1.4 Solvability of the Mixed Formulation of the Poisson Problem

The following theorem is important for imposing boundary conditions.

Theorem 1.27. The trace operator γν : H(div; Ω) → H −1/2 (∂Ω) is surjective and has a
continuous right inverse, i.e., for any g ∈ H −1/2 (∂Ω) there is v ∈ H(div; Ω) such that γν v = g
and ∥v∥H(div;Ω) ≤ c∥g∥H −1/2 (∂Ω) .

In addition to the primal mixed formulation (1.31), using

−(v, ∇q)L2 (Ω) = (div v, q)L2 (Ω) for v ∈ H(div; Ω) and q ∈ H01 (Ω),

the dual mixed formulation of the Poisson problem on H(div; Ω) can be introduced: find
(u, p) ∈ H(div; Ω) × L2 (Ω) such that

(u, v)L2 (Ω) + (div v, p)L2 (Ω) = 0 for all v ∈ H(div; Ω), (1.32a)
(div u, q)L2 (Ω) = −(f, q)L2 (Ω) for all q ∈ L2 (Ω). (1.32b)

Compared to the primal formulation, the dual formulation leads to a more efficient dis-
cretization, while the former offers some advantages when applied to, e.g., problems from
linear elasticity.

Interestingly the Dirichlet condition p|∂Ω = 0 does not appear explicitly. The following
theorem shows that it is contained implicitly in the formulation.

Theorem 1.28. The dual mixed formulation (1.32) has a unique solution (u, p) ∈
H(div; Ω) × L2 (Ω) for every f ∈ H −1 (Ω). We even have p ∈ H01 (Ω).

Proof. The problem (1.32) is a saddle point problem for the Hilbert spaces V = H(div; Ω),
W = L2 (Ω) and the bilinear forms

a(v, w) := (v, w)L2 (Ω) , b(v, q) := (div v, q)L2 (Ω) .

The kernel of b
Z := {v ∈ V : (div v, q)L2 (Ω) = 0 for all q ∈ L2 (Ω)}
coincides with the space of divergence-free functions. The bilinear form a is continuous
on V × V and coercive on Z, because for v ∈ Z we have

a(v, v) = ∥v∥2L2 (Ω) = ∥v∥2L2 (Ω) + ∥div v∥2L2 (Ω) = ∥v∥2H(div;Ω) .

It remains to show the LBB condition for b. Let 0 ̸= q ∈ L2 (Ω) be arbitrary. Since C0∞ (Ω)
is dense in L2 (Ω), there is 0 ̸= q̃ ∈ C0∞ (Ω) such that
1
∥q − q̃∥2L2 (Ω) ≤ ∥q∥2L2 (Ω) .
2
We define v having the components
Z x1
v1 (x1 , . . . , xd ) := q̃(t, x2 , . . . , xd ) dt
−∞

23
1 Non-Conforming and Mixed Methods

and vi = 0, i ≥ 2. Then v satisfies pointwise

∂v1
div v = = q̃.
∂x1

Hence v ∈ V \ {0} and in the same way as in the proof of the Poincaré-Friedrichs inequality
we obtain ∥v∥L2 (Ω) ≤ c∥q̃∥L2 (Ω) , which implies

∥v∥2H(div;Ω) = ∥v∥2L2 (Ω) + ∥div v∥2L2 (Ω) ≤ (c2 + 1) ∥q̃∥2L2 (Ω) .

From
1 2  1 
(q̃, q)L2 (Ω) = ∥q̃∥L2 (Ω) + ∥q∥2L2 (Ω) − ∥q − q̃∥2L2 (Ω) ≥ ∥q̃∥2L2 (Ω) + ∥q∥2L2 (Ω)
2 4
1
≥ ∥q̃∥L2 (Ω) ∥q∥L2 (Ω)
2
(using Young’s inequality 2ab ≤ a2 + b2 ⇔ (a − b)2 ≥ 0 for a, b ∈ R) it follows that

b(v, q) (div v, q)L2 (Ω) (q̃, q)L2 (Ω) ∥q∥L2 (Ω)

≥√ =√ ≥ √ .
∥v∥H(div;Ω) c2 + 1 ∥q̃∥L2 (Ω) c2 + 1 ∥q̃∥L2 (Ω) 2 c2 + 1

Theorem 1.19 shows that there is a unique solution (u, p) ∈ H(div; Ω) × L2 (Ω).

We show that p ∈ H01 (Ω). Due to [C0∞ (Ω)]3 ⊂ V , (1.32) implies

(u, φ)L2 (Ω) = −(p, div φ)L2 (Ω) for all φ ∈ [C0∞ (Ω)]3 .

The definition of the weak derivative shows u = ∇p and thus p ∈ H 1 (Ω). Furthermore,
from (1.32) we obtain for φ ∈ [C ∞ (Ω)]3 ⊂ V
Z
0 = (u, φ)L2 (Ω) + (div φ, p)L2 (Ω) = (∇p, φ)L2 (Ω) + (div φ, p)L2 (Ω) = φ · ν p ds
∂Ω

and thus p ∈ H01 (Ω).

1.4.1 The Raviart-Thomas Element

Let Ω = Ω1 ∪ Ω2 with Σ := Ω1 ∩ Ω2 and u|Ωi ∈ H(div; Ωi ). In the exercises it will be shown

that u ∈ H(div; Ω) iff the normal component u · ν is continuous across the interface Σ. In
the following, we present a conforming discretization of V = H(div; Ω) in two dimensions.
To this end, let Th be a triangulation of Ω and
   
aτ x
2 2
Vh := {vh ∈ [L (Ω)] : vh (x, y)|τ = + cτ with aτ , bτ , cτ ∈ R
bτ y
and vh · ν is continuous at ∂τ for all τ ∈ Th }.

Raviart-Thomas functions vh ∈ Vh obviously satisfy vh |τ ∈ [Π1 ]2 (each component is a

univariate linear polynomial) and has three degrees of freedom per element τ ∈ Th . The

24
 1.4 Solvability of the Mixed Formulation of the Poisson Problem

normal component of Raviart-Thomas functions is constant on every edge e ∈ ∂τ . To see

this, let e be given by {α + tβ, t ∈ [0, 1]}. Then νe = [−β2 , β1 ]T /∥β∥2 and
     
aτ aτ
νe · vh (α + tβ) = νe · + cτ (α + tβ) = νe · + cτ α .
bτ bτ

The three degrees of freedom aτ , bτ , and cτ of vh on τ are determined by the normal com-
ponents νe · vh on the three edges e ∈ ∂τ . To see this, the two components of vh (xi ) in each
of the three vertices xi , i = 1, 2, 3, of τ can be computed from the normal components of the
two adjacent edges. For the six values vh (xi ), i = 1, 2, 3, a unique linear bivariate polynomial
vh |τ ∈ [Π21 ]2 can be found. As an element of [Π21 ]2 , all normal components of vh |τ are linear
on the edges. They are even constant because their values coincide on the two endpoints.
The linear ansatz    
aτ x
vh (x, y)|τ = + Bτ
bτ y
with Bτ ∈ R2×2 together with νe · vh (α + tβ) = ce yields that Bτ = cτ I.

Figure 1.2: Raviart-Thomas element

Interpolation operator. Using Theorem 1.23, the finite element error can be estimated
by the approximation error. This, in turn, can be estimated as before by the interpolation
error. We define Iτ : H 1 (τ ) → Vh |τ by
Z
(v − Iτ v) · νe ds = 0 for all e ∈ ∂τ
e

and we define Ih : H 1 (Ω) → Vh by

(Ih v)|τ = Iτ (v|τ ) for all τ ∈ Th .

Hence, the interpolation is defined by the mean values of the normal components of v ∈
H 1 (Ω) over the edges. The Gauß divergence theorem shows
Z XZ
div(v − Iτ v) dx = (v − Iτ v) · νe ds = 0. (1.33)
τ e∈∂τ e

Since Iτ v is linear and thus div(Iτ v) is constant, we obtain the property

Z
1
div(Iτ v) = div v dx. (1.34)
|τ | τ

25
1 Non-Conforming and Mixed Methods

Discrete problem. According to Theorem 1.28, the dual mixed formulation of the Poisson
problem is uniquely solvable in H(div; Ω) × L2 (Ω). It remains to check whether the unique
solvability is valid also for the discrete space Vh × Wh , where

Wh := {q ∈ L2 (Ω) : q|τ ∈ Π0 for all τ ∈ Th }.

Employing the spaces Vh and Wh , we can investigate the discrete variational problem asso-
ciated with (1.32): find (uh , ph ) ∈ Vh × Wh such that

(uh , vh )L2 (Ω) + (div vh , ph )L2 (Ω) = 0 for all vh ∈ Vh , (1.35a)

(div uh , qh )L2 (Ω) = −(f, qh )L2 (Ω) for all qh ∈ Wh . (1.35b)

Theorem 1.29. The discretization (1.35) using Raviart-Thomas elements has a unique
solution (uh , ph ) ∈ Vh × Wh for every right-hand side f ∈ L2 (Ω).

Proof. For vh ∈ Zh = {vh ∈ Vh : (div vh , qh )L2 (Ω) = 0 for all qh ∈ Wh } it holds that
div vh dx = 0 for all τ ∈ Th . Furthermore, (div vh )|τ ∈ Π0 shows that div vh = 0 in τ
R
τ
and Vh ⊂ V implies div vh = 0. Hence, ∥vh ∥2H(div;Ω) = ∥vh ∥2L2 (Ω) = a(vh , vh ) and a is coer-
cive on Zh . It remains to show the LBB condition for b. According to the Fortin criterion
(Lemma 1.24), we have to construct a projector Πh : V → Vh with h-independent continuity
constant and

b(v − Πh v, qh ) = (div(v − Πh v), qh )L2 (Ω) = 0 for all v ∈ V and qh ∈ Wh .

The last condition is satisfied by Πh := Ih due to (1.33).

By transformation to the reference element it is readily shown that

∥Iτ v∥L2 (τ ) ≤ c∥v∥L2 (τ ) for all τ ∈ Th .

The continuity of Πh follows from (1.34) as

∥Πh v∥2H(div;τ ) = ∥Ih v∥2L2 (τ ) + ∥div Ih v∥2L2 (τ ) ≤ c∥v∥2L2 (τ ) + ∥div v∥2L2 (τ ) ≤ c∥v∥2H(div;τ ) .

The assertion follows from summation over τ ∈ Th .

After the solvability of both the continuous and the discrete problem have been investigated,
error estimates between the two solutions are to be derived. To this end, we need estimates
on the quality of the approximation in Vh .

Lemma 1.30. Let Th be a non-degenerate triangulation of Ω. Then for v ∈ H 1 (Ω) it holds

that
∥v − Ih v∥H(div;Ω) ≤ ch|v|H 1 (Ω) + inf ∥div v − qh ∥L2 (Ω) .
qh ∈Wh

26
 1.5 Computing the Solution of the Discrete Problem

Proof. Let τ ∈ Th . The trace theorem 1.26 shows that the functional v 7→ e v · νe ds, e ∈ ∂τ ,
R

is bounded on [H 1 (τ )]2 . Due to [Π0 ]2 ⊂ Vh it holds that Iτ v = v for v ∈ [Π0 ]2 . The

Bramble-Hilbert lemma and the transformation to the reference element show
∥v − Ih v∥L2 (Ω) ≤ ch|v|H 1 (Ω) .
The estimate for ∥div(v − Ih v)∥L2 (Ω) is a consequence of (1.34).

From Poincaré’s inequality we obtain

inf ∥p − qh ∥L2 (Ω) ≤ ch|p|H 1 (Ω) .
qh ∈Wh

Hence, the previous lemma together with Theorem 1.23 using −div u = f yields the error
estimate
 
∥u − uh ∥H(div;Ω) + ∥p − ph ∥L2 (Ω) ≤ c h|u|H 1 (Ω) + h|p|H 1 (Ω) + inf ∥f − fh ∥L2 (Ω) .
fh ∈Wh

If the right-hand side satisfies f ∈ H 1 (Ω) then also the third term leads to a power of h.

1.5 Computing the Solution of the Discrete Problem

The discretization of saddle point problems (1.27) results in linear systems of equations of
the form     
A B x f
T = . (1.36)
B −C y g
Due to the coercivity of the bilinear form a, the matrix A ∈ Rm×m having the entries
aij := a(φj , φi ), i, j = 1, . . . , m,
is (symmetric) positive definite. Here, we assume that
Vh = span{φ1 , . . . , φm } and Wh = span{ψ1 , . . . , ψn }.
From the discrete LBB condition (1.30) it follows that B ∈ Rm×n with
bij = b(φi , ψj ), i = 1, . . . , m, j = 1, . . . , n,
has rank n; see the remark on page 14. Here, we assume that m ≥ n. As a generalization
of the saddle point structure we consider an additional (symmetric) positive semi-definite
matrix C ∈ Rn×n . The coefficient matrix
 
A B
 :=
B T −C
is thus symmetric. However, due to
 T    T  
x x 0 0
 = xT Ax > 0, x ̸= 0,  = −y T Cy ≤ 0,
0 0 y y
it is indefinite.

27
1 Non-Conforming and Mixed Methods

Lemma 1.31. The linear system (1.36) is uniquely solvable and the Schur complement
S := C + B T A−1 B is (symmetric) positive definite.

Proof. The matrix S ∈ Rn×n is obviously positive definite. The decomposition

    
I 0 A B A B
=
−B T A−1 I B T −C 0 −S

shows that  is non-singular.

For indefinite problems we have mentioned the MINRES method in the lecture series Intro-
duction to Iterative Numerical Methods. In the following we introduce two methods which
are tailored to saddle point problems.

1.5.1 The Uzawa Method

The following classical iteration method for the solution of (1.36) in the case C = 0 is called
Uzawa method. Let (x0 , y0 ) ∈ Rm+n . For k = 0, 1, . . . let

xk+1 = A−1 (f − Byk ),

yk+1 = yk + ω(B T xk+1 − g)

with a parameter ω > 0. Eliminating xk+1 , we obtain

yk+1 = yk + ω(B T A−1 (f − Byk ) − g) = yk + ω(B T A−1 f − g − Syk ),

which is the Richardson method applied to

Sy = B T A−1 f − g.

From the lecture series Introduction to Iterative Numerical Methods we obtain

Theorem 1.32. The Uzawa method converges iff ω ∈ (0, 2/λmax (S)). The optimum relax-
ation parameter is
2
ωopt := .
λmin (S) + λmax (S)

In each step of the Uzawa method, the linear system

Axk+1 = f − Byk (1.37)

has to be solved, which is equivalent to minimizing the functional

1
hk (x) := xT Ax − xT (f − Byk ).
2

28
 1.5 Computing the Solution of the Discrete Problem

In order to reduce the numerical effort, xk+1 is not found as the exact solution of (1.37)
but an iterative method (e.g., the cg method) is used. In the Arrow-Hurwicz method
only one step with stepsize ε is taken in the direction of the negative gradient −∇hk (xk ) =
f − Byk − Axk of hk in xk : Let (x0 , y0 ) ∈ Rm+n . For k = 0, 1, . . . let

xk+1 = xk + ε(f − Axk − Byk ),

yk+1 = yk + ω(B T xk+1 − g).

1.5.2 The Bramble-Pasciak CG Method

Although  is indefinite, it is possible to introduce a transformation T such that T  is

positive definite. This allows to exploit the improved convergence properties of the conjugate
gradients (CG) method.

Let CA ∈ Rm×m be a positive definite matrix with

0 < xT CA x < xT Ax, x ̸= 0. (1.38)

Multiplication of (1.36) from the left by

ACA−1 − I 0
 
T := ∈ R(m+n)×(m+n)
B T CA−1 −I

leads to
ACA−1 A − A (ACA−1 − I)B x ACA−1 f − f
      
x
T Â = = . (1.39)
y B T (CA−1 A − I) C + B T CA−1 B y B T CA−1 f − g

Theorem 1.33. The matrix T Â is (symmetric) positive definite.

Proof. T Â is obviously symmetric. Due to the assumptions, A − CA is positive definite.

From
 T    T 
(ACA−1 − I)Ax + (ACA−1 − I)By

u x u
T Â =
v y v B T (CA−1 A − I)x + (C + B T CA−1 B)y
= uT (A − CA )CA−1 (Ax + By) + v T [B T CA−1 (A − CA )x + (C + B T CA−1 B)y]

we obtain in particular
 T  
x x
T Â = [(A − CA )x + By]T CA−1 [(A − CA )x + By] + xT (A − CA )x + y T Cy.
y y

The decomposition      
x x x
= 0 + 1 ,
y 0 y
where x1 is the unique solution of Ax1 = −By and x0 := x − x1 , shows

xT1 Ax1 = y T B T A−1 By.

29
1 Non-Conforming and Mixed Methods

From this we obtain

 T  
x0 x
T Â 0 = xT0 (A − CA )CA−1 (A − CA )x0 + xT0 (A − CA )x0
0 0
≥ λmin ((A − CA )CA−1 (A − CA )) ∥x0 ∥22
| {z }
=:c1 >0

and
 T  
x1 x 1 1
T Â 1 = xT1 CA x1 + xT1 (A − CA )x1 + y T Cy = xT1 Ax1 + y T (C + B T A−1 B)y
y y 2 2
1 1
≥ λmin (A) ∥x1 ∥22 + λmin (S) ∥y∥22 .
|2 {z } |2 {z }
=:c2 >0 =:c3 >0

Due to  T  
x0 x
T Â 1 = xT0 (A − CA )CA−1 (Ax1 + By) = 0
0 y
we have
 T    T    T  
x x x0 x0 x1 x
T Â = T Â + T Â 1 ≥ c1 ∥x0 ∥22 + c2 ∥x1 ∥22 + c3 ∥y∥22 .
y y 0 0 y y

The estimate ∥x∥22 = ∥x0 + x1 ∥22 ≤ 2∥x0 ∥22 + 2∥x1 ∥22 leads to
 T  
x x 1
T Â ≥ min{c1 , c2 }∥x∥22 + c3 ∥y∥22
y y 2
and thus the assertion.

We show that T Â is spectrally equivalent with

 
A − CA
Ĉ := ,
CS

where CS ∈ Rn×n is matrix that is spectrally equivalent to the Schur complement S. In

addition to the assumption (1.38), let CA futhermore be spectrally equivalent to A.

Theorem 1.34. If CA is spectrally equivalent to A, i.e.

αA xT CA x ≤ xT Ax ≤ βA xT CA x for all x ∈ Rm

and if CS is spectrally equivalent to S, i.e.

αS xT CS x ≤ xT Sx ≤ βS xT CS x for all x ∈ Rn

with constants 1 < αA ≤ βA and 0 < αS ≤ βS independent of m and n, then Ĉ is spectrally

equivalent to T Â, i.e.

α̂ xT T Âx ≤ xT Ĉx ≤ β̂ xT T Âx for all x ∈ Rm+n .

30
 1.5 Computing the Solution of the Discrete Problem

Proof. We have to estimate the smallest and the largest eigenvalue of Ĉ −1 T Â. To this end,
consider the eigenvalue problem
ACA−1 A − A (ACA−1 − I)B x1
     
A − CA 0 x1
T −1 T −1 =λ .
B (CA A − I) C + B CA B x2 0 C S x2
The first row reads

(ACA−1 − I)Ax1 + (ACA−1 − I)Bx2 = λ(A − CA )x1 = λ(ACA−1 − I)CA x1

or equivalently Bx2 = (λCA − A)x1 .

If λ ∈ [1, βA ] then nothing has to be proved. Hence, we may consider the case λ ̸∈ [1, βA ].
Then λCA − A is invertible and we obtain

x1 = (λCA − A)−1 Bx2 .

Plugging this into the second equality, we obtain

B T CA−1 (A − CA )(λCA − A)−1 Bx2 + [C + B T CA−1 B]x2 = λCS x2 .

Due to

CA−1 (A − CA )(λCA − A)−1 = CA−1 [A − λCA + (λ − 1)CA ](λCA − A)−1

= (λ − 1)(λCA − A)−1 − CA−1 ,

this is equivalent with

(λ − 1)B T (λCA − A)−1 Bx2 + Cx2 = λCS x2 . (1.40)

For λ > βA we have that λCA − A is positive definite. The spectral equivalence of A and CA
implies
λ − βA T
x1 Ax1 ≤ xT1 (λCA − A)x1 , x1 ∈ Rm ,
βA
and thus
βA
xT1 (λCA − A)−1 x1 ≤ xT A−1 x1 , x1 ∈ Rm .
λ − βA 1
Using the spectral equivalence of S and CS , from (1.40) and βA (λ − 1) ≥ λ − βA we obtain
λ T
x Sx2 ≤ λxT2 CS x2 = xT2 Cx2 + (λ − 1)(Bx2 )T (λCA − A)−1 Bx2
βS 2
βA
≤ xT2 Cx2 + (λ − 1) (Bx2 )T A−1 Bx2
λ − βA
λ−1 T
≤ βA x Sx2
λ − βA 2
and thus
λ λ−1
≤ βA ,
βS λ − βA
i.e. λ2 − βA (βS + 1)λ + βA βS ≤ 0. From this it follows that λ− ≤ λ ≤ λ+ , where
r
1 1
λ± = βA (βS + 1) ± [βA (βS + 1)]2 − βA βS .
2 4

31
1 Non-Conforming and Mixed Methods

Hence, we have
βA < λ ≤ λ+ =: β̂.
It remains to investigate the case λ < 1. Then A − λCA is positive definite. The spectral
equivalence of A and CA yields
βA − λ T
xT1 (A − λCA )x1 ≤ x1 Ax1
βA
and thus
βA
xT1 (A − λCA )−1 x1 ≥ xT1 A−1 x1 , x1 ∈ Rm .
βA − λ
Using the spectral equivalence of S and CS , again from (1.40) and βA /(βA − λ) = 1/(1 −
λ/βA ) ≤ 1/(1 − λ) we obtain
λ T
x2 Sx2 ≥ λxT2 CS x2 = xT2 Cx2 + (1 − λ)(Bx2 )T (A − λCA )−1 Bx2
αS
βA
≥ xT2 Cx2 + (1 − λ) (Bx2 )T A−1 Bx2
βA − λ
1−λ T
≥ βA x Sx2
βA − λ 2
and thus
1−λ λ
βA ≤ ,
βA − λ αS
which is equivalent with λ2 − βA (αS + 1)λ + αS βA ≤ 0, i.e. λ− ≤ λ ≤ λ+ , where
r
1 1
λ± = βA (αS + 1) ± [βA (αS + 1)]2 − αS βA .
2 4
In summary we obtain
1 > λ ≥ λ− =: α̂.

The following method is obtained from applying the cg method w.r.t. the scalar product
(·, ·)Ĉ −1 to (1.39). In the following, all vectors x ∈ Rm+n are split (if necessary) into the
upper components x1 ∈ Rm and the lower components x2 ∈ Rn . Note that due to r(k+1) =
r(k) − αk w(k) and w(k) = T ŵ(k) it holds that
(k+1) (k+1) (k) (k) (k) (k)
v1 = (A−CA )−1 r1 = (A−CA )−1 r1 −αk (A−CA )−1 (A−CA )CA−1 ŵ1 = v1 −αk CA−1 ŵ1 .
(0) (0)
Similarly, with r(0) = T r̂(0) it follows that v1 = CA−1 r̂1 .

Algorithm 1.35 (Bramble-Pasciak cg method).

(0) (0) (0) (0) (0) (0)
r̂1 = f − Ax1 − Bx2 , r̂2 = g − B T x1 + Cx2 ;
(0) (0) (0) (0) (0) (0) (0) (0) (0) (0)
v1 = CA−1 r̂1 , r1 = Av1 − r̂1 , r2 = B T v1 − r̂2 ; v2 = CS−1 r2 ;
p(0) = v (0) and ρ0 = (v (0) , r(0) );
k = 0;
do {

32
 1.5 Computing the Solution of the Discrete Problem

(k) (k) (k) (k) (k) (k)

ŵ1 = Ap1 + Bp2 ; ŵ2 = B T p1 − Cp2 ;
(k) (k) (k) (k) (k)
s(k) = CA−1 ŵ1 ; w1 = As(k) − ŵ1 ; w2 = B T s(k) − ŵ2 ;
x(k+1) = x(k) + αk p(k) , where αk = ρk /(p(k) , w(k) );
r(k+1) = r(k) − αk w(k) ;
(k+1) (k) (k+1) (k+1)
v1 = v1 − αk s(k) ; v2 = CS−1 r2 ;
(k+1) (k+1)
ρk+1 = (v ,r );
p(k+1) = v (k+1) + γk p(k) , where γk = ρk+1 /ρk ;
k = k + 1;
} while (ρk > ε2 ρ0 );

33
2 Stokes’ Equations

To describe the stationary flow of an incompressible and extremely viscous liquid (e.g. honey)
in a bounded domain Ω ⊂ Rd , d = 2, 3, the system of Stokes’ equations is used

⃗ + ∇p = f
−∆u in Ω, (2.1a)
div u = 0 in Ω, (2.1b)
u = u0 auf ∂Ω (2.1c)

with the vector-Laplacian ∆ ⃗ = [∆, . . . , ∆]T . Here, u : Ω → Rd is the unknown velocity

and p : Ω → R denotes the unknown pressure. In contrast to less viscous fluids, no mi-
croscale vortices occur here. The system is linear and therefore relatively easy to analyze.
Nevertheless fundamental numerical difficulties can be observed.

The condition (2.1b) guarantees the preservation of mass. This can be seen from the Gauss
divergence theorem for arbitrary ω ⊂ Ω
Z Z
u · ν ds = div u dx = 0.
∂ω ω

Thus, on average, as much liquid flows out through the boundary ∂ω as flows in. In partic-
ular, the compatibility condition Z
u0 · ν ds = 0
∂Ω

has to hold for the Dirichlet data u0 . Furthermore, it is obvious that the pressure p cannot
be unique because only its gradient appears. In order to guarantee that p is uniquely
determined, we thus need a further scaling condition. A possible choice is
Z
p dx = 0.
Ω

2.1 Variational Formulation

Using the integration by parts formulas

⃗ ϕ)L2 (Ω) = (∇u, ∇ϕ)L2 (Ω)

−(∆u, and (∇p, ϕ)L2 (Ω) = −(p, div ϕ)L2 (Ω) , ϕ ∈ [C0∞ (Ω)]d ,

we obtain from (2.1a)

(∇u, ∇ϕ)L2 (Ω) − (div ϕ, p)L2 (Ω) = (f, ϕ)L2 (Ω) .

34
 2.1 Variational Formulation

As usual, we assume homogeneous Dirichlet boundary conditions u0 = 0 on ∂Ω. With the

bilinear forms a : [H01 (Ω)]d × [H01 (Ω)]d → R,
d Z
X
a(v, w) := (∇v, ∇w)L2 (Ω) = ∂j vi ∂j wi dx,
i,j=1 Ω

and b : [H01 (Ω)]d × L20 (Ω) → R,

b(v, q) := −(div v, q)L2 (Ω) ,

where we set L20 (Ω) := {q ∈ L2 (Ω) : Ω q dx = 0}, the variational formulation of Stokes’
R

equations reads: find u ∈ [H01 (Ω)]d and p ∈ L20 (Ω) such that

a(u, v) + b(v, p) = (f, v)L2 (Ω) , v ∈ [H01 (Ω)]d , (2.2a)

b(u, q) = 0, q ∈ L20 (Ω). (2.2b)

Notice that these are d + 1 equations. For d = 2 they read in componentwise notation
Z Z Z
∂1 u1 ∂1 v1 + ∂2 u1 ∂2 v1 dx − p ∂1 v1 = f1 v1 dx, v1 ∈ H01 (Ω),
ZΩ ZΩ ZΩ
∂1 u2 ∂1 v2 + ∂2 u2 ∂2 v2 dx − p ∂2 v2 = f2 v2 dx, v2 ∈ H01 (Ω),
Ω Z Ω Ω

(∂1 u1 + ∂2 u2 )q dx = 0, q ∈ L20 (Ω).

Ω

Theorem 2.1. Weak solutions (u, p) ∈ [H01 (Ω)]d × L20 (Ω) of (2.2) with u ∈ [C 2 (Ω)]d and
p ∈ C 1 (Ω) are classical solutions (2.1) for all right-hand sides f ∈ [L2 (Ω)]d .

Proof. Let q := div u ∈ L2 (Ω). Then there is a constant c such that q − c ∈ L20 (Ω).
Using (2.2b), we obtain
Z
∥div u∥L2 (Ω) = (div u, q)L2 (Ω) = (div u, c)L2 (Ω) = c div u dx.
2
Ω

The boundary condition u|∂Ω = 0 and the Gauss divergence theorem show
Z Z
div u dx = u · ν ds = 0.
Ω ∂Ω

Hence, ∥div u∥L2 (Ω) = 0 and u ∈ [C 2 (Ω)]d yield that div u(x) = 0 for every x ∈ Ω.

The first equation (2.1a) is obtained by returning to a Poisson problem. The weak solution
u ∈ [H01 (Ω)]d satisfies

(∇u, ∇v)L2 (Ω) = g(v), g(v) := (f, v)L2 (Ω) + (p, div v)L2 (Ω) .

Hence, u is the classical solution of the Poisson problem

−∆u = g = f − ∇p in Ω

with homogeneous Dirichlet boundary conditions, which is the desired equation.

35
2 Stokes’ Equations

Equation (2.2) appears to be a saddle point problem with a bilinear form a coercive on [H01 (Ω)]d .
In order to guarantee the existence of solutions of (2.2), according to Theorem 1.19 only the
LBB condition for the bilinear form b must be proved. This is done via the following estimate,
the proof of which is beyond the scope of this lecture.

Lemma 2.2. Let Ω be a bounded Lipschitz domain. Then

(i) The image of the linear mapping ∇ : L2 (Ω) → [H −1 (Ω)]d is closed.

(ii) There is a constant c = c(Ω) such that

∥p∥L2 (Ω) ≤ c(∥∇p∥H −1 (Ω) + ∥p∥H −1 (Ω) ) for all p ∈ L2 (Ω),

∥p∥L2 (Ω) ≤ c∥∇p∥H −1 (Ω) for all p ∈ L20 (Ω). (2.3)

Theorem 2.3. Let Ω be a bounded Lipschitz domain. Then Stokes’ equations (2.2) have a
unique solution (u, p) ∈ [H01 (Ω)]d × L20 (Ω) for all f ∈ ([H01 (Ω)]d )′ .

Proof. In order to be able to apply Theorem 1.19, we have to check the LBB condition for b.
For p ∈ L20 (Ω) we obtain from (2.3)
1
∥∇p∥H −1 (Ω) ≥ ∥p∥L2 (Ω) .
c
The definition of the H −1 (Ω)-norm shows that there is v ∈ [H01 (Ω)]d with ∥v∥H 1 (Ω) = 1 and
1 1
(v, ∇p)L2 (Ω) ≥ ∥∇p∥H −1 (Ω) ≥ ∥p∥L2 (Ω) .
2 2c
Since b(v, p) = −(div v, p)L2 (Ω) = (v, ∇p)L2 (Ω) , this shows
1
b(v, p) = (v, ∇p)L2 (Ω) ≥ ∥p∥L2 (Ω)
2c
and thus the LBB condition.

2.2 The Discrete LBB Condition for Stokes’ Equations

As we have already noticed, the LBB condition (1.22) does not automatically imply the
discrete LBB condition (1.30) for discrete spaces Vh × Wh . In the case of Stokes’ equations
it is possible that a non-trivial pressure q ∈ Wh exists in the kernel of b, i.e. we have

(div v, q)L2 = 0 for all v ∈ Vh .

In this case, the discretization Vh × Wh is called unstable. A weaker form of instability is

present if no q ∈ Wh \ {0} is in the kernel of b but the LBB constant depends on the grid
size h.

36
 2.2 The Discrete LBB Condition for Stokes’ Equations

Instable Stokes elements

The Q1 /P0 rectangle element is quite popular due to its simplicity. This element discretizes
the flow with piecewise bilinear functions, the pressure is approximated with piecewise con-
stant functions:

Vh := {v ∈ [C(Ω)]2 : v|∂Ω = 0 and v|τ bilinear for all τ ∈ Th },

Wh := {q ∈ L20 (Ω) : q|τ ∈ Π0 for all τ ∈ Th }.

This discretization is unstable because the kernel of Bh′ : Wh → Vh′ , Bh′ q := b(·, q), is non-
trivial. We number the vertices in the element τij as shown in the figure.

(i, j + 1) (i + 1, j + 1)

(i + 12 , j + 12 )

(i, j) (i + 1, j)

Since div v is linear and q is constant on each element, we obtain

Z
q div v dx = h2 q(i+ 1 ,j+ 1 ) div v(i+ 1 ,j+ 1 )
2 2 2 2
τij
1
= h2 q(i+ 1 ,j+ 1 ) v1,(i+1,j+1) + v1,(i+1,j) − v1,(i,j+1) − v1,(i,j) +
2 2 2h


+v2,(i+1,j+1) + v2,(i,j+1) − v2,(i+1,j) − v2,(i,j) .

Summing over the elements and sorting the terms w.r.t. the vertices yields
Z X
q div v dx = h2

v1,(i,j) (∇1 q)ij + v2,(i,j) (∇2 q)ij
Ω i,j

with the finite differences

1  
(∇1 q)ij := q(i+ 1 ,j+ 1 ) + q(i+ 1 ,j− 1 ) − q(i− 1 ,j+ 1 ) − q(i− 1 ,j− 1 ) ,
2h 2 2 2 2 2 2 2 2

1  
(∇2 q)ij := q(i+ 1 ,j+ 1 ) + q(i− 1 ,j+ 1 ) − q(i+ 1 ,j− 1 ) − q(i− 1 ,j− 1 ) .
2h 2 2 2 2 2 2 2 2

Since v ∈ [H01 (Ω)]2 , the summation takes into account only inner vertices. We have q ∈
Ker Bh′ if Z
q div v dx = 0 for all v ∈ Vh ,
Ω

i.e. (∇1 q)ij and (∇2 q)ij vanish at all inner vertices. This happens if

q(i+ 1 ,j+ 1 ) = q(i− 1 ,j− 1 ) , q(i+ 1 ,j− 1 ) = q(i− 1 ,j+ 1 ) .

2 2 2 2 2 2 2 2

37
2 Stokes’ Equations

To satisfy these two equations, it is not required that q is constant. It is sufficient that
(
a, i + j even,
q(i+ 1 ,j+ 1 ) =
2 2
b, i + j odd,

with values a, b such that Ω q dx = 0. In particular, a and b have opposite signs, i.e.
R

a checkerboard pattern appears. Hence, this kind of instability is called checkerboard

instability.

Remark. A conceivable remedy for the problem would be the restriction of Vh to the or-
thogonal complement of Yh = Ker Bh′ . In this case the kernel is trivial, but the constant in
the discrete LBB condition (1.29) depends on h.

The Taylor-Hood Element

The commonly used Taylor-Hood element employs polynomials for the velocity of a degree
higher than for the pressure. Compared to the unstable Q1 /P0 elements, the pressure is
assumed to be continuous:

Vh := {v ∈ [C(Ω)]d : v|∂Ω = 0 and v|τ ∈ [Π2 ]d for all τ ∈ Th },

Wh := {q ∈ C(Ω) ∩ L20 (Ω) : q|τ ∈ Π1 for all τ ∈ Th }.

In the following theorem, we investigate discretizations Th of Ω ⊂ R2 using triangles.

Theorem 2.4. Let Th be a quasi-uniform triangulation of the polygonal domain Ω ⊂ R2 .

Then the Taylor-Hood element satisfies the discrete LBB condition (1.30).

Proof. We show that for all qh ∈ Wh there is 0 ̸= vh ∈ Vh such that

b(vh , qh )
≥ c̃LBB ∥qh ∥L2 (Ω) .
∥vh ∥H 1 (Ω)

38
 2.2 The Discrete LBB Condition for Stokes’ Equations

The proof will be done in two steps.

(i) For given K > 0 we introduce the set

WhK := {qh ∈ Wh : ∥qh ∥L2 (Ω) > Kh|qh |H 1 (Ω) }.

Let qh ∈ WhK . Due to the continuous LBB condition (1.22) there is v ∈ [H01 (Ω)]2 such that
∥v∥H 1 (Ω) = 1 and
b(v, qh ) ≥ cLBB ∥qh ∥L2 (Ω) .
Furthermore, let 0 ̸= vh ∈ Vh the Clément approximation of v, which satisfies

∥v − vh ∥L2 (Ω) ≤ ch∥v∥H 1 (Ω) ,

where the constant does not depend on qh and thus not on K. In summary it follows that

b(vh , qh ) = b(v, qh ) + b(vh − v, qh ) ≥ cLBB ∥qh ∥L2 (Ω) + (vh − v, ∇qh )L2 (Ω)
≥ cLBB ∥qh ∥L2 (Ω) − ∥v − vh ∥L2 (Ω) |qh |H 1 (Ω)
≥ cLBB ∥qh ∥L2 (Ω) − ch∥v∥H 1 (Ω) |qh |H 1 (Ω) .

Due to the boundedness of the Clément operator we have ∥vh ∥H 1 (Ω) ≤ c′ ∥v∥H 1 (Ω) = c′ .
Furthermore, qh ∈ WhK implies

b(vh , qh ) cLBB − c/K

≥ ∥qh ∥L2 (Ω) .
∥vh ∥H 1 (Ω) c′

Hence, the LBB condition is valid for sufficiently large K.

(ii) We consider the case qh ̸∈ WhK . For every edge e in the interior of Ω denote by te its
tangential vector and by νe its normal. At the midpoint me of e we set

te · vh (me ) = te · (∇qh )(me ), νe · vh (me ) = 0.

Furthermore, let vh (x) := 0 in the vertices of all triangles and in the midpoints of edges
on ∂Ω. This construction shows

∥vh ∥L2 (Ω) ≤ |qh |H 1 (Ω) . (2.4)

Since (∇qh )|τ is constant and (vh · ∇qh )|τ is quadratic, it will be proved in the exercises1 that
3 3
|τ | X |τ | X
Z
vh · ∇qh dx = vh (mei ) · (∇qh )|τ = [tei · (∇qh )|τ ]2
τ 3 i=1 3 i=1
|τ |
Z
c
≥ c∥(∇qh )|τ ∥22 = ∥(∇qh )|τ ∥22 dx
3 3 τ
with a constant c > 0 depending on the degeneracy of the grid Th . Summation over all
τ ∈ Th yields Z
c X c
b(vh , qh ) ≥ ∥(∇qh )|τ ∥22 dx = |qh |2H 1 (Ω) .
3 τ ∈T τ 3
h

1
Let T ∈ Rm×d be a full-rank matrix and T = QR be a QR-decomposition of T . Then ∥T T x∥22 =
xT QRRT QT x ≥ λmin (RRT )∥QT x∥22 = λmin (RRT )∥x∥22 .

39
2 Stokes’ Equations

Then we obtain
2
b(vh , qh ) c |qh |H 1 (Ω) c |qh |H 1 (Ω) ∥qh ∥L2 (Ω) c
≥ ≥ ≥ ∥qh ∥L2 (Ω) .
∥vh ∥H 1 (Ω) 3 ∥vh ∥H 1 (Ω) 3Kh ∥vh ∥H 1 (Ω) 3Kcinv

Here, we have used (2.4) and the inverse inequality ∥vh ∥H 1 (Ω) ≤ cinv h−1 ∥vh ∥L2 (Ω) .

Theorem 2.5. The Taylor-Hood element satisfies the error estimate

∥u − uh ∥H 1 (Ω) + ∥p − ph ∥L2 (Ω) ≤ chk ∥u∥H k+1 (Ω) + ∥p∥H k (Ω)



with k ∈ {1, 2} provided that u ∈ H k+1 (Ω) and p ∈ H k (Ω).

Proof. Due to the approximation estimates

inf ∥u − vh ∥H 1 (Ω) ≤ chk ∥u∥H k+1 (Ω) , inf ∥p − qh ∥L2 (Ω) ≤ chk ∥p∥H k (Ω)
vh ∈Vh qh ∈Wh

for k ∈ {1, 2}, the assertion follows from Theorem 1.23.

Remark. Instead of using quadratric ansatz functions for the velocity, also piecewise linear
functions on a refined triangulation can be employed:

Vh := {v ∈ [C(Ω)]d : v|∂Ω = 0 and v|τ ∈ [Π1 ]d for all τ ∈ Th/2 },

Wh := {q ∈ C(Ω) ∩ L20 (Ω) : q|τ ∈ Π1 for all τ ∈ Th }.

This choice also yields a stable discretization, and the number of degrees of freedom is the
same as for Taylor-Hood element. This variant is also referred to as Taylor-Hood discretiza-
tion.

The MINI element

A disadvantage of the Taylor-Hood element is that the values for pressure and velocity refer
to different grids or different polynomial bases. The following MINI element is a P1 /P1
element with a velocity space that is enriched with bubble functions

bτ (x) := 9λτ,1 λτ,2 λτ,3 ,

i.e.

Vh := {v ∈ [C(Ω)]2 : v|∂Ω = 0 and v|τ ∈ [Π1 ⊕ span{bτ }]2 for all τ ∈ Th },

= [S01,0 (Th ) ⊕ B3 ]2

with B3 := {v ∈ C(Ω) : v|τ ∈ span{bτ } for all τ ∈ Th } and as before

Wh := {q ∈ C(Ω) ∩ L20 (Ω) : q|τ ∈ Π1 for all τ ∈ Th }.

40
 2.2 The Discrete LBB Condition for Stokes’ Equations

Figure 2.1: Bubble functions for the unit triangle.

Figure 2.2: MINI element (v is prescribed at the dots, q at the crosses).

Here, (λτ,1 , λτ,2 , λτ,3 ) denote the barycentric coordinates of x in the triangle τ ∈ Th . Let
x1 , x2 , x3 be the vertices of τ . If for x ∈ τ it holds that

(λ1 + λ2 + λ3 ) x = λ1 x1 + λ2 x2 + λ3 x3

and at least one of the coefficients λ1 , λ2 , λ3 ∈ R+0 does not vanish then we say that the
coefficients (λ1 , λ2 , λ3 ) are barycentric coordinates of x with respect to τ . The vertices
themselves have the coordinates (1, 0, 0), (0, 1, 0) and (0, 0, 1). Notice that
P3with λi also αλi
are barycentric coordinates of x for any α > 0. The additional condition i=1 λi = 1 makes
λi unique. On the unit triangle these coincide with x1 , x2 and 1 − x1 − x2 . Notice that
bτ |∂τ = 0.

Theorem 2.6. Let Ω be a polygonal domain in R2 . Then the MINI element satisfies the
discrete LBB condition.

Proof. We prove the validity of the Fortin criterion, i.e., we show the existence of a projection
Πh : [H01 (Ω)]2 → Vh satisfying

(div(v − Πh v), qh )L2 (Ω) = 0 for all v ∈ [H01 (Ω)]2 and qh ∈ Wh

(0)
and its boundedness with an h-independent constant cΠ > 0. To this end, let Πh : H01 (Ω) →
S01,0 (Th ) be the Clément operator. This choice satisfies
(0)
∥v − Πh v∥L2 (Ω) ≤ c1 h∥v∥H 1 (Ω)

41
2 Stokes’ Equations

(0) (1)
and ∥Πh v∥H 1 (Ω) ≤ c2 ∥v∥H 1 (Ω) . Furthermore, we define Πh : L2 (Ω) → B3 as
X Z −1 Z
(1)
Πh v = βτ bτ with βτ := bτ dx v dx.
τ ∈Th τ τ

(1)
The projection Πh can be understood as a two-step process. As first, an L2 -projection is
done onto the piecewise constant functions. Then the constants are replaced with bubble
functions with the same integral. This implies
(1)
∥Πh v∥L2 (Ω) ≤ c3 ∥v∥L2 (Ω) .
(0) (1) (0) (1)
We set Πh v := Πh v + Πh (v − Πh v). Since (Πh v)|τ preserves the mean value of v|τ , we
(0)
obtain with ṽ := v − Πh v
Z Z
(1)
v − Πh v dx = ṽ − Πh ṽ dx = 0 for all τ ∈ Th .
τ τ

Integration by parts yields

X
(div(v − Πh v), qh )L2 (Ω) = (div(v − Πh v), qh )L2 (τ )
τ ∈Th
X Z
=− (v − Πh v, ∇qh )L2 (τ ) + (v − Πh v) · ν qh ds.
τ ∈Th ∂τ

The integrals over inner edges vanish due to changing signs of the normals. The integrals
along the edges on the boundary vanish as v|∂Ω = 0 = (Πh v)|∂Ω . Since ∇qh is piecewise
constant, we obtain in summary
X Z
b(v − Πh v, qh ) = (div(v − Πh v), qh )L2 (Ω) = − (∇qh )|τ v − Πh v dx = 0.
τ ∈Th τ

We extend the defintion of Πh to vector-valued functions v ∈ [H01 (Ω)]2 by applying it com-

(0) (1)
ponentwise. Using the inverse inequality, it follows from the continuity of Πh and Πh
that
(0) (1) (0)
∥Πh v∥H 1 (Ω) ≤ ∥Πh v∥H 1 (Ω) + cinv h−1 ∥Πh (v − Πh v)∥L2 (Ω)
(0)
≤ c2 ∥v∥H 1 (Ω) + c3 cinv h−1 ∥v − Πh v∥L2 (Ω)
≤ (c2 + c1 c3 cinv ) ∥v∥H 1 (Ω) .

Lemma 1.24 gives the assertion.

Remark. From the previous proof it is visible that only one additional degree of freedom
is sufficient to guarantee the desired stability.

Theorem 2.7. The MINI element satisfies the error estimate


∥u − uh ∥H 1 (Ω) + ∥p − ph ∥L2 (Ω) ≤ ch ∥u∥H 2 (Ω) + ∥p∥H 1 (Ω) ,

provided that u ∈ H 2 (Ω) and p ∈ H 1 (Ω).

42
 2.3 Solution of the Discrete Problems

Proof. From the approximation estimates

inf ∥u − vh ∥H 1 (Ω) ≤ ch∥u∥H 2 (Ω) , inf ∥p − qh ∥L2 (Ω) ≤ ch∥p∥H 1 (Ω)

vh ∈Vh qh ∈Wh

the assertion follows from Theorem 1.23.

The MINI element has a smaller number of degrees of freedom than the Taylor-Hood element.
However, the latter can benefit from a possible higher regularity of the solution.

2.3 Solution of the Discrete Problems

The variational formulation of the Stokes’ equations as a saddle point problem leads to
systems of linear equations of the form (1.36). For the solution of these we have introduced
the Uzawa and the Bramble-Pasciak cg method. The Uzawa method was equivalent to the
Richardson iteration for the Schur complement, and for the Bramble-Pasciak cg method
(given a preconditioner for A) the condition of the Schur complement is decisive for the
convergence.

The following lemma shows that the condition number of the Schur complement S =
B T A−1 B in the case of Stokes’ equations is bounded. To this end, let Vh be an m-dimensional
and Wh an n-dimensional space for the discretization of V and W , respectively. Denote by
Jh : Rn → Wh the natural injection
n
X
Jh x := xi φi
i=1

with the finite element basis {φi } of Wh . We already know that the mass matrix M = Jh∗ Jh ∈
Rn×n for Wh is well-conditioned. The following result for M −1 S avoids the discussion about
the scaling of the entries of S w.r.t. h and holds in d-dimensional space.

Lemma 2.8. Let the discrete LBB condition (1.30) be valid. Then the maximum and
minimum eigenvalue of M −1 S satisfy

0 < c̃2LBB ≤ λmin ≤ λmax ≤ d.

Proof. The smallest eigenvalue of M −1 S can be bounded from below2 by

∥M −1 Sy∥2M (Sy, M −1 Sy) (Sy, M −1 Sy) (Sy, y)

λ2min ≥ min = min = min
y̸=0 ∥y∥2M y̸=0 ∥y∥2M y̸=0 (Sy, y) ∥y∥2M
(y, M −1 Sy)S (Sy, y) (Sy, y)
≥ min 2
min 2
= λmin min .
y̸=0 ∥y∥S y̸=0 ∥y∥M y̸=0 ∥y∥2 M

2
choose y the eigenvalue corresponding to λmin

43
2 Stokes’ Equations

With x := A−1 By the nominator satisfies due to ∥x∥A = maxz̸=0 (x, z)A /∥z∥A

(Sy, y) = (B T A−1 By, y) = (A−1 By, By) = (x, Ax) = ∥x∥2A

(Ax, z)2 (By, z)2
= max = max .
z̸=0 ∥z∥2A z̸=0 ∥z∥2A

In summary we obtain
2
(By, z)2 (y, B T z)

λmin ≥ min max = min max .
y̸=0 z̸=0 ∥z∥2 2
A ∥y∥M y̸=0 z̸=0 ∥z∥A ∥y∥M

For vh = Jh z and qh = Jh y it holds that

∥z∥A = ∥∇vh ∥L2 (Ω) , ∥y∥M = ∥qh ∥L2 (Ω) , and (y, B T z) = (div vh , qh )L2 (Ω) .

Hence, using the discrete LBB condition (1.30) we obtain

 2
(div vh , qh )L2 (Ω)
λmin ≥ min max ≥ c̃2LBB .
qh ∈Wh vh ∈Vh ∥∇vh ∥L2 (Ω) ∥qh ∥L2 (Ω)

Similarly, we obtain
2 2
∥div vh ∥L2 (Ω)
 
(div vh , qh )L2 (Ω)
λmax ≤ max max ≤ max .
qh ∈Wh vh ∈Vh ∥∇vh ∥L2 (Ω) ∥qh ∥L2 (Ω) vh ∈Vh ∥∇vh ∥L2 (Ω)

The assertion follows from the Cauchy-Schwarz inequality as

Z d
!2 Z d
X X
∥div vh ∥2L2 (Ω) = ∂i vh dx ≤ d (∂i vh )2 dx ≤ d∥∇vh ∥2L2 (Ω) .
Ω i=1 Ω i=1

44
3 Finite Element Methods for
Electro-Magnetic Problems

The classical electromagnetic field is described by the four vectors E, D, H, and B which are
functions of the position x ∈ R3 and the time t ∈ R. The vectors E and H are referred to as
the electric and magnetic field, while D and B are the electric and magnetic displacements,
respectively. They are governed by Maxwell’s equations

div B = 0, (Gauss’s law for magnetism)

∂t B + curl E = 0, (Faraday’s law of induction)
div D = ρ, (Gauss’s law)
curl H − ∂t D = J. (Ampère’s circuital law)

Here and in the following,

 
∂2 v3 − ∂3 v2
curl v := ∂3 v1 − ∂1 v3  = ∇ × v
∂1 v2 − ∂2 v1

is the curl operator of v = [v1 , v2 , v3 ]T . ρ denotes the given electric charge density and
J the electric current density. In the presence of conductors Ω moving with velocity v, we
adopt Ohm’s law
J = σ(E + v × B) in R3 ,
where σ is the electric conductivity of the given conductor occupying Ω, and σ = 0 outside
the conductor.

Remark. From the identity

0 = div curl H = div J + ∂t div D = div J + ∂t ρ, (3.1)

it follows that div J = −∂t ρ, which is the charge conservation law.

The relation between displacement field D and electric field E, as well as the magnetic
field H and the magnetic displacement B is specified by the material. For linear materials
D = εE and µH = B hold, where ε denotes the permittivity and µ is the permeability of
the material. In the following we shall assume that

0 < εm ≤ ε ≤ εM , 0 < µm ≤ µ ≤ µ M , (3.2)

with εm , εM , µm , µM ∈ R.

45
3 Finite Element Methods for Electro-Magnetic Problems

Wave propagation

The Maxwell’s equations are of hyperbolic type. To see this, let us first consider the simple
case of a non-moving homogeneous isotropic medium, i.e. the case where ε, µ, and σ are
constant and v = 0. Then
µ ∂t H + curl E = 0,
curl H − ε ∂t E = J.
Taking the time derivative of the first equation and the curl of the second one, we obtain by
adding
εµ ∂t2 H + curl curl H = curl J.
Since div H = µ−1 div B = 0, we obtain from the vector identity −∆ ⃗ = curl curl − ∇div
⃗ = curl J.
εµ ∂t2 H − ∆H
If J is given, we obtain a hyperbolic problem that describes propagation of electromagnetic
waves in the space with speed (εµ)−1/2 . If Ohm’s law is assumed in Ω, we obtain
⃗ = 0 in Ω,
εµ ∂t2 H + µσ ∂t H − ∆H
⃗ = 0 in R3 \ Ω
εµ ∂t2 H − ∆H
with a condition at infinity and suitable interface conditions. Here, we also have wave
propagation but waves are damped in the conductor, the damping being proportional to the
electric conductivity σ.

Time-harmonic Maxwell’s equations

The time-harmonic formulation of Maxwell’s equations is based on the following assumptions

E(x, t) = Re (e−iωt E(x)),
D(x, t) = Re (e−iωt D(x)),
H(x, t) = Re (e−iωt H(x)),
B(x, t) = Re (e−iωt B(x)).
These assumptions make sense, for instance, when studying the propagation of electromag-
netic waves at a given frequency ω ̸= 0. With this assumption the Maxwell’s equations
read
div(µH) = 0,
−iωµH + curl E = 0,
div(εE) = r,
curl H + iωεE = j,
where ρ(x, t) = Re (e−iωt r(x)) and J(x, t) = Re (e−iωt j(x)). Eliminating H, the previous
system can be reformulated as a second order equation
curl(µ−1 curl E) − ω 2 εE = F,

46
and the divergence condition
−ω 2 div(εE) = div F,
where F = iωj.

If the electric field is to be found in a domain Ω, a perfect conducting surface implies the
boundary condition E × ν = 0 on ∂Ω.

Remark. Let ε, µ be constant. Then, eliminating E instead of H, one obtains for the
magnetic field
curl curl H − ω 2 εµH = curl j.
Using div H = 0, we obtain the Helmholtz equation
⃗ − ω 2 εµH = curl j.
−∆H

Time-independent problems

In the static case, we have ∂t B = 0 = ∂t D. This situation enables decoupling electricity

and magnetism in the following way. Due to curl E = 0 in R3 , we deduce the existence of a
scalar potential ϕ : R3 → C.

Theorem 3.1. Let Ω be a simply connected Lipschitz domain. E ∈ [L2 (Ω)]3 satisfies
curl E = 0 iff there is ϕ ∈ H 1 (Ω)/R such that E = −∇ϕ.

Assuming Ohm’s law in the static conductor Ω, i.e. J = σE in Ω, from (3.1) we obtain the
electrostatic problem
−div σ∇ϕ = 0 in Ω.
This is an elliptic equation that requires suitable boundary conditions. If a part of the
boundary ∂Ω can be connected to an electricity generator, we prescribe the potential ϕ
when we have a voltage generator (Dirichlet condition) or the normal derivative of ϕ by

σ ∂ν ϕ = J · ν

when we have a current generator (Neumann condition).

In the static case, the displacement current term ∂t D can be neglected and Maxwell’s equa-
tions for the magnetic field simplifies to curl H = J. The first equation div B = 0 implies
the existence of a vector potential A.

Theorem 3.2. Let Ω be a simply connected Lipschitz domain. B ∈ [L2 (Ω)]3 satisfies
Z
div B = 0 in Ω, B · ν = 0,
∂Ω

iff there is A ∈ [H 1 (Ω)]3 such that B = curl A. A is unique if we impose the gauge condition
div A = 0.

47
3 Finite Element Methods for Electro-Magnetic Problems

Then the so-called magnetostatic problem reads

curl(µ−1 curl A) = J.

Due to ν · B = ν · (∇ × A) = ∇ · (A × ν), the boundary condition A × ν = 0 means that the

normal component of B vanishes.

Summary

The electrostatic problem is a Poisson-type problem, which has already been discussed. We
focus on the problems involving the curl curl operator. Both, the magnetostatic and the
time-harmonic problem (formulation with E and div F = 0) lead to the following unified
type
1
curl curl u − ω 2 u = f in Ω, (3.3a)
µ
div u = 0 in Ω, (3.3b)
u × ν = 0 on ∂Ω (3.3c)

with given f such that div f = 0.

Remark. Notice that by taking the divergence of the first equation we obtain that div u = 0,
so that the second equation is redundant provided ω ̸= 0.

It is clear that the solvability of the previous system is strictly related to the frequency ω.
It is ill-posed for ω 2 = λ, where λ solves
1
curl curl u = λu in Ω, (3.4a)
µ
div u = 0 in Ω, (3.4b)
u × ν = 0 on ∂Ω (3.4c)

for some non-vanishing vector field u. The previous problem is the so-called cavity prob-
lem for Maxwell’s equations, which is also known as interior Maxwell’s eigenvalue
problem.

3.1 The Space H(curl; Ω)

For Maxwell’s equations the space

H(curl; Ω) := {v ∈ [L2 (Ω)]3 : curl v ∈ [L2 (Ω)]3 }

is of great importance. It is equipped with the norm

∥v∥2H(curl;Ω) := ∥v∥2L2 (Ω) + ∥curl v∥2L2 (Ω) .

48
 3.1 The Space H(curl; Ω)

Theorem 3.3. Let Ω ⊂ R3 be a bounded Lipschitz domain with unit outer normal ν. Then
the mapping γt : [C ∞ (Ω)]3 → [C ∞ (∂Ω)]3 with γt v = v|∂Ω ×ν can be extended to a continuous
linear operator γt : H(curl; Ω) → [H −1/2 (∂Ω)]3 and the Green’s identity for v ∈ H(curl; Ω)
and w ∈ [H 1 (Ω)]3
Z Z
⟨γt v, γ0 w⟩∂Ω = v · curl w dx − w · curl v dx. (3.5)
Ω Ω

Note that γt : H(curl; Ω) → [H −1/2 (∂Ω)]3 is not surjective since γt v has vanishing normal
component.

With the trace γt we can define H0 (curl; Ω) to be the subspace of H(curl; Ω) consisting of
vector fields v with vanishing tangential component v × ν on the boundary of Ω.

Using the terminology of differential forms, the space H(curl; Ω) can be identified with an
element of the following de Rham complex

∇ curl div
H 1 (Ω) −→ H(curl; Ω) −→ H(div; Ω) −→ L2 (Ω) (3.6)

The de Rham complex shows the sequence of spaces of differential forms An under the
exterior derivative dn : An → An+1 . In particular, we have dn+1 dn = 0. If we consider
homogeneous boundary conditions then the de Rham complex takes the following form

∇ curl div
H01 (Ω) −→ H0 (curl; Ω) −→ H0 (div; Ω) −→ L2 (Ω) (3.7)

Here, H0 (div; Ω) denotes the subspace of H(div; Ω) consisting of functions v with vanishing
normal component v · ν = 0 on the boundary ∂Ω.

If Ω is simply connected then the two sequences are exact, i.e., the range of each operator
equals the kernel of the next one. For instance, every divergence-free vector field of H(div; Ω)
is the curl of an element of H(curl; Ω) and a curl-free vector field of H0 (curl; Ω) is the gradient
of an element of H01 (Ω); see Thms. 3.1 and 3.2.

The previous sequences are related with an important tool of vector calculus: the so-called
Helmholtz decomposition. The following theorem is a result of

Theorem 3.4. Let Ω be a simply connected Lipschitz domain. Then v ∈ [L2 (Ω)]3 can
be decomposed into the sum of an irrotational (curl-free) and a solenoidal (divergence-free)
part, which refer to a unique scalar potential ϕ ∈ H 1 (Ω)/R and a unique vector potential
A ∈ [H 1 (Ω)]3 with div A = 0 such that

v = ∇ϕ + curl A.

The Helmholtz decomposition is orthogonal in the sense that ∇ϕ · curl A = 0.

49
3 Finite Element Methods for Electro-Magnetic Problems

Proof. Let ϕ ∈ H01 (Ω) by defined by

(∇ϕ, ∇w)L2 (Ω) = (v, ∇w)L2 (Ω) for all w ∈ H01 (Ω).
Then v0 := v − ∇ϕ is divergence-free and ∂Ω v0 · ν ds = 0 follows from the Gauss divergence
R

theorem. Theorem 3.2 shows the existence of a vector potential A such that v0 = curl A
with div A = 0. The orthogonality of the decomposition follows from
∇ϕ · (∇ × A) = −A · (∇ × ∇ϕ) = 0.

In order to deal with physically more relevant multiply connected domains Ω, a natural
approach consists in assuming the existence of a finite number of mutually disjoint sets
Λ1 , . . . , Λn such that Ω \ ni=1 Λi is simply connected. In this case, the sequences (3.6) and
S
(3.7) are no longer exact. However, the range of each operator is a closed subset of the kernel
of the next one. A standard procedure is to consider quotient spaces, which are then called
cohomology spaces. For instance, looking at the kernel of the curl operator H0 (curl0 ; Ω)
thought as a subset of H0 (curl; Ω), we can consider the orthogonal complement K of ∇H01 (Ω)
in H0 (curl0 ; Ω), i.e.
H0 (curl0 ; Ω) = ∇H01 (Ω) ⊕ K.
While on simply connected domains K = {0}, it turns out that on multiply connected
domains the dimension of K equals n. The space K consists of vector fields v which satisfy
curl v = 0, div v = 0 (via integration by parts) in Ω and v × ν = 0 on ∂Ω. It is clear that the
presence of a non-trivial cohomology space K has the effect of changing the construction of
Helmholtz’s decomposition which then reads
v = ∇ϕ + curl A + ζ
with ζ ∈ K.

The following Friedrichs inequality for the space H0 (curl; Ω) is a central tool for the
coercivity of Maxwell problems. It plays the same role as the Friedrichs inequality for the
space H01 (Ω). Let
Z = H0 (curl; Ω) ∩ H(div0 ; Ω)
consisting of vector fields in H(curl; Ω) with vanishing tangential component along the
boundary and with zero divergence in Ω. Then there is α > 0 such that
∥curl v∥L2 (Ω) ≥ α∥v∥L2 (Ω) for all v ∈ Z. (3.8)

3.2 Variational Formulation

The variational formulation of the problem (3.3) with ω ̸= 0 can be naturally obtained in
the space H0 (curl; Ω). After multiplying the first equation by a test function v ∈ H0 (curl; Ω)
and integrating by parts using (3.5) we obtain
Z Z Z
1
curl u · curl v dx − ω 2
u · v dx = f · v dx
Ω µ Ω Ω
| {z } | {z }
a(u,v) ℓ(v)

50
 3.2 Variational Formulation

with given f ∈ H(div0 ; Ω). Hence, the variational problem reads: find u ∈ H0 (curl; Ω) such
that
a(u, v) = ℓ(v) for all v ∈ H0 (curl; Ω). (3.9)
Notice that we did not take into account the condition div u = 0 due to the remark after (3.3).
For small ω the condition div u = 0 will be a significant source of trouble in the numerical
approximation. Therefore, we enforce the divergence-free condition by requiring that u ⊥
∇H01 (Ω). Introducing the bilinear form

b(q, v) = (∇q, v)L2 (Ω) ,

the corresponding mixed formulation reads: find u ∈ H0 (curl; Ω) and p ∈ H01 (Ω) such that

a(u, v) + b(p, v) = ℓ(v) for all v ∈ H0 (curl; Ω), (3.10a)

b(q, u) = 0 for all q ∈ H01 (Ω). (3.10b)

The variational problem (3.9) and its mixed formulation (3.10) are equivalent in the case
ω ̸= 0. This is due to ∇H01 (Ω) ⊂ H0 (curl; Ω). Choosing v = ∇q with q ∈ H01 (Ω) in (3.9)
shows that any solution u of (3.9) is a solution of the mixed problem (3.10) with p = 0. Vice
versa, we may choose v = ∇p in the mixed problem. As a consequence, we have a(u, v) = 0
and ℓ(v) = 0, which shows that 0 = b(p, v) = (∇p, v)L2 (Ω) = ∥∇p∥2L2 (Ω) . Together with the
boundary conditions, this implies p = 0. Notice that this holds also for ω = 0. Hence, (3.10)
is more general than (3.9).

The following theorem gives the conditions for the well-posedness of problem (3.10).

Theorem 3.5. Assume that ω 2 is not an interior Maxwell eigenvalue (see problem (3.4)).
Then (3.10) has a unique solution u ∈ H0 (curl; Ω), which satisfies

∥u∥H(curl;Ω) ≤ C∥f ∥L2 (Ω) ,

2 1+µM λi
where C = c max{µ−1m + ω , |λi −ω 2 | , i = 1, 2, . . . } and λi , i = 1, 2, . . . , denote the interior
Maxwell eigenvalues and µm , µM were defined in (3.2).

Proof. It is clear that the bilinear forms a : H0 (curl; Ω) × H0 (curl; Ω) → R and b : H01 (Ω) ×
H0 (curl; Ω) → R are continuous. The continuity constant of a is µ−1 m +ω .
2

The inf-sup condition for b is trivial: given q ∈ H01 (Ω), we can take v = ∇q ∈ H0 (curl; Ω)
and obtain from ∥v∥H(curl;Ω) = ∥∇q∥L2 (Ω) ≤ ∥q∥H 1 (Ω)

b(q, v) ∥∇q∥2L2 (Ω)

sup ≥ ≥ c∥q∥H 1 (Ω) ,
0̸=v∈H0 (curl;Ω) ∥v∥H(curl;Ω) ∥q∥H 1 (Ω)

where c is related to Friedrichs’ inequality.

The coercivity of the bilinear form a is less immediate. It is clear that if ω = 0, the bilinear
form a is coercive on the kernel

Z = {v ∈ H0 (curl; Ω) : b(q, v) = 0 for all q ∈ H01 (Ω)}

51
3 Finite Element Methods for Electro-Magnetic Problems

consisting of divergence-free vector fields in H0 (curl; Ω), which is a consequence of Friedrichs’

inequality (3.8). In view of Lemma 1.13, for ω ̸= 0 we need to prove that
a(u, v)
sup ≥ α∥u∥H(curl;Ω) for all u ∈ Z. (3.11)
0̸=v∈Z ∥v∥H(curl;Ω)
Observe that any element of Z can be represented as a Fourier series in terms of the interior
Maxwell eigenfunctions {ui }, i.e.
(µ−1 curl ui , curl v)L2 (Ω) = λi (ui , v)L2 (Ω) for all v ∈ H0 (curl; Ω),
(∇q, ui )L2 (Ω) = 0 for all q ∈ H01 (Ω),
(ui , uj )L2 (Ω) = 0, i ̸= j, ∥ui ∥H(curl;Ω) = 1,
with eigenvalues λi ≥ 0. As a consequence of the orthogonality, we have
(µ−1 curl ui , curl uj )L2 (Ω) = λi (ui , uj )L2 (Ω) = 0, i ̸= j.
We also observe that eigenfunctions ui satisfy
1
(1 + µM λi ) (ui , ui )L2 (Ω) = (ui , ui )L2 (Ω) + µM ( curl ui , curl ui )L2 (Ω)
µ
≥ (ui , ui )L2 (Ω) + (curl ui , curl ui )L2 (Ω) = ∥ui ∥2H(curl;Ω) = 1.
Given u ∈ Z and its representation u = ∞
P∞
i=1 αi ui , let us construct v = i=1 βi ui ∈ Z in
P
order to prove (3.11). Define βi = ±αi such that
αi βi (λi − ω 2 ) = αi2 |λi − ω 2 | ≥ 0
for every i. In particular, it turns out that ∥v∥2H(curl;Ω) = i βi2 = i αi2 = ∥u∥2H(curl;Ω) .
P P
With this v it follows
∞ ∞ ∞ ∞
1X X X X
a(u, v) = ( αi curl ui , βj curl uj )L2 (Ω) − ω (
2
αi ui , βj uj )L2 (Ω)
µ i=1 j=1 i=1 j=1
∞ ∞
X 1 X
= αi βi ( curl ui , curl ui )L2 (Ω) − ω 2 αi βi (ui , ui )L2 (Ω)
i=1
µ i=1
∞
X
= αi βi (λi − ω 2 ) (ui , ui )L2 (Ω)
i=1
∞
X |λi − ω 2 |
≥ αi2
i=1
1 + µ M λi
|λi − ω 2 | 2 |λi − ω 2 |
≥ min ∥u∥H(curl;Ω) = min ∥u∥H(curl;Ω) ∥v∥H(curl;Ω) .
i 1 + µM λi i 1 + µM λi

The assertion follows from Theorem 1.19.

An alternative approach

Problem (3.3) can be reformulated as the following problem: find u ∈ H0 (curl; Ω)∩H(div; Ω)
such that
(µ−1 curl u, curl v)L2 (Ω) + γ(div u, div v)L2 (Ω) + α(u, v)L2 (Ω) = (f, v)L2 (Ω) (3.12)

52
 3.3 Approximation of the Variational Problem

for all v ∈ H0 (curl) ∩ H(div; Ω). Here, γ > 0 and α ∈ R denote constants. When div f = 0,
any solution u of (3.12) belongs to the space H(div0 ; Ω) and is also a solution of the following
curl-curl problem: find u ∈ H0 (curl; Ω) such that
(µ−1 curl u, curl v)L2 (Ω) + α(u, v)L2 (Ω) = (f, v)L2 (Ω)
for all v ∈ H0 (curl; Ω). This can be seen from using the test function v = ∇ϕ, where
ϕ ∈ H01 (Ω) is the solution of γ∆ϕ + αϕ = ψ with given ψ ∈ L2 (Ω). Then (3.12) becomes
γ(div u, ∆ϕ)L2 (Ω) + α(div u, ϕ)L2 (Ω) = 0,
which shows that (div u, ψ)L2 (Ω) = 0 and thus div u = 0.

For positive α, (3.12) is uniquely solvable by the Riesz representation theorem applied to
the Hilbert space
XN := H0 (curl; Ω) ∩ H(div; Ω)
with the inner product
(v, w)XN := (curl v, curl w)L2 (Ω) + (div v, div w)L2 (Ω) + (v, w)L2 (Ω) .
For α ≤ 0, the problem (3.12) is well-posed as long as α ̸= −λγ,j for j ≥ 1, where λγ,j denote
the eigenvalues of the modified problem. In particular, in the case where α = 0 and Ω is
simply connected, it is uniquely solvable due to Friedrichs’ inequality
∥v∥L2 (Ω) ≤ c(∥curl v∥L2 (Ω) + ∥div v∥L2 (Ω) ), v ∈ XN .

Since it is difficult to construct finite-element subspaces for XN , the problem (3.12) used to be
discretized by H 1 -conforming nodal finite elements. However, the subspace [H 1 (Ω)]3 ∩ XN
turns out to be a closed subspace of XN . Therefore, any H 1 -conforming finite element
method for (3.12) must fail if the solution u does not belong to [H 1 (Ω)]3 , which happens
when Ω is non-convex.

Example 3.6. Consider the weak solution u ∈ H 1 (Ω) of

−∆u = f in Ω,
u = 0 on ∂Ω
with f ∈ L2 (Ω). Define w = ∇u ∈ [L2 (Ω)]3 . Then curl w = 0 and div w = ∆u = −f ∈
L2 (Ω). Hence, w ∈ H(curl; Ω) ∩ H(div; Ω). But w ̸∈ [H 1 (Ω)]3 unless u ∈ H 2 (Ω).

Hence, the solutions obtained by H 1 -conforming finite element methods in such situations
converge to the wrong solution (the projection of u in [H 1 (Ω)]3 ∩ XN ).

3.3 Approximation of the Variational Problem

In the following, the finite-element approximation of the variational problem (3.9) is done
by considering a finite-dimensional conforming subspace Vh ⊂ H0 (curl; Ω) and by looking for
uh ∈ Vh such that
a(uh , vh ) = ℓ(vh ) for all vh ∈ Vh . (3.13)

53
3 Finite Element Methods for Electro-Magnetic Problems

The discussion of the previous section shows that a discrete counterpart of the mixed for-
mulation (3.10) will play an important role for the analysis of (3.13). Let Qh ⊂ H01 (Ω) be a
conforming finite-dimensional approximation of H01 (Ω). We can consider the discretization
of (3.10): find (uh , ph ) ∈ Vh × Qh such that

a(uh , vh ) + b(ph , vh ) = ℓ(vh ) for all vh ∈ Vh , (3.14a)

b(qh , uh ) = 0 for all qh ∈ Qh . (3.14b)

Remark. The finite-element approximation (3.13) is equivalent with (3.14) if ∇Qh ⊂ Vh .

This can be proved similar to the discussion before Theorem 3.5. The required inclusion
∇Qh ⊂ Vh can be guaranteed by choosing Qh to be made of standard Lagrangian elements
and Vh as edge finite elements, which will be defined in Sect. 3.5. Although (3.13) and (3.14)
are equivalent in exact arithmetics, the problem (3.14) is more stable in practice.

In order to prove the discrete version of Theorem 3.5, we have to consider the discrete interior
Maxwell eigenvalues λi,h , i = 1, . . . , Nh , satisfying

uh ∈ Vh : (µ−1 curl uh , curl vh )L2 (Ω) = λh (uh , vh )L2 (Ω) for all vh ∈ Vh .

Theorem 3.7. Assume that ∇Qh ⊂ Vh and that ω 2 is not an interior Maxwell eigenvalue
and let u be the unique solution of (3.9). Let us assume that h is such that ω 2 does
not coincide with any discrete interior Maxwell’s eigenvalues λi,h , i = 1, . . . , Nh . Then,
problem (3.14) has a unique solution uh which satisfies

∥u − uh ∥H(curl;Ω) ≤ C inf ∥u − vh ∥H(curl;Ω) ,

vh ∈Vh

1+µM λi,h
where C = c max{µ−1 2
m +ω , |λi,h −ω 2 |
, i = 1, 2, . . . , Nh }.

Proof. The proof is analogous to the one of Theorem 3.5. In particular, we shall show the
well-posedness of the mixed problem (3.14). This is done via proving the discrete inf-sup
conditions for the bilinear forms a and b.

The inf-sup condition for b is immediate due to ∇Qh ⊂ Vh . The inf-sup condition for a
can be proved as in the continuous case by considering a basis of Vh consisting of discrete
Maxwell eigenfunctions ui,h . The discrete Friedrichs inequality in this case reads

λi,h
(µ−1 curl ui,h , curl ui,h )L2 (Ω) ≥ ∥ui,h ∥2H(curl;Ω) .
1 + µM λi,h

Arguing as in the proof of Theorem 3.5, we get the inf-sup constant for a

|λi,h − ω 2 |
min .
i 1 + µM λi,h

The final estimate is then obtained from Theorem 1.23 by observing that we do not need to
estimate the terms ∥p − ph ∥H 1 (Ω) and inf qh ∈Qh ∥p − qh ∥H 1 (Ω) since both p and ph are zero.

54
 3.4 Conforming Finite Elements for H(curl; Ω)

3.4 Conforming Finite Elements for H(curl; Ω)

Partitioning Ω ⊂ R3 into sub-domains is an essential feature for discrete methods. Continuity

properties at interfaces between such sub-domains are a crucial part in the definition of finite-
element approximations.

Let Th be a partition of Ω into elements τ (tetrahedra or hexahedra). We only deal with

compatible meshes in the sense that the intersection between two elements is a common face,
edge or vertex.

Let
W (Ω) = {v ∈ [L2 (Ω)]3 : v|τ ∈ H(curl; τ ) for all τ ∈ Th }
the set of functions which are in H(curl; τ ) for each element τ . The following lemma shows
necessary and sufficient conditions for a function in W (Ω) to be in H0 (curl; Ω). Notice first
that due to (3.5) for v ∈ H(curl; Ω)
X XZ Z Z Z
⟨γt v, γ0 ϕ⟩∂τ = v · curl ϕ dx − ϕ · curl v dx = v · curl ϕ dx − ϕ · curl v dx
τ ∈Th τ ∈Th τ τ Ω Ω

= ⟨γt v, γ0 ϕ⟩∂Ω , ϕ ∈ [H 1 (Ω)]3 . (3.15)

P
Lemma 3.8. It holds that H0 (curl; Ω) = {v ∈ W (Ω) : τ ∈Th ⟨γt v, γ0 ϕ⟩∂τ = 0 for all ϕ ∈
[H 1 (Ω)]3 }.

Proof. If v ∈ H0 (curl; Ω), then by (3.15) we have that τ ∈Th ⟨γt v, ϕ⟩∂τ = 0 for all ϕ ∈
P
[H (Ω)] . Let us prove the converse. Using integration by parts (3.5) on each element τ , we
1 3

obtain
Z XZ XZ
v · curl ϕ dx = v · curl ϕ dx = ϕ · curl v dx, ϕ ∈ [H 1 (Ω)]3 .
Ω τ ∈Th τ τ ∈Th τ

Hence, ℓ(ϕ) := Ω v · curl ϕ dx ≤ ∥ϕ∥L2 (Ω) τ ∈Th |v|2H(curl;τ ) for all ϕ ∈ [C0∞ (Ω)]3 . The Riesz
R P

representation theorem shows that there is w ∈ [L2 (Ω)]3 such that

Z
ℓ(ϕ) = w · ϕ dx, ϕ ∈ [C0∞ (Ω)]3 .
Ω

Hence, w = curl v is the weak derivative, which shows that v ∈ H(curl; Ω). With (3.15) we
have v ∈ H0 (curl; Ω).

Hence, functions in H(curl; Ω) have “continuous” tangential traces at the interfaces.

Construction of Elements

The following general definition of finite elements is due to Ciarlet.

55
3 Finite Element Methods for Electro-Magnetic Problems

Definition 3.9. A finite element is a tripel (τ, Π, Σ) satisfying the following properties:

(i) τ ⊂ Rd is a polyhedron.

(ii) Π ⊂ C(τ ) is a subspace of finite dimension s. The functions in Π are called form
functions.

(iii) Σ is a set of s linearly independent functionals defined over Π. Every p ∈ Π is uniquely

determined by the values of the s functionals in Σ. These are the so-called generalized
interpolation conditions.

Let τ ⊂ R3 be a tetrahedron and let ei , i = 1, . . . , 6, denote its edges and fi , i = 1, . . . , 4, its

faces. Given k ∈ N0 we define the Nédélec elements

Nk (τ ) = [Πk (τ )]3 ⊕ {x × p(x), p ∈ [Πk,h (τ )]3 },

where Πk,h denotes the space of homogeneous polynomials of degree k, i.e. polynomials p
satisfying p(αx) = αk p(x) for all x ∈ R3 and α ∈ R.

Lemma 3.10. For k ≥ 0 and any v ∈ Nk (τ ) the degrees of freedom can be defined by
Z
v · t pk ds, pk ∈ Πk (ei ), i = 1, . . . , 6, (3.16a)
ei
Z
v · pk−1 dσ, pk−1 ∈ [Πk−1 (fj )]3 with pk−1 · νfj = 0, j = 1, . . . , 4, (3.16b)
f
Zj
v · pk−2 dx, pk−2 ∈ [Πk−2 (τ )]3 , (3.16c)
τ

where t denotes the tangential direction. The space Πk denotes the space of polynomials of
degree (the sum of the univariate degrees) at most k.

Let us check that the number of conditions
involved in the previous
lemma is the same
as the dimension of Nk . As dim Πk = k+d k
, (3.16a) imposes k+1
k
= k + 1 conditions
on each edge, (3.16b) gives 2 k−1 = k(k + 1) conditions on each face, and (3.16c) adds
k+1

3 k+1 = (k − 1)k(k + 1)/2 more conditions in each of the three components of the vector


k−2
field. As it can be easily seen, the sum of the three contributions is equal to

dim Nk = (k + 1)(k + 3)(k + 4)/2.

Lemma 3.8 states that in order to construct an approximation of H(curl; Ω) we need to

ensure continuity of the tangential components across neighboring elements. The following
lemma implies that the degrees of freedom (3.16a) and (3.16b) guarantee such continuity.
Therefore elements in Nk are called edge elements.

56
 3.4 Conforming Finite Elements for H(curl; Ω)

Lemma 3.11. Assume that v ∈ Nk (τ ) is such that (3.16a) and (3.16b) vanish on a given
face f ⊂ τ and on all edges e of f . Then v × ν = 0 on f .

Our last comment is about the characterization of the kernel

Nk0 (τ ) := {v ∈ Nk (τ ) : curl v = 0}

of the curl operator in Nk (τ ). It is easily seen that Nk0 (τ ) = ∇Πk+1 (τ ). This will be of great
importance to satisfy the assumption ∇Qh ⊂ Vh made in Sect. 3.3. Let us define

Nk+ (τ ) = Nk (τ )/∇Πk+1 (τ ).

Obviously, curl Nk (τ ) = curl Nk+ (τ ).

The following lemma establishes a connection with Raviart-Thomas spaces (see Sect. 1.4.1),
which are used for a conforming discretization of H(div; Ω). With

RTk (τ ) := {p ∈ [Πk (τ )]3 + xΠk (τ ) : p · ν ∈ Rk (∂τ )},

where Rk (∂τ ) := {v ∈ L2 (∂τ ) : v|e ∈ Πk (e) for all e ⊂ ∂τ }, and

RT0k (τ ) := {v ∈ RTk (τ ) : div v = 0}

we have

Lemma 3.12. The curl operator is bijective from Nk+ (τ ) onto RT0k (τ ).

Proof. Applying curl to Nk (τ ) and using that curl(x × p) = x · div p − 3p ∈ [Πk (τ )]3 + xΠk (τ ),
we obtain that curl Nk+ (τ ) ⊂ RT0k (τ ). A simple count shows that the dimension of Nk+ (τ )
and the dimension of RT0k (τ ) are equal.

Interpolation operators for H(curl; Ω)

The first step to derive error estimates is the definition of an interpolation operator. The
main question for the definition of the interpolation operator is the regularity of the function
to be interpolated. Remember that H 1 -regularity is not sufficient for a nodal interpolant.

In order to define the degrees of freedom in (3.16a) and (3.16b), we assume that v belongs
to the space

X(τ ) = {v ∈ [Lp (τ )]3 , curl v ∈ [Lp (τ )]3 , ν × v|∂τ ∈ [Lp (∂τ )]3 }

for some p > 2.

57
3 Finite Element Methods for Electro-Magnetic Problems

Theorem 3.13. Let τ be an affine element and Iτ : X(τ ) → Nk (τ ) be the interpolation

operator. Then there is a constant c depending only on k and on the shape of τ such that

∥ curl(v − Iτ v)∥L2 (τ ) ≤ chm

τ |curl v|H m (τ ) , v ∈ [H m (τ )]3 ,

for 1 ≤ m ≤ k + 1.

The spaces Nk defined in the previous section can be used to define conforming approxima-
tions of H(curl; Ω). With the element spaces Nk (τ ) we define

Nk (Th ) := {v ∈ H(curl; Ω) : v|τ ∈ Nk (τ ), τ ∈ Th }.

Given ε > 0, let us consider the space

X := H(curl; Ω) ∩ [H 1/2+ε (Ω)]3 ,

so that the global interpolation operator Ih : X → Nk (Th ), Ih v|τ := Iτ v|τ , can be defined.

Theorem 3.14. Let Th be a non-degenerate family of decompositions of Ω. Then there is

a constant c independent of h such that for 1 < m ≤ k + 1

∥v − Ih v∥L2 (Ω) ≤ chm |v|H m (Ω) , v ∈ [H m (τ )]3 .

Moreover, for 1/2 < s = 1/2 + ε ≤ 1 and p > 2

∥v − Ih v∥L2 (Ω) ≤ chs (|v|H s (Ω) + ∥curl v∥Lp (Ω) ), v ∈ X.

Finally,
∥ curl(v − Ih v)∥L2 (Ω) ≤ cht |curl v|H t (Ω) , v ∈ [H t (τ )]3 ,
where 1 ≤ t ≤ k + 1.

3.5 Assembling Finite Element Stiffness Matrices

We discuss the construction of two families N of lowest order elements in three space dimen-
sions. Let τ be a tetrahedron and let τ consist of four vertices xi , xj , xk , xl with associated
barycentric coordinates λi , λj , λk , λl ; see Sect. 2.2. Let eij = xj − xi be the edge connecting
xi and xj . The face fijk is defined by the vertices xi , xj , xk . Its outward normal to fijk is
denoted by νijk . The height of a tetrahedron from xl to the opposite face fijk will be denoted
by hl . We note that the gradient of the barycentric coordinates is related to the normals,
i.e.
1
∇λl (x) = − νijk , x ∈ τ.
hl

58
 3.5 Assembling Finite Element Stiffness Matrices

Lowest order Nédélec elements N

Example 3.15 (Lowest order edge elements). We investigate the easiest case k = 0,
which is probably most often used in practice. The corresponding elements are referred to as
Whitney elements. The only meaningful degrees of freedom are those presented in (3.16a).
The space N0 is given by
N0 (τ ) = [Π0 (τ )]3 ⊕ (x × [Π0 (τ )]3 )
           
 1 0 0 0 −x3 x2 
= span 0 , 1 , 0 ,  x3  ,  0  , −x1  .
0 0 1 −x2 x1 0
 

The six degrees of freedom are ei v · t ds, i = 1, . . . , 6, and it can be checked that the
R

quantity v · t for v ∈ N0 (τ ) is constant along the edges.

The form function φij associated with edge eij and the corresponding linear functions lij are
φij = λi ∇λj − λj ∇λi ,
Z
1
lij (v) = v · t ds ≈ [v(xi ) + v(xj )] · eij ,
eij 2
where the approximation is exact when v · t is linear. Due to
Z
φij (xi ) · eij = ∇λj · eij = ∇λj · t ds = λj (xj ) − λj (xi ) = 1,
eij
Z
φij (xj ) · eij = −∇λi · eij = − ∇λi · t ds = −λi (xj ) + λi (xi ) = 1,
eij

we have lij (φij ) = 1. Let emn be another (closed) edge such that emn ∩ eij = ∅. Then
λi |emn = 0 = λj |emn . Hence, without loss of generality we may consider m = i and n ̸∈ {i, j}.
Then ∇λj · emn = 0 and λj |emn = 0 and therefore φij |emn · emn = 0. This shows that the
vector field φij is orthogonal to other edges and verifies lmn (φij ) = 0 for (i, j) ̸= (m, n).

The lowest order edge element is

Nτ = span{φij , eij ⊂ ∂τ } = N0 (τ ),
which is a linear polynomial space. It can be shown that Nτ = {vh ∈ [Π1 (τ )]3 : vh =
a + b × x, a, b ∈ R3 }. In particular, we have ∇λl ∈ Nτ .

For the ease of computing the stiffness matrix, we also calculate the curl of φij
2 ek ′ l ′
curl φij = 2∇λi × ∇λj = sin α,
hi hj ∥ek′ l′ ∥
where ek′ l′ is the conjugate edge, i.e. the edge connecting the other two vertices xk′ and xl′
in τ , and α denotes the angle between νik′ l′ and νjk′ l′ .

The mass matrix can be computed using the integral formula of barycentric coordinates
Z
α1 !α2 !α3 !α4 !
λα1 1 λα2 2 λα3 3 λα4 4 dx = 6 |τ |.
(3 + 4i=1 αi )!
P
τ

59
4 Generalized Finite Elements

We consider Dirichlet boundary value problems

−div(α∇u) = f in Ω, (4.1a)
u = 0 on ∂Ω, (4.1b)

with given f and α. The function α : Ω → R could be nonsmooth, i.e. fail to have continuous
derivatives, but is assumed to satisfy

0 < α0 ≤ α(x) ≤ α1 < ∞, x ∈ Ω.

The variational formulation, find u ∈ E0 (Ω) such that a(u, v) = f (v) for all v ∈ E0 (Ω), of
this problem employs the bilinear form
Z
a(u, v) = α ∇u · ∇v dx,
Ω

which induces the energy norm ∥v∥E(Ω) := a1/2 (v, v) over

E(Ω) = {v : ∥v∥E(Ω) < ∞} and E0 (Ω) := {v ∈ E(Ω) : v|∂Ω = 0}.

The FEM uses a polyhedrization of the domain to construct piecewise polynomial approx-
imating functions. The supports of the shape functions are unions of polyhedra relative to
the polyhedrization. But for domains in 3d, it is quite difficult to generate a good mesh.
One of the important aspects of the Generalized Finite Element Method (GFEM) is that it
permits the use of partition of unity functions, whose support may not depend on any mesh
or may depend on a simple mesh that does not conform to the geometry of Ω. In this sense,
the GFEM is a meshless method and this feature allows to avoid the use of sophisticated
mesh generators. Another important aspect of GFEM is that local approximation spaces
can have functions other than polynomials which locally approximate the unknown solution
well.

4.1 Partition of Unity

Let ωi ⊂ Ω, i = 1, . . . , N , be open sets such that Ω = N

i=1 ωi . We assume that every x ∈ Ω
S
belongs to at most κ of the subdomains ωi . Then let {ϕi }i=1,...,N be a family of functions
defined on Ω such that ϕi satisfies the following properties:

(i) ϕi (x) = 0, x ∈ Ω \ ωi , i = 1, . . . , N ,

(ii)
PN
i=1 ϕi (x) = 1, x ∈ Ω,

60
 4.1 Partition of Unity

(iii) ∥ϕi ∥L∞ (Ω) ≤ c1 , i = 1, . . . , N ,

(iv) ∥∇ϕi ∥L∞ (Ω) ≤ c2

diam ωi
, i = 1, . . . , N ,

where 0 < c1 , c2 < ∞. The property (ii) states that {ϕi } is a partition of unity of Ω.

Example 4.1. As an example, consider the classical FEM with nodal points pi and trian-
gular elements satisfying the minimal angle condition. Let ωi be the patch associated with
node pi , i.e., the union of triangles with pi as one of their vertices. Let ϕi be the piecewise
linear functions with (
1, i = j,
ϕi (pj ) =
0, i ̸= j.
Then it is easily seen that {ϕi } satisfies the previous properties with c1 = 1 and c2 depending
on the minimal angle condition.

Example 4.2. Let Ω = (0, 1)2 and let pij = (ih, jh), i, j = 0, . . . , m, h := 1/m. Set
ωij = BRh (pij ) ∩ Ω,
where R is chosen such that {ωij } covers Ω. With a function ϕ(r), 0 ≤ r < ∞, having
bounded first derivative and ϕ(r) > 0 for 0 ≤ r < R and ϕ(r) = 0 for r ≥ R, define
s 
 2  2
x − ih y − jh 
ϕ̃ij (x, y) = ϕ  + .
h h

The family {ϕ̃ij } satisfies (i), (iii), and (iv) but in general not (ii). If we define the Shepard
functions
ϕ̃ij (x, y)
ϕij (x, y) = Pm ,
k,l=0 ϕ̃kl (x, y)

then the family {ϕij } satisfies all conditions (i)–(iv). To see (iii) and (iv) for this family, we
use the fact that m
X
ϕ̃ij (x, y) ≥ τ > 0, (x, y) ∈ Ω.
i,j=0

To every ωi we associate an mi -dimensional space of functions defined on ωi :

( mi
)
X (i) (i) (i) (i) (i)
Vi = ξi = αj ξj , αj ∈ R, with ξj ∈ E(ωi ) and ξj = 0 on ω i ∩ ∂Ω .
j=1

The space Vi is called a local approximation space. Note that the (essential) Dirichlet
boundary condition is built into Vi . Then we let
( N
)
X
SGFEM = ψ = ϕi ξi , ξi ∈ Vi
i=1
= span{ηij , j = 1, . . . , mi , i = 1, . . . , N },
(i)
where ηij := ϕi ξj are the shape functions for SGFEM . The space SGFEM is called the global
approximation space.

61
4 Generalized Finite Elements

Theorem 4.3. We have SGFEM ⊂ E0 (Ω).

(i) (i)
Proof. Using (i) we see that (ϕi ξj )(x) = 0 for x ∈ ∂ωi ∩ Ω. Hence, ϕi ξj can be extended
(i)
by zero to all of Ω and will be in E(Ω). Furthermore, since ξj = 0 on ω i ∩ ∂Ω, we see that
(i) (i)
ϕi ξj |∂Ω = 0. So, for all i and j, ϕi ξj ∈ E0 (Ω), and hence the span of these functions is
in E0 (Ω).

The Generalized Finite Element Method (GFEM) is now defined to be the Galerkin
method with the conforming discretization S = SGFEM . We denote the approximate solution
by uS = uGFEM ∈ SGFEM . Galerkin orthogonality implies Céa’s lemma

∥u − uS ∥E(Ω) = inf ∥u − ξ∥E(Ω) . (4.2)

ξ∈S

Hence, in order to prove error estimates for ∥u − uS ∥E(Ω) , we have to make sure that ξ u ∈ S
exists such that ∥u − ξ u ∥E(Ω) is small.

For each i we assume that the exact solution u ∈ E(Ω) can be approximated on ωi by a
function ξiu ∈ Vi . The corresponding local errors are
Z 1/2
u
δi := ∥u − ξi ∥L2α (ωi ) = α|u − ξi | dx
u 2
(4.3)
ωi

and Z 1/2
εi := ∥u − ξiu ∥E(ωi ) = α|∇(u − ξiu )|2 dx . (4.4)
ωi
The local approximation is ensured by the appropriate selection of the spaces Vi . We define
the global approximation
XN
u
ξ := ϕi ξiu ∈ SGFEM .
i=1

Theorem 4.4. Suppose u ∈ E0 (Ω). Then

N
!1/2
X
∥u − ξ u ∥L2α (Ω) ≤ c1 κ1/2 δi2
i=1

and !1/2
N N
u 1/2
X δi2 X
∥u − ξ ∥E(Ω) ≤ (2κ) c22 2
+ c 2
1 ε2i .
i=1
(diam ωi ) i=1

Proof. Recalling the definition of ξi and using the fact that {ϕi } is a partition of unity on Ω,
we have
Z Z N
X
u 2
∥u − ξ ∥L2α (Ω) = α|u − ξ | dx =
u 2
α| ϕi (u − ξiu )|2 dx.
Ω Ω i=1

62
 4.1 Partition of Unity

Using the fact that any x ∈ Ω is in at most κ subdomains ωi , we see that the sum N
P
i=1 ϕi (u−
ξiu ) has at most κ terms for any x ∈ Ω. Hence, using the Cauchy-Schwarz inequality, we
have
N
X N
X
u 2
| ϕi (u − ξi )| ≤ κ |ϕi (u − ξiu )|2 .
i=1 i=1
Thus, we have
Z N
X N Z
X N
X
∥u − ξ u ∥2L2α (Ω) ≤κ α |ϕi (u − ξiu )|2 dx ≤ c21 κ α|u − ξiu |2 dx = c21 κ δi2 . (4.5)
Ω i=1 i=1 ωi i=1

The second part follows similarly. Proceeding as above, we have

Z
u 2
∥u − ξ ∥E(Ω) = α|∇(u − ξ u )|2 dx
Ω
Z N
X
= α|∇ ϕi (u − ξiu )|2 dx
Ω i=1
Z N
X
= α| (u − ξiu )∇ϕi + ϕi ∇(u − ξiu )|2 dx
Ω i=1
Z N
!2 Z N
!2
X X
≤2 α (u − ξiu )∇ϕi dx + 2 α ϕi ∇(u − ξiu ) dx
Ω i=1 Ω i=1
N Z
X N Z
X
≤ 2κ α|(u − ξiu )∇ϕi |2 dx + 2κ α|ϕi ∇(u − ξiu )|2 dx.
i=1 ωi i=1 ωi

Hence, using the assumptions (4.3) and (4.4) we obtain

N N
!
X δi2 X
∥u − ξ u ∥2E(Ω) ≤ 2κ c22 + c 2
1 ε2i .
i=1
(diam ωi )2 i=1

Since εi is usually proportional to δi /diam ωi , the terms in the previous estimate are in some
sense balanced. The following theorem gives conditions to ensure this balance.

Theorem 4.5. Let u ∈ E0 (Ω). Suppose, the patches ωi and the local approximation
spaces Vi satisfy the following assumptions:

(i) For all i for which |ω i ∩ ∂Ω| = 0, Vi contains the constant functions and

∥v∥L2α (ωi ) ≤ c3 diam ωi ∥v∥E(ωi )

αv dx = 0.
R
for all v ∈ E(ωi ) satisfying ωi

(ii) For all i for which |ω i ∩ ∂Ω| > 0

∥v∥L2α (ωi ) ≤ c4 diam ωi ∥v∥E(ωi )

for all v ∈ E(ωi ) satisfying v|ωi ∩∂Ω = 0.

63
4 Generalized Finite Elements

Then there is ξ˜iu ∈ Vi so that the corresponding global approximation

N
X
ξ˜u = ϕi ξ˜iu
i=1

satisfies !1/2
N
X
∥u − ξ˜u ∥L2α (Ω) ≤ c5 (diam ωi )2 ε2i ,
i=1
√
where c5 := κ c1 (c23 + c24 )1/2 and

N
!1/2
X
∥u − ξ˜u ∥E(Ω) ≤ c6 ε2i , (4.6)
i=1

where c6 := (2κ(c21 + c22 (c23 + c24 )))1/2 .

Proof. Let ξiu satisfy (4.3) and (4.4). Set Aint = {i ∈ {1, . . . , N } : |ω i ∩ ∂Ω| = 0} and
Abnd = {i ∈ {1, . . . , N } : |ω i ∩ ∂Ω| > 0}.

For i ∈ Aint let ξ˜iu = ξiu + ri , where ri is a constant chosen so that u − ξ˜iu has a zero
α-average on ωi . By assumption (i), ξ˜iu ∈ Vi . Then with v = u − ξ˜iu and noting that
∇(u − ξ˜iu ) = ∇(u − ξiu ), from (i) and (4.4) we have

Z
∥u − ξ˜iu ∥2L2α (ωi ) ≤ c23 (diam ωi )2 α|∇(u − ξ˜iu )|2 dx
Zωi
= c23 (diam ωi )2 α|∇(u − ξiu )|2 dx
ωi
2 2 2
≤ c3 (diam ωi ) εi .

We also have
Z
∥u − ξ˜iu ∥2E(ωi ) = α|∇(u − ξiu )|2 dx = ε2i .
ωi

For i ∈ Abnd let ξ˜iu = ξiu . Now u|ωi ∩∂Ω = 0, and we know that ξ˜iu |ωi ∩∂Ω = 0. Thus, using
assumption (ii) with v = u − ξ˜iu and (4.4), we have

∥u − ξ˜iu ∥L2α (ωi ) = ∥u − ξiu ∥L2α (ωi ) ≤ c4 diam ωi ∥u − ξiu ∥E(ωi ) = c4 diam ωi εi .

Also, ∥u − ξ˜iu ∥2E(ωi ) ≤ ε2i .

64
 4.1 Partition of Unity

Following the steps leading to (4.5), we get

N
X
∥u − ξ˜u ∥2L2α (Ω) ≤ c21 κ ∥u − ξ˜iu ∥2L2α (ωi )
i=1
!
X X
= c21 κ ∥u − ξ˜iu ∥2L2α (ωi ) + ∥u − ξ˜iu ∥2L2α (ωi )
i∈Aint i∈Abnd
N
X
≤ c21 κ(c23 + c24 ) (diam ωi )2 ε2i .
i=1

Similarly, we obtain
N
X
∥u − ξ˜u ∥2E(Ω) ≤ 2κ(c21 + c22 (c23 + c24 )) ε2i .
i=1

Remark. In the previous theorems the conditions on the patch ωi are considerably weak.
The ωi can, in particular, be multiply connected.

We return to the GFEM. Suppose the assumptions of the previous theorem is satisfied and
suppose that u is the exact solution and uGFEM its Galerkin approximation. It follows
from (4.2) and (4.6) that

N
!1/2
X
∥u − uGFEM ∥E(Ω) = inf ∥u − ξ∥E(Ω) ≤ ∥u − ξ˜u ∥E(Ω) ≤ c ε2i .
ξ∈SGFEM
i=1

We can write uGFEM as

X mi
X
N
uGFEM = i=1 cij ηij ,
j=1

where c = [cij ] is the solution of the linear system

mk
N X
X
a(ηkl , ηij )ckl = f (ηij ), 1 ≤ j ≤ mi , 1 ≤ i ≤ N,
k=1 l=1

or Ac = f , where A is the stiffness matrix, whose elements are

Z
Aij,kl = a(ηkl , ηij ) = α ∇ηkl · ∇ηij dx
ωi ∩ωk

and f is the load vector, whose components are

Z
f (ηij ) = f ηij dx.
ωi

65
4 Generalized Finite Elements

4.2 GFEM and Classical FEM

The GFEM is a very general method. It is based on the generalization of the ideas of classical
FEMs and covers many standard FEMs. Using polynomial functions together with other
special functions we get the XFEM, which is a special case of the GFEM. The specific choice
of ϕi and Vi lead to the methods referred to in the literature by different names.

Example 4.6. The classical FEM in 1d based on continuous, piecewise polynomials of

degree k, is the same as a suitably chosen GFEM. We show this for k = 2 by proving that
the finite-dimensional approximation space used in the GFEM is the same as the classical
FEM space.

Suppose Ω = I = (0, 1) and for fixed n let xi = ih, 0 ≤ i ≤ n, with h = 1/n, be uniformly
distributed nodes in I. We consider the subivision of I by the intervals Ii = (xi , xi+1 ). The
standard FEM space is given by

SFEM = {v ∈ C[0, 1] : v|Ii ∈ Πk , i = 0, . . . , n − 1}.

We construct a GFEM space as follows: To each node xi , we associate a function ϕi , which

is the usual piecewise linear continuous hat function centered at xi such that ϕi (xj ) = δij .
We let ω i = supp ϕi = [xi−1 , xi+1 ], 0 < i < n. For i = 0, n we define ω 0 = supp ϕ0 = [x0 , x1 ]
and ω n = supp ϕn = [xn−1 , xn ]. Clearly, {ϕi }ni=0 form a partition of unity in I and satisfy
the assumptions (i)–(iv). For the local approximation spaces Vi , 0 ≤ i ≤ n, we consider

Vi = span{1, x − xi }.

We then define the GFEM space as

n
X
SGFEM = {ψ = ϕi li },
i=0

where li (x) = αi + βi (x − xi ) ∈ Vi with αi , βi ∈ R. The functions li are only defined in ωi ,

but since ϕi (x) = 0 at x = xi−1 and x = xi+1 , ϕi li has a continuous zero-extension to I.

We will show that SFEM = SGFEM . Since the function ϕi li are continuous on I, it is clear
that functions in SGFEM are also continuous on I. Also, since ϕi and li are piecewise linear, it
is clear that every ψ ∈ SGFEM is a continuous, piecewise quadratic function. Thus SGFEM ⊂
SFEM .

We next show that SFEM ⊂ SGFEM , i.e., for a given q ∈ SFEM we can find constants αi , βi
and hence li for 0 ≤ i ≤ n such that
n
X
q(x) = ϕi (x)li (x), x ∈ I. (4.7)
i=0

We first note that equality of q and ϕi li at the nodes xj , 0 ≤ j ≤ n, implies

Pn
i=0

n
X
q(xj ) = ϕi (xj )li (xj ) = ϕj (xj )lj (xj ) = αj . (4.8)
i=0

66
 4.2 GFEM and Classical FEM

We now consider the function ni=0 ϕi li with these αi . Since q and ni=0 ϕi li are both con-
P P
tinuous, have the same values at the nodes xj and their restrictions to the Ii are quadratics,
they will be equal for all x ∈ I if they are equal at the mid points of Ii , i.e.
n
X
q(xj+1/2 ) = ϕi (xj+1/2 )li (xj+1/2 ), 0 ≤ j < n,
i=0

where xj+1/2 = xj + h/2. Imposing these conditions yields

n
X
q(xj+1/2 ) = ϕi (xj+1/2 )li (xj+1/2 )
i=0
= ϕj (xj+1/2 )lj (xj+1/2 ) + ϕj+1 (xj+1/2 )lj+1 (xj+1/2 )
1
= [αj + βj (xj+1/2 − xj ) + αj+1 + βj+1 (xj+1/2 − xj+1 )]
2 
1 h
= αj + αj+1 + (βj − βj+1 ) ,
2 2
which can be written as
2
βj − βj+1 = [2q(xj+1/2 ) − (αj + αj+1 )], 0 ≤ j < n. (4.9)
h
For an arbitrary given value of β0 , we can solve for βi , 1 ≤ i ≤ n, uniquely in terms of β0 .
Using these βi and the αi as given in (4.8), we have constructed li , 0 ≤ i ≤ n, such that (4.7)
holds. Thus SFEM ⊂ SGFEM and using the fact that SGFEM ⊂ SFEM , we have SFEM = SGFEM .

It is well known that for k = 2, a basis of SFEM consists of nodal hat functions ϕi , 0 ≤ i ≤ N ,
and quadratic bubble functions Bi , 0 ≤ i < n, given by
(
h−2 (x − xi )(xi+1 − x), xi ≤ x ≤ xi+1 ,
Bi (x) = (4.10)
0, otherwise.

It will be useful later in this section to have an expression for Bi of the form (4.7). From (4.8)
with q = Bi , it is clear that
(i)
αj = Bi (xj ) = 0, 0 ≤ j ≤ n. (4.11)

Also, since (
1
4
, j = i,
Bi (xj+1/2 ) =
0, j ̸= i,
from (4.9) with q = Bi we have
(
1
(i) (i) h
, j = i,
βj − βj+1 =
0, j ̸= i.
(i) (i)
We can solve this system uniquely in terms of β0 . If we take β0 = 1/h, the solution of this
system is (
1
(i) , 0 ≤ j ≤ i,
βj = h (4.12)
0, i < j ≤ n.

67
4 Generalized Finite Elements

Thus, using (4.11) and (4.12) in (4.7), we get

i
1X
Bi (x) = ϕj (x)(x − xj ). (4.13)
h j=0

The above expression for Bi is of the form (4.8) and thus Bi is a linear combination of the
shape functions ηij of SGFEM .

Remark. The shape functions of SGFEM may be linearly dependent giving rise to a singular
linear system. This can easily be seen in Example 4.6 where there are 2(n + 1) shape
functions in SGFEM given by ηij = ϕi (x)(x − xi )j , j = 0, 1, 0 ≤ i ≤ n. But the dimension of
SFEM = SGFEM is 2n + 1. Thus, the number of shape functions in SGFEM is greater than its
dimension. The solution of singular or ill-conditioned linear systems obtained from GFEM
will be considered in Sect. 4.4.

Remark. We recall that the functions in the local approximation space Vi for i, for which
ω i ∩ ∂Ω ̸= ∅, must satisfy the Dirichlet boundary conditions on ω i ∩ ∂Ω. In the 1d-setting,
if the exact solution u of a BVP satisfies the boundary conditions u(0) = 0 at x = 0, we
take α0 = 0 and the functions in V0 satisfy the boundary condition at x = 0. Likewise, if
u(1) = 0 is the specified boundary condition at x = 1, we take αn = 0 and

Vn = span{x − 1};

the functions in Vn satisfy the boundary condition at x = 1. A minor modification of the

above analysis shows that SGFEM = SFEM in this case also.

Example 4.7. Consider the domain Ω = (0, 1)2 and for n ∈ N let xi = ih, yj = jh, where
h = 1/n and 0 ≤ i, j ≤ n. We subdivide Ω into the squares Ωij = (xi , xi+1 ) × (yj , yj+1 ),
0 ≤ i, j < n. The nodes of this polyhedrization Ωh are pij = (xi , yj ), 0 ≤ i, j ≤ n. A
standard FEM space with respect to Ωh is

SFEM = {v ∈ C(Ω) : v|Ωij ∈ Qk },

where Qk = span{xl y m }kl,m=0 denotes the space of polynomials of degree at most k in each
variable. It is possible to find a GFEM space SGFEM with suitably chosen partition of unity
functions {ϕij } and local approximation spaces Vij so that SGFEM = SFEM . We again do this
for k = 2. In this case, the functions in SFEM are continuous piecewise biquadratics. We
construct a GFEM space as follows: To each node pij we associate a function

ϕij (x, y) = ϕi (x)ϕj (y), (4.14)

where ϕi and ϕj are one-dimensional hat functions centered at xi and yj , respectively, as

discussed in Example 4.6. ϕij is the standard piecewise bilinear hat function centered at pij
satisfying ϕij (pij ) = 1 and ϕij is zero at every other node. We let ω ij = supp ϕij =
[xi−1 , xi+1 ] × [yj−1 , yj+1 ]. We note that when i = 0 or j = 0, we replace xi−1 by xi or
yi−1 by yi , accordingly, in the definition of ωij . Similarly, when i = n or j = n, we replace
xi+1 by xi or yi+1 by yi , accordingly. Then {ϕij } satisfies (i)–(iv), in particular, they are a
partition of unity on Ω. For local approximation spaces Vij , 0 ≤ i, j ≤ n, we take

Vij = span{(x − xi )l (y − yj )m , l = 0, 1, m = 0, 1}.

68
 4.2 GFEM and Classical FEM

Thus Vij is the space of bilinear functions defined on ωij . We now define the GFEM space
as n
X
SGFEM = {ψ = ϕij lij }, (4.15)
i,j=0

where lij (x, y) = aij +bij (x−xi )+cij (y −yj )+dij (x−xi )(y −yj ) ∈ Vij , aij , bij , cij , dij ∈ R. We
note that lij is defined only on ωij , but since ϕij |∂ωij = 0, ϕij lij has a continuous extension
to Ω. Thus SGFEM is equivalently given by
SGFEM = span{ϕij , (x − xi )ϕij , (y − yj )ϕij , (x − xi )(y − yj )ϕij , i, j = 0, . . . , n}. (4.16)

We now show that SFEM = SGFEM . Since the functions ϕij lij are continuous in Ω, it is clear
from (4.15) that the functions in SGFEM are continuous in Ω. Also since ϕij , lij are bilinear
on each rectangle of Ωh , the functions in SGFEM are continuous piecewise biquadratics, and
hence SGFEM ⊂ SFEM .

It remains to show that SFEM ⊂ SGFEM . We will do this by proving that every element of a
basis of SFEM is contained in SGFEM . For k = 2 a well-known basis of SFEM consists of the
following functions, which can be grouped into four categories:

(a) The functions ϕij corresponding to the nodes pij , 0 ≤ i, j ≤ n.

(1)
(b) The functions Sij corresponding to the line segments (pij , pi+1,j ), 0 ≤ i < n, 0 ≤ j ≤ n,
defined by
(1)
Sij (x, y) = Bi (x)ϕj (y). (4.17)
Here, Bi is the one-dimensional quadratic bubble defined in (4.10). We note that for
(1)
j ≥ 1, supp Sij = [xi , xi+1 ] × [yj−1 , yj+1 ]. For j = 0, the support is [xi , xi+1 ] × [y0 , y1 ].
(2)
(c) The functions Sij corresponding to the line segments (pij , pi,j+1 ), 0 ≤ i ≤ n, 0 ≤ j < n,
defined by
(2)
Sij (x, y) = ϕi (x)Bj (y). (4.18)
(2)
We note that for i ≥ 1, supp Sij = [xi−1 , xi+1 ] × [yj , yj+1 ]. For i = 0, the support is
[x0 , x1 ] × [yi , yi+1 ].

(d) The functions Bij corresponding to the rectangles Ωij , 0 ≤ i, j < n, defined by
Bij (x, y) = Bi (x)Bj (y). (4.19)
We note that supp Bij = [xi , xi+1 ] × [yj , yj+1 ].

It is immediate from (4.16) that ϕij ∈ SGFEM for 0 ≤ i, j ≤ n. Using (4.17), (4.13), and
(4.14), we have
(1)
Sij (x, y) = Bi (x)ϕj (y)
i
1X
= ϕl (x)(x − xl )ϕj (y)
h l=0
i
1X
= (x − xl )ϕlj (x, y)
h l=0

69
4 Generalized Finite Elements

(1)
and therefore from (4.16), we have Sij ∈ SGFEM , 0 ≤ i < n, 0 ≤ j ≤ n. Similarly, using
(4.18), (4.13), (4.14), and (4.16), we have

j
(2) 1X
Sij (x, y) = (y − yl )ϕil (x, y) ∈ SGFEM
h l=0

for 0 ≤ i ≤ n and 0 ≤ j < n. Finally, from (4.19), (4.13), (4.14), and (4.16) we have

Bij (x, y) = Bi (x)Bj (y)

j
" i
#" #
1 X 1 X
= ϕl (x)(x − xl ) ϕm (y)(y − ym )
h l=0 h m=0
i j
1 XX
= 2 (x − xl )(y − ym )ϕlm (x, y) ∈ SGFEM
h l=0 m=0

for 0 ≤ i, j < n. Thus we have shown that all the basis elements for SFEM belong to SGFEM .
Therefore, SFEM = SGFEM .

Note that not for any domain Ω, every FEM relative to every polyhedrization of Ω can be
cast in the framework of GFEM.

In the case that SFEM = SGFEM , Theorems 4.4 and 4.5 can also be seen in the context of
FEM. However, application of these theorems does not yield the well-known error estimates
for FEM. In Example 4.6, SFEM is the space of continuous piecewise quadratic polynomials.
It is well known that
∥u − uFEM ∥E(Ω) ≲ h2 ∥u∥H 3 (Ω) ,
where uFEM denotes the Galerkin projection onto SFEM . Here, u is a smooth solution of
an elliptic linear Dirichlet BVP posed on Ω = (0, 1) with u(0) = u(1) = 0. We apply
Theorem 4.4 or Theorem 4.5 to obtain an error estimate. To this end, we choose ξiu ∈ Vi ,
0 ≤ i ≤ n, such that
∥u − ξiu ∥E(ωi ) ≲ h∥u∥H 2 (ωi ) =: εi .
RecallPthat Vi = span{1, x − xi }, 0 < i < n, V0 = span{x}, and Vn = span{x − 1}. Let
ξ u = ni=0 ϕi ξiu . It is easy to check that the assumptions of Theorem 4.5 hold in this example.
Thus, we obtain
n
X n
X
∥u − ξ˜u ∥2E(Ω) ≲ ε2i ≲h 2
∥u∥2H 2 (ωi ) ≲ h2 ∥u∥2H 2 (Ω) .
i=0 i=0

Thus, using ∥u − uGFEM ∥E(Ω) = inf ξ∈SGFEM ∥u − ξ∥E(Ω) , we have

∥u − uGFEM ∥E(Ω) ≤ ∥u − ξ˜u ∥E(Ω) ≲ h∥u∥H 2 (Ω) .

This loss of power of h can be explained as follows: The only assumption on the partition of
unity functions
Pn {ϕi } are (i)–(iv). It was not assumed that {ϕi } reproduce linear polynomials,
i.e., that i=0 xi ϕi (x) = x for all x ∈ Ω. Assuming the latter property, one can prove the
desired O(h2 ) error estimate.

70
 4.3 Choice of Local Approximation Spaces

4.3 Choice of Local Approximation Spaces

The choice of local approximation spaces Vi can be governed by available information on the
exact solution u. In this section we discuss some types of available information and show
how it can be used in the process of selecting Vi .

Information in terms of Sobolev spaces

We assume that the only available information on u is that it lies in H m (ωi ) and
 1/2
Z X (m)
∥u∥H m (ωi ) =  |Dk u|2 dx ≤ Ki , m = 0, 1, . . . , i = 1, . . . , N. (4.20)
ωi |k|≤m

Here, k = (k1 , k2 ), ki ≥ 0, and |k| = k1 + k2 . Our aim is to choose the spaces Vi so that

sup inf ∥u − ξi ∥E(ωi )

(m) ξi ∈Vi
u∈E(ωi ), ∥u∥H m (ωi ) ≤Ki

is small. In this situation, the space of polynomials of degree at most p is a good choice
(p)
for Vi . Denote this space by Vi = Wi . Then for m ≥ 1
(m)
Ki min{p,m−1}
εi ≤ c h ,
pm−1 i

where hi = diam ωi and c is independent of u, h, p, and m.

Information in terms of the BVP

In addition to the smoothness property (4.20), we often know more about the boundary
value problem. For example, if u is the solution of (4.1) with α = 1 and f = 0, then u
is a harmonic function. In this situation, we use harmonic polynomials instead of general
polynomials. Let
(p),H (p)
Wi = {v ∈ Wi : v is harmonic in ωi }.
Suppose ωi is star-shaped with respect to a ball1 and ∂ωi is piecewise analytic with internal
(p),H
angles αi = βi π, where 0 ≤ βi < 2 − λ and λ > 0. Then with shape functions in Wi it is
known that  (2−λ)(m−1)
(m) m−1 log p
εi ≤ cKi hi , p ≥ m − 1, m ≥ 1,
p
where c is independent of u but depends on the shape of ωi . We note that the rate of
convergence w.r.t. p depends on the angle of corners of the boundary.
1
A domain is called star-shaped (with respect to a point x0 ) if the segment [x, x0 ] ⊂ Ω for all x ∈ Ω and
star-shaped with respect to a ball if it is star-shaped with respect to each point of the ball.

71
4 Generalized Finite Elements

(p),H (p)
Remark. The dimension of Wi is 2p+1, whereas the dimension of Wi is 12 (p+1)(p+2).
Hence, for a given asymptotic rate of convergence, the space of harmonic polynomials has a
smaller number of degrees of freedom than the space of standard polynomials.

In order to handle general (analytic) diffusion coefficients α, one can rely on generalized
harmonic polynomials; see the monograph by Vekua, 1967.

If the right-hand side is nonzero, then we have to add additional shape functions. For
example, if f = 1, we add the shape function ξ = x2 + y 2 .

Selection of Vi when u has singularities

The solution u of (4.1) can be singular due to one or more of the following reasons:

1. the boundary ∂Ω has corners;

2. the boundary conditions changes, e.g., from Dirichlet to Neumann;

3. the coefficient α is rough, i.e., it is piecewise constant;

4. the right-hand side is not smooth.

We address only the first two cases. The character of the singular behavior of the solution
is well known. We will assume that ∂Ω has a corner A located at the origin and that the
boundary ∂Ω near A consists of two straight lines. This assumption is only for the sake of
simplicity. If f is sufficiently smooth, then in the neighborhood of the origin
s
X
u(r, φ) = ak rλk (log r)µk ψk (φ) + ζ(r, φ), (4.21)
k=0

where λk+1 ≥ λk , µk+1 ≥ µk , ψk is a smooth function of φ and ζ(r, φ) is smoother than any
of the terms in the sum. Here, (r, φ) are polar coordinates. We note that (4.21) is also true
if A is a point where a Dirichlet condition changes to Neumann.

With this information we select the shape functions in Vi to be the functions rλk (log r)µk ψk (φ),
k = 0, . . . , s, together with polynomials. Then the error of the approximation of u by func-
tions in Vi is only the approximation of ζ(r, φ) by polynomials.

Remark. There is a large literature on expansions of the form (4.21). In particular, an

expansion similar to (4.21) is valid for elasticity and Maxwell problems.

Construction of these singular functions may not be easy, especially in elasticity problems.
Hence, a numerical treatment is unavoidable. Either we solve the associated (local) prob-
lem numerically (with the GFEM) or use analytic formulae with numerically determined
parameters.

The use of special shape functions in ωi for which A ∈ ω i is crucial. Also, we have to use some
of the special shape functions in patches ωi which A ̸∈ ω i but ωi is close to A. The number

72
 4.4 Implementational Issues in the GFEM

of special shape functions needed depends on the accuracy requirement. Determining the
optimal number of terms as well as in which elements special shape functions are needed is
not simple. Usually, two terms in patches ωi for which A ∈ ω i and one term in all ωi that
are the direct neighbors of these patches is sufficient.

4.4 Implementational Issues in the GFEM

Implementation of the GFEM consists of five major parts, namely:

(a) the selection of local approximating functions;

(b) the selection of partition of unity (PU) functions;

(c) the construction of the stiffness matrix;

(d) the solution of the linear system; and,

(e) the computation of data of interest.

(a) We have already discussed the selection of local approximating functions in the previ-
ous section, which depends on the available information on the unknown solution u of the
problem (4.1).

(b) The primary role of PU functions {ϕi } in GFEM is to paste together the local approx-
imation functions to form global approximation functions that are conforming, i.e., global
approximation functions that are in E0 (Ω). In theory, any PU satisfying (i)–(iv) will do; we
may consider Shepard functions with disks as their supports or finite element hat functions,
or any family of particle shape functions used in meshless methods.

But the choice of patches {ωi } and the associated PU functions {ϕi } affects many aspects of
the implementation of GFEM, e.g. (c) and (d). We first discuss the effect of patches and of
PU functions in the construction of the stiffness matrix. A typical element of the stiffness
matrix is of the form Z
α ∇ηkl · ∇ηij dx. (4.22)
ωi ∩ωk

Since these integrals are evaluated by numerical integration, it is important to choose {ωi }
such that {ωi ∩ ωk }, 1 ≤ i, k ≤ N , are simple domains, in which numerical integration can be
performed efficiently. For example, if the ωi are discs (in R2 ) or balls (in R3 ), then ωi ∩ ωk is
a lens-shaped domain, and accurate numerical integration over such domains is known to be
difficult. If the ωi are chosen to be rectangles, then ωi ∩ωk is also a rectangle. It is much easier
to perform numerical integration on rectangular domains. Moreover, since ηij = ϕi ξij , the
integrand in (4.22) has terms involving {ϕi } and {∇ϕi }, and thus the numerical evaluation
of (4.22) depends also on the smoothness of the PU function {ϕi } and their derivatives.

The choice of the PU functions {ϕi } also affects the linear system Ac = f . We have mentioned
that the shape functions of SGFEM could be linearly dependent or independent, depending
on the choice of PU functions. This, in turn, leads to either a singular or a non-singular

73
4 Generalized Finite Elements

system. We further note that the condition number of the stiffness matrix depends on the
choice of the PU functions.

(d) We now comment on solving the linear system Ac = f . We have mentioned before that
the stiffness matrix A could be positive semi-definite or severely ill-conditioned. When A is
positive semi-definite, the system has non-unique solutions. The lack of unique solvability
of Ac = f does not imply that the GFEM has non-unique solutions.

A solution of Ac = f can be obtained either with a specialized direct solver or an iterative

solver. We describe an iterative scheme to solve Ac = f . We first perturb the matrix A
by εI, where ε > 0 is small. Let
Aε = A + εI.
Clearly, Aε is positive definite. We first compute
c0 = A−1
ε f,
r0 = f − Ac0 ,
z0 = A−1
ε r0 ,
v0 = Az0 .
Then, for i = 1, 2, . . . , we compute
ci = ci−1 + zi−1 ,
ri = ri−1 − vi−1 ,
zi = A−1
ε ri ,
vi = Azi ,
until the ratio
|ziT Azi |
|cTi Aci |
is sufficiently small, which is attained, say, for i = I. Then cI is considered a solution
of Ac = f . In practice, the above ratio becomes sufficiently small in one or two steps.

Lemma 4.8. We have that

 i  i
ε ε
∥zi ∥A ≤ ∥z0 ∥A and ∥f − Aci ∥ ≤ ∥r0 ∥,
λ2 + ε λ2 + ε

where 0 = λ1 < λ2 ≤ . . . denote the eigenvalues of A.

Proof. We observe
ri+1 = ri − vi = ri − AA−1 −1 −1
ε ri = (Aε − A)Aε ri = εAε ri . (4.23)
Similarly, zi+1 = (I − A−1
ε A)zi = εAε zi . Let (λi , ηi ) be the eigenpairs of A ∈ R
−1 n×n
such
that P
0 = λ1 < λ2 ≤ . . . and ηi are pairwise orthonormal.
Pn εβjIf c is a solution of Ac = f , then
f = nj=2 βj ηj . As a consequence, r0 = εA−1
ε f = j=2 λj +ε ηj and with (4.23)

n
X (i)
ri = αj ηj .
j=2

74
 4.4 Implementational Issues in the GFEM

Hence, we obtain
n (i) n
X (αj )2 ε2 X (i) 2 ε2
T
ri+1 ri+1 = ε2 riT A−1 −1
ε Aε ri =ε 2
2
≤ 2
(αj ) = rT r
2 i i
j=2
(λj + ε) (λ2 + ε) j=2
(λ2 + ε)

and similarly zi+1

T
Azi+1 ≤ ε2 /(λ2 + ε)2 ziT Azi . Furthermore, it is proved by induction that
f − Aci = ri . We obtain
 i
ε ε
∥f − Aci ∥ = ∥ri ∥ ≤ ∥ri−1 ∥ ≤ ∥r0 ∥
λ2 + ε λ2 + ε

and ziT Azi ≤ [ε/(λ2 + ε)]2i z0T Az0 .

(e) With the vector c, the solution at a particular point x ∈ Ω is obtained as

X mi
N X mi
N X
X
uGFEM (x) = cij ηij (x) and ∇uGFEM (x) = cij ∇ηij (x).
i=1 j=1 i=1 j=1

We note that the computation of ηij and ∇ηij involves computation of ϕi , ξij , ∇ϕi , and ∇ξij .

We show that P thePnon-uniqueness of c does not mean that uGFEM is non-unique. Let
u′GFEM (x) = N i=1
mi
c ′
j=1 ij ∇η ij (x) be a second solution resulting from another solution c′
auf Ac = f . Define J : RM → SGF EM with M := N i=1 mi by Jx := j=1 xij ηij . Then
P PN Pmi
i=1
the kernel of J coincides with the kernel of A. Since c′ − c ∈ Ker A = Ker J, we have

uGFEM (x) = Jc(x) = Jc′ (x) = u′GFEM (x).

75