Solving PDEs numerically: the finite element method (1)

The study of partial differential equations is a fascinating field of mathematics with many concrete applications, be it in physics, mechanics, meteorology, medicine, urbanism… Mathematicians model the world as equations to better understand it and, if possible, compute a theoretical solution to a physical problem. However, we can’t always get an analytical solution – for example if the geometry is too complex – and sometimes must rely on numerical methods to get an approximation.


Partial differential equations can model a lot of phenomena we see in our daily lives: wave propagation (such as sound propagation or light diffusion), fluid dynamics, mechanical deformation, etc.. They are a powerful tool to analyze and predict the behavior of physical systems both in space and time.

The finite element method is a numerical method that helps us solve a lot of engineering problems like heat transfer, fluid flow… It is complex and really interesting, so I plan on devoting a few articles to it. I won’t go into the mathematical details too much; instead, I’d rather present the core ideas of the method and try to sum up its main goals. But I will try and give you as many references as possible if you’re interested in reading more – and these articles may still be more technical than the last ones!

This first article explains the basic notions of the finite element method. The second one will focus on some implementation details. Finally, the third one will talk about advanced techniques and further developments.

Modeling of the diffusion of a Wi-Fi wave in a very simple apartment with the source at the center, thanks to the finite element method

The idea for this series of articles comes from the course I took at my school of Engineering this year, taught by Bertrand Thierry, that really fascinated me! I dived more into the method during the final project and took it as a chance to start to develop a small library in Python, called SimpleFEMPy, to solve PDEs with the finite element method. It is far from being complete and not yet available, but it already shows some promise and I hope to put it online for free in the months to come.

It takes a lot after FreeFem++, an open-source finite element software developed by researchers from the french university Paris Sorbonne. It is very efficient and compatible with many usual tools in the domain: mesh generators like GMSH, scientific visualization software like Paraview

So, what is the finite element method?

The main idea of the finite element method is to find an approximation of the solution to a static (i.e.: time-independent) physical problem. To do so, we discretize our domain spatially in order to get an equivalent of our world with a finite number of points so we are able to ask a computer to work on the problem with us (I already talked about discretization in a previous article, don’t hesitate to check it out!). We then model our problem as a linear system and solve it to compute the solution at each point of our discretized domain.

This part is a very crude overview of the basic ideas of the method.

Spatial discretization

So, basically, we consider that the part of the world our phenomenon lives in can be – approximately – ‘cut down’ into small construction blocks (the ‘elements’) that, when pieced together, more or less reform the whole portion of space. If you are in 2D, the elements are flat polygons: triangles or squares, usually. If you are in 3D, they become tetrahedrons or cubes.

Let’s say we take the ‘unit square’ or the ‘unit cube’, meaning the 1×1 square or 1x1x1 cube centered around the origin of our spatial system. Here is a little embedded visualization tool to better understand how these reference spaces can be cut into ‘elements’:

Warning: This app is not completely responsive, so to get the best experience, you should take a look at it on a computer…

[code_fem_spacedisc]

The 2D/3D buttons set the number of dimensions in your reference space (2 for the unit square, 3 for the unit cube). The Regular/Triangulation buttons set the type of discretization.

If you want to better approximate complex geometries (in particular, regular domains like circles), you must refine your elements and choose new ones like curved triangles… but this greatly complicates mathematical expressions!

Modeling the problem

Once you have defined your elements, you need to ‘put your actual physical problem into words’. Or, to be precise, into mathematical equations. Typically, a physical model solved by the finite element method is defined by two things:

  • one or more equations that represent(s) the behavior of the quantity we are studying: the solution we are searching for (so, the physical quantity we are interested in) can be linked to its derivatives and to the physical characteristics of the problem thanks to one or more mathematical equations

For example, if you study heat diffusion, then you look at the heat source in your room and the space derivatives of the temperature.

  • the boundary conditions: these define a specific behavior of the solution or its derivative at the limit of the domain (usually called its ‘border’ or its ‘boundary’)

Indeed, because we study static problems, we don’t have an initial condition to insure a unique solution. Instead, we say that we have a boundary value problem and the specific values of our solution (or its derivatives) at the border of our domain are the only way to – maybe – mathematically prove we have only one possible solution.

‘Weakening’ the equations

Even if I won’t go into the details here, because it would require quite a bit of mathematical background, it is important to note that the finite element method doesn’t solve the actual physical problem per se. In fact, due to some complex math stuff (I recommend taking a look at Sobolev spaces if you want to learn more), searching for a solution that satisfies our aforementioned equations exactly can be too hard. Thus, we ‘weaken’ these equations and replace them with a ‘weak’, or ‘variational formulation‘.

Our goal isn’t to find the exact solution to the exact physical problem but the exact solution to an approximate problem that is globally equivalent. This is a way to deal with situations where the solution isn’t ‘smooth’ enough to be computed directly.

Checking how far off we are from the real solution is therefore necessary, and this is the whole point of finite element analysis, which I’ll talk a bit about in the third article of the series.

Solving the equations

Now that we have both our discrete domain and our equations, we need to actually solve the problem. To do so, the first step is to convert our equations to a linear system; in other words, we state that our solution x verifies a formula like: A.x = b, where:

  • A is a square matrix with as many rows as there are points in the discrete domain; its coefficients are extracted from the parts of our equation(s) that depend on the solution or its derivatives
  • b is a column vector with the same number of rows; its coefficients come from the parts of our equation(s) that do not depend on the solution or its derivatives but rather on purely physical data (a heat source, an incident acoustic wave…)

Given this system, we can the second step is now to solve it for each point of our discrete domain. So, in the end, we have an approximate solution of our ‘weakened problem’, which is itself an approximate problem derived from the initial model of the physical situation.

Nowadays, a lot of very powerful libraries are dedicated to solving this kind of system efficiently. For example, in Python, you can use NumPy and SciPy.

One of the reasons these libraries are so quick is because they rely on the most recent BLAS implementations, a set of routines designed to optimize linear algebra computations.

Usually, the A matrix is actually decomposed into two matrices: the ‘mass matrix’ M and the ‘stiffness matrix’ K. The first one considers the terms that depend on the solution itself and the second one considers the terms that depend on the solution’s first derivatives. The names ‘mass’ and ‘stiffness’ were coined because of the applications of the finite element method in the field of solid mechanics.

Possible applications

As I’ve said before, the finite element method is a very powerful tool that can be applied to many real-life examples. In particular, it can be combined with other numerical methods like the finite difference method to also model and solve time-dependent problems.

I’m showing a few (simple) examples of FEM applications below!

Acoustic diffraction

Here is a visualization of the diffraction of sound on the middle of a ring depending on the angle of the incident wave: imagine you throw a sound wave from outside the ring and study how it collides with the inner border. With this simple shape, it is quite obvious what will happen. But if you replace the ring with a more complex geometry like a plane (in 2D top-view), you can get an intuition of how well it is designed for aerodynamics.

Diffraction of sound on the middle of a ring depending on the angle of the incident wave

Same phenomenon on a 2D plane geometry: we see how the shape of the plane is carefully designed to diffract waves

Wave interference

When two waves produced by identical sources meet, a phenomenon called wave interference occurs. Based on the position of the sources, when these waves meet, they can either combine to make an even bigger wave (this is a constructive interference) or cancel out (this is a destructive interference).

Schematic of the wave interference concept (by Haade, from Wikipedia: https://commons.wikimedia.org/wiki/File:Interference_of_two_waves.svg)

This small animation shows this phenomenon: consider you have a fixed source on the left and one that moves around on the right – then, sometimes, waves meet up and combine, or sometimes they destroy each other.

Wave interference between two identical sources depending on the position of the second source

Wave diffusion and absorption

Suppose you have a very simple – and unrealistic – apartment composed of 3 rooms. You want to know where to put your Wi-Fi router to optimize the path of the waves – for example, have it spread to the kitchen but avoid your baby’s bedroom. Physics provide you with a set of equations that precisely model the behavior of the Wi-Fi wave emitted from your source and take into account the difference between the ‘open spaces’ and the walls.

With this very basic model, we already observe the diffraction phenomenon when the wave passes through the door on the right and the almost complete absorption when it tries to traverse a (plaster) wall instead of air.

Wi-Fi wave diffusion in a small apartment (with plaster walls) depending on the position of the router

In the following chapter…

In the second article, I will talk more about the actual implementation of the thing. By presenting some objects of my new Python library, SimpleFEMPy, I will try and explain how these mathematical concepts are represented and coded on a computer.

The reference section at the end of this article gives you a few links to great courses about the finite element method. Most are in English and greatly helped me understand the secrets of the FEM…

References
    1. My teacher’s website (the course notes are in French, but really check it out if you’re interested in the topic!): https://www.ljll.math.upmc.fr/bthierry/
    2. Joseph Flaherty’s course on Finite Element Analysis (see ‘Course Notes’): http://www.cs.rpi.edu/~flaherje/
    3. Pascal Frey’s course on Finite Element Approximation: https://www.ljll.math.upmc.fr/frey/cours/UdC/ma691/ma691_ch7.pdf
    4. Dr. Olivier De Weck and Dr. Il Yong Kim’s course on FEM for solid mechanics and CAD/CAE: http://web.mit.edu/16.810/www/16.810_L4_CAE.pdf
    5. Aude Rondepierre and Adeline Rouchon’s course on PDE theory (in French): https://www.math.univ-toulouse.fr/~rondep/CoursTD/polyic2.pdf
    6. FreeFem++’s website: https://freefem.org/
    7. GMSH’s website: http://www.gmsh.info/
    8. Paraview’s website: https://www.paraview.org/
    9. I. Wikimedia Foundation, “Finite element method” (https://en.wikipedia.org/wiki/Finite_element_method), December 2018. [Online; last access 09-January-2019].
    10. I. Wikimedia Foundation, “Boundary value problem” (https://en.wikipedia.org/wiki/Boundary_value_problem), February 2018. [Online; last access 13-January-2019].
    11. I. Wikimedia Foundation, “Weak formulation” (https://en.wikipedia.org/wiki/Weak_formulation), August 2018. [Online; last access 13-January-2019].
    12. I. Wikimedia Foundation, “Sobolev space” (https://en.wikipedia.org/wiki/Sobolev_space), December 2018. [Online; last access 13-January-2019].
    13. I. Wikimedia Foundation, “Finite difference method” (https://en.wikipedia.org/wiki/Finite_difference_method), December 2018. [Online; last access 14-January-2019].
    14. I. Wikimedia Foundation, “Basic Linear Algebra Subprograms” (https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms), November 2018. [Online; last access 13-January-2019].
    15. I. Wikimedia Foundation, “Wave interference” (https://en.wikipedia.org/wiki/Wave_interference), January 2019. [Online; last access 14-January-2019].

Leave a Reply

Your email address will not be published. Required fields are marked *