Concentration Phenomenon in Some Non-Local Equation

We are interested in the long time behaviour of the positive solutions of the Cauchy problem involving the following integro-differential equation $$\partial\_t u(t, x) = \left(a(x) -- \int\_{\Omega} k(x, y)u(t, y) dy\right ) u(t, x) + \int\_{\Omega} m(x, y)[u(t, y) -- u(t, x)] dy\quad \text{ for}\quad (t, x) $\in$ \mathbb{R}\_{+} \times \Omega,$$ together with the initial condition $u(0, \cdot) = u0 \quad \text{ in }\quad \Omega$. Such a problem is used in population dynamics models to capture the evolution of a clonal population structured with respect to a phenotypic trait. In this context, the function u represents the density of individuals characterized by the trait, the domain of trait values $\Omega$ is a bounded subset of $\mathbb{R}^N$ , the kernels $k$ and $m$ respectively account for the competition between individuals and the mutations occurring in every generation, and the function a represents a growth rate. When the competition is independent of the trait, we construct a positive stationary solution which belongs to the space of Radon measures on $\Omega$. Moreover, when this '' stationary '' measure is regular and bounded, we prove its uniqueness and show that, for any non negative initial datum in $L^{\infty} (\Omega) \cap L^1 (\Omega)$, the solution of the Cauchy problem converges to this limit measure in $L^2 (\Omega)$. We also construct an example for which the measure is singular and non-unique, and investigate numerically the long time behaviour of the solution in such a situation. These numerical simulations seem to reveal some dependence of the limit measure with respect to the initial datum.

1. Introduction. In this paper, we are interested in the evolution of a clonal population structured with respect to a phenotypic trait and essentially subjected to three processes: mutation, growth, and competition. As an example, one can think of a virus population structured by its virulence, as this trait can be easily quantified from experimental data. For such type of population, a mathematical model commonly used (see for instance [11,10,22,24,12,13,31,15,14,32,33]) is based on the following partial integro-differential equations: ∀x ∈ Ω, u(0, x) = u 0 (x), (2) in which the non-negative function u is the density of individuals of the considered population characterized by the trait x, the functions k and a are respectively a kernel accounting for interaction between individuals through competition and a growth rate, the set Ω is a bounded domain of R N , and M denotes a linear diffusion operator modelling the mutation process. Depending on the context, several kinds of diffusion operators have been considered in the literature, see [11,10,24,12,3,15,21,27,32,30,28] among others. In the present work, our attention is focused on models for which the operator M is integral and of the form the function m being a positive kernel satisfying some integrability conditions. Lately, this type of integro-differential equation has attracted a lot of interest and much effort has been made to analyse the solutions of (1). In particular, existence of global solutions in C 1 (R + , L 1 (Ω) ∩ L ∞ (Ω)) for any non-negative initial datum in L ∞ (Ω) and fairly general assumptions on the functions k, m, a and the domain Ω was established in [12,13,21,14], and the local stability of bounded continuous stationary solutions in unidimensional domains Ω ⊂ R was investigated in [32,14,33]. However, the analysis of stationary solutions of (1) in higher dimension remains to be done, whereas the long time behaviour of positive solutions of problem (1)-(2) is still not fully understood.
When mutations are neglected (that is, when m ≡ 0), equation (1) is reduced to and, for any generic positive initial datum u 0 , the solution to (4)-(2) is known to converge weakly to a positive Radon measure µ [23,21,26]. This measure is, in some sense, a stationary solution of (4) representing an evolutionarily stable strategy for the system. For example, when the kernel k is positive and does not depend on the trait, the support of µ lies in the set A := {x ∈ Ω | ∀y ∈ Ω, a(y) ≤ a(x)}, and one may check in such a situation that the sum of Dirac masses is indeed a stationary solution of the equation. Moreover, when the measure µ is unique, any positive solution u(t, ·) of (4)-(2) converges weakly to it as t tends to infinity (see [26] for a proof).
Since the mutation process can be seen as a diffusion operator on the trait space, it is expected that the long time behaviour of a positive solution to problem (1)-(2) is simple and that concentration phenomena do not occur. This is verified when the mutation operator M is a classical elliptic operator [27,17]. For an integral operator, as in the present case, the existence of bounded equilibria when the domain Ω is unidimensional seems to give credit to such a conjecture. However, we prove that it is false in higher dimension, by exhibiting a class of situations for which a positive singular measure µ, solution of (1), can be constructed and by investigating numerically the long time behaviour of positive solutions of the corresponding Cauchy problem.
1.1. Main results. Let us state precisely the assumptions on the domain Ω, the kernels k and m and the function a under which the results are obtained. First, we suppose that the domain Ω is an open bounded connected set of R N with a Lipschitz boundary, that the function a is such that a is continuous on Ω and positive, (5) and that the function m is a non-negative symmetric Carathéodory kernel function, that is, Finally, we assume that the kernel k is independent of the trait, that is, for any (x, y) in Ω × Ω, k(x, y) = k(y), and that it satisfies the following condition: there exist positive constants c 0 and C 0 such that c 0 1 Ω ≤ k ≤ C 0 1 Ω , where 1 Ω denotes the characteristic function of the set Ω.
We start by considering stationary solutions of (1), that is, satisfying the equation Under the above assumptions, we prove that there exists a positive Radon measure solving (8) in a weak sense.
Theorem 1.1. Assume that the functions a, m and k respectively satisfy conditions (5) to (7). Then, there exists a positive Radon measure µ such that ∀ϕ ∈ C c (Ω), Let λ p be the principal eigenvalue of the operator M + a defined by We have the following characterisation for the measure µ.
• If λ p is associated with an eigenfunction ϕ p which belongs to L 1 (Ω), then µ is an absolutely continuous measure, that is, there existsū in L 1 (Ω) such that dµ(x) =ū(x)dx, and the unique strong solution to equation (8). Moreover, the functionū is in L ∞ (Ω) when the principal eigenfunction ϕ p belongs to L ∞ (Ω). • Otherwise, µ is a singular measure.
As a consequence of this theorem, the existence of a singular measure solution to (9) is strongly related to the non-existence of an integrable principal eigenfunction for the non-local operator M + a, a result which was proved recently [16,35,18].
Next, the global stability of the measure µ and the long time behaviour of positive solutions of the problem (1)-(2) are analysed. When the measure is regular, we have the following convergence result. Theorem 1.2. Assume that the functions a, m and k respectively satisfy conditions (5) to (7) and that there exists an absolutely continuous positive Radon measure µ, such that dµ(x) =ū(x) dx, solution to (8). Assume further that the functionū belongs to L ∞ (Ω). Then, for any non-negative initial datum u 0 in L 1 (Ω) ∩ L ∞ (Ω), the positive solution u of problem (1)-(2) satisfies lim t→∞ u(t, ·) −ū L 2 (Ω) = 0.
Note that this last global stability result implies the uniqueness of the measure solution of (8). On the contrary, when no regular positive Radon measure exists, the convergence of a positive solution of (1) to a stationary one is delicate to analyse. To shed light on the possible dynamics prevailing in such a situation, we have conducted a few numerical experiments.

Numerical simulations.
In order to gain some insight on the long time behaviour of solutions of problem (1)-(2), we compute approximate solutions, for different choices of growth function a and initial datum u 0 , using a numerical method. Limiting ourselves to preliminary computations, the retained configuration is such that the two-dimensional domain Ω = B 1/4 (0) the open ball of radius 0.25 centred at the origin, and the competition and mutation kernels are two constant functions, such that k ≡ 1 and m ≡ ρ, with ρ a positive constant, respectively. The problem to be solved thus reduces to the equation ∂ t u(t, x) = ρˆΩ(u(t, y) − u(t, x)) dy + a(x) −ˆΩ u(t, y) dy u(t, x), for any (t, x) in R * + × Ω, completed by the initial condition ∀x ∈ Ω, u(0, x) = u 0 (x).
1.2.1. A simple growth rate. We first consider a configuration in which the growth rate achieves its maximum at a single point, a case for which the uniqueness of the stationary solution can be shown.
Proposition 1. For any positive value ρ, there exist a unique positive measure µ which is a stationary solution of (10), and a critical value ρ * such that this measure is singular for ρ < ρ * or regular otherwise. In addition, for any non-negative initial datum u 0 in L 1 (Ω) ∩ L ∞ (Ω), the solution of problem (10)-(11) converges weakly to µ.
This proposition is a direct consequence of Theorem 1.1 and of the uniform L 1 estimates obtained in Section 3. To illustrate its conclusions, we take a(x) = 1 − x 2 , where · 2 denotes the Euclidean norm in R 2 , and solve numerically the problem. The obtained results, presented in Figure 1, provide a clear picture of the dynamics of the solution with respect to the constant value of the mutation rate. configurations, in which only the mutation rate differs. The competition rate is constant and set to 1 and the growth rate function achieves its maximum only at the origin while the initial datum u 0 is uniform with value 1. We have set ρ = 1 for the first simulation (subfigures (A) to (D)), and ρ = 0.1 for the second one (subfigures (E) to (H)). In both situations, we observe the convergence to a stationary solution, either to a regular measure (see subfigure (D)) or to a singular measure with one Dirac mass at the origin (see subfigure (H)), the latter being characteristic of a concentration phenomenon. In the regular case (subfigures (A) to (D)), the stationarity being attained numerically around t = 590 times for two configurations, which differ only in their initial datum. The mutation and competition rates are constant and set respectively to 2 and 1, the growth rate function achieves its maximum at four points and the initial datum u 0 is such that it vanishes on three (subfigures (A) to (D)) or two (subfigures (E) to (H)) of these points. In both cases, rapid convergence of the approximate solution towards an identical regular stationary state is observed, the numerical stationarity being attained around t = 85. The mutation and competition rates are constant and set respectively to 0.01 and 1, the growth rate function achieves its maximum at four points and the initial datum u 0 is such that it vanishes on three of these four points. We observe a slow convergence of the numerical solution towards the approximation of a singular stationary measure containing a single Dirac mass. The approximate solution continues to take large increasing values in a single element at t = 1000).

1.2.2.
A complex growth rate. We next explore a situation for which the growth rate reaches its maximum at multiple points. In such a setting, we expect the stationary measure to be non-unique. In order to numerically investigate this conjecture, we consider a growth rate function of the form . With this choice, for a sufficiently small constant value ρ of the mutation rate, we can show that there are at least four different positive Radon measures which are stationary solutions of the problem. The impact of the non-uniqueness can be seen in the simulations presented in Figures 3 and 4. Indeed, as soon as the regime allows several singular stationary measures to exist, we observe that the outcome of the simulation may drastically differ depending on the initial datum. In contrast, when the mutation rate is such that the stationary measure is regular, the stationary state is a global attractor (see Figure 2).  11) at different times. The mutation and competition rates are constant and set respectively to 0.01 and 1, the growth rate function achieves its maximum at four points and the initial datum u 0 is such that it vanishes on two of these four points. We observe a slow convergence of the numerical solution towards the approximation of a singular stationary measure containing two Dirac masses.
1.3. Outline. The paper is organised as follows. We start by recalling important facts about the spectral properties of the class of non-local operators considered in Section 2. Uniform estimates are derived by means of non-linear relative entropy formulas in Section 3, leading to proofs of Theorems 1.1 and 1.2 in Section 4. Finally, the numerical method used for the simulations is briefly described in an appendix.

Spectral properties of non-local operators. Consider the eigenproblem
in which M is the integral operator defined by (3), whose kernel satisfies condition (6). When the function a is not constant, neither the operator M + a + λ nor its inverse are compact, and the Krein-Rutman theory fails in providing existence of the principal eigenvalue of M + a. However, a variational formula, introduced in [7] to characterise the principal eigenvalue of elliptic operators, can be transposed to the operator M + a. Namely, the quantity (13) is well defined and called the generalised principal eigenvalue of M + a. It is also well-known (see [16,35,20,4,5,34]) that λ p (M + a) is not always an eigenvalue of M + a in a reasonable Banach space, which means there is not always a positive continuous eigenfunction associated with it. Nevertheless, as shown in [18], there always exists an associated positive Radon measure.  (3), whose kernel m satisfies condition (6), and define Then, there exists a positive Radon measure µ p such that, for any function ϕ in In addition, we have the following dichotomy: • or there exists g p in C(Ω), g p > 0, and ν a positive singular measure with respect to the Lebesgue measure, whose support lies in the set A, such that The measure µ p can be characterised more precisely and there exists a simple criterion guaranteeing its regularity. 16,20,5,19]). Under the assumptions of the preceding theorem, one has dµ p (x) = ϕ p (x) dx with ϕ p in C(Ω), ϕ p > 0, if and only if λ p (M+a) < −σ.
We conclude with an alternate characterisation of λ p (M + a). In the spirit of the works in [8,6,9], let us define the quantity As in the case of elliptic operators, we have the following result. 20,5,19]). Let Ω ⊂ R N be a bounded domain, a be a continuous function over Ω, M be the integral operator defined by (3), whose kernel m satisfies condition (6). Then, we have

3.
A priori estimates. In this section, given a non-negative initial datum u 0 in L 1 (Ω) ∩ L ∞ (Ω), we establish some uniform in time a priori estimates on the solution of problem (1)- (2). To do so, we start by proving a non-linear relative entropy identity satisfied by any solution of (1).

Proposition 3 (general identity).
Let Ω ⊂ R N be a bounded domain and assume that the functions a, m and k respectively satisfy conditions (5) to (7). Let H be a smooth (of class C 1 at least) function. Finally, let the functionū in L 1 (Ω) ∩ L ∞ (Ω) be a positive stationary solution of (1), and the function u in C 1 ((0, +∞), L ∞ (Ω)) be a solution of (1). Then we have where H H,ū [u], D H,ū (u) and Γ are respectively defined by where ν is the measure on Ω × Ω defined by ν := m(x, y)ū(x)ū(y) dx dy.
proof. Since the kernel k satisfies condition (7) and defining Γ as above, we have from equation (1) that, for all t in R * + and almost every x in Ω, The functionū being a positive stationary solution of (1), we also have that, for almost every x in x in Ω, and we can rewrite the above equation as follows: Multiplying this identity by the functionū(x)H u(t,x) and integrating over Ω, we find that, for all t in R * + , Using definition (3) and rearranging the terms, we get, for all t in R * + , where ν stands for the positive measure ν := m(x, y)ū(x)ū(y)dxdy.
Due to the symmetry property of the kernel m, we straightforwardly see that and, by combining the above equalities and using the remaining definitions, we finally reach Remark 1. When k ≡ 0, equation (1) is linear and relative entropy formulas are well-known in this case [29].
Remark 2. When the function H is non-decreasing andū is only assumed to be a stationary super-solution of (1), we clearly see from the above proof that Similarly, ifū is a positive stationary sub-solution of (1), we have Equipped with this general relative entropy identity, we are in a position to derive useful differential inequalities.

Proposition 4.
Let Ω ⊂ R N be a bounded domain and assume that the functions a, k and m respectively satisfy conditions (5) to (7). Let the functionsū and u be defined as in Proposition 3, q in [1, +∞) and H q be the smooth convex function Moreover, we have for all t ∈ R * + , where ν is the positive measure in Ω × Ω defined by ν := m(x, y)ū(x)ū(y)dxdy.
Remark 3. When q = 2, we observe that H H2,ū [u](t) = u(t, ·) 2 L 2 (Ω) , yielding a Lyapunov functional which involves the L 2 (Ω)-norm of u(t, ·) instead of a weighted L q (Ω)-norm. Indeed, we have, for all t in R * + , d dt log Proof of Proposition 4. Owing to Proposition 3, we have, for H = H q ,  (18) and, combining these two relations, we end up with Equality (16) then follows from direct computations, by using symmetry and an obvious change of variables.
We may now obtain uniform in time a priori bounds on the L 1 (Ω)-norm of a solution to problem (1)-(2). Lemma 3.1. Let Ω ⊂ R N be a bounded domain and assume that the functions a, m and k respectively satisfy conditions (5) to (7). Let the function u in C 1 (R * + , L 1 (Ω)∩ L ∞ (Ω)) be a non-negative solution of the Cauchy problem (1)-(2) with non-negative initial datum u 0 in L 1 (Ω) ∩ L ∞ (Ω), such that u 0 L 1 (Ω) = 0. Then, there exist two positive constants c 1 and C 1 , depending on u 0 , such that Proof. We first observe that large (respectively small) constants are stationary super-solutions (respectively sub-solutions) of equation (1). Indeed, given a constant Therefore, from Proposition 3 and Remark 2, replacing the stationary solutionū by a large, respectively small, constant and considering the convex function H = H 1 , we obtain, for all t in R * + and C ≥ sup x∈Ω a(x) respectively, for all t in R * + and c ≤ inf x∈Ω a(x) Using the fact that the kernel k satisfies condition (7), we finally get, for all t in Let us c 1 and C 1 be the following positive constants, From the logistic character of these differential inequalities and since u 0 L 1 (Ω) > 0, we deduce that, for all t in R + , Remark 4. One may observe that the above proof relies on the fact that the function k is bounded above and below by positive constants. As a consequence, such an uniform L 1 (Ω)-norm estimate holds more generally if the kernel k depends on the trait x.

Proofs.
We are now in a position to prove Theorems 1.1 and 1.2.

Construction of a stationary state.
To construct a solution to stationary equation (8) in the space of Radon measures (and thus prove Theorem 1.1), we look for a measure µ satisfying (9). Owing to Theorem 2.1, consider a positive measure µ p associated with λ p (M + a) and normalised in order to have´Ω dµ p (x) = 1. We are reduced to proving the following assertion.
Claim 4.1. There exists a unique positive real number θ such that θ µ p satisfies (9).
. We have, for all ϕ in C c (Ω), so that θ µ p is a solution to (9). To conclude, it remains to show that −λ p (M+a) > 0. We first notice that the couple λ = − inf x∈Ω a(x) and ϕ ≡ 1 satisfies It then follows directly from the alternate characterisation of the principal eigenvalue, given by Theorem 2.2, that Remark 5. From the above computation, it is clearly seen that the uniqueness of the stationary state follows from the uniqueness of the measure associated with the principal eigenvalue λ p (M + a).

4.2.
Long time behaviour. This final subsection is devoted to the proof of Theorem 1.2. We now assume that the positive measure µ constructed above is absolutely continuous and bounded, that is, there exists a functionū in L 1 (Ω) ∩ L ∞ (Ω) such that dµ(x) =ū(x) dx. This measure being linked to the principal eigenvalue λ p (M + a) through Claim 4.1, we have in addition thatū = θ ϕ p , with ϕ p a function in L 1 (Ω) ∩ L ∞ (Ω). From the regularity of ϕ p , we infer that the function u is a strong solution of equation (8). Thus, knowing that a positive continuous stationary solution exists, we may derive further a priori estimates on the solution u of problem (1)-(2).
On the other hand, the uniform upper bound is a straightforward application of Proposition 4. More precisely, the functionū in L 1 (Ω) ∩ L ∞ (Ω) being a positive stationary solution of (1), the quantity F(t) = log is a decreasing function of t. Therefore, we have and, using Lemma 3.1, we finally reach To show that the solution u converges to a stationary state, we introduce the decomposition ∀t ∈ R * + , u(t, ·) = λ(t)ū + h(t, ·), (19) in which the function h is such that Proof. Let us denote by ·, · L 2 the standard inner product over L 2 (Ω). It stems from the decomposition of the solution u that ∀t ∈ R + , u(t, ·),ū L 2 = θ u(t, ·), ϕ p L 2 = λ(t)θ 2 ϕ p , ϕ p L 2 .
Since the function ϕ p is positive and bounded over Ω, we have from Lemma 3.1 which gives We next observe that an obvious upper bound for h(t, ·) L 2 (Ω) can be derived from Lemma 4.2 since, by construction, one has Substituting to the solution u its decomposition into equation (1), we moreover have As H H2,ū [h](t) = h(t, ·) 2 L 2 (Ω) by definition, we obtain, by multiplying the above equation by h and integrating over the set Ω, where we have used the definition ofū and the fact that ū, h(t, ·) L 2 (Ω) = 0 at any time t. By proceeding as in the proof of Proposition 3 with the choice H = H 2 , we see that with Γ(t) = ˆΩ k(y)ū(y) dy −ˆΩ k(y)u(t, y) dy = −λ p −ˆΩ k(y)u(t, y) dy .
The function H H2,ū [h] being non-negative, we either have that it is positive at all times or there exists t 0 in R + such that H H2,ū [h](t 0 ) = 0. In the latter case, this means that u(t 0 , x) = λ(t 0 )θϕ p (x) for almost every x in Ω. Setting w(t, ·) := η(t)θϕ p , with the function η solution of the initial value problem we can check that w is a solution of (1) for all times larger than or equal to t 0 . Since w(t 0 , ·) = u(t 0 , ·), we have, by uniqueness of the solution of the Cauchy problem (1), that u(t, · = w(t, ·) for all t ≥ t 0 , which yields that h(t, ·) = 0 and λ(t) = η(t) for all t ≥ t 0 and therefore lim t→+∞ λ(t) = lim t→+∞ η(t) = 1. In the former case, we resort to the following result. Indeed, let us assume for the moment that Claim 4.4 holds. It follows from (19) and we deduce from Proposition 4 that The Claim implying that h(t, ·) 2 L 2 (Ω) = H H2,ū [h](t) vanishes as t tends to infinity, we then have that , and an elementary analysis of this ordinary differential equation finally shows that λ(t) tends to 1 as t tends to infinity.
Let us now conclude by giving the postponed proof of the above Claim. Proof of Claim 4.4: The function H H2,ū [h] being positive, we see, using (22) and following the proof of Proposition 4, that the function F(t) = log and is thus both smooth and non-increasing. To prove the Claim, it is thus sufficient, due to this monotonicity property, to exhibit a sequence (t n ) n∈N tending to infinity such that the corresponding sequence (H H2,ū [h](t n )) n∈N has a vanishing limit, which amounts to proving that inf t∈R+ H H2,ū [h](t) = 0. Now, let us assume that inf t∈R+ H H2,ū [h](t) = κ > 0. Then, using (21) and (20) respectively, we infer there exist three positive constants α, β and γ such that ∀t ∈ R * + , 0 < κ ≤ H H2,ū [h](t) ≤ α and 0 < β ≤ H H1,ū [u](t) ≤ γ. As a consequence from the above, there exists a constant c 0 such that lim t→+∞ F(t) = c 0 and lim t→+∞ d dt Let (t n ) n∈N be a sequence that tends to infinity. On the one hand, it follows from (23) and (24) that lim On the other hand, the sequence (h(t n , ·)) n∈N being bounded in L 2 (Ω), it has a weakly convergent subsequence in L 2 (Ω), which we still denote (h(t n , ·)) n∈N , with limith. We have where the function g(x) :=ˆΩ m(x, y)ū (y) u(x) dy is well defined over Ω and bounded from below by a positive constant C, the functionsū and m being both positive and bounded. By Fatou's lemma and the weak convergence of (h(t n , ·)) n∈N in L 2 (Ω), we thus getˆΩ and lim n→+∞ˆΩˆΩ m(x, y)h(t n , x)h(t n , y) dx dy =ˆΩˆΩ m(x, y)h(x)h(y) dx dy.
Therefore, one has which enforces thath = τū for some real constant τ . Recalling that that the function h(t n , ·) is orthogonal toū for any integer n, we then find that implying that τ = 0. Using (25) and (26) with the fact thath ≡ 0 thus leads to the following contradiction Hence, one has inf t∈R+ H H2,ū [h](t) = 0, which implies that lim inf t→+∞ H H2,ū [h](t) = 0, the function being non-negative. The proof is complete.
Appendix A. The numerical method. To investigate numerically the asymptotic behaviour of solutions to problem (1)-(2), we assume that we are in the setting of Lemma 4.2, so that solutions belong to C 1 ((0, T ], L 2 (Ω)) for all positive times T . We use a spatial discretisation based on the finite element method and solve the resulting ordinary differential equation with an implicit-explicit (IMEX) timediscretisation scheme. In this way, the linear terms in the equation are treated implicitly in time, while the non-linear reaction term is dealt with explicitly in time.
The first step in deriving the method is to write the problem in variational form. To this end, we multiply equation (1) by a test function w in L 2 (Ω) and integrate over the domain Ω to obtain ∀t ∈ R * + , ∀w ∈ L 2 (Ω), m(x, y)u(t, y) dy w(x) dx+ Ω a(x) −ˆΩ m(x, y) dy u(t, x)w(x)dx−ˆΩ ˆΩ k(x, y)u(t, y) dy u(t, x)w(x)dx.
Next, given a mesh T h of the domain Ω, that is a collection of geometrically simple elements (we used triangles in the simulations) partitioning 1 the set Ω, the positive real parameter h being related to the size of the elements in T h , an approximation space V h , which is a finite dimensional subspace of L 2 (Ω) composed of functions which are constant on each element K of the partition T h , is constructed. Setting N = dim V h and considering a basis {w i } i=1,...,N of V h , the restriction of the above variational problem to the set V h gives rise to a system of N ordinary differential equations, of the form where the components of the vector u h are the coefficients at time t of the semidiscrete approximate solution u h (t) in the basis of V h , i.e. u h (t, x) = The initial value function of this semi-discrete problem is the projection of the function u 0 onto the chosen finite element space. While the mass matrices are sparse, that is most of their entries are zero, the matrix K(u) associated with the non-linear term is a priori dense, which prevents the use of fully implicit time-integration schemes to solve the previous differential system when the integer N is large. A common strategy is thus to consider a combination of implicit and explicit (linear multistep or Runge-Kutta) schemes (see [1,2] for instance), the most simple case of being that of the forward and backward Euler methods, yielding the system of algebraic equations where ∆t is the length of the fixed time step (which has to be chosen small enough due to stability restrictions) and u (n) h is the approximation of u h at time n ∆t, n ∈ 1, . . . , T ∆t , to be solved at each iteration.
This method was implemented using FreeFem++ [25] to obtain the numerical simulations presented in Subsection 1.2. The same mesh was used in all the simulations and composed of approximately 30000 triangles. The length of the time step was chosen constant and equal to 0.1. Note that the piecewise constant numerical solutions were interpolated continuously for the visual representations in Figures 1  to 2.