Lecture AIMS Notes 2006
Peter J. Olver
14. Finite Elements
In this part, we introduce the powerful finite element method for finding numerical approximations to the solutions to boundary value problems involving both ordinary and partial differential equations can be solved by direct integration. The method relies on the characterization of the solution as the minimizer of a suitable quadratic functional. The innovative idea is to restrict the infinite-dimensional minimization principle characterizing the exact solution to a suitably chosen finite-dimensional subspace of the function space. When properly formulated, the solution to the resulting finite-dimensional minimization problem approximates the true minimizer. The finite-dimensional minimizer is found by solving the induced linear algebraic system, using either direct or iterative methods. We begin with one-dimensional boundary value problems involving ordinary differential equations, and, in the final section, show how to adapt the finite element analysis to partial differential equations, specifically the two-dimensional Laplace and Poisson equations.
14.1. Finite Elements for Ordinary Differential Equations.
The characterization of the solution to a linear boundary value problem via a quadratic minimization principle inspires a very powerful and widely used numerical solution scheme, known as the finite element method . In this final section, we give a brief introduction to the finite element method in the context of one-dimensional boundary value problems involving ordinary differential equations. The underlying idea is strikingly simple. We are trying to find the solution to a boundary value problem by minimizing a quadratic functional P[ u ] on an infinite-dimensional vector space U . The solution u U to this minimization problem is found by solving a differential equation subject to specified boundary conditions. However, minimizing the functional on a finite-dimensional subspace W U is a problem in linear algebra, and, moreover, one that we already know how to solve! Of course, restricting the functional P[ u ] to the subspace W will not, barring luck, lead to the exact minimizer. Nevertheless, if we choose W to be a sufficiently "large" subspace, the resulting minimizer w W may very well provide a reasonable approximation to the actual solution u U . A rigorous justification of this process, under appropriate hypotheses, requires a full analysis of the finite element method, and we refer the interested reader to [45, 49]. Here we shall concentrate on trying to understand how to apply the method in practice. 3/15/06 233
c 2006 Peter J. Olver
To be a bit more explicit, consider the minimization principle P[ u ] = for the linear system K[ u ] = f, where K = L L,
1 2
L[ u ]
2
- f ,u
(14.1)
representing our boundary value problem. The norm in (14.1) is typically based on some form of weighted inner product v , v on the space of strains v = L[ u ] V , while the inner product term f , u is typically (although not necessarily) unweighted on the space of displacements u U . The linear operator takes the self-adjoint form K = L L, and must be positive definite -- which requires ker L = {0}. Without the positivity assumption, the boundary value problem has either no solutions, or infinitely many; in either event, the basic finite element method will not apply. Rather than try to minimize P[ u ] on the entire function space U , we now seek to minimize it on a suitably chosen finite-dimensional subspace W U . We begin by selecting a basis 1 , . . . , n of the subspace W . The general element of W is a (uniquely determined) linear combination (x) = c1 1 (x) + + cn n (x) (14.2) of the basis functions. Our goal, then, is to determine the coefficients c 1 , . . . , cn such that (x) minimizes P[ ] among all such functions. Substituting (14.2) into (14.1) and expanding we find 1 m c c - P[ ] = b c = 2 i,j = 1 ij i j i = 1 i i where (a) c = ( c1 , c2 , . . . , cn ) is the vector of unknown coefficients in (14.2), (b) M = (mij ) is the symmetric n n matrix with entries mij =
T T n n 1 2
cT M c - cT b,
(14.3)
L[ i ] , L[ j ] ,
i, j = 1, . . . , n,
(14.4)
(c) b = ( b1 , b2 , . . . , bn ) is the vector with entries bi = f , i , i = 1, . . . , n. (14.5)
Observe that, once we specify the basis functions i , the coefficients mij and bi are all known quantities. Therefore, we have reduced our original problem to a finite-dimensional problem of minimizing the quadratic function (14.3) over all possible vectors c R n . The coefficient matrix M is, in fact, positive definite, since, by the preceding computation,
n
c Mc =
i,j = 1
T
mij ci cj = L[ c1 1 (x) + + cn n ]
2
= L[ ]
2
> 0
(14.6)
In this case, an orthonormal basis is not of any particular help.
3/15/06
234
c 2006
Peter J. Olver
as long as L[ ] = 0. Moreover, our positivity assumption implies that L[ ] = 0 if and only if 0, and hence (14.6) is indeed positive for all c = 0. We can now invoke the original finite-dimensional minimization Theorem 12.12 to conclude that the unique minimizer to (14.3) is obtained by solving the associated linear system M c = b. (14.7)
Solving (14.7) relies on some form of Gaussian Elimination, or, alternatively, an iterative linear system solver, e.g., GaussSeidel or SOR. This constitutes the basic abstract setting for the finite element method. The main issue, then, is how to effectively choose the finite-dimensional subspace W . Two candidates that might spring to mind are the space P (n) of polynomials of degree n, or the space T (n) of trigonometric polynomials of degree n. However, for a variety of reasons, neither is well suited to the finite element method. One criterion 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, the linear algebraic system (14.7) will typically be rather large, and so the coefficient matrix M should be as sparse as possible, i.e., have lots of zero entries. Otherwise, computing the solution will be too timeconsuming to be of much practical value. Such considerations prove to be of absolutely crucial importance when applying the method to solve boundary value problems for partial differential equations in higher dimensions. The really 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 P[ u ]. The governing differential equation requires its solutions to have a certain degree of smoothness, whereas the associated minimization principle typically requires only half as many derivatives. Thus, for second order boundary value problems, P[ u ] only involves first order derivatives. It can be rigorously shown that the functional has the same minimizing solution, even if one allows (reasonable) functions that fail to have enough derivatives to satisfy the differential equation. Thus, one can try minimizing over subspaces containing fairly "rough" functions. Again, the justification of this method requires some deeper analysis, which lies beyond the scope of this introductory treatment. For second order boundary value problems, a popular and effective choice of the finitedimensional subspace is to use continuous, piecewise affine functions. Recall that a function is affine, f (x) = a x + b, if and only if its graph is a straight line. The function is piecewise affine if its graph consists of a finite number of straight line segments; a typical example is plotted in Figure 14.1. Continuity requires that the individual line 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 mesh points a = x0 < x1 < x2 < < xn-1 < xn = b. The formulas simplify if one uses equally spaced mesh points, but this is not necessary for the method to apply. Let W denote the vector space consisting of all continuous, piecewise affine functions, with corners at the nodes, that satisfy the homogeneous boundary 3/15/06 235
c 2006 Peter J. Olver
0.8 0.6 0.4 0.2 0.2 0.4 0.6 0.8 1
Figure 14.1.
A Continuous Piecewise Affine Function.
conditions. To be specific, let us treat the case of Dirichlet (fixed) boundary conditions (a) = (b) = 0. Thus, on each subinterval (x) = cj + bj (x - xj ), Continuity of (x) requires cj = (x+ ) = (x- ) = cj-1 + bj-1 hj-1 , j j j = 1, . . . , n - 1, (14.9) for xj x xj+1 , j = 0, . . . , n - 1. (14.8)
where hj-1 = xj -xj-1 denotes the length of the j th subinterval. The boundary conditions (14.8) require (a) = c0 = 0, (b) = cn-1 + hn-1 bn-1 = 0. (14.10) The function (x) involves a total of 2 n unspecified coefficients c0 , . . . , cn-1 , b0 , . . . , bn-1 . The continuity conditions (14.9) and the second boundary condition (14.10) uniquely determine the bj . The first boundary condition specifies c0 , while the remaining n - 1 coefficients c1 = (x1 ), . . . , cn-1 = (xn-1 ) are arbitrary. We conclude that the finite element subspace W has dimension n - 1, which is the number of interior mesh points. Remark : Every function (x) in our subspace has piecewise constant first derivative w (x). However, the jump discontinuities in (x) imply that its second derivative (x) has a delta function impulse at each mesh point, and is therefore far from being a solution to the differential equation. Nevertheless, the finite element minimizer (x) will, in practice, provide a reasonable approximation to the actual solution u (x). The most convenient basis for W consists of the hat functions, which are continuous, piecewise affine functions that interpolate the same basis data as the Lagrange polynomials (13.22), namely j (xk ) = 3/15/06 1, 0, j = k, j = k, for 236 j = 1, . . . , n - 1, k = 0, . . . , n.
c 2006 Peter J. Olver
1.2 1 0.8 0.6 0.4 0.2 1 -0.2 2 3 4 5 6 7
Figure 14.2.
A Hat Function.
The graph of a typical hat function appears in Figure 14.2. The explicit formula is easily established: x-x j-1 , xj-1 x xj , x -x j j-1 xj+1 - x j (x) = j = 1, . . . , n - 1. (14.11) , xj x xj+1 , x j+1 - xj 0, x xj-1 or x xj+1 , An advantage of using these basis elements is that the resulting coefficient matrix (14.4) turns out to be tridiagonal. Therefore, the tridiagonal Gaussian Elimination algorithm in (4.47) will rapidly produce the solution to the linear system (14.7). Since the accuracy of the finite element solution increases with the number of mesh points, this solution scheme allows us to easily compute very accurate numerical approximations. Example 14.1. Consider the equilibrium equations K[ u ] = - d dx c(x) du dx = f (x), 0<x< ,
for a non-uniform bar subject to homogeneous Dirichlet boundary conditions. In order to formulate a finite element approximation scheme, we begin with the minimization principle based on the quadratic functional P[ u ] =
1 2
u
2
- f ,u =
0
1 2
c(x) u (x)2 - f (x) u(x) dx.
We divide the interval [ 0, ] into n equal subintervals, each of length h = /n. The resulting uniform mesh has j xj = j h = , j = 0, . . . , n. n 3/15/06 237
c 2006 Peter J. Olver
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, The associated linear system (14.7) has coefficient matrix entries mij = i , j =
0
(14.12)
i (x) j (x) c(x) dx,
i, j = 1, . . . , n - 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, j = 1, . . . , n - 1. j (x) = - 1/h, xj x xj+1 , 0, otherwise, Therefore, the coefficient matrix has the tridiagonal form s + s -s
0 1 1
where
- s1 1 M= 2 h
s1 + s 2 - s2
..
- s2 s2 + s 3 - sn-3
.
..
.. . . sn-3 + sn-2 - sn-2 - sn-2 sn-2 + sn-1
- s3
,
(14.13)
xj+1
sj =
c(x) dx
xj
(14.14)
is the total stiffness of the j th subinterval. For example, in the homogeneous case c(x) 1, the coefficient matrix (14.13) reduces to the very special form 2 -1 -1 2 -1 -1 2 -1 1 .. .. .. M= (14.15) . . . . h -1 2 -1 -1 2 The corresponding right hand side has entries bj = f , j = f (x) j (x) dx
xj xj-1 xj+1
0
1 = h 3/15/06
(14.16) (xj+1 - x)f (x) dx .
c 2006 Peter J. Olver
(x - xj-1 )f (x) dx + 238
xj
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 14.3.
Finite Element Solution to (14.19).
In this manner, we have assembled the basic ingredients for determining the finite element approximation to the solution. In practice, we do not have to explicitly evaluate the integrals (14.14, 16), but may replace them by a suitably close numerical approximation. When h 1 is small, then the integrals are taken over small intervals, and we can use the trapezoid rule , [5, 7], to approximate them: sj h c(xj ) + c(xj+1 ) , 2 bj h f (xj ). (14.17)
Remark : The j th entry of the resulting finite element system M c = b is, upon dividing by h, given by - u(xj+1 ) - 2 u(xj ) + u(xj-1 ) cj+1 - 2 cj + cj-1 = - = -f (xj ). 2 h h2 (14.18)
The left hand side coincides with the standard finite difference approximation to minus the second derivative - u (xj ) at the mesh point xj . As a result, for this particular differential equation, the finite element and finite difference numerical solution methods happen to coincide. Example 14.2. Consider the boundary value problem - d du (x + 1) = 1, dx dx u(0) = 0, u(1) = 0. (14.19)
One might be tempted use more accurate numerical integration procedures, but the improvement in accuracy of the final answer is not very significant, particularly if the step size h is small.
3/15/06
239
c 2006
Peter J. Olver
The explicit solution is easily found by direct integration: u(x) = - x + log(x + 1) . log 2 (14.20)
It minimizes the associated quadratic functional P[ u ] =
0 1 2
(x + 1) u (x)2 - u(x) dx
(14.21)
over all possible functions u C1 that satisfy the given boundary conditions. The finite element system (14.7) has coefficient matrix given by (14.13) and right hand side (14.16), where
xj+1
sj =
xj
(1 + x) dx = h (1 + xj ) +
1 2
h2 = h + h 2
j+
1 2
xj+1
,
bj =
1 dx = h.
xj
The resulting solution is plotted in Figure 14.3. The first three graphs contain, respectively, 5, 10, 20 points in the mesh, so that h = .2, .1, .05, while the last plots the exact solution (14.20). Even when computed on rather coarse meshes, the finite element approximation is quite respectable. Example 14.3. Consider the SturmLiouville boundary value problem - u + (x + 1)u = x ex , u(0) = 0, u(1) = 0. (14.22)
The solution minimizes the quadratic functional
1
P[ u ] =
0
1 2
1 u (x)2 + 2 (x + 1) u(x)2 - ex u(x) dx,
(14.23)
while
over all functions u(x) that satisfy the boundary conditions. We lay out a uniform mesh of step size h = 1/n. The corresponding finite element basis hat functions as in (14.12). The matrix entries are given by 2 2h + (xi + 1), i = j, h 3 1 h mij = i (x) j (x) + (x + 1) i (x) j (x) dx 1 - + (xi + 1), | i - j | = 1, 0 h 6 0, otherwise,
1
bi = x e x , i =
0
x ex i (x) dx xi exi h.
The integration is made easier by noting that the integrand is zero except on a small subinterval. Since the function x + 1 (but not i or j ) does not vary significantly on this subinterval, it can be approximated by its value 1 + xi at a mesh point. A similar simplification is used in the ensuing integral for bi .
3/15/06
240
c 2006
Peter J. Olver
0.1 0.08 0.06 0.04 0.02
0.1 0.08 0.06 0.04 0.02
0.2
0.4
0.6
0.8
1
0.2
0.4
0.6
0.8
1
0.1 0.08 0.06 0.04 0.02
0.1 0.08 0.06 0.04 0.02
0.2
0.4
0.6
0.8
1
0.2
0.4
0.6
0.8
1
Figure 14.4.
Finite Element Solution to (14.22).
The resulting solution is plotted in Figure 14.4. As in the previous figure, the first three graphs contain, respectively, 5, 10, 20 points in the mesh, while the last plots the exact solution, which can be expressed in terms of Airy functions, cf. [36]. So far, we have only treated homogeneous boundary conditions. An inhomogeneous boundary value problem does not immediately fit into our framework since the set of functions satisfying the boundary conditions does not form a vector space. One way to get around this problem is to replace u(x) by u(x) = u(x) - h(x), where h(x) is any convenient function that satisfies the boundary conditions. For example, for the inhomogeneous Dirichlet conditions u(a) = , u(b) = , we can subtract off the affine function h(x) = ( - )x + b - a . b-a
Another option is to choose an appropriate combination of elements at the endpoints: h(x) = 0 (x) + n (x). Linearity implies that the difference u(x) = u(x) - h(x) satisfies the amended differential equation K[ u ] = f , where f = f - K[ h ], now supplemented by homogeneous boundary conditions. The modified boundary value problem can then be solved by the standard finite element method. Further details are left as a project for the motivated student. Finally, one can employ other functions beyond the piecewise affine hat functions (14.11) to span finite element subspace. Another popular choice, which is essential for higher order boundary value problems such as beams, is to use splines. Thus, once we have chosen our mesh points, we can let j (x) be the basis Bsplines. Since j (x) = 0 3/15/06 241
c 2006 Peter J. Olver
for x xj-2 or x xj+2 , the resulting coefficient matrix (14.4) is pentadiagonal , which means mij = 0 whenever | i - j | > 2. Pentadiagonal matrices are not quite as pleasant as their tridiagonal cousins, but are still rather sparse. Positive definiteness of M implies that an iterative solution technique, e.g., SOR, can effectively and rapidly solve the linear system, and thereby produce the finite element spline approximation to the boundary value problem.
14.2. Finite Elements for the Laplace and Poisson Equations.
Finite element methods are also effectively employed to solving boundary value problems for elliptic partial differential equations. In this section, we concentrate on applying these ideas to the two-dimensional Poisson equation. For specificity, we concentrate on the homogeneous Dirichlet boundary value problem. Theorem 14.4. The function u(x, y) that minimizes the Dirichlet integral
1 2
u
2
- u,f =
1 2
u2 + x
1 2
u2 - f u dx dy y
(14.24)
among all C1 functions that satisfy the prescribed homogeneous Dirichlet boundary conditions is the solution to the boundary value problem - u = f in u=0 on . (14.25)
In the 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 possible solutions to the partial differential equation. Nevertheless, they do provide good approximations to the actual solution. An important practical consideration, impacting the speed of the calculation, is to employ functions with small support. The resulting finite element matrix will then be sparse and the solution to the linear system can be relatively rapidly calculate, usually by application of an iterative numerical scheme such as the GaussSeidel or SOR methods discussed in Section 7.4. Finite Elements and Triangulation For one-dimensional boundary value problems, the finite element construction rests on the introduction of a mesh a = x0 < x1 < < xn = b on the interval of definition. The mesh nodes xk break the interval into a collection of small subintervals. In two-dimensional problems, a mesh consists of a finite number of points xk = (xk , yk ), k = 1, . . . , m, known as nodes, usually lying inside the domain R 2 . As such, there is considerable freedom in the choice of mesh nodes, and completely uniform spacing is often not possible. We regard the nodes as forming the vertices of a triangulation of the domain , consisting of a finite number of small triangles, which we denote by T1 , . . . , TN . The nodes are split into two categories -- interior nodes and boundary nodes, the latter lying on or close to the boundary of the domain. A curved boundary is approximated by the polygon through the boundary nodes formed by the sides of the triangles lying on the edge of the domain; see Figure 14.5 for a typical example. Thus, in computer implementations of the finite element 3/15/06 242
c 2006 Peter J. Olver
Figure 14.5.
Triangulation of a Planar Domain.
method, the first module is a routine that will automatically triangulate a specified domain in some reasonable manner; see below for details on what "reasonable" entails. 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. "Piecewise affine" means that, on each triangle, the graph of w is flat, and so has the formula w(x, y) = + x + y, for (x, y) T . (14.26)
Continuity of w requires that its values on a common edge between two triangles must agree, and this will impose certain compatibility conditions on the coefficients , , and , , associated with adjacent pairs of triangles T , T . The graph of z = w(x, y) forms a connected polyhedral surface whose triangular faces lie above the triangles in the domain; see Figure 14.5 for an illustration. The next step is to choose a basis of the subspace of piecewise affine functions for the given triangulation. As in the one-dimensional version, the most convenient basis consists of pyramid functions k (x, y) which assume the value 1 at a single node xk , and are zero at all the other nodes; thus 1, i = k, (14.27) 0, i = k. Note that k will be nonzero only on those triangles which have the node xk as one of their vertices, and hence the graph of k looks like a pyramid of unit height sitting on a flat plane, as illustrated in Figure 14.7. k (xi , yi ) =
Here and subsequently, the index is a superscript, not a power!
3/15/06
243
c 2006
Peter J. Olver
Figure 14.6.
Piecewise Affine Function.
Figure 14.7.
Finite Element Pyramid Function.
The pyramid functions k (x, y) corresponding to 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, which is supposed to be a good approximation to the curved boundary of the original domain . Thus, the finite-dimensional finite element subspace W is the span of the interior node pyramid functions, and so general element w W is a linear combination thereof:
n
w(x, y) =
k=1
ck k (x, y),
(14.28)
where the sum ranges over the n interior nodes of the triangulation. Owing to the original specification (14.27) of the pyramid functions, the coefficients ck = w(xk , yk ) u(xk , yk ), 3/15/06 244 k = 1, . . . , n,
c 2006
(14.29)
Peter J. Olver
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 c 1 = = cn = 0. Thus, the interior node pyramid functions 1 , . . . n form a basis for finite element subspace W , which therefore has dimension equal to n, the number of interior nodes. Determining the explicit formulae for the finite element basis 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 (14.26) that takes the value 1 at the vertex xk and 0 at its other two vertices xl and xm . Thus, we are in need of a formula for an affine function or element
k (x, y) = k + k x + k y,
(x, y) T ,
(14.30)
that takes the prescribed values
k (xi , yi ) = k (xj , yj ) = 0, k (xk , yk ) = 1,
at three distinct points. These three conditions lead to the linear system
k (xi , yi ) = k + k xi + k yi = 0, k (xj , yj ) = k + k xj + k yj = 0, k (xk , yk ) = k + k xk + k yk = 1.
(14.31)
The solution produces the explicit formulae
k =
x i yj - x j yi ,
k =
yi - y j ,
k =
xj x - i ,
(14.32)
for the coefficients; the denominator
is, up to sign, twice the area of the triangle T .
1 = det 1 1
xi xj xk
yi yj = 2 area T yk
(14.33)
Example 14.5. Consider an isoceles right triangle T with vertices x1 = (0, 0), x2 = (1, 0), x3 = (0, 1).
Using (14.3233) (or solving the linear systems (14.31) directly), we immediately produce the three affine elements 1 (x, y) = 1 - x - y, 2 (x, y) = x, 3 (x, y) = y. (14.34)
As required, each k equals 1 at the vertex xk and is zero at the other two vertices. The finite element pyramid function is then obtained by piecing together the individual affine elements, whence k (x, y) = 3/15/06
k (x, y),
if (x, y) T which has xk as a vertex, otherwise. 245
c 2006
0,
(14.35)
Peter J. Olver
Figure 14.8.
Vertex Polygons.
Figure 14.9.
Square Mesh Triangulations.
Continuity of k (x, y) is assured since the constituent affine elements have the same values at common vertices. The support of the pyramid function (14.35) is the polygon supp k = Pk =
[
T
(14.36)
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 those that are connected to xk by a single edge of the triangulation. In Figure 14.8 the shaded regions indicate two of the vertex polygons for the triangulation in Figure 14.5. Example 14.6. The simplest, and most common triangulations are based on regular meshes. Suppose that the nodes lie on a square grid, and so are of the form x i,j = (i h + a, j h + b) where h > 0 is the inter-node spacing, and (a, b) represents an overall offset. If we choose the triangles to all have the same orientation, as in the first picture in Figure 14.9, then the vertex polygons all have the same shape, consisting of 6 triangles 3/15/06 246
c 2006 Peter J. Olver
of total area 3 h2 -- the shaded region. On the other hand, if we choose an alternating, perhaps more sthetically pleasing triangulation as in the second picture, then there are two types of vertex polygons. The first, consisting of four triangles, has area 2 h 2 , 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, in order to ensure convergence of the finite element solution to the true minimizer, one should choose a triangulation with the following properties: (a) The triangles are not too long and skinny. In other words, their sides should have comparable lengths. In particular, obtuse triangles should be avoided. (b) The areas of nearby triangles T should not vary too much. (c) The areas of nearby vertex polygons Pk should also not vary too much. For adaptive or variable meshes, one might very well have wide variations in area over the entire grid, with small triangles in regions of rapid change in the solution, and large ones in less interesting regions. But, overall, the sizes of the triangles and vertex polygons should not dramatically vary as one moves across the domain. The Finite Element Equations We now seek to approximate the solution to the homogeneous Dirichlet boundary value problem by restricting the Dirichlet functional to the selected finite element subspace W . Substituting the formula (14.28) for a general element of W into the quadratic Dirichlet functional (14.24) and expanding, we find
n n 2 n
P[ w ] = P
c i i
=
i=1
=
1 k c c - b c = 2 i,j = 1 ij i j i = 1 i i
n
ci
i
- f (x, y)
1 2
c i i
i=1
i=1
n
cT K c - bT c.
T
dx dy
Here, K = (kij ) is the symmetric n n matrix, while b = ( b1 , b2 , . . . , bn ) is the vector that have the respective entries kij = i , j = i j dx dy, (14.37)
bi = f , i =
f i dx dy.
Thus, to determine the finite element approximation, we need to minimize the quadratic function P (c) = 1 cT K c - bT c (14.38) 2 over all possible choices of coefficients c = ( c1 , c2 , . . . , cn ) R n , i.e., over all possible function values at the interior nodes. Restricting to the finite element subspace has reduced us to a standard finite-dimensional quadratic minimization problem. First, the coefficient matrix K > 0 is positive definite due to the positive definiteness of the original functional; 3/15/06 247
c 2006 Peter J. Olver T
the proof in Section 14.1 is easily adapted to the present situation. The minimizer is obtained by solving the associated linear system K c = b. (14.39)
The solution to (14.39) can be effected by either Gaussian elimination or an iterative technique. To find explicit formulae for the matrix coefficients kij in (14.37), we begin by noting that the gradient of the affine element (14.30) is equal to
k (x, y) = a = k k k
=
1
yi - y j , xj - x i
(x, y) T ,
(14.40)
which is a constant vector inside the triangle T , while outside k (x, y) =
k = a , k
k = 0. Therefore,
if (x, y) T which has xk as a vertex, otherwise,
0,
(14.41)
reduces to a piecewise constant function on the triangulation. Actually, (14.41) is not quite correct since if (x, y) lies on the boundary of a triangle T , then the gradient does not exist. However, this technicality will not cause any difficulty in evaluating the ensuing integral. We will approximate integrals over the domain by integrals over the 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 .
(14.42)
Now, according to (14.41), 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, 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 (14.41),
kij =
T
a a dx dy = a a area T = i j i j
1 2
a a | | . i j
Let T have vertices xi , xj , xk . Then, by (14.40, 41, 33),
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 | | = . (14.43) kii = 2 ( )2 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
3/15/06
248
c 2006
Peter J. Olver
Figure 14.10.
Right and Equilateral Triangles.
elemental stiffnesses depend only on the angles of the triangle and not on its size. Thus, similar triangles have the same elemental stiffnesses. Indeed, if we denote the angle in T at the vertex xk by k , then
kii = 1 2 cot k + cot j ,
while
kij = kji = - 1 cot k , 2
i = j,
(14.44)
depend only upon the cotangents of the angles. Example 14.7. 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. (14.45)
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 elemental stiffnesses are k11 = k22 = k33 =
1 3
.577350,
1 k12 = k21 = k13 = k31 = k23 = k32 = - 23 -.288675.
(14.46)
Assembling the Elements The elemental stiffnesses of each triangle will contribute, through the summation (14.42), 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, K =
N
K ,
=1
(14.47)
to produce the full finite element matrix, in accordance with (14.42). 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 (14.39) only refers to the n 3/15/06 249
c 2006 Peter J. Olver
Figure 14.11.
The Oval Plate.
5 4 13 12 1 7 8 9 10 2 3 11
1 2 3
4 5 9 8 6 10
7
12 13
6
11
14
Triangles Figure 14.12.
Nodes A Coarse Triangulation of the Oval Plate.
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 k ij when both xi and xj are interior nodes. For the homogeneous boundary value problem, this is all we require. As we shall 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 14.8. 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 14.11. 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 Poisson equation with Dirichlet boundary conditions, for the equilibrium temperature u(x, y). Let us describe how to set up the finite element approximation to such a boundary value problem. 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 14.12. There are 13 nodes in all, 3/15/06 250
c 2006 Peter J. Olver
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 1313, its rows and columns labeled by all the nodes, while the reduced matrix K appearing in the finite element equations (14.39) 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 T 1 is equilateral, and so has elemental stiffnesses (14.46). Its vertices are labeled 1, 5, and 6, and therefore we place the stiffnesses (14.46) in the rows and columns numbered 1, 5, 6 to form the summand
.577350 0 0 0 -.288675 K1 = -.288675 0 0 . . .
0 0 0 0 0 0 0 0 . . .
0 0 -.288675 0 0 0 0 0 0 0 0 0 0 0 .577350 0 0 -.288675 0 0 0 0 0 0 . . . . . . . . .
-.288675 0 0 0 -.288675 .577350 0 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 T 2 has the same equilateral elemental stiffness matrix (14.46), but now its vertices are 1, 6, 7, and so it will contribute
.577350 0 0 0 0 K2 = -.288675 -.288675 0 . . .
0 0 0 0 0 0 0 0 . . .
0 0 0 -.288675 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 .577350 0 0 0 -.288675 0 0 0 0 . . . . . . . . . . . .
-.288675 0 0 0 0 -.2886750 .5773500 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 (14.45). Its vertices are labeled 1, 4, and 5, with 3/15/06 251
c 2006 Peter J. Olver
Figure 14.13.
A Square Mesh for the Oval Plate.
vertex 5 at the right angle. Therefore, its contribution is .5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 .5 -.5 0 0 -.5 K4 = 0 0 0 0 0 0 0 0 0 0 0 0 . . . . . . . . . . . . -.5 0 0 0 0 0 -.5 0 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 ... . . .. . . . . .
.
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
K = K1 + K2 + 0 3.732 -1 B -1 4 B B B 0 -1 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 2006
(14.48)
3/15/06
252
Peter J. Olver
-.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 C -.7887 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
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 . (14.49) 0 -1 3.732 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 14.13 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 (14.45). Summing the corresponding matrices K over all the triangles, as in (14.47), the 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 entries in (14.48) are k22 = 4, and k2j = -1 for the four adjacent nodes j = 1, 3, 4, 9. 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 3/15/06 253
c 2006 Peter J. Olver
finite difference solution to the Laplace or Poisson equation. The finite element approach has the advantage of applying to much more general triangulations. In general, while the finite element matrix K for a two-dimensional boundary value problem is not as nice as the tridiagonal matrices we obtained in our one-dimensional problems, it is still very sparse and, on regular grids, highly structured. This makes solution of the resulting linear system particularly amenable to an iterative matrix solver such as GaussSeidel, Jacobi, or, for even faster convergence, successive over-relaxation (SOR). 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 (14.39). According to (14.37), the entries b i 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 i dx dy
f i dx dy
T
b . i
(14.50)
Typically, the exact computation of the various triangular integrals is not convenient, and so we resort to a numerical approximation. Since we are assuming that the individual triangles are small, we can adopt 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 (x, y) T by a constant. The i integral (14.50) is then approximated by b = i
f i dx dy c i i (x, y) dx dy = 1 3
T
T
c area T = i
1 6
c | |. i
(14.51)
The formula for the integral of the affine element i (x, y) follows from solid geometry. Indeed, it equals the volume under its graph, a tetrahedron of height 1 and base T , as illustrated in Figure 14.14. 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 (14.50) becomes
bi
1 3
f (xi , yi ) area T =
1 3
f (xi , yi ) area Pi ,
(14.52)
where Pi is the vertex polygon (14.36) corresponding to the node xi . In particular, for the square mesh with the uniform choice of triangles, as in Example 14.6, area Pi = 3 h2 for all i, and so bi f (xi , yi ) h2 (14.53)
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. 3/15/06 254
c 2006 Peter J. Olver
Figure 14.14.
Finite Element Tetrahedron.
Example 14.9. For the coarsely triangulated oval plate, the reduced stiffness matrix is (14.49). 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 (14.49), while b=
4 3
2+
3 3 4 ,
2, 2 +
3 3 4
T
= ( 4.39872, 2.66667, 4.39872 ) .
T
The entries of b are, by (14.52), equal to 4 = f (xi , yi ) times one third the area of the corresponding vertex polygon, which for node 2 is the square consisting of 4 right triangles, 1 each of area 1 , whereas for nodes 1 and 3 it consists of 4 right triangles of area 2 plus 2 three equilateral triangles, each of area 43 ; see Figure 14.12. The solution to the final linear system is easily found: c = ( 1.56724, 1.45028, 1.56724 ) . 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 14.15. A more accurate solution, based on a square grid triangulation of size h = .1 is plotted in the second figure. Inhomogeneous Boundary Conditions So far, we have restricted our attention to problems with homogeneous Dirichlet boundary conditions. The solution to the inhomogeneous Dirichlet problem - u = f in , u=h on ,
T
is also obtained by minimizing the Dirichlet functional (14.24). However, now the minimization takes place over the affine subspace consisting of all functions that satisfy the 3/15/06 255
c 2006 Peter J. Olver
Figure 14.15.
Finite Element Solutions to Poisson's Equation for an Oval Plate.
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 the interpolation condition (14.27), and so has the same piecewise affine form (14.35). The corresponding finite element approximation
m
w(x, y) =
i=1
ci i (x, y),
(14.54)
has the same form as before, (14.28), 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. (14.55)
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 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 (14.48), and the upper 3 3 subblock is the reduced matrix (14.49). The remaining entries of the 3/15/06 256
c 2006 Peter J. Olver
Figure 14.16.
Solution to the Dirichlet Problem for the Oval Plate.
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 (14.56) We similarly split the coefficients ci of the finite element function (14.54) 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 (14.55). The solution to the finite element approximation (14.54) is obtained by solving the associated linear system K c + K h = b, or K c = f = b - K h. (14.57)
Example 14.10. For the oval plate discussed in Example 14.8, 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 14.16. 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, -1, -1, -1, -1, 0, 1, 1, 1, 1, 0 ) . Using the previously computed formulae (14.49, 56) for the interior coefficient matrix K and boundary coefficient matrix K, we approximate the solution to the Laplace equation by solving (14.57). We are assuming that there is no external forcing function, f (x, y) 0, and so the right hand side is b = 0, and so we must solve K c = f = - K h = T ( 2.18564, 3.6, 7.64974 ) . The finite element function corresponding to the solution c = T ( 1.06795, 1.8, 2.53205 ) is plotted in the first illustration in Figure 14.16. 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. 3/15/06 257
c 2006 Peter J. Olver T T