domain: math info phys spi

text:
Franz Chouly
email: franz.chouly@u-bourgogne.fr

An introductive course to some numerical approximation methods for ordinary and partial differential equations

Introduction

These lectures notes were written for students of the first year of the international Master 'Maths for Physics' at the University of Burgundy in 2021, as a part of the course called 'Numerical Methods'. It was designed for a total of 22 hours of classes. The main objective of this course is to provide a first insight into numerical methods to solve mathematical problems inspired from physics and engineering. It focuses both on the mathematical fundations and justifications of numerical methods and on practical aspects related to their implementation. We start to study ordinary differential equations (ODEs), and present some methods for their discretization, i.e., to approximate the timederivative. Then we focus on partial differential equations (PDEs), and present a well-known method, the Finite Element Method (FEM), to approximate the solution of elliptic PDEs. The last part will be devoted to time-dependent PDEs, such as the heat equation, where we combine both ideas coming from the two previous chapters. I thank deeply all the students for their remarks associated to this course, and for their enthusiasm despite the very special context, and I thank also my colleagues at the IMB and ICB for their support.

Discretization of ODEs

This part concerns ordinary differential equations and their numerical approximation. Ordinary differential equations are omnipresent to describe the behavior of many systems in physics and applied science. For instance (see [START_REF] Hairer | Solving ordinary differential equations[END_REF]), if you are interested in computing the velocity v of a small ball of mass m when it is dropped from a given height, you need to solve:

v (t) = g - k m v(t), t > 0, v(0) = v 0 ,
where t is the time variable, g is the gravity constant, k is the resistance of the air, and v 0 is the initial value of the velocity. If you are interested in predicting the motion of a spherical pendulum, you need to solve: 0 = -sin(θ(t)) + sin(θ(t)) cos(θ(t))( φ(t)) 2 θ(t) 0 = -sin 2 (θ(t)) φ(t) -2 sin(θ(t)) cos(θ(t)) φ(t) θ(t), t > 0 ϕ(0) = ϕ 0 , φ(0) = φ0 , θ(0) = θ 0 , θ(0) = θ0 , where (θ, ϕ) represent the Euler angles of a pendulum of lenght 1 (with mass and gravity constant 1, also). Initial conditions are given, respectively the initial position (θ 0 , ϕ 0 ) and velocity ( θ0 , φ0 ) It is well known that a lot of differential equations modelling physical problems do not admit, up to what is known today, closed-form solutions. For instance the first ODE we presented may have one, but the second one surely not. So we need some techniques to approximate their solution efficiently and accurately.

Introduction and simplest schemes

Let us start to study this simple initial-value problem:

y (t) = f (t, y(t)), y(0) = y 0 . (1)
Here y : R + → R is the (unknown) solution, a function that depends on the (time) variable t, y is the first-order derivative of y, and f : R + × R → R a given (known) function of two variables. The (real) scalar y 0 is known. The first equation of ( 1) is a scalar first-order ordinary differential equation (ODE), and the second one is the initial condition, at time t = 0. We suppose that f is kind enough, for instance globally Lipschitz, so that (1) admits one unique solution y, well-defined and continuously differentiable on the whole interval R + . This setting is simplistic, of course, but will be sufficient to introduce almost all the ideas related to numerical approximation. Most of the techniques and proofs that we will see in the reminder extend easily to more general setting (ODE well-posed only within a bounded interval of R + , second-order ODE, systems of ODE, etc). Suppose that f , though kind enough to make (1) is well-posed, is not so kind, however, so that it is easy to find an explicit, closed-form solution. So if we are really interested in solving this, one realistic option is to solve it approximately using a computer. The main problem, at this stage, is that a computer is very efficient to perform quickly basic arithmetical operations, but not clever enough to understand the concepts of function and derivative. So we can not give it directly [START_REF] Allaire | An introduction to mathematical modelling and numerical simulation[END_REF] and tell him : 'Please, solve this !'. We need to preprocess it for him, so that he has to deal only with sequences of numbers. The simplest way to do this is to slice the half-real line into small slabs, of size τ , where τ > 0 is a chosen parameter, called the time-step. As a result we define the sequence of real values t n := nτ , for any n ∈ N: 0 = t 0 < t 1 (= τ ) < . . . < t n (= nτ ) < . . .

To each value of n, we also associate a real unknown y n , and we will work enough so that at the end, the computer provides an explicit value of y n such that y n y(t n ), and we will precise later on the meaning of . The idea is that, if τ is small enough, the sequence (y n ) is fine enough to represent approximately the function y. Of course, plenty of values of y are missing, but could be, if necessary, recovered, using, for instance, linear interpolation between two successive values y n and y n+1 . Now, we need to approximate the value of the derivative y using the sequence (y n ). The basic idea is to make use of the finite difference:

y (t n ) y n+1 -y n τ .
Let us just remark for the moment, that, if (y n ) represents accurately the function y, and if τ is very small, we should approximate correctly the derivative. The initial condition is simply transposed as

y 0 = y 0 .
From now on, three possible options arise naturally for the discretization of [START_REF] Allaire | An introduction to mathematical modelling and numerical simulation[END_REF], that depend each one of the treatment of the right-and side of the first equation (f (t, y(t))):

1. The first option: we take f (t, y(t)) at time-step n, which leads to the equation, for all n ≥ 0:

y n+1 -y n τ = f (t n , y n ). ( 2
)
This is called the forward Euler scheme.

2. The second option: we take f (t, y(t)) at time-step n + 1, which leads to the equation, for all n ≥ 0:

y n+1 -y n τ = f (t n+1 , y n+1 ). (3)
This is called the backward Euler scheme.

3. The third option: we take a barycentric combination of time-steps n and n + 1. This leads to the following equation, for all n ≥ 0:

y n+1 -y n τ = 1 2 (f (t n , y n ) + f (t n+1 , y n+1 )). ( 4
)
This is called the Crank-Nicolson scheme.

Before entering the mathematical analysis of these procedures, let us just give a first insight into their implementation. We illustrate the implementation aspects using Octave, which is an open-source clone of Matlab, that has a graphical user interface and a script language very convenient for scientific computing. It can be downloaded here: https://www.gnu.org/software/octave/index There is also an online version of Octave, available here: https://octave-online.net

We recommend to download the package on your computer, which makes its use more comfortable. Nevertheless, the above online version can be helpful in case of emergency for instance. Of course, other alternatives can be chosen if preferred, for instance the commercial software Matlab (if you are rich enough), the open-source french software Scilab (from the INRIA, https://www.scilab. org) or Python-based solutions such as SciPy, with nice interfaces such as Spyder or Jupyter notebooks (https://www.scipy.org/about.html).

Let us start with the implementation of the first scheme (2), forward Euler. Note that the relationship (2) can be rewritten:

y n+1 = y n + τ f (t n , y n ). (5)
This provides a simple recurrence formula to compute the approximation y n+1 , at time t n+1 , for the previous approximated solution y n at time t n . There is nothing more to do. Since the approximated solution at a given discrete time is inferred explicitely from the previous solution, this scheme is called explicit (and is sometimes referred as explicit Euler time-marching scheme). For the backward Euler scheme (3), it can be also rewritten as follows:

y n+1 -τ f (t n+1 , y n+1 ) = y n . (6)
This provides also a recurrence formula, but this is not as simple as for forward Euler. Indeed, the approximation y n+1 is not obtained explicitely, and, instead, needs to be found as the zero of the function I -τ f (t n+1 , •) -y n (I denotes the identity). This operation can be done by many ways, for instance using Newton's method, but at the price of an extra effort, and an extra computational expense. For this reason, this scheme is called implicit (and is often referred to as implicit Euler time-marching scheme). Despite of this difficulty, this scheme remains interesting in many situations. We will see the precise reason for this later on. The Crank-Nicolson scheme (4) can be rewritten as:

y n+1 - 1 2 τ f (t n+1 , y n+1 ) = y n + 1 2 τ f (t n , y n ). (7)
Remark that this scheme remains still implicit. Remark also that any other convex combination between the time-steps n and n+1 could be meaningful, but Crank-Nicolson refers to the values ( 1 2 , 1 2 ) (otherwise we talk about θ-schemes). At first glance, this scheme seems to have no advantage above the backward Euler scheme, but we will see later on it has some.

A final remark about implementation is that, though all these schemes are very simple to implement, their main drawback is that, to obtain the solution at a given final time, you need to compute iteratively all the intermediate solutions, which can be long for some problems. In fact, for ODEs, there is no straightforward way to solve them quickly using parallel architectures, and some special schemes are needed if you want to achieve this (see for instance [START_REF] Gander | 50 years of time parallel time integration[END_REF] or [START_REF]Time parallel time integration[END_REF]). Now, let us provide another viewpoint on the above schemes. Let us integrate the first equation of (1) between time t n and time t n+1 :

y(t n+1 ) -y(t n ) = tn+1 tn f (t, y(t))dt. (8)
In some sense, we transform our initial problem to a problem of another kind: the computation of an integral. Now let us try, instead of approximating the derivative y in (1), to approximate the integral in the right-hand side of [START_REF] Gander | Numerical analysis of partial differential equations using Maple and MATLAB[END_REF].

There are many quadrature rules that can be used, and let us try with the simplest ones. For instance, if we use the left rectangle rule

b a ϕ(t)dt (b -a)ϕ(a)
we get from (8):

y(t n+1 ) -y(t n ) τ f (t n , y(t n )) (9)
from which we recover the forward Euler scheme, by substituting y n to y(t n ) for each value of n. You can also show, as an exercice, that the backward Euler scheme, resp. the Crank-Nicolson scheme, is recovered identically from the right rectangle rule, resp. the trapezoidal rule. Now, to conclude this first part, let us give some very first notions about numerical stability of the above schemes. To this purpose, we consider the following variant of (1):

y (t) = -α y(t), y(0) = 1. ( 10
)
Here α > 0 is a real parameter. The above equation has the following explicit solution:

y(t) = exp(-αt).
This solution remains bounded, in the sense that its modulus remains smaller than 1:

|y(t)| ≤ 1, ∀t ∈ R + .
Let us approximate now this solution using the forward Euler scheme. Following equation ( 5), we obtain:

y n+1 = y n -ατ y n = (1 -ατ ) y n . (11)
Here, there is no need to use a computer to obtain the discrete solution, which is simply given by the formula below, for any value of n:

y n = (1 -ατ ) n .
So the modulus of the discrete solution is:

|y n | = |1 -ατ | n .
We expect, that, at least, the discrete solution remains bounded, but we observe here that this is not always the case. A necessary and sufficent condition, in this situation, is that:

|1 -ατ | ≤ 1
or, equivalently:

-1 ≤ 1 -ατ ≤ 1 which can be reformulated as τ ≤ 2 α ,
(where we used that α and τ are positive). So the discrete solution (y n ) remains bounded only if the time-step τ is small enough. We say in this case that the forward Euler scheme is conditionally stable: stable if τ is small enough, unstable otherwise. Of course, we want to avoid numerical unstability, since, in this situation, the discrete solution approximates very poorly the exact one. Using the same argument, one can assess that the backward Euler scheme and the Crank-Nicolson scheme are unconditionnaly stable (these are good exercices).

One-step methods

Now that we have a first insight into numerical approximation of ODEs, and before going into more complex (and widely used) methods, let us be more precise and go into the mathematical analysis of the simplest methods. To this purpose, let us introduce a general class of approximation methods, called one-step methods, as follows:

y n+1 = y n + τ Φ τ (t n , y n ), (12)
for n ∈ N, and where Φ τ : R + × R → R is a given function, that may depend on the time-step τ , and that is supposed at least continuous. The simplest example of one-step method is the forward Euler scheme [START_REF] Ern | Theory and practice of finite elements[END_REF], and this case corresponds to Φ τ (t, y) = f (t, y).

We have seen previously that y n is only an approximation of y(t n ), and we will try now to say something general and precise about this approximation. For this purpose, we introduce the global error :

e n := y(t n ) -y n at each step n. We introduce also the truncation error T n as

T n := y(t n+1 ) -y(t n ) τ -Φ τ (t n , y(t n )),
where y is the solution to Problem [START_REF] Allaire | An introduction to mathematical modelling and numerical simulation[END_REF]. The truncation error is the basic notion, and represents roughly how much the one-step method approximates the ODE.

The first interesting point is that we can infer a bound on the global error from the truncation error (see [START_REF] Süli | An introduction to numerical analysis[END_REF]):

Theorem 1 Suppose that, for the one-step method (12), the function Φ τ is Lipschitz: there exists L > 0 such that, for all t ∈ R + and u, v ∈ R, there holds:

|Φ τ (t, u) -Φ τ (t, v)| ≤ L|u -v|, (13)
then there holds, for all n ∈ N:

|e n | ≤ 1 L e Ltn -1 max k=0,...,n |T k |. (14)
Proof. First let us rewrite differently the definition of the truncation error:

y(t n+1 ) = y(t n ) + τ Φ τ (t n , y(t n )) + τ T n
and then we substract [START_REF] Hairer | Solving ordinary differential equations[END_REF]:

y(t n+1 ) -y n+1 = y(t n ) -y n + τ Φ τ (t n , y(t n )) -τ Φ τ (t n , y n ) + τ T n
so we obtain an equation on the global error:

e n+1 = e n + τ (Φ τ (t n , y(t n )) -Φ τ (t n , y n )) + τ T n .
Now we use the triangle inequality and the Lipschitz condition [START_REF] Hairer | Solving ordinary differential equations[END_REF] to bound the left hand-side above:

|e n+1 | ≤ |e n | + τ L |y(t n ) -y n | + τ |T n |,
which can be rewritten:

|e n+1 | ≤ (1 + τ L)|e n | + τ |T n |.
Now, let us show by induction that, for all n ≥ 0, there holds:

|e n | ≤ Θ n L ((1 + τ L) n -1) (15)
where Θ n := max k=0,...,n |T k |. Indeed, for n = 0, there holds: e 0 = y(0) -y 0 = y 0 -y 0 = 0, so that (15) holds for n = 0. Suppose now that ( 15) is true for n. From this relationship and the bound above, we deduce:

|e n+1 | ≤ (1 + τ L)|e n | + τ |T n | ≤ (1 + τ L) Θ n L ((1 + τ L) n -1) + τ Θ n ≤ Θ n L (1 + τ L) n+1 -(1 + τ L) Θ n L + τ Θ n ≤ Θ n+1 L (1 + τ L) n+1 -1)
where we used the (obvious) bound Θ n ≤ Θ n+1 . We then use the relationship

1 + τ L ≤ exp(τ L),
(remember exp(•) is a convex function) to infer from [START_REF] Sayas | Variational techniques for elliptic partial differential equations[END_REF]:

|e n | ≤ Θ n L (exp(nτ L) -1) ,
and remarking that τ n = t n , we obtain [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF]. This ends the proof.

Let us now apply this result to the forward Euler scheme:

Corollary 1 Let us consider the forward Euler method (2) to discretize Problem (1), and make the assumption that f is globally Lipschitz: there exists L > 0 such that, for all t ∈ R + and u, v ∈ R, there holds:

|f (t, u) -f (t, v)| ≤ L|u -v|, (16)
then there holds, for all n ∈ N:

|e n | ≤ τ 2L e Ltn -1 max t∈[0;tn] |y (t)|. ( 17
)
Proof. It suffices to apply Theorem 1 to the forward Euler scheme, recalling that Φ τ = f , so L is the Lipschitz constant of f . Let us estimate the consistency error, which expression is, for the forward Euler method:

T n = y(t n+1 ) -y(t n ) τ -f (t n , y(t n )).
Then we use the first equation from (2):

y (t n ) = f (t n , y(t n ))
and the Taylor-Lagrange expansion of y at time t n :

y(t n+1 ) = y(t n ) + τ y (t n ) + τ 2 2 y (ξ n ), with ξ n ∈ [t n ; t n+1 ].
We then deduce:

T n = τ y (t n ) + τ 2 2 y (ξ n ) τ -y (t n ) = y (t n ) + τ 2 y (ξ n ) -y (t n ).
Therefore:

T n = τ 2 y (ξ n ).
This ends the proof.

As a result, the forward Euler scheme (2) converges in O(τ ). We say this is a first-order scheme, or a scheme of order one. More generally, we say that a scheme is of order p (p ≥ 0) when convergence is in O(τ p ). The above statement is important. Not only it tolds us that the approximation error vanishes when τ decreases, but it provides us an information with respect to the convergence rate (or speed). For a scheme of order one, notably, this means that dividing by two the time step results in a global error divided by two. This is not so bad, and this makes this scheme already applicable for numerical solving of ODEs. Nonetheless, in terms of accuracy, it remains still among the worst schemes one can think about. This is the reason why we will see in the next part of this chapter some schemes that allow to increase the accuracy. These schemes are not as simple as forward Euler, but they are not that much complicated and the gain in accuracy is sometimes huge.

Before seing other schemes, let us state a few additional remarks corresponding to Theorem 1 and Corollary 1. Particularly, note that, as expected, the bound on the global error increases, exponentially, with time. As a result, the numerical solution gets less and less accurate as time is evolving. This is something you can observe generally in practice. This is a problem for long-time simulations (think about simulation in celestial dynamics, astrophysics, or for spatial technology), which has motivated the need of very accurate schemes. Note also that the bound on the global error depends on the solution itself and its variations: in Corollary 1, it depends on y , so, as it can be expected 'genuinely', an ODE for which the solution has great/irregular variations with time can be more difficult to approximate correctly.

Runge-Kutta (RK) methods

Runge-Kutta (RK) methods are still one-step methods, but for which f is evaluated at some intermediate points. If chosen correctly, these evaluation allow to recover a better accuracy than the forward Euler method. For instance the modified Euler method is defined by the recurrence formula

y n+1 = y n + τ f t n + 1 2 τ, y n + 1 2 τ f (t n , y n ) .
It can be shown, as an exercice, that this method is second-order accurate: the bound in the global error is in O(h 2 ), this time. The gain in accuracy is substantial: dividing by two the time step, the approximation error is divided by four. Plenty of RK methods exist (see for instance the books of Butcher [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF], and Hairer, Norsett & Wanner [START_REF] Hairer | Solving ordinary differential equations[END_REF][START_REF] Hairer | Solving ordinary differential equations[END_REF] for an extensive study of these methods).

Let us just mention the most popular one, the classical fourth-order RK method (see [START_REF] Süli | An introduction to numerical analysis[END_REF]):

y n+1 = y n + 1 6 ỹ1 n + 2ỹ 2 n + 2ỹ 3 n + ỹ4 n ,
where

ỹ1 n = f (t n , y n ), ỹ2 n = f (t n + 1 2 τ, y n + 1 2 τ ỹ1 n ), ỹ3 n = f (t n + 1 2 τ, y n + 1 2 τ ỹ2 n ), ỹ4 n = f (t n + τ, y n + τ ỹ3 n ).
An (awkward) calculation allows to assess that this method is fourth-order: the approximation error is in O(h 4 ). This is very good, still with a method that is quite simple. Remark that the above RK methods are both explicit: y n+1 is obtained directly, without the need to solve an equation involving f .

Linear multistep methods

A drawback of above RK methods is that they may require too many evaluations of f to improve the accuracy. Therefore, another strategy to improve the accuracy is to use previous computations y n-1 , y n-2 , . . . , and not only y n , to predict y n+1 . The interesting point, here, is that these previous computations have already been made, and stored somewhere on your favorite computer. This is the main idea of linear multistep methods, which are, for instance, behind the function lsode in Octave. Following [START_REF] Süli | An introduction to numerical analysis[END_REF], a linear k-step method is defined according to the formula below:

k j=0 α j y n+j = τ k j=0 β j f (t n+j , y n+j ). (18)
In (18) the coefficients α 0 , . . . , α k and β 0 , . . . , β k are real values. Each different choice corresponds to a different method. We will assume that: 1) either α 0 or β 0 is different from 0, 2) α k = 0, so as to avoid degenerate cases. Remark that if β k = 0, the method is explicit: y n+k is obtained directly from previous values of the approximate solution. If β k = 0, the method is implicit: an algebraic equation needs to be solved at each time-step to recover y n+k . To alleviate the notation, we will use the same convention as in the book of Süli & Mayers [START_REF] Süli | An introduction to numerical analysis[END_REF] and set:

f n := f (t n , y n ).
Let us present first some popular examples of linear multistep methods:

• For the choice

β 0 = 1, β 1 = 0, α 0 = -1, α 1 = 1, we get y n+1 -y n = τ f n ,
which is the (explicit) forward Euler scheme (2).

• For the choice

β 0 = 0, β 1 = 1, α 0 = -1, α 1 = 1, we get y n+1 -y n = τ f n+1 ,
which is the (implicit) backward Euler scheme (3).

• For the choice

β 0 = 1 2 , β 1 = 1 2 , α 0 = -1, α 1 = 1, we get y n+1 -y n = τ 2 (f n + f n+1 ),
which is the (implicit) Crank-Nicolson scheme (4).

• For the choice

β 0 = - 9 24 , β 1 = 37 24 , β 2 = - 59 24 , β 3 = 55 24 , β 4 = 0,
and

α 0 = 0, α 1 = 0, α 2 = 0, α 3 = -1, α 4 = 1,
we get

y n+4 -y n+3 = τ 24 (-9f n + 37f n+1 -59f n+2 + 55f n+3 ),
which is the (explicit) Adams-Bashforth four-step scheme.

• For the choice

β 0 = 1 24 , β 1 = - 5 24 , β 2 = 19 24 , β 3 = 9 24 ,
and

α 0 = 0, α 1 = 0, α 2 = -1, α 3 = 1, we get y n+3 -y n+2 = τ 24 (f n -5f n+1 + 19f n+2 + 9f n+3 ),
which is the (implicit) Adams-Moulton three-step scheme.

Note first that this formalism allows to recover the three different schemes presented in 1.1. The last two schemes are new. You already have guessed that the difficulty in the design of sophisticated linear multistep methods is to set appropriately the values of the coefficients (α j ) and (β j ). They must satisfy conditions that make the scheme stable and accurate (consistent and convergent). These last properties are vital: if they are not satisfied, the method is of no use in practice. Let us say a little bit more about these properties.

Stability and accuracy

Let us consider a linear multistep scheme to solve the initial-value problem [START_REF] Allaire | An introduction to mathematical modelling and numerical simulation[END_REF].

Let us now suppose that this problem is well-posed on a bounded time-interval [0, T ], and let N be the natural number such that T = N τ (we suppose that the time-step divides exactly T ). We need k values y 0 , . . . , y k-1 to initiate the computation of the approximate solution using the recurrence formula (18). The first thing to note is that, due to (limited) machine precision, at least, the initial values y 0 , . . . , y k-1 will not be used exactly to compute the approximated solution (y n ), but instead we will only have approximate values: let us call them ŷ0 , . . . , ŷk-1 . They should satisfy, of course:

ŷ0 y 0 , ŷk-1 y k-1 ,
with a precision of, for instance, ten digits (or more). The precision for the initial condition can be very good, but the final approximate solution y n at a time t n may be obtained from these initial values after thousands (or more) calculations. This means that the final solution can be very inaccurate, due to the accumulation of round-off errors. To limit such effects, we want to design linear multistep schemes that verify the following property, known as zero-stability:

The linear k-step method (18) is said zero-stable if there exists a constant K > 0 such that, for any small enough value of the time-step τ , and for any sequence (y n ) generated from (18) and initial conditions y 0 , . . . , y k-1 , and another sequence (ŷ n ), still generated using (18), and initial conditions ŷ0 , . . . , ŷk-1 , there holds

|y n -ŷn | ≤ K max 0≤j≤k-1 |y j -ŷj |, (19)
for any natural number n smaller than N .

Let us introduce now a very useful, and nice, result to assess in practice if a multistep method is zero-stable (or not). This is in fact an algebraic equivalent of zero-stability and is called root condition. For this purpose, let us define the first and second characteristic polynomials associated to (18), as, respectively:

P α (z) = k j=0 α j z j , P β (z) = k j=0 β j z j ,
for z ∈ C. To start with, we need the following lemma:

Lemma 1 Let us consider the recurrence relation

α k u n+k + α k-1 u n+k-1 + . . . + α 1 u n+1 + α 0 u n = 0, ( 20
)
for n a natural number, with the associated characteristic polynomial

P (z) = α k z k + α k-1 z k-1 + . . . + α 1 z + α 0 ,
for z ∈ C. We suppose α 0 = 0 and α k = 0. Let us write P in factorized form

P (z) = α k (z -z 1 ) m1 (z -z 2 ) m2 . . . (z -z l-1 ) m l-1 (z -z l ) m l ,
where z 1 , z 2 , . . . , z l-1 , z l are the l distinct (complex) roots of P , of associated multiplicities m 1 , m 2 , . . . , m l-1 , m l . Then, if a sequence (u n ) satisfies (20), it can be written as follows

u n = l j=1 P j (n)z n j , (21)
for any natural number n. In the above formula, the P j are polynomials of degree m j -1.

We can admit this Lemma for the moment. The proof can be found in the book of Süli & Mayers [START_REF] Süli | An introduction to numerical analysis[END_REF]Lemma 12.1]. We provide now our main result known as the Root Condition, which is a necessary and sufficent condition for zero-stability [START_REF] Süli | An introduction to numerical analysis[END_REF].

Theorem 2 (Root Condition) Let us consider the linear k-step method (18) to solve the initial-value problem (1). Then the method (18) is zero-stable if, and only if, the roots of its first characteristic polynomial P α are all within the domain

S 1 := {z ∈ C | |z| ≤ 1}
with, additionally, each root z j that satisfies |z j | = 1 being simple.

Proof. Let us prove that the Root Condition is a necessary condition.

Let us take the particular case f = 0 in (1), so that we solve y = 0 with initial condition y 0 , which solution is a constant one: y(t) = y 0 for t ≥ 0. In this case, the method (18) reads

k j=0 α j y n+j = 0. ( 22
)
We then apply the above Lemma 1 to infer that the approximate solution (y n ) can be written as

y n = l j=1 P j (n)z n j .
Therefore, if there exists a root z j of the first characteristic polynomial P α that satisfies z j / ∈ S 1 , this means that |z j | > 1 and therefore

|z j | n → +∞
where n tends to +∞. As a result, there exists a set of initial values y 0 , . . . , y k-1 such that |y n | tends to +∞ when n tends to +∞ (i.e., when τ → 0, since n ≤ N and N τ = T ). The other case we need to consider is the case of a root z j on the boundary of S 1 , (|z j | = 1) of multiplicity greater than 1. Since its associated polynomial P j is of multiplicity m j -1 ≥ 1, then it grows with n as n mj -1 , and once again, there exists a set of initial values y 0 , . . . , y k-1 such that |y n | tends to +∞ when n tends to +∞. Let us take now another sequence ŷ0 = 0, ŷ1 = 0, . . . , ŷk-1 = 0 such that the sequence (ŷ n ) generated from these initial values and (22) be identically equal to 0. It appears then that there is no way that the property (19) be satisfied for the two sequences (y n ) and (ŷ n ), when the time-step τ is small enough. Therefore we proved that if the root condition is not satisfied, zero-stability does not hold any more.

The proof that the Root Condition is a sufficent condition is much more technical, and can be admitted. For a detailed proof, see for instance the book of Gautschi, Section 6.3 [START_REF] Gautschi | Numerical analysis[END_REF].

Let us now apply this result to assess that some of the presented linear multistep methods are zero-stable. First, let us consider the forward Euler method [START_REF] Allaire | Numerical linear algebra[END_REF]. The associated first characteristic polynomial is

P α (z) = z -1,
which has one unique simple root z 1 = 1, that verifies |z 1 | = 1. So Theorem 2 ensures that the forward Euler method is zero-stable. The same property holds for the backward Euler scheme (3) and the Crank-Nicolson scheme (4). For the Adams-Bashforth four-step scheme, the first characteristic polynomial is

P α (z) = z 4 -z 3 = z 3 (z -1),
with a triple root z 1 = 0 and a simple root z 2 = 1, so the Root Condition holds, and Theorem 2 ensures that it is zero-stable. You can treat the case of the Adams-Moulton three-step scheme as an exercice.

Let us provide also a counter-example: the following linear multi-step method

y n+3 + y n+2 -y n+1 -y n = 2τ (f n+2 + f n+1 ),
is not zero-stable, though it can look nice at first glance. Indeed, the first characteristic polynomial is

P α = z 3 + z 2 -z -1 = z(z 2 -1) + z 2 -1 = (z + 1)(z 2 -1) = (z + 1) 2 (z -1),
has a double root z 1 = -1 on the boundary of S 1 , so the Root Condition is not satisfied, and Theorem 2 ensures that it can not be zero-stable.

In fact, as you may have noticed in the proof of Theorem 2, the Root Condition is linked with solving the ODE y = 0, with zero right-hand side. This is why we talk about 'zero'-stability. In some sense, a multistep scheme is zerostable if small perturbations in the initial data imply small perturbations in the whole (discrete) trajectory, and this notion of zero-stability can be viewed as a sensibility to initial small perturbations.

For a multistep method to be useful, it not only needs to be (zero-)stable, it should also be accurate. So the other crucial notion for this purpose is that of consistency. Let us take y the unique solution to the initial-value problem (1) and the recurrence formula (18) for a linear k-step scheme. We define the corresponding truncation error as follows:

T n := k j=0 α j y(t n+j ) -τ k j=0 β j f (t n+j , y(t n+j )) τ k j=0 β j . ( 23
)
Though it looks more complicated, one can realize it is a quite straightforward adaptation of the definition for one-step scheme, with still the same idea behind: the exact solution y should satisfy approximately the recurrence relation (18). Note however the normalization by the quantity

τ k j=0 β j (= P β (1)),
that should be supposed different from 0. We now define the consistency property as follow: the multistep scheme (18) is consistent if, for any ε > 0, there exists τ ε > 0, such that, for any value of the time-step τ that satisfies τ ≤ τ , and for any natural number n, there holds

|T n | ≤ ε.
Let us now use some Taylor expansions to rewrite the numerator of (18):

k j=0 α j y(t n+j ) -τ k j=0 β j f (t n+j , y(t n+j )) = k j=0 α j y(t n+j ) -τ k j=0 β j y (t n+j ) = k j=0 α j y(t n ) + (jτ )y (t n ) + O(τ 2 ) -τ k j=0 β j y (t n ) + (jτ )y (t n ) + O(τ 2 ) = y(t n ) k j=0 α j + τ   y (t n ) k j=0 (jα j -β j )   + O(τ 2 ).
where we used first that y is the solution to (1), and then we applied Taylor expansions of y and y at time t n (supposing that y is smooth enough). Finally, we just put together all the terms of the same order in τ . So (necessary and sufficent) conditions for consistency are

k j=0 α j = 0, k j=0 jα j = k j=0 β j ,
which can be rewritten also

P α (1) = 0, P α (1) = P β (1)( = 0).
Remark particularly that any consistent method is associated with a simple root of the first characteristic polynomial, that belongs to the boundary of S 1 (this is, luckily, still compatible with zero-stability).

As an example, still we consider the forward Euler method. It is indeed consistent, since

k j=0 α j = -1 + 1 = 0, k j=0 jα j = 0 + 1 = 1 = 1 + 0 = k j=0 β j .
The same calculations can be done as exercices for the other linear multistep methods provided previously. In fact, we can go further and find the largest natural number p such that T n = O(τ p ). Algebraic conditions on the coefficients (α j ) and (β j ) can be derived following the same path as previously. It only requires Taylor expansions to order higher than 2 when we reformulate the truncation error (23) (see the book of Süli & Mayers [START_REF] Süli | An introduction to numerical analysis[END_REF], Section 12.8, for more details). As for one-step methods, we then talk about methods of order p, and once again, the interest to derive methods of high order is to increase the accuracy of the discrete solution (y n ), with only a slight increase in the computational cost.

We end this section with an important result that allows us to understand why both notions of zero-stability and consistence are fundamental (and sufficent) to assess that the approximated solution (y n ) obtained with the multistep scheme (18) is accurate, in the sense it gets closer to y(t n ) as τ becomes smaller. This is a well-known result in the theory of numerical ODEs, due to Germund Dahlquist. The proof is difficult and will be omitted for this introductory course.

Theorem 3 (Dahlquist's Equivalence Theorem) Let us consider the initial-value problem (1), with f globally Lipschitz, and the multistep scheme (18), that is supposed to be consistent. Let us suppose that the initial values y 0 . . . y k-1 are chosen consistently. Therefore, the multistep scheme (18) is zero-stable, if, and only if, it converges. Moreover, if the exact solution y has continuous derivatives up to order p and the truncation error satisfies T n = O(τ p ), for p ≥ 1, then the global error e n = y(t n ) -y n satisfies

e n = O(τ p ).
The above Theorem 3 ensures the convergence of the five multistep methods we have presented. For the forward Euler method (2) especially, we recover the result obtained thanks to the Corollary 1. We are done now with this first tour about numerical approximation of ODEs. Before entering in the (wonderful) world of numerical PDEs, we just give a few extra comments on this topic.

Extra remarks

This first part about numerical approximation of ODEs is just a small introduction, but can be enough to have a first idea on the topic, and for the reminding part of this course. We will revisit this topic at the end, when we will study PDEs with evolution in time. There are many issues we have not discussed, and if you are interested, there are many good books that deal with the topic (among many others, let us quote [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF][START_REF] Gander | An introduction using Maple and MATLAB[END_REF][START_REF] Hairer | Solving ordinary differential equations[END_REF][START_REF] Hairer | Solving ordinary differential equations[END_REF][START_REF] Süli | An introduction to numerical analysis[END_REF]). In fact, numerical analysis of ODEs has been object of intensive research along the twentieth century, and many methods have been invented and analysed. This is still an active research topic, though.

To give you a few hints about other issues, first, there is also an important result by Germund Dahlquist, the so-called Dalquist barrier, which tells you, roughly speaking, that the condition of zero-stability limits the multistep methods in terms of accuracy [START_REF] Süli | An introduction to numerical analysis[END_REF]. Also, another issue are stiff systems of ODEs, which present some specific difficulties in terms of numerical approximation, and for which we need a stronger notion of stability (A-stability) and specific multistep schemes (BDF schemes for instance). A good and simple introduction to this issue can be found, still, in the book of Süli & Mayers [START_REF] Süli | An introduction to numerical analysis[END_REF]Section 12.11], and one can have a look at [START_REF] Hairer | Solving ordinary differential equations[END_REF] for an extensive treatment. Also, a fascinating topic is about geometric numerical integration and numerical schemes that preserve some important properties for some classes of dynamical systems (symplectic numerical integration for instance). An introduction to these issues can be found in the book of 'Geometric Numerical Integration' of Ernst Hairer, Gerhard Wanner and Christian Lubich [START_REF] Hairer | Geometric numerical integration[END_REF].

FEM for elliptic PDEs

Let us move now from ODEs to Partial Differential Equations (PDEs). PDEs allow to model a great range of physical phenomena very accurately, since common physical laws, such as the conservation of mass, momentum, energy, can be expressed using PDEs. Basically, and as you know, PDEs relate some unknown quantities (the temperature, the motion of material points, etc) to their derivative in time and space. The most famous PDEs are known sometimes since the eighteen or nineteen centuries, as well as closed-form solutions for very specific configurations. Various approximation techniques have been invented from the end of the nineteen century, so as to calculate effective approximate solutions. Now, the most popular methods are the finite difference method, the finite element method, the finite volume method and the spectral method. An introductive presentation to such methods is given for instance in the book of Martin Gander and Felix Kwok [START_REF] Gander | Numerical analysis of partial differential equations using Maple and MATLAB[END_REF].

As an introduction to this topic, we will focus here on the finite element method, and on the most simple elliptic partial differential equation (a linear, scalar, PDE). We first present it, and then reformulate it to find what is called its weak formulation. Indeed, the finite element method does not discretize directly the differential operators. In fact, it relies on the construction of an appropriate finite element space, which is a finite dimensional vector space, that approximates the true vector space in which we seek the solution (which is an infinite dimensional vector space, endowed with a Hilbert (or Banach) structure). At the end of the whole process, we will simply obtain an algebraic linear system, that can be implemented and solved using a computer. Finally, we will answer the question : does the finite element method provide an accurate approximation of the exact solution ?

Model problem: reaction-diffusion

Let us introduce a first model problem, as simple as possible. A PDE provides a relationship between the different derivatives, in time and space, of an unknown. Let us forget about time for the moment and consider only PDEs in space, that correspond to stationnary problems (the solution does not evolve in time). In general (but not always) the solution of the PDE is defined only on a subdomain of the whole space. Let us denote by Ω this subdomain. We will suppose that Ω is an open subset of R d , d ≥ 1 (in practice, d = 1, 2 or 3). We will suppose that Ω is non-empty and bounded. We will denote by Γ := ∂Ω the boundary of Ω. For d = 1 we will take Ω as a finite union of open intervals. For d = 2 (respectively d = 3) the boundary Γ will be supposed polygonal (respectively polyhedral). This last assumption can be alleviated in fact, but will avoid some technical difficulties at a specific moment. We do not make any extra assumption about Ω, and, for instance, it does not need to be connected or convex.

We want to solve the following problem:

c u -∆u = f in Ω, u = 0 on Γ. (24)
In (24) u : Ω → R is the unknown solution, c ≥ 0 is a constant, f : Ω → R is a (known) source term. The operator ∆ is the Laplace operator, defined as

∆v := div (∇v) = ∂ 2 ∂x 2 1 + • • • + ∂ 2 ∂x 2 d v,
for any v : Ω → R smooth enough. Note that we used the notation x 1 , . . . , x d for the coordinate system in R d . The above equation cu -∆u = 0 is known as the reaction-diffusion equation, and is among the simplest partial differential equations. It is a linear, scalar, elliptic partial differential equation. For c = 0, we recover the Laplace equation (or Poisson's problem), that has many physical interpretations: for instance the unknown u can represent the (steady) temperature induced by the source term f , or the vertical displacement of a thin membrane induced by f , which in this case represents a surfacic density of loads. For c > 0 there is a reaction term cu can have a special meaning for some problems in chemistry or biology, and also will be useful when we will study time-evolving PDEs.

Remark also that in Problem (24) there is a second equation on the boundary Γ. This is a boundary condition. Since it involves the unknown u directly, it is called a Dirichlet (or essential) boundary condition. There are other types of boundary conditions (Neumann, Robin), that we should study later on. This condition is in general motivated by physical considerations : for instance, if u represents the displacement of a membrane, it means that the membrane is fixed at its boundary (u = 0 means 'no vertical displacement'), or if u represents a temperature, it means that the temperature is fixed (to 0 K) on the boundary.

Weak formulation

The first important idea of the FEM is to reformulate Problem (24) as a weak problem, which is equivalent in some specific sense. To this purpose we pick an arbitrary test function v : Ω → R, and multiply it to the first equation of (24):

c uv -(∆u)v = f v,
and we integrate over the whole domain Ω:

c Ω uv - Ω (∆u)v = Ω f v.
We then apply the Green formula to the second term:

c Ω uv + Ω ∇u • ∇v - Γ (∇u • n)v = Ω f v.
Here n denotes the outward unit normal to the boundary Γ. Finally, we impose an extra condition on the test function v. It should vanish on the boundary Γ. This causes the boundary term above to vanish, and we get finally:

c Ω uv + Ω ∇u • ∇v = Ω f v.
Before we go further, I guess that now, many questions should arise. Why did we do that ? Are the above calculations mathematically meaningful ? Why such a complicated way to proceed ? We will try to provide answers to these questions in the sequel.

First of all, just suppose that u and v are smooth enough so that above calculations are valid. Then, let us try to provide, at least, a precise meaning to the last equation. For this purpose, we need two functional spaces, first let us introduce the Sobolev space

H 1 (Ω) := {v ∈ L 2 (Ω) | (∇v) ∈ L 2 (Ω; R d )},
where L 2 (Ω) is the (Lebesgue) space of square-integrable (real-valued) functions, and L 2 (Ω; R d ) is the space of square-integrable vector-valued functions. Take care that we need to take the gradient of v, ∇v in the sense of distributions (or in the weak sense) since functions in L 2 (Ω) are not necessarily differentiable (in the classic sense). We introduce for functions v, w in H 1 (Ω) a scalar product, defined as (v, w) 1,Ω := (v, w) 0,Ω + (∇v, ∇w) 0,Ω , where (v, w) 0,Ω := Ω vw denotes the standard scalar product in L 2 (Ω). The corresponding norm will be denoted by • 1,Ω :

v 1,Ω := (v, v) 1,Ω
for v ∈ H 1 (Ω). It can be easily proven that, endowed with this scalar product, H 1 (Ω) is a Hilbert space, i.e., a vector space with a scalar product, in which every Cauchy sequence converges (completeness). This is in fact a property inherited from the structure of L 2 (Ω) [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF][START_REF] Sayas | Variational techniques for elliptic partial differential equations[END_REF] We define now H 1 0 (Ω) which is the closure of the space of test functions D(Ω) in H 1 (Ω), for the norm • 1,Ω . It can be proven that

H 1 0 (Ω) := {v ∈ H 1 (Ω) | v = 0 on Γ}.
This is a non-trivial result, and an important consequence of the trace theorem (see as well the references cited above), since it is not at all obvious at first glance that functions in H 1 (Ω) are properly defined on the boundary Γ. Particularly, the trace theorem needs regularity assumptions on the boundary, and if Γ is a Lipschitz (or strictly) polygonal boundary, this theorem can be applied. The above characterization of H 1 0 (Ω) is much more convenient than its original definition, and will serve to treat the Dirichlet boundary conditions. We can now define the bilinear form

a(v, w) := c Ω vw + Ω ∇v • ∇w,
for every functions v and w in H 1 (Ω). We check that indeed, a(•, •) is meaningful when the functions are restricted to this space: applying the Cauchy-Schwarz inequality, we assess indeed the Lebesgue integrals are finite. As well, we suppose that f ∈ L 2 (Ω), and for every w in H 1 (Ω), we define the linear form L(w) := Ω f w, and we check once again that this definition is meaningful.

We have worked enough now to define the weak formulation (or variational formulation) of Problem (24), which is

Find u ∈ H 1 0 (Ω) such that a(u, v) = L(v), for any v ∈ H 1 0 (Ω). ( 25
)
Remark that the boundary condition u = 0 on Γ has been integrated in the definition of the space H 1 0 (Ω) (we say it is an essential boundary condition). Since the bilinear form a(•, •) is symmetric, this weak formulation is the firstorder optimality condition associated to the minimisation problem

Find u ∈ H 1 0 (Ω) that minimizes J over H 1 0 (Ω), (26)
where

J (v) := 1 2 a(v, v) -L(v).

Lagrange finite element space on simplices

Of course, Problem (24) (or (25)) can be solved analytically only for very special cases, where for instance Ω is square-shaped and the expression of f is simple. Otherwise, there is little hope to find an exact solution (and the situation is even worse for more complex PDEs, with more complex boundary conditions). So we introduce an approximation technique known as the Lagrange Finite Element Method. The first step, as in the previous chapter about ODEs, is to decompose the domain Ω into small, simple cells. This process, in the FEM terminology, is called, meshing. So we introduce a simplicial mesh as a collection of simplices in R d that partition the domain Ω. We call each simplex K and K h the collection of all the simplices in Ω (we will precise in a few minutes the meaning of the index h). The mesh needs to satisfy two conditions:

1. It should cover exactly the domain Ω:

Ω := K∈K h K,
where the simplices K are supposed to be closed.

2. The intersection of two distinct simplices K and K is either empty, or a simplex of dimension lower than d, that is common to both K and K .

For each simplex K in the mesh, h K denotes the maximal distance between two points inside K (it is the 'size' of K), and

h := max K∈K h h K
is the global size of the mesh. What we expect now is that, if the mesh size is small enough, we will recover an accurate enough approximated solution (as for the time-step τ when we studied ODEs). Before going further, let us be more explicit about these notions for problems in one or two dimensions. First in one dimension (d = 1) and for Ω that is an open interval, the different simplices K 1 , . . . , K N +1 are simply closed intervals, that cover Ω, and intersect each other only at both ends. We denote

x i = K i ∩ K i+1
the nodes of the mesh. Remark that if all the simplices are of the same size, we recover a situation similar as in the previous chapter. Now, for problems in two dimensions, d = 2, simplices K are triangles. The intersection between two triangles is either empty, or a common vertex, or a common edge. The vertices of the triangles are also called the nodes of the mesh. Now, from the mesh K h , we will build a finite dimensional subspace of H 1 0 (Ω), that we will call V h , as follows. We define, now for every value of the dimension d, the nodes of the mesh as the vertices of the simplices. The interior nodes are the nodes that do not belong to the boundary Γ. Let us denote them by x 1 , . . . , x N . Let us define a family of basis functions ϕ 1 , . . . , ϕ N such that 1. Each ϕ i is globally continuous on Ω and its restriction to each simplex K is a piecewise affine function (a polynomial function of degree at most 1):

ϕ i | K ∈ P 1 (K),
where P 1 (K) is the space of (multivariate) polynomials of degree 1.

2. Each ϕ i is equal to 1 at node x i , and equal to 0 for the other nodes:

ϕ i (x j ) = δ ij ,
where δ ij is the Kronecker symbol.

From the two above properties, one can check directly that each ϕ i has a compact support located on the patch ω i of simplices that share the same vertex x i . So this function vanishes everywhere except in simplices that are in the neighbourhood of the node x i .

One can also check that the functions (ϕ i ) are linearly independent. Then we set

V h := span(ϕ 1 , . . . , ϕ N ).

In dimension 1, (ϕ i ) is the so-called family of hat functions, and in dimension 2, each ϕ i is a pyramid of basis ω i and height 1. One can also check that V h is the space of continuous functions that are affine on each simplex K, and that this is a proper subspace of H 1 0 (Ω): because of this property, we talk about conforming finite elements. The terminology 'Lagrange' finite elements comes from the fact that it relies on piecewise polynomial Lagrange interpolation between the nodes of the mesh. Remark finally that any function v h in V h satisfies automatically the boundary condition v h = 0 on Γ.

Once the space V h is built, the discrete problem that serves to approximate Problem (25) simply reads:

Find u h ∈ V h such that a(u h , v h ) = L(v h ), for any v h ∈ V h . ( 27
)
This is our Finite Element Method. You see, this is easy, we just took Problem (25) and added the subscript h everywhere. Seriously speaking, the main difference with ( 25) is that now, we find an approximate solution in a simple, finite dimensional, vector space V h . The main interest of ( 27) is the following: it can be recasted in matrix form. Indeed, first write

u h = N j=1 U j ϕ j
and set U the column vector of size N that contains all the scalar unknowns (U j ). Remark, that, due to the definition of the basis (ϕ i ), we get the relationship

U j = u h (x j ),
for any index j. We say that the unknowns are nodal unknowns, since they simply represent the value of u h at interior nodes. Take v h = ϕ i ∈ V h for an arbitrary index i, and we get from (27):

a(u h , ϕ i ) = L(ϕ i ).
Then use the above expression for u h :

a   N j=1 U j ϕ j , ϕ i   = L(ϕ i ).
By linearity

N j=1 U j a (ϕ j , ϕ i ) = L(ϕ i ). ( 28
)
There remains to introduce the (symmetric) stiffness matrix K = (K ij ) with

K ij = Ω ∇ϕ i • ∇ϕ j = Ω ∇ϕ j • ∇ϕ i = K ji
as well as the (symmetric) mass matrix

M ij = Ω ϕ i ϕ j = Ω ϕ j ϕ i = M ji
for i, j = 1, . . . , N . We define the load vector F = (F i ) with the relationship

F i = L(ϕ i ) for i = 1, . . . , N . Since there holds, by definition, a(ϕ i , ϕ j ) = cM ij + K ij , formula (28) becomes N j=1 (c M ij + K ij ) U j = F i .
So we finally reformulate the discrete Problem (27) in matrix form:

(c M + K) U = F. ( 29
)
At this point, we get something that can be treated by a computer. So the whole story about FEM, in practice, consists, in 1) building the matrices M and K as well as the vector F (assembly procedure), 2) solving the linear system (29), 3) recovering u h from the vector U of nodal unknowns (postprocessing).

Practical implementation

There have been plenty of available finite element codes for decades. They differ in various aspects, but are more or less based on the same structure and ideas.

To start with, and to stay within the (now friendly) Octave world, we will play a little bit with the Octave scripts of Martin Gander and Felix Kwok, available here:

https://www.unige.ch/ ~gander/siamprograms.php

If you have a look at all the scripts, and at the detailed explanations, you will realize that solving a PDE with FEM is rather more complicated than solving an ODE. Let us start just giving the main ideas.

1st step: mesh generation. The first thing is to get a mesh K h of the domain Ω. A simple and convenient way to handle it is to store it as a matrix of mesh nodes N, which contains the coordinates of all the nodes, ordered by columns, and another matrix of simplices T that contains the information about connectivity. In fact, each line of T contains the indices of the nodes that belong to the same simplex. It can also contain additional information, such as special physical properties, or information about the boundary. To generate a mesh that corresponds to various shapes, simply use the function NewMesh that returns, at least, N and T. There are also functions to plot, refine and smooth a mesh.

2nd step: assembly of the linear system (29). This is the most important part and, in some sense, the very core of every finite element program. In fact, the most efficient, and now widespread way to compute the matrices involved in (29) is to go through all the simplices, and add the corresponding contribution to the global matrices. This way of thinking comes from the following observation. Recall that, for the stiffness matrix, we have

K ij = Ω ∇ϕ i • ∇ϕ j = K∈K h K ∇ϕ i • ∇ϕ j .
So it suffices to compute only integrals on simplices K that involve restrictions on K of the basis functions. Remark that for a wide range of values for i and j, these integral vanish. The precise way to do it is described in details in the corresponding chapter in the book of Martin Gander and Felix Kwok [START_REF] Gander | Numerical analysis of partial differential equations using Maple and MATLAB[END_REF].

The function FEPoisson does the whole job, and calls subroutines Compute-ElementMassMatrix and ComputeElementStiffnessMatrix, that compute elementary matrices (associated to each simplex) that will contribute into the global matrix.

3rd step: solving of the linear system (29). In general, one just call a library for solving linear systems. In Octave, this is done automatically and efficiently through the instruction \. In general, direct methods such as LU or Cholesky decomposition or variants are used for small systems, and iterative methods such as CG or GMRES are used for large systems [START_REF] Allaire | Numerical linear algebra[END_REF][START_REF] Gander | An introduction using Maple and MATLAB[END_REF]. Still the function FEPoisson does this part.

4th step: postprocessing. Once U is recovered, one can provide the user some various information about the solution. At least, it should be visualized. This is the object of the function PlotSolution.

Existence of solutions

Now that we have a first insight into the method and its behavior, let us try to justify it, particularly to justify that u h is a good approximation of u if h is small enough and if the simplices have a nice shape (we will put a precise definition behind these fuzzy words soon). The first step is to make sure that the exact solution u exists, and is unique, so that we approximate something well-defined. For ODEs, we had the Cauchy-Lipschitz (or Picard) Theorem that ensured this fundamental point. Here, we need an equivalent. Before we prove this main theorem of existence, we need a preliminary result, which is known as Poincaré's inequality.

Theorem 4

Let Ω be an open, connected and bounded set in R d with polygonal boundary. There exists c P > 0 such that, for every v ∈ H 1 0 (Ω):

c P v 0,Ω ≤ ∇v 0,Ω . ( 30
)
As a consequence, we deduce that ∇ • 0,Ω is a norm on H 1 0 (Ω), equivalent to the H 1 -norm.

Proof. We proceed by contradiction. Suppose that, for every n ≥ 1, there exists v n ∈ H 1 0 (Ω) such that 1 n v n 0,Ω ≥ ∇v n 0,Ω .

Define w n := v n / v n 0,Ω , so there holds

1 n ≥ ∇w n 0,Ω ,
with w n 0,Ω = 1. As a result, there holds

w n 2 1,Ω = w n 2 0,Ω + ∇w n 2 0,Ω ≤ 1 + 1 n 2 ≤ 2,
so that the sequence (w n ) is bounded in the H 1 -norm. As a result, we use the Rellich-Kondrachov Theorem, and deduce there is a subsequence (w ϕ(n) ) that converges strongly in L 2 (Ω). Moreover, since

∇w ϕ(n) 0,Ω ≤ 1 ϕ(n) ,
we deduce that ∇w ϕ(n) tends to 0 for n → +∞, and thus (w ϕ(n) ) is a Cauchy sequence in H 1 0 (Ω). Since H 1 0 (Ω) is a Hilbert (complete) space, we deduce that there exists a limit w such that

w ϕ(n) → w (strongly) in H 1 0 (Ω).
There holds ∇w 0,Ω = lim n→+∞ ∇w ϕ(n) 0,Ω = 0, from which we deduce that w is a constant on Ω. Furthermore, since w = 0 on ∂Ω, we conclude that w = 0. On the other hand, we have

w 0,Ω = lim n→+∞ w ϕ(n) 0,Ω = 1,
which contradicts w = 0. To prove the last assertion, note first that we have, by definition, for v ∈ H 1 0 (Ω):

∇v 0,Ω ≤ v 1,Ω .
Then we use (30) to prove that

v 2 1,Ω = v 2 0,Ω + ∇v 2 0,Ω ≤ 1 c 2 P ∇v 2 0,Ω + ∇v 2 0,Ω ≤ 1 c 2 P + 1 ∇v 2 0,Ω .
This ends the proof.

Beware that the above Theorem makes use of three non-trivial results in functional analysis:

1. The Rellich-Kondrachov Theorem, that ensures, below the appropriate assumptions on Ω, that the injection H 1 (Ω) → L 2 (Ω) is compact (particularly, Ω needs to be bounded).

2. The kernel of the (weak) gradient operator ∇ consists of the constant distributions on a connected domain (a result that is difficult to establish when d ≥ 1).

3. The Trace Theorem, which ensures that H 1 0 (Ω) is the kernel of the trace operator.

There is a direct (simplest ?) proof of Poincaré's inequality, which is obtained thanks to explicit calculation of the Sobolev norms [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF][START_REF] Sayas | Variational techniques for elliptic partial differential equations[END_REF]. Nevertheless, the above proof has the advantage that it can be generalized easily for more complex boundary conditions (not only homogeneous Dirichlet). Some generalizations can be found in the litterature with the name of Poincaré-Friedrichs inequalities, Peetre-Tartar lemma, etc.

We now state our main theorem of existence:

Theorem 5 Let Ω be an open bounded, connected, set in R d with polygonal boundary. Let f ∈ L 2 (Ω), for any c ≥ 0, there exists one unique solution u ∈ H 1 0 (Ω) to Problem (25), that is also the unique minimizer of J on H 1 0 (Ω). Moreover, the following a priori estimate holds: there exists C > 0 such that

u 1,Ω ≤ C f 0,Ω . (31)
Proof. We use the previous result, Theorem 4, which ensures that the bilinear form a(•, •) defines a scalar product on H 1 0 (Ω), even for c = 0, equivalent to the canonical scalar product (•, •) 1,Ω . Since H 1 0 (Ω) is a Hilbert space, we can apply the Riesz-Fréchet representation theorem, to assess that there is one unique

vector u ∈ H 1 0 (Ω) such that a(u, •) = L(•).
Moreover, there holds, for any v in H 1 0 (Ω):

J (v) -J (u) = 1 2 a(v -u, v -u).
This ensures that u is the unique minimizer of J on H 1 0 (Ω). Let us derive finally the a priori error estimate. Remember first that

u 2 1,Ω ≤ 1 c 2 P + 1 ∇u 2 0,Ω ≤ 1 c 2 P + 1 a(u, u) = 1 c 2 P + 1 L(u).
Using the Cauchy-Schwarz inequality, we have also

L(u) ≤ f 0,Ω u 1,Ω .
Therefore the estimate (31) is deduced. Remark that the constant C does not depend of c.

In the above proof, we have established that there exists α > 0 such that, for any v ∈ H 1 0 (Ω), there holds:

α v 2 1,Ω ≤ a(v, v).
When such property is satisfied, we say that the bilinear form a(•, •) is elliptic (or coercive) on H 1 0 (Ω). Having such a property (or a variant) is a fundamental point to prove existence and uniqueness when one looks at PDEs. From the viewpoint of convex analysis, using the ellipticity of a(•, •), we deduce that the functional J is strongly convex on H 1 0 (Ω) (and we can recover the existence and uniqueness of solutions following this path).

Let us now establish the strong-weak equivalence, in other terms let us prove that Problem (24) and Problem (25) are equivalent in some sense.

Proposition 1 Any solution u ∈ H 1 0 (Ω) to Problem (25) satisfies also the Problem (24), where the first equation has to be understood in the distributional sense. Reciprocally, if u ∈ H 1 0 (Ω) satisfies (24) (in the sense of distributions), it is also a solution to (25).

∇u, ∇ϕ = Ω f ϕ -c Ω uϕ so (cu -∆u -f ) , ϕ = 0, which means that cu -∆u -f = 0 in D (Ω). Moreover, since u ∈ H 1 0 (Ω)
, there also holds u = 0 on the boundary Γ. As a result, u is a solution to Problem (24). The other assertion can also be proven using the density of D(Ω) into H 1 0 (Ω), as we will show now. Let us suppose that u ∈ H 1 0 (Ω) satisfies

cu -∆u = f in D (Ω). Then pick v ∈ H 1 0
(Ω) and a sequence of test functions (ϕ n ) that tends to v in the H 1 -norm. First, note that, for any n ≥ 0, we can bound

|a(u, v) -a(u, ϕ n )| ≤ (c + 1) u 1,Ω v -ϕ n 1,Ω
where we used the Cauchy-Schwarz inequality. Similarly, we bound

|L(v) -L(ϕ n )| ≤ f 0,Ω v -ϕ n 1,Ω .
There just remains to write

a(u, v) -L(v) = (a(u, v) -a(u, ϕ n )) + (a(u, ϕ n ) -L(ϕ n )) + (-L(v) + L(ϕ n )).
Since cu -∆u = f in D (Ω), the second term a(u, ϕ n ) -L(ϕ n ) vanishes, and we conclude using the triangle inequality, the above bounds, and letting n → +∞.

Following the same proof as above, we deduce this next result for the FE solution u h :

Theorem 6 Let Ω be an open bounded, connected, set in R d with polygonal boundary. Let V h a finite dimensional space such that V h ⊂ H 1 0 (Ω). Let f ∈ L 2 (Ω)
, for any c ≥ 0, there exists one unique solution u h ∈ V h to Problem (27), that is also the unique minimizer of J on V h . Moreover, the following a priori estimate holds: there exists C > 0 such that

u h 1,Ω ≤ C f 0,Ω . (32)
Moreover, from the ellipticity of a(•, •) on V h (⊂ V ), we deduce that, in (29), the global matrix (c M + K) is symmetric, positive and definite (even for c = 0).

Accuracy of the FE approximation

Now we are ensured that both solutions u and u h exist, and are not multiple, we can say something about the accuracy of the FE approximation, and we will show that the norm u -u h 1,Ω will decrease when the size h of the mesh decreases, provided some specific assumption on the geometry of the simplices. We will start with an abstract result, known as Cea's Lemma.

Lemma 2 Let V be a Hilbert space and V h a closed vector space such that V h ⊂ V . Let a(•, •) be a continuous bilinear form on V that satisfies the following ellipticity condition α v 2 1,Ω ≤ a(v, v). Let L(•) be a continuous linear form on V . We define u the unique solution of the weak problem

a(u, v) = L(v), ∀v ∈ V,
and u h the unique solution of the approximate weak problem

a(u h , v h ) = L(v h ), ∀v h ∈ V h .
Then there exists a constant C > 0 that depends solely on the continuity and ellipticity constants of a(•, •) such that:

u -u h V ≤ C inf v h ∈V h u -v h V (33)
where • V is the natural norm on V .

Before we study the proof, let us comment a little bit on this lemma. The first point is that there are very few assumptions on V h , except that it must be a subspace of V . It has not to be something specific like a finite element space (it does not even need to be a finite dimensional vector space). The second point is that you bound the difference between the continuous and discrete errors by an approximation error that quantifies how accurately the space V h approximates V (the richer the space is, the better the approximation should be).

Proof. Let u and u h be the continuous and discrete weak solutions. First, pick

v h ∈ V h and, noting that v h -u h ∈ V h ⊂ V , we can write a(u, v h -u h ) -a(u h , v h -u h ) = L(v h -u h ) -L(v h -u h ) = 0 and thus a(u -u h , v h -u h ) = 0.
Using the above equality we can write

a(u -u h , u -u h ) = a(u -u h , u -v h + v h -u h ) = a(u -u h , u -v h ).
We then use the ellipticity of a(•, •), combined with the above identity:

α u -u h 2 V ≤ a(u -u h , u -u h ) = a(u -u h , u -v h ).
There remains to use the continuity of a(•, •):

α u -u h 2 V ≤ C a u -u h V u -v h V ,
with C a > 0 the continuity constant of a(•, •). As a result we get

u -u h V ≤ C a α u -v h V ,
and we conclude by setting C := Ca α and taking the infimum over V h . Now to say something more precise about the (concrete) Lagrange finite element space V h we defined and used in the previous sections, we need some technical results that allow to bound the error for the Lagrange interpolation in appropriate Sobolev norms. Before this, we need now to introduce a geometric assumption on the mesh, and more precisely on the simplices of the mesh.

Let us consider a mesh sequence (K h ) h>0 index by h (the size of the simplices). We say that this sequence is shape-regular if there is a constant C R > 0, independent of h, such that, for any value of h, and for every simplex K ∈ K h , there holds

h K ρ K ≤ C R ,
where h K is the size of K, as defined previously, and ρ K is the size of the largest ball included in K.

Remark first that we consider a mesh sequence, and not only one fixed mesh, such that the maximum size of the simplices tend to 0. This definition is traditionnaly attributed to Philippe G. Ciarlet (one of the pioneers of this method). It means, roughly speaking, that the geometry of the simplices can not degenerate when the mesh is refined. For d = 2, when simplices are triangles, it can be shown that this is equivalent to the minimal angle condition: the smallest angle in each triangle is bounded by below (the triangles can not be "too flat"). For structured mesh, such a geometric condition is always satisfied, but, for unstructured mesh, and complicated geometries, this condition may be hard to achieve in practice (this is still the object of current research). Of course, for d = 1, this geometric condition is always satisfied (h K = ρ K . . . ).

We now introduce the Lagrange interpolation operator as follows:

I h (v) := N i=1 v(x i )ϕ i ,
for any v ∈ H 1 0 (Ω) ∩ C 0 (Ω). We impose this condition on v since continuity is required to define nodal (ponctual) values of v, and functions in H 1 0 (Ω) are not always continuous (at least when d ≥ 2). Remark that I h (v) belongs to the Lagrange finite element space V h : this function is simply obtained from v by setting its value at interior nodes, and then doing piecewise polynomial interpolation between nodes.

The following estimate can be proven for the linear Lagrange interpolation: Lemma 3 Let d ≤ 3. Let V h be the linear Lagrange finite element space built from a shape-regular mesh K h of the domain Ω. Then, there is a constant C > 0, independent of h (but dependent of the shape regularity constant C R ), such that, for every v ∈ H 2 (Ω) ∩ H 1 0 (Ω), there holds:

v -I h (v) 1,Ω ≤ Ch v 2,Ω .
The Sobolev space H 2 (Ω) denotes the subspace of L 2 (Ω) with partial weak derivatives up to order two that belong to L 2 (Ω). It can be proven that functions in H 2 (Ω) are always continuous for d ≤ 3 (this is a consequence of the Sobolev embedding Theorems). The proof of this result in the general case (d ≥ 2) is a bit harsh and can be admitted at this stage (for an excellent detailed proof, see for instance [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF]). In the case d = 1, the proof is much more simple [START_REF] Gander | Numerical analysis of partial differential equations using Maple and MATLAB[END_REF].

We end this section with the following result, that bounds the error related to the Finite Element solution u h , as a function of the mesh size h: Theorem 7 Let d ≤ 3, let (K h ) h>0 be a shape regular sequence of meshes, let V h be the Lagrange finite element space built upon these meshes. We still denote by u the solution to Problem (25) and by u h the solution to Problem (27). We suppose that u ∈ H 2 (Ω) ∩ H 1 0 (Ω). Then there exists a constant C > 0, that depends only on the continuity and ellipticity constants of a(•, •), and on the shape-regularity constant C R , such that

u -u h 1,Ω ≤ Ch u 2,Ω
for every h > 0.

Proof. We set V = H 1 0 (Ω) and V h ⊂ V the Lagrange finite element space. Then the result is obtained using Cea's Lemma 2, setting v h = I h (u), and applying the interpolation estimate in Lemma 3.

The above result is of fundamental importance. First, it allows to assess, that, for shape-regular meshes, the approximate FE solution u h tends to the solution u when the mesh size h tends to 0 (when the simplices get smaller and smaller). Moreover, the H 1 -norm is a quite strong one, since it takes also into account the spatial derivatives of the solution. Also, it quantifies precisely how much the quality of the approximation is improved when the mesh is refined. If h is divided by two, so is the approximation error u -u h 1,Ω . We talk of linear convergence in this situation. This behavior can be illustrated very accurately in practice in many situations.

Extra remarks

This presentation was a tiny introduction to this wide topic, and many aspects of the method were hidden, mostly to simplify the presentation, and also because we do not have infinite time to discuss about this topic.

First of all, the reaction-diffusion, or even Poisson's problem, is a good, and common choice to start with, but there are many more partial differential equations, with a great variety of boundary conditions. Sometimes, these more complex problems can be treated by rather straightforward adaptations of the present framework, sometimes not, and in this case, it has been object of intensive research activity.

Also, many type of finite elements, with specific attractive features, have been invented since 1950. There is even a Table of Finite Elements designed by Douglas Arnold: http://www-users.math.umn.edu/ ~arnold/femtable/ The Lagrange finite element is just the simplest one (and still among the most widely used). There are also many methods inspired from the finite element method (just to quote a few recent ones: the eXtended Finite Element Method, the Isogeometric Method, the Virtual Element Method, the Hybrid High Order Method, the cut Finite Element method, etc, etc).

Finally, one can go much further into the mathematical analysis, and prove much more things about the properties of the numerical approximation: estimates in various norms can be derived, as well as, for instance, a posteriori error estimates that allow to determine on which cells of the mesh the error is the larger. This allows for instance adaptive mesh refinement, which is a powerful technique to improve the accuracy of the solution.

Concerning all these different topics, there are very good monographs that review most of them in details, particularly [START_REF] Allaire | An introduction to mathematical modelling and numerical simulation[END_REF][START_REF] Brenner | The mathematical theory of finite element methods[END_REF][START_REF] Ern | Theory and practice of finite elements[END_REF][START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF].

Time-dependent PDEs

In this final chapter, we will illustrate how the notions we have studied in the two previous chapters can be combined to deal with time-dependent PDEs, that allow to model unsteady phenomena.

The heat equation

Let us introduce a simple model of PDE with an evolution in time: the heat equation. As in the previous chapter, we need a domain Ω, which is still an open subset of R d , d ≥ 1 (in practice, d = 1, 2 or 3). We will suppose again that Ω is non-empty and bounded. We will denote by Γ := ∂Ω the boundary of Ω. For d = 1 we will take Ω as a finite union of open intervals. For d = 2 (respectively d = 3) the boundary Γ will be supposed polygonal (respectively polyhedral). We consider the evolution of the temperature u in the domain Ω, so we want to solve the following heat equation:

∂u ∂t -κ ∆u = f in R + * × Ω, u = 0 on R + * × Γ, u(0) = u 0 in Ω. (34)
In the above equation ( 34)

u : R + × Ω → R
is the unknown solution (the temperature). We will denote by u(t, x) its (scalar) value at time t ≥ 0 and at point x ∈ Ω. The Laplace operator is still denoted by ∆, whereas the notation ∂ ∂t denotes the partial derivative in time. As for the reaction-diffusion equation, f : R + * × Ω → R is a given source term, that models the local production, or absorption, of heat. The coefficient κ > 0 in front of the Laplace operator is a diffusion coefficient, that models the conductive properties of the material under consideration. Still we impose a Dirichlet condition u = 0 on the boundary, for the sake of simplicity. It means that we know the value of the temperature u on Γ. This needs to be changed if a different situation is considered. Conversely to the reaction-diffusion equation, and in the same way as for ODEs, we need an initial condition, and we suppose the temperature is known at the initial time t = 0, and equal to u 0 , a given function.

Remark that if we seek steady (stationnary) solutions to this equation (if we set u(t, x) = u(x)) the partial derivative in time vanishes, and we recover the Poisson's equation we have seen before (reaction-diffusion with c=0). This is the reason why we told you earlier that Poisson's equation can also model the diffusion of heat in a continuous medium.

As for most of the other PDEs, and as for the reaction-diffusion equation, closed-form solutions are only known for very particular domains and/or values of the source term, initial conditions and boundary conditions. Notably, a well-known solution, due to Joseph Fourier, can be obtained using the Fourier transform (this motivated the introduction of the Fourier transform). Nevertheless, this solution is valid only for Ω = R d . So for most practical situations, we need to discretize Problem (34) and to find an approximate solution with our (now) favorite computer. Before doing this, and as for the reaction-diffusion equation earlier, we need first to derive a weak form of this equation.

For this purpose, it will be convenient to adopt the following convention, that will allow us to see the space-time field u as a function of time, only:

u : t → u(t) ∈ Ω so we can say that u(t)(x) = (u(t))(x) = u(t, x).
This said and done, we proceed exactly as for the heat equation, at a fixed time t > 0. We, once again, pick an arbitrary test function v : Ω → R, and multiply it to the first equation of (34):

∂u ∂t v -κ (∆u)v = f v
and we integrate over the whole domain Ω:

Ω ∂u ∂t v -κ Ω (∆u)v = Ω f v.
We then apply the Green formula to the second term:

Ω ∂u ∂t v + κ Ω ∇u • ∇v -κ Γ (∇u • n)v = Ω f v.
Still n denotes the outward unit normal to the boundary Γ. Finally, we impose the test function v to vanish on the boundary Γ:

Ω ∂u ∂t v + κ Ω ∇u • ∇v = Ω f v.
We use the identity

Ω ∂u ∂t v = d dt Ω uv
and get the weak formulation corresponding to (34):

d dt (u(t), v) 0,Ω + κ(∇u(t), ∇v) 0,Ω = (f (t), v) 0,Ω , t > 0, u(0) = u 0 , (35)
where u(t), v ∈ H 1 0 (Ω). Remark that it integrates all the three equations present in (34) (the boundary condition has been included in H 1 0 (Ω)). We also made explicit the dependency of each term with t. Now, let us first establish a result of existence and uniqueness for this equation, as we did in the previous section. This result will allow moreover to precise the functional spaces in which the solution is defined. There are plenty of ways to establish it, for instance using a Galerkin approximation, using the semi-group theory or even using the Ladyzenskaya-Babuska-Brezzi (LBB) theory. See for instance [START_REF] Allaire | An introduction to mathematical modelling and numerical simulation[END_REF][START_REF] Ern | Theory and practice of finite elements[END_REF][START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF]. Here, we follow the path of Grégoire Allaire and state the result proven in his book. We skip the proof, that is a bit long (but not really difficult), and we recommand you study it in the book.

Theorem 8 Let us suppose that that the initial condition u 0 belongs to L 2 (Ω) and the source term f belongs to L 2 ((0, +∞); L 2 (Ω)). Then Problem (35) admits one unique solution u ∈ L 2 ((0, +∞);

H 1 0 (Ω)) ∩ C([0, +∞); L 2 (Ω)).
The above notation L 2 ((0, +∞); L 2 (Ω)) means functions square-integrable in the time-interval (0, +∞), with values in L 2 (Ω) (L 2 ((0, +∞); H 1 0 (Ω)) means the same, except that the functions take their value in H 1 0 (Ω)). The notation C([0, +∞); L 2 (Ω)) means functions of time with values in L 2 (Ω), that are continuous in time. For the detailed proof, see [1, Theorem 8.2.3] (I just adapted here for a time domain that is unbounded).

The heat equation satisfies many interesting mathematical properties, that have also a physical meaning. We will focus especially here on the following stability estimate: Proposition 2 Suppose that u 0 ∈ L 2 (Ω), and f is identically equal to zero, then the following stability estimate holds, for all t ≥ 0:

1 2 u(t) 2 0,Ω + κ t 0 ∇u(s) 2 0,Ω ds = 1 2 u 0 2 0,Ω , (36)
where u is the solution to (35).

Proof. Let us take u 0 ∈ L 2 (Ω), f = 0, and t > 0. We choose

v = u(t) ∈ H 1 0 (Ω)
in (35) and get:

d dt u(t), u(t) 0,Ω + κ (∇u(t), ∇u(t)) 0,Ω = 0.
This can be reformulated as

1 2 d dt u(t) 2 0,Ω + κ ∇u(t) 2 0,Ω . = 0.
Then we obtain (36) by integration in time.

We deduce from (36) that the solution u remains bounded in the L 2 (Ω)-norm. This norm decreases when t increases. The term in κ is a dissipation term, that depends on the physical properties of the material. There is dissipation meanwhile the temperature is not homogeneous in the material, and dissipation stops when the temperature gets uniform. We will try to mimic this behavior numerically.

Finite elements

To discretize Problem (35), we use what is called the method of lines: we first semi-discretize the space variables using finite elements, ang we will obtain a system of ODEs. Then we apply a time-marching scheme to compute the ultimate numerical approximation of these equations, as we did in the first chapter. Other strategies can be designed, for instance one could think about discretizing globally the space-time cylinder in which the equation lives. Nevertheless, the method of lines is among the simplest one and has some advantages in terms of implementation, computation cost, etc.

So we just use the Lagrange finite element space V h , as defined in the previous chapter, and built from a mesh K h of the domain Ω. The semi-discrete (FE) weak form associated to (35) consists in finding u h : t → u h (t) ∈ V h for any t > 0 that is solution to d dt (u h (t), v h ) 0,Ω + κ(∇u h (t), ∇v h ) 0,Ω = (f (t), v h ) 0,Ω , t > 0, u h (0) = u h0 , (37)

for any v h ∈ V h . The initial condition u h0 ∈ V h is an approximation of the exact initial condition u 0 , that can be obtained, for instance, using Lagrange interpolation (u h0 = I h (u 0 )). Still using the same matrix notations as in previous chapter, this problem can be recasted as :

M dU dt (t) + κKU(t) = F(t). (38)
As a result, we simply obtain a system of (linear) ODEs, and we find ourselves in the situation of the first chapter. Since K is symmetric definite positive, Cauchy-Lipschitz Theorem (for instance) ensures the existence and uniqueness of a solution U ∈ C([0; +∞); R N ). Now there just remain to apply our favorite time-marching scheme for this system of ODEs. Before this, let us state the following result, which is obtained exactly the same way as in the continuous case:

Proposition 3 Suppose that u h0 ∈ V h , and f is identically equal to zero, then the following stability estimate holds, for all t ≥ 0:

1 2 u h (t) 2 0,Ω + κ t 0 ∇u h (s) 2 0,Ω ds = 1 2 u h0 2 0,Ω , (39)
where u h is the solution to (37).

Backward Euler

Let us use the backward Euler scheme. We keep the same conventions as in the previous chapter, and, for a time-step τ > 0, and introduce a sequence (u n h ), such that, for each n :

u n h u h (t n ),
with u 0 h = u h0 , and, for n ≥ 0:

1 τ (u n+1 h -u n h , v h ) 0,Ω + κ(∇u n+1 h , ∇v h ) 0,Ω = (f (t n+1 ), v h ) 0,Ω , (40)
still with v h ∈ V h . In matrix form, we can write it as

M U n+1 -U n τ + κKU n+1 = F(t n+1 ). (41)
In this form you may (should) recognize Euler scheme as it was introduced at the very beginning of this course. Remark that Problem (40) can be recasted as

1 τ (u n+1 h , v h ) 0,Ω + κ(∇u n+1 h , ∇v h ) 0,Ω = (f (t n+1 ), v h ) 0,Ω + 1 τ (u n h , v h ) 0,Ω , (42)
so we realize that u n+1 h is solution to the discretized reaction-diffusion equation presented in the previous chapter (and you understand now the interest of the constant c, with here c = 1 τ ). In the same way, the linear system that needs to be solved at each time-step is

(M + τ κK) U n+1 = τ F(t n+1 ) + MU n . (43)
This is a system identical to what we obtained for the reaction-diffusion equation. Of course, since the matrices involved are the same at each time-step, they can be computed before we enter in the loop that computes the sequence U 1 , U 2 , . . . . We can even precompute, for instance, their Cholesky decomposition before we enter the loop. Theorems of the previous section ensure existence and uniqueness of the fully discrete solution u n h for any value of n ≥ 0. Theory about ODEs also ensures stability and accuracy of this time-marching scheme. Let us, for the end, provide another argument inspired from the stability estimates we presented previously. Proposition 4 Suppose that u 0 h ∈ V h , and f is identically equal to zero, then the following stability estimate holds, for all n ≥ 1:

where u n h is the solution to (40).

Proof. Let us take u 0 h ∈ V h , f = 0, and n ≥ 1. We choose We combine this with the first inequality and get: Then we obtain (44) by summation.

v h = u n+1 h ∈ V h in (
1 2τ u n+1 h -u n
There are some comments about this result: then it allows to check that the backward Euler scheme is unconditionnaly stable, since the L 2 -norm of u n h remains bounded by the initial condition whatever is the value of n and whatever is the value of τ . We also recover a "physical dissipation" that mimics the integral term in the estimate (36). There is also an extra dissipation term, which is purely numerical, and comes from the scheme : we say that the scheme is dissipative. It dissipates artificial energy. In some situations (let us say, for large κ and τ small enough) this dissipation is neglectible, but sometimes it is not really desired since it deteriorates the accuracy of the numerical solution.

Extra remarks

There are still many aspect to discuss about time-evolution PDEs: the general process to solve other PDEs such as the wave equation or the Schrödinger equation is the same, but care has to be taken to chose the appropriate finite element method and time discretization, so as to recover appropriately some interesting properties, such as energy conservation. Also, some PDEs such as the transport equation or some more complex first-order hyperbolic PDEs need a very specific treatment in order to recover appropriate solutions.

Proof.

  Let u ∈ H 1 0 (Ω) be the solution to Problem (25), and pick ϕ ∈ D(Ω). There holds ∆u, ϕ = -∇u, ∇ϕ by definition of the Laplace operator, in the sense of distributions (the brackets •, • denote the duality product in D (Ω)). But, since u ∈ H 1 0 (Ω), ∇u is in fact a regular distribution, which means that ∇u, ∇ϕ = Ω ∇u • ∇ϕ Now, from Problem (25), we deduce

Perspectives

This class is just an introduction to some basic facts when one needs to solve numerically some ODEs or PDEs related to physical models. Many topics have not been discussed. For instance, I did not talk about non-linear PDEs, which require some specific treatment, and specific care, to be solved appropriately. Another interesting aspect is to use some modern finite element software, such as FreeFEM++ (http://www3.freefem.org/), GetFEM++ (http: //getfem.org/) or FEniCS (https://fenicsproject.org/), that include a high-level language related to weak formulations. For instance, you can try to use FreeFEM++ here: https://freefem.org/tryit


Yu-Lin Huang
email: ylin.huang@univ-artois.fr

Gildas Morvan

Frédéric Pichon

David Mercier

SPSC: a new execution policy for exploring discrete-time stochastic simulations

Keywords: stochastic simulation, multi-agent-based simulation, solution space exploration

In this paper, we introduce a new method called SPSC (Simulation, Partitioning, Selection, Cloning) to estimate efficiently the probability of possible solutions in stochastic simulations. This method can be applied to any type of simulation, however it is particularly suitable for multi-agent-based simulations (MABS). Therefore, its performance is evaluated on a well-known MABS and compared to the classical approach, i.e., Monte Carlo.

Introduction

Multi-agent-based simulations (MABS) are widely used in various fields to study complex systems [START_REF] Railsback | Agent-Based and Individual-Based Modeling: A Practical Introduction[END_REF]. Most of them are combined with stochasticity to represent non fully controllable phenomena and use a discrete-time approach to facilitate model construction. Such model can generally be described as taking some initial conditions and some parameter set as inputs, in order to return outputs at each time step (c.f. Figure 1).

Before running into exploration of the parameter set or the initial condition space, we must first analyze outcomes from a fixed parameter set and initial conditions. Let us denote a stochastic simulation outputs (called observables in the following) at a final time step T as a random vector X T . Then a key question to address is: what is the probability P(X T ∈ S) = θ S of a specific solution S? The classical method to handle this question is Monte Carlo simulation (MC) [START_REF] Rubinstein | Simulation and the Monte Carlo method[END_REF]. It consists in simulating a number n of replications and building an estimator θS of θ S defined as:

θS = 1 n n i=1 1 S (X i T ) (1)
where 1 S is the indicator function of the set S and X i T the value of observables in the i th replication. The issue with this approach is that for the estimator to be good, the number n has generally to be large, as illustrated in Section 2.

Some methods have been developed to speed-up the computation of such simulations, such as splitting [START_REF] L'ecuyer | Rare Event Simulation using Monte Carlo Methods, chap. Splitting Techniques[END_REF] or polyagent [START_REF] Parunak | Pheromones, probabilities and multiple futures[END_REF]. However, they look for specific solutions (rare or mean), assume a particular modeling approach (Markov chains or agent-based) and require some low-level manipulations of the model.

In this paper, we propose a policy that simulates an authorized number N of replications and is as generic as the MC approach yet provides a better estimator when computational resources are limited (i.e. small N ).

The paper is organized as follows: Section 2 recalls and illustrates a standard approach to determine the required number of replications in Monte Carlo simulation for a single observable. The design principles and the approach proposed to answer the above-mentioned issues are presented in Section 3 and then applied to a classical MABS in Section 4. Section 5 concludes the paper.

Monte Carlo simulation, how many replications?

We recall in this section a standard approach to determine the number n(X T,i ) of replications to obtain a good estimator θSi of the probability P(X T,i ∈ S i ) = θ Si of some solution S i where X T,i is one observable of the vector X T . Suppose a desired relative error for the estimator θSi at confidence level 1 -α :

P( | θSi -θ Si | θ Si ≤ ) ≥ 1 -α. (2)
The minimal value for n(X T,i ) to verify (2) can be determined by applying the following algorithm [1, p. 449]:

1. Simulate n 0 replications. (n 0 observations X 1 T,i , . . . , X n0 T,i ) 2. Compute n(X T,i ) = ( Z 1-( α 2 ) • s i • X T,i ) 2 (3)
where Z 1-(α/2) is the 100(1 -α/2) quantile of the normal distribution, s i stands for the sample standard deviation over the n 0 observations and X T,i is the sample mean value over the n 0 observations. The conventional values for n 0 , and α are respectively 150, 0.05 and 0.05.

Afterward, we can then deduce the necessary number n satisfying every observable as:

n = max X T ,i ∈XT n(X T,i ) (4)
To illustrate this algorithm, let us take an academic example. We consider an environment containing vegetation and 2 types of agents: preys consuming the vegetation and predators hunting preys for food. Both preys and predators can move without restriction in the environment. This model has been implemented on the Similar platform [START_REF] Morvan | On time and consistency in multi-agent-based simulations[END_REF] and is based on the NetLogo wolf sheep predation model [START_REF] Wilensky | NetLogo wolf sheep predation model[END_REF]. The set of observables here consists of the populations of different species at each time step. The necessary number n for some arbitrarily chosen parameter set and initial state of the simulation, using the conventional values for n 0 , and α, is 3600 (c.f. Table 1). However, if we want a more precise estimation, the necessary number n of replications increases drastically: for example, considering a relative error = 0.005 yields a necessary number of replications n = 7249285.

A new execution policy for stochastic simulations

In this section, we introduce a new execution policy for stochastic simulations called SPSC (Simulation, Partitioning, Selection, Cloning). This approach relies on a decomposition of the probability of interest that we explain first.

Decomposition of the probability of interest

The probability P(X T ∈ S) concerns the observables with respect to a specific solution S at some final time step T . Thanks to the law of total probability, considering some intermediate time step j before T , we can write

P(X T ∈ S) = Sj ∈Pj P(X T ∈ S|X j ∈ S j )P(X j ∈ S j ) (5)
where P j is a partition of the state space of the random vector X j . More generally, considering all time steps before T , we can obtain the following decomposition by assuming a discrete-time system where X i depends only on X i-1 :

P(X T ∈ S) = S T -1 ∈P T -1 . . . S1∈P1 T -1 i=0 P(X i+1 ∈ S i+1 |X i ∈ S i ) (6)
where ] the following steps are applied (c.f. Figures 2 and3):

P i , i = 1, ..., T -1,
Simulation Simulate N replications from t (i) to t (i+1) , i ∈ {0, . . . , m -1},
where N corresponds to the number of replications we authorize for the simulation. Partitioning At time t (i+1) , form a partition of the space of observables of these N replications. This can be done by applying a clustering algorithm. Selection Choose one or multiple representative replications (which we call delegates) from each partition and discard the other replications. Cloning Clone the selected delegates to obtain N replications in total. Once the iterations are finished, we have created for each t (i) a partition P (i) for the state space of X t (i) . For an element S (i) ∈ P (i) of the partition at time step t (i) , after the selection and cloning steps, it has n i cloned replications. Besides, among these n i cloned replications, after evolving to the next time step t (i+1) , some of them (n i+1 replications) belong to some element S (i+1) ∈ P (i+1) . We propose to use the numbers n i and n i+1 to approximate the conditional probability P(X t (i+1) ∈ S (i+1) |X t (i) ∈ S (i) ) by:

P (X t (i+1) ∈ S (i+1) |X t (i) ∈ S (i) ) = n i+1 n i (7)
Finally, we define an estimator θS for P(X T ∈ S) using a similar decomposition as that of Equation ( 6), based on time steps t (i) and ( 7):

θS = S (m-1) ∈P (m-1)
. . .

S (1) ∈P (1) m-1 i=0 P (X t (i+1) ∈ S (i+1) |X t (i) ∈ S (i) ) (8)
where S (m) = S and S (0) = S 0 .

t(1) t(0)

Implementation

We describe here a simple implementation of SPSC used in the experiment in Section 4:

Simulation No special action is taken in this step. Partitioning N replications provide N instances of observables. To form a partition in the space of observables, we can take advantage of existing unsupervised learning algorithm which can separate instances by multiple subgroups.

The well-known clustering process kmeans has been chosen to fulfill the task. The number of cluster k is preset to 15.

Selection From any element S (i) ∈ P (i) of an intermediate partition, we select the replication which is the nearest to the center using euclidean distance on the space of observables. Cloning After partitioning and selection, k delegates are obtained to be cloned.

To come back to N replications in total, we clone each delegate N k times. If k does not divide N , we select randomly the remainder number N -k * N k of delegates to produce one more clone per selected delegate.

The time interval [0, T ] is homogeneously split into m = 5 pieces (i.e. ∀i ∈ {0, 1, . . . , 5}, t (i) = i × T 5 ).

Experiment

Let us take the prey predator model mentioned previously in Section 2 as an example. As this model is well-known and well-studied, we can give some possible solutions before launching simulations:

S 1 : Extinction of preys and predators, only vegetation remains. S 2 : Predators go extinct, preys live without nature enemy's harass. S 3 : All species survive and form a stable ecosystem. Now the question is, for an arbitrary parameter set and initial condition, what is the probability of these solutions at a given time step (e.g. T = 1000)? To answer this question, the MC approach recalled in Section 2 is generally used.

In the following, we compare the performances of MC and SPSC. The validation is done by comparing the outputs of these methods with the same limited number of replications N = 50. By repeating the simulations 1000 times, we will be able to compare statistically the results obtained by MC and SPSC. Two performance measures are considered here : 1) The detection rate of a specific solution S. 2) The precision of the probability estimator for a specific solution. Before evaluating these performance measures, we have done 30000 replications using MC in order to provide reference values for the comparisons:

P ref (S 1 ) ≈ 0.0065, P ref (S 2 ) ≈ 0.0126, P ref (S 3 ) ≈ 0.981 (9)
The detection rates obtained with MC and SPSC policies, i.e, the capacity of identifying a specific solution, from N = 50 replications are summed up in Table 2. The first three columns indicate the detection rate of single solutions and the last column indicate the detection rate for the three solutions simultaneously. We can then deduce that SPSC explores more efficiently the solution space.

To evaluate the precision of the probability estimator, the absolute error between the probability estimator outcomes and the references is computed:

Err(S i ) = | P (S i ) -P ref (S i )|. ( 10
)
Furthermore, to gain an entire vision on the three solutions simultaneously, we consider also the mean of three solutions relative errors: We can notice that the distribution of errors are not normal nor symmetric. Thus, to compare the errors from MC and SPSC policies, the Wilcoxon-Mann-Whitney test is applied with the threshold α = 0.05. The test results are presented in Table 3 with p-value and alternative hypotheses, we can then conclude that SPSC yields better probability estimates for each solution than MC.

Err = 1≤i≤3 | P (S i ) -P ref (S i ) P ref (S i ) |. (11)

Conclusions and perspectives

We have introduced a generic policy called SPSC for executing stochastic simulations that deals with the weakness of MC when the number of replications is limited. It treats simulations as black boxes and therefore, does not rely upon a priori knowledge. We have also presented a simple implementation of SPSC and run it on a classic stochastic MABS model. By comparing the results obtained with SPSC and with MC, we can conclude that SPSC gives a better solution probability estimation and can reveal more different solutions than MC.

The first perspectives of this work are related to the impact of the parameters (N , k, etc.), the partitioning algorithm as well as the selection and cloning strategies on the performance. For instance, instead of having the same number of clones for each delegate, we could clone more the delegates from small partitions and less the delegates from large partitions. Theoretical properties of the proposed solution as well as its interest for multimodal transport simulation will also be investigated.

Moreover, since we deal with small sample size at each intermediate time step, we could take advantage of modern tools [START_REF] Kanjanatarakul | Prediction of future observations using belief functions: A likelihood-based approach[END_REF] for statistical inference to compute the estimator θS .

Fig. 1 .

 1 Fig. 1. Illustration of a discrete-time stochastic simulation

Fig. 2 .

 2 Fig. 2. SPSC process diagram

Fig. 3 .

 3 Fig. 3. Illustration of the first and second iterations of SPSC, starting from a single initial state.

Fig. 4 .

 4 Fig. 4. Comparison of errors

Table 1 .

 1 Determination of the necessary number of replications for the prey predator model implemented on the Similar platform. The parameters applied are n0 = 150, = 0.05 and α = 0.05.

	XT,i	si	XT,i	n(XT,i)
	number of preys	697.83 783.77	1219
	number of predators 196.95 128.67	3600

  is a partition of the state space of X i , S T = S and S 0 is the initial state of the simulation.

	3.2 SPSC: Simulation, Partitioning, Selection, Cloning

Inspired by the decomposition (6), we split the time interval [0, T ] into m pieces: [t (0) , t (1) ], [t (1) , t (2) ], ..., [t (m-1) , t (m) ] where t (0) = 0 < t (1) < ... < t (m-1) < t (m) = T . Then, for each interval [t (i) , t (i+1)

Table 2 .

 2 Detection rates obtained with MC and SPSC when launching 50 replications.Histograms of these errors for the policies SPSC and MC are shown in Fig.4.

								S1		S2		S3 S1, S2 and S3
					MC 0.236 0.455	1			0.102
					SPSC 0.332 0.617	1			0.203
		1000	S1 absolute errors	1000	S2 absolute errors	1000	S3 absolute errors	1000		Mean relative errors
		800					800					800					800
	SPSC	400 600					400 600					400 600					400 600
		200					200					200					200
		0.00 0	0.02	0.04	0.06	0.08	0.00 0	0.02	0.04	0.06	0.08	0.00 0	0.02	0.04	0.06	0.08	0	0	1	2	3	4
		1000				1000
		800					800
	MC	400 600					400 600
		200					200
		0.00 0	0.02	0.04	0.06	0.08	0.00 0	0.02	0.04	0.06	0.08	0.00 0	0.02	0.04	0.06	0.08

Table 3 .

 3 Hypothesis and p-value given by Wilcoxon-Mann-Whitney test.

	Target solution Alternative hypotheses	p-value	Conclusion
	S1	ErrSP SC < ErrMC	2.2e-16 ErrSP SC < ErrMC
	S2	ErrSP SC < ErrMC 3.163e-16 ErrSP SC < ErrMC
	S3	ErrSP SC < ErrMC 9.849e-09 ErrSP SC < ErrMC
	S1, S2 and S3	ErrSP SC < ErrMC	2.2e-16 ErrSP SC < ErrMC

Acknowledgement

The ELSAT2020 project is co-financed by the European Union with the European Regional Development Fund, the French state and the Hauts de France Region Council.