PE281 Finite Element Method Course Notes

PE281 Finite Element Method Course Notessummarized by Tara LaForceStanford, CA 23rd May 20061Derivation of the MethodIn order to derive the fundamental concepts of FEM we will start by lookingat an extremely simple ODE and approximate it using FEM.1.1The Model ProblemThe model problem is: u′′ u x 0 x 1u (0) 0u (1) 0(1)and this problem can be solved analytically: u (x) x sinhx/sinh1. Thepurpose of starting with this problem is to demonstrate the fundamentalconcepts and pitfalls in FEM in a situation where we know the correct answer,so that we will know where our approximation is good and where it is poor.In cases of practical interest we will look at ODEs and PDEs that are toocomplex to be solved analytically.FEM doesn’t actually approximate the original equation, but rather the weakform of the original equation. The purpose of the weak form is to satisfythe equation in the "average sense," so that we can approximate solutionsthat are discontinuous or otherwise poorly behaved. If a function u(x) is asolution to the original form of the ODE, then it also satisfies the weak formof the ODE. The weak form of Eq. 1 is1

Z1′′( u u) vdx 0Z1xvdx(2)0The function v(x) is called the weight function or test function. v(x) can beany function of x that is sufficiently well behaved for the integrals to exist.The set of all functions v that also have v(0) 0, v(1) 0 are denoted byH. (We will put many more constraints on v shortly.)The new problem is to find u so thatR10( u′′ u x) vdx 0 f orallv Hu (0) 0u (1) 0(3)Once the problem is written in this way we can say that the solution u belongsto the class of trial functions which are denoted H̃. When the problem iswritten in this way the classes of test functions H and trial functions H̃are not the same.R 1 ′′For example, u must be twice differentiable and have theproperty that 0 u vdx , while v doesn’t even have to be continuous aslong as the integral in Eq. 3 exists and is finite. It is possible to approximateu in this way, but having to work with two different classes of functionsunnecessarily complicates the problem. In order to make sure that H and H̃are the same we can observe that if v is sufficiently smooth thenZ1′′ u vdx 0Z01u′ v ′ dx u′ v 10(4)This formulation must be valid since u must be twice differentiable and v wasarbitrary. This puts another constraint on v that it must be differentiableand thatRthose derivatives must be well-enough behaved to ensure that the1integral 0 u′ v ′ dx exists. Moreover, since we decided from the outset thatv(0) 0 and v(1) 0, the second term in Eq. 4 is zero regardless of thebehavior of u′ at these points. The new problem isZ10(u′ v ′ uv xv) dx 0.(5)Notice that by performing the integration by parts we restricted the class oftest functions H by introducing v ′ into the equation. We have simultaneouslyexpanded the class of trial functions H̃, since u is no longer required to have2

a second derivative in Eq. 5. The weak formulation defined in Eq. 5 is calleda variational boundary-value problem.In Eq. 5 u and v have exactly the same constraints on them:R1R11. u and v must be square integrable, that is: 0 uvdx 0 u2 dx 2. ThederivativesR 1 ′ firstR 1 ′ 2 of u and v must be square integrable, that is:′uvdx (u ) dx (this actually guarantees the first property)003. We had already assumed that v(0) 0 and v(1) 0 and we know fromthe original statement of the problem that u(0) 0 and u(1) 0.Now we have that H̃ H H01 . Any function w is a member of H01R1if 0 (u′ )2 dx and w(0) w(1) 0. H01 is the space of admissiblefunctions for the variational boundary-value problem (ie. all admissible testand trial functions are in H01 )We will consider the variational form Eq. 5 to be the equation that we wouldlike to approximate, rather than the original statement in Eq. 1. Once wehave found a solution to Eq. 5 in this way we can ask the question whetherthis formulation is also a solution to Eq. 1: That is, whether this solution isa function satisfying Eq. 1 at every x in 0 x 1, or whether we have founda solution that satisfies only the weak form of the equation. In the case thatwe can only find a solution to the weak form, no "classical" solution exists.1.2Galerkin ApproximationsWe now have the problem re-stated so that we are looking for u H01 suchthatZ1′ ′(u v uv) dx 0Z1xvdx(6)0for all v H01 . In order to narrow down the number of functions we willconsider in our approximate solutions we will make two more assumptionsabout H01 . First, we will assume that H01 is a linear space of functions (thatis if v1 , v2 H01 and a, b are constants then av1 bv2 H01 .)1The second assumption is that H0 is infinite dimensional. For example1 if wehave the sine series ψn (x) 2sin(nπx) for n 1, 2, 3, . and v H0 then3

Pv can be represented by v (x) n 1 an ψn (x). The scalar coefficients anR1are given by an 0 v (x) ψn (x) dx, just like usual. Hence infinititely manycoefficients an must be found to define v exactly. As in Fourier analysis,many of these coefficients will be zero. We will also truncate the series inorder to have managable length series, just like in discrete Fourier analysis.Unlike in Fourier analysis, though the basis functions do not have to be sinesand cosines, much less smooth functions can be used. In fact our set of basisfunctions do not even have to be smooth and can contain discontinuities inthe derivatives, but they must be continuous. We will assume that the infiniteseries converges so that we can consider only the first N basis functions andget a good approximation vN of the original test (or trial) function:v vN XNi 1(7)βi φi (x)where φi are as-yet unspecified basis functions. This subspace of functions is(N )denoted H0 and is a subspace of H01 . Galerkin’s method consists of finding(N )an approximate solution to Eq. 6 in a finite-dimensional subspace H0 of(1H0 of admissible functionsrather than in the whole space H01 . Now we arePNlooking for uN i 1 αi φi (x). The new approximate problem we have is(N )to find uN H0 such thatZ1′(u′N vN uN vN ) dx 0Z1(8)xvN dx0(N )for all vN H0 . Since the φi are known (in principle) uN will be completelydetermined once the coefficients αi have been found.PPNIn order to find that αn we put Ni 1 αi φi (x) andi 1 βi φi (x) into Eq. 8.hPi hPi NNdd Z 1 i 1 βi φi (x) dxj 1 αj φi (x) hdxi hPiPNNβφ(x)αφ(x) iijii 1j 10 PNx i 1 βi φi (x)for all N independent sets of βi .This can be expanded and factored to give4 dx 0(9)

XNi 1βi XNj 1 Z01 φ′j(x) φ′i Z(x) φj (x) φi (x) dx αj 01 xφi (x) dx 0(10)for all N independent sets of βi . The structure of Eq. 10 is easier to see if itis re-written asXNi 1βi XNj 1 Kij αj Fi 0(11)for all βi . WhereKij R1 0φ′j(x) φ′i (x) φj (x) φi (x) dx F Z1xφi (x) dx(12)0and where i, j 1, ., N. The N N matrix of Kij is called the stiffnessmatrix and the vector F is the load vector. Since the βi are known Kij andF can be calculated directly. But the βi were arbitrary so we can choose eachelement βi for eachP equation. For the first equation choose β1 1 and βn 0for n 6 1. Now Nj 1 K1j αj F1 . Similarly for the second equation choosePβ2 1 and βn 0 for n 6 2 so that Nj 1 K2j αj F2 . In this way we havechosen N independent equations that can be used to find thePNN unknownsαi . Moreover the N coefficients αi can be found from αj j 1 (K 1 )jiFiwhere (K 1 )ji are the elements of the inverse of K.The stiffness matrix K is symmetric for this simple problem, which makesthe computation of the matrix faster since we don’t have to compute all ofthe elements, symmetric matricies are also much faster to invert.1.3Finite Elements Basis FunctionsNow we have done a great deal of work, but it may not seem like we aremuch closer to finding a solution to the original ODE since we still knownothing about φi . The purpose of using such a general formulation is thatany set of linearly independent functions will work to solve the ODE. Nowwe are finally going to talk about what kind of functions we will want to useas basis functions. The finite element method is a general and systematictechnique for constructing basis functions for Galerkin approximations. In5

FEM the basis functions φi are defined piecewise over subregions. Over anysubdomain the φi will be chosen to be polynomials of low degree, thoughother possibilities do exist. finite elements are the subregions of the domain over which each basisfunction is defined. Hence each basis function has compact support overan element. Each element has length h. The lengths of the elementsdo NOT need to be the same (but generally we will assume that theyare.) nodes or nodal points are defined within each element. In Figure 1 thefive nodes are the endpoints of each element (numbered 0 to 4). the finite element mesh is the collection of elements and nodal pointsthat make up the domain and is shown in Figure 1. An element i isdenoted by Ωi .Now we need to construct the actual basis functions using the three criteriadefined before: 1) The basis functions are simple functions defined piecewiseover the finite element mesh, 2) the basis functions must be in the class of testfunctions H01 , and 3) The basis functions are chosen so that the parametersαi are the values of uN (x) at the nodal points.The simplest set of basis functions are the “hat functions” on elements i 1, 2, 3.φi (x) x xi 1hixi 1 xhi 10 forxi 1 x xi forxi x xi 1 for x xi 1 , x xi 1(13)where hi xi xi 1 is the length of element i. The derivatives areφ′i (x) 1hi 1hi 10 forxi 1 x xi forxi x xi 1 for x xi 1 , x xi 1(14)The equations for elements 0 and 4 have been left out since we decided thatu(0) u(1) 1, so no basis functions are required. In general the basisfunctions for the first and last elements are half of the functions since there6

is no i 1 or i 1 node, respectively. The hat functions are shown in Figure 2.The mathematical term for hat functions is piece-wise linear basis functionsΩ1h1h20Ω4Ω3Ω2h31h4234Figure 1: Four finite elements on the interval [0 1].1φ1Ω11φ3φ2Ω4Ω3Ω2φ 1Ω10φ 3φ 2Ω4Ω3Ω20h1h2h3h1h4h2h3h4 1 10123401234Figure 2: Four hat functions (top) and their derivatives (bottom) on the interval [0 1].Looking at the three criteria above, clearly the functions in Eq. 13 are simpleand defined element-wise. It is easy to show that they are in H01 , since theyhave square-integrable first derivatives. They also satisfy the third criteriasince φi (xj ) 1 if i j and 0 otherwise. Hence each function contributes tothe value of uN at exactly one node and αi uN (xi ).It is less clear that the hat functions will give a continuous representationof vN and uN . Let v be the sine function with period 2 shown in Figure7

3. At the nodes (0, 1, 2, 3, 4) sine has the values (0,0.7071,1,0.7071,0).The representation vN on the finite element mesh is vN 0.7071φ1 (x) φ2 (x) 0.7071φ3 (x). When the elements are summed up the sine waveis approximated by piecewise linear functions between each of the nodes,and is exactly represented at each node. When more nodes are used theapproximation improves and in the limit of N the sine wave would beexactly represented. In FEM we will never proceed all the way to the limit,so the interval size h will always have finite size h. This is why the termfinite elements is used.1vNsin(pi x)φ20.7071φ30.7071φ1Ω1Ω4Ω3Ω20h1h2h3h4 100.511.522.533.54Figure 3: The finite element approximation of sin(πx) using five nodes on the interval[0 1].1.4The Stiffness Matrix K and the Load Vector F forHat FunctionsRecall from Eq. 12 that each element of the stiffness matrix K is given by R1Kij 0R φ′i (x) φ′j (x) φi (x) φj (x) dx P 4e 1 Ωe φ′i (x) φ′j (x) φi (x) φj (x) dxP 4e 1 Kije8(15)