What is the difference between Finite Difference Methods, Finite Element Methods and Finite Volume Methods for solving PDEs?

1. Can you help me explain the basic difference between FDM, FEM and FVM?

2. What is the best method and why?

3. Advantage and disadvantage of them?

Solution 1:

This is a difficult question to answer.

"The FDM is the oldest and is based upon the application of a local Taylor expansion to approximate the differential equations. The FDM uses a topologically square network of lines to construct the discretization of the PDE. This is a potential bottleneck of the method when handling complex geometries in multiple dimensions. This issue motivated the use of an integral form of the PDEs and subsequently the development of the finite element and finite volume techniques." (http://www2.imperial.ac.uk/ssherw/spectralhp/papers/HandBook.pdf)

Here are two references to review so you can get a better feel for these methods.

• http://files.campus.edublogs.org/blog.nus.edu.sg/dist/4/1978/files/2012/01/CN4118R_Final_Report_U080118W_OliverYeo-1r6dfjw.pdf (see page 10 for a very nice comparison in the types of problems they were interested in - computational fluid dynamics)

• There are some nice references for these methods at http://www2.imperial.ac.uk/ssherw/spectralhp/papers/HandBook.pdf (See section 7 for very nice references)

Solution 2:

Here is an old scicomp.SE question that answered some of your question: What are criteria to choose between finite-differences and finite-elements?

In my humble opinion, FEM is the most flexible one in terms of dealing with complex geometry and complicated boundary conditions. FEM also allows the adaptive/local procedure to get higher order local approximation or battling singularities. FEM's basis can be discontinuous and not well-defined pointwisely, which is a nice heritage from the Hilbert space framework. For computational fluid dynamics and electromagnetism, FEM is the way to incorporate the intrinsic geometrical properties of the solutions.

For FVM: partly you can refer to my answer here: How should a numerical solver treat conserved quantities? It is also worth noting that FVM can only have lower order of approximation.

In some recently development in FEM addresses the problem I mentioned in the answer above. For example, for convection-dominated pde, tradition continuous Galerkin framework for FEM doesn't work well, which introduces dissapation over time and oscillation over material-layers for the numerical solution. Now there are Discontinuous Galerkin FEM (higher order FVM) and hybrized DGFEM (see here: Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems) to remedy these two effects.

FDM and FVM are easy to implement, but you get trade-off from this convenience of implementation for limited usage for different PDEs.

Solution 3:

It shall be argued in this post that the whole idea of one Numerical Method being superior to another is merely a prejudice that rests upon insufficient in-depth analysis of the real thing.

The argumentation will proceed at hand of a two-dimensional example.
The reader is invited not to skip through but take notice of the details.

Numerical Analysis of Diffusion starts with a well known Partial Differential Equation (PDE). The problem will be restricted here to the simpler case of two space dimensions: $$ \frac{\partial Q_x}{\partial x} + \frac{\partial Q_y}{\partial y} = 0 $$ $ (x,y) = $ Cartesian coordinates. A possible interpretation of the vector $ (Q_x,Q_y) $ is the heat flux. The differential equation then follows from the law of conservation of energy. In case of pure diffusion of heat, also known as conduction, the components of the heat flux are related to temperature $T$ as follows: $$ Q_x = - \lambda \frac{\partial T}{\partial x} \qquad \qquad Q_y = - \lambda \frac{\partial T}{\partial y} $$ Where $ \lambda = $ thermal conductivity. Hence the final differential equation for the temperature field is actually of the second degree. In order to make the PDE amenable for numerical treatment, an integration procedure has to be resorted to. At this point, there occurs a splitting into several distinct roads, all leading to a numerical solution, more or less efficiently.

Triangle isoparametrics

The simplest Finite Element in two dimensions - and my absolute favorite - is the linear triangle:

Let's summarize the isoparametrics (= affine transformation) in the first place: $$ \begin{cases} x = x_1 + (x_2-x_1)\xi + (x_3-x_1)\eta \\ y = y_1 + (y_2-y_1)\xi + (y_3-y_1)\eta \\ f = f_1 + (f_2-f_1)\xi + (f_3-f_1)\eta \end{cases} $$ It will be shown next how partial differentiation at such a linear triangle takes place. First do the chain rules with global $(x,y)$ and local $(\xi,\eta)$ coordinates: $$ \begin{cases} \Large \frac{\partial f}{\partial \xi} = \frac{\partial f}{\partial x}\frac{\partial x}{\partial \xi} + \frac{\partial f}{\partial y}\frac{\partial y}{\partial \xi} \\ \Large \frac{\partial f}{\partial \eta} = \frac{\partial f}{\partial x}\frac{\partial x}{\partial \eta} + \frac{\partial f}{\partial y}\frac{\partial y}{\partial \eta} \end{cases} \quad \Longleftrightarrow \quad \Large \begin{bmatrix} \frac{\partial f}{\partial \xi} \\ \frac{\partial f}{\partial \eta} \end{bmatrix} = \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{bmatrix} \begin{bmatrix} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \end{bmatrix} $$ But what we need is the inverse, with determinant / Jacobian $\;\Delta = (\partial x / \partial \xi)(\partial y / \partial \eta) - (\partial x / \partial \eta)(\partial y / \partial \xi)$ : $$ \Large \begin{bmatrix} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \end{bmatrix} = \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{bmatrix}^{-1} \begin{bmatrix} \frac{\partial f}{\partial \xi} \\ \frac{\partial f}{\partial \eta} \end{bmatrix} = \begin{bmatrix} \frac{\partial y}{\partial \eta} & -\frac{\partial y}{\partial \xi} \\ -\frac{\partial x}{\partial \eta} & \frac{\partial x}{\partial \xi} \end{bmatrix} / \Delta \begin{bmatrix} \frac{\partial f}{\partial \xi} \\ \frac{\partial f}{\partial \eta} \end{bmatrix} $$ Giving full discretization of the function derivatives, with determinant / Jacobian $\;\Delta = (x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1)$ : $$ \Delta \begin{bmatrix} \partial f / \partial x \\ \partial f / \partial y \end{bmatrix} = \begin{bmatrix} (y_3-y_1) & -(y_2-y_1) \\ -(x_3-x_1) & (x_2-x_1) \end{bmatrix} \begin{bmatrix} f_2-f_1 \\ f_3-f_1 \end{bmatrix} $$ Here $\Delta$ is the area of a vector parallelogram, which is twice the area of the triangle. The above can also be written as: $$ \Delta \left[ \begin{array}{c} \partial f / \partial x \\ \partial f / \partial y \end{array} \right] = \left[ \begin{array}{ccc} +(y_2 - y_3) & +(y_3 - y_1) & +(y_1 - y_2) \\ -(x_2 - x_3) & -(x_3 - x_1) & -(x_1 - x_2) \end{array} \right] \left[ \begin{array}{c} f_1 \\ f_2 \\ f_3 \end{array} \right] $$ The matrix in this last formula should be memorized; it is called a differentiation matrix.

Finite Element Method

When using a Finite Element method, the differential equation may be multiplied at first with an arbitrary (test)function. Subsequently the PDE is integrated over the domain of interest. Let the test function be called $f$, then: $$ \iint f . \left[ \frac{\partial Q_x}{\partial x} + \frac{\partial Q_y}{\partial y} \right] \, dx dy = 0 $$ It can be shown that this integral formulation is (more or less) equivalent with the original partial differential equation. This is due to the fact that $f$ is an arbitrary function. It should be non-zero, continuous and integrable, though.
Partial integration, or applying Green's theorem (which is the same), results in an expression with line-integrals over the boundaries and an area integral over the bulk field. The latter is given by: $$ - \iint \left[ \frac{\partial f}{\partial x}.Q_x + \frac{\partial f}{\partial y}.Q_y \right] \, dx dy $$ Mind the minus sign. The advantage accomplished herewith is a reduction of the difficulty of the problem: only derivatives of the first degree are left. As a next step, the domain of interest is split up into "elements" $E$. Due to this, also the integral will split up into separate contributions, each contribution corresponding with an element: $$ - \sum_E \iint \left[ \frac{\partial f}{\partial x}.Q_x + \frac{\partial f}{\partial y}.Q_y \right] \, dx dy $$ It is clear that $\partial f / \partial x$ and $\partial f / \partial y$ are constants. While considering only 2-D diffusion, $Q_x$ and $Q_y$ are also partial derivatives of the first degree, hence constants. Herewith the bulk Finite Element formulation, for one triangle, is given by: $$ - \left[ \frac{\partial f}{\partial x}.Q_x + \frac{\partial f}{\partial y}.Q_y \right] \iint dx dy = - \left[ \frac{\partial f}{\partial x}.Q_x + \frac{\partial f}{\partial y}.Q_y \right] \Delta/2 $$ The remaining integral is equal, namely, to de area of the triangle. Applying now the differentiation matrix, we find: $$ = - \frac{1}{2} \left[ \begin{array}{ccc} f_1 & f_2 & f_3 \end{array} \right] \left[ \begin{array}{cc} y_2 - y_3 & x_3 - x_2 \\ y_3 - y_1 & x_1 - x_3 \\ y_1 - y_2 & x_2 - x_1 \end{array} \right] \left[ \begin{array}{c} Q_x \\ Q_y \end{array} \right] = $$ $$ = \frac{1}{2} \left[ \begin{array}{ccc} f_1 & f_2 & f_3 \end{array} \right] \left[ \begin{array}{c} (y_3 - y_2) Q_x - (x_3 - x_2) Q_y \\ (y_1 - y_3) Q_x - (x_1 - x_3) Q_y \\ (y_2 - y_1) Q_x - (x_2 - x_1) Q_y \end{array} \right] $$ Actually, we don't want to subdivide the Finite Element domain into triangular elements, but rather into quadrilateral elements. However, any quad element, in turn, can be subdivided yet into triangles, even in two different ways:

In addition, what we want is a configuration in which all quad vertices play an equally important role. In order to accomplish this, all of the four triangles must be present in our formulation, simultaneously. For just one quadrilateral, it boils down to renumbering vertices in the formulation for a single triangle, according to the following permutations:

   1  2  3       2  4  1       3  1  4       4  3  2

Also an upper label (not a power) will be attached to the values $(Q_x,Q_y)$, because it must be denoted at which triangle the discretization takes place. Any contributions are summed now over the four triangles (and the whole is divided by a factor two again): $$ \frac{1}{4} \left[ \begin{array}{ccc} f_1 & f_2 & f_3 \end{array} \right] \left[ \begin{array}{c} (y_3 - y_2) Q^1_x - (x_3 - x_2) Q^1_y \\ (y_1 - y_3) Q^1_x - (x_1 - x_3) Q^1_y \\ (y_2 - y_1) Q^1_x - (x_2 - x_1) Q^1_y \end{array} \right] + $$ $$ \frac{1}{4} \left[ \begin{array}{ccc} f_2 & f_4 & f_1 \end{array} \right] \left[ \begin{array}{c} (y_1 - y_4) Q^2_x - (x_1 - x_4) Q^2_y \\ (y_2 - y_1) Q^2_x - (x_2 - x_1) Q^2_y \\ (y_4 - y_2) Q^2_x - (x_4 - x_2) Q^2_y \end{array} \right] + $$ $$ \frac{1}{4} \left[ \begin{array}{ccc} f_3 & f_1 & f_4 \end{array} \right] \left[ \begin{array}{c} (y_4 - y_1) Q^3_x - (x_4 - x_1) Q^3_y \\ (y_3 - y_4) Q^3_x - (x_3 - x_4) Q^3_y \\ (y_1 - y_3) Q^3_x - (x_1 - x_3) Q^3_y \end{array} \right] + $$ $$ \frac{1}{4} \left[ \begin{array}{ccc} f_4 & f_3 & f_2 \end{array} \right] \left[ \begin{array}{c} (y_2 - y_3) Q^4_x - (x_2 - x_3) Q^4_y \\ (y_4 - y_2) Q^4_x - (x_4 - x_2) Q^4_y \\ (y_3 - y_4) Q^4_x - (x_3 - x_4) Q^4_y \end{array} \right] \mbox{ } $$ The four overlapping triangles at the corners of the quadrilateral are a small but significant twist to the standard Finite Element Method, which is motivated by the end-result.

Finite Volume Method

In order to save unnecessary paperwork, the following shorthand notation has been adopted. It may be interpreted as an outer product: $$ r_{ij} \times Q_k = (y_i - y_j) Q^k_x - (x_i - x_j) Q^k_y = (x_j - x_i) Q^k_y - (y_j - y_i) Q^k_x $$ Terms belonging to $f_k, k=1 ... 4$ are collected together. By doing so, the standard Finite Element assembly procedure is demonstrated at a small scale. What else is the Finite Element matrix than just an incomplete system of equations? $$ \frac{1}{4} \left[ \begin{array}{cccc} f_1 & f_2 & f_3 & f_4 \end{array} \right] \left[ \begin{array}{c} r_{32} \times Q_1 + r_{42} \times Q_2 + r_{34} \times Q_3 + 0 \\ r_{13} \times Q_1 + r_{14} \times Q_2 + 0 + r_{34} \times Q_4 \\ r_{21} \times Q_1 + 0 + r_{41} \times Q_3 + r_{42} \times Q_4 \\ 0 + r_{21} \times Q_2 + r_{13} \times Q_3 + r_{23} \times Q_4 \end{array} \right] $$ Subsequently use: $$ r_{32} = r_{34} + r_{42} \qquad r_{14} = r_{13} + r_{34} \qquad r_{41} = r_{42} + r_{21} \qquad r_{23} = r_{21} + r_{13} $$ To put the above in a more handsome form: $$ \left[ \begin{array}{cccc} f_1 & f_2 & f_3 & f_4 \end{array} \right] \left[ \begin{array}{c} \frac{1}{2} r_{42} \times \frac{1}{2} (Q_1+Q_2) + \frac{1}{2} r_{34} \times \frac{1}{2} (Q_1+Q_3) \\ \frac{1}{2} r_{13} \times \frac{1}{2} (Q_1+Q_2) + \frac{1}{2} r_{34} \times \frac{1}{2} (Q_2+Q_4) \\ \frac{1}{2} r_{21} \times \frac{1}{2} (Q_1+Q_3) + \frac{1}{2} r_{42} \times \frac{1}{2} (Q_3+Q_4) \\ \frac{1}{2} r_{21} \times \frac{1}{2} (Q_2+Q_4) + \frac{1}{2} r_{13} \times \frac{1}{2} (Q_3+Q_4) \end{array} \right] $$ It's a triviality, but nevertheless: a picture says more that a thousand words.


It is seen that the four pieces-of-equations correspond with four pieces of line-integrals, each of them belonging to one of the vertices. Midpoints of triangle sides are connected by lines at which the integration takes place. The heat flux at a midpoint is the average of values at the vertices.
Let's adopt another point of view now and no longer concentrate on elements but on vertices. Instead of arranging vertices around an element, elements are arranged around a vertex. Label triangle side midpoints as $a,b,c,d,e,f,g,h$.
It is immediately noted that the lines connecting the midsides of the triangles around a vertex, when tied together, neatly delineate a closed area, which can be interpreted as a kind of 2-D Finite Volume. Expressed in the outer product formalism, we find: $$ r_{ba} \times Q_a + r_{cb} \times Q_c + r_{dc} \times Q_c + r_{ed} \times Q_e + r_{fe} \times Q_e + r_{gf} \times Q_g + r_{hg} \times Q_g + r_{ah} \times Q_a $$ Which is the content of one equation in the Finite Element global matrix. All terms together represent a discretization of the following circular integral: $$ \sum r \times Q = \oint Q_y dx - Q_x dy $$ With help of Green's theorem, however, such a circular integral can be converted into a "volume" integral, over the area indicated in the above figure: $$ \oint Q_y dx - Q_x dy = + \iint \left[ \frac{\partial Q_x}{\partial x} + \frac{\partial Q_y}{\partial y} \right] \, dx dy $$ Conservation of heat is integrated over a finite volume, which is wrapped around a vertex. So we have arrived at sort of a Finite Difference method. To be more precise: at a Finite Volume Method. It is remarked that this F.V. procedure has been applicable for curvilinear grids from the start.
Apply a Finite Element (Galerkin) method to a mesh of quadrilaterals. Subdivide each of the quads into four (overlapping) triangles, in the two ways that are possible. Then such a method is equivalent to a Finite Volume method: midsides of the triangles, around the vertex of interest, are neatly connected together, to form the boundary of a 2-D finite volume, and the conservation law is integrated over this volume.
Unification of a Finite Element and a Finite Volume method has been accomplished herewith, for a restricted class of 2-D diffusion problems.