The Finite Element Method FEM

This is a specific Galerkin method, and as such the trial and test spaces are the same and the aim is to make the residual orthogonal to all elements in the space. First, some references. A good beginning book is Claes Johnson's ``Numerical Solution of PDE's'' (1987), other books are Phillippe Ciarlet's ``FEM for Elliptic Problems'' (1978) and Brenner and Scott's ``FEM for PDE's'' (1999). Most of these books emphasize the mathematical aspects. However, a lot of work has been put on the implementation side, which is a crucial component of the method and is unfortunately glossed over in these more mathematical books. Some good references on implementation issues are: .

As the name implies, finite elements use basis sets that are compactly supported. There are some nice advantages of FEM:

(1) mathematical- you are forced to think from the outset what function spaces you are using and in what sense, then, is a computed solution ``close'' to the exact solution, i.e. the spaces have built-in norms. A second mathematical aspect is that in casting the problem mathematically, one is forced to seriously consider whether the partial differential equation, or ordinary differential equation, is well-posed. That is not to say that you can use the analysis of Galerkin and the proceed to implement the approximate solution by some non-Galerkin technique.

(2) computationally - the technique allows us to obtain the solution everywhere in the domain, not just at grid locations. This, in contrast to a spectral (a special collocation case) or finite-difference technique, where you are approximating the solution at a finite set of grid locations. That is not to say that you can circumvent this issue by careful interpolation, but this is already provided by the Galerkin technique (and more importantly, you know an error estimate everywhere in the domain). Perhaps the most important advantage is that the technique lends itself very naturally to tiling very complex domains with elements that fit more naturally than simple uniform or rectangular lattices. It is thus very popular in boundary value problems (and eigenvalue problems) that come from very complex structures, such as bridges, buildings,structural components of vehicles, etc. This has been made considerably easier to do lately, with the availability of very smart automated mesh generation packages (well, in 2d...but things are getting better in 3d), which take care of the most tedious and difficult part of the implementation of the method.

Again, here we focus only on elliptic problems with simple boundary conditions.

Dirichlet Model Problem Consider

$$- (a u) = f, \Omega(x,y)$$

$$u = 0, \partial\Omega.$$ (112)

Here $\Omega \subseteq \mathbb{R}^2$, $a, u, f: \Omega \rightarrow \mathbb{R}$. We assume that $0 \le \underline{a} \leq a(x,y) \leq \overline a$. We call this the Strong Formulation (SF).

A reminder:

$$u \equiv \frac{\partial u}{\partial x} \hat x + \frac{\partial u}{\partial y} \hat y,$$

$$- (a u) \equiv - \frac{\partial}{\partial x} (a \frac{\partial u}{\partial x}) - \frac{\partial }{\partial y} (a \frac{\partial u}{\partial y}).$$

Remark: note that when $a=1$, we get the standard Poisson Equation problem we've studied before.

Weak Formulation (WF): Find $u$, vanishing on $\delta \Omega$, such that

$$\int_\Omega a u \cdot v dx dy = \int_\Omega f v dx dy, \qquad \forall v \mbox{which vanish on $\delta \Omega$} .$$ (113)

Derivation: Use the divergence theorem (integration by parts in multi-dimensions). For ${\bf w}$ and $v$,

$$\int_\Omega ($$div$${\bf w}) v dx dy = -\int_\Omega {\bf w} \cdot$$   grad$$v dx dy + \int_{\delta \Omega} {\bf w}\cdot {\bf\hat n} v ds,$$

where ${\bf\hat n}$ is the outward normal to $\Omega$

Let ${\bf w} = a$   grad$u$ yields (113), using $v=0$ on $\delta \Omega$.

More precisely: for the SF we require $a(x,y)$ differentiable, and $u$ twice differentiable, i.e. $a \in C^1$ and $u \in C^2$. For the WF we require $a$ integrable and grad$u \cdot$   grad$v$ integrable.

A little detour: we remind ourselves that a function space $L^2(\Omega)$ is comprised of the set of functions $u$ for which

$$\int_\Omega \vert u\vert^2 d \Omega < \infty.$$

Example: Take the function $u=x^\beta$ and ask for what values of $\beta$ is $u$ an element of $L^2$ over the whole real line. If $x$ and $\beta$ are real,

$$\int_{-\infty}^\infty x^{2\beta} dx < \infty$$

if $\beta > -1/2$.

The Sobolev Space is a function space defined by the integrability of its elements. The Sobolev space $H^1(\Omega)$ is defined as

$$H^1(\Omega) = \left\{ u \in L^2(\Omega) \vert \qquad \partial_{x_i} u \in L^2(\Omega) \right\}$$

i.e. the function and all of its first derivatives are $L^2(\Omega)$. So, looking at spaces, $L^2 \supset H^1 \supset C^2 \supset C^1$. Another space, which we will make use of:

$$H_0^1(\Omega) =\left\{ u \in H^1(\Omega) \vert \qquad u = 0 \qquad \mbox{ on } \delta \Omega \right\}$$

Theorem: If $u \in C^2$ and $u$ satisfies the SF then $u$ satisfies the WF. The proof is above.

Theorem: The WF has a unique solution.

The Variational Problem (V): Find $u \in H^1_0(\Omega)$ such that

$$E(u) = \min_{v \in H^1_0} E(v),$$

where

$$E(v) = \frac{1}{2} \int_{\Omega} a \vert$$   grad$$v\vert^2 d \Omega - \int_\Omega f v d \Omega,$$

so that the functional $E:H^1_0 \rightarrow \mathbf{R}$.

Theorem $u$ satisfies the WF if and only if $u$ satisfies the variational problem.

Proof: suppose $u$ satisfies the WF. Let's evaluate

$$E(v) = E(u+(v-u)) = \frac{1}{2} \int a \vert$$   grad$$u +$$   grad$$(v-u)\vert^2 - \int f(u+(v-u)).$$

Note that the integral sign assumes the integration is over the whole domain $\Omega$. Expanding the first integral

$$E(v) = E(u) + \frac{1}{2} \int a \vert$$   grad$$(v-u)\vert^2 + \int a$$   grad$$u \cdot$$   grad$$(v-u) - \int f(v-u).$$

The last two terms are zero since WF holds. Then

$$E(v) = E(u) + \frac{1}{2} \int a \vert$$grad$$(v-u)\vert^2 \geq E(u)$$

which implies the VF. $\Box$

Remark: For non-self-adjoint problems, can have WF but no VF. Note, however, that any linear second order SF problem can be put in self adjoint form (see Sturm-Liouville theory).

The problem

$$-$$   div$$\cdot (a$$   grad$$u) + {\bf b}\cdot$$   grad$$u + c u = f$$

with zero Dirichlet boundary condtions has a WF but does not have a VF unless ${\bf b}$ is zero.

Exercise Let

$$b(u,v) = \int [a$$   grad$$u \cdot$$   grad$$v + {\bf b}\cdot ($$grad$$u) v + c u v ]$$

Let $E(w) = \frac{1}{2} b(w,w) - F(w)$. Use $w=u+tv$, i.e. thinking of $tv$ as a perturbation from $u$. Find the stationary condition for $E$ about $t=0$. What can you conclude from this calculation?

$\Box$

Theorem: If $u$ satisfies the VF, then $u$ satisfies the WF.

Proof: Choose any $v\in H^1_0$, and consider a real quantity $t$ and $\rho: \mathbf{R} \rightarrow \mathbf{R}$, defined as

$$\rho(t) = E(u+t v) = \frac{1}{2} \int a \vert$$   grad$$u + t$$   grad$$v\vert^2 - \int f(u+tv)$$

if $\rho(0) = \min \rho(t)$, the $\rho'(0) = 0$.

$$0 = \rho'(0) = \int$$   grad$$u \cdot$$   grad$$v - \int f v,$$

which is the WF.

In what follows we will define

$$b(u,v) \equiv \int a u \cdot v$$

$$F(v) \equiv \int f v,$$

so that $E(u) = \frac{1}{2} b(u,u) - F(u)$.

The numerical approximation of the VF problem is called the Rayleigh-Ritz Method and it can be summarized as follows: Choose a subspace $S_h \in H^1_0$, where $S_h$ is finite dimensional subspace. Then, finding the $u_h \in S_h$ that minimizes $E(v)$ yields an approximation to $u_h \approx u$, if $S_h$ is ``sufficiently large.'' This could be called the VF$_h$ problem.

The numerical approximation of the WF problem can be done using Galerkin techniques. Briefly, use a finite dimensional set of basis for $S_h$. Supposing dim$S_h = N$, then, every function $v$ in $S_h$ can be written as

$$v = \sum_{i=1}^N \alpha_i \phi_i.$$

where $\phi$ are the bases. Substituting in WF gives

$$b(u_h,v) = F(v)$$

or

$$b(\sum_{i=1}^N \alpha_i \phi_i,\phi_j) = F(\phi_j) \qquad i,j = 1,2,\ldots N$$

which is thus

$$\sum_{i=1}^N \alpha_i b(\phi_i,\phi_j) = F(\phi_j) \qquad i,j = 1,2,\ldots N.$$

$b(\phi_i,\phi_j)$ with $i,j = 1,2,\ldots N$ is a matrix, we shall call it the stiffness matrix $M$, and $F(\phi_j)$, with $j = 1,2,\ldots N$ is a vector we shall call the load vector ${\bf F}$. So in matrix notation, the task is to solve

$$M {\bf\alpha} = {\bf F}$$

for the unknown vector ${\bf\alpha}$.

What happens if the problem is instead nonlinear? for example, suppose we want to solve

$$-$$div$$\cdot (a(x,u)$$   grad$$u(x)) = f(x) \qquad x \in \Omega$$

subject to the zero Dirichlet boundary conditions, with the usual assumptions on $a(x,u)$. In that case we are lead to a nonlinear set of equations, given by

$$\int a(x,\sum_j \alpha_j \phi_j)$$   grad$$\phi_j \cdot$$   grad$$\phi_i = \int f \phi_i$$

with $i,j=1,2,\ldots,N$. We know how to solve these types of systems using a root finding algorithm (For Newton's method, see).

Exercise Set up the WF problem and its Galerkin approximation for

$$-u''(x) = f \qquad 0 < x < 1,$$ (114)

with boundary conditions $u(0)=u(1) = 0$. Let $S_h = {\cal P}^9$, where the basis set is $\phi_i = x^i(1-x)$, $i = 0,1,\ldots 9$. Prove that the underlying linear algebraic problem has a stiffness matrix that is symmetric and nonsingular.

$\Box$

Model Problem with Neumann Boundary Conditions

$$- (a u) = f \Omega$$ (115)

$$a \frac{\partial u}{\partial n} = \partial \Omega$$ (116)

If $u \in H^m$ then the partial derivatives up to order $m-1$ are defined on $\partial \Omega$ and are $L^2$ functions there. However, partials of order $m$ cannot be sensibly defined on $\partial \Omega$. So to define a space of functions that is $H^1$ such that $a \partial/\partial n$ of these functions are zero on the boundary is nonsensical.

Instead, let $v \in H^1$.We show that this is a good function space for the trial functions:

$$-\int_\Omega$$   div$$\cdot (a$$   grad$$u) v = \int_\Omega a$$   grad$$u \cdot$$   grad$$v - \int_{\partial \Omega} (a$$   grad$$u ) \cdot {\hat n} v = \int_\Omega f v$$

which is

$$\int_\Omega a$$   grad$$u \cdot$$   grad$$v - \int_{\partial \Omega} a \frac{\partial u}{\partial n} = \int_\Omega f v.$$

If $u$ solves (116), then the second term on the left hand side is zero. Hence, the WF of the problem is

$$\int_\Omega a$$   grad$$u \cdot$$   grad$$v = \int_\Omega f v, \qquad \forall v \in H^1.$$

Remark: we say that $a \frac{\partial u}{\partial n} = 0$ is a natural boundary condition, whereas $u=0$ is an essential boundary condition.

Exercise Set up the WF problem and its Galerkin approximation for

$$-u''(x) + a u'(x) + b u(x) = f \qquad 0 < x < 1,$$

$$u'(x) = c_0 \qquad x=0$$

$$u'(x) = c_1 \qquad x=1$$

Let $S_h \subset H^1$ be the set of $N$ hat piece-wise linear functions defined on the distinct nodes $x_i$, $i=1,2,\ldots N$ defined on the interval $0 \leq x \leq 1$. Note: you should get that $c_1$ and $c_2$ make contributions to the load vector.

$\Box$

An Error Estimate: for the problem posed in (114): take $S_h \subset H^1_0$, then

$$b(u_h,v) = (f,v) \qquad \forall v \in S_h.$$

and

$$b(u,v) = (f,v) \qquad \forall v \in S_h.$$

where

$$u_h = \sum_{i=1}^{N} \alpha_i \phi_i$$

where the set $\{\phi_i \}_{i=1}^N$ span $S_h$. Then

$$b(u-u_h,v) = 0 \qquad \forall v \in S_h.$$

Theorem For any $v \in S_h$, as above,

$$\vert\vert(u-u_h)'\vert\vert \leq \vert\vert(u-v)'\vert\vert.$$

Proof: Let $v,u_h \in S_h$, and $w=u_h-v$. Then

$$\vert\vert (u - u_h)'\vert\vert^2 = b(u-u_h,u-u_h) + b(u-u_h,w)$$

$$= b(u-u_h,u-u_h+w) = b(u-u_h,u-v) \leq \vert\vert (u-u_h)' \vert\vert \vert\vert(u-v)'\vert\vert$$

by Cauchy Schwarz. $\Box$

Furthermore, if $\phi_i$ are piecelinear hat functions,

$$\vert u - u_h\vert \leq \frac{h^2}{8} \max_{y \in \Omega} \vert u''(y)\vert$$

where $h$ is the largest distance between the nodes. This is the same estimate you would get for the finite difference approximation of the equation using center difference approximation for the second derivative.

$\Box$

Worked Two-Dimensional Example Calculation Here we will work through a two-dimensional calculation. We will assume that you will be doing all of the work. We'll provide you with some answers for what you should get along the way.

We'll solve

$$\nabla^2 u = f(x,y)\qquad (x,y) \in \Omega$$

$$u(x,y)=\alpha(x,y) \qquad (x,y)\in\partial\Omega_1$$

$$\frac{\partial u}{\partial n(x,y)} = \beta(x,y) \qquad (x,y)\in\partial\Omega_2$$

See Figure 18.

Figure 18: (a) Domain for worked example; (b) node and element assignment.

We will assume in this example that $f=0$, $\alpha = 0$, and $\beta(\theta) = \sin 2 \theta$. We will take the inner radius to be $r=1$ and the outer one to be $r=2$. Here $x=r \cos \theta$, and $y= r \sin \theta$. The Laplacian in polar coordinates, applied to some function $w$, reads

$$\frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial w}{\partial r} \right) + \frac{1}{r^2} \frac{\partial^2 w}{\partial \theta^2}.$$

We will tile the domain with 8 triangular elements. Each node point in the figure is numbered.

To label node points we will use a local and a global scheme. A single index will denote a node in the global scheme as given in the figure. A double index will denote the local scheme. In this local scheme the node points are numbered $k1$, $k2$, $k3$ going in counterclockwise direction where $k$ is the element number. For example, in element V the local nodes 51, 52, 53 correspond to the global nodes 2, 7, 8, respectively. We will use piecewise bilinear elements. hence our approximation will have $C^0$ global continuity.

Construction of the elements Let

$$\phi_{ki}(x,y)= \gamma_1^i + \gamma_2^i x + \gamma_3^i y$$

$$\phi_{ki}(x_{ki},y_{ki})=1\qquad \phi_{ki}(x_{kj},y_{kj})=0\quad i\neq j$$

Take $i=1$

$$\begin{split}\phi_{k1}(x_1,y_1)&=\gamma_1^1+\gamma_2^1x_1+\ga... ...,y_3)&=\gamma_1^1+\gamma_2^1x_3+\gamma_3^1y_3 = 0 \end{split}$$

Solving gives

$$\begin{split}\gamma_1^1 &= \frac{x_2y_3-x_3y_2}{2A} \gamma_... ...c{y_2-y_3}{2A} \gamma_1^3 &= \frac{x_3-x_2}{2A} \end{split}$$

where

$$A=\frac{1}{2} \det \left\vert \begin{array}{ccc} 1 & 1 & 1 \... ..._1 & x_2 & x_3 \\ y_1 & y_2 & y_3 \\ \end{array}\right\vert$$

i.e. the area of the triangle.

The other two basis functions that are nonzero in the $k$-th element are

$$\phi_{ki}(x,y)=\frac{1}{2A}\big[ (x_{k(i+1)}y_{k(i+2)}-x_{k(i+2)}y_{k(i+1)}+(y_{k(i+1)}-y_{k(i+2)})x + (x_{k(i+2)}-x_{k(i+1)})y \big]$$

where $i=2,3$ and the indicial operations are taken $\mod 3$ with counterclockwise local numbering assumed.

We want an approximation of the solution in the form

$$v(x,y)=\sum_{j=1}^{N} v_j \phi_j(x,y)$$ (117)

where $N$ is the number of node points. Again we multiply the given differential equation by $\phi_i$ and integrate. This way we get

$$\int_{\Omega}\nabla^2 v \phi_i(x,y) \mathrm{d}x \mathrm{d}y = \int_{\Omega} f\phi_i(x,y) \mathrm{d}x \mathrm{d}y$$ (118)

Integration by parts produces

$$\int_{\Omega}\nabla^2 v \phi_i(x,y) \mathrm{d}x \mathrm{d}y = ... ... \mathrm{d}s - \int_{\Omega}\nabla v \nabla \phi_i \mathrm{d}x \mathrm{d}y$$ (119)

where $\hat{n}$ is the outward unit normal on each segment of $\partial \Omega$ and $\hat{n}\cdot\nabla v=\frac{\partial v}{\partial n}$.

Using (118) and (120) we rewrite (119) as follows

$$\int_{\partial\Omega}\frac{\partial v}{\partial n} \phi_i \mathr... ... \mathrm{d}x \mathrm{d}y = \int_{\Omega} f\phi_i \mathrm{d}x \mathrm{d}y$$ (120)

since $f=0$, the RHS vanishes. The

$$\int_{\Omega} \nabla \phi_j \cdot \nabla \phi_i \mathrm{d}x \mathrm{d}y$$

terms give entries into what is commonly called the stiffness matrix $M$, a $N\times N$ square matrix. The quantities

$$\int_{\partial\Omega}\frac{\partial v}{\partial n} \phi_i \mathrm{d}s - \int_{\Omega} f\phi_i \mathrm{d}x \mathrm{d}y$$

are called load vector $F$. So if $V=[v_1,v_2,\ldots,v_N]^T$ then (121) is nothing more than

$$MV=F$$

We most efficiently compute the stiffness matrix by using the relation

$$\int_{\Omega} = \sum \int_{\Omega_k}$$

where $\Omega_k$ is the domain of each element. The plan is finding the contributions for the stiffness matrix from every element and then adding up all these contributions. So we need to evaluate

$$\int_{\Omega^e} \nabla \phi_l \cdot \nabla \phi_m \mathrm{d}x \... ...la x} +\frac{\partial\phi_l}{\nabla y}\frac{\partial\phi_m}{\nabla y} \right)$$

$l,m=1,2,3$. You should be able to show that

$$\begin{multline*} \int_{\Omega^e} \nabla \phi_l \cdot \nabla \phi_m \mathrm{d}... ...+2})(y_{m+1}-y_{m+2}) +(x_{l+2}-x_{l+1})(x_{m+2}-x_{m+1}) \Big] \end{multline*}$$

The load vector $\int \partial v/\partial n \mathrm{d}s$ can be found as follows: It is easy to see that it vanishes for all interior nodes because those $\phi_i$ are zero along the boundary $\partial \Omega$. Let us look at node 2. Assume that $\partial v/\partial n = 1$ and use the fact that $\phi_2$ is linear on $\partial \Omega$. Then

$$\begin{split}\int \frac{\partial v}{\partial n} \phi_2 \mat... ...\mathrm{d}s &=\frac{l_{1-2}}{2}+\frac{l_{2-3}}{2} \end{split}$$

where $l_{i-j}$ denotes the distance from node $i$ to node $j$.

To solve the given problem, follow these steps:

1. Find the coordinates of all the node points.
   

   node number	coordinates
   1	(0, 2)
   2	(1.4142, 1.4142)
   3	(2, 0)
   4	(1, 0)
   5	(0.7071, 0.7071)
   6	(0, 1)
   7	(0.75, 1.299)
   8	(1.299, 0.75)

2. Make a table assigning to each element its corresponding node points in counterclockwise motion.

   element number	nodes
   I	1,7,2
   II	1,6,7
   III	6,5,7
   IV	2,7,8
   V	5,8,7
   VI	3,2,8
   VII	5,4,8
   VIII	4,3,8

3. Find the local stiffness matrices for each element. Use symmetries and a matlab script to make the task less tedious.

   As all stiffness matrices are symmetric only the upperdiagonal part is given. Element I Element II
   

   0.4116	-0.7897	0.3781
   	2.1224	-1.3327
   		0.9546

   0.4346	-0.2353	-0.1993
   	0.7026	-0.4673
   		0.6667

   
   
   

   Element III Element IV
   

   0.4086	-0.2426	-0.1659
   	0.7563	-0.5136
   		0.6796

   0.7045	-0.3523	-0.3523
   	0.5311	-0.1789
   		0.5311

   
   
   

   Element V Element VI
   

   0.8646	-0.4323	-0.4323
   	0.5051	-0.0728
   		0.5051

   0.4116	0.3781	-0.7897
   	0.9546	-1.3327
   		2.1224

   
   
   

   Element VII Element VIII
   

   0.7563	-0.2426	-0.5136
   	0.4086	-0.1659
   		0.6796

   0.7026	-0.2353	-0.4673
   	0.4346	-0.1993
   		0.6667

   
   
   

4. Add up the element stiffness matrices to obtain the global stiffness matrix.
   
   

   0.8462	0.3781	0	0	0	-0.2353	0.989	0
   	2.6137	0.3781	0	0	0	-1.685	-1.685
   		0.8462	-0.2353	0	0	0	0.989
   			1.1112	-0.2426	0	0	-0.6332
   				2.3772	-0.2426	-0.9459	-0.9549
   					1.1112	-0.6332	0
   						4.5059	-0.2517
   							4.5049

   
   
   You should do this on your own. Click here for a possible computer program solution. e) Find the force vector.

   The only contribution comes from the second node. We will integrate the boundary terms over the arc of the circle. If many points on the arc are used the difference between integration along the arc to integration along the triangle edges is small and the integrals are much more convenient along the arcs. Let us just consider element VI which will give us $F_2/2$ because of symmetry. Using polar coordinates we have $\mathrm{d}s = 2 \mathrm{d}\theta$ and we can interpolate $\phi_2$ linearly which will gives us

   $$\frac{F_2}{2}=\int_{\text{node 2}}^{\text{node 3}} 1\cdot\phi_2 ... ...rac{\pi}{4}} \left(\frac{\pi}{4}-\theta\right)\sin(2\theta) \mathrm{d} \theta$$

5. Impose the Dirichlet boundary conditons.

   All nodes except 2, 7, 8 have zero function value. Therefore the stiffness matrix reduces to a $3\times3$ matrix and we end up with the following system of equations:

   $$\left(\begin{array}{ccc} 2.6137 & -1.685 & -1.685\\ -1.685 & 4.5... ...end{array}\right) = \left(\begin{array}{c} 0.5708 \\ 0 \\ 0 \end{array}\right)$$

6. Solve the resulting system of equations.

   We obtain $v_1 = 0.4464$, $v_2=0.1769$, $v_3=0.1769$.

7. Compare the results with the exact solution. The exact solution is $u(r,\theta) = 4/17 (r^2 - 1/r^2) \sin 2 \theta$.

Our approximation is, not surprisingly, poor. However, only three free nodes were included in the calculation. More elements are needed for a better result. $\Box$

$\Box$