Finite Element Method Basics

The core Partial Differential Equation Toolbox™ algorithm uses the Finite Element Method (FEM) for problems defined on bounded domains in 2-D or 3-D space. In most cases, elementary functions cannot express the solutions of even simple PDEs on complicated geometries. The finite element method describes a complicated geometry as a collection of subdomains by generating a mesh on the geometry. For example, you can approximate the computational domain Ω with a union of triangles (2-D geometry) or tetrahedra (3-D geometry). The subdomains form a mesh, and each vertex is called a node. The next step is to approximate the original PDE problem on each subdomain by using simpler equations.

For example, consider the basic elliptic equation.

$$-\nabla \cdot \left(c\nabla u\right)+au=f\text{ on domain }\Omega$$

Suppose that this equation is a subject to the Dirichlet boundary condition $$u=r$$ on $$\partial {\Omega }_D$$ and Neumann boundary conditions on $$\partial {\Omega }_N$$. Here, $$\partial \Omega =\partial {\Omega }_D\cup \partial {\Omega }_N$$ is the boundary of Ω.

The first step in FEM is to convert the original differential (strong) form of the PDE into an integral (weak) form by multiplying with test function $$v$$ and integrating over the domain Ω.

$${\displaystyle \underset{\Omega }{\int }\left(-\nabla ·\left(c\nabla u\right)+au-f\right)vd\Omega }=0\forall v$$

The test functions are chosen from a collection of functions (functional space) that vanish on the Dirichlet portion of the boundary, $$v=0$$ on $$\partial {\Omega }_D$$. Above equation can be thought of as weighted averaging of the residue using all possible weighting functions $$v$$. The collection of functions that are admissible solutions, u, of the weak form of PDE are chosen so that they satisfy the Dirichlet BC, u = r on $$\partial {\Omega }_D$$.

Integrating by parts (Green’s formula) the second-order term results in:

$${\displaystyle \underset{\Omega }{\int }\left(c\nabla u\nabla v+auv\right)d\Omega }-{\displaystyle \underset{\partial {\Omega }_N}{\int }\overrightarrow{n}·\left(c\nabla u\right)}vd\partial {\Omega }_N+{\displaystyle \underset{\partial {\Omega }_D}{\int }\overrightarrow{n}·\left(c\nabla u\right)}vd\partial {\Omega }_D={\displaystyle \underset{\Omega }{\int }fvd\Omega }\forall v$$

Use the Neumann boundary condition to substitute for second term on the left side of the equation. Also, note that $$v=0$$ on $$\partial {\Omega }_D$$ nullifies the third term. The resulting equation is:

$${\displaystyle \underset{\Omega }{\int }\left(c\nabla u\nabla v+auv\right)d\Omega }+{\displaystyle \underset{\partial {\Omega }_N}{\int }quv}d\partial {\Omega }_N={\displaystyle \underset{\partial {\Omega }_N}{\int }gv}d\partial {\Omega }_N+{\displaystyle \underset{\Omega }{\int }fvd\Omega }\forall v$$

Note that all manipulations up to this stage are performed on continuum Ω, the global domain of the problem. Therefore, the collection of admissible functions and trial functions span infinite-dimensional functional spaces. Next step is to discretize the weak form by subdividing Ω into smaller subdomains or elements $${\Omega }^e$$, where $$\Omega =\cup {\Omega }^e$$. This step is equivalent to projection of the weak form of PDEs onto a finite-dimensional subspace. Using the notations $$u_h$$ and $$v_h$$ to represent the finite-dimensional equivalent of admissible and trial functions defined on $${\Omega }^e$$, you can write the discretized weak form of the PDE as:

$${\displaystyle \underset{{\Omega }^e}{\int }\left(c\nabla u_h\nabla v_h+a{u}_h{v}_h\right)d{\Omega }^e}+{\displaystyle \underset{\partial {\Omega }_N^e}{\int }q{u}_{h}v{}_{h}d\partial {\Omega }_N^e}={\displaystyle \underset{\partial {\Omega }_N^e}{\int }gv{}_{h}d\partial {\Omega }_N^e}+{\displaystyle \underset{{\Omega }^e}{\int }f{v}_{h}d{\Omega }^e}\forall v_h$$

Next, let ϕ_i, with i = 1, 2, ... , N_p, be the piecewise polynomial basis functions for the subspace containing the collections $$u_h$$ and $$v_h$$, then any particular $$u_h$$ can be expressed as a linear combination of basis functions:

$$u_h={\displaystyle \sum _1^{N_p}{U}_i{\varphi }_i}$$

Here U_i are yet undetermined scalar coefficients. Substituting $$u_h$$ into to the discretized weak form of PDE and using each $$v_h={\phi }_i$$ as test functions and performing integration over element yields a system of N_p equations in terms of N_p unknowns U_i.

Note that finite element method approximates a solution by minimizing the associated error function. The minimizing process automatically finds the linear combination of basis functions which is closest to the solution u.

FEM yields a system KU = F where the matrix K and the right side F contain integrals in terms of the test functions ϕ_i, ϕ_j, and the coefficients c, a, f, q, and g defining the problem. The solution vector U contains the expansion coefficients of u_h, which are also the values of u_h at each node x_k (k = 1,2 for a 2-D problem or k = 1,2,3 for a 3-D problem) since u_h(x_k) = U_i.

FEM techniques are also used to solve more general problems, such as:

• Time-dependent problems. The solution u(x,t) of the equation

  $$d\frac{\partial u}{\partial t}-\nabla \cdot \left(c\nabla u\right)+au=f$$

  can be approximated by

  $$u_h(x,t)={\displaystyle \sum _{i=1}^N{U}_i(t){\varphi }_i}(x)$$

  The result is a system of ordinary differential equations (ODEs)

  $$M\frac{dU}{dt}+KU=F$$

  Two time derivatives result in a second-order ODE

  $$M\frac{d^{2}U}{d{t}^2}+KU=F$$

• Eigenvalue problems. Solve

  $$-\nabla \cdot \left(c\nabla u\right)+au=\lambda du$$

  for the unknowns u and λ, where λ is a complex number. Using the FEM discretization, you solve the algebraic eigenvalue problem KU = λMU to find u_h as an approximation to u. To solve eigenvalue problems, use solvepdeeig.

• Nonlinear problems. If the coefficients c, a, f, q, or g are functions of u or ∇u, the PDE is called nonlinear and FEM yields a nonlinear system K(U)U = F(U).

To summarize, the FEM approach:

1. Represents the original domain of the problem as a collection of elements.

2. For each element, substitutes the original PDE problem by a set of simple equations that locally approximate the original equations. Applies boundary conditions for boundaries of each element. For stationary linear problems where the coefficients do not depend on the solution or its gradient, the result is a linear system of equations. For stationary problems where the coefficients depend on the solution or its gradient, the result is a system of nonlinear equations. For time-dependent problems, the result is a set of ODEs.

3. Assembles the resulting equations and boundary conditions into a global system of equations that models the entire problem.

4. Solves the resulting system of algebraic equations or ODEs using linear solvers or numerical integration, respectively. The toolbox internally calls appropriate MATLAB^® solvers for this task.

References

[1] Cook, Robert D., David S. Malkus, and Michael E. Plesha. Concepts and Applications of Finite Element Analysis. 3rd edition. New York, NY: John Wiley & Sons, 1989.

[2] Gilbert Strang and George Fix. An Analysis of the Finite Element Method. 2nd edition. Wellesley, MA: Wellesley-Cambridge Press, 2008.