11 Chapter Numerical Methods: Finite Elements
In Chapter 5, we introduced the first, the oldest, and in many ways the simplest class of numerical algorithms for approximating the solutions to partial differential equations: finite differences. In the present chapter, we study the second of the two major numerical paradigms: the finite element method. Finite elements are of more recent vintage, having first appeared soon after the Second World War; historical details can be found in [134]. As a consequence of their flexibility and ability to readily handle complicated geometries, finite elements have, in many situations, become the basic method of choice for solving equilibrium boundary value problems governed by elliptic partial differential equations. Extensions to dynamical problems are also extensively utilized, but these will not be treated here due to lack of space. Finite elements rely on a more sophisticated understanding of the partial differential equation, in that, unlike finite differences, they are not obtained by simply replacing derivatives by their numerical approximations. Rather, they can be based on an associated minimization principle which, as we learned in Chapter 10, can be used to characterize the unique solution to a positive definite boundary value problem. Restricting the minimizing functional to an appropriately chosen finite-dimensional subspace of functions reduces it to a finite-dimensional minimization problem, that can then be solved by numerical linear algebra. When properly formulated, the resulting finite-dimensional minimization problem has a solution that well approximates the true minimizer. After setting up the basic framework, we illustrate the basic constructions in the context of boundary value problems for ordinary differential equations. The following section then extends finite element analysis to boundary value problems associated with the two-dimensional Laplace and Poisson equations. A rigorous justification and proof of convergence of the finite element approximations requires further analysis, and we refer the interested reader to [134, 151]. In this chapter, we shall concentrate on understanding how to formulate and implement the finite element method in practical situations.
11.1. Minimization and Finite Elements.
To explain the key ideas behind the finite element method, we return to the abstract framework that was developed in Chapter 10. Recall that we can characterize the solution to a positive definite boundary value problem as the function u U , belonging to a certain infinite-dimensional function space, that minimizes an associated quadratic functional Q: U R. This sets the stage for the first key idea of the finite element method. Instead of trying to minimize the functional Q over the entire infinite-dimensional function space, 1/19/12 394
c 2012 Peter J. Olver
we will seek to minimize it over a finite-dimensional subspace W U . The effect is to reduce a problem in analysis -- a boundary value problem for a differential equation -- to a problem in linear algebra, and hence one that a computer is capable of solving. On the surface, the idea seems crazy: how could one expect to come close to finding the minimizer in a gigantic infinite-dimensional function space by restricting the search to a paltry finite-dimensional subspace? But this is where the magic of infinite dimensions comes into play. One can, in fact, approximate all (reasonable) functions arbitrarily closely by functions belonging to finite-dimensional subspaces. Indeed, you are already familiar with two examples: Fourier series, where one approximates rather general periodic functions by trigonometric polynomials, and interpolation theory, in which one approximates functions by ordinary polynomials, or, more sophisticatedly, by splines. Thus, the finite element idea perhaps is not as outlandish as it might initially seem. To be a bit more explicit, let us begin with a linear operator L: U V between real inner product spaces, where, as in Section 10.1, u ; u is used to denote the inner product in U , and v ; v the inner product in V . To ensure uniqueness of solutions, we always assume that L has trivial kernel: ker L = {0}. According to Theorem 10.25, the element u U that minimizes the quadratic function(al) Q[ u ] = where
1 2
L[ u ]
2
- f ;u ,
(11.1)
denotes the norm in V , is the solution to the linear system S[ u ] = f, where S = L L, (11.2)
with L : V U denoting the adjoint operator. The hypothesis that L have trivial kernel implies that S is a self-adjoint, positive definite linear operator, and so the solution to (11.2) is unique. In our applications, L is a linear differential operator between function spaces, e.g., the gradient, while Q[ u ] represents a quadratic functional, e.g., the Dirichlet principle, and the associated linear system (11.2) forms a positive definite boundary value problem, e.g., the Poisson equation along with suitable self-adjoint boundary conditions. To form a finite element approximation to the solution u U , rather than try to minimize Q[ u ] on the entire function space U , we now seek to minimize it on a suitably chosen finite-dimensional subspace W U . We will specify W by selecting a set of linearly independent functions 1 , . . . , n U , and letting W be their span. Thus, 1 , . . . , n form a basis of W , whereby dim W = n, and the general element of W is a (uniquely determined) linear combination w(x) = c1 1 (x) + + cn n (x) (11.3) of the basis functions. Our goal is to minimize Q[ w ] over all possible w W ; in other words, we need to determine the coefficients c1 , . . . , cn such that Q[ w ] = Q[ c1 1 + + cn n ] (11.4)
is as small as possible. Substituting (11.3) back into (11.1), and then expanding, using the linearity of L and then the bilinearity of the inner product, the resulting expression is the 1/19/12 395
c 2012 Peter J. Olver
quadratic function 1 P (c) = 2
n n
kij ci cj -
i,j = 1 i=1
bi ci =
1 2
cT K c - cT b,
(11.5)
in which T c = ( c1 , c2 , . . . , cn ) R n is the vector of unknown coefficients in (11.3); K = (kij ) is the symmetric n n matrix with entries kij =
T
L[ i ] ; L[ j ] ,
i, j = 1, . . . , n;
(11.6)
b = ( b1 , b2 , . . . , bn ) is the vector with entries bi = f ; i , i = 1, . . . , n. (11.7)
Note that formula (11.6) uses the inner product on the target space V , whereas (11.7) relies on the inner product on the domain space U . Thus, once we've specified the basis functions i , the coefficients kij and bi are all known quantities. We have effectively reduced our original problem to the finitedimensional problem of minimizing the quadratic function (11.5) over all possible vectors c R n . The symmetric matrix K is, in fact, positive definite, since, by the preceding computation,
n
cT K c =
i,j = 1
kij ci cj = L[ c1 1 (x) + + cn n ]
2
= L[ w ]
2
> 0
(11.8)
as long as L[ w ] = 0. Moreover, our initial assumption tells us that L[ w ] = 0 if and only if w = 0, and hence, by linear independence, (11.8) is indeed positive for all c = 0. We can now invoke the finite-dimensional minimization result contained in Example 10.24 to conclude that the unique minimizer to (11.5) is obtained by solving the associated linear system K c = b, and hence c = K -1 b. (11.9) Remark : When of moderate size, the linear system (11.9) can be solved by basic Gaussian Elimination. When the size (i.e., the dimension, n, of the subspace W ) becomes too large, as is often the case in dealing with partial differential equations, it is better to rely on an iterative linear system solver, e.g., the GaussSeidel method or Successive Over-Relaxation (SOR); see [108] for details. This summarizes the basic abstract setting for the finite element method. The key issue, then, is how to effectively choose the finite-dimensional subspace W . Two candidates that might spring to mind are the space of polynomials of degree n, or the space of trigonometric polynomials (truncated Fourier series) of degree n. However, for a variety of reasons, neither is well suited to the finite element method. One constraint is that the functions in W must satisfy the relevant boundary conditions -- otherwise W would not be a subspace of U . More importantly, in order to obtain sufficient accuracy of the approximate solution, the linear algebraic system (11.9) will typically -- especially when 1/19/12 396
c 2012 Peter J. Olver
0.8 0.6 0.4 0.2 0.2 0.4 0.6 0.8 1
Figure 11.1.
A Continuous Piecewise Affine Function.
dealing with partial differential equations -- be quite large, and hence it is desirable that the coefficient matrix K be as sparse as possible, i.e., have lots of zero entries. Otherwise, computing the solution may well be too time-consuming to be of much practical value. With this in mind, the second innovative contribution of the finite element method is to first (paradoxically) enlarge the space U of allowable functions upon which to minimize the quadratic functional Q[ u ]. The governing differential equation requires its (classical) solutions to have a certain degree of smoothness, whereas the associated minimization principle typically requires only half as many continuous derivatives. Thus, for second order boundary value problems, the differential equation requires continuous second order derivatives, while the quadratic functional Q[ u ] only involves first order derivatives. It fact, it can be rigorously shown that, under rather mild hypotheses, the functional retains the same minimizing solution, even when one allows functions that fail to qualify as classical solutions to the differential equation.
11.2. Finite Elements for Ordinary Differential Equations.
To make the preceding discussion concrete, let us focus our attention on a boundary value problem governed by a second order ordinary differential equation. For example, we might consider a SturmLiouville problem (10.69) with homogeneous Dirichlet boundary conditions. Once we understand how the finite element constructions work in this simplest context, we will be in a good position to develop finite element approximations for boundary value problems governed by elliptic partial differential equations. For such one-dimensional boundary value problems, a popular and effective choice of the finite-dimensional subspace W is to use continuous, piecewise affine functions. Recall that a function is affine if and only if its graph is a straight line: f (x) = a x + b. (The function is linear , in accordance with the general Definition B.32, if and only if b = 0.) A function is called piecewise affine if its graph consists of a finite number of straight line segments; a typical example is plotted in Figure 11.1. Continuity requires that the individual segments be connected together end to end. Given a boundary value problem on a bounded interval [ a, b ], let us fix a finite collection of nodes a = x0 < x1 < x2 < < xn-1 < xn = b. 1/19/12 397
c 2012 Peter J. Olver
The formulas simplify if one uses equally spaced nodes, but this is not necessary for the method to be carried out. Let W denote the vector space consisting of all continuous, piecewise affine functions w(x) on the interval a x b that satisfy the homogeneous boundary conditions and are affine when restricted to each subinterval [ xj , xj+1 ]. Thus, in the case of homogeneous Dirichlet boundary conditions, we require that w(a) = w(b) = 0. On each subinterval, w(x) = cj + bj (x - xj ), for xj x xj+1 , j = 0, . . . , n - 1, (11.10)
for certain constants cj , bj . Continuity of w(x) requires cj = w(x+ ) = w(x- ) = cj-1 + bj-1 hj-1 , j j j = 1, . . . , n - 1, (11.11)
where hj-1 = xj -xj-1 denotes the length of the j th subinterval. The boundary conditions (11.10) require w(a) = c0 = 0, w(b) = cn-1 + hn-1 bn-1 = 0. (11.12)
The function w(x) involves a total of 2 n unspecified coefficients c0 , . . . , cn-1 , b0 , . . . , bn-1 . The continuity conditions (11.11) and the second boundary condition (11.12) uniquely determine the bj . The first boundary condition specifies c0 , while the remaining n - 1 coefficients c1 = w(x1 ), . . . , cn-1 = w(xn-1 ) are arbitrary. We conclude that the finite element subspace W has dimension n - 1, which is the number of interior nodes. Remark : Every function w(x) in our subspace has piecewise constant first derivative w (x). However, the jump discontinuities in w (x) imply that its second derivative w (x) may include delta function impulses at the nodes, and hence w(x) is far from being a solution to the differential equation. Nevertheless, in practice, the finite element minimizer w W will (under suitable assumptions) provide a reasonable approximation to the actual solution u . The most convenient basis for W consists of the hat functions, which are continuous, piecewise affine functions satisfying j (xk ) = 1, 0, j = k, j = k, for j = 1, . . . , n - 1, k = 0, . . . , n. (11.13)
The graph of a typical hat function appears in Figure 11.2. The explicit formula is easily established: x-x j-1 , xj-1 x xj , x -x j j-1 xj+1 - x j = 1, . . . , n - 1. (11.14) j (x) = , xj x xj+1 , x j+1 - xj 0, x xj-1 or x xj+1 , 1/19/12 398
c 2012 Peter J. Olver
1.2 1 0.8 0.6 0.4 0.2 1 -0.2 2 3 4 5 6 7
Figure 11.2.
A Hat Function.
One advantage of using these basis functions is that, thanks to (11.13), the coefficients in the linear combination (11.3) coincide with the value of the sum at the nodes: cj = w(xj ), j = 1, . . . , n. (11.15)
An additional benefit is that the resulting coefficient matrix (11.6) turns out to be tridiagonal, and hence the tridiagonal Gaussian Elimination algorithm, [108], will rapidly produce the solution to the associated linear system (11.9). Since the accuracy of the finite element solution increases with the number of nodes, this numerical scheme allows us to easily compute very accurate approximations to the solution to the boundary value problem. The sparse tridiagonal nature of the finite element matrix is a consequence of the fact that the basis functions are zero on much of the interval, or, in more mathematical language, that they have small support, in the following sense. Definition 11.1. The support of a function f (x), written supp f , is the closure of the set where f (x) = 0. Thus, a point x will belong to the support provided f is not zero there, or at least is not zero at nearby points. For example, the support of the hat function (11.14) is the (small) interval [ xj-1 , xj+1 ]. Example 11.2. Let (x) > 0 for 0 x . Consider the equilibrium equations S[ u ] = - d dx (x) du dx = f (x), 0 < x < , u(0) = u() = 0,
for a non-uniform bar with stiffness and fixed ends, when subject to external forcing f . In order to formulate a finite element approximation scheme, we begin with the minimization principle based on the quadratic functional
Q[ u ] =
0
1 2
(x) u (x)2 - f (x) u(x) dx,
which is a special case of the general SturmLiouville minimization principle (10.72). 1/19/12 399
c 2012 Peter J. Olver
We divide the interval [ 0, ] into n equal subintervals, each of length h = /n. The resulting uniform mesh has nodes j , j = 0, . . . , n. n The corresponding finite element basis hat functions are explicitly given by xj-1 x xj , (x - xj-1 )/h, j (x) = j = 1, . . . , n - 1. (x - x)/h, xj x xj+1 , j+1 0, otherwise, xj = j h =
(11.16)
The associated linear system (11.9) has coefficient matrix entries kij = ; i j =
0
(x) (x) (x) dx, i j
i, j = 1, . . . , n - 1.
Therefore, due to the small support of the basis hat functions, the finite element coefficient matrix assumes the tridiagonal form s +s -s
0 1 1
Since the function i (x) vanishes except on the interval xi-1 < x < xi+1 , while j (x) vanishes outside xj-1 < x < xj+1 , the integral will vanish unless i = j or i = j 1. Moreover, xj-1 x xj , 1/h, i (x) = j = 1, . . . , n - 1. - 1/h, xj x xj+1 , 0, otherwise,
where
- s1 1 K= 2 h
s1 + s2 - s2
- s2 s2 + s3 - sn-3
- s3
..
.
..
.. . . sn-3 + sn-2 - sn-2 - sn-2 sn-2 + sn-1
,
(11.17)
xj+1
sj =
(x) dx
xj
(11.18)
is the total stiffness of the j th subinterval. The corresponding right hand side has entries
bj = f ; j =
0
f (x) j (x) dx
xj xj-1 xj+1
1 = h
(11.19) (xj+1 - x)f (x) dx .
(x - xj-1 )f (x) dx +
xj
In practice, we do not have to explicitly evaluate the integrals (11.18, 19), but may replace them by a suitably close numerical approximation. When the step size h 1 is small, 1/19/12 400
c 2012 Peter J. Olver
0.08
0.08
0.06
0.06
0.04
0.04
0.02
0.02
0.2
0.4
0.6
0.8
1
0.2
0.4
0.6
0.8
1
0.08
0.08
0.06
0.06
0.04
0.04
0.02
0.02
0.2
0.4
0.6
0.8
1
0.2
0.4
0.6
0.8
1
Figure 11.3.
Finite Element Solution to (11.23).
then the integrals are taken over small intervals, and so the elementary trapezoid rule, [28, 129], produces sufficiently accurate approximations: sj h [ (xj ) + (xj+1 ) ], 2 bj h f (xj ). (11.20)
In this case, the j th equation in the finite element linear system is, upon dividing by h, - cj+1 - 2 cj + cj-1 = f (xj ). h2
The finite element system K c = b is then solved for c, whose entries, according to (11.15), coincide with the values of the finite element approximation to the solution at the nodes: cj = w(xj ) u(xj ). In particular, in the homogeneous case (x) 1, the coefficient matrix (11.17) reduces to the very special form 2 -1 -1 2 -1 -1 2 -1 1 .. .. .. K= (11.21) . . . . h -1 2 -1 -1 2 (11.22)
The left hand side coincides with the standard finite difference approximation (5.5) to minus the second derivative - u (xj ) at the node xj . As a result, in this particular case the finite element and finite difference numerical solution schemes happen to coincide. Example 11.3. Consider the boundary value problem - 1/19/12 du d (x + 1) = 1, dx dx 401 u(0) = 0, u(1) = 0.
c 2012
(11.23)
Peter J. Olver
The explicit solution is easily found by direct integration: u(x) = - x + log(x + 1) . log 2 (11.24)
It minimizes the associated quadratic functional
Q[ u ] =
0
1 2
(x + 1) u (x)2 - u(x) dx
(11.25)
over all C2 functions u(x) that satisfy the given boundary conditions. The finite element system (11.9) has coefficient matrix given by (11.17) and right hand side (11.19), where
xj+1
sj =
xj
(1 + x) dx = h (1 + xj ) +
1 2
h2 = h + h2
j+
1 2
xj+1
,
bj =
1 dx = h.
xj
The resulting piecewise affine approximation to the solution is plotted in Figure 11.3. The first three graphs contain, respectively, 5, 10, 20 nodes, so that h = .2, .1, .05, while the last plots the exact solution (11.24). Even when computed on rather coarse meshes, the finite element approximation is quite respectable. Remark : One can obtain a smoother, and hence more realistic approximation to the solution by smoothly interpolating the finite element approximations cj u(xj ) at the nodes, e.g., by use of cubic splines, [108]. Alternatively, one can require that the finite element functions themselves be smoother, e.g., by making the finite element subspace consist of piecewise cubic splines that satisfy the boundary conditions.
11.3. Finite Elements in Two Dimensions.
The same basic strategy underlies the development of finite element methods for approximating the solution to boundary value problems governed by elliptic partial differential equations. In this section, we concentrate on the simplest case: the two-dimensional Poisson equation. Having mastered this, the reader will be well-equipped to carry over the method to more general partial differential equations and higher dimensional problems. As before, we concentrate on the practical design of the finite element procedure, and refer the reader to a more advanced text, e.g., [134], for the analytical details and proofs of convergence. Most of the multi-dimensional complications lie not in the underlying theory, but rather in the arena of data management and organization. For specificity, consider the homogeneous Dirichlet boundary value problem - u = f in , u=0 on , (11.26)
on a bounded domain R 2 . According to Theorem 10.30, the solution u = u is characterized as the unique minimizing function for the Dirichlet functional Q[ u ] =
1 2
u
2
- u;f =
1 2
1 u2 + 2 u2 - f u dx dy x y
(11.27)
among all C2 functions u(x, y) that satisfy the prescribed boundary conditions. 1/19/12 402
c 2012 Peter J. Olver
Figure 11.4.
Triangulation of a Planar Domain.
To construct a finite element approximation, we restrict the Dirichlet functional to a suitably chosen finite-dimensional subspace. As in the one-dimensional situation, the most convenient finite-dimensional subspaces consist of functions that may lack the requisite degree of smoothness that qualifies them as candidate solutions to the partial differential equation. Nevertheless, they will provide good approximations to the actual classical solution. As before, an important practical consideration is to employ functions with small support, meaning that they vanish on most of the domain. This strategy ensures a sparse finite element matrix, the benefit being that the approximate solution to the linear finite element system can be relatively rapidly calculated, usually by application of an iterative numerical scheme. Triangulation The first step is to introduce a mesh consisting of a finite number of nodes xk = (xk , yk ), k = 1, . . . , m, usually lying inside the domain R 2 . Unlike the finite difference approach, we are not restricted to a rectangular mesh, giving the finite element method considerably more flexibility in its discretization of the domain. We regard the nodes as forming the vertices of a triangulation of the domain, consisting of a collection of nonS overlapping small triangles, which we denote by T1 , . . . , TN , whose union, T = T , approximates ; see Figure 11.4 for a typical example. The nodes are split into two categories -- interior nodes and boundary nodes, the latter lying on or close to . A curved boundary will thus be approximated by the polygon given by the outer boundary of the triangulation, T , whose vertices are the boundary nodes. Thus, in any programmed implementation of the finite element method, the first module is a routine that will automatically triangulate a specified domain in some "reasonable" manner, as explained below. As in our one-dimensional finite element construction, the functions w(x, y) in the finite-dimensional subspace W will be continuous and piecewise affine, which means that, on each triangle, the graph of w is a flat plane and hence has the formula w(x, y) = + x + y, for (x, y) T , (11.28)
Here and subsequently, the index is a superscript, not a power.
1/19/12
403
c 2012
Peter J. Olver
Figure 11.5.
Piecewise Affine Function.
for certain constants , , . Continuity of w requires that its values on a common edge between two triangles must agree, and this will impose compatibility constraints on the coefficients , , and , , associated with adjacent pairs of triangles T and T . The full graph of the piecewise affine function z = w(x, y) forms a connected polyhedral surface whose triangular faces lie above the triangles in the domain; see Figure 11.5 for an illustration. In addition, we require that the piecewise affine function w(x, y) vanish at the boundary nodes, which implies that it vanishes on the entire polygonal boundary of the triangulation, T , and hence approximately satisfies the homogeneous Dirichlet boundary conditions on the curved boundary of the original domain, . The next step is to choose a basis of the subspace of piecewise affine functions associated with the given triangulation and subject to the imposed homogeneous Dirichlet boundary conditions. The analog of the one-dimensional hat function (11.13) is the pyramid function k (x, y), which has the value 1 at a single node xk = (xk , yk ), and vanishes at all the other nodes: 1, i = k, k (xi , yi ) = (11.29) 0, i = k. Because, on any triangle, the pyramid function k (x, y) is uniquely determined by its values at the vertices, it will be nonzero only on those triangles which have the node xk as one of their vertices. Hence, as its name implies, the graph of k forms a pyramid of unit height sitting on a flat plane; a typical example appears in Figure 11.6. The pyramid functions k (x, y) associated with the interior nodes xk automatically satisfy the homogeneous Dirichlet boundary conditions on the boundary of the domain -- or, more correctly, on the polygonal boundary of the triangulated domain. Thus, the finite element subspace W is the span of the interior node pyramid functions, and so a general 1/19/12 404
c 2012 Peter J. Olver
Figure 11.6.
Finite Element Pyramid Function.
element w W is a linear combination thereof:
n
w(x, y) =
k=1
ck k (x, y),
(11.30)
where the sum ranges over the n interior nodes of triangulation. the Owing to the original specification (11.29) of the pyramid functions, the coefficients ck = w(xk , yk ) u(xk , yk ), k = 1, . . . , n, (11.31)
are the same as the values of the finite element approximation w(x, y) at the interior nodes. This immediately implies linear independence of the pyramid functions, since the only linear combination that vanishes at all nodes is the trivial one c1 = = cn = 0. Determining the explicit formulae for the pyramid functions is not difficult. On one of the triangles T that has xk as a vertex, k (x, y) will be the unique affine function (11.28) that takes the value 1 at the vertex xk and 0 at its other two vertices xi and xj . Thus, we seek a formula for an affine function or element
k (x, y) = + k x + k y, k
(x, y) T ,
(11.32)
that takes the prescribed values
k (xi , yi ) = k (xj , yj ) = 0, k (xk , yk ) = 1,
or, explicitly,
k (xi , yi ) = + k xi + k yi = 0, k k (xj , yj ) = + k xj + k yj = 0, k k (xk , yk )
(11.33)
=
k
+
k
xk +
k yk
= 1.
An easy exercise in linear algebra, using either Cramer's Rule or direct Gaussian elimination, produces the explicit formulae = k 1/19/12 xi y j - xj y i ,
k =
yi - yj ,
k =
xj - xi ,
c 2012
(11.34)
Peter J. Olver
405
Figure 11.7. where the denominator 1 xi = det 1 xj 1 xk
Vertex Polygons.
is, up to sign, twice the area of the triangle T ; see Exercise .
yi yj = 2 area T yk
(11.35)
Example 11.4. Consider an isoceles right triangle T with vertices x1 = (0, 0), x2 = (1, 0), x3 = (0, 1).
Using (11.3435) (or solving the linear system (11.33) directly), we immediately produce the three corresponding affine elements 1 (x, y) = 1 - x - y, 2 (x, y) = x, 3 (x, y) = y. (11.36)
As required, each k equals 1 at the vertex xk and is zero at the other two vertices. A pyramid function is then obtained by piecing together the individual affine elements: k (x, y) =
k (x, y),
if (x, y) T and xk is a vertex of T , otherwise.
0,
(11.37)
Continuity of k (x, y) is assured since the constituent affine elements have the same values at common vertices, and hence also along common edges. The support of the pyramid function (11.37) is the polygon supp k = Pk =
[
T
(11.38)
consisting of all the triangles T that have the node xk as a vertex. In other words, k (x, y) = 0 whenever (x, y) Pk . We will call Pk the k th vertex polygon. The node xk lies on the interior of its vertex polygon Pk , while the vertices of Pk are all the nodes 1/19/12 406
c 2012 Peter J. Olver
Figure 11.8.
Square Mesh Triangulations.
connected to xk by a single edge of the triangulation. In Figure 11.7, the shaded regions indicate two of the vertex polygons for the triangulation in Figure 11.4. Example 11.5. The simplest, and most common triangulations are based on regular meshes. For example, suppose that the nodes lie on a square grid. If we choose the triangles to all have the same orientation, as in the first picture in Figure 11.8, then the vertex polygons all have the same shape, consisting of 6 triangles of total area 3 h2 -- the shaded region. On the other hand, if we choose an alternating triangulation, as in the second picture, then there are two types of vertex polygons. The first, consisting of four triangles, has area 2 h2 , while the second, containing 8 triangles, has twice the area, 4 h2 . In practice, there are good reasons to prefer the former triangulation. In general, to ensure convergence of the finite element solution to the true minimizer, one should choose a triangulation with the following properties: The individual side lengths of any triangle should be comparable, and so long, skinny triangles and obtuse triangles should be avoided. The areas of nearby triangles T should not vary too much. The areas of nearby vertex polygons Pk should also not vary too much. While the nearby triangles should be of comparable size, one might very well allow wide variations over the entire domain, with small triangles in regions where the solution is changing rapidly, and large triangles in regions of less variability. The Finite Element Equations We now seek to approximate the solution to the homogeneous Dirichlet boundary value problem by restricting the Dirichlet functional (11.27) to the selected finite element subspace W . Using the general framework of Section 11.1, we substitute the formula (11.30) for a general element of W into the quadratic Dirichlet functional (10.79). Expanding, we 1/19/12 407
c 2012 Peter J. Olver
find
n n 2 n
Q[ w ] = Q
i=1
ci i
=
n i=1
ci i
n
- f (x, y)
i=1 1 2
ci i
dx dy
=
1 b c = k c c - 2 i,j = 1 ij i j i = 1 i i
cT K c - bT c.
T
Here, K = (kij ) is an symmetric n n matrix while b = ( b1 , b2 , . . . , bn ) R n , with respective entries kij = i ; j =
is a vector in
i j dx dy, (11.39) f i dx dy,
bi = f ; i =
which also follow directly from the general formulae (11.67). Thus, the finite element approximation (11.30) will minimize the quadratic function P (c) =
1 2
cT K c - b T c
T
(11.40)
over all possible choices of coefficients c = ( c1 , c2 , . . . , cn ) R n , i.e., over all possible function values at the interior nodes. As above, the minimizer is obtained by solving the associated linear system K c = b, (11.41) whose solution can be produced by either Gaussian elimination or a suitable iterative linear systems solver. To find explicit formulae for the matrix coefficients kij in (11.39), we begin by noting that the gradient of the affine element (11.32) is equal to
gk = k (x, y) = k /x k /y
=
k k
=
1
yi - yj , xj - xi
(x, y) T ,
(11.42)
which is a constant vector inside the triangle T , while k = 0 outside. Therefore,
k (x, y) =
gk ,
if (x, y) T which has xk as a vertex, otherwise,
0,
(11.43)
which reduces to a piecewise constant function on the triangulation. Actually, (11.43) is not quite correct since if (x, y) lies on the boundary of a triangle T , then the gradient does not exist, but this will not cause us any difficulty evaluating the ensuing integrals. We will approximate integrals over the domain by summing the corresponding integrals over the individual triangles -- which relies on our assumption that the polygonal boundary of the triangulation is a reasonably close approximation to the true boundary . In particular, kij
T
i j dx dy
kij .
(11.44)
c 2012
1/19/12
408
Peter J. Olver
Figure 11.9.
Right and Equilateral Triangles.
Now, according to (11.43), one or the other gradient in the integrand will vanish on the entire triangle T unless both xi and xj are vertices. Therefore, the only terms contributing to the sum are those triangles T that have both xi and xj as vertices. If i = j there are only two such triangles, having a common edge, while if i = j every triangle in the ith vertex polygon Pi contributes. The individual summands are easily evaluated, since the gradients are constant on the triangles, and so, by (11.43),
kij = gi gj dx dy = gi gj area T = 1 2 gi gj | | .
T
Let T have vertices xi , xj , xk . Then, by (11.42, 43, 35),
kij =
(xi - xk ) (xj - xk ) 1 (yj - yk )(yk - yi ) + (xk - xj )(xi - xk ) | | = - , i = j, 2 2 ( ) 2 | | xj - xk 2 1 (yj - yk )2 + (xk - xj )2 kii = | | = (11.45) 2 ( )2 2 | | (xi - xk ) (xi - xj ) + (xi - xk ) (xj - xk ) = - kij - kik . =- 2
In this manner, each triangle T specifies a collection of 6 different coefficients, kij = kji , indexed by its vertices, and known as the elemental stiffnesses of T . Interestingly, the elemental stiffnesses depend only on the three vertex angles in the triangle and not on its size. Thus, similar triangles have the same elemental stiffnesses. Indeed, according to Exercise , kii = 1 2 cot j + cot k ,
while
1 kij = kji = - 2 cot k ,
i = j,
(11.46)
where 0 < k < denotes the angle in T at the vertex xk .
Example 11.6. The right triangle with vertices x1 = (0, 0), x2 = (1, 0), x3 = (0, 1) has elemental stiffnesses k11 = 1, k22 = k33 =
1 2
,
k12 = k21 = k13 = k31 = - 1 , 2
k23 = k32 = 0.
(11.47)
The same holds for any other isoceles right triangle, as long as we chose the first vertex to be at the right angle. Similarly, an equilateral triangle has all 60 angles, and so its 1/19/12 409
c 2012 Peter J. Olver
Figure 11.10. elemental stiffnesses are k11 = k22 = k33 =
1 3
The Oval Plate.
.5774,
1 k12 = k21 = k13 = k31 = k23 = k32 = - 23 -.2887.
(11.48)
Assembling the Elements The elemental stiffnesses of each triangle will contribute, through the summation (11.44), to the finite element coefficient matrix K. We begin by constructing a larger matrix K, which we call the full finite element matrix , of size m m where m is the total number of nodes in our triangulation, including both interior and boundary nodes. The rows and columns of K are labeled by the nodes xi . Let K = (kij ) be the corresponding mm matrix containing the elemental stiffnesses kij of T in the rows and columns indexed by its vertices, and all other entries equal to 0. Thus, K will have (at most) 9 nonzero entries. The resulting m m matrices are all summed together over all the triangles,
N
K=
=1
K ,
(11.49)
to produce the full finite element matrix, in accordance with (11.44). The full finite element matrix K is too large, since its rows and columns include all the nodes, whereas the finite element matrix K appearing in (11.41) only refers to the n interior nodes. The reduced n n finite element matrix K is simply obtained from K by deleting all rows and columns indexed by boundary nodes, retaining only the elements kij when both xi and xj are interior nodes. For the homogeneous boundary value problem, this is all we require. As we will subsequently see, inhomogeneous boundary conditions are most easily handled by retaining (part of) the full matrix K. The easiest way to digest the construction is by working through a particular example. Example 11.7. A metal plate has the shape of an oval running track, consisting of a rectangle, with side lengths 1 m by 2 m, and two semicircular disks glued onto its shorter ends, as sketched in Figure 11.10. The plate is subject to a heat source while its edges are held at a fixed temperature. The problem is to find the equilibrium temperature distribution within the plate. Mathematically, we must solve the planar Poisson equation, subject to Dirichlet boundary conditions, for the equilibrium temperature u(x, y). 1/19/12 410
c 2012 Peter J. Olver
5 1 2 3 9 8 10 11 14 4 5 6 7 12 13 7 8 6 1
4
13 12
2
3 11
9
10
Triangles Figure 11.11.
Nodes A Coarse Triangulation of the Oval Plate.
Let us describe how to set up the finite element approximation. We begin with a very coarse triangulation of the plate, which will not give particularly accurate results, but does serve to illustrate how to go about assembling the finite element matrix. We divide the rectangular part of the plate into 8 right triangles, while each semicircular end will be approximated by three equilateral triangles. The triangles are numbered from 1 to 14 as indicated in Figure 11.11. There are 13 nodes in all, numbered as in the second figure. Only nodes 1, 2, 3 are interior, while the boundary nodes are labeled 4 through 13, going counterclockwise around the boundary starting at the top. The full finite element matrix K will have size 13 13, its rows and columns labeled by all the nodes, while the reduced matrix K appearing in the finite element equations (11.41) consists of the upper left 3 3 submatrix of K corresponding to the three interior nodes. Each triangle T will contribute the summand K whose values are its elemental stiffnesses, as indexed by its vertices. For example, the first triangle T1 is equilateral, and so has elemental stiffnesses (11.48). Its vertices are labeled 1, 5, and 6, and therefore we place the stiffnesses (11.48) in the rows and columns numbered 1, 5, 6 to form the summand .5774 0 0 0 -.2887 K1 = -.2887 0 0 . . . 0 0 0 0 0 0 0 0 . . . 0 0 0 0 0 0 0 0 . . . 0 -.2887 -.2887 0 0 0 0 0 0 0 0 0 0 0 0 0 0 .5774 -.2887 0 0 -.2887 .5774 0 0 0 0 0 0 0 0 0 . . . . . . . . . . . . 0 0 0 0 0 0 0 0 . . . ... ... ... ... ... ... ... ... .. .
,
where all the undisplayed entries in the full 13 13 matrix are 0. The next triangle T2 has the same equilateral elemental stiffness matrix (11.48), but now its vertices are 1, 6, 7, and 1/19/12 411
c 2012 Peter J. Olver
so it will contribute .5774 0 0 0 0 K2 = -.2887 -.2887 0 . . . 0 0 0 0 0 0 0 0 . . . 0 0 0 0 0 0 0 0 . . . 0 0 0 0 0 0 0 0 . . . 0 -.2887 -.2887 0 0 0 0 0 0 0 0 0 0 0 0 0 .5774 -.2887 0 -.2887 .5774 0 0 0 . . . . . . . . . 0 0 0 0 0 0 0 0 . . . ... ... ... ... ... ... ... ... .. .
.
Similarly for K3 , with vertices 1, 7, 8. On the other hand, triangle T4 is an isoceles right triangle, and so has elemental stiffnesses (11.47). Its vertices are labeled 1, 4, and 5, with vertex 5 at the right angle. Therefore, its contribution is .5 0 0 0 -.5 K4 = 0 0 0 . . . 0 0 0 0 0 0 0 0 . . . 0 0 0 0 0 0 0 .5 0 -.5 0 0 0 0 0 0 . . . . . . -.5 0 0 -.5 1.0 0 0 0 . . . 0 0 0 0 0 0 0 0 . . . 0 0 0 0 0 0 0 0 . . . 0 0 0 0 0 0 0 0 . . . ... ... ... ... ... ... ... ... .. .
.
Continuing in this manner, we assemble 14 contributions K1 , . . . , K14 , each with (at most) 9 nonzero entries. The full finite element matrix is the sum
c K = K1 + K2 + 0 3.732 -1 B -1 4 B B 0 -1 B B B 0 -1 B B -.7887 0 B B B -.5774 0 B = B -.5774 0 B B 0 B -.7887 B B 0 -1 B B 0 0 B B B 0 0 B @ 0 0 0 0
+ K14 0 -1 3.732 0 0 0 0 0 0 -.7887 -.5774 -.5774 -.7887
0 -1 0 2 -.5 0 0 0 0 0 0 0 -.5
-.7887 0 0 -.5 1.577 -.2887 0 0 0 0 0 0 0
-.5774 0 0 0 -.2887 1.155 -.2887 0 0 0 0 0 0
-.5774 0 0 0 0 -.2887 1.155 -.2887 0 0 0 0 0
c 2012
(11.50)
1/19/12
412
Peter J. Olver
Figure 11.12.
A Square Mesh for the Oval Plate.
-.7887 0 0 0 0 0 -.2887 1.577 -.5 0 0 0 0 0 -1 0 0 0 0 0 -.5 2 -.5 0 0 0 0 0 -.7887 0 0 0 0 0 -.5 1.577 -.2887 0 0 0 0 -.5774 0 0 0 0 0 0 -.2887 1.155 -.2887 0 0 0 -.5774 0 0 0 0 0 0 0 -.2887 1.155 -.2887 0 0 C C -.7887 C C C -.5 C C 0 C C C 0 C C 0 C. C 0 C C C 0 C C 0 C C C 0 C C -.2887 A 1.577
1
It is not difficult to directly construct K, bypassing K entirely. For a finer triangulation, the construction is similar, but the matrices become much larger. The procedure can, of course, be automated. Fortunately, if we choose a very regular triangulation, then we do not need to be nearly as meticulous in assembling the stiffness matrices, since many of the entries are the same. The simplest case is when we use a uniform square mesh, and so triangulate the domain into isoceles right triangles. This is accomplished by laying out a relatively dense square grid over the domain R 2 . The interior nodes are the grid points that fall inside the oval domain, while the boundary nodes are all those grid points lying adjacent to one or more of the interior nodes, and are near but not necessarily precisely on the boundary . Figure 11.12 shows the nodes in a square grid with intermesh spacing h = .2. While a bit crude in its approximation of the boundary of the domain, this procedure does have the advantage of making the construction of the associated finite element matrix relatively painless. For such a mesh, all the triangles are isoceles right triangles, with elemental stiffnesses (11.47). Summing the corresponding matrices K over all the triangles, as in (11.49), the 1/19/12 413
c 2012 Peter J. Olver
Since only nodes 1, 2, 3 are interior nodes, the reduced finite element matrix only uses the upper left 3 3 block of K, so 3.732 -1 0 K = -1 4 -1 . (11.51) 0 -1 3.732
rows and columns of K corresponding to the interior nodes are seen to all have the same form. Namely, if i labels an interior node, then the corresponding diagonal entry is kii = 4, while the off-diagonal entries kij = kji , i = j, are equal to either -1 when node i is adjacent to node j on the grid, and is equal to 0 in all other cases. Node j is allowed to be a boundary node. (Interestingly, the result does not depend on how one orients the pair of triangles making up each square of the grid, which only plays a role in the computation of the right hand side of the finite element equation.) Observe that the same computation applies even to our coarse triangulation. The interior node 2 belongs to all right isoceles triangles, and the corresponding nonzero entries in (11.50) are k22 = 4 and k21 = k23 = k24 = k29 = -1, indicating the four adjacent nodes. Remark : Interestingly, the coefficient matrix arising from the finite element method on a square (or even rectangular) grid is the same as the coefficient matrix arising from a finite difference solution to the Laplace or Poisson equation, as described in Example 5.7. The finite element approach has the advantage of applying to much more general triangulations. The Coefficient Vector and the Boundary Conditions So far, we have been concentrating on assembling the finite element coefficient matrix T K. We also need to compute the forcing vector b = ( b1 , b2 , . . . , bn ) appearing on the right hand side of the fundamental linear equation (11.41). According to (11.39), the entries bi are found by integrating the product of the forcing function and the finite element basis function. As before, we will approximate the integral over the domain by an integral over the triangles, and so bi =
f (x, y) i(x, y) dx dy
T
f (x, y) i (x, y) dx dy
b . i
(11.52)
Typically, the exact computation of the various triangular double integrals is not so convenient, and we resort to a numerical approximation. Since we are assuming that the individual triangles are small, we can get away with a very crude numerical integration scheme. If the function f (x, y) does not vary much over the triangle T -- which will certainly be the case if T is sufficiently small -- we may approximate f (x, y) c for i (x, y) T by a constant. The integral (11.52) is then approximated by b = i
f (x, y) i (x, y) dx dy c i i (x, y) dx dy = 1 3
T
T
c area T = i
1 6
c | |. i
(11.53) The formula for the integral of the affine element i (x, y) follows from solid geometry: it equals the volume under its graph, a tetrahedron of height 1 and base T , as illustrated in Figure 11.13. How to choose the constant c ? In practice, the simplest choice is to let c = f (xi , yi ) i i be the value of the function at the ith vertex. With this choice, the sum in (11.52) becomes bi
1 3
f (xi , yi ) area T = 414
1 3
f (xi , yi ) area Pi ,
c 2012
(11.54)
1/19/12
Peter J. Olver
Figure 11.13.
Finite Element Tetrahedron.
where Pi is the vertex polygon (11.38) corresponding to the node xi . In particular, for the square mesh with the uniform choice of triangles, as in Example 11.5, area Pi = 3 h2 for all i, and so bi f (xi , yi ) h2 (11.55)
is well approximated by just h2 times the value of the forcing function at the node. This is the underlying reason to choose the uniform triangulation for the square mesh; the alternating version would give unequal values for the bi over adjacent nodes, and this would introduce unnecessary errors into the final approximation. Example 11.8. For the coarsely triangulated oval plate, the reduced stiffness matrix is (11.51). The Poisson equation - u = 4 models a constant external heat source of magnitude 4 over the entire plate. If we keep the edges of the plate fixed at 0 , then we need to solve the finite element equation K c = b, where K is the coefficient matrix (11.51). The entries of b are, by (11.54), equal to 4 = f (xi , yi ) times one third the area of the corresponding vertex polygon, which for node 1 2 is the square consisting of 4 right triangles, each of area 2 , whereas for nodes 1 and 3 it 1 consists of 4 right triangles of area 2 plus three equilateral triangles, each of area 43 ; see Figure 11.11. Thus, b=
4 3
2+
3 3 4 ,
2, 2 +
3 3 4
T
= ( 4.3987, 2.6667, 4.3987 ) .
T
The solution to the final linear system is easily found: c = ( 1.5672, 1.4503, 1.5672 ) . Its entries are the values of the finite element approximation at the three interior nodes. The finite element solution is plotted in the first illustration in Figure 11.14. A more accurate solution, based on a square grid triangulation of size h = .1 is plotted in the second figure. 1/19/12 415
c 2012 Peter J. Olver T
Figure 11.14.
Finite Element Solutions to Poisson's Equation for an Oval Plate.
Inhomogeneous Boundary Conditions So far, we have restricted our attention to problems with homogeneous Dirichlet boundary conditions. According to Theorem 10.31, the solution to the inhomogeneous Dirichlet problem - u = f in , u = h on , is also obtained by minimizing the Dirichlet functional (10.79). However, now the minimization takes place over the set of functions that satisfy the inhomogeneous boundary conditions. It is not difficult to fit this problem into the finite element scheme. The elements corresponding to the interior nodes of our triangulation remain as before, but now we need to include additional elements to ensure that our approximation satisfies the boundary conditions. Note that if xk is a boundary node, then the corresponding boundary element k (x, y) satisfies (11.29), and so has the same piecewise affine form (11.37). The corresponding finite element approximation
m
w(x, y) =
i=1
ci i (x, y),
(11.56)
has the same form as before, (11.30), but now the sum is over all nodes, both interior and boundary. As before, the coefficients ci = w(xi , yi ) u(xi , yi ) are the values of the finite element approximation at the nodes. Therefore, in order to satisfy the boundary conditions, we require cj = h(xj , yj ) whenever xj = (xj , yj ) is a boundary node. (11.57)
Remark : If the boundary node xj does not lie precisely on the boundary , we need to approximate the value h(xj , yj ) appropriately, e.g., by using the value of h(x, y) at the nearest boundary point (x, y) . The derivation of the finite element equations proceeds as before, but now there are additional terms arising from the nonzero boundary values. Leaving the intervening details to the reader, the final outcome can be written as follows. Let K denote the full m m finite element matrix constructed as above. The reduced coefficient matrix K is obtained 1/19/12 416
c 2012 Peter J. Olver
Figure 11.15.
Solution to the Dirichlet Problem for the Oval Plate.
by retaining the rows and columns corresponding to only interior nodes, and so will have size n n , where n is the number of interior nodes. The boundary coefficient matrix K is the n (m - n) matrix consisting of the entries of the interior rows that do not appear in K, i.e., those lying in the columns indexed by the boundary nodes. For instance, in the the coarse triangulation of the oval plate, the full finite element matrix is given in (11.50), and the upper 3 3 subblock is the reduced matrix (11.51). The remaining entries of the first three rows form the boundary coefficient matrix 0 -.7887 -.5774 -.5774 -.7887 0 0 0 0 0 K = -1 0 0 0 0 -1 0 0 0 0 . 0 0 0 0 0 0 -.7887 -.5774 -.5774 -.7887 (11.58) We similarly split the coefficients ci of the finite element function (11.56) into two groups. We let c R n denote the as yet unknown coefficients ci corresponding to the values of the approximation at the interior nodes xi , while h R m-n will be the vector of boundary values (11.57). The solution to the finite element approximation (11.56) is obtained by solving the associated linear system K c + K h = b, or K c = f = b - K h. (11.59)
Example 11.9. For the oval plate discussed in Example 11.7, suppose the right hand semicircular edge is held at 10 , the left hand semicircular edge at -10 , while the two straight edges have a linearly varying temperature distribution ranging from -10 at the left to 10 at the right, as illustrated in Figure 11.15. Our task is to compute its equilibrium temperature, assuming no internal heat source. Thus, for the coarse triangulation we have the boundary nodes values h = ( h4 , . . . , h13 ) = ( 0, -10, -10, -10, -10, 0, 10, 10, 10, 10 ) . Using the previously computed formulae (11.51, 58) for the interior coefficient matrix K and boundary coefficient matrix K, we approximate the solution to the Laplace equation 1/19/12 417
c 2012 Peter J. Olver T T
by solving (11.59). We are assuming that there is no external forcing function, f (x, y) 0, T and hence b = 0, and so we must solve K c = f = - K h = ( 2.1856, 3.6, 7.6497 ) . The T finite element function corresponding to the solution c = ( 1.0679, 1.8, 2.5320 ) is plotted in the first illustration in Figure 11.15. Even on such a coarse mesh, the approximation is not too bad, as evidenced by the second illustration, which plots the finite element solution for a square mesh with spacing h = .2 between nodes.
1/19/12
418
c 2012
Peter J. Olver