15.1
A One-Dimensional Problem: Bending of a Beam
Consider a beam of unit length supported at its ends in 0 and 1, stretched along its axis by
a force P , and subjected to a transverse load f (x)dx per element dx, as illustrated in Figure
15.1.
0
1
dx
−P
P
f (x)dx
Figure 15.1: Vertical deflection of a beam
The bending moment u(x) at the absissa x is the solution of a boundary problem (BP)
of the form
−u (x) + c(x)u(x) = f(x), 0 < x < 1
u(0) = α
u(1) = β,
401
402
CHAPTER 15. INTRODUCTION TO THE FINITE ELEMENTS METHOD
where c(x) = P/(EI(x)), where E is the Young’s modulus of the material of which the beam
is made and I(x) is the principal moment of inertia of the cross-section of the beam at the
abcissa x, and with α = β = 0. For this problem, we may assume that c(x) ≥ 0 for all
x ∈ [0, 1].
Remark: The vertical deflection w(x) of the beam and the bending moment u(x) are related
by the equation
d2w
u(x) = −EI
.
dx2
If we seek a solution u ∈ C2([0, 1]), that is, a function whose first and second derivatives
exist and are continuous, then it can be shown that the problem has a unique solution
(assuming c and f to be continuous functions on [0, 1]).
Except in very rare situations, this problem has no closed-form solution, so we are led to
seek approximations of the solutions.
One one way to proceed is to use the finite difference method , where we discretize the
problem and replace derivatives by differences. Another way is to use a variational approach.
In this approach, we follow a somewhat surprising path in which we come up with a so-called
“weak formulation” of the problem, by using a trick based on integrating by parts!
First, let us observe that we can always assume that α = β = 0, by looking for a solution
of the form u(x) − (α(1 − x) + βx). This turns out to be crucial when we integrate by parts.
There are a lot of subtle mathematical details involved to make what follows rigorous, but
we here, we will take a “relaxed” approach.
First, we need to specify the space of “weak solutions.” This will be the vector space V of
continuous functions f on [0, 1], with f (0) = f (1) = 0, and which are piecewise continuously
differentiable on [0, 1]. This means that there is a finite number of points x0, . . . , xN+1 with
x0 = 0 and xN+1 = 1, such that f (xi) is undefined for i = 1, . . . , N, but otherwise f is
defined and continuous on each interval (xi, xi+1) for i = 0, . . . , N.1 The space V becomes a
Euclidean vector space under the inner product
1
f, g V =
(f (x)g(x) + f (x)g (x))dx,
0
for all f, g ∈ V . The associated norm is
1
1/2
f
=
(f (x)2 + f (x)2)dx
.
V
0
Assume that u is a solution of our original boundary problem (BP), so that
−u (x) + c(x)u(x) = f(x), 0 < x < 1
u(0) = 0
u(1) = 0.
1We also assume that f (x) has a limit when x tends to a boundary of (xi, xi+1).
15.1. A ONE-DIMENSIONAL PROBLEM: BENDING OF A BEAM
403
Multiply the differential equation by any arbitrary test function v ∈ V , obtaining
−u (x)v(x) + c(x)u(x)v(x) = f(x)v(x),
(∗)
and integrate this equation! We get
1
1
1
−
u (x)v(x)dx +
c(x)u(x)v(x)dx =
f (x)v(x)dx.
(†)
0
0
0
Now, the trick is to use integration by parts on the first term. Recall that
(u v) = u v + u v ,
and to be careful about discontinuities, write
1
N
xi+1
u (x)v(x)dx =
u (x)v(x)dx.
0
i=0
xi
Using integration by parts, we have
xi+1
xi+1
xi+1
u (x)v(x)dx =
(u (x)v(x)) dx −
u (x)v (x)dx
xi
xi
xi
xi+1
= [u (x)v(x)]x=xi+1
u (x)v (x)dx
x=x
−
i
xi
xi+1
= u (xi+1)v(xi+1) − u (xi)v(xi) −
u (x)v (x)dx.
xi
It follows that
1
N
xi+1
u (x)v(x)dx =
u (x)v(x)dx
0
i=0
xi
N
xi+1
=
u (xi+1)v(xi+1) − u (xi)v(xi) −
u (x)v (x)dx
i=0
xi
1
= u (1)v(1) − u (0)v(0) −
u (x)v (x)dx.
0
However, the test function v satisfies the boundary conditions v(0) = v(1) = 0 (recall that
v ∈ V ), so we get
1
1
u (x)v(x)dx = −
u (x)v (x)dx.
0
0
Consequently, the equation (†) becomes
1
1
1
u (x)v (x)dx +
c(x)u(x)v(x)dx =
f (x)v(x)dx,
0
0
0
404
CHAPTER 15. INTRODUCTION TO THE FINITE ELEMENTS METHOD
or
1
1
(u v + cuv)dx =
f vdx,
for all v ∈ V.
(∗∗)
0
0
Thus, it is natural to introduce the bilinear form a : V × V → R given by
1
a(u, v) =
(u v + cuv)dx,
for all u, v ∈ V ,
0
and the linear form f : V → R given by
1
f (v) =
f (x)v(x)dx,
for all v ∈ V .
0
Then, (∗∗) becomes
a(u, v) = f (v),
for all v ∈ V.
We also introduce the energy function J given by
1
J(v) = a(v, v) − f(v) v ∈ V.
2
Then, we have the following theorem.
Theorem 15.1. Let u be any solution of the boundary problem (BP).
(1) Then we have
a(u, v) = f (v),
for all v ∈ V,
(WF)
where
1
a(u, v) =
(u v + cuv)dx,
for all u, v ∈ V ,
0
and
1
f (v) =
f (x)v(x)dx,
for all v ∈ V .
0
(2) If c(x) ≥ 0 for all x ∈ [0, 1], then a function u ∈ V is a solution of (WF) iff u
minimizes J(v), that is,
J(u) = inf J(v),
v∈V
with
1
J(v) = a(v, v) − f(v) v ∈ V.
2
Furthermore, u is unique.
15.1. A ONE-DIMENSIONAL PROBLEM: BENDING OF A BEAM
405
Proof. We already proved (1).
To prove (2), first we show that
v 2V ≤ 2a(v, v), for all v ∈ V.
For this, it suffices to prove that
1
v 2
(f (x))2dx,
for all v
V ≤ 2
∈ V.
0
However, by Cauchy-Schwarz for functions, for every x ∈ [0, 1], we have
x
1
1
1/2
|v(x)| =
v (t)dt ≤
|v (t)|dt ≤
|v (t)|2dt
,
0
0
0
and so
1
1
v 2 =
((v(x))2 + (v (x))2)dx
(v (x))2dx
V
≤ 2
≤ 2a(v, v),
0
0
since
1
a(v, v) =
((v )2 + cv2)dx.
0
Next, it is easy to check that
1
J(u + v) − J(u) = a(u, v) − f(v) + a(v, v), for all u, v ∈ V .
2
Then, if u is a solution of (WF), we deduce that
1
1
J(u + v) − J(u) = a(v, v) ≥
v
2
4
V ≥ 0
for all v ∈ V.
since a(u, v) − f(v) = 0 for all v ∈ V . Therefore, J achieves a minimun for u.
We also have
θ2
J(u + θv) − J(u) = θ(a(u, v) − f(v)) +
a(v, v) for all θ ∈
2
R,
and so J(u + θv) − J(u) ≥ 0 for all θ ∈ R. Consequently, if J achieves a minimum for u,
then a(u, v) = f (v), which means that u is a solution of (WF).
Finally, assuming that c(x) ≥ 0, we claim that if v ∈ V and v = 0, then a(v, v) > 0. This
is because if a(v, v) = 0, since
v 2V ≤ 2a(v, v) for all v ∈ V,
we would have v
= 0, that is, v = 0. Then, if v = 0, from
V
1
J(u + v) − J(u) = a(v, v) for all v ∈ V
2
we see that J(u + v) > J(u), so the minimum u is unique
406
CHAPTER 15. INTRODUCTION TO THE FINITE ELEMENTS METHOD
Theorem 15.1 shows that every solution u of our boundary problem (BP) is a solution
(in fact, unique) of the equation (WF).
The equation (WF) is called the weak form or variational equation associated with the
boundary problem. This idea to derive these equations is due to Ritz and Galerkin.
Now, the natural question is whether the variational equation (WF) has a solution, and
whether this solution, if it exists, is also a solution of the boundary problem (it must belong
to C2([0, 1]), which is far from obvious). Then, (BP) and (WF) would be equivalent.
Some fancy tools of analysis can be used to prove these assertions. The first difficulty is
that the vector space V is not the right space of solutions, because in order for the variational
problem to have a solution, it must be complete. So, we must construct a completion of the
vector space V . This can be done and we get the Sobolev space H10(0, 1). Then, the question
of the regularity of the “weak solution” can also be tackled.
We will not worry about all this. Instead, let us find approximations of the problem (WF).
Instead of using the infinite-dimensional vector space V , we consider finite-dimensional sub-
spaces Va (with dim(Va) = n) of V , and we consider the discrete problem:
Find a function u(a) ∈ Va, such that
a(u(a), v) = f (v),
for all v ∈ Va.
(DWF)
Since Va is finite dimensional (of dimension n), let us pick a basis of functions (w1, . . . , wn)
in Va, so that every function u ∈ Va can we written as
u = u1w1 + · · · + unwn.
Then, the equation (DWF) holds iff
a(u, wj) = f (wj),
j = 1, . . . , n,
and by plugging u1w1 + · · · + unwn for u, we get a system of k linear equations
n
a(wi, wj)ui = f (wj),
1 ≤ j ≤ n.
i=1
Because a(v, v) ≥ 1 v
, the bilinear form a is symmetric positive definite, and thus
2
Va
the matrix (a(wi, wj)) is symmetric positive definite, and thus invertible. Therefore, (DWF)
has a solution given by a linear system!
From a practical point of view, we have to compute the integrals
1
aij = a(wi, wj) =
(wiwj + cwiwj)dx,
0
and
1
bj = f (wj) =
f (x)wj(x)dx.
0
15.1. A ONE-DIMENSIONAL PROBLEM: BENDING OF A BEAM
407
However, if the basis functions are simple enough, this can be done “by hand.” Otherwise,
numerical integration methods must be used, but there are some good ones.
Let us also remark that the proof of Theorem 15.1 also shows that the unique solution of
(DWF) is the unique minimizer of J over all functions in Va. It is also possible to compare
the approximate solution u(a) ∈ Va with the exact solution u ∈ V .
Theorem 15.2. Suppose c(x) ≥ 0 for all x ∈ [0, 1]. For every finite-dimensional subspace
Va (dim(Va) = n) of V , for every basis (w1, . . . , wn) of Va, the following properties hold:
(1) There is a unique function u(a) ∈ Va such that
a(u(a), v) = f (v),
for all v ∈ Va,
(DWF)
and if u(a) = u1w1 + · · · + unwn, then u = (u1, . . . , un) is the solution of the linear
system
Au = b,
(∗)
with A = (aij) = (a(wi, wj)) and bj = f (wj), 1 ≤ i, j ≤ n. Furthermore, the matrix
A = (aij) is symmetric positive definite.
(2) The unique solution u(a) ∈ Va of (DWF) is the unique minimizer of J over Va, that is,
J(u(a)) = inf J(v),
v∈Va
(3) There is a constant C independent of Va and of the unique solution u ∈ V of (WF),
such that
u − u(a)
≤ C inf u − v
.
V
v∈V
V
a
We proved (1) and (2), but we will omit the proof of (3) which can be found in Ciarlet
[22].
Let us now give examples of the subspaces Va used in practice. They usually consist of
piecewise polynomial functions.
Pick an integer N ≥ 1 and subdivide [0, 1] into N + 1 intervals [xi, xi+1], where
1
xi = hi,
h =
,
i = 0, . . . , N + 1.
N + 1
We will use the following fact: every polynomial P (x) of degree 2m + 1 (m ≥ 0) is
completely determined by its values as well as the values of its first m derivatives at two
distinct points α, β ∈ R.
408
CHAPTER 15. INTRODUCTION TO THE FINITE ELEMENTS METHOD
There are various ways to prove this. One way is to use the Bernstein basis, because
the kth derivative of a polynomial is given by a formula in terms of its control points. For
example, for m = 1, every degree 3 polynomial can be written as
P (x) = (1 − x)3 b0 + 3(1 − x)2x b1 + 3(1 − x)x2 b2 + x3 b3,
with b0, b1, b2, b3 ∈ R, and we showed that
P (0) = 3(b1 − b0)
P (1) = 3(b3 − b2).
Given P (0) and P (1), we determine b0 and b3, and from P (0) and P (1), we determine b1
and b2.
In general, for a polynomial of degree m written as
m
P (x) =
bjBm
j (x)
j=0
in terms of the Bernstein basis (Bm
0 (x), . . . , Bm
m (x)) with
m
Bm
j (x) =
(1 − x)m−jxj,
j
it can be shown that the kth derivative of P at zero is given by
k
k
P (k)(0) = m(m − 1) · · · (m − k + 1)
(−1)k−i b
i
i
,
i=0
and there is a similar formula for P (k)(1).
Actually, we need to use the Bernstein basis of polynomials Bm[r, s], where
k
m
s − x m−j x − r j
Bm
j [r, s](x) =
,
j
s − r
s − r
with r < s, in which case
m(m − 1) · · · (m − k + 1)
k
k
P (k)(0) =
(−1)k−i b
(s − r)k
i
i
,
i=0
with a similar formula for P (k)(1). In our case, we set r = xi, s = xi+1.
Now, if the 2m + 2 values
P (0), P (1)(0), . . . , P (m)(0), P (1), P (1)(1), . . . , P (m)(1)
15.1. A ONE-DIMENSIONAL PROBLEM: BENDING OF A BEAM
409
are given, we obtain a triangular system that determines uniquely the 2m + 2 control points
b0, . . . , b2m+1.
Recall that Cm([0, 1]) denotes the set of Cm functions f on [0, 1], which means that
f, f (1), . . . , f (m) exist are are continuous on [0, 1].
We define the vector space V m as the subspace of Cm([0, 1]) consisting of all functions f
N
such that
1. f (0) = f (1) = 0.
2. The restriction of f to [xi, xi+1] is a polynomial of degree 2m + 1, for i = 0, . . . , N.
Observe that the functions in V 0 are the piecewise affine functions f with f (0) = f (1) =
N
0; an example is shown in Figure 15.2.
y
x
0
ih
1
Figure 15.2: A piecewise affine function
This space has dimension N , and a basis consists of the “hat functions” wi, where the
only two nonflat parts of the graph of wi are the line segments from (xi−1, 0) to (xi, 1), and
from (xi, 1) to (xi+1, 0), for i = 1, . . . , N, see Figure 15.3.
The basis functions wi have a small support, which is good because in computing the
integrals giving a(wi, wj), we find that we get a tridiagonal matrix. They also have the nice
property that every function v ∈ V 0 has the following expression on the basis (w
N
i):
N
v(x) =
v(ih)wi(x),
x ∈ [0, 1].
i=1
410
CHAPTER 15. INTRODUCTION TO THE FINITE ELEMENTS METHOD
y
wi
x
(i − 1)h ih (i + 1)h
Figure 15.3: A basis “hat function”
In general, it it not hard to see that V m has dimension mN + 2(m
N
− 1).
Going back to our problem (the bending of a beam), assuming that c and f are constant
functions, it is not hard to show that the linear system (∗) becomes
2 + 2ch2
−1 + c h2
u
f
3
6
1
−1 + c h2 2 + 2ch2 −1 + c h2
u
f
6
3
6
2
1
. .
..
..
.
. ..
. ..
= h
.
h
.
.
−1 + c h2 2 + 2ch2 −1 + c h2
u
f
6
3
6
N −1
−1 + c h2 2 + 2ch2
u
f
6
3
N
We can also find a basis of 2N + 2 cubic functions for V 1 consisting of functions with
N
small support. This basis consists of the N functions w0i and of the N + 2 functions w1i
15.1. A ONE-DIMENSIONAL PROBLEM: BENDING OF A BEAM
411
uniquely determined by the following conditions:
w0i(xj) = δij, 1 ≤ j ≤ N
(w0i) (xj) = 0, 0 ≤ j ≤ N + 1
w1i(xj) = 0, 1 ≤ j ≤ N
(w1i) (xj) = δij, 0 ≤ j ≤ N + 1,
with δij = 1 iff i = j and δij = 0 if i = j. Some of these functions are displayed in Figure
15.4. The function w0i is given explicitly by
1
w0i(x) =
(x − (i − 1)h)2((2i + 1)h − 2x), (i − 1)h ≤ x ≤ ih,
h3
1
w0i(x) =
((i + 1)h − x)2(2x − (2i − 1)h), ih ≤ x ≤ (i + 1)h,
h3
for i = 1, . . . , N . The function w1j is given explicitly by
1
w1j(x) = − (ih − x)(x − (i − 1)h)2, (i − 1)h ≤ x ≤ ih,
h2
and
1
w1j(x) =
((i + 1)h − x)2(x − ih), ih ≤ x ≤ (i + 1)h,
h2
for j = 0, . . . , N + 1. Furthermore, for every function v ∈ V 1, we have
N
N
N +1
v(x) =
v(ih)w0i(x) +
v jih)w1j(x), x ∈ [0, 1].
i=1
j=0
If we order these basis functions as
w10, w01, w11, w02, w12, . . . , w0N, w1N, w1N+1,
we find that if c = 0, the matrix A of the system (∗) is tridiagonal by blocks, where the blocks
are 2 × 2, 2 × 1, or 1 × 2 matrices, and with single entries in the top left and bottom right
corner. A different order of the basis vectors would mess up the tridiagonal block structure
of A. We leave the details as an exercise.
Let us now take a quick look at a two-dimensional problem, the bending of an elastic
membrane.
412
CHAPTER 15. INTRODUCTION TO THE FINITE ELEMENTS METHOD
y
w0i
w10
w1j
w1N+1 x
0
ih
jh
1
Figure 15.4: The basis functions w0i and w1j
15.2
A Two-Dimensional Problem: An Elastic
Membrane
Consider an elastic membrane attached to a round contour whose projection on the (x1, x2)-
plane is the boundary Γ of an open, connected, bounded region Ω in the (x1, x2)-plane, as
illustrated in Figure 15.5. In other words, we view the membrane as a surface consisting of
the set of points (x, z) given by an equation of the form
z = u(x),
with x = (x1, x2) ∈ Ω, where u: Ω → R is some sufficiently regular function, and we think
of u(x) as the vertical displacement of this membrane.
We assume that this membrane is under the action of a vertical force τ f (x)dx per surface
element in the horizontal plane (where τ is the tension of the membrane). The problem is
15.2. A TWO-DIMENSIONAL PROBLEM: AN ELASTIC MEMBRANE
413
τ f (x)dx
g(y)
u(x)
x2
Ω
dx
x
y
x
Γ
1
Figure 15.5: An elastic membrane
to find the vertical displacement u as a function of x, for x ∈ Ω. It can be shown (under
some assumptions on Ω, Γ, and f ), that u(x) is given by a PDE with boundary condition,
of the form
−∆u(x) = f(x), x ∈ Ω
u(x) = g(x),
x ∈ Γ,
where g : Γ → R represents the height of the contour of the membrane. We are looking for
a function u in C2(Ω) ∩ C1(Ω). The operator ∆ is the Laplacian, and it is given by
∂2u
∂2u
∆u(x) =
(x) +
(x).
∂x21
∂x22
This is an example of a boundary problem, since the solution u of the PDE must satisfy the
co