Solving PDEs numerically: the finite element method (3)

In this last article from the FEM series, I present some additional concepts related to the finite element method and further developments. In particular, I focus on finite element analysis, some methods for solving time-dependent problems and optimization ideas.

The goal here is not to go into details about all of these topics – because they would each require at least one entire article! – but just to present them roughly and spike your interest.

Finite Element Analysis

As we’ve seen before, the finite element method relies on approximations: we discretize our domain, we use variational formulations instead of the initial physical equations and we then turn these formulations into linear systems; in the end, we have an approximated solution of an approximated problem!

What we might wonder about at this point is: how truly viable is this technique? In other words, how far are we from the actual solution?

This is the focus of Finite Element Analysis (FEA): this subfield of mathematics is devoted to applying the FEM to solve physical problems and estimating the approximation errors we make when computing our solution. ‘FEA’ refers to the whole process of meshing the domain, solving the problem and analyzing the results.

For example, FEA studies the ‘mesh convergence’ of the model, i.e. how accurate the solution is depending on how refined the discretization of our domain is. Of course, the more points you have in your discretized domain, the more precise you will be. However, this implies more computation and a longer solving time. Thus, we need to make a trade-off between the meshing and solving steps. Sometimes, it is more interesting to only refine part of our domains: by using these so-called ‘adaptative meshes’, we can keep a coarse discretization where we know accuracy can’t be improved but have a fine one where necessary. Therefore, we increase the precision of our solution without increasing too much the size of the problem itself.

Measuring convergence (which roughly corresponds to accuracy) is a purely mathematical problem I won’t dive in here – it uses various norms and analysis tool; if you are interested, I encourage you to read this article: “What is Convergence in Finite Element Analysis?”.

Applying the FEM for time-dependent problems

In my previous articles, I explained that the finite element method only applies to static problems, meaning problems that do not depend on time. On the other hand, you may be aware that most of the famous physical partial differential equations we want to solve involve time. So, does this mean that the FEM is restricted to a few problems nobody ever deals with?

Not really! The trick is to think of time-dependent problems as a series of time-independent problems… at each timestep, you can think of your current state as a new subproblem to solve, completely independent of time, and therefore you can apply the finite element method on these new subproblems.

In a lot of examples, the domain does not change over time. This means that we can mix the FEM with a finite difference method, for example, to solve our time-dependent problem:

  • first, we can discretize the domain once and for all at the beginning
  • now, we can apply a time-discretization scheme like a Euler scheme or the Crank-Nicolson scheme: these gradually compute the evolution of our solution at each timestep based on an initial solution and the result at the previous iteration
  • finally, we can use the finite element method at each iteration to get the current solution

This technique is quite useful because it is not too hard to implement and it allows us to compute quite accurate results for time-dependent physical problems.

Still, you have to solve a problem with the FEM at each iteration. If the domain is not modified over time, it basically reduces to multiplying matrices and vectors or, at most, solving a linear system (depending on the time-discretization scheme you chose). Both these can be done efficiently nowadays thanks to BLAS routines. But if the domain changes, the mix technique asks for a lot of computation because you need to recompute the matrices given by your discretization whenever it does…

To sum up, solving time-dependent problems with the finite element method requires you to select:

  • a time-discretization scheme: Euler explicit, Euler implicit, Crank-Nicolson… and the associated time step
  • a discretization of your domain for the FEM part (with, among other things, a refinement step)

It is worth noting that time-discretization and space-discretization steps must verify a few criteria for the process to work: it may seem counter-intuitive but choosing a refinement that is too precise may endanger the stability of the time-discretization scheme!

Implementation optimization tricks

As I was programming my small FEM library, SimpleFEMPy, I came across a few interesting optimization tricks. Even if I didn’t implement all of them, I thought I’d mention them here nonetheless.

What if your problem changes just a little?

Let’s assume you have taken an equation and turned it into a variational formulation. You should get a left-hand side with integrals that depend on the solution and a right-hand side with others that do not. Then, we know that the finite element method eventually relies on the solving of a linear system of the form: A.x = b (where A is a square matrix that represents the left-hand side of our weak formulation, x is a column vector that holds our solution and b is a column vector that represents the right-hand side of our weak formulation).

Now, let’s say you want to solve the same equation for 2 different source terms. Since the source term doesn’t depend on the solution, this means that only the right-hand side will change from one case to the other.

This may not seem like much at first, but it is actually really interesting. Indeed, under the hood, ‘solving the system’ usually means ‘factorizing A in a nice form’ first and then solving a simpler system. But, here, look: A does not change when we take another source term. So, we only need to factorize it once, and then we can reuse this ‘nice form’ the second time to bypass the first step and only solve the second simple system. And, of course, this is still true if we take a third problem configuration where only solution-independent variables change… So, as long as you don’t modify your left-hand side, you can get your new solution very quickly if you have stored your factorized A matrix.

As usual, this optimization trick involves a trade-off between storage and computation time… Be aware that if you are on a memory-limited device and if your domain has many points, given the size of A, this may not be the best idea…

Sparsity is the key!

Last time, we talked about the assembly of the A matrix. In particular, we saw that it is better to iterate over our triangles than over our points because there are fewer so the loop will be shorter, and that we just need to hold a table somewhere to match local and global indices. Basically, for each triangle in our mesh, we simply compute the values of our matrix on its 3 vertices, and in the end we sum all these ‘elementary contributions’ into our global matrix.

It is easy to see that, with that process, we will get null values more often than not: as soon as a vertex is not in a triangle (so, almost all the time since we only consider 3 points out of every other for one triangle), there will be no contribution.

This means that our A matrix is filled with a lot of zeros: we say that it is ‘sparse‘. To be more precise, a sparse matrix is a matrix where more than half of the coefficients are zeros; so, here, it is a very sparse matrix. Just to give you some idea of how empty the left-hand side matrix can be, for a medium-sized domain and a basic Helmholtz equation like the one I showed in the second article, there are only about 12-15% of non-null coefficients in A!

The upside is that sparse matrices can be stored efficiently in memory to limit the required memory size: in fact, there is no need to store all the zeros; instead, it is more efficient to remember the position and the value of the non-null coefficients, since we can assume all the others are null. This is the idea of the COO (‘coordinate list’) format, for example, where the matrix is stored as 3 lists: the rows, the columns and the values of the non-null coefficients. The LIL (‘list of lists’) format is quite similar but has one list of column index and value for each rows of the matrix. Both these formats are interesting for incremental construction since you just need to append the new coefficients without having to worry about redundancy.

However, if you want to get rid of the duplicate positions afterwards, CSR or CSC formats (for ‘compressed sparse row’ and ‘compressed sparse column’) are a good way to go: with these formats, you store how many values there are in the current row or column, and then only the column or row index and value of the coefficient. These formats are optimized for matrix-vector multiplication and for row or column access. Be careful: if you want to change the sparsity of a matrix in the CSR or CSC format, you’ll need to modify all the compressed indices!

This small example shows you how a small matrix would be stored in the various formats:

COO, LIL, CSR and CSC formats for a small example matrix

Here, we see that the various formats require more or less data; still, they all store less than 16 values as a dense matrix would. The CSR and CSC formats are really interesting as soon as the (m rows x n columns) matrix you want to store contains less than (m (n-1) – 1) / 2 non-null values.

Using the right format at the right time in our process allows us to avoid storage issues and to perform some computations more efficiently.

Relying on tried and tested libraries

I’ve talked about BLAS routines several times in this series, but I’ll say it one last time: linear algebra is complex and some people have already worked on it a lot, which resulted in these optimized algorithms for classical linear algebra operations. Whenever you aim at efficiency when multiplying big matrices and vectors, factorizing matrices or solving linear systems – just what we need for the FEM -, you are probably better off finding a well-written library with optimized objects to make the best of BLAS rather than trying to write your own method.

In Python, NumPy and SciPy are well-known references. These libraries allow you to use multidimensional arrays (namely: matrices) and apply all the useful operations efficiently. They make use of the BLAS routines to solve linear systems quickly and provide easy-to-use sparse matrices.

In C++, many linear algebra libraries are available such as Eigen or Armadillo. Many of those libraries actually rely on the Fortran LAPACK library, but they provide a nice interface to this linear algebra package. For the really crazy, a few Javascript libraries or Github projects have been developed to do maths in a browser (mathjs, deeplearn.js…), but scientific computation online is hardly as efficient.

However, learning how these libraries actually work is worthwhile. To me, knowing the basics of high-performance computation is essential for engineers today. Even if you don’t immerse yourself into all the technical details, it is beneficial to grasp the main idea of asymptotic performance bounds, LU/QR decomposition, direct and iterative solvers…

Just be aware that you probably won’t be able to match the optimized running time of the aforementioned libraries, sadly!

Final conclusion

In this series, I’ve scratched the surface of the field of Finite Element Analysis. Although quite young (it really started in the 1960s), this branch of mathematics has already found many applications in the industry and everyday-life problems solving.

For example, the Laplace equation appears in many problems in physics, mechanics, fluid dynamics… When we take funky parameters, we can get fun-looking results like this one!

While the common finite difference method is a good way of computing PDEs numerically, the FEM allows us to work on more complex geometries and gives us a more accurate spatial discretization. It usually gives a better solution – even if it is extremely problem-dependent – and has already proven how powerful it can be on many engineering problems.

There are various other developments of the finite element method like the stochastic FEM to solve probability models numerically or problem-specific domain refinement (with adaptative meshing or with additional points in our triangles, like the P1-bubble polynomials in mechanics).

I’ve enjoyed my course on finite elements a lot and I am currently working on SimpleFEMPy but also on a free online FEM solver project with Bertrand Thierry: Hopefully, it will go live soon!

    1. My teacher’s website:
    2. Joseph Flaherty’s course on Finite Element Analysis (see ‘Course Notes’):
    3. I. Wikimedia Foundation, “Finite element method” (, December 2018. [Online; last access 09-January-2019].
    4. Simwiki, “What is FEA | Finite Element Analysis?” (
    5. Ajay Harish, “How to measure convergence in FEA” (, Nov. 2018
    6. I. Wikimedia Foundation, “Finite difference method” (, December 2018. [Online; last access 24-February-2019].
    7. I. Wikimedia Foundation, “Sparse matrix” (, February 2019. [Online; last access 24-February-2019].
    8. I. Wikimedia Foundation, “Basic Linear Algebra Subprograms” (, November 2018. [Online; last access 13-January-2019].
    9. NumPy’s website:
    10. SciPy’s website:
    11. Eigen’s website:
    12. Armadillo’s website:
    13. LAPACK’s website:
    14. mathjs’ website:
    15. deeplearn.js’ website:

Leave a Reply

Your email address will not be published.