Volker John
Numerical Methods for Partial Differential Equations
1. Some Partial Differential Equations From Physics...................................3
2. Finite Difference Methods for Elliptic Equations...................................17
3. Introduction to Sobolev Spaces.....................................................47
4. The Ritz Method and the Galerkin Method........................................67
5. Finite Element Methods...........................................................77
6. Interpolation.......................................................................97
7. Finite Element Methods for Second Order Elliptic Equations......................109
References...........................................................................121
1
Chapter 1
Some Partial Differential Equations From Physics
Remark 1.1. Contents. This chapter introduces some partial differential equations (pde's) from physics to show the importance of this kind of equations and to motivate the application of numerical methods for their solution. □
1.1 The Heat Equation
Remark 1.2. Derivation. The derivation follows (Wladimirow, 1972, p. 39). Let x = (x1,x2,x3)T G Q C K3, where Q is a domain, t G M, and consider the following physical quantities
• u(t,x) - temperature at time t and at the point x with unit [K],
• p(t, x) - density of the considered species with unit [kg/m3],
• c(t, x) - specific heat capacity of the species with unit [J/kg k] = [w s/kg k] ,
• k(t,x) - thermal conductivity of the species with unit [w/m k],
• F(t,x) - intensity of heat sources or sinks with unit [w/m3]. Consider the heat equilibrium in an arbitrary volume V C Q and in an
arbitrary time interval (t,t + At). First, there are sources or sinks of heat: heat can enter or leave V through the boundary dV, or heat can be produced or absorbed in V. Let n(x) be the unit outer normal at x G dV. Due to Fourier's1 law, one finds that the heat
enters through dV into V. One obtains with integration by parts (Gaussian
theorem)
■t+At
Q
i
1 Jean Baptisté Joseph Fourier (1768 - 1830)
3
4
1 Some Partial Differential Equations From Physics
In addition, the heat
r-t+At r
Q2 = /        / F{t,x) dx dt, [W s] = [J], Jt Jv
is produced in V.
Second, a law for the change of the temperature in V has to be derived. Using a Taylor series expansion, on gets that the temperature at x changes in (t, t + At) by
u{t + At,x) -u{t,x) = ^(i, x)At + 0((At)2).
Now, a linear ansatz is utilized, i.e.,
u(t + At,x) -u(t,x) = ^-{t,x)At.
With this ansatz, one has that for the change of the temperature in V and for arbitrary sufficiently small At, the heat
u(t + At, x) — u(t, x)
cp—--r-^--       dx dt
Jv At
r-t+At
r+Zit f du J       J cp—(t,x) dx dt, [J],
is needed. This heat has to be equal to the heat sources, i.e., it holds Q3 = Q2 + Qi, from what follows that
cp— - V • {kVu) - F
(t, x) dx dt = 0.
Since the volume V was chosen to be arbitrary and At was arbitrary as well, the term in the integral has to vanish. One obtains the so-called heat equation
cp— - V • (kVu) = F in (0, T) x Q.
At this point of modeling one should check if the equation is dimensionally correct. One finds that all terms have the unit [w/m3].
For a homogeneous species, c, p, and k are positive constants. Then, the heat equation simplifies to
^--e2Au = f m(0,T)xO, (1.1)
with e2 = k/(cp), [m%] and / = F/(cp), [K/s]. To obtain a well-posed problem, (1.1) has to be equipped with an initial condition u(0,x) and appropriate boundary conditions on (0,T) x dQ. □
1.1 The Heat Equation
5
Remark 1.3. Boundary conditions. For the theory and the numerical simulation of partial differential equations, the choice of boundary conditions is of utmost importance. For the heat equation (1.1), one can prescribe the following types of boundary conditions:
• Dirichlet2 condition: The temperature u(t,x) at a part of the boundary is prescribed
u = gi on (0,T) x dQo
with dQo C dQ. In the context of the heat equation, the Dirichlet condition is also called essential boundary conditions.
• Neumann3 condition: The heat flux is prescribed at a part of the boundary
du
-k— =g2 on (0,T) x dQN
with OOm C dQ. This boundary condition is a so-called natural boundary condition for the heat equation.
• Mixed boundary condition, Robin4 boundary condition: At the boundary, there is a heat exchange according to Newton's5 law
du
k— + h(u - uenv) = 0 on (0,T) x dQm, on
with dQm C dQ, the heat exchange coefficient h, [w/m2k2], and the temperature of the environment uenv.
□
Remark 1.4- The stationary case. An important special case is that the temperature is constant in time u(t, x) = u{x). Then, one obtains the stationary heat equation
-e2Au = f in J?. (1.2)
This equation is called Poisson6 equation. Its homogeneous form, i.e., with f{x) = 0, is called Laplace7 equation. Solution of the Laplace equation are called harmonic functions. The Poisson equation is the simplest partial differential equation. The most part of this lecture will consider numerical methods for solving this equation. □
Remark 1.5. Another application of the Poisson equation. The stationary distribution of an electric field with charge distribution f(x) satisfies also the Poisson equation (1.2). □
2 Johann Peter Gustav Lejeune Dirichlet (1805 —1859)
3 Carl Gottfried Neumann (1832 - 1925)
4 Gustave Robin (1855 - 1897)
5 Isaac Newton (1642 - 1727)
6 Simeon Denis Poisson (1781 - 1840)
7 Pierre Simon Laplace (1749 - 1829)
6
1 Some Partial Differential Equations From Physics
Remark 1.6. Non-dimensional equations. The mathematical analysis as well as the application of numerical methods relies on equations for functions without physical units, the so-called non-dimensional equations. Let
• L - a characteristic length scale of the problem, [m],
• U - a characteristic temperature scale of the problem, [K],
• T* - a characteristic time scale of the problem, [s].
If the new coordinates and functions are denoted with a prime, one gets with the transformations
from (1.1) the non-dimensional equation
9 ,tt ^9t' 9   ( 9 m ,,9<\9<     ,    . (tT\
T* dti       L2  ^d{xrf      J \    T* J
Usually, one denotes the non-dimensional functions like the dimensional functions, leading to
For the analysis, one sets L = 1 m, U = 1 K, and T* = 1 s which yields
^-e2Au = f in (0,T) x Q, (1.3)
with a non-dimensional temperature diffusion e2 and a non-dimensional right-hand side f(t, x).
The same approach can be applied to the stationary equation (1.2) and one gets
-s2Au = / in O, (1.4)
with the non-dimensional temperature diffusion e2 and the non-dimensional right-hand side f(x). □
Remark 1.7. A standard approach for solving the instationary equation. The heat equation (1.3) is an initial value problem with respect to time and a boundary value problem with respect to space. Numerical methods for solving initial value problems were topic of Numerical Mathematics 2.
A standard approach for solving the instationary problem consists in using a so-called one-step (9-scheme for discretizing the temporal derivative. Consider two consecutive discrete times tn and in+i with r = in+i — tn. Then,
1.1 The Heat Equation
7
Fig. 1.1 Solution of the two-dimensional example of Example 1.8.
the application of a one-step ^-scheme yields for the solution at £n+i Un+1 ~ Un - 8e2Aun+1 - (1 - 8)e2Aun = 8fn+1 + (1 - 8)fn,
T
where the subscript at the functions denotes the time level. This equation is equivalent to
un+1 - T6e2Aun+1 =un + r(l - 6)e2Aun + rdfn+1 + r(l - 6>)/„. (1.5)
For 8 = 0, one obtains the forward Euler scheme, for 8 = 0.5 the Crank-Nicolson scheme (trapezoidal rule), and for 8 = 1 the backward Euler scheme.
Given un, (1.5) is a boundary value problem for un+i. That means, one has to solve in each discrete time a boundary value problem. For this reason, this lecture will concentrate on the numerical solution of boundary value problems. □
Example 1.8. Demonstrations with the code MooNMD JOHN & MATTHIES (2004).
• Consider the Poisson equation (1.4) in Q = (0, l)2 with e = 1. The right-hand side and the Dirichlet boundary conditions are chosen such that u(x,y) = sin(7ra;) sin(7rj/) is the prescribed solution, see Figure 1.1 Hence, this solution satisfies homogeneous Dirichlet boundary conditions. Denote by Uh{x,y) the computed solution, where h indicates the refinement of a mesh in Q. The errors obtained on successively refined meshes with the simplest finite element method are presented in Table 1.1. One can observe in Table 1.1 that \\u — u^W^^^ converges with second order and |V(m — Mfe)II.l2(i?) converges with first order. A main topic of the numerical analysis of discretizations for partial differential equations consists in showing that the computed solution converges to the solution
8
1 Some Partial Differential Equations From Physics
Table 1.1 Example 1.8, two-dimensional example.
h degrees	of freedom	II'" - Uh\\L2(Q)	\\V(u-uh)\\L2(n)
1/4	25	8.522e-2	8.391e-l
1/8	81	2.256e-2	4.318e-l
1/16	289	5.726e-3	2.175e-l
1/32	1089	1.437e-3	1.089e-l
1/64	4225	3.596e-4	5.451e-2
1/128	16641	8.993e-5	2.726e-2
1/256	66049	2.248e-5	1.363e-2
1/512	263169	5.621e-6	6.815e-3
of an appropriate continuous problem in appropriate norms. In addition, to prove a certain order of convergence (in the asymptotic regime) is of interest.
• Consider the Poisson equation (1.4) in Q = (0, l)3 with e = 1 and / = 0. At z = 1 the temperature profile should be u(x,y, 1) = 16x(l—x)y(l—y) and at the opposite wall should be cooled u(x,y,0) = 0. At all other walls, there should be an undisturbed temperature flux §^(x,y, z) = 0. A approximation of the solution computed with a finite element method is presented in Figure 1.2.
The analytical solution is not known in this example (or it maybe hard to compute). It is important for applications that one obtains, e.g., good visualizations of the solution or approximate values for quantities of interest. One knows by the general theory that the computed solution converges to the solution of the continuous problem in appropriate norms and one hopes that the computed solution is already sufficiently close.
□
Fig. 1.2 Contour lines of the solution of the three-dimensional example of Example 1.8.
1.3 The Navier—Stokes Equations
1.2 The Diffusion Equation
9
Remark 1.9. Derivation. Diffusion is the transport of a species caused by the movement of particles. Instead of Fourier's law, Newton's law for the particle flux through dV per time unit is used
dQ = —DVu ■ n ds
with
• u(t,x) - particle density, concentration with unit [mol/m3],
• D(t,x) - diffusion coefficient with unit [m2/s]-
The derivation of the diffusion equation proceeds in the same way as for the heat equation. It has the form
c— - V • (DVu) +qu = F in (0, T) x Q, (1.6)
where
• c(t,x) - is the porosity of the species, [•],
• q{t,x) - is the absorption coefficient of the species with unit p/s],
• F(t,x) - describes sources and sinks, [mol/s m3].
The porosity and the absorption coefficient are positive functions. To obtain a well posed problem, an initial condition and boundary conditions are necessary.
If the concentration is constant in time, one obtains
- V • {DVu) +qu = F in O. (1.7)
Hence, the diffusion equation possesses a similar form as the heat equation.
□
1.3 The Navier—Stokes Equations
This section was not presented in the course. It is included in the lecture notes for students who are interested in.
Remark 1.10. Generalities. The Navier8-Stokes9 equations are the fundamental equations of fluid dynamics. In this section, a viscous fluid (with internal friction) with constant density (incompressible) will be considered. □
Remark 1.11. Conservation of mass. The first basic principle of the flow of an incompressible fluid is the conservation of mass. Let V be an arbitrary
8 Claude Louis Marie Henri Navier (1785 - 1836)
9 George Gabriel Stokes (1819 - 1903)
10
1 Some Partial Differential Equations From Physics
volume. Then, the change of fluid in V satisfies d
J pv ■ n ds = J V • {pv) dx,
dV V change flux through the boundary of V
where
• v(t,x) - velocity (vi,V2,vs)T at time t and at point x with unit [m/s],
• p - density of the fluid, [ks/m3].
Since V is arbitrary, the terms in the volume integrals have to be the same. One gets the so-called continuity equation
^ + V • {pv) = 0 in (0, T) x Q.
Since p is constant, one obtains the first equation of the Navier-Stokes equation, the so-called incompressibility constraint,
V -v = 0 in (0,T) x Í2. (1.8)
□
Remark 1.12. Conservation of linear momentum. The second equation of the Navier-Stokes equations represents Newton's second law of motion
net force = mass x acceleration.
It states that the rate of change of the linear momentum must be equal to the net force acting on a collection of fluid particles.
The forces acting on an arbitrary volume V are given by
Fv = J —Pn ds + J S'n ds + J pg dx,
ev ev v
outer pressure        friction gravitation
where
• S'{t,x) - stress tensor with unit [N/m2],
• P{t,x) - the pressure with unit [N/m2],
• g{t,x) - standard gravity (directed), [m/s2].
The pressure possesses a negative sign since it is directed into V, whereas the stress acts outwardly.
The integral on dV can be transformed into an integral on V with integration by parts. One obtains the force per unit volume
-VP + V • §' + pg.
1.3 The Navier—Stokes Equations
11
On the basis of physical considerations (Landau & Lifschitz, 1966, p. 53) or (John, 2016, Chapter 2), one uses the following ansatz for the stress tensor
§' = j](vv + VvT -       •       + C(V • v)I,
where
• rj - first order viscosity of the fluid, [ks/m s],
• ( - second order viscosity of the fluid, [ks/m s],
• I - unit tensor.
For Newton's second law of motion one considers the movement of particles with velocity v(t,x(t)). One obtains the following equation
\/>.\.p dv(t'xit))
dt
force per unit volume        mass per unit volume
acceleration
= p (dtv + (v ■ V)v).
The second formula was obtained with the chain rule. The detailed form of the second term is
/ VidxVi + v2dyv1 + «3dz«i \ (v ■ V)v =    vi_dxv2 + v2dyv2 + v3dzv2 \vidxv3 + v2dvv3 + v3dzv3 J
If both viscosities are constant, one gets
9v „.       VP I/77     A „,„
— - vAv + (v ■ V)v +-= g + - [{ + C V V • v ,
at p p V3 /
where v = r]/p, [m2/s] is the kinematic viscosity. The second term on the right-hand side vanishes because of the incompressibility constraint (1.8). One obtains the dimensional Navier-Stokes equations
dv VP
— - vAv + (v ■ V)v H--=g,   V • v = 0 in (0, T) x Q.
Remark 1.13. Non-dimensional Navier-Stokes equations. The final step in the modeling process is the derivation of non-dimensional equations. Let
• L - a characteristic length scale of the problem, [m],
• U - a characteristic velocity scale of the problem, [m/s],
• T* - a characteristic time scale of the problem, [s].
Denoting here the old coordinates with a prime, one obtains with the transformations
x' v t'
X = —,    u = —,    t = -
V U T*
12
1 Some Partial Differential Equations From Physics
the non-dimensional equations
-^-dtu--^-Au+(u- V)« + Vp = /,   V-w = 0 in(0,T)xJ?, U 1 U Li
with the redefined pressure and the new right-hand side
P Lg
P{t,x) =       (t,x),   f(t,x) = —(t,x).
This equation has two dimensionless characteristic parameters: the Strouhal10 number St and the Reynolds 11 number Re
U1 v
Setting T* = L/U, one obtains the form of the incompressible Navier-Stokes equations which can be found in the literature
^ - Re^Au + (u ■ V)w + Vp = / in (0, T) x Q, V-w = 0 in (0,T) x Q.
□
Remark 1.14- About the incompressible Navier-Stokes equations. The Navier-Stokes equations are not yet understood completely. For instance, the existence of an appropriately defined classical solution for Q C M3 is not clear. This problem is among the so-called millennium problems of mathematics Fefferman (2000) and its answer is worth one million dollar. Also the numerical methods for solving the Navier-Stokes equations are by far not developed sufficiently well as it is required by many applications, e.g. for turbulent flows in weather prediction. □
Remark 1.15. Slow flows. Am important special case is the case of slow flows which lead to a stationary (independent of time) flow field. In this case, the first term in the in the momentum balance equation vanish. In addition, if the flow is very slow, the nonlinear term can be neglected as well. One gets the so-called Stokes equations
-Re'^Au + Vp = / in O, V • u = 0 in O.
□
10 Čeněk Strouhal (1850 - 1923)
11 Osborne Reynolds (1842 - 1912)
1.4 Classification of Second Order Partial Differential Equations
13
1.4 Classification of Second Order Partial Differential Equations
Definition 1.16. Quasi-linear and linear second order partial differential equation. Let Q C Kd, d G N. A quasi-linear second order partial differential equation denned on Q has the form
d
ajk{x)djdkU + F (x,u,diu,... ,ddu) = 0 (1.9)
or in nabla notation
V • (A(x)Vu) +F(x,u,d!U,..., ddu) = 0. A linear second order partial differential equation has the form
d
a,jk{x)djdkU + b(x) ■ Vm + c(x)u = F(x).
Remark 1.17. The matrix of the second order operator. If u[x) is sufficiently regular, then the application of the Theorem of Schwarz12 yields djdku(x) = dkdju(x). It follows that equation (1.9) contains the coefficient d3dku(x) twice, namely in a3k{x) and a,kj(x). For defmiteness, one requires that
a-jk(x) = akj{x).
Now, one can write the coefficient of the second order derivative with the symmetric matrix
/au{x) ■ ■ ■ aid(a;)\ A(x)=       :     •.. :
\adl(x) ■ ■ ■ add(x)J
All eigenvalues of this matrix are real and the classification of quasi-linear second order partial differential equations is based on these eigenvalues. □
Definition 1.18. Classification of quasi-linear second order partial differential equation. On a subset Q C Q let a be the number of positive eigenvalues of A(x), (3 be the number of negative eigenvalues, and 7 be the multiplicity of the eigenvalue zero. The quasi-linear second order partial differential equation (1.9) is said to be of type (a, j3,7) on Q. It is called to be
• elliptic on Q if it is of type (d, 0,0) = (0, d, 0),
Hermann Amandus Schwarz (1843 — 1921)
14
1 Some Partial Differential Equations From Physics
• hyperbolic on Q, if its type is (d — 1,1,0) = (1, d — 1,0),
• parabolic on Q, if it is of type (d — 1,0,1) = (0, d — 1,1).
In the case of linear partial differential equations, one speaks of a parabolic equation if in addition to the requirement from above it holds that
ra,nk(A(x), b(x)) = d
in fi. □
Remark 1.19. Other cases. Definition 1.18 does not cover all possible cases. However, the other cases are only of little interest in practice. □
Example 1.20. Types of second order partial differential equations.
• For the Poisson equation (1.4) one has a„ = — e2 < 0 and al3 = 0 for i 7^ j. It follows that all eigenvalues of A are negative and the Poisson equation is an elliptic partial differential equation. The same reasoning can be applied to the stationary diffusion equation (1.7).
• In the heat equation (1.3) there is besides the spatial derivatives also the temporal derivative. The derivative in time has to be taken into account in the definition of the matrix A. Since this derivative is only of first order, one obtains in A a zero row and a zero column. One has, e.g., a„ = —e2 < 0, i = 2,..., d + 1, an = 0, and al3 = 0 for i 7^ j. It follows that one eigenvalue is zero and the others have the same sign. The vector of the first order term has the form b = (1,0, ...,0)T G Md+1, where the one comes from dtu(t,x). Now, one can see immediately that (A, b) possesses full column rank. Hence, (1.3) is a parabolic partial differential equation.
• An example for a hyperbolic partial differential equation is the wave equation
dttu -e2Au = f in (0, T) x n.
□
1.5 Literature
Remark 1.21. Some books about the topic of this class. Books about finite difference methods are
• Samarskij (1984), classic book, the English version is Samarskii (2001)
• LeVeque (2007)
Much more books can be found about finite element methods
• Ciarlet (2002), classic text,
• Strang & Fix (2008), classic text,
• Braess (2001), very popular book in Germany, English version available,
1.5 Literature
15
• Brenner & Scott (2008), rather abstract treatment, from the point of view of functional analysis,
• Ern & Guermond (2004), modern comprehensive book,
• Grossmann & Roos (2007)
• Solin (2006), written by somebody who worked a lot in the implementation of the methods,
• Goering et al. (2010), introductory text, good for beginners,
• Deuflhard & Weiser (2012), strong emphasis on adaptive methods
• Dziuk (2010).
These lists are not complete.
These lectures notes are based in some parts on lecture notes from Sergej Rjasanow (Saarbrücken) and Manfred Dobrowolski (Würzburg). □
Chapter 2
Finite Difference Methods for Elliptic Equations
Remark 2.1. Model problem. The model problem in this chapter is the Poisson equation with Dirichlet boundary conditions
-Au = finfl, u = g on oil,
where il C M2. This chapter follows in wide parts Samarskij (1984). □
2.1 Basics on Finite Differences
Remark 2.2. Grid. This section considers the one-dimensional case. Consider the interval [0,1] that is decomposed by an equidistant grid
x% = ih,    i = 0,. .. ,n,    h = 1/n, - nodes, ujh = {xi : i = 0,.. ., n} - grid.
□
Definition 2.3. Grid function. A vector uh = (uq, ..., un)T G Mn+1 that assigns every grid point a function value is called grid function. □
Definition 2.4. Finite differences. Let v(x) be a sufficiently smooth function and denote by vt = v(xi), where xt are the nodes of the grid. The following quotients are called
17
2 Finite Difference Methods for Elliptic Equations
tangent
X^—\ Xi X^\
Fig. 2.1 Illustration of the finite differences.
_ Vj+1 - Vj
Vx'1 - h
forward difference,
- backward difference,
vs., =   1+1„,—- central difference,
2h
vl+1 - 2vl
second order difference,
see Figure 2.1.
Remark 2.5. Some properties of the finite differences. It is [exercise)
Vx,i       ~i^(.Vx,i ^^c,i)i     Vx~x^ (,Vx,i)x,i-
Using the Taylor series expansion for v(x) at the node x%, one gets (exer-
vx,i = v'{xi) + -hv"{xi) + O (h2)
ve,í = v'{xi) - ^hv"{Xi) + O (h2) v i,i =v'(xi) + 0(h2) , v^=v"(x,) + 0(h2).
2.1 Basics on Finite Differences 19
Definition 2.6. Consistent difference operator. Let I be a differential operator. The difference operator : Mn+1 —> Mn+1 is called consistent with L of order k if
max \(Lu)(xi) - (Lhuh)i\ = \\Lu - Lhuh\\oo,uh = O (hk)
0<i<n
for all sufficiently smooth functions u{x). □
Example 2.7. Consistency orders. The order of consistency measures the quality of approximation of L by L^.
The difference operators vx^,v^^,v%^ are consistent to L = ^ with order
1,1, and 2, respectively. The operator VxX,i is consistent of second order to ,2
L = ^3, see Remark 2.5. □
Example 2.8. Approximation of a more complicated differential operator by difference operators. Consider the differential operator
Lu = Tx{k{x)Tx)'
where k(x) is assumed to be continuously differentiable. Define the difference operator      as follows
(Lhuh)% = (aiij,,)^, = ^[a(x%+1)ux^(x%+1) - a(xi)ux,i(zj))
where a is a grid function that has to be determined appropriately. One gets with the product rule
(Lu)i = k'{xi){u')i + k{xi){u")i and with a Taylor series expansion for ut+i, which is inserted in (2.2),
{Lhuh)i =---(m )i H----(u )i H----(u )i + O (h ) .
Thus, the difference of the differential operator and the difference operator is (Lu)i - [Lhuh)i = (hHxi) - a^ha^        + (k(Xi) _ ai+12+a'j (u'% -h{ai+1^ai)(n"'h + 0(h2). (2.3)
In order to define Lh such that it is consistent of second order to L, one has to satisfy the following two conditions
20 2 Finite Difference Methods for Elliptic Equations
From the first requirement, it follows that al+i — a% = O (h). Hence, the third term in the consistency error equation (2.3) is of order O [h2). Possible choices for the grid function are [exercise)
Note that the 'natural' choice, at = kt, leads only to first order consistency.
(exercise) □
2.2 Finite Difference Approximation of the Laplacian in Two Dimensions
Remark 2.9. The five point stencil. The Laplacian in two dimensions is defined by
d2u d2u dx2 dy2
The simplest approximation uses for both second order derivatives the second order differences. One obtains the so-called five point stencil and the approximation
LA(l ^ J\_a — axx -\- LLyy — r ^2 '
(2-4)
see Figure 2.2. From the consistency order of the second order difference, it follows immediately that Au approximates the Laplacian of order 0{hl + h%). □
Remark 2.10. The five point stencil on curvilinear boundaries. There is a difficulty if the five point stencil is used in domains with curvilinear boundaries. The approximation of the second derivative requires three function values in each coordinate direction
[x - h~,y), (x,y), [x + h+,y), (x,y-h~), (x, y), (x, y + h+),
see Figure 2.3. A guideline of defining the approximation is that the five point stencil is recovered in the case h~ = A possible approximation of this type is
(-^i y h-)
Fig. 2.3 Sketch to Remark 2.10.
d2u ^ 1  f u(x + h+,y) -u(x,y)     u(x,y) - u(x - hx ,y)\
dx2     hx V ht hx
with hx = (hx + h~)/2. Using a Taylor series expansion, one finds that the error of this approximation is
d2u      1   fu{x + h+,y)-u{x,y)     u{x,y) - u{x - hx ,y)
dx2    hx V ht hx
For /i+ 7^ /i^. , this approximation is of first order. A different way consists in using
d2u      1  (u{x + h+,y) - u{x,y)     u{x,y) - u{x - hx ,y)
dx2     hx \ hi hx
with hx = max{/i+, hx}. However, this approximation possesses only the order zero, i.e., there is actually no approximation.
Altogether, there is a loss of order of consistency at curvilinear boundaries.
□
22
2 Finite Difference Methods for Elliptic Equations
Example 2.11. The Dirichlet problem. Consider the Poisson equation that is equipped with Dirichlet boundary conditions (2.1). First, K2 is decomposed by a grid with rectangular mesh cells xt = ihx,yj = jhy, hx,hy > 0, i, j G Z. Denote by
= {o}        inner nodes, five point stencil does not contain any boundary node,
= {*}        inner nodes that are close to the boundary, five point stencil contains boundary nodes, 7h = {*}        boundary nodes, LUh =      U to£ inner nodes, U Jh grid,
see Figure 2.4.
The finite difference approximation of problem (2.1) that will be studied in the following consists in finding a mesh function u{x) such that
—Au(x) = 4>(x) x G w^, -A*u(x) = (p(x) x G w*h, (2.6) u(x) = g(x) x G jh,
where <f>(x) is a grid function that approximates f(x) and A* is an approximation of the Laplacian for nodes that are close to the boundary, e.g., defined by (2.5). The discrete problem is a large sparse linear system of equations. The most important questions are:
• Which properties possesses the solution of (2.6)?
• Converges the solution of (2.6) to the solution of the Poisson problem and if yes, with which order?
□
2.3 The Discrete Maximum Principle for a Finite Difference Approximation
23
2.3 The Discrete Maximum Principle for a Finite Difference Approximation
Remark 2.12. Contents of this section. Solutions of the Laplace equation, i.e., of (2.1) with f(x) = 0, fulfill so-called maximum principles. This section shows, that the finite difference approximation of an operator, where the five point stencil of the Laplacian is a special case, satisfies a discrete analog of one of the maximum principles. □
Theorem 2.13. Maximum principles for harmonic functions. Let Q C
Md be a bounded domain and u G C2(Q)C\C(Q) be harmonic in Q, i.e., u(x) solves the Laplace equation —Au = 0 in Q.
• Weak maximum principle. It holds
maxti(a;) = max u{x).
That means, u{x) takes its maximal value at the boundary.
• Strong maximum principle. If Q is connected and if the maximum is taken in Q (note that Q is open), i.e., u(xg) = max^,^ u{x) for a point Xg G Q, then u{x) is constant
u(x) = maxti(a;) = u(xq)   V x G Q.
Proof. See the literature, e.g., (Evans, 2010, p. 27, Theorem 4) or the course on the theory of partial differential equations. ■
Remark 2.14- Interpretation of the maximum principle.
• The Laplace equation models the temperature distribution of a heated body without heat sources in Q. Then, the weak maximum principle just states that the temperature in the interior of the body cannot be higher than the highest temperature at the boundary.
• There are maximum principles also for more complicated operators than the Laplacian, e.g., see Evans (2010).
• Since the solution of the partial differential equation will be only approximated by a discretization like a finite difference method, one has to expect that basic physical properties are satisfied by the numerical solution also only approximately. However, in applications, it is often very important that such properties are satisfied exactly.
□
Remark 2.15. The difference equation. In this section, a difference equation of the form
a(x)u(x)=   ^  b(x,y)u(y) + F(x), x G ujh (J-yh, (2.7)
24 2 Finite Difference Methods for Elliptic Equations
Fig. 2.5 Grid that is not allowed in Section 2.3.
will be considered. In (2.7), for each node x, the set S(x) is the set of all nodes on which the sum has to be performed, but x G" S(x). That means, a(x) describes the contribution of the finite difference scheme of a node x to itself and b(x,y) describes the contributions from the neighbors.
It will be assumed that the grid u>h of inner nodes is connected, i.e., for all xa,xe G u>h exist x\,..., xm G u>h with x\ G S(xa), X2 G S(xi),..., xe G S(xm). E.g., the situation depicted in Figure 2.5 is not allowed.
It will be assumed that the coefficients a(x) and b(x, y) satisfy the following conditions:
a(x) > 0, b(x,y) > 0, V x G wfe,V y G S(x),
a(x) = 1, b(x,y) = 0 V x G "fh (Dirichlet boundary condition).
The values of the Dirichlet boundary condition are incorporated in (2.7) in the function F(x). □
Example 2.16. Five point stencil for approximating the Laplacian. Inserting the approximation of the Laplacian with the five point stencil (2.4) for x = (x,y) G     in scheme (2.7) gives
hlh%
—u(x,y) =
j^u(x + hx,y) + j^u(x - hx,y)
+ —^u{x,y + hy) + -r^u{x, y - hy)
-4>{x,y).
It follows that
a(x)
2(h2x + hp hlh?y '
b(x,y) G {hx2,hy2}, S{x) = {(a; - hx,y), {x + hx,y), {x,y - hy), {x,y + hy)}.
For inner nodes that are close to the boundary, only the one-dimensional case (2.5) will be considered for simplicity. Let x + hx G "fh, then it follows
2.3 The Discrete Maximum Principle for a Finite Difference Approximation 25 by inserting (2.5) in (2.7)
T~ )u(x'y> = —n~=--1--rr+--^(x), (2.8)
ft^   \/la;        /la; / hxtlx hxhj
on 7fc—>F(i)
such that a(z) = J-       +       , 6(x,y) = =±=, and 5(a;) = {(a; - h~,y)}.
X     \     X XX Q
Remark 2.17. Reformulation of the difference scheme. Scheme (2.7) can be reformulated in the form
d(x)u(x) =   J2  b(x,y)(u(y)-u(x))+F(x) (2.9)
withd(x)=a(x)-J2yes(a:)b(x,y)- D
Example 2.18. Five point stencil for approximating the Laplacian. Using the five point stencil for approximating the Laplacian, form (2.9) of the scheme is obtained with
, N     2{h2x + hl)     2      2 ,
nxny nx ny
for x G u°h.
The coefficients a(x) and b(x,y) are the weights of the finite difference stencil for approximating the Laplacian. A minimal condition for consistency is that this approximation vanishes for constant functions since the derivatives of constant functions vanish. It follows that also for the nodes x G ui^, it is a(x) = X^es(a;) b{x,y), compare (2.8). However, as it could be seen in Example 2.16, in this case the contributions from the neighbors on jh are included in the scheme (2.7) in F(x). Hence, one obtains for nodes that are close to the boundary
d(x)= J2 b(x,y)-    J2    b(x,y)=    J2    h^x^- (2-n)
yeS(a:) yeS(a:) ,ygfh yeS(a:),ye-fh
In the one-dimensional case, one has, by the definition of hx and with h~ =
hx > h+,
JM      1/1 1
d(x) =
hx J     hxhx      hxhx      hxhx -\- hx hx
2 1 > T-,---- = -r~r- > 0.
26
2 Finite Difference Methods for Elliptic Equations
Lemma 2.19. Discrete maximum principle (DMP) for inner nodes.
Let u{x)    const on cu^ and d(x) > 0 for all x € cu^. Then, it follows from
Lhu(x) := d(x)u(x) -   ^2  b(x>y){u(y) -u(x))< 0 (2.12)
(or Lhii(x) > 0, respectively) on cu^ that u{x) does not possess a positive maximum (or negative minimum, respectively) on cu^.
Proof The proof is performed by contradiction. Let L^u(x) < 0 for all x 6 oj^ and assume that u{x) has a positive maximum on uj^ at x, i.e., u(x) = max^^ u{x) > 0. For the node x, it holds that
Lhu(x) = d(x) u(x) -   ^2   H^,y)    (u(y)-u(x))     >d(x)u(x)>0. (2.13)
>0     >0       y£S(as)     >Q     ^ ^ definition of S
Hence, it follows that L^uix) = 0 and, in particular, that all terms of L^uix) have to vanish. For the first term, it follows that d(x) = 0. For the terms in the sum to vanish, it must hold
u(y) = u(x)   Vt/G S(x). (2.14)
From the assumption u{x) ^ const, it follows that there exists a node x 6 uj^ with u(x) > u{x). Because the grid is connected, there is a path x, x\, . . . ,xm,x such that, using (2.14) for all nodes of this path,
x\ 6 S(x),   u{x\) = u(x),
X2 G S(xi), u(x2) = u(x\) = u(x),
x e S(xrn), u(xrn) = u(xrn-i) = . . . = u(x) > u(x).
The last inequality is a contradiction to (2.14) for xm. ■
Corollary 2.20. DMP for boundary value problem. Let u(x) > 0 for x € jh and Lhu(x) < 0 (or Lhu(x) > 0, respectively) on u>h- Then, the grid function u{x) is non-positive (or non-negative, respectively) for all x € U 7fe.
Proof. Let L^uix) < 0 on oj^ u 7^. Assume that there is a node x 6 oj^ with u{x) > 0. Then, the grid function has either a positive maximum on uj^, which is a contradiction to the DMP for the inner nodes, Lemma 2.19, or u{x) has to be constant, i.e., u{x) = u(x) > 0 for all x 6 uj^. For the second case, consider a boundary-connected inner node x 6 uj^. Using the same calculations as in (2.13) and taking into account that the values of u at the boundary are non-positive, one obtains
Lhu(x) = d(x)u(x) - b(x,y) (u(y) - u(x))
b(y,y){u{y)-u{x))>o,
which is a contradiction to the assumption on Lyt. ■
Corollary 2.21. Unique solution of the discrete Laplace equation with homogeneous Dirichlet boundary conditions. The discrete Laplace
2.3 The Discrete Maximum Principle for a Finite Difference Approximation
27
equation Lhu(x) = 0 for x G Lüh U jh possesses only the trivial solution u(x) = 0.
Proof. The statement of the corollary follows by applying Corollary 2.20 and its analog for the non-positivity of the grid function if u{x) < 0 for x 6 7^ and Lhu{x) < 0 on uj^.
Corollary 2.22. Comparison lemma. Let
Lhu(x) = f[x) for x G lüh; u(x) = g(x) for x G jh, Lhu(x) = f(x) for x G Lüh;    u(x) = g(x) for x G jh,
with \f(x)\ < f(x) and \g(x)\ < g(x). Then, it is \u(x)\ < u(x) for all x G u>h U jh- The function u[x) is called majorizing function.
Proof. Exercise. ■
Remark 2.23. Remainder of this section. The remaining corollaries presented in this section will be applied in the stability proof in Section 2.4. In this proof, the homogeneous problem (right-hand side vanishes) and the problem with homogeneous Dirichlet boundary conditions will be analyzed separately.
□
Corollary 2.24. Homogeneous problem. For the solution of the problem
Lhu(x) =0,        X G Lüh,
u(x) = g(x), x G 7fe, with d(x) = 0 for all x G       it holds that
Wu\\l™(üjhu-,h) ^ Il9ll(°°(7fc) •
Proof. Consider the problem
Lhü{x) = 0, x 6 Lüh,
ü(x) = g(x) = const = ||ff||,oo(7fc) , x E ~fh-
It will be shown that u(x) = |(/||^oo= const by inserting this function in the problem.1 For inner nodes that are not close to the boundary, it holds that
Lhü(x)=    d(x)    u(x) - b(x,y) («(3/) - «0=)) = 0.
, yeS(ai) "---'
=0, (2.10) "      K   ' =0
With the same arguments as in Example 2.18, one can derive the representation (2.11) for inner nodes that are close to the boundary. Inserting (2.11) in (2.12) and using in addition ü{x) = u{y) yields
1 The corresponding continuous problem is —Au = 0 in Q, u = const = |«?11^00(^.^) on
It is clear that u = ||g||£oo(7/j is the solution of this problem. It is shown that the discrete
analog holds, too.
28
2 Finite Difference Methods for Elliptic Equations
Lhu(x) = d(x)u(x) -   ^   Kx,y) (u(y) -u(x)) =        ^ b(x,y)u(x) =        ^ b(x,y)ü(y).
y(ES(x),y(E-Yh
This expression is exactly the contribution of the nodes on 7^ that is included in F(x) in scheme (2.7), see also Example 2.16. That means, the finite difference equation is also satisfied by the nodes that are close to the boundary.
Now, the application of Corollary 2.22 gives Tt(x) > |'ti(a;)| for all x £ uj^ u 7^, such that
IMIi~(^u7fc) <ü(x) = NI,=o(7fc),
which is the statement of the corollary. ■
Corollary 2.25. Problem with homogeneous boundary condition.
For the solution of the problem
Lhu{x) = f(x), x G uh, u(x) = 0,     x £ jh,
with d(x) > 0 for all x € cu^, it is
IMIi°°(cfcU7fc) ^ \\D 1f\\l<*>(uh)
with D = diag(d(x)) for x € u>h-Proof Consider the grid function
J(x) = \f(x)\ > f(x)Vxeuh.
From the discrete maximum principle, it follows that the solution of the problem
Lhü(x) = f(x), x e ujh, ü(x) = 0,      x e lh,
is non-negative, i.e., it holds u(x) > 0 for x £ uj^ u7/J,. Define the node x by the condition
«(<=)= INIi-(lJfcu7fc)-
In x, it is
Lhü(x) = d(x)ü(x) - b(^, y) (ü(y) -ü(x)) = \f(x)\,
from what follows that
. 1/(^)1 . 1/(^)1
u(x) < —       < max —-—— = max d(x)      x^ojh d(x) x^ojh
= ||D-if|i
d(x)       " J»l^(^h)
Since u{x) < u{x) for all x 6 uj^ u 7^ because of Corollary 2.22, the statement of the corollary is proved. ■
Corollary 2.26. Problem with homogeneous boundary condition and inhomogeneous right-hand side close to the boundary. Consider
2.4 Stability and Convergence of the FD Approximation of the Poisson Problem 29
Lhu(x) = f(x), x G uih, u(x) = 0,     x € jh,
with f(x) = 0 for all leuj. With respect to the finite difference scheme, it will be assumed that d(x) = 0 for all x G cu^, and d(x) > 0 for all x G cu^. Then, the following estimate is valid
with D+ = diag(0,d(x)^1). The zero entries appear for x G lu^ and the entries d(x)^1 for lEuJ.
Proof Let fix) = \f(x)\, x 6 oj^, and g(x) = 0,x 6 7^. The solution u{x) is non-negative, u{x) > 0 for all x 6 uj^ u 7^, see the DMP for the boundary value problem, Corollary 2.25. Define x by
ü(x)= |N|,oo(lJfcU7fc).
One can choose x 6 uj^, because if x 6 uj^, then it holds that
d(x)ü(x) - b(ß,y) (ü(y) - u(x)) = fix) = 0,
^ yeS(W)^-^'---'
i.e., Tt(x) = u(y) for all y £ S^Ic). Let £ 6 a;^ and x,x\, . . . ,xm,x be a connection with Xi ^ o;^, 'i = 1, . .. , m. For xm, it holds analogously that
ü{xrn) = ||^||;oo(^hU7h) = ü(y) V y e S{xrn).
Hence, it follow in particular that u(x) = H^ll^oo^^y^^^ such that one can choose x = x. It follows that
d{x)        ü{x) -   ^   b(x, y) (ü(y) - ü(x)) = J(x).
^(*)^^-Z-'
Since all terms in the sum over y 6 uj^ are non-negative, it follows, using also Corollary 2.22, that _ _
fix) fix)
IMI;co,    ,,   1 < ||«||,00,    ,,   1 < :L^-L < max <   D+f I    ,    . .
2.4 Stability and Convergence of the Finite Difference Approximation of the Poisson Problem with Dirichlet Boundary Conditions
Remark 2.27. Decomposition of the solution. A short form to write (2.6) is
Lhu(x) = f(x), x G uih,   u(x) = g(x), x G jh. The solution of (2.6) can be decomposed into
30
2 Finite Difference Methods for Elliptic Equations
u(x) = Ui(x) + u2(x),
with
LfrUi{x) = f(x), x G LUh, Ui(x) = 0, xE)), (homogeneous boundary cond.), LhU2(x) = 0, x G a;^,   U2(x) = g(a;), x e 71 (homogeneous right-hand side).
□
Stability with Respect to the Boundary Condition
Remark 2.28. Stability with respect to the boundary condition. From Corollary 2.24, it follows that
ll«2||,»(ll,h)<||9||,»(7h). (2.15)
□
Stability with Respect to the Right-Hand Side
Remark 2.29. Decomposition of the right-hand side. The right-hand side will be decomposed into
f(x) = r(x) + r(x)
with
fW = { f0{Xh I I f*(x) = f{x) - f°(x).
Since the considered finite difference scheme is linear, also the function ui(x) can be decomposed into
U\ (x) = u\ (x) + u\ (x)
with
Lhu°{x) = f°(x), x G uh, u\(x) = 0, x G jh, Lhu[(x) = f*{x), x G uih,   u[(x) = 0, x G 7^.
□
Remark 2.30. Estimate for the inner nodes. Let S((0, 0),R) be a circle with center (0,0) and radius R, which is chosen such that R > \\x\\2 for all x G Q. Consider the function
u(x) = a [R2 - x2 - y2)   with a > 0,
2.4 Stability and Convergence of the FD Approximation of the Poisson Problem 31
that takes only non-negative values for (x,y) G Q. Applying the definition of the five point stencil, it follows that
Au(x) = -aA(x2 +y2 - R2)
(x + hx)2 - 2x2 + (x - hx)2     (y + hy)2 - 2y2 + (y - hy)2
hi
h2
and
-4a =: -f{x), x G ui°h,
1  ((x + K)2-x2 x2-(x-h-)2
A*u(x) = —a
hx
l ({y + Kf-y2   y2 - {y-Kf
hV
-f(x), x G ui*h.
Hence, u(x) is the solution of the problem
Lhü{x)=J{x), x€ujh, Ti(x) = a [R2 — x2 — y2) > 0, x G jh-
It is u(x) > 0 for all x G "fh- Choosing a = j \\f°\\i<x>(uh)i one obtains
7(*) = 4a=|iril^fc)>|r0z)|, X€u°h,
J(x)>0=\r(x)\, XGLÜ*h.
Now, Corollary 2.22 (Comparison Lemma) can be applied, which leads to
R2
hW^M ^ ll"lli»(Uh) < aR = — ll/°lli»(Uh) '
(2.16)
One gets the last lower or equal estimate because (0,0) does not need to belong to Q or uih- □
Remark 2.31. Estimate for the nodes that are close to the boundary. Corollary 2.26 can be applied to estimate u*(x). For x G ui%, it is d(x) = 0, see Example 2.18. For x G ui£, one has
d(x) = a (a;) -      ^      b(x, y) > j-,
h?
with h = ma,x{hx, hy}, since all terms in the sum are of the form
32 2 Finite Difference Methods for Elliptic Equations
1111
hxhx     hxhx     foy^y ^y^y
see Example 2.18. One obtains
Kii«»(Uh)<l^+r||,»(11,h)<ft2iirii,»(Uh). (2.i7)
□
Lemma 2.32. Stability estimate The solution of the discrete Dirichlet problem (2.6) satisfies
r2 ..... -2
i~(Uhu7h) < Il9ll,»(7h) + — + V IMIi~K) (2-18)
with R > \\x\\2 for all x € Q and h = max{hx, hy}, i.e., the solution u{x) can be bounded in the norm \\'\\i^^hUlh^ by the data of the problem.
Proof. The statement of the lemma is obtained by combining the estimates (2.15), (2.16), and (2.17). ■
Convergence
Theorem 2.33. Convergence. Let u(x) be the solution of the Poisson equation (2.1) and Uh{x) be the finite difference approximation given by the solution of (2.6). Then, it is
\\u-uh\\^(u]hUlh)<Ch2
with h = max{hx,hy}.
Proof. The error in the node (a^, yj) is defined by eij = u(xi, yj) — U}1(xi, yj). With
-Au(xi,yj) = -Au(xi,yj) + O (h2) = f{xi,yj) + 0(h2) ,
one obtains by subtracting the finite difference equation, the following problem for the error
-Ae(x) = ip(x), x e w°, ip(x) = O (h2) , -A*e(x) = if>(x), x e w*, ip(x) = 0(1), e(x) =0,       x e lh:
where ip(x) is the consistency error, see Section 2.2. Applying the stability estimate (2.18) to this problem, one obtains immediately
2.5 An Efficient Solver for the Dirichlet Problem in the Rectangle
33
1---
hx
0 1 n,
Fig. 2.6 Grid for the Dirichlet problem in the rectangular domain.
2.5 An Efficient Solver for the Dirichlet Problem in the Rectangle
Remark 2.34- Contents of this section. This section considers the Poisson equation (2.1) in the special case Q = (0,lx) x (0,lv). In this case, a modification of the difference stencil in a neighborhood of the boundary of the domain is not needed. The convergence of the finite difference approximation was already established in Theorem 2.33. Applying this approximation results in a large linear system of equations Au = / which has to be solved. This section discusses some properties of the matrix A and it presents an approach for solving this system in the case of a rectangular domain in an almost optimal way.
A number of result obtained here will be need also in Section 2.6. □
Remark 2.35. The considered problem and its approximation. The considered continuous problem consists in solving
-Au = fmO = (0,lx)x (0,y, u = g on dO,
and the corresponding discrete problem in solving
—Au(x) = (f>(x), x G uih, u(x) = g(x), x G 7fe,
where the discrete Laplacian is of the form (for simplicity of notation, the subscript h is omitted)
Au = ~ 2UJ + ^ + ^+1 ~ ^ + ^ =: Axu + Ayu, (2.19)
with hx = lx/nx,hy = ly/ny, i = 0,... ,nx, j = 0,...,%, see Figure 2.6. □
Remark 2.36. The linear system of equations. The difference scheme (2.19) is equivalent to a linear system of equations Au = /.
For assembling the matrix and the right-hand side of the system, usually a lexicographical enumeration of the nodes of the grid is used. The nodes are
34
2 Finite Difference Methods for Elliptic Equations
called enumerated lexicographically if the node (ii,ji) has a smaller number than the node {12,32)1 if f°r the corresponding coordinates, it is
Vi < s/2 or (j/j = j/2) a (zi < x2).
Using this lexicographical enumeration of the nodes, one obtains for the inner nodes a system of the form
A = BlockTriDiag(C, B, C) G K("x-i)K-i)x(nI-i)(n!)-i))
B = ™i««(-4'4 + 4'"4)eR(B""1)X(B""1)'
C = Diag (y-J^) e K^"1^^-1),
' (f>(x), lEw
,/ n , g{x±hx,y)
<p(a;) H---, a; G a>£, close to lower
or upper boundary,
,. . g[x,y±hy)
<p(x) H--———, x G 0}^,, close to left
h2
or right boundary,
,. .     g(x±hx,y) g(x,yx±hv)
<p(x) H--7.--1--7.-—, x G Lit, corner of inner nodes.
In this approach, the known Dirichlet boundary values are already substituted into the system and they appear in the right-hand side vector. The matrices B and C possess some modifications for nodes that have a neighbor on the boundary.
The linear system of equations has the following properties:
• high dimension: N = (nx — 1)(% — 1) ~ 103 • • • 107,
• sparse: per row and column of the matrix there are only 3, 4, or 5 non-zero entries,
• symmetric: hence, all eigenvalues are real,
• positive definite: all eigenvalues are positive. It holds that
Amin = A(M) ~ ^ (L + ^=o(l),
Amax = \nj!-i,ny-i) ~ vr2 (jp + ^ = O {h-2) , (2.20)
with h = m.ax.{hx,hy\, see Remark 2.37 below.
• high condition number: For the spectral condition number of a symmetric and positive definite matrix, it is
«2 (A) = ^ = O {h-2) .
2.5 An Efficient Solver for the Dirichlet Problem in the Rectangle
35
Since the dimension of the matrix is large and the matrix is sparse, iterative solvers are an appropriate approach for solving the linear system of equations. The main costs for iterative solvers are the matrix-vector multiplications (often one per iteration). The cost of one matrix-vector multiplication is for sparse matrices proportional to the number of unknowns. Hence, an optimal solver with respect to the number of floating point operations is given if the number of operations for solving the linear system of equations is proportional to the number of unknowns. It is known that the number of iterations of many iterative solvers depends on the condition number of the matrix:
• (damped) Jacobi method, SOR, SSOR. The number of iteration is proportional to k2(A). That means, if the grid is refined once, h —> h/2, then the number of unknowns is increased by around the factor 4 in two dimensions and also the number of iterations increases by a factor of around 4. Altogether, for one refinement step, the total costs increase by a factor of around 16.
• (preconditioned) conjugate gradient (PCG) method. The number of iterations is proportional to y^K^A), see the corresponding theorem from the class Numerical Mathematics II. Hence, the total costs increase by a factor of around 8 if the grid is refined once.
• multigrid methods. For multigrid methods, the number of iterations on each grid is bounded by a constant that is independent of the grid. Hence, the total costs are proportional to the number of unknowns and these methods are optimal. However, the implementation of multigrid methods is involved.
Remark 2.37. An eigenvalue problem. The derivation of an alternative direct solver is based on the eigenvalues and eigenvectors of the discrete Laplacian. It is possible to computed these quantities only in special situations, e.g., if the Poisson problem with Dirichlet boundary conditions is considered, the domain is rectangular, and the Laplacian is approximated with the five point stencil.
Consider the following eigenvalue problem
Denote the node x = (xl, yj) by x%] and grid functions in a similar way. The solution of this problem is sought in product form (separation of variables)
□
Av(x) = \v(x), x G ujh, v{x) = 0,       x £ 7fe.
V,
(fc) _    <M,x (ky),y
It is
36
2 Finite Difference Methods for Elliptic Equations
where i = 0,. .., nx, j = 0,.. ., ny refers to the nodes and kx = 1,... , nx — 1, ky = l,...,ny — 1 refers to the eigenvalues. Note that the number of eigenvalues is equal to the number of inner nodes, i.e., it is [nx — l)(ny — 1). In this ansatz, also a splitting of the eigenvalues in a contribution from the x coordinate and a contribution from the y coordinate is included. From the boundary condition, it follows that
(kx),x _    (k*),x _    <ky),y _    lky),y _ n
Dividing by v^^'xv^y^'y arKJ rearranging terms, the eigenvalue problem can be split
AA        , x W =    avvj _ vl v)v>
with Afc = AJ^ + X^K Both sides of this equation have to be constant since one of them depends only on i, i.e., on x, and the other one only on j, i.e., on y. The splitting of Afc can be chosen such that the constant is zero. Then, one gets
Axv^x + A<f>f^* = 0,   Ayvf*hy + \^vf*hy = 0.
The solution of these eigenvalue problems is known (exercise)
(k*),x       /~2~ .   ( kxiri\        (x)      4   . 2(kxir V lx      V nx J'      k*      h2x \2nx
v<M.y = 112 sin (bill],  A(»> = ±sin2 ( V[
(2.21)
with i = 0,.. .,nx, j = 0,. .. ,ny and kx = 1,. .., nx — 1, ky = 1,.. .,ny — 1. Using a Taylor series expansion, one obtains now the asymptotic behavior of the eigenvalues as given in (2.20). Note that because of the splitting of the eigenvalues into the directional contributions, the number of individual terms for computing the eigenvalues is only proportional to [nx + ny). □
Remark 2.38. On the eigenvectors, weighted Euclidean inner product. Since the matrix corresponding to A is symmetric, the eigenvectors are orthogonal with respect to the Euclidean vector product. They become orthonormal with
>' ly       \ lly J ny It follows that the solution of the full eigenvalue problem is
2.5 An Efficient Solver for the Dirichlet Problem in the Rectangle 37
respect to the weighted Euclidean vector product
{u,v) = hxhy    ^2   u(x)v(x) = hxhy'^2'^2uijVij, (2.22)
x£ujhU~/h i=0 j=0
with
h - — h - ^-
Tix fly
i.e., then it is
This property can be checked by using the relation
E. 2 /iir\ n SmUJ = 2' n>1-
i=0 v '
The norm induced by the weighted Euclidean vector product is given by
(\ 1/2 M'.. VW )     . (2.23)
The weights are such that this norm can be bounded for constants independently of the mesh, i.e.,
1/2
||l||h = {hxhy(nx + 1)(% + I))'7' = (ljynx + lny + 1)     < 2 (lxly)1/2 .
(2.24)
□
Remark 2.39. Solver based on the eigenvalues and eigenvectors. One uses the ansatz
^(aO = X>«(fc)(aO (2-25) fc
with the Fourier coefficients fk = </,««> = ^EE/.sin sin (**SL) ,    fc = (kx,ky),
\l"Trill 11 11 ~'
x y i=o j=o
with /jj = f(xij). The solution u(a;) of (2.19) is sought as a linear combination of the eigenfunctions
u{x) = ^2 UkV<k^ (x)
38 2 Finite Difference Methods for Elliptic Equations
with unknown coefficients u^. With this ansatz, one obtains for the finite difference operator
Au = ^Mfcyli)(fc) = y~\feAfet>(fe). fc fc
Since the eigenfunctions form a basis of the space of the grid functions, a comparison of the coefficients with the right-hand side (2.25) gives
or, for each component, using (2.21),
= _ E h(k) = _ 2^ £ sin (^) sin (^) ,
i = 0,.. . ,nx, j = 0,. .. ,nv.
It is possible to implement this approach with the Fast Fourier Transform (FFT) with
O (nxny\og2nx +nxny\og2ny) = 0(N\og2N) , N = (nx — l)(nv - 1),
operations. Hence, this method is almost, up to a logarithmic factor, optimal.
□
2.6 A Higher Order Discretizations
Remark 2.4-0. Contents. The five point stencil is a second order discretization of the Laplacian. In this section, a discretization of higher order will be studied. In these studies, only the case of a rectangular domain Q = (0, lx) x (0, ly) and Dirichlet boundary conditions will be considered. □
Remark 2.4-1- Derivation of a fourth order approximation. Let u{x) be the solution of the Poisson equation (2.1) and assume that u(x) is sufficiently smooth. It is
d2u
Luix) = Auix) = Lxu(x) + Lyu(x),   Lau := —-.
dx2a
Let the five point stencil be represented by the following operator
Au = Axu + Ayu. Applying a Taylor series expansion, one finds that
2.6 A Higher Order Discretizations
39
(hj)
Fig. 2.7 The nine point stencil.
Au-Au^^Llu + ^Llu + O^). (2.26)
From the equation — Lu = /, it follows with differentiation that
Lxu = —Lxf — LxLyu,   Lyu = —Lyf — LyLxu. Inserting these expressions in (2.26) gives
An - An = ~^Lxf - ^Lyf - ^-^LxLyU + O (h4) . (2.27) The operator LxLy = dx%dyi can be approximated as follows
LxLyU ^ AxAyU 'Ojixyy
The difference operator in this approximation requires nine points, see Figure 2.7,
AxAyu = —jj^ (ui+1j+1 - 2uij+1 + - 2ui+1j + Auij
Therefore it is called nine point stencil.
One checks, as usual by using a Taylor series expansion, that this approximation is of second order
LxLyu — AxAyu = O (h2) ■
Inserting this expansion in (2.27) and using the partial differential equation shows that the difference equation
h^ÄxAy)u=   f+hlLxf + % 12 v \       12 12
40 2 Finite Difference Methods for Elliptic Equations
is a fourth order approximation of the differential equation (2.1). In addition, one can replace the derivatives of f(x) also by finite differences
Lxf = Axf + 0 (h2x) ,   Lyf =Ayf + 0 (hi) .
Finally, one obtains a finite difference equation —A'u = <j> with
A' = Ax+Ay + ^-^AxAy,   ^ = f + !^Axf + ^Ayf.
Remark 2.4-2. On the convergence of the fourth order approximation. The finite difference problem with the higher order approximation property can be written with the help of the second order differences. Since the convergence proof is based on the five point stencil, the following lemma considers this stencil. It will be proved that one can estimate the values of the grid function by the second order differences. This result will be used in the convergence proof for the fourth order approximation. □
Lemma 2.43. Stability estimate. Let
<^h = {(ihx,jhy) : i = 1,. ..,nx- l,j = 1,. .. ,ny - 1},
and let y be a grid function on cu^ U jh with y(x) = 0 for x € j^. Then, the following estimate holds
\\y\\i^h^h)<M\\Ay\\h,
max{/2 I2 }
with the mesh-independent constant M = " 7 ^ ^S ^e ma^x obtained
by using the five point stencil A = Ax + Ay for approximating the second derivatives, and the norm on the right-hand side is defined in (2.23). Proof. Let {vfj}, h = (kx, ky), be the orthonormal basis with
/ kx TTi \        f kv 7tj - sin | - sin 1 -
f]"x\y
which was derived in Remark 2.37. Then, there is a unique representation of the grid function y = Y2k ykyk and it holds with (2.22)
Ay = J2 y^kVk,    \\Ayfh = J2vlXl (2-28)
fe fe
It follows for x £ ujh, because of |sin(x)| < 1 for all x 6 M, that
\y(x)\=   ^2vkVk(x)   < J2\Vk\ \'"k(v)\ < —,= J2\Vk\-
k k V'x'y k
Using this estimate, applying the Cauchy—Schwarz inequality for sums, and utilizing (2.28)
2.6 A Higher Order Discretizations
41
°'8.0    0.2    0.4    0.6    0.8    1.0    1.2 1.4
Fig. 2.8 The function sin(0)/0.
Now, one has to estimate the last sum. It is already known that 4      2 (kx'ft\  ,4   . 2 (^y77
Afe = —r sin     —— J + — sin -—   ,    kx = 1, . . . , nx - 1, ky = 1, . . . , ny - 1.
fr£        \ 2nx I     M        \ 2nv
k27r2 / sin 4>x \ 2     fcy7T2 f sin tfty
ll    \   4>x   ) l\    \ 4>y
)2>4(f + |)>,4^+^
In performing this estimate, it was used that the function sin (</>)/</> is monotonically decreasing on (0,7r/2), see Figure 2.8, and that
sin (b     sin("7r/2)      2 , ,
(ft 7T/2 7T
The estimate will be continued by constructing a function that majorizes (k2 + fc2) 2 and that can be easily integrated. Let G = {(x,y) : x > 0, y > 0, x2 + y2 > 1} be the first quadrant of the complex plane without the part that belongs to the unit circle, see Figure 2.9. The function {k2 + k2^ has its smallest value in the square [kx — 1, kx] X [ky — 1, ky] in the point (kx, ky). Using the lower estimate of       one obtains
42
2 Finite Difference Methods for Elliptic Equations
A — 16
fe,fe^(i,i)   k fe,fe^(i,i)
E ^  ^  ^ E (^+fer2
fe,fe^(l,l)
smallest value in square v = 1
L j2 /   /   (fc? + feS)"2
16 fc.fc^tl,!)7"--17^-1
;4 _ rfc^ rfcy
16 E /   /   (*2 + ^:
16 k,k*(i,i)Jk*-lJk«-1
I4    f 2
— J   (x2 + y2) cfcccfa/
rd.      p r 16 h Jo
2^   2 dxdy
polar coord.   I      I        I p
— 11       — dd>dp =
p4 16 2 \      2    p=1 J 64
For performing this computation, one has to exclude p —> 0.
4	.   2 (   n \		7T	- 4
~ hi	sin - + \2nxJ		2ny	
_ 7T2	f 2lx \ 2 . 2	V2Z:J	-f	'3k.)2
~ z2	8         7T2 8 ^r2 + T2"^2 -	16 ¥'		
\(-[ t\ = —7- sin   1 - I -f- -rr sm   1 - 1 = -r- sin   1 - I H--7r sin^
ft y 7T ^
(2.30)
For this estimate, the following relations and the monotonicity of sin(:e)/a;, see Figure 2.S, were used
la haT\     7r      / sin <pa \2     (sin(7r/4)\2 8
Collecting all estimates gives
2.6 A Higher Order Discretizations 43
1   _ .-2 1        [4      ^4 [4
2^AT-A(i,i)+    2_v Al-256+"64-16-
k      k fe,fe^(l,l) fe
Inserting this estimate in (2.29), the final bound has the form
ll»lll»(Uhu7h) < -J= WML J =: M \\Ay\\h .
Theorem 2.44. Convergence of the higher order finite difference scheme. Let Q = (0,/^) x (0,ly). The finite difference scheme
—A'u(x) = (f>(x), x G uih, u(x) = g(x), x G 7fe,
with
A' = Ax + Ay + ^^AxAy,   ^ = f + ^LAxf+'^Ayf, converges of fourth order.
Proof. Analogously as in the proof of Theorem 2.33, one finds that the following equation holds for the error e = u(xi,yj) — Uij:
-A'e{x) = ip(x), ip = O (h4) ,x e ojh, e(x) = 0,       x e lh-
Let Q}-L be the vector space of grid functions, which are non-zero only in the interior, i.e., at the nodes from uj^, and which vanish on 7^. Let Aay = —Aay, y 6 Qh, a G {x,y}. The operators Aa : Qh     Qh are linear and they have the following properties:
• They are symmetric and positive definite, i.e., Aa = j4* > 0, where .A* is the adjoint (transposed) of Aa, and (Aau,v) = (u,Aav), Vw,uG Qh-
• They are elliptic, i.e., (Aau,u) > \^*\u, u), \/u 6 J?^, with
4     ■   2 (Kha\ 8
At    = —- sin - >
h%        \ 2,la j 1%, see (2.30).
• They are bounded, i.e., it holds (Aau:u) < X^_1(u,u) with
^-1 = ST-'(£) * £   ^   (^.«) (»D
and ||j4q;||2 < since the spectral norm of a symmetric positive definite matrix
equals the largest eigenvalue.
• They are commutative, i.e., AxAy = AyAx.
• It holds AxAy = (AxAy)*.
The error equation on uj^ is given by
h?
Axe + Aye — (kx + k,y)AxAye = A'e = tp   with   ka = -y^-. (2.32)
44
2 Finite Difference Methods for Elliptic Equations
Using the boundedness of the operators, one finds with (2.31) for all v 6 that
(ltXAXAyV + KyAXAyV,v)    =   ((K-XAX)AyV,v)  + ((KyAy)AXV,v)
12 feg
Now, it follows for all v 6 i?^ that
(yl'v,     = ((Ax + Ay) v, v) — (itxAxAyv + kyAxAyv, v) 2
> ^ ((^ + Ay)v,v) > 0.
The matrices on both sides of this inequality are symmetric and because the matrix on the lower estimate is positive definite, also the matrix at the upper estimate is positive definite. The matrices commute since the order of applying the finite differences in x and y direction does not matter. Using these properties, one gets (exercise)
'■(Ax + Ay)e
<\\A'e\\h = U\\h,
(2.33)
where the last equality follows from (2.32). The application of Lemma 2.43 to the error, (2.33), (2.32), and (2.24) yields
Hell,«
I2 31 < -=\\(^+Ay)e\\h <
4 -\J Ix^y 4-^/ Ix^y
2 \J Ix^y
' = (hxhy(nx + l)(ny + l))1/2 ||^||IOO({JhU7h)
xiy
W\\h
al2 ( nx + 1 ny + 1 4   \   nx ny
1/2
IIV-II
I°°(ufcU7fc)
Remark 2.4-5. On the discrete maximum principle. Reformulation of the finite difference scheme —A'u = <j> in the form studied for the discrete maximum principle gives for the node
a(x)u(x) = ^2 b(xyy)u(y) + <t>{x),
2
a(x) b(x,y)
b(x,y) b(x,y)
hi + hi   12    +hy> hihy- 3 [hi + hy)> u'
1      1 /, 2    , 2\    2       1 I 5      1 \    . , .
K + /l»)l2l2 =fi     12 -12     ' l±1'3,
h% 12
1/1 5
hlhl     6 \h2x hy
6
hi hi
1  ( 1
1
\2\hl h\
(left, right node) i,j±l, (bottom, top node)
i ± 1, j ± 1, (other neighbors)
2.7 Summary 45
Hence, the assumptions for the discrete maximum principle, see Remark 2.15, are satisfied only if
1        hx r-
V5 «y
Consequently, the ratio of the grid widths has to be bounded and it has to be of order one. In this case, one speaks of an isotropic grid. □
2.7 Summary
Remark 2.4-6. Summary.
• Finite difference methods are the simplest approach for discretizing partial differential equations. The derivatives are just approximated by difference quotients.
• They are very popular in the engineering community.
• One large drawback are the difficulties in approximating domains that are not of tensor-product type. However, in the engineering communities, a number of strategies have been developed to deal with this issue in practice.
• Another drawback arises from the point of view of numerical analysis. The numerical analysis of finite difference methods is mainly based on Taylor series expansions. For this tool to be applicable, one has to assume a high regularity of the solution. These assumptions are generally not realistic.
• In Numerical Mathematics, one considers often other schemes then finite difference methods. However, there are problems, where finite difference methods can compete with other discretizations.
□
Chapter 3
Introduction to Sobolev Spaces
Remark 3.1. Contents. Sobolev spaces are the basis of the theory of weak or variational forms of partial differential equations. A very popular approach for discretizing partial differential equations, the finite element method, is based on variational forms. In this chapter, a short introduction into Sobolev spaces will be given. Recommended literature are the books Adams (1975); Adams & Fournier (2003), and Evans (2010). □
3.1 Elementary Inequalities
Lemma 3.2. Inequality for strictly monotonically increasing function. Let f : M+ U {0} —> K be a continuous and strictly monotonically increasing function with f(0) = 0 and f(x) —> oo for x —> oo. Then, for all a,teI+U {0}, it is
ab< [ f(x) dx+ [ f-\y) dy, Jo Jo
where f~1{y) is the inverse of f(x).
Proof. Since f(x) is strictly monotonically increasing, the inverse function exists. The proof is based on a geometric argument, see Figure 3.1.
Consider the interval (0, a) on the x-axis and the interval (0, b) on the y-axis. Then, the area of the corresponding rectangle is given by ab, f^ f(x) dx is the area below the curve, and Jq f~1(y) dy is the area between the positive y-axis and the curve. From Figure 3.1, the inequality follows immediately. The equal sign holds only iff f(a) = b. M
Remark 3.3. Young's1 inequality. Young's inequality
ab < |a2 + yb2   V a, b £ R+ U {0},£ £ R+, (3.1) 1 William Henry Young (1863 - 1942)
47
48
3 Introduction to Sobolev Spaces
follows from Lemma 3.2 with f(x) = ex, f^1{y) = £~1y- It is also possible to derive this inequality from the binomial theorem. For proving the generalized Young inequality
ab<^-ap + -^bq,   Va,6eK+U{0},£GK+, (3.2)
withp_1+g_1 = l,p,q € (1, oo), one chooses f(x) = xv~x, f^1{y) = y1/(p~1'1 and applies Lemma 3.2 with intervals where the upper bounds are given by ea and e~1b. □
Remark 3.4- Cauchy-Schwarz inequality.
• The Cauchy2-Schwarz3 inequality (for vectors, for sums)
|fej/)| < Uhkh V^EeR"' (3-3)
where (•, •) is the Euclidean product and ||-||2 the Euclidean norm, is well known.
• One can prove this inequality with the help of Young's inequality. First, it is clear that the Cauchy-Schwarz inequality is correct if one of the vectors is the zero vector. Now, let x,y with ||a;||2 = \\y\\2 = 1- One obtains with the triangle inequality and Young's inequality (3.1)
| (£>2/) | =
n 1    n 1 n
2 ^ 1 " 2
Hence, the Cauchy-Schwarz inequality is correct for x, y. Last, one considers arbitrary vectors x 0, y 0. Now, one can utilize the homogeneity of the Cauchy-Schwarz inequality. From the validity of the Cauchy-Schwarz inequality for x and y, one obtains by a scaling argument
2 Augustin Louis Cauchy (1789 - 1857)
3 Hermann Amandus Schwarz (1843 — 1921)
3.1 Elementary Inequalities 49
Kllill^Vi' II1II2^1
£ y
Both vectors x,y have the Euclidean norm 1, hence
n-ii Vii Ife^l^1 \&y)\<m\2\\y\\2-
\\x\\2\\y\\2> -' "-112
• The generalized Cauchy-Schwarz inequality or Holder's4 inequality
V=i       /      \i=i J
with p^1 + q^1 = q € (l,oo), can be proved in the same way with the help of the generalized Young inequality.
□
Definition 3.5. Lebesgue spaces. The space of functions that are Lebesgue5 integrable on Q to the power of p € [1, oo) is denoted by
which is equipped with the norm
i/p
ll/lli>(fi)= (J \f(x)\'dx
For p = oo, this space is given by
= {/ : |/(a;)| < oo almost everywhere in Q}
with the norm
=ess supxea\f(x)\.
□
Lemma 3.6. Holder's inequality. Let p^1 +      = l,p,q € [l,oo]. If u G
LP(Q) and v € Lq(Q), then it is uv € L1(f2) and it holds that
IIHIli(r) < \\u\\Lp(n) IMIw(r) ■ (3-4) If P = 9 = 2, then this inequality is also known as Cauchy-Schwarz inequality
4 Otto Holder (1859 - 1937)
5 Henri Lebesgue (1875 - 1941)
50
3 Introduction to Sobolev Spaces
IIHIli(r) < MIl^q) ■ (3-5)
Proof, i) p,q 6 (1, oo). First, one has to show that -ut;(a3) can be estimated from above by an integrable function. Setting in the generalized Young inequality (3.2) e = 1, a = |u(a3)|, and b = -u(£c) gives
\u(x)v(x)\ < - \u(x)\p + - \v(x)\q . (3.6)
v q
Since the right-hand side of this inequality is integrable, by assumption, it follows that uv 6 L1(i7). Integrating (3.6), Holder's inequality is proved for the case = IMIw(fi) = 1
/  \u(x)v(x)   dx < — I  \u(x)\p dx H--/   -u(ic) |g dx = 1.
Jq pJq qJo
The general inequality follows, for the case that both functions do not vanish almost everywhere, with the same homogeneity argument as used for proving the Cauchy—Schwarz inequality of sums. In the case that one of the functions vanishes almost everywhere, (3.4) is trivially satisfied.
ii) p = 1, q = oo. It is
j \u(x)v(x)\ dx < j |u(a;)|ess suPa!efi\v(x) dx = |M|Li(fi) IMILco(fi) .
3.2 Weak Derivative and Distributions
Remark 3.7. Contents. This section introduces a generalization of the derivative which is needed for the definition of weak or variational problems. For an introduction to the topic of this section, e.g., see Haroske & Triebel (2008)
Let Q C Rd be a domain with boundary _T = dO, d G N, Q 7^ 0. A domain is always an open set. □
Definition 3.8. The space C^°(J?). The space of infinitely often differen-tiable real functions with compact (closed and bounded) support in Q is denoted by C0°°(i7)
C*0°°(J?) = {v : v £ C°°(0), supp(v) C O},
where
supp(ij) = {x G Q : v(x) 7^ 0}.
In particular, functions from Cq°(J?) vanish in a neighborhood of the boundary. □
Definition 3.9. Convergence in Cg°(J?). The sequence {(f>n(x)}^=1, (f>n G Cg°(J?) for all n, is said to convergence to the zero functions if and only if a) 3K C Q,K compact with supp(</>„) C K for all n,
3.2 Weak Derivative and Distributions
51
b) Da(f>n(x) —> 0 for n —> oo on K for all multi-indices a = (a.\,.. ., ad),
|a| = Qi + ... + ad. It is
lim 4>n = 4> lim (4>n - 4>) = 0.
□
Definition 3.10. Weak derivative. Let f,F G L\oc(Q). A function m belongs to L\oc(Q) if for each compact subset Q' C Q, it holds
/   |u(a;)| dx < oo.
If for all functions g G C^°(J?), it holds that
f F(x)g(x) dx = (-l)^ [ f(x)Dag(x) dx,
then F(x) is called weak derivative of f(x) with respect to the multi-index a. □
Remark 3.11. On the weak derivative.
• The notion 'weak' is used in mathematics if something holds for all appropriate test elements (test functions).
• One uses the same notations for the derivative as in the classical case : F(x)=Daf(x).
• If f{x) is classically differentiate on Q, then the classical derivative is also the weak derivative.
• The assumptions on f(x) and F(x) are such that the integrals in the definition of the weak derivative are well defined. In particular, since the test functions vanish in a neighborhood of the boundary, the behavior of f(x) and F(x) if x approaches the boundary is not of importance.
• The main aspect of the weak derivative is due to the fact that the (Lebesgue) integral is not influenced from the values of the functions on a set of (Lebesgue) measure zero. Hence, the weak derivative is defined only up to a set of measure zero. It follows that f(x) might be not classically differentiate on a set of measure zero, e.g., in a point, but it can still be weakly differentiate.
• The weak derivative is uniquely determined, in the sense described above.
□
Example 3.12. Weak derivative. The weak derivative of the function f(x) = \x\ is
C -1 x < 0, f'(x) = \  0 x = 0, 1 x > 0.
52 3 Introduction to Sobolev Spaces
In x = 0, one can use also any other real number. The proof of this statement follows directly from the definition and it is left as an exercise. □
Definition 3.13. Distribution. A continuous linear functional defined on Cg°(J?) is called distribution. The set of all distributions is denoted by (C0°°(rt))'.
Let u G Co°(n) and if> G (C^°(J?)) , then the following notations are used for the application of the distribution to the function
ip(u) = (ip,u) G R.
Remark 3.14- On distributions. Distributions are a generalization of functions. They assign each function from C^°(0) a real number. □
Example 3.15. Regular distribution. Let u G L\oc(Q). Then, a distribution is defined by
[ u(x)4>(x) dx = {i>,4>)   V <j> G C"0°°(J?). Jn
This distribution will be identified with u G L\oc(Q).
Distributions with such an integral representation are called regular, otherwise they are called singular. □
Example 3.16. Dirac distribution. Let £ G Q be fixed, then
defines a singular distribution, the so-called Dirac6 distribution or 5-distribu-tion. It is denoted by 5^ = S(x — £). □
Definition 3.17. Derivatives of distributions. Let <j) G (C^°(J?))' be a distribution. The distribution ip G (C^°(J?))' is called derivative in the sense of distributions or distributional derivative of <f> if
^,«) = (-l)H^)JDtt«) V^eCHfl),
a = (qi, ... ,ad), QLj > 0, j = 1,... ,d, \a\ = qi + ... + ad. □
Remark 3.18. On derivatives of distributions.
• Each distribution has derivatives in the sense of distributions of arbitrary order.
• If the derivative in the sense of distributions Dau(x) with u G L\oc(Q) is a regular distribution, then also the weak derivative of u{x) exists and both derivatives are identified.
□
6 Paul Adrien Maurice Dirac (1902 - 1984)
3.3 Lebesgue Spaces and Sobolev Spaces
53
3.3 Lebesgue Spaces and Sobolev Spaces
Remark 3.19. On the spaces LP{Q). These spaces were introduced in Definition 3.5.
• The elements of LP(Q) are, strictly speaking, equivalence classes of functions that are different only on a set of Lebesgue measure zero.
• The spaces LP(Q) are Banach7 spaces (complete normed spaces). A space X is complete, if each so-called Cauchy sequence {un}^=0 G X, i.e., for all e > 0 there is an index no(e) such that for all i,j > rao(e)
IK -uj\\x <£'
converges and the limit is an element of X.
• The space L2(Q) becomes a Hilbert8 spaces with the inner product
(f,g)= f f(x)g(x)dx,    ||/|U2 = (/,/)1/2, f,g€L2{Q). J n
• The dual space of a space X is the space of all bounded linear functionals defined on X. Let Q be a domain with sufficiently smooth boundary _T and consider the Lebesgue space LP(Q), p G [1, oo], then
(Lp(Q))' = Lq(Q) with p,ge(l,oo), 1 + 1 = 1,
p q
{lW)'= l°°(n), {l°°{n))'? l\n),
where the prime symbolizes the dual space. The spaces L1(i7), L°°(Q) are not reflexive, i.e., the dual space of the dual space is not the original space again.
□
Definition 3.20. Sobolev9 spaces. Let k G N U {0} and p G [l,oo], then the Sobolev space Wk'p{Q) is defined by
Wk'p{Q) := {u G Lp(fi) : Dau G Lp(n) V a with \a\ < k}.
This space is equipped with the norm
\\u\\w".p(n) :=        Wdau\\lp(n)- (3-7)
\a\<k
□
7 Stefan Banach (1892 - 1945)
8 David Hilbert (1862 - 1943)
9 Sergei Lvovich Sobolev (1908 - 1989)
54
3 Introduction to Sobolev Spaces
Remark 3.21. On the spaces Wk'p{0).
• Definition 3.20 has the following meaning. From u G LP(Q), p G [l,oo), it follows in particular that u G Lloc(Q), such that u defines (represents) a distribution. Then, all derivatives Dau exist in the sense of distributions. The statement Dau G LP(Q) means that the distribution Dau G (Cg°(Q))' can be represented by a function from LP(Q).
• One can add elements from Wk'p(Q) and one can multiply them with real numbers. The result is again a function from Wk'p{Q). With this property, the space Wk'p(Q) becomes a vector space (linear space). It is straightforward to check that (3.7) is a norm, [exercise)
• It is Dau(x) = u(x) for a = (0,..., 0) and W°'P{Q) = LP{Q).
• The spaces Wk'p(Q) are Banach spaces.
• Sobolev spaces have for p G [l,oo) a countable basis {<pn(x)}^=1 (Schau-der10 basis), i.e., each element u(x) can be written in the form
• Sobolev spaces are reflexive for p G (1, oo).
• The subspace C°°{Q) n Wk'p{Q) is dense in Wk'p{Q)1 see (Gilbarg & Trudinger, 1983, p. 154). Under a certain condition on the smoothness of the boundary of a bounded domain Q, one can show that Cg°(J?) is dense in Wk'p(Q), p G [1, oo), with respect to the norm (3.7), e.g., (Adams, 1975, Thm. 3.18). With this property, one can characterize the Sobolev spaces Wk'p(Q) as completion of the functions from Cg°(J?) with respect to the norm (3.7). It follows that Ck(H) is dense in Wk'p(Q), p G [l,oo).
• The Sobolev space Hk(Q) = Wk'2(Q) is a Hilbert space with the inner product
Definition 3.22. The space Wq,p{Q). The Sobolev space Wq,p{Q) is defined as the completion of Cg°(J?) in the norm of Wk'p(Q)
and the induced norm = {u,u)
1/2
□
Wk'p(Q) = C0°°(ß)
□
10 Juliusz Pawel Schauder (1899 - 1943)
3.4 The Trace of a Function from a Sobolev Space
55
3.4 The Trace of a Function from a Sobolev Space
Remark 3.23. Motivation. This class considers boundary value problems for partial differential equations. In the theory of weak or variational solutions, the solution of the partial differential equation is searched in an appropriate Sobolev space. Then, for the boundary value problem, this solution has to satisfy the boundary condition. However, since the boundary of a domain is a manifold of dimension [d — 1), and consequently it has Lebesgue measure zero, one has to clarify how a function from a Sobolev space is defined on this manifold. This definition will be presented in this section. □
Definition 3.24. Lipschitz boundary, Lipschitz domain, (Grisvard, 1985, Def. 1.2.1.1). Let Q be a bounded domain in Rd, then Q is called Lipschitz11 domain, respectively the boundary _T of Q is called Lipschitz boundary, if for every x g _T there exists a neighborhood U of x in Md and new orthogonal coordinates (j/j,. .., y^) such that
1) U is a hypercube in the new coordinates
U = {(yi, ■ ■ ■ ,yd) ■ -a-i < Vi < a,i, i = 1,... ,d} .
2) There exists a Lipschitz continuous function <f>, defined in
U'= {(yi,. ■ ■ ,yd-i) ■ -a-i < Vi < au i = 1,... ,d - 1} , such that
\4>{y')\ < y for eveir v' = (yi, - •• ,yd-i) e u',
ODU = {y = {y',yn)€V : yn < 4>{y')} , rnU = {y = (y',yn)€V : yn = cb(y')}.
□
Remark 3.25. Lipschitz boundary.
• In a neighborhood of y, Q is below the graph of <f> und the boundary _T is the graph of <f>.
• The domain Q is not on both sides of the boundary at any point of r.
• The outer normal vector is defined almost everywhere at the boundary and it is almost everywhere continuous.
□
Example 3.26. On Lipschitz domains.
• Domains with Lipschitz boundary are, for example, balls or polygonal domains in two dimensions where the domain is always on one side of the boundary.
11 Rudolf Otto Sigismund Lipschitz (1832 - 1903)
56
3 Introduction to Sobolev Spaces
Fig. 3.2 Polyhedral domain in three dimensions that is not Lipschitz continuous (at the corner where the arrow points to).
• A domain that is not a Lipschitz domain is a circle with a slit
At the slit, the domain is on both sides of the boundary. • In three dimension, a polyhedral domain is not not necessarily a Lipschitz domain. For instance, if the domain is build of two bricks that are laying on each other like in Figure 3.2, then the boundary is not Lipschitz continuous where the edge of one brick meets the edge of the other brick.
Theorem 3.27. Trace theorem. Let Q C Md, d > 2, with a Lipschitz boundary. Then, there is exactly one linear and continuous operator 7 : WAl'p(i7) Lp(r), p G [1, 00), that gives for functions u G C{Q) l~l WAl'p(i7) the classical boundary values
i.e., ju(x) = u(x)\xer-
Proof. The proof can be found in the literature, e.g., in Adams (1975); Adams & Fournier
Remark 3.28. On the trace.
• The operator 7 is called trace or trace operator.
• By definition of the trace, one gets for u G C(Q) the classical boundary values. By the density of C°°{0) C C(JTj in W1'p{0) for domains with smooth boundary that for all u G W1'P(Q) there is a sequence {un}^=1 G C°°(Q) with un —v u in W1'P(Q). Then, the trace of u is defined to be ju = limn^00(7M„).
• Since a linear and continuous operator is bounded, there is a constant C > 0 with
0 = {(x,y) : x2 +y2 <l}\{(x,y) : x>0,y = 0}.
□
ju{x) = u(x), x€T, Vue C{f2) n W1'p{ú),
(2003).
\ju\\
Lr(r)
<C\\u\\
W1-r(n)
Vue W1'p{ú)
3.5 Sobolev Spaces with Non-Integer and Negative Exponents
57
IM
c(w1-p(n),LP(r))
< C.
It is
jDau(x) = 0   V« G Wo'p{ü),\a\ < k-1.
(3.8)
3.5 Sobolev Spaces with Non-Integer and Negative Exponents
Remark 3.29. Motivation. Sobolev spaces with non-integer and negative exponents are important in the theory of variational formulations of partial differential equations.
Let Q C Rd be a domain and p G (1, oo) with p^1 + q^1 = 1. □
Definition 3.30. The space W~k'q{n). The space W~k'q{n), k G NU{0}, contains distributions that are defined on Wk'p{Q)
W-k'q(0) = {^€(Cr(0))' : IHw-*.,(fi) < oo}
with
li   li (v'u)
Remark 3.31. On the spaces W-k'p(Q).
• It is W-k'q{0) = [wofc'p(«)]', i.e., W-k'q{0) can be identified with the
dual space of W^'P{Q). In particular, it is H~1{Q) = (H%{Q)) ■
• It is
... c W2'p(0) c W1'P{Q) c Lp(0) c W'^iO) c W~2'q{Q)...
□
Definition 3.32. Sobolev—Slobodeckij space. Let s G M, then the Sobo-lev-Slobodeckij12 or Sobolev space HS(Q) is defined as follows:
• s G Z. ifs(J?) = W"'2(J?).
• s > 0 with s = fc + a, fc G N U {0}, a G (0,1). The space Hs(0) contains all functions u for which the following norm is finite:
L. N. Slobodeckij
58
3 Introduction to Sobolev Spaces
with
0, v)Hs(n) = (u, v)Hk + (u, v)k+a,    \u\2k+a = (u, u)k+a
and
(u,v)k+a = ^
(Dau(x) - Dau(y)) (Dav(x) - Dav(y))
II \\d+2a \\X-V\\2
dxdy,
s < 0. ifs(J?) = [if^J?)]' with Hos{H) = C^°(«)"'"H °<!7)
3.6 Theorem on Equivalent Norms
Definition 3.33. Equivalent norms. Two norms 11-11-j^ and |-||2 on the linear space X are said to be equivalent if there are constants C\ and C2 such that
Ci IMIj < |H|2 < C2 IMIj Vtiel.
□
Remark 3.34- On equivalent norms.
• Many important properties, like continuity or convergence, do not change if an equivalent norm is considered.
• In finite-dimensional spaces, all norms are equivalent.
□
Theorem 3.35. Equivalent norms in Wk'p{Q) (Smirnow, 1967, § 114, Satz 3). Let Q c Md be a bounded domain with Lipschitz boundary _T, p g [l,oo], and fceN. Let {fi}\=i be a system of Junctionals with the following properties:
1) fi : Wk'p(Q) ->• K+ U {0} is a seminorm,
2) boundedness: 3Cl > 0 with 0 < < Cl \\v\\wk,p(n), V v g Wk'p{Q),
3) fi is a norm on the polynomials of degree k — 1, i.e., if for v g Pk-i = |S|o;|<fc-i Ca^" |> it holds that fi(v) = 0, i = 1,. .., I, then it is v = 0.
Then, the norm \\'\\wk<p(n) defined in (3.7) and the norm
3.6 Theorem on Equivalent Norms
59
\\Wk-r(n)
'\Wk-r(n)
1/p
with
1/p
= E / i^wr ^
are equivalent.
Remark 3.36. On seminorms. For a seminorm /i(-), one cannot conclude from fi{v) = 0 that v = 0. The third assumptions however states, that this conclusion can be drawn for all polynomials up to a certain degree. □
Example 3.37. Equivalent norms in Sobolev spaces. • The following norms are equivalent to the standard norm (3.7) in W1'P(Q):
u dx
u ds
a) MlW1-r(Q) ~ ^
b) Mlw1-r(n) = ^
c) \\u\\'w^(n) = (/,IMIP ds + Mw-1-^)) In W^J?), it is
1/p
1/p
1/p
lwfc.P(17)
Sn1
+ \u\wk.p(n)
i/p
equivalent to the standard norm. Here, n denotes the outer normal on _T
with ||n||2 = 1.
In the case Wq'p(0), one does not need the regularity of the boundary. It is
Mlw0fcj(i?) = Mwfc-i=(i?) >
i.e., in the spaces Wq'p{0) the standard seminorm is equivalent to the standard norm.
In particular, it is for u G Hq(Q) (k = l,p = 2)
Ci \\u\\m(a) < IIVmIIl2(r) < C"2 \\u\\HHa). It follows that there is a constant C > 0 such that
\H\mn)<C\\Vn\\LHn)    Vu(=H*(n). (3.9)
□
60
3 Introduction to Sobolev Spaces
3.7 Some Inequalities in Sobolev Spaces
Remark 3.38. Motivation. This section presents a generalization of the last part of Example 3.37. It will be shown that for inequalities of type (3.9), it is not necessary that the trace vanishes on the complete boundary.
Let Q C Md be a bounded domain with Lipschitz boundary A and let i"i C A with measRd-i (Ai) = Jr  ds > 0.
One considers the space
V0 = {v€ WAl'p(i7) : v\n = 0} C WAl'p(i7) if A C A,
v0 = w01,p(fl) if A =r,
with p£ [1, oo). □
Lemma 3.39. Friedrichs13 inequality, Poincare14 inequality, Poincare— Friedrichs inequality. Let p G [l,oo) and measVd-i (Ai) > 0. Then, it is for all u G Vo
/" |u(a;)|p da; < CP [ \\Vu(x)f2 dx, (3.10) where \\-\\2 is the Euclidean vector norm.
Proof. The inequality will be proved with the theorem on equivalent norms, Theorem 3.35. Let : W1'"^) -> m+ u {0} with
h(u) = ^ |«(S)|" ds)
This functional has the following properties:
1) fl(u) is a seminorm.
2) It is bounded, since
0 < fi(u) = i^j   \u(s)\p ds^j /P < i^j \u(s)\p ds^j
= \\u\\LP(r) = \hu\\LP(r) < ^ IMIw1^^) ■
The last estimate follows from the continuity of the trace operator. 3) Let v £ Pq, i.e., v is a constant. Then, one obtains from
0 = h(v) = \v(s)\v ds^j      = \v\ (measg^x (A))1/" ,
that \v\ = 0.
Hence, all assumptions of Theorem 3.35 are satisfied. That means, there are two constants C*i and c2 with
13 Kurt Otto Friedrichs (1901 - 1982)
14 Henri Poincare (1854 - 1912)
3.8 The Gaussian Theorem
61
Cl(yj   \u(s)\p ds + j^\\Vu(x)\\p dx^j      < M\wi,p(Q) < C2 \\u\\'w^p(Q)
1 ]lW1>P(f2)
for all u 6 W1,P{Q). In particular, it follows that
f \u(x)\p dx + f \\Vu(x)\\p dx<c$( [ \u(s)\p ds + f \\Vu(x)\\p dx] Jq Jq \Jri_ Jq J
or, neglecting the non-negative term on the left-hand side,
/ \u(x)\p dx<CP( / \u(s)\p ds+ I \\Vu(x)\\p2 dx) Jq \Jri_ Jq J
with Cp = Cp. Since u G Vq vanishes on Fi, the statement of the lemma is proved. ■
Remark 3.4-0. On the Poincare-Friedrichs inequality. In the space Vq becomes H^/i.p a norm that is equivalent to |HIvi/i,p(j?)- The classical Poincare-Friedrichs inequality is given for _T\ = _T and p = 2
IMIl2(R) < Cp \\Vu\\L2(n) Vm G Hl(Q),
where the constant depends only on the diameter of the domain O, e.g., see (Galdi, 2011, Theorem II. 5.1). □
3.8 The Gaussian Theorem
Remark 3.4-1- Motivation. The Gaussian theorem is the generalization of the integration by parts from calculus. This operation is very important for the theory of weak or variational solutions of partial differential equations. One has to study, under which conditions on the regularity of the domain and of the functions it is well defined. □
Theorem 3.42. Gaussian theorem. Let Q C Md,d > 2, be a bounded domain with Lipschitz boundary r. Then, the following identity holds for all u G W1'1(Q)
I diu{x) dx =     u(s)rii(s) ds, (3-11) Jn J r
where n is the unit outer normal vector on T.
Proof. It is referred to the literature. ■
Corollary 3.43. Vector field. Let the conditions of Theorem 3.42 on the domain Q be satisfied and let u G (WAl'I(i7))   be a vector field. Then, it is
V • u(x) dx = / u(s) ■ n(s) ds. n Jr
62 3 Introduction to Sobolev Spaces
Proof. The statement follows by adding (3.11) from i = 1 to i = d. ■
Corollary 3.44. Integration by parts. Let the conditions of Theorem 3.42 on the domain Q be satisfied. Consider u G WAl'p(i7) and v G WAl'9(i7) with p € (1, oo) and ^ + ^ = 1. Then, it is
/ dlu(x)v(x) dx = / m(s)i)(s)tIi(s) ds — I u(x)dlv(x) dx. Jn Jr Jn
Proof,  exercise. ■
Corollary 3.45. First Green15's formula. Let the conditions of Theorem 3.42 on the domain Q be satisfied, then it is
Vu(x) ■ X7v(x) dx = f ^—(s)v(s) ds — f Au(x)v(x) dx n Jr an Jn
for all u G H2(Q) and v G if1 (J?).
Proof. From the definition of the Sobolev spaces, it follows that the integrals are well denned. Now, the proof follows the proof of Corollary 3.44, where one has to sum over the components and one has to take di'v instead of v. ■
Remark 3.46. On the first Green's formula. The first Green's formula is the formula of integrating by parts once. The boundary integral can be equiva-lently written in the form
J Vu(s) ■ n(s)v(s) ds.
The formula of integrating by parts twice is called second Green's formula.
□
Corollary 3.47. Second Green's formula. Let the conditions of Theorem 3.42 on the domain Q be satisfied, then one has
/ (Au(x)v(x) — Av(x)u(x)) dx = / (—(s)v(s)--(s)u(s) I ds
Jn Jr \dn On J
for all u,v G H2(0).
3.9 Sobolev Imbedding Theorems
Remark 3.48. Motivation. This section studies the question which (Sobolev) spaces are subspaces of other Sobolev spaces. With this property, called
15 Georg Green (1793 - 1841)
3.9 Sobolev Imbedding Theorems
63
imbedding, it is possible to estimate the norm of a function in the subspace by the norm in the larger space, compare (3.12). □
Lemma 3.49. Imbedding of Sobolev spaces with same integration power p and different orders of the derivative. Let Q c Md be a domain, p g [l,oo], and k<m, then it is Wm'p{0) c Wk'p(Q).
Proof. The statement of this lemma follows directly from the definition of Sobolev spaces, see Definition 3.20. B
Lemma 3.50. Imbedding of Sobolev spaces with the same order of the derivative k and different integration powers. Let Q c Md be a
bounded domain, k > 0, and p,q g [1, oo] with q > p. Then, it is Wk'q(Q) c Wk'p{Q).
Proof,  exercise. ■
Remark 3.51. Imbedding of Sobolev spaces with the same order of the derivative k and the same integration power p in imbedded domains. Let Q c Md be a domain with sufficiently smooth boundary E, k > 0, and p g [1, oo]. Then, there is a map E : Wk'p{Q) ->• Wk'p(Rd), the so-called (simple) extension, with
• Ev\n = v,
• l-E^llw^^) - C IMIw"=.p(r)' with C > 0 independent of v,
e.g., see (Adams, 1975, Chapter IV) for details. Likewise, the natural restriction e   :   Wk'p{Rd) ->• Wk'p(Q) can be defined and it is \\ev\\wk,T{n) <
\\v\\wk<P(Rd)- D
Theorem 3.52. A Sobolev inequality. Let Q c Md be a bounded domain with Lipschitz boundary E, k > 0, and p g [1, oo) with
k > d for p = 1, k > d/p    for p > 1.
Then, there is a constant C such that for all u g Wk'p{Q), it follows that u g Cb(^), where
Cb(^) = {v £ C{Q) : v is bounded},
and it is
Wu\\cB(n) = IMIl~(r) < C IMIw^fi) • (3-12)
Proof. See literature, e.g., Adams (1975); Adams & Fournier (2003). ■
Remark 3.53. On the Sobolev inequality.
• The Sobolev inequality states that each function with sufficiently many weak derivatives (the number depends on the dimension of Q and the integration power) can be considered as a continuous and bounded function in Q, i.e., there is such a representative in the equivalence class where this function belongs to. One says that Wk'p(Q) is imbedded in Cb{&)-
64
3 Introduction to Sobolev Spaces
Fig. 3.3 The function f(x) of Example 3.55 for d = 2.
• It is
C(Q) c Cb(O) c C(O).
Consider Q = (0,1) and f\(x) = 1/rr and ^(z) = sin(l/a;). Then, /j G C(i?), /i £ CB(Í2) and h e Cb(í2), /2 Č C(77). • Of course, it is possible to apply this theorem to weak derivatives of functions. Then, one obtains imbeddings like Wk'p{Q) CSB(Q) for (k-s)p > d,p > 1. A comprehensive overview on imbeddings can be found in Adams (1975); Adams & Fournier (2003).
□
Example 3.54- H1(f2) in one dimension. Let d = 1 and Í2 be a bounded interval. Then, each function from if1 (J?) (fc = l,p = 2) is continuous and bounded in Í2. □
Example 3.55. H1(f2) in higher dimensions. The functions from if1(i?) are in general not continuous for d > 2. This property will be shown with the following example.
Let Í2 = {x G
|a;||2 < 1/2} and f(x) = In |ln |a;||21, see Figure 3.3.
For ||a;||2 < 1 /2, it is |ln ||a;||2| = — In ||a;||2 and one gets for x 7^ 0
dp í    \ 1 1 -Ei -Ei
In llxll
Ila;ll2lnlla;ll2
For p < d, one obtains
dx.
	p	1
\\xh		la; L In
		
		
<1		
1^112 111 11^112
The estimate of the second factor can be obtained, e.g., with a discussion
of the curve. Using now spherical coordinates, p ■ sphere, yields
and Sd 1 is the unit
3.9 Sobolev Imbedding Theorems
domain without Lipschitz boundary in (0,0)
0.0       0.2       0.4       0.6       0.8 1.0
Fig. 3.4 Domain of Example 3.56.
\dif{x)\pdx<
dx
n \\x\\d2\\a\\x\\2\d
1/2 pd-!
,(sd
1/2
S*-1 Jo pd\\np\ dP    _ _____iod-i
dpdcu
p\\np\a
;(Sd
dt
< oo,
because of d > 2.
It follows that dif G LP(Q) with p < d. Analogously, one proves that / G Lp(Q) with p < d. Altogether, one has / G W1'p(Q) with p < d. However, it is / G" L°°(Q) and consequently / G" C_a(i7). This example shows that the condition k > d/p for p > 1 is sharp.
In particular, it was proved for p = 2 that from / G if1 (J?) in general it does not follow that / G C(Q). □
Example 3.56. The assumption of a Lipschitz boundary. Also the assumption that Q is a Lipschitz domain is of importance.
Consider Q = {(x,y) el2 : 0 < x < 1, \y\ < xr, r > 1}, see Figure 3.4 for r = 2. The boundary is not Lipschitz in (0,0).
For u(x, y) = x~E/p with 0 < £ < r, it is
dxu = x-EJp-x (-^ = C(e,p)x-£/p-\ dyu = 0.
Using the same notation for the constant, which might take different values at different occasions, it follows that
66
3 Introduction to Sobolev Spaces
= C{s,p) [ x-£-p+r dx. Jo
This value is finite for — e — p + r > — lor for p < 1 + r — e, respectively. If one chooses r > e > 0, then it is u G W1'P(Q). But for e > 0, the function u(x) is not bounded in Q, i.e., u G" L°°(Q) and consequently u G" C_a(i7).
The unbounded values of the function are compensated in the integration by the fact that the neighborhood of the singular point (0,0) possesses a small measure. □
Chapter 4
The Ritz Method and the Galerkin Method
Remark 4-1- Contents. This chapter studies variational or weak formulations of boundary value problems of partial differential equations in Hilbert spaces. The existence and uniqueness of an appropriately defined weak solution will be discussed. The approximation of this solution with the help of finite-dimensional spaces is called Ritz method or Galerkin method. Some basic properties of this method will be proved.
In this chapter, a Hilbert space V will be considered with inner product a(-,-) : V x V —> R and norm \\v\\v = a(?j, tj)1-72. □
4.1 The Theorems of Riesz and Lax—Milgram
Theorem 4.2. Representation theorem of Riesz1.   Let f G V' be a
continuous and linear functional, then there is a uniquely determined u G V with
a{u,v) = f{v)   VtjGF. (4.1) In addition, u is the unique solution of the variational problem
F(v) = ^a(v,v) - f(v)     min V v G V. (4.2)
Proof. First, the existence of a solution u of the variational problem will be proved. Since / is continuous, it holds
\f(v)\<c\\v\\v Vvev,
from what follows that
F(v)>\\\v\\2v-C\\v\\v>-l-C2, 1 Frigyes Riesz (1880 - 1956)
67
68
4 The Ritz Method and the Galerkin Method
where in the last estimate the necessary criterion for a local minimum of the expression of the first estimate,
|       - c = o   ^   IMIv = C,
is used. Hence, the function F(-) is bounded from below and
k = inf F(v)
exists.
Let {^fcjfc^n be a sequence with F(vk) k for k oo. A straightforward calculation (parallelogram identity in Hilbert spaces) gives
hk -vi\\v + K +vi\\v = 2Kllv + 2|h||v ■ Using the linearity of /(■) and k < F(v) for all v 6 V, one obtains
\\vk ~ vi\\y
= 2 K \\2V + 2 \\vtWl - 4 (J -*±^L fv - 4f(vk) - 4f(Vl) +8/(^±^)
= 4F(vk) + 4F(Vl) - 8F
< 4F(vk) + 4F(vi) -8k^0
for k,l oo. Hence, {v/cjfc^n is a Cauchy sequence. Because V is a complete space, there exists a limit u of this sequence with u £ V. Because F(-) is continuous, it is F(u) = k and u is a solution of the variational problem.
In the next step, it will be shown that each solution of the variational problem (4.2) is also a solution of (4.1). It is for arbitrary v £ V
&(s) = F{u + sv) = — a{u + sv, u + sv) — f(u + sv)
1 e2 = — a('ii, u) + £a('ii, t?) + —a(v, v) — f(u) — sf(v).
If u is a minimum of the variational problem, then the function <P(s) has in particular a local minimum at s = 0. The necessary condition for a local minimum leads to
0 = #'(0) = a(u, v) - f(v)   for all v e V.
Finally, the uniqueness of the solution will be proved. It is sufficient to prove the uniqueness of the solution of the equation (4.1). If the solution of (4.1) is unique, then the existence of two solutions of the variational problem (4.2) would be a contradiction to the fact proved in the previous step. Let u\ and u2 be two solutions of the equation (4.1). Computing the difference of both equations gives
a ('til — u2,v) =0   for all v £ V.
This equation holds, in particular, for v = u\ — u2- Hence, ||'tii — 'ti2||v = 0, such that 'til = u2- ■
Definition 4.3. Bounded bilinear form, coercive bilinear form, V-elliptic bilinear form. Let &(•, •) : V x V —> M. be a bilinear form on the Banach space V. Then, it is bounded if
\b{u,v)\ <M\\u\\v \\v\\v   \Ju,v£ V.M > 0, (4.3)
4.1 The Theorems of Riesz and Lax—Milgram
69
where the constant M is independent of u and v. The bilinear form is coercive or ^-elliptic if
b(u,u) > m \\u\\v   VueF,m > 0, (4.4) where the constant m is independent of u. □
Remark 4-4- Application to an inner product. Let V be a Hilbert space. Then, the inner product a(-, •) is a bounded and coercive bilinear form, since by the Cauchy-Schwarz inequality
\a(u,v)\ < \\u\\v \\v\\v    V u,v G V,
and obviously a(u,u) = Hence, the constants can be chosen to be
M = 1 and m = 1.
Next, the representation theorem of Riesz will be generalized to the case of coercive and bounded bilinear forms. □
Theorem 4.5. Theorem of Lax2-Milgram3. Let &(•, •) : V x V —> R be
a bounded and coercive bilinear form on the Hilbert space V. Then, for each bounded linear functional f G V there is exactly one u G V with
b{u,v)=f{v)   VtjGF. (4.5)
Proof. One defines operators T,T'  : V -> V by
a(Tu, v) = b(u, v) V v e V,    a(T'u, v) = b(v, u) V v e V. (4.6)
These operators are linear, e.g., using that b{-, ■) is a bilinear form, one gets
a {T{a\ui + ct2U2), v) = aib(ui ,v) + ct2b{u2, v) = a {a\Tu\ + ct2Tu2, v)    V v 6 V.
Because this relation holds for all v G V, it is T{a\u\ + ct2U2) = a\Tu\ + ct2Tu2. Since b(u, ■) and b(-,u) are continuous linear functionals on V, it follows from Theorem 4.2 that the elements Tu and T'u exist and they are denned uniquely. Because the operators satisfy the relation
a(Tu, v) = b(u, v) = a(T'v, u) = a(u, T'v), (4.7)
T' is called adjoint operator of T. Setting v = Tu in (4.6) and using the boundedness of 6(-, ■) yields
\\Tu\\y = a(Tu,Tu) = b(u,Tu) < M \\u\\v \\Tu\\v => ||T«||V < M ||«||v
for all u 6 V. Hence, T is bounded. Since T is linear, it follows that T is continuous. Using the same argument, one shows that T' is also bounded and continuous. Define the bilinear form
d(u, v) := a(TT'u, v) = a(T'u, T'v)   V u, v e V, (4.8)
where (4.7) was used. Hence, this bilinear form is symmetric. Using the coercivity of 6(-, ■), the Cauchy—Schwarz inequality, the definition of and (4.8) gives
2 Peter Lax, born 1926
3 Arthur Norton Milgram (1912 - 1961)
70
4 The Ritz Method and the Galerkin Method
™2 \Mv ^ b(v,vf = a(T'v,vf < \\v\\v \\T'v\\^ = \\v\\y a(T'v,T'v) = \\v\\v d(v, v).
Applying now the boundedness of a(-, ■) and of T' yields
m2 \\v\\y < d(v,v) = a(T'v,T'v) = \\T'v\\v < M \\v\\y . (4.9)
Hence, d(-, ■) is also coercive and, since it is symmetric, it defines an inner product on V. From (4.9), one has that the norm induced by d^c^v)1/*2 is equivalent to the norm From Theorem 4.2, it follows that there is a exactly one w 6 V with
d(w,v) = f(v) ViEV.
Now, inserting u = T'w in 6(-, ■) gives with (4.6)
b(T'w, v) = a(TT'w, v) = d(w, v) = f(v)   V v e V,
hence u = T'w is a solution of (4.5).
The uniqueness of the solution is proved analogously as in the symmetric case. ■
4.2 Weak Formulation of Boundary Value Problems
Remark 4-6. Model problem. Consider the Poisson equation with homogeneous Dirichlet boundary conditions
^ = nin^crd'
ii = 0 on oil. v '
Definition 4.7. Weak formulation of (4.10). Let / G L2(Q). A weak formulation of (4.10) consists in finding u G V = ifg(i?) such that
a{u,v) = {f,v)   VtieF (4.11)
with
a(u, v) = (Vu, Vv) = / Vu(x)-Vv(x) J n
Vu(x) ■ Vv(x) dx
in
and (•, •) is the inner product in L2(i7).
Remark 4-8. On the weak formulation.
• The weak formulation is also called variational formulation.
• As usual in mathematics, 'weak' means that something holds for all appropriately chosen test functions.
• Formally, one obtains the weak formulation by multiplying the strong form of the equation (4.10) with the test function, by integrating the equation on O, and applying integration by parts. Because of the Dirichlet boundary condition, one can use as test space ifg(J?) and therefore the integral on the boundary vanishes.
Weak Formulation of Boundary Value Problems
71
The ansatz space for the solution and the test space are defined such that the arising integrals are well defined.
The weak formulation reduces the necessary regularity assumptions for the solution by the integration and the transfer of derivatives to the test function. Whereas the solution of (4.10) has to be in C2(Q) l~l C(Q), the solution of (4.11) has to be only in if1 (J?). The latter assumption is much more realistic for problems coming from applications. The regularity assumption on the right-hand side can be relaxed to / € if_1(J?). Then, the right-hand side of the weak formulation has the form
/W = </,«>/f-i(fi),«01(«)'
where the symbol {•, •,)H-1(n),H1(n) denotes the dual pairing of the spaces Hr](Q) and H'1^).
□
72
4 The Ritz Method and the Galerkin Method
Theorem 4.9. Existence and uniqueness of the weak solution. Let
f G L2(Q). There is exactly one solution of (4.11).
Proof. Because of the Poincare inequality (3.10), there is a constant C with It follows for v e Hl(n) C H1^) that
M\HHa) = (lMl'2(fi) + llvHl'2(fi))1/2 < (c \\VvfLHn) + \\VvfLHn)y/2
< C\\Vv\\L2(Q)<C\\v\\HHQ). Hence, a(-, ■) is an inner product on Hq{Q) with the induced norm
which is equivalent to the norm IHI^i^)-Define for / £ L2(J?) the linear functional
f(v):= / f(x)v(x) dx   \/ v e Hq{Q).
Using the Cauchy—Schwarz inequality (3.5) and the Poincare inequality (3.10) shows that this functional is continuous on Hq(Q)
\}{V)\ = \{f,v)\  <  \\f\\L2(Q) \ML2W  <C\\f\\L2(Q) \\Vv\\L2(n) = C\\f\\L2(Q) \MHl(n) .
Applying the representation theorem of Riesz, Theorem 4.2, gives the existence and uniqueness of the weak solution of (4.11). In addition, u(x) solves the variational problem
?(v) = - \\Vv\\l - / f(x)v(x) dx -> min    for all v e 2 JQ
(Q).
Example 4-10. A more general elliptic problem. Consider the problem
-V • (A(x)Vu) + c(x)u = / in Q C Rd, . ,
m = 0on<9J?, ( '
with A(x) G Rdxd for each point x G Q. It will be assumed that the coefficients a,ij(x) and c(x) > 0 are bounded, / G L2(Q), and that the matrix (tensor) A(x) is for all x G Q uniformly elliptic, i.e., there are positive constants m and M independent of x such that
m\\y\\22 <yTA{x)y< M\\y\\22   V y G Rd, V x G Q.
The weak form of (4.12) is obtained in the usual way by multiplying (4.12) with test functions v G Hq(O), integrating on J?, and applying integration by parts: Find u G Hr\(Q), such that
a(u,v) = (f,v) Vv£Hl{n)
4.3 The Ritz Method and the Galerkin Method 73 with
Jn
This bilinear form is bounded (exercise). The coercivity of the bilinear form is proved by using the uniform ellipticity of A(x) and the non-negativity of c(x):
Applying the Theorem of Lax-Milgram, Theorem 4.5, gives the existence and uniqueness of a weak solution of (4.12).
If the tensor is not symmetric, alJ{x) 7^ a.ji{x) for one pair then the solution cannot be characterized as the solution of a variational problem. □
4.3 The Ritz Method and the Galerkin Method
Remark J^.ll. Idea of the Ritz method. Let V be a Hilbert space with the inner product a(-, •). Consider the problem
F(v) = ^a(v,v) - f(v)     min, (4-13)
where / : V —> K is a bounded linear functional. As already proved in Theorem 4.2, there is a unique solution u G V of this variational problem which is also the unique solution of the equation
a{u,v) = f{v)   VueK (4.14)
For approximating the solution of (4.13) or (4.14) with a numerical method, it will be assumed that V has a countable orthonormal basis (Schauder basis). Then, there are finite-dimensional subspaces V"i, V2,... C V with dimVfc = k, which have the following property: for each u € V and each £ > 0 there is a K G N and a     G 14 with
\\u-uk\\v<e   Vk>K. (4.15)
Note that it is not required that there holds an inclusion of the form Vfc C V/b+i.
The Ritz approximation of (4.13) and (4.14) is defined by: Find G Vfc with
a(uk,vk) = f{vk)   ViJfcGVfc. (4.16)
74
4 The Ritz Method and the Galerkin Method
□
Lemma 4.12. Existence and uniqueness of a solution of (4.16). There exists exactly one solution of (4.16).
Proof. Finite-dimensional subspaces of Hilbert spaces are Hilbert spaces as well. For this reason, one can apply the representation theorem of Riesz, Theorem 4.2, to (4.16) which gives the statement of the lemma. In addition, the solution of (4.16) solves a minimization problem on Vk ■ H
Lemma 4.13. Best approximation property. The solution of (4.16) is the best approximation of u in Vk, i.e., it is
\\u-uk\\v=  inf \\u-vk\\v. (4-17)
Proof. Since Vk C V, one can use the test functions from Vk in the weak equation (4.14). Then, the difference of (4.14) and (4.16) gives the orthogonality, the so-called Galerkin orthogonality,
a(u-uk,vk)=0   VvkeVk- (4.18)
Hence, the error u — Uk is orthogonal to the space Vk'- u — Uk -L Vk - That means, Uk is the orthogonal projection of u onto Vk with respect of the inner product of V.
Let now Wk 6 Vk be an arbitrary element, then it follows with the Galerkin orthogonality (4.18) and the Cauchy—Schwarz inequality that
II" - "fcllv = a(u — uk,u — uk) = a(u — uk,u - (uk - wk)) = a(u - uk,u - vk)
< II'" - "fcllv II'" _ "fcllv ■
Since wk 6 Vk was arbitrary, also Vk 6 Vk is arbitrary. If \\u — Uk\\v > 0, division by I" ~ "fcllv giyes the statement of the lemma, since the error cannot be smaller than the best approximation error. If \\u — Uk\\v = 0, the statement of the lemma is trivially true.
■
Theorem 4.14. Convergence of the Ritz approximation. The Ritz approximation converges
lim \\u — Uk\\v =0.
k—>oo
Proof. The best approximation property (4.17) and property (4.15) give
ll«-«fcllv=   inf   II'" ~ "fcllv ^ £
for each e > 0 and k > K{e). Hence, the convergence is proved. ■
Remark 4-15. Formulation of the Ritz method as linear system of equations. One can use an arbitrary basis {4>i}k=i of Vk for the computation of Uk-First of all, the equation for the Ritz approximation (4.16) is satisfied for all vk G V/c if and only if it is satisfied for each basis function (j>l. This statement follows from the linearity of both sides of the equation with respect to the test function and from the fact that each function Vk G Vk can be represented
4.3 The Ritz Method and the Galerkin Method
75
as linear combination of the basis functions. Let vk = '*l2l=lOLi4>i, then from (4.16), it follows that
fc fc a(uk,Vk) = ^2 ata{uk, 4>t) = ^ a>if(<j>i) = f(vk). fc=i fc=i
This equation is satisfied if a(uk,(j>i) = f{(j>i), i = l,...,fc. On the other hand, if (4.16) holds then it holds in particular for each basis function (j>l.
Now, one uses as ansatz for the solution also a linear combination of the basis functions
fc
Uk = ^h^j
with unknown coefficients u] € K. Using as test functions the basis functions yields
fc fc
^"(^fc^.) = s^Ja{4>],4>t)u3 = /(&)>  « = l, • • • ,k.
This equation is equivalent to the linear system of equations Au = /, where
A = K)^=1 = a(<^J,^)lfcjJ=1
is called stiffness matrix. Note that the order of the indices is different for the entries of the matrix and the arguments of the inner product. The right-hand side is a vector of length fc with the entries f% = f(<f>i), i = 1,..., fc.
Using the one-to-one mapping between the coefficient vector (v1,... ,vk)T and the element vk = X«=i one can show that the matrix A is symmetric and positive definite (exercise)
A = AT <^=> a(v, w) = a(w, v)   V v, w € Vk, xTAx>0ioix^0-^^a(v,v)>0 Vti6l4,!)/0.
□
Remark 4-16. The case of a bounded and coercive bilinear form. If &(•,•) is bounded and coercive, but not symmetric, it is possible to approximate the solution of (4.5) with the same idea as for the Ritz method. In this case, it is called Galerkin method. The discrete problem consists in finding € Vfc such that
b{uk, vk) = f{vk)   Vwfce yfc. (4.19)
□
Lemma 4.17. Existence and uniqueness of a solution of (4.19). There is exactly one solution of (4.19).
76
4 The Ritz Method and the Galerkin Method
Proof. The statement of the lemma follows directly from the Theorem of Lax—Milgram, Theorem 4.5. B
Remark 4-18. On the discrete solution. The discrete solution is not the orthogonal projection into Vk in the case of a bounded and coercive bilinear form, which is not the inner product of V. □
Lemma 4.19. Lemma of Cea4, error estimate. Let b :  V x V —> M be
a bounded and coercive bilinear form on the Hilbert space V and let f € V be a bounded linear functional. Let u be the solution of (4.5) and Uk be the solution of (4.19), then the following error estimate holds
M
\\u-uk\\v<— inf \\u-vk\\v, (4.20) m vkevk
where the constants M and m are given in (4.3) and (4.4).
Proof. Considering the difference of the continuous equation (4.5) and the discrete equation (4.19), one obtains the error equation
b(u - uk,vk) = 0 Vi>fceVfc,
i.e., Galerkin orthogonality holds. With (4.4), the Galerkin orthogonality, and (4.3), it follows that
, 1 1
\\u ~ uk II v — —b(u — uk:u — uk) = —b(u — uk :u — vk) m m
M
< - \W - Uk\\y   \\U - Vk\\y  , VvkeVk,
m
from what the statement of the lemma follows immediately. ■
Remark 4-20. On the best approximation error. It follows from estimate (4.20) that the error is bounded by a multiple of the best approximation error, where the factor depends on properties of the bilinear form &(•, •). Thus, concerning error estimates for concrete finite-dimensional spaces, the study of the best approximation error will be of importance. □
Remark 4-21. The corresponding linear system of equations. The corresponding linear system of equations is derived analogously to the symmetric case. The system matrix is still positive definite but not symmetric. □
Remark 4-22. Choice of the basis. The most important issue of the Ritz and Galerkin method is the choice of the spaces 14, or more concretely, the choice of an appropriate basis {4>t}f=1 that spans the space Vk- From the point of view of numerics, there are the requirements that:
• it should be possible to compute the entries a%] of the stiffness matrix efficiently,
• and that the matrix A should be sparse.
□
4 Jean Cea, born 1932
Chapter 5
Finite Element Methods
5.1 Finite Element Spaces
Remark 5.1. Mesh cells, faces, edges, vertices. A mesh cell if is a compact polyhedron in Kd, d G {2, 3}, whose interior is not empty. The boundary dK of K consists of m-dimensional linear manifolds (points, pieces of straight lines, pieces of planes), 0 < m < d— 1, which are called m-faces. The 0-faces are the vertices of the mesh cell, the 1-faces are the edges, and the (d—l)-faces are just called faces. □
Remark 5.2. Finite-dimensional spaces defined on K. Let SEN. Finite element methods use finite-dimensional spaces P(K) C CB(K) that are defined on K. In general, P(K) consists of polynomials. The dimension of P(K) will be denoted by dimP(K) = NK. □
Example 5.3. The space P(K) = Pi (if).The space consisting of linear polynomials on a mesh cell K is denoted by P1(K):
There are d + 1 unknown coefficients a%, i = 0,..., d, such that dim Pi (K) =
Remark 5.4- Linear functionals defined on P(K), nodal functionals. For the definition of finite elements, linear functional that are defined on P(K) are of importance. These functionals are called nodal functionals.
Consider linear and continuous functionals <L>k,i, ■ ■ ■ ,@k,Nk '■ C"(K) —> M which are linearly independent. There are different types of functionals that can be utilized in finite element methods:
• point values: <P(v) = v(x), x G K,
• point values of a first partial derivative: <P(v) = dtv(x), x G K,
NK =d + l.
□
77
78
5 Finite Element Methods
• point values of the normal derivative on a face E of K: <P(v) = \7v(x)-ue, rig is the outward pointing unit normal vector on E,
• integral mean values on K: <P(v) = jj^ jK v(x) dx,
• integral mean values on faces E: $(v) = p^y JEv(s) ds.
The smoothness parameter s has to be chosen in such a way that the func-tionals <Pff,i) • • • j &k,nk are continuous. If, e.g., a functional requires the evaluation of a partial derivative or a normal derivative, then one has to choose at least s = 1. For the other functional given above, s = 0 is sufficient. □
Definition 5.5. Unisolvence of P(K) with respect to the functionals
<L>k,i, ■ ■ ■ ,&k,nk• The space P(K) is called unisolvent with respect to the functionals <L>k,i, ■ ■ ■ ,®k,nk if there is for each a € ~M.Nk , a = (a1,... ,aNK)T, exactly one p G P(K) with
$Kti{p)=ai, \<i<NK.
□
Remark 5.6. Local basis. Unisolvence means that for each vector a € M.Nk , a = (ai,... ,a^tK)T, there is exactly one element in P(K) such that al is the image of the i-th functional, i = 1,..., N^.
Choosing in particular the Cartesian1 unit vectors for a, then it follows from the unisolvence that a set {(J>k,í}íÍ^i exists with <j>k,í £ P{K) and
®K,i(<i>K,i) = i,j = l,...,NK.
Consequently, the set {(pK^}1^ forms a basis of P(K). This basis is called local basis. □
Remark 5.7. Transform of an arbitrary basis to the local basis. If an arbitrary basis {ft}^ of P{K) is known, then the local basis can be computed by solving a linear system of equations. To this end, represent the local basis in terms of the known basis
nk
4>k,j = ^2 CjkPk,   cjk EK, j = 1,..., NK, fc=i
with unknown coefficients Cjk- Applying the definition of the local basis leads to the linear system of equations
nk
&k,i{4>k,j) = ^CjfcOifc = Sij,    i,j = l,...,NK,    aik = $k,i{pk)-fc=l
Because of the unisolvence, the matrix A = (alJ) is non-singular and the coefficients Cjk are determined uniquely. □
1 René Descartes (1596 - 1650)
5.1 Finite Element Spaces
79
Example 5.8. Local basis for the space of linear functions on the reference triangle. Consider the reference triangle K with the vertices (0,0), (1,0), and (0,1). A linear space on K is spanned by the functions l,x,y. Let the functionals be defined by the values of the functions in the vertices of the reference triangle. Then, the given basis is not a local basis because the function 1 does not vanish at the vertices.
Consider first the vertex (0,0). A linear basis function ax + by + c that has the value 1 in (0,0) and that vanishes in the other vertices has to satisfy the following set of equations
The solution is a = —1,6 = —l,c = 1. The two other basis functions of the local basis are x and y, such that the local basis has the form {1 — x — y, x, y}.
□
Remark 5.9. Triangulation, grid, mesh, grid cell. For the definition of global finite element spaces, a decomposition of the domain Q into polyhedra K is needed. This decomposition is called triangulation Th and the polyhedra K are called mesh cells. The union of the polyhedra is called grid or mesh.
A triangulation is called admissible, see the definition in (Ciarlet, 1978, p. 38, p. 51), if:
• It holds Q = \JKe-rhK.
• Each mesh cell K G Th is closed and the interior K is non-empty.
• For distinct mesh cells K\ and K2 there holds K\ l~l K2 = 0-
• For each K G Th, the boundary dK is Lipschitz continuous.
• The intersection of two mesh cells is either empty or a common m-face, m G {0,... ,d - 1}.
be continuous linear functionals of the same types as given in Remark 5.4, where for each K, v\K G P(K) has to be understood in the sense that the polynomial in K is extended continuously to the boundary of K. The restriction of the functionals to CS(K) defines a set of local functionals &k,i-i ■ ■ ■ 1 &k,nk j where it is assumed that the local functionals are unisol-vent on P(K). The union of all mesh cells Kj, for which there is a p G P(Kj) with $i(p) 7^ 0, will be denoted by ul. □
Example 5.11. On subdomains cul. Consider the two-dimensional case and let $l be defined as nodal value of a function in x G K. If x G K, then oj, = K. In the case that a; is on a face of K but not in a vertex, then oj, is the union
□
Remark 5.10. Global and local functionals. Let
<PU...,<PN : {v€L°°(Ü) : v\K€P(K)}^R
80 5 Finite Element Methods
Fig. 5.1 Subdomains UJi.
of K and the other mesh cell whose boundary contains this face. Last, if x is a vertex of K, then oj, is the union of all mesh cells that possess this vertex, see Figure 5.1. □
Definition 5.12. Finite element space, global basis. A function v(x) defined on Q with v\k G P{K) for all K G Th is called continuous with respect to a global functional $l defined in Remark 5.10 if
$i{v\Kl) = $i{v\K2), yK1,K2euii.
The space
S = |« G L°°(Q) : v\k G P{K) and v is continuous with respect to
$i,i = l,...,N}
is called finite element space.
The global basis {4>j}f=i of S is defined by the condition
4>j & S,   <?l(<^J) = 5lJ,   i,j = 1,...,JV.
□
Example 5.13. Piecewise linear global basis function. Figure 5.2 shows a piece-wise linear global basis function in two dimensions. Because of its form, such a function is called hat function. □
Remark 5.14- On global basis functions. A global basis function coincides on each mesh cell with a local basis function. This property implies the uniqueness of the global basis functions.
Whether the continuity with respect to {^>i}|^=1 implies the continuity of the finite element functions depends on the functionals that define the finite element space. □
Definition 5.15. Parametric finite elements. Let if be a reference mesh cell with the local space P(K), the local functionals <Pi,... and a class of bijective mappings {Fx '■ K —> K}. A finite element space is called a parametric finite element space if:
• The images {K} of {FK} form the set of mesh cells.
• The local spaces are given by
P(K) = \p : p = p0FK\p€P(K)y (5.1)
• The local functionals are defined by
&k,Mx)) = ^ (*(*)) = ^ . I5-2)
where x = (xi,...,id)T are the coordinates of the reference mesh cell and it holds x = Fk(x), v = v o Fx-
□
Remark 5.16. Motivations for using parametric finite elements. Definition 5.12 of finite elements spaces is very general. For instance, different types of mesh cells are allowed. However, as well the finite element theory as the implementation of finite element methods become much simpler if only parametric finite elements are considered. □
82
5 Finite Element Methods
5.2 Finite Elements on Simplices
Definition 5.17. d-simplex. A d-simplex K C Md is the convex hull of
(d + 1) points ai,
■ , CLd+l G
which form the vertices of K.
Remark 5.18. On d-simplices. It will be always assumed that the simplex is not degenerated, i.e., its d-dimensional measure is positive. This property is equivalent to the non-singularity of the matrix (exercise)
A-
( an	ai2 ■		-i\
«21	«22 •	• a,2,d-\	-l
adi	ad2 •	■ a,d,d-\	-l
\ 1	1 .	. 1	
where a* = (au,a2l,.
j = l,
,d + l.
For d = 2, the simplices are the triangles and for d = 3 they are the tetrahedra. □
Definition 5.19. Barycentric coordinates. Since K is the convex hull of the points {fflij-fij1, the parametrization of K with a convex combination of the vertices reads as follows
K -
C d+1 d+1 ^
: x = ^2 Xidi, 0 < Aj < 1, ^ Aj = 1 I .
I i=i i=i J
The coefficients Ai,..., Xd+i are called barycentric coordinates of x G K. □
Remark 5.20. On barycentric coordinates.
• From the definition, it follows that the barycentric coordinates are the solution of the linear system of equations
d+1 d+1
^2 ajiK = Xj,    1 < j < d,   ^ Aj = 1.
i=i i=i
Since the system matrix is non-singular, see Remark 5.18, the barycentric coordinates are determined uniquely.
• The barycentric coordinates of the vertex aI5 i = 1,..., d + 1, of the simplex are X% = 1 and X3• = 0 if i j. Since Xl(aJ) = 8%], the barycentric coordinate X% can be identified with the linear function that has the value 1 in the vertex at and that vanishes in all other vertices a3 with j i.
• The barycenter of the simplex is given by
j     d+1 d+1 j
5.2 Finite Elements on Simplices
83
Fig. 5.3 The unit simplices in two and three dimensions.
Hence, its barycentric coordinates are \ =       + 1), i = 1,..., d + 1.
□
Remark 5.21. Simplicial reference mesh cells. A commonly used reference mesh cell for triangles and tetrahedra is the unit simplex
K =
J* G Rd : < 1, £j > 0, i = 1,... ,d| ,
see Figure 5.3. The class {Fk} of admissible mappings are the bijective affine mappings
FKx = BKx + b,     BK€ Rdxd, det (BK) / 0, f> € Rd. (5.3)
The images of these mappings generate the set of the non-degenerated simplices {K} C Rd. □
Definition 5.22. Affine family of simplicial finite elements. Given a simplicial reference mesh cell K, affine mappings {Fk}, and an unisolvent set of functionals on K. Using (5.1) and (5.2), one obtains a local finite element space on each non-degenerated simplex. The set of these local spaces is called affine family of simplicial finite elements. □
Definition 5.23. Polynomial space Pk. Let x = [x\,... ,x^)T, k G NU
{0}, and a = (ai,.. ., c<d)T■ Then, the polynomial space Pk is given by
(■ d d \
Pk = span I Yl XT = xa : a.£NU{0} for i = 1,..., d, ^ ai < k \ ■ U=l 1=1 J
Remark 5.24- Lagrangian2 finite elements. In many examples given below, the linear functionals on the reference mesh cell K are the values of the
Joseph-Louis de Lagrange (1736 — 1813)
84
5 Finite Element Methods
Fig. 5.4 The finite element Pq(K).
a2
Fig. 5.5 The finite element P\{K).
polynomials with the same barycentric coordinates as on the general mesh cell K. Finite elements whose linear functionals are values of the polynomials on certain points in K are called Lagrangian finite elements. □
Example 5.25. Pq : piecewise constant finite element. The piecewise constant finite element space consists of discontinuous functions. The linear functional is the value of the polynomial in the barycenter of the mesh cell, see Figure 5.4. It is dim P0 (K) = 1. □
Example 5.26. Pi : conforming piecewise linear finite element. This finite element space is a subspace of C(Q). The linear functionals are the values of the function in the vertices of the mesh cells, see Figure 5.5. It follows that AhnP^K) = d+ 1.
The local basis for the functionals {(Pl(v) = v(at), i = l,...,d+ 1} is {^i}i=i smce ^i(-^j) = compare Remark 5.20. Since a local basis exists, the functionals are unisolvent with respect to the polynomial space Pi(K).
Now, it will be shown that the corresponding finite element space consists of continuous functions. Let K\,K2 be two mesh cells with the common face E and let u £ Pi(= S). The restriction of on E is a linear function on E as well as the restriction of vK2 on E. It has to be shown that both linear functions are identical. A linear function on the (d— l)-dimensional face E is
5.2 Finite Elements on Simplices
85
Fig. 5.6 The finite element P2(K).
uniquely determined with d linearly independent functionals that are defined on E. These functionals can be chosen to be the values of the function in the d vertices of E. The functionals in S are continuous by the definition of S. Thus, it must hold that both restrictions on E have the same values in the vertices of E. Hence, it is \g = v^2\E and the functions from Pj are continuous. □
Example 5.27. P2 : conforming piecewise quadratic finite element. This finite element space is also a subspace of C(Q). It consists of piecewise quadratic functions. The functionals are the values of the functions in the d+1 vertices of the mesh cell and the values of the functions in the centers of the edges, see Figure 5.6. Since each vertex is connected to each other vertex, there are J2i=1 i = d(d+l)/2 edges. Hence, it follows that dim P2(K) = (d+l)(d+2)/2.
The part of the local basis that belongs to the functionals {&i(v) = v(at), i = l,...,d+l}, is given by
{&(A) = Ai(2Ai-l),   i = l,...,d+l}.
Denote the center of the edges between the vertices at and a3 by al3. The corresponding part of the local basis is given by
{<j>ij = 4AlAJ,   i,j = l,...,d+l, i <j}.
The unisolvence follows from the fact that there exists a local basis. The continuity of the corresponding finite element space is shown in the same way as for the Pi finite element. The restriction of a quadratic function defined in a mesh cell to a face E is a quadratic function on that face. Hence, the function on E is determined uniquely with d(d + l)/2 linearly independent functionals on E.
The functions (j>lJ are called in two dimensions edge bubble functions. □
Example 5.28. P3 : conforming piecewise cubic finite element. This finite element space consists of continuous piecewise cubic functions. It is a subspace of C(Q). The functionals in a mesh cell K are defined to be the values in
86
5 Finite Element Methods
the vertices ((d + 1) values), two values on each edge (dividing the edge in three parts of equal length) (2^i=1i = d(d + 1) values), and the values in the barycenter of the 2-faces of K, see Figure 5.7. Each 2-face of K is defined by three vertices. If one considers for each vertex all possible pairs with other vertices, then each 2-face is counted three times. Hence, there are (d + l)(d — l)d/6 2-faces. The dimension of PsiK) is given by
dimP3 (A) = (<*+!)+ d(d + 1) + (___W±_1 = (*+l)(* + 2)(<* + 3)
6 6
For the functional
v(a,i), i = 1,..., d + 1, (vertex), v(auj), i, j = 1,... ,d+ l,i 7^ j,     (point on edge), v(aijk), i = l,...,d+l,i<j<k, (point on 2-face) |,
the local basis is given by
= ^(3Ai - l)(3Ai - 2),   4>uj(A) = ^XiXjtfXi - 1),
<?>ijfc(A) = 27AiAjAfc
In two dimensions, the function </>ijfc(A) is called cell bubble function. □
Example 5.29. Cubic Hermite3 element. The finite element space is a sub-space of C(Q), its dimension is (d + l)(d + 2)(d + 3)/6 and the functional are the values of the function in the vertices of the mesh cell ((d +1) values), the value of the barycenter at the 2-faces of K ((d+ l)(d — l)d/6 values), and the partial derivatives at the vertices (d(d+ 1) values), see Figure 5.8. The
3 Charles Hermite (1822 - 1901)
5.2 Finite Elements on Simplices
87
Fig. 5.8 The cubic Hermite element.
dimension is the same as for the P3 element. Hence, the local polynomials can be defined to be cubic.
This finite element does not define an affine family in the strict sense, because partial derivatives on the reference cell are mapped to directional derivatives on the physical cell. Concretely, the functionals for the partial derivatives <Pl(v) = div(0) on the reference cell are mapped to the functionals &i{v) = dt%v(a), where a = Fk{0) and tl are the directions of edges which are adjacent to a, i.e., a is an end point of this edge. This property suffices to control all first derivatives. One has to take care of this property in the implementation of this finite element.
Because of this property, one can use the derivatives in the direction of the edges as functionals
v(fli), (vertices) V«(tij) • [a,j — cii), i,j = l,...,d — l,ij^j, (directional deriv.) v(a,ijk), i < j < k, (2-faces)
with the corresponding local basis
&(A) = -2A? + 3A? - 7Xi £J<fc,^,Ml A,-Afc) <?>ij(A) = AjAj(2Aj - \j - 1), <^ijfc(A) = 27AjAjAfc.
The proof of the unisolvence can be found in the literature.
Here, the continuity of the functions will be shown only for d = 2. Let K\, K2 be two mesh cells with the common edge E and the unit tangential vector t. Let V\,V2 be the end points of E. The restrictions v\k1,v\k2 to E satisfy four conditions
v\Kl(Vi)=v\K2(Vl),   dtv\Kl(Vi) = dtv\K2(Vi), t = l,2.
Since both restrictions are cubic polynomials and four conditions have to be satisfied, their values coincide on E.
®ijk{v) =
88
5 Finite Element Methods
The cubic Hermite finite element possesses an advantage in comparison with the P3 finite element. For d = 2, it holds for a regular triangulation Th that
where #(•) denotes the number of triangles, nodes, and edges, respectively. Hence, the dimension of P3 is approximately #{V) + 2#(P) + #(if) ~ 7#(V), whereas the dimension of the cubic Hermite element is approximately 3#(y)+#(if) ~ 5#(F). This difference comes from the fact that both spaces are different proper subspaces of the space of all continuous piecewise cubic functions. The elements of both spaces are continuous functions, but for the functions of the cubic Hermite finite element, in addition, the first derivatives are continuous at the nodes. That means, these two spaces are different finite element spaces whose degree of the local polynomial space is the same (cubic). One can see at this example the importance of the functionals for the definition of the global finite element space. □
Example 5.30. Pfc : non-conforming linear finite element, Crouzeix-Raviart finite element, Crouzeix & Raviart (1973). This finite element consists of piecewise linear but discontinuous functions. The functionals are given by the values of the functions in the barycenters of the faces such that dimP"c(K) = (d + 1). It follows from the definition of the finite element space, Definition 5.12, that the functions from Pfc are continuous in the barycenter of the faces
P"c = {v G L2(Q) : v\k £ Pi (if), v(x) is continuous at the barycenter
Equivalently, the functionals can be defined to be the integral mean values on the faces and then the global space is defined to be
where £(K) is the set of all (d — l)-dimensional faces of K.
For the description of this finite element, one defines the functionals by
$i(v) = «(aj_i]i+i) for d = 2,   <t>t(v) = ^(al_2,I-i,I+i) for d = 3,
where the points are the barycenters of the faces with the vertices that correspond to the indices, see Figure 5.9. This system is unisolvent with the local basis
#(k) «2#(n #(p)«2#(n
of all faces}.
(5.4)
(5.5)
<f>i(\) = 1 — d\i,   i = 1,..., d + 1.
5.3 Finite Elements on Parallelepipeds and Quadrilaterals
89
□
5.3 Finite Elements on Parallelepipeds and Quadrilaterals
Remark 5.31. Reference mesh cells, reference map to parallelepipeds. One can find in the literature two reference cells: the unit cube [0,l]d and the large unit cube [—1, l]d. It does not matter which reference cell is chosen. Here, the large unit cube will be used: K = [—1, l]d. The class of admissible reference maps {Fk\ to parallelepipeds consists of bijective affine mappings of the form
FKx = BKx + b,   BK G Rdxd, b G Rd.
If Bx is a diagonal matrix, then K is mapped to d-rectangles.
The class of mesh cells that is obtained in this way is not sufficient to triangulate general domains. If one wants to use more general mesh cells than parallelepipeds, then the class of admissible reference maps has to be enlarged, see Remark 5.40. □
Definition 5.32. Polynomial space Qk. Let x = (x\,... ,Xd)T and denote by a = (qi, ..., ay)T a multi-index. Then, the polynomial space Qk is given by
Qk = span < Y[ XT = xa '■ 0 < Qj < k for i = 1,..., d > .
□
Example 5.33. Qi vs. P1. The space Q1 consists of all polynomials that are d-linear. Let d = 2, then it is
Qi = span{l,a;, y, xy},
whereas
90
5 Finite Element Methods
01 o2 Fig. 5.10 The finite element Qi(K).
Pi = span{l, x, y}.
□
Remark 5.34- Finite elements on d-rectangles. For simplicity of presentation, the examples below consider d-rectangles. In this case, the finite elements are just tensor products of one-dimensional finite elements. In particular, the basis functions can be written as products of one-dimensional basis functions.
□
Example 5.35. Qo : piecewise constant finite element. Similarly to the Pg space, the space Qo consists of piecewise constant, discontinuous functions. The functional is the value of the function in the barycenter of the mesh cell K and it holds dimQ0(K) = 1. □
Example 5.36. Qi : conforming piecewise d-linear finite element. This finite element space is a subspace of C(Q). The functional are the values of the function in the vertices of the mesh cell, see Figure 5.10. Hence, it is dimQi(ir) =2d.
The one-dimensional local basis functions, which will be used for the tensor product, are given by
<i>i{x)
1
(1
02 (£)
1
(1
With these functions, e.g., the basis functions in two dimensions are computed
by
<M*)<My), <M£)<Ms/), <M£)<Ms/), foi^foiy)-
The continuity of the functions of the finite element space Qi is proved in the same way as for simplicial finite elements. It is used that the restriction of a function from Qk(K) to a face E is a function from the space Qk(E), k > 1. □
5.3 Finite Elements on Parallelepipeds and Quadrilaterals
91
Fig. 5.11 The finite element Q2(K).
Fig. 5.12 The finite element Q3(K).
Example 5.37. Q2 •' conforming piecewise d-quadratic finite element. It holds that Q2 C C(Q). The functionals in one dimension are the values of the function at both ends of the interval and in the center of the interval, see Figure 5.11. In d dimensions, they are the corresponding values of the tensor product of the intervals. It follows that dimQ2(if) = 3d.
The one-dimensional basis function on the reference interval are defined
by
<£i(£) =-£),    4>2(x) = (1 - x)(\ +£),   4>3(x) = ^(1 + x)x.
The basis function Yli=i <t>2{xi) is called cell bubble function. □
Example 5.38. Q3 : conforming piecewise d-cubic finite element. This finite element space is a subspace oiC(O). The functionals on the reference interval are given by the values at the end of the interval and the values at the points x = —1/3, x = 1/3. In multiple dimensions, it is the corresponding tensor product, see Figure 5.12. The dimension of the local space is dim Q3(K) = 4d. The one-dimensional basis functions in the reference interval are given by
92
5 Finite Element Methods
Fig. 5.13 The finite element Q;ot(_ff).
Mx) = ~{ßx + l)(3z - l)(x - 1),   fa{x) = ^-(x + l)(3z - l)(x - 1), 16 16
16 16
Example 5.39. Q\ot : rotated non-conforming element of lowest order, Ran-nacher-Turek element, Rannacher & Turek (1992): This finite element space is a generalization of the Pfc finite element to quadrilateral and hexahedral mesh cells. It consists of discontinuous functions that are continuous at the barycenter of the faces. The dimension of the local finite element space is dim Q\ot (K) = 2d. The space on the reference mesh cell is defined by
Qiot (k) = {P ■ V e span{l,£,y,£2 - y2}} for d = 2,
Q\ot (k) = {p : p G span{l,£ ,y, z,x2 - y2,y2 - z2}} for d = 3.
Note that the transformed space
Q\ot(K) = {P = poFK\p€ Q\ot(K)}
contains polynomials of the form ax2 — by2, where a,b depend on Fk-For d = 2, the local basis on the reference cell is given by
4>i(x,y) = -'^(x2 -y2)        +       Mx,y) = ^(x2 - y2) + if + i,
<M£,y) = -|(*2 -y2) + \y + 7,   <M£,y) = l(x2 -y2) -\x + \.
(5.6)
Analogously to the Crouzeix-Raviart finite element, the functionals can be defined as point values of the functions in the barycenters of the faces, see Figure 5.13, or as integral mean values of the functions at the faces.
5.4 Transform of Integrals
93
Consequently, the finite element spaces are defined in the same way as (5.4) or (5.5), with Pfc(K) replaced by Q\ot(K).
In the code ParMooN Wilbrandt et al. (2017), the mean value oriented Q\ot finite element space is implemented for two dimensions and the point value oriented Qjot finite element space for three dimensions. For d = 3, the integrals on the faces of mesh cells, whose equality is required in the mean value oriented Q5ot finite element space, involve a weighting function which depends on the particular mesh cell K. The computation of these weighting functions for all mesh cells is an additional computational overhead. For this reason, it was suggested in (Schieweck, 1997, p. 21) to use for d = 3 the simpler point value oriented form of the Q\ot finite element. □
Remark 5.4-0. Parametric mappings. The image of an affine mapping of the reference mesh cell K = [—1, l]d, d G {2,3}, is a parallelepiped. If one wants to consider finite elements on general d-quadrilaterals, then the class of admissible reference maps has to be enlarged.
The simplest non-affine parametric finite element on quadrilaterals in two dimensions uses bilinear mappings. Let K = [—1, l]2 and let
Fk(x) = ( FM\) = ( ^ J ^ J ai3f( J ™ ), F>K G QU i = 1, 2, \Fk\x) j     \a21 + a22x + a23y + a24xy j
be a bilinear mapping from K on the class of admissible quadrilaterals. A quadrilateral K is called admissible if
• the length of all edges of K is larger than zero,
• the interior angles of K are smaller than ir, i.e., K is convex.
This class contains, e.g., trapezoids and rhombi. □
Remark 5-41- Parametric finite element functions. The functions of the local space P[K) on the mesh cell K are defined by p = p o F^1. These functions are in general rational functions. However, using d-linear mappings, then the restriction of Fk on an edge of K is an affine map. For instance, in the case of the Qi finite element, the functions on K are linear functions on each edge of K. It follows that the functions of the corresponding finite element space are continuous, compare Example 5.26. □
5.4 Transform of Integrals
Remark 5.42. Motivation. The transformation of integrals from the reference mesh cell to mesh cells of the grid and vice versa is used as well for the analysis as for the implementation of finite element methods. This section provides an overview of the most important formulae for transformations.
Let K c Rd be the reference mesh cell, K be an arbitrary mesh cell, and FK : K —> K with x = FK(x) be the reference map. It is assumed that the
94 5 Finite Element Methods
reference map is a continuous differentiable one-to-one map. The inverse map is denoted by F^1 : K —> K. For the integral transforms, the derivatives (Jacobians) of FK and F^1 are needed
dx' dx'
DFK(x)ij =DFK1{x)l3 = i,j = l,...,d.
Remark 5.4-3. Integral with a function without derivatives. This integral transforms with the standard rule of integral transforms
/ v(x) dx = / v{x)\detDFK{x)\ dx, (5.7) J k J k
where v(x) = v(FK(x)). □ Remark 5.44- Transform of derivatives. Using the chain rule, one obtains
= V&v(x)-[{DF^1(FK(x)))T)i
= ([[DFK1(FK(x)))T>)i • V*v(x), (5.8) dv dv , .dx,     „ . .
^-(-) = E^-(-)^ = v^)-((^(a;)f
= Vv(x)-{{DFK(FK1(x)))T)i. (5.9)
The index i denotes the i-th row of a matrix. Derivatives on the reference mesh cell are marked with a symbol on the operator. □
Remark 5.45. Integrals with a gradients. Using the rule for transforming integrals and (5.8) gives
/ b(x) ■ Vv(x) dx
J k
= Jb{FK{x))- [(DFKr)T (FK(x))} V&v{x)\detDFK{x)\ dx. (5.10) Similarly, one obtains
5.4 Transform of Integrals 95
Vv(x) • Vw(x) dx
= J [{DFKr)T (FK(x))] V&v(x) ■ [{DFKr)T (FK(x))] V&w(x) x \det DFK(x)\ dx. (5.11)
□
Example 5.4-6. Affine transform. The most important class of reference maps are affine transforms (5.3), where the invertible matrix Bk and the vector b are constants. It follows that
x = BKX (x-b)= BKxx - BKxb.
In this case, there are
DFK = BK,   DFKX = B'1,   det DFK = det (BK).
One obtains for the integral transforms from (5.7), (5.10), and (5.11)
/ v{x) dx = \det{BK)\ / v{x) dx, (5.12) J k J k
l b(x) ■ Vv{x) dx = \det{BK)\ / b(FK(x)) ■ BKTVAv(x) dx, (5.13) J k J k
Vv{x) ■ Vw{x) dx = |det (BK)\ I BKTV&v(x) ■ BKTV&w(x) dx.
k
a
k k
(5.14)
Setting v(x) = 1 in (5.12) yields
\det(BK)\ = \^\. (5.15) \K\
Chapter 6 Interpolation
Remark 6.1. Motivation. Variational forms of partial differential equations use functions in Sobolev spaces. The solution of these equations shall be approximated with the Ritz method in finite-dimensional spaces, the finite element spaces. The best possible approximation of an arbitrary function from the Sobolev space by a finite element function is a factor in the upper bound for the finite element error, e.g., see the Lemma of Cea, estimate (4.20).
This section studies the approximation quality of finite element spaces. Estimates are proved for interpolants of functions. Interpolation estimates are of course upper bounds of the best approximation error and they can serve as factors in finite element error estimates. □
6.1 Interpolation in Sobolev Spaces by Polynomials
Lemma 6.2. Unique determination of a polynomial with integral conditions. Let Q be a bounded domain in Md with Lipschitz boundary. Let m € NU{0} be given and let for all derivatives with multi-index a, \a\ < m, a value aa € K be prescribed. Then, there is a uniquely determined polynomial p € Pm{Q) such that
dap(x) dx = aa,    \a\ < m. (6-1)
Proof. Let p 6 Prn{Q) be an arbitrary polynomial. It has the form
v(&) = ^2 bßxß-
\ß\<m
Inserting this representation in (6.1) leads to a linear system of equations Mb = a with
M = (Maß), Maß = / daxß dx, b = (bß), a= (aa), Ja
97
98
6 Interpolation
for \a\ , \j3\ < m. Since M is a squared matrix, the linear system of equations possesses a unique solution if and only if M is non-singular.
The proof is performed by contradiction. Assume that M is singular. Then, there exists a non-trivial solution of the homogeneous system. That means, there is a polynomial q e Pm(O) \ {0} with
/  dctq(x) dx = 0 for all \a\ < m. J Q
The polynomial q{x) has the representation q{x) = X]|/3|<m cf3x^• Now, one can choose a Cj3 ^ 0 with maximal value \j3\. Then, it is dpq{x) = Ccp = const ^ 0, where C > 0 comes from the differentiation rule for polynomials, which is a contradiction to the vanishing of the integral for dpq(x). ■
Remark 6.3. To Lemma 6.2. Lemma 6.2 states that a polynomial is uniquely determined if a condition on the integral on Q is prescribed for each derivative. □
Lemma 6.4. Poincare-type inequality. Denote by Dkv(x), k G NU {0}, the total derivative of order k of a function v(x), e.g., for k = 1 the gradient ofv(x). Let Q be convex and be included into a ball of radius R. Let I G NU{0} with k < I and let p G K with p G [1, oo). Assume that v G Wl'p(Q) satisfies
I dav(x) dx = 0 for all \a\ < I — 1, Jn
in
then it holds the estimate
\\Dk4L,{n)<CRl'k\\Dl4L,{nv
where the constant C does not depend on Q and on v(x).
Proof. There is nothing to prove if k = /. In addition, it suffices to prove the lemma for k = 0 and 1 = 1, since the general case follows by applying the result to da_v{x).
Since Q is assumed to be convex, the integral mean value theorem can be written in the form
v{x) — v{y) = I   S7v{tx + (1 — t)y) ■ (x — y) dt,    x,y £ j?. Jo
Integration with respect to y yields
Jo y yields
v(x) /    dy -      v(y) dy =       /   Vv(tx + (1 - t)y) ■ (x - y) • Jn        Jn JnJo
tion that the se
-y) dt dy.
It follows from the assumption that the second integral on the left-hand side vanishes that
v(x) = — J    J   Vv(tx + (1 - t)y) ■(x-y) dt ,
Now, taking the absolute value on both sides, using that the absolute value of an integral is estimated from above by the integral of the absolute value, applying the Cauchy—Schwarz inequality for vectors (3.3), and the estimate \\x — y\\2 < 2R yields
6.1 Interpolation in Sobolev Spaces by Polynomials 1
99
f f Vv(tx + (1 - t)y) ■ (x - y) dt d J q jo
i- f f \Vv(tx + (1 -t)y) ■ (x -y)\ dt < u\ Jq Jo
f  f  ||V«(ta!+(l-t)3/)||2 dtdy. J n Jo
2R
(6.2)
Then, (6.2) is raised to the power p and integrated with respect to x. One obtains with Holder's inequality (3.4), with p~1 + q~1 = 1        p/q — p = p(l/q — 1) = —1, that
fnWx)\" dx < JQ (jQjQ  WVv(tx + f1 -*)3/)ll2 dt dv^j dx
■iVif
Jq \JqJo
CUP
\q dt dy
\q\rli
CK"
\Q\
j^l^  \\Vv(tx + (l-t)y)\\l dtdy^
f (f f \\Vv(tx+(l-t)y)\\% dtdy)dx. Jq \Jn Jo /
Applying the theorem of Fubini allows the commutation of the integration
f \v(x)\p dx < f   f (f \\Vv(tx+ (1 -t)y)\\% dy] dx dt.
Jq \ft\ Jo Jq \Jq J
Using the integral mean value theorem in one dimension gives that there is a to 6 [0, 1] such that
J \v(x)\p dx < ^ J (^j \\Wv(t0x + (1- to)y)\\l dy^jdx.
The function ||Vt/(a;)||2 will be extended to m,d by zero and the extension will be also denoted by ||Vv(ic)||2- Then, it is
f \v(x)\p dx<^- ( (f \\Vv(tQx + (1 - tQ)y)\{p dy) dx. (6.3) Jq w Jq Vr** J
Let to £ [0,1/2]. Since the domain of integration is Md, a substitution of variables tox + (1 — to)y = z can be applied and leads to
/   \\Vv(t0x + (1 - to)y)\\l dy = —L- f   \\Vv(z)\\p dz < 2 \\Vv\\PLP(Q) , Jmd I - to JRd ^ !
since 1/(1 — to) < 2. Inserting this expression in (6.3) gives
j \v(x)\p dx<2CRp\\Vv\\PLP(Q).
If to > 1/2 then one changes the roles of x and y, applies the theorem of Fubini to change the sequence of integration, and uses the same arguments. ■
100
6 Interpolation
Remark 6.5. On Lemma 6.4- Lemma 6.4 proves an inequality of Poincaré-type. It says that it is possible to estimate the Lp(f2) norm of a lower derivative of a function v(x) by the same norm of a higher derivative if the integral mean values of some lower derivatives vanish.
An important application of Lemma 6.4 is in the proof of the Bramble1-Hilbert2 lemma. The Bramble-Hilbert lemma considers a continuous linear functional that is defined on a Sobolev space and that vanishes for all polynomials of degree less than or equal to m. It states that the value of the functional can be estimated by the Lebesgue norm of the (m + l)th total derivative of the functions from this Sobolev space. □
Theorem 6.6. Bramble-Hilbert lemma. Let m G N U {0}, p G [l,oo], and F : Wrn+1'P(Q) —> K be a continuous linear functional, and let the conditions of Lemma 6.2 and Lemma 6.4 be satisfied. Let
F(p) = 0 VpGPm(i2),
then there is a constant C{Q), which is independent of v and F, such that
\F[v)\ < C{Q)\\Dm+1v\\LT(n)   V v G Wm+1-P(i7).
Proof. Let v e Wm+1'p(f2). It follows from Lemma 6.2 that there is a polynomial from Pm(Q) with
da(v + p)(x) dx = 0 for \a\ < m.
J q
Lemma 6.4 gives, with Z = m + 1 and considering each term in IHIj^m+i.p^j individually, the estimate
\\v + v\\wm+i,HQ) < C(Q) \\Dm+Hv + P)\\LP(Q) = C(Q) \\Dm+1v\\LHn) . From the vanishing of F for p £ Pm(J?) and the continuity of _F, it follows that \F(v)\ = \F(v + p)\<C\\v + p\\wm+1,p(Q) ^CWWD^vW^^.
Remark 6.7. Strategy for estimating the interpolation error. The Bramble-Hilbert lemma, more precisely Lemma 6.4, will be used for estimating the interpolation error for finite elements. The strategy is as follows:
• Show first the estimate on the reference mesh cell k.
• Transform the estimate on an arbitrary mesh cell k to the reference mesh cell k.
• Apply the estimate on k.
• Transform back to k.
One has to study what happens if the transforms are applied to the estimate.
□
1 James H. Bramble, born 1930
2 Stephen R. Hilbert
6.1 Interpolation in Sobolev Spaces by Polynomials
101
Remark 6.8. Assumptions, definition of the interpolant. Let K C Md,d € {2,3}, be a reference mesh cell (compact polyhedron), P(K) a polynomial space of dimension N, and <P1,..., <PN : CS(K) —> M continuous linear functionals. It will be assumed that the space P(K) is unisolvent with respect to these functionals. Then, there is a local basis 4>i, ■ ■ ■ ,4>n £ P{K). Consider v G CS(K), then the interpolant Ip/b G P(K) is defined by
n
Iki>{&) = ^^)k{x).
The operator Ig, is a continuous and linear operator from CS(K) to P(K). From the linearity, it follows that I^ is the identity on P(K)
Ikp = P VpGP(ii).
□
Example 6.9. Interpolation operators.
• Let if C Md be an arbitrary reference cell, P(K) = Po(-řf), and
é(v) = j-í-r /" v(x) dx.
\k\ Jk
The functional <I> is bounded, and hence continuous, on C°(K) since
i- i   i r \&\
\${v)\< -J—i /  \v{x)\ dx < \—pmax|í)(á)| = ||i)||Co/^-N. 1      1     \k\ Jk \k\ ňek
For the constant function 1 G P0(K), it is $(1) = 1^0. Hence, {</>} = {1} is the local basis and the space is unisolvent with respect to <P. The operator
I^v(x) = S(v)4>(x) = -j—r / v(x) dx \k\ Jk
is an integral mean value operator, i.e., each continuous function on K will be approximated by a constant function whose value equals the integral mean value, see Figure 6.1
• It is possible to define <P(v) = v(xq) for an arbitrary point Xq £ K. This functional is also linear and continuous in C°(K). The interpolation operator 1^ defined in this way interpolates each continuous function by a constant function whose value is equal to the value of the function at xq, see also Figure 6.1.
102
6 Interpolation
Interpolation operators which are defined by using values of functions, are
called Lagrangian interpolation operators. This example demonstrates that the interpolation operator 1^ depends on P(K) and on the functionals <Pt. □
Theorem 6.10. Interpolation error estimate on a reference mesh
cell. Let Pm(K) C P(K), let p G [l,oo), and let S G N U {0} such that (m + 1 — š)p > d > (m — s)p and š > s, where s appears in the definition of the interpolation operator. Then there is a constant C that is independent of v(x) such that
ll*-^*L™+1.P(K)^c'll£,m+li)lli>(K)   VfiGWm+^(if). (6.4)
Proof. Since K is bounded, one has the Sobolev imbedding, Theorem 3.52,
Wm+1'P(K) = jy(™+i-s)+s,P(_K-) _> Cg(K).
Because K is convex, the imbedding CS(K) CS(K) is compact, see (Adams, 1975, Theorem 1.31), such that the interpolation operator is well defined in WrnJrl'p(K). From the identity of the interpolation operator in Pm(^)j the triangle inequality, the bound-edness of the interpolation operator (it is a linear and continuous operator mapping CS(K) -> P(K) C Wm+1'P(K)), and the Sobolev imbedding, one obtains for q e Pm(K)
II5-/k5IIm/-+i.p(k) = \\'d + 9 - Ik('s + 9)\\wm+i.P(ii)
< \\'" + S\\wm+i,T{k) + \\lkt" + ")\\w™+i.v(k)
< 11« + <í\\Wm+l,p(k) + c II5 + 9llc"(k)
< C\\v + q\\wm+1,Pik).
Now, q(x) is chosen such that
da(v + q) dx = 0   V \a\ < m
L
IK
holds. Hence, the assumptions of Lemma 6.4 are satisfied. It follows that
6.1 Interpolation in Sobolev Spaces by Polynomials 103
\\v + q\\wm+,,r{k) < C \\D™ + Hv + q)\\LHk) = C ||Dm + 1«||„(it) ■
Definition 6.11. Quasi-uniform and regular family of triangulations,
(Brenner & Scott, 2008, Def. 4.4.13). Let {Th} with 0 < h < 1, be a family of triangulations such that
max fix < h diam(J?),
KeTh
where hx is the diameter of K = Fk(K), i.e., the largest distance of two points that are contained in K. The family is called to be quasi-uniform, if there exists a C > 0 such that
min pk > Ch diam(J?) (6.5)
KeTh
for all h G (0,1], where px is the diameter of the largest ball contained in K.
The family is called to be regular, if there is exists a C > 0 such that for all K G Th and for all h G (0,1]
pK > ChK.
Remark 6.12. Assumptions on the reference mapping and the triangulation. For deriving the interpolation error estimate for arbitrary mesh cells K, and finally for the finite element space, one has to study the properties of the mapping from K to K and of the inverse mapping. Here, only the case of an affine family of finite elements whose mesh cells are generated by affine mappings
FK x = BK x + b,
will be considered, see (5.3), where Bk is a non-singular d x d matrix and b is a d vector. For the global estimate, a quasi-uniform family of triangulations will be considered. □
Lemma 6.13. Estimates of matrix norms. For each matrix norm \\-\\, one has the estimates
\\BK\\<ChK,    WB^W <ChK\ (6.6) where the constants depend on the matrix norm.
Proof. Since K is a Lipschitz domain with polyhedral boundary, it contains a ball B(xq, r) with xq £ K and some r > 0. Hence, xq + y £ K for all ||y||2 = r- It follows that the images
xq = Bkxq -\-b,    x = Bk(xq + y) + b = xq + BKy are contained in K. Hence, one obtains for all y
104
6 Interpolation
\\BKy\\2 = \\x ~ xo\\2 < hK-
Now, it holds for the spectral norm that
||BK||2 =   sup   ||BKz||2 = -   sup   ||BKz||2 < —.
Isll2 — 1 Ilsll2 — '
A bound of this form, with a possible different constant, holds also for all other matrix norms since all matrix norms are equivalent, see Remark 3.34.
The estimate for JB^1] proceeds in the same way with interchanging the roles of K and K. M
Theorem 6.14. Local interpolation estimate. Let an affine family of finite elements be given by its reference cell K, the functionals {^i}, and a space of polynomials P(K). Let all assumptions of Theorem 6.10 be satisfied. Then, for all v € Wm+1'P(K), p € [l,oo), there is a constant C, which is independent of v, such that
\\D\v-IKv)\\LT(K)<Ch™+1-k\\Dm+\\\LT(K),   0<fe<m + l. (6.7)
Proof. The idea of the proof consists in transforming the left-hand side of (6.7) to the reference cell, using the interpolation estimate on the reference cell, and transforming back.
i) . Denote the elements of the matrices Bx and B^} by bij and 6^ 1\ respectively. Since H-B^Hoo = max^j \bij \ is also a matrix norm, it holds that
\bij\<ChK,    Jfc^1'! < Ch'1. (6.8)
Using element-wise estimates for the matrix Br (Leibniz formula for determinants), one obtains
\detBK\ < ChdK,    jdetB"1! < Ch~d. (6.9)
ii) . The next step consists in proving that the transformed interpolation operator is equal to the natural interpolation operator on K. The latter one is given by
N
lKV = ^$K,i(V)(i)K,i: (6-10) i = l
where {(pK,i} is the basis of the space
P(K) = {p : K^R : p = p o F'1 ,p e P(K)} , which satisfies &kk,j) =       The functionals are defined by
=&i(voFK) =&i(v). (6.11) Hence, it follows for v = <pj o F^-1 from the condition on the local basis on K that
&k,i(h ofk1) = ®i(4>j) = 5ih i.e., the local basis on K is given by (pK,j = 4>j ° F^-1. Using (6.11) and (6.10), one gets
6.1 Interpolation in Sobolev Spaces by Polynomials 105
N N I N \
IfrV = ^2$i(v)(j>i = ^2<PK>i(vo F'1) (f>K>i oFK = i J oFK
i = l i = l %-v-' \i = l /
—v
= IKv oFK.
Consequently, I^-v is transformed correctly. in). One obtains with the chain rule
dv(x) dv(x) ^(-l)     dv(x) dv(x) ^
dxi 4  dx-;   -i"1    ' 4  dx-; J%
j=i       J j=i J
It follows with (6.8) that (with each derivative one obtains an additional factor of or B^1, respectively)
||D^(:z)||2 < ChKk \\D%v(A)\\2 ,    ||D|«(£)||2 < ChkK \\Dkv(x)\\2. One gets with (6.9)
/ \\Dkv(x)\\P2 dx < ChKkp \detBK\ f \\D%.v{x)\\P dx<ChKkp+d f \\D%.v{x)\\P dx J k J k J k
(6.12)
and
/  \\D%H£)\\l dx<Chkp\detBK1\ ( \\Dkv(x)\\P2 dx < Chk^-d f \\D*v(x)\\p dx.
J k J k J k
(6.13)
Using now the interpolation estimate on the reference cell (6.4) yields
ll^(*-J**)|IL(*)^CllD*+1*IIL.(*)>    0<fe<m + l. (6.14) It follows with (6.12), (6.14), and (6.13) that
\\Dk(v-IKv)\\PLHK) < ChKkp+d\\D%(v-Ikv)\\PLHk)
<ChKkp+d\\D^v\\PLp(k)
< rh(-m+1^k'>p llnm+1,illp
s ^nK v\\LP(K) ■
Taking the p-th root proves the statement of the theorem. ■
Remark 6.15. On estimate (6.7).
• Note that the power of      does not depend on p and d.
• Consider a quasi-uniform triangulation and define
h = max \ Hk\-
KeTh
Then, one obtains by summing over all mesh cells an interpolation estimate for the global finite element space
106 6 Interpolation
/ \1/p
\\Dk(v-Ihv)\\LP(ü)= (  £ \\Dk(v-IKv)\\PLP{K) \KeTh ,
< ( £ Ch^-k)*\\D™+^LT(K)
i/p
\KeTh
<Ch^m+1-^\\Dm+1v\\LT{n). (6.15)
□
Corollary 6.16. Finite element error estimate. Let u(x) be the solution of the model problem (4.10) with u € Hm+1{Q) and let uh{x) be the solution of the corresponding finite element problem. Consider a family of quasi-uniform triangulations and let the finite element spaces Vh contain polynomials of degree m. Then, the following finite element error estimate holds
||V(« - uh)\\L2(n) < Chm\\Dm+1u\\LHn) = Chm \u\Hm+1{n). (6.16)
Proof. The statement follows by combining Lemma 4.13 (for v = Hq(Q)) and (6.15) P(u-uh)\\L.{Q) = jfvhp(u-v*)\\L2{Q)
< ||V(«-/h«)||L2(fi) <Chm\u\Hm+1(n).
Remark 6.17. To (6.16). Note that Lemma 4.13 provides only information about the error in the norm on the left-hand side of (6.16), but not in other norms. □
6.2 Inverse Estimate
Remark 6.18. On inverse estimates. The approach for proving interpolation error estimates can be used also to prove so-called inverse estimates. With inverse estimates, a norm of a higher order derivative of a finite element function is estimated by a norm of a lower order derivative of this function. Likewise, norms in different Lebesgue spaces are estimated. One obtains as penalty a factor with negative powers of the diameter of the mesh cell. □
Theorem 6.19. Inverse estimate. Let 0 < k < I be natural numbers and let p,q € [l,oo]. Then there is a constant C;nv, which depends only on k,l,p,q,K, P(K), such that
\\Dlvh\\        < Cimh(tl)-d{p~1-q~1) \\Dkvh\\ V vh G P{K). (6.17)
6.2 Inverse Estimate
107
Proof. In the first step, (6.17) is shown for = 1 and k = 0 on the reference mesh cell. Since all norms are equivalent in finite-dimensional spaces, one obtains
\\Dl^\LHK)^\\ih\\w^(K)<CPh\\LP{K)     V0heP(#). (6.18)
If k > 0, then one sets
P(K) = {dcvh  : vh e P(K), \a\ = k} , which is also a space consisting of polynomials. The application of (6.18) to P(K) gives
llD'«1L,(*) = E \\Dl-k(d^h)\\Lq(k)<c J2 pavh\\LP(k)
| ol\ — k \ol | — k
This estimate is transformed to an arbitrary mesh cell K analogously as for the interpolation error estimates, compare the proof of Theorem 6.14. From the estimates (6.12) and (6.13) for the transformations, one obtains
\\nl„h\\        < ch~l+d/q IInlt,h\\        < ch~l+d/q IInkt,h\\
\\u " \\h1(K) ^ unK " \\Li(k) ^ unK " \\Lv(k)
r,     , k-l+d/q-d/p II „fc kii = ^mv«K \\U   0   II LP (if) ■
Remark 6.20. On the proof. The crucial point in the proof is the equivalence of all norms in finite-dimensional spaces. Such a property does not hold in infinite-dimensional spaces. □
Corollary 6.21. Global inverse estimate. Let p = q and let {Th} be a
quasi-uniform family of triangulations of Q, then
\\DlA\^n)^C^hk~l\\DkA\^ny (6-19)
where
i/p
Remark 6.22. On W'Wrv.h^y The cell-wise definition of the norm is important for k > 2 or / > 2 since in these cases finite element functions generally do not possess the regularity for the global norm to be well defined. It is also important for / > 1 and non-conforming finite element functions. □
Chapter 7
Finite Element Methods for Second Order Elliptic Equations
7.1 General Convergence Theorems
Remark 7.1. Motivation. In Section 5.1, non-conforming finite element methods have been introduced, i.e., methods where the finite element space Vh is not a subspace of V, which is the space in the definition of the continuous variational problem. The property Vh (£ V is given for the Crouzeix-Raviart and the Rannacher-Turek element. Another case of non-conformity is given if the domain does not possess a polyhedral boundary and one has to apply some approximation of the boundary.
For non-conforming methods, the finite element approach is not longer a Ritz method. Hence, the convergence proof from Theorem 4.14 cannot be applied in this case. In addition, in practice, one is interested also in the order of convergence in other norms than I'H^ or one has to take into account that the values of the bilinear or linear form need to be approximated numerically. The abstract convergence theorem, which will be proved in this section, allows the numerical analysis of complex finite element methods. □
Remark 7.2. Notations, Assumptions. Let {h > 0} be a set of mesh widths and let Sh,Vh normed spaces of functions which are defined on domains {Qh C Rd}. It will be assumed that the space Sh has a finite dimension and that Sh and Vh possess a common norm ||-||fe. In the application of the abstract theory, Sh will be a finite element space and Vh is defined such that the restriction and/or extension of the solution of the continuous problem to Qh is contained in Vh. The index h indicates that Vh might depend on h but not that Vh is finite-dimensional. Strictly speaking, the modified solution of the continuous problem does not solve the given problem any longer. Hence, it is consequent that the continuous problem does not appear explicitly in the abstract theory.
Given the bilinear forms
109
110 7 Finite Element Methods for Second Order Elliptic Equations
ah : Sh xSh ^ R,
ah : {Sh + Vh) x (Sh + Vh) R.
Let the bilinear form ah be regular in the sense that there is a constant m > 0, which is independent of h, such that for each vh g Sh there is a wh g 5* with 11 wh 11, =1 such that
ii Wh
m\\vh\\h<ah(vh,wh). (7.1)
This condition is equivalent to the requirement that the stiffness matrix A with the entries al3 = ah(4>3,4>i), where {(f>t} is a basis of Sh, is uniformly non-singular, i.e., its regularity is independent of h. For the second bilinear form, only its boundedness will be assumed
~ah(u,v) <M\\u\\h\\v\\h   Vu,v£Sh+Vh. (7.2)
Let the linear functionals {fh{-)} : Sh —> R be given. Then, the following discrete problems will be considered: Find uh g Sh with
ah(uh,vh) = fh{vh)   \/vheSh. (7.3)
Because the stiffness matrix is assumed to be non-singular, there is a unique solution of (7.3). □
Theorem 7.3. Abstract finite element error estimate. Let the conditions (7.1) and (7.2) be satisfied and let uh be the solution of (7.3). Then, the following error estimate holds for each u g Vh
\\~      fell   ^n  -t   J||-      fell   ^ \ah(vh,wh)-ah(vh,wh)\
\\u — u   ,<C   mi   < \\u — v    ,+  sup--r—rv,-
" «"es"\" wh/sh \\wh\\h
\ah(u,wh)-fh(wh)\ +C sup--—f--(7.4)
with C = C(m,M).
Proof. Because of (7.1), there is for each vh e Sh a wh e Sh with ||«Jh||h = 1 and
m\\uh -vh\\h < ah(uh -vh,wh). Using the definition of uh from (7.3), one obtains
m\\uh -vh\\h < fh(wh) -ah(vh,wh)+ah(vh,wh) +ah(u-vh,wh) -ah(u,wh). From (7.2) and ||«flh||h = 1, it follows that
ah(u-vh,wh) <M\\u-vh\\h. Rearranging the terms appropriately and using ||tt;ft/ H^^H^H^ = 1 gives
7.1 General Convergence Theorems 111
m\\uh - vh\\i < M\\ü - u'1 IL + sup i lift —      n lift ,,
whĚSh
\hh(vh,wh) - ah(vh,wh)\
\ah(ü,wh) - fh(wh)\ +   sup   J-j,—,-'-. (7.5)
»'es* II™ IL
Applying the triangle inequality
and inserting the estimate (7.5) gives (7.4).
Remark 7.4- To Theorem 7.3.
• An important special case of this theorem is the case that the stiffness matrix is uniformly positive definite, i.e., the condition
m\\vh\\l<ah(vh,vh)   VvheSh (7.6)
is satisfied. Dividing (7.6) by ||^fe||fe reveals that condition (7.1) is implied by (7.6).
• If the continuous problem is also defined with the bilinear form ah(-,-), then
\ah(vh,wh) -ah(vh,wh)\ whesh \\w Ik
can be considered as consistency error of the bilinear forms and the term
\ah(u,wh) - fh(wh)\
sup -n~hi\-
™"es" II™ IL
as consistency error of the right-hand sides.
Theorem 7.5. First Strang1 lemma Let Sh be a conforming finite element space, i.e., Sh c V, with \\-\\h = \\-\\v and let the space Vh be independent of h. Consider a continuous problem of the form
dh(u,v) = f(v) VtjGF,
then the following error estimate holds
ah(vh,wh)-ah(vh,wh) v - ~ „hTbh 1 ir    " \\v     T* \uňh\\
\u — u\\<C  inf   < \\u — v   L+ sup
vhesh I 11 v    ,„fceSfc \\'w"\\v
_ \f(wh)-fh(™h)\
+C   Sup   J-—rj-'-.
Proof. The statement of this theorem follows directly from Theorem 7.3.
1 Gilbert Strang, born 1934
112 7 Finite Element Methods for Second Order Elliptic Equations
7.2 Finite Element Method with the Non-conforming Crouzeix—Raviart Element
Remark 7.6. The continuous problem. Let Q C Md, d G {2, 3}, be a bounded domain with polygonal Lipschitz boundary. Let
Lu = f   in n,   u = 0   on dO, (7.7)
where the operator is given by
Lu = -V • (AVu)
with
A(x) = (ai](x))lJ=1,   0^ G W^n^p > d, (7.8) It will be assumed that there are two positive real numbers m, M such that
m\\£\\l<£,TA(x)£<M\\£\\l   V£GKd,a;GJ?. (7.9)
From the Sobolev inequality, Theorem 3.52, it follows that al3 G L°°(Q). With
a{u,v) = / (A(x)Vu(x)) ■ Vv(x) dx J n
and the Cauchy-Schwarz inequality, one obtains
\a{u,v)\ < \\A\\L^{n) jjVu(x) ■ Vv(x)\ dx < C\\X7u\\L2(n) \\X7v\\L2(n)
for all u,v G Hq(Q). In addition, it follows from (7.9) that m || Vm||^2(j7) < a(u, u) V?iGifg(i7).
Hence, the bilinear form is bounded and elliptic. Using the Theorem of Lax-Milgram, Theorem 4.5, it follows that there es a unique weak solution u G Hl(n) of (7.7) with
a{u,v) = f{v)   Vu G Hl(Q).
□
Remark 7.7. Assumptions and the discrete problem. The non-conforming Crouzeix-Raviart finite element Pfc was introduced in Example 5.30. To simplify the presentation, it will be restricted here on the two-dimensional case. In addition, to avoid the estimate of the error coming from approximating the domain, it will be assumed that Q is a convex domain with polygonal boundary. It can be shown that in this case the boundary is Lipschitz.
7.2 Finite Element Method with the Non-conforming Crouzeix—Raviart Element
113
Fig. 7.1 Function from Pfc.
Let Th be a regular triangulation of Q with triangles. Let P"c (nc -non-conforming) be denote the finite element space of piecewise linear functions which are continuous at the midpoints of the edges. This space is nonconforming if it is applied for the discretization of a second order elliptic equation since the continuous problem is given in Hq(O) and the functions of Hq(Q) do not possess jumps. The functions of Pfc have generally jumps, see Figure 7.1, and they are not weakly differentiable. In addition, the space is also non-conforming with respect to the boundary condition, which is not satisfied exactly. The functions from Pf c vanish in the midpoint of the edges at the boundary. However, in the other points at the boundary, their value is generally not equal to zero.
The bilinear form
ah(u,v)= / (A(x)Vu(x))-Vv(x) dx   Vu,o e Hq (J?) + P"c.
Then, the non-conforming finite element method is given by: Find uh G Pj1 with
The goal of this section consists in proving the linear convergence with
the solution of the continuous problem (7.7) is smooth, i.e., that u G H2(Q), that / G L2(Q), and that the coefficients alJ(x) are weakly differentiable
will be extended to H^(O) + Pfc by
ah(uh,vh) = (f,vh) ^vh€P1c.
with bounded derivatives.
□
114
7 Finite Element Methods for Second Order Elliptic Equations
Remark 7.8. The error equation. The first step of proving an error estimate consists in deriving an equation for the error. To this end, multiply the continuous problem (7.7) with a test function from vh G Pfc, integrate the product on O, and apply integration by parts on each triangle. This approach gives
(f,vh) = - £  / V-(A(x)Vu(x))vh(x)dx
= J2   [ (A(x)Vu(x))-Vvh(x) dx KeTh K
-   J2    f     (AS)^U(S))-nK(s)vh(s) ds
KeT"JaK
= ah(u,vh)-Y   Í   {A{s)Vu{s))-nK{s)vh{s) ds, rrrtk J ok
Kerh '
where is the unit outer normal at the edges of the triangles. Subtracting the finite element equation, one obtains
ah{u-uh,vh) = Y   [   {A{s)X7u{s))-nK{s)vh{s)ds   V vh G P^. (7.10) ,77^^ JdK
KeTh
Lemma 7.9. Estimate of the right-hand side of the error equation (7.10). Assume that u G H2(0) and al3 G W1'00{0), i, j = 1, 2, then it is
V   [   A{s)Vu{s)-nK{s)vh{s) ds
KeTh
< Ch\\u\\H2(n) \\vh\\h.
Proof. Every edge of the triangulation which is in Q appears exactly twice in the boundary integrals on dK. The corresponding unit normals possess opposite signs. One can choose for each edge one unit normal and then one can write the integrals in the form
J2 [ [\(A(s)Vu(s))-nE(s)vh(s)\]E ds = J2 [ (A(s)Vu(s)) ■ nE(s) [\vh\]E(s) ds,
e e
where the sum is taken over all edges {E}. Here, denotes the jump of vh
i] (q\ _ í vh\Kl(s) -vh\K2(s) s UeW     \vh(s) s
e E CÍ2,
e E c da,
where tie is directed from K\ to K2 or it is the outer normal on dfl. For writing the integrals in this form, it was used that V'ti(s), A(s), and ti_e(s) are almost everywhere continuous, such that these functions can be written as factor in front of the jumps. Because of the continuity condition for the functions from P^c and the homogeneous Dirichlet boundary condition, it is for all vh £ P^c that (P) = 0 for the midpoints P of all
edges. From the linearity of the functions on the edges, it follows that
7.2 Finite Element Method with the Non-conforming Crouzeix—Raviart Element
115
V
		1	
		ex ^\ E Ki	
-1		1	
Fig. 7.2 Reference configuration.
[|^|]s(s) ds = 0   VS. (7.11)
Let E be an arbitrary edge in Q which belongs to the triangles K\ and Ki- The next goal consists in proving the estimate
/ (A(s)Vu(s))-nE(s) [\vh\]    (s) ds
J E
< Ch\MH2(Ki) (||V,"||L2(Ki) + ||V,"||L2(K2)) . (7.12)
To this end, one uses a reference configuration [K\, K2, EJ, where K\ is the unit triangle and K2 is the triangle that is obtained by reflecting the unit triangle at the y-axis. The common edge E is the interval (0,1) on the y-axis. The unit normal on E will be chosen to be the Cartesian unit vector ex, see Figure 7.2. This choice is the other way around than in the definition of the jump, but it is just for simplicity of notation and it does not influence the estimate. The reference configuration can be transformed to (Ki, K2, E) by a map which is continuous and on both triangles _£Q affine. For this map one, can prove the same properties for the transform as proved in Chapter 6.
Using (7.11), the Cauchy—Schwarz inequality, and the trace theorem, one obtains for an arbitrary constant a £ M
f (A(s)Vu(s)) -ex [|^|]^ ds=  f ((A(s)Vu(s)) ■ ex - a) [|^|]^ ds
j E j E
<C||(^Vfl)-e,-a||Hl(Ai)||[|«h|]^||i.a(fi)(7.13) In particular, one can choose a such that
((i(s)V'ti(s)) ■ ex - a) ds = 0.
Using first that (a2 + 62)1/2 < a + b for a, b > 0, then the L2(J?) term in the first factor of the right-hand side of (7.13) can be bounded using the estimate from Lemma 6.4 for k = 0 and 1 = 1 and the choice of a
116
7 Finite Element Methods for Second Order Elliptic Equations
|| (iw) -ex -a\\Hlriti)
< (|| (AVu) ■ ex -a\\LHki) + || V ((AVu) ■ ex - a)\\L,(ki))
<C\\v((AVu)-ex-a)\\L2(i<:i)
= C\\V((AVu)-ex)\\L2(I(i).
To estimate the second factor, the trace theorem and the equivalence of norms in finite-dimensional spaces are applied
lltKIUL^) <c(\\vh\\HH&l) + \\vh\\HH&2))
<C(llV51L2(K1)+llV5'*ll^(K2))- (7'14)
To apply the norm equivalence, one has to prove that the terms in the last line are in fact norms. Let the terms in the last line be zero, then it follows that vh = c\ in K\ and vh = C2 in K<i - Because vh is continuous in the midpoint of E, one finds that c\ = C2 and consequently that = 0. Hence, also the left-hand side of the estimate is zero. It
follows that the right-hand side of estimate (7.14) defines a norm in the quotient space of Pfc with respect to = 0.
Altogether, one obtains for the reference configuration
I (A(s)Vu(s)) ■ ex [\vh\]6 ds
J E
< C ||V ((iVu) ■ ex)\\L2(ki) (||V«h||La(Ai) + ||V«h||La(Aa)) .
This estimate has to be transformed to the triple (Ki, K2, E). In this step, one gets for the integral on the edge the factor C (Ch for V and C'h~1 for ds). For the product of the norms on the right-hand side, one obtains the factor Ch (Ch for the first factor and C for the second factor). In addition, one uses that A(s) and all first order derivatives of A(s) are bounded to estimated the first term on the right-hand side (exercise). In summary, (7.12) is proved.
The statement of the lemma follows by summing over all edges and by applying on the right-hand side the Cauchy—Schwarz inequality. ■
Theorem 7.10. Finite element error estimate. Let the assumptions of Lemma 7.9 be satisfied, then it holds the following error estimate
\\u ~ uh\\h - Ch WuWhz({?) \\u ~ uh\\h + ch2 \\u\\2h*(n) •
Proof. Applying Lemma 7.9, it follows from the error equation (7.10) that
\ah(u-uh,vh)\ < Ch\\u\\H2(n) \\vh\\h    V/e P?c.
Let Ih : Hq(Q) P^c be an interpolation operator with optimal interpolation order in H'll^. Then, one obtains with the Cauchy—Schwarz inequality, the triangle inequality, and the interpolation estimate
7.3 L2(J?) Error Estimate 117
||« - uh\\l = ah(u - uh,u - uh) = ah(u - uh,u - Ihu) + ah(u - uh,Ihu - uh)
< \ah(u-uh,u-Ihu)\ +Ch\\u\\H2(Q) \\lhu-uh\\h
< \\u-uh\\h \\u-Ihu\\h + Ch\\u\\H2(Q) (\\lhu-u\\h + ||«-«h||J
< Ch\\u-uh\\h \\u\\H2{n) + Ch\\u\\H2{n) (h\\u\\H2{n) + ||« - ■
■
Remark 7.11. To the error estimate. If h is sufficiently small, than the second term of the error estimate is of higher order and this term can be absorbed in the constant of the first term. One obtains the asymptotic error estimate
i.e., first order convergence. □
7.3 L2(f2) Error Estimate
Remark 7.12. Motivation. A method is called quasi-optimal in a given norm, if the order of the method is the same as the optimal approximation order. Already for one dimension, one can show that at most linear convergence in H1(Q) can be achieved for the best approximation in Pj. This statement can be already verified with the function v(x) = x2. Hence, all considered methods so far are quasi-optimal in the energy norm.
However, the best approximation error in L2(Q) is of one order higher than the best approximation error in H1(Q). A natural question is whether finite element methods converge also of higher order with respect to the error in L2(Q) than with respect to the error in the energy norm.
In this section, it will be shown that one can obtain for finite element methods a higher order of convergence in L2{Q) than in if1 (12). However, there are more restrictive assumptions to prove this property in comparison with the convergence proof for the energy norm. □
Remark 7.13. Model problem. Let Q c Md, d g {2, 3}, be a convex polyhedral domain with Lipschitz boundary. The model problem has the form
-Au = f in J2,   u = 0 ondO. (7.15)
For proving an error estimate in P2(i7), the regularity of the solution of (7.15) plays an essential role. □
Definition 7.14. m-regular differential operator. Let L be a second order differential operator. This operator is called m-regular, m > 2, if for all / g Hrn~2(Q) the solutions of Lu = / in Q, u = 0 on dQ, are in the space Hm(Q) and the following estimate holds
118 7 Finite Element Methods for Second Order Elliptic Equations
IMIff-(R) < c + c\\u\\H^n) • (7-16)
□
Remark 7.15. On the m-regularity.
• The definition is formulated in a way that it can be applied also if the solution of the problem is not unique.
• For the Laplacian, the term I^H^-i,-^ can be estimated by l/H^^-) such that with (7.16) one obtains (exercise)
IMIff2(R) ^ C ll/lli2(fi) •
• Many regularity results can be found in the literature. Loosely speaking, they say that regularity is given if the data of the problem (coefficients of the operator, boundary of the domain) are sufficiently regular. For instance, an elliptic operator in divergence form (A = V • (AV)) is 2-regular if the coefficients are from W1'P(Q), p > 1, and if dQ is a C2 boundary. Another important result is the 2-regularity of the Laplacian on a convex domain. A comprehensive overview on regularity results can be found in Grisvard (1985).
□
Remark 7.16. Variational form and finite element formulation of the model problem. The variational form of (7.15) is: Find u G Hq(Q) with
(V«,V«) = (/,«)  Vu effort).
The Pi finite element space, with zero boundary conditions, will be used for the discretization. Then, the finite element problem reads as follows: Find uh G Pi such that
(W\ Vvh) = (f,vh)   V/ePi. (7.17)
□
Theorem 7.17. Finite element error estimates. Letu(x) be the solution of (7.15), let (7.15) be 2-regular, and let uh(x) be the solution of (7.17). Then, the following error estimates hold
P{u-uh)\\L2(n)<Ch\\f\\LHn),
\\u-uh\\L*(n) < Ch2 Wfh^n) ■
Proof. With the error estimate in H1(f2), Corollary 6.16, and the 2-regularity, one obtains
llV('"-'"h)IL2(fi) <Ch\\u\\H2in) <Ch\\f\\L2(Q). (7.18)
For proving the L2(J?) error estimate, let w £ Hq{Q) be the unique solution of the so-called dual problem
7.4 Outlook
119
(Vv, Vw) = (u-uh,v) V«effo(i7).
For a symmetric differential operator, the dual problem has the same form like the original (primal) problem. Hence, the dual problem is also 2-regular and it holds the estimate
IMIji2(fi) <C\\u-uh\\L2(Qy
For performing the error estimate, the Galerkin orthogonality of the error is utilized
(V(« - uh), Vvh) = (Vu, Vvh) - (Vuh, Vvh) = (/, vh) - (/, vh) = 0
for all vh 6 P\. Now, the error u — uh is used as test function v in the dual problem. Let Ihw be the interpolant of w in Pi. Using the Galerkin orthogonality, the interpolation estimate, and the regularity of w, one obtains
= (V(u-uh),V™) = (V(u-uh),V(w-Ihw))
< \\V(u-uh)\\L2(Q)\\V(W-IhW)\\L2(Q)
< Ch\\w\\H2(Q) \\V(u-uh)\\L2(Q) <Ch\\u-uh\\L2{n)\\V(u-uh)\\L2{n).
Finally, division by \\u — '"^H^^) an<^ ^ne application of the already known error estimate (7.18) for ||V('ti — 'u^)||^2^^ are used f°r completing the proof of the theorem. ■
7.4 Outlook
Remark 7.18. Outlook to forthcoming classes. This class provided an introduction to numerical methods for solving partial differential equations and the numerical analysis of these methods. There are many further aspects that will be covered in forthcoming classes.
Further aspects for elliptic problems.
• Adaptive methods and a posteriori error estimators. It will be shown how it is possible to estimate the error of the computed solution only using known quantities and in which ways one can decide where it makes sense to refine the mesh and where not.
• Multigrid methods. Multigrid methods are for certain classes of problems optimal solvers.
• Numerical analysis of problems with other boundary conditions or taking into account quadrature rules.
Time-dependent problems. As mentioned in Remark 1.7, standard approaches for the numerical solution of time-dependent problems are based on solving stationary problems in each discrete time.
• The numerical analysis of discretizations of time-dependent problems has some new aspects, but also many tools from the analysis of steady-state problems are used.
120
7 Finite Element Methods for Second Order Elliptic Equations
Convection-diffusion equations. Convection-diffusion equations are of importance in many applications. Generally, the convection (first order differential operator) dominates the diffusion (second order differential operator).
• In the convection-dominated regime, the Galerkin method as presented in this class does not work. One needs new ideas for discretizations and these new discretizations create new challenges for the numerical analysis.
Problems with more than one unknown function. The fundamental equation of fluid dynamics, the Navier-Stokes equations, Section 1.3, belong to this class.
• It will turn out that the discretization of the Navier-Stokes equations requires special care in the choice of the finite element spaces. The numerical analysis becomes rather involved.
□
References
Adams, R. A. (1975) Sobolev spaces. Academic Press [A subsidiary of Harcourt Brace Jovanovich, Publishers], New York-London, pp. xviii+268. Pure and Applied Mathematics, Vol. 65.
Adams, R. A. & Fournier, J. J. F. (2003) Sobolev spaces. Pure and Applied Mathematics (Amsterdam), vol. 140, second edn. Elsevier/Academic Press, Amsterdam, pp. xiv+305.
Alt, H. (1999) Lineare Funktionalanalysis. Eine anwendungsorientierte Einführung, 3. edn. Springer Berlin.
BRAESS, D. (2001) Finite elements, second edn. Cambridge: Cambridge University Press, pp. xviii+352. Theory, fast solvers, and applications in solid mechanics, Translated from the 1992 German edition by Larry L. Schumaker.
BRENNER, S. C. & scott, L. R. (2008) The mathematical theory of finite element methods. Texts in Applied Mathematics, vol. 15, third edn. New York: Springer, pp. xviii+397.
ClARLET, P. G. (1978) The finite element method for elliptic problems. Amsterdam: North-Holland Publishing Co., pp. xix+530. Studies in Mathematics and its Applications, Vol. 4.
ClARLET, P. G. (2002) The finite element method for elliptic problems. Classics in Applied Mathematics, vol. 40. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM), pp. xxviii+530. Reprint of the 1978 original [North-Holland, Amsterdam; MR0520174 (58 #25001)].
CROUZEIX, M. & RAVIART, P.-A. (1973) Conforming and nonconforming finite element methods for solving the stationary Stokes equations. I. Rev. Frangaise Automat. Informat. Recherche Opérationnelle Sér. Rouge, 7, 33-75.
DEUFLHARD, P. & WEISER, M. (2012) Adaptive numerical solution of PDEs. de Gruyter Textbook. Berlin: Walter de Gruyter & Co., pp. xii+421.
DziUK, G. (2010) Theorie und Numerik partieller Differentialgleichungen. Walter de Gruyter GmbH & Co. KG, Berlin, pp. x+319.
121
122
References
Ern, A. & GUERMOND, j.-L. (2004) Theory and practice of finite elements. Applied Mathematical Sciences, vol. 159. New York: Springer-Verlag, pp. xiv+524.
EVANS, L. C. (2010) Partial differential equations. Graduate Studies in Mathematics, vol. 19, second edn. Providence, RI: American Mathematical Society, pp. xxii+749.
FEFFERMAN, C. (2000). http://www.claymath.org/millennium/Navier-Stokes_Equations/.
galdi, G. P. (2011) An introduction to the mathematical theory of the Navier-Stokes equations. Springer Monographs in Mathematics, second edn. Springer, New York, pp. xiv+1018. Steady-state problems.
GOERING, H., HANS-GÖRG, R. & TOBISKA, L. (2010) Die Finite-Elemente-Methode fr Anfänger, fourth edn. Wiley-VCH, Berlin, pp. ix + 219.
GRISVARD, P. (1985) Elliptic problems in nonsmooth domains. Monographs and Studies in Mathematics, vol. 24. Boston, MA: Pitman (Advanced Publishing Program), pp. xiv+410.
GROSSMANN, C. & ROOS, H.-G. (2007) Numerical treatment of partial differential equations. Universitext. Berlin: Springer, pp. xii+591. Translated and revised from the 3rd (2005) German edition by Martin Stynes.
HAROSKE, D. D. & TRIEBEL, H. (2008) Distributions, Sobolev spaces, elliptic equations. EMS Textbooks in Mathematics. European Mathematical Society (EMS), Zürich, pp. x+294.
john, V. (2016) Finite element methods for incompressible flow problems. Springer Series in Computational Mathematics, vol. 51. Springer, Cham, pp. xiii+812.
John, V. & Matthies, G. (2004) MooNMD—a program package based on mapped finite element methods. Comput. Vis. Sei., 6, 163-169.
LANDAU, L. & LifSCHITZ, E. (1966) Lehrbuch der theoretischen Physik, vol. VI, Hydrodynamik. Akademie-Verlag Berlin.
LeVeQUE, R. j. (2007) Finite difference methods for ordinary and partial differential equations. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM), pp. xvi+341. Steady-state and time-dependent problems.
RANNACHER, R. & TUREK, S. (1992) Simple nonconforming quadrilateral Stokes element. Numer. Methods Partial Differential Equations, 8, 97-111.
samarskii, A. A. (2001) The theory of difference schemes. Monographs and Textbooks in Pure and Applied Mathematics, vol. 240. New York: Marcel Dekker Inc., pp. xviii+761.
SamARSKIJ, A. (1984) Theorie der Differenzenverfahren. Mathematik und ihre Anwendungen in Physik und Technik, vol. 40. Akademische Verlagsgesellschaft Geest & Portig K.-G, Leipzig.
schieweck, F. (1997) Parallele Lösung der stationären inkompressiblen Navier-Stokes Gleichungen. Otto-von-Guericke-Universität Magdeburg, Fakultät für Mathematik. Habilitation.
References
123
SMIRNOW, W. I. (1967) Lehrgang der höheren Mathematik. Teil V. VEB Deutscher Verlag der Wissenschaften, Berlin, pp. xiii+570. Zweite, berichtigte Auflage, Übersetzung aus dem Russischen von Renate Helle und Brigitte Mai, Hochschulbücher für Mathematik, Band 6.
SOLIN, P. (2006) Partial differential equations and the finite element method. Pure and Applied Mathematics (New York). Hoboken, NJ: Wiley-Interscience [John Wiley & Sons], pp. xviii+472.
STRANG, G. & Fix, G. (2008) An analysis of the finite element method, second edn. Wellesley, MA: Wellesley-Cambridge Press, pp. x+402.
WlLBRANDT,   U.,   bartsch,   C,   ahmed,   N.,   ALIA,   N.,   ANKER, F.,
Blank, L., Caiazzo, A., Ganesan, S., Giere, S., Matthies, G., Meesala, R., Shamim, A., Venkatesan, J. & John, V. (2017) ParMooN—A modernized program package based on mapped finite elements. Corn-put. Math. Appl., 74, 74-88. wladimirow, W. S. (1972) Gleichungen der mathematischen Physik. VEB Deutscher Verlag der Wissenschaften, Berlin.