A structure-preserving discontinuous Galerkin scheme for the Fischer-KPP equation

An implicit Euler discontinuous Galerkin scheme for the Fisher-Kolmogorov-Petrovsky-Piscounov (Fisher-KPP) equation for population densities with no-flux boundary conditions is suggested and analyzed. Using an exponential variable transformation, the numerical scheme automatically preserves the positivity of the discrete solution. A discrete entropy inequality is derived, and the exponential time decay of the discrete density to the stable steady state in the L1 norm is proved if the initial entropy is smaller than the measure of the domain. The discrete solution is proved to converge in the L2 norm to the unique strong solution to the time-discrete Fisher-KPP equation as the mesh size tends to zero. Numerical experiments in one space dimension illustrate the theoretical results.


Introduction
The preservation of the structure of nonlinear diffusion equations on the discrete level is of paramount importance in applications. While there has been an enormous progress on structure-preserving schemes for ordinary differential equations (see, e.g., [16]), the development of structure-preserving numerical techniques for nonlinear diffusion equations is still an ongoing quest, in particular for higher-order methods. In this paper, we analyze a toy problem, the Fisher-Kolmogorov-Petrovsky-Piscounov (Fisher-KPP) equation with no-flux boundary conditions, to devise an implicit Euler discontinuous Galerkin scheme which preserves the positivity of the solution, the entropy structure, and the exponential equilibration on the discrete level. In a future work, we aim to extend the scheme to diffusion systems.
The Fisher-KPP equation [12] is the reaction-diffusion equation ∇u¨n " 0 on BΩ, up0q " u 0 in Ω, (2) where D ą 0 is the diffusion coefficient, Ω Ă R d a bounded domain, and n the exterior unit normal vector on the boundary BΩ. The variable upx, tq models a population density or chemical concentration, influenced by diffusion and logistic growth. The Fisher-KPP equation admits traveling-wave solutions upx, tq " φpx´ctq, which switch between the unstable steady state u˚" 0 and the stable steady state u˚" 1. By the maxiumum principle, the density stays nonnegative if it does so initially, and it satisfies the entropy inequality (3) d dt If there are no reaction terms, we have conservation of the total mass, and the logarithmic Sobolev inequality implies the exponential decay of the (mathematical) entropy Sptq " ş Ω puptqplog uptq´1q`1qdx (see, e.g., [20,Chapter 2]). When reaction terms are present, the situation is more delicate, since there are two steady states, u˚" 0 and u˚" 1. If the initial entropy Sp0q is smaller than the measure of Ω, then uptq converges exponentially fast to u˚" 1 in the L 1 pΩq norm. Our objective is to preserve the aforementioned properties on the discrete level.
It is well known that the preservation of the positivity or nonnegativity of discrete solutions for (1) may fail in standard (finite-element) schemes, in particular when the solution vanishes in some region; see Section 5 for an example. Our key idea to preserve the positivity is to employ the exponential transformation u " e λ . Such a transformation or a variant is used, for instance, in the Il'in scheme [19] and in the existence analysis of drift-diffusion equations [13]. Moreover, it allows for the preservation of L 8 pΩq bounds in volume-filling cross-diffusion systems [8,20]. The implicit Euler scheme for (1)- (2) in the exponential variable then reads as 1 △t`e λ k´e λ k´1˘" divpe λ k ∇e λ k q`e λ k p1´e λ k q in Ω, (4) ∇λ k¨n " 0 on BΩ, (5) where here and in the following, we set D " 1 for simplicity and we choose 0 ă △t ă 1. At first glance, one may think that this formulation unnecessarily complicates the problem, but we will show that it enjoys some useful properties.
We propose a discontinuous Galerkin (DG) discretization for problem (4)-(5) with variable λ k h , where h ą 0 is the maximal diameter of the mesh elements. The nonlinear diffusion term is discretized by an interior penalty DG method. By construction, the discrete densities exppλ k h q are positive, and the scheme also preserves the entropy structure and large-time asymptotics. Our main results can be sketched as follows: ‚ Existence of a solution λ k h to the implicit Euler DG scheme (8), given a function λ k´1 h (Proposition 6). This result is based on the Leray-Schauder fixed-point theorem and a coercivity estimate. ‚ Discrete entropy inequality (Lemma 7). The inequality follows from scheme (8) using the test function λ k h and the convexity of u Þ Ñ uplog u´1q`1.
‚ Exponential decay of the discrete entropy S k h :" ż Ω`e λ k h pe λ k h´1 q`1˘dx ď S 0 h e´κ k△t (Proposition 9) and of the L 1 norm of e λ k h´1 (Theorem 11). The result holds if S 0 h ă |Ω|. This condition implies a positive lower bound for the total mass ş Ω e λ k h dx, which is needed to guarantee that the discrete solution converges to the stable steady state u˚" 1 and not to the steady state u˚" 0. The case S 0 h ě |Ω| is discussed in Remark 12. ‚ Convergence of the scheme (Theorem 14): There exists a unique strong solution u k P H 2 n pΩq to the implicit Euler discretization associated to (1)- (2) such that e λ k h Ñ u k strongly in L 2 pΩq as h Ñ 0. The result is based on a compactness property, which is a consequence of the gradient estimate from the entropy inequality and a coercivity estimate. This yields a very weak semi-discrete solution, which turns out to be a strong solution thanks to a duality argument. Let us put our results into context and review the state of the art of stucture preservation in DG methods. The DG scheme was introduced in the early 1970s for first-order hyperbolic problems in [22,31]. The development of discontinuous finite-element schemes for secondorder elliptic problems can be traced back to [27] with similar approaches in, for instance, [2,5,29,34]; see also [3].
The design of structure-preserving DG methods is a rather recent topic. Positivitypreserving DG schemes for parabolic equations were developed in, e.g., [9,15,23,33,36]. The positivity preservation is ensured by using a special slope limiter (as in [9,15]), together with a strong stability preserving Runge-Kutta time discretization (as in [33,36]), while in [23], the positivity of the discrete solution is enforced through a reconstruction algorithm, based on positive cell averages. As far as we know, the use of an exponential transformation to ensure the positivity of the discrete solutions within a DG scheme is new. Positivitypreserving schemes for the Fisher-KPP equation were already studied in the literature, but only for finite-difference approximations [17,24], without a convergence analysis, and for continuous finite-element discretizations [35].
Other important properties are entropy stability (the entropy is bounded for all times) and entropy monotonicity (the entropy is nonincreasing). Entropy-stable DG schemes for the compressible Euler and Navier-Stokes equation were studied in [14,28], while a discrete version of the entropy inequality (and hence entropy monotonicity) was proved in [33] for Fokker-Planck-type equations and aggregation models. We are not aware of results in the literature regarding the preservation of the entropy structure of the Fisher-KPP equation on the discrete level.
The paper is organized as follows. We state our notation and some auxiliary results related to the DG method in Section 2. The DG scheme is introduced and studied in Section 3: The existence of a solution to the DG scheme, the discrete entropy inequality, and the exponential decay of the entropy are proved. The convergence of the numerical scheme is proved in Section 4. Finally, Section 5 is devoted to some numerical experiments in one space dimension.

Notation and auxiliary results
We start with some notation. Let T h " tK i : i " 1, . . . , N h u be a family of simplicial partitions of the bounded domain Ω Ă R d for d " 1, 2, 3. The mesh parameter h is defined by h " max KPT h h K , where h K " diampKq. The elements may be tetrahedra in three space dimensions, triangles in two dimensions, and intervals in one dimension. In two and three dimensions, we suppose that T h is shape regular (see, e.g., [30,Section 2.1]) and, for simplicity, without hanging nodes. Our analysis actually extends also to k-irregular meshes [18]. We denote by E h the set of interior faces or edges of the elements in T h .
On the partition T h , we define the broken Sobolev space The traces of functions in H 1 pΩ, T h q belong to the space T pΓ h q " ś KPT h L 2 pBKq, where Γ h is the union of all boundaries BK for all K P T h . The functions in T pΓ h q are single-valued on BΩ and double-valued on Γ h zBΩ.
Let q be a piecewise smooth function and q be a piecewise smooth vector field on T h . We write K´and K`for the two elements sharing the face f , i.e. f " BK´X BK`, and nf or the unit normal vector pointing to the exterior of K˘. Furthermore, we set q˘" q| Kȃ nd φ˘" φ| K˘. Then we define averages: tqu " 1 2 pq´`q`q, tφu " 1 2 pφ´`φ`q, jumps: vqw " q´n´`q`n`, vφw " φ´¨n´`φ`¨n`.
Note that the jump of a scalar function is a vector which is normal to f , and the jump of a vector-valued function is a scalar. The mesh size function h P L 8 pΓ h q is defined by Furthermore, we introduce the finite-element space of degree p P N associated to the partition T h : where P p pKq is the set of polynomials on K with degree at most p, and the space of test functions H 2 n pΩq " tφ P H 2 pΩq : ∇φ¨n " 0 on BΩu. Next, we recall some auxiliary results.
Lemma 1 (Inverse trace inequality; Lemma 2.1 in [32]). Let K P R d (d " 2, 3) be an element with diameter h K , let f be an edge or face of K, and let n f be a unit normal vector normal to f . Then for all polynomials ξ P P p pKq of degree p, there exists a constant C inv ą 0, independent of h K and p, such that Lemma 2 (Multiplicative trace inequality; Lemma A.2 in [30]). Let K be a shape-regular element. Then there exists a constant C ą 0 such that for all ξ P H 1 pKq, Lemma 3 (Discrete Poincaré-Wirtinger inequality; Theorem 4.1 in [7]). There exists a constant C PW ą 0 such that for all ξ P H 1 pΩ, T h q, We also need a compactness result for functions ξ P H 1 pΩ, T h q. For this, we define the DG norm (7) }ξ} DG "ˆ}ξ} 2 Lemma 4 (DG compact embedding; Lemma 8 in [7]). Let pξ h q Ă H 1 pΩ, T h q be a sequence such that }ξ h } DG ď C for all h P p0, 1q and some C ą 0. Then there exists a subsequence ph i q with h i Ñ 0 as i Ñ 8 and a function ξ P H 1 pΩq such that where 1 ď q ă q˚and q˚" 4 for d " 3, q˚" 8 for d " 1, 2.
3. Analysis of the DG scheme: existence and structure preservation The DG discretization of the weak formulation of (4)-(5) reads as follows. Let ε ě 0 and λ 0 ż The form B : H 1 pΩ, T h q 3 Ñ R represents the interior penalty DG discretization of the nonlinear diffusion term. It is linear in the second and third argument and is defined by where αpuq is a stabilization function, given by (10) αpuq " 3 2 C 2 inv`m axtpe u q´, pe u q`u˘2 max expp}u} L 8 pK´q q, expp}u} L 8 pK`q q ( . We recall that the constant C inv is defined in Lemma 1. The third term on the left-hand side of (8) is a regularization term (only) needed for the existence analysis to derive a uniform (but ε-depending) bound for the fixed-point argument. For linear elements p " 1, we may allow for ε " 0; see Appendix A.
3.1. Existence of a discrete solution. We show that problem (8) possesses a solution. First, we prove a coercivity property for the form B.
Lemma 5 (Coercivity of B). The form B, defined in (9), satisfies for all v P H 1 pΩ, T h q, Proof. Definition (9) gives for v P H 1 pΩ, T h q: We estimate the second integral by using Young's inequality: where β f ą 0 is a parameter which will be defined below. The first integral on the righthand side is estimated according to To proceed, we set Taking into account the inverse trace inequality (6), we infer that Consequently, we obtain Inserting this estimate into (11), it follows that With the definitions of αpvq (see (10)) and β f as well as the property h ď h K˘, the difference in the bracket can be computed as This shows that By the definition of the jumps and the mean-value theorem for x P f , We use definition (10) and insert the previous estimate into (12): Since pe v q˘ě expp´}v} L 8 pK˘q q, we have This finishes the proof.
The idea is to apply the Leray-Schauder fixed-point theorem. We define the fixedpoint operator Φ : The left-hand side defines the bilinear form apw, φq, which is coercive, apw, wq " ε}w} 2 L 2 pΩq . The right-hand side defines a linear form which is continuous on L 2 pΩq (using the fact that in finite dimensions, all norms are equivalent). Thus, Φ is well defined by the Lax-Milgram lemma. As the right-hand side of (13) is continuous with respect to w, standard arguments show that Φ is continuous. Furthermore, Φpw, 0q " 0. It remains to prove that there exists a uniform bound for all fixed points of Φ. To this end, let v P V h and σ P r0, 1s such that Φpv, σq " v.
Let spvq " vplog v´1q`1 ě 0. The convexity of s implies that Then, using the test function φ " v in (13) gives, because of the properties Bpv; v, vq ě 0 (Lemma 5) and e v p1´e v qv ď 0, This is the desired uniform bound. We infer the existence of a solution to (8) by the Leray-Schauder fixed-point theorem.

3.2.
Discrete entropy inequality and exponential decay. Let λ k h P V h be a solution to (8). We show that the entropy is nonincreasing with respect to k P N.
Lemma 7 (Discrete entropy inequality). Let ε ě 0 and let λ k h P V h be a solution to (8). Then where the constant C 0 ą 0 only depends on C inv and C PW from Lemmas 1 and 3.
Proof. We take φ h " λ k h as a test function in (8) and use inequality (14) to find that It remains to estimate the first term on the right-hand side. For this, we use the coercivity estimate of Lemma 5 and the discrete Poincaré-Wirtinger inequality from Lemma 3: Setting C 0 " 2 mint1, C 2 inv uC´2 PW finishes the proof. We wish to bound the total mass ş Ω exppλ k h qdx from below and above. Since spuq " uplog u´1q`1 is invertible only on r0, 1s and on r1, 8q but not globally on r0, 8q, we introduce the following functions: In particular, σ´˝s " id on r0, 1s and σ`˝s " id on r1, 8q.
Lemma 8 (Bounds for the total mass). Let ε ě 0 and let λ k h be a solution to (8). Then Observe that if S 0 h ă |Ω|, the lower bound σ´pS 0 h {|Ω|q is positive. Thus, the total mass can never vanish, which excludes the case of solutions converging for k Ñ 8 to the zero solution. The reason for the difference between S 0 h ă |Ω| and S 0 h ě |Ω| lies in the fact that (4)-(5) admits two steady states, λ k h " 0 (corresponding to u k " e λ k h " 1) and λ k h "´8 (corresponding to u k " 0). The assumption S 0 h ă |Ω| will be crucial to prove the decay estimate for the entropy; see Proposition 9. We discuss the case S 0 h ě |Ω| in Remark 12. Proof of Lemma 8. First, we show the lower bound. If S 0 h ě |Ω|, we have σ´pS 0 h {|Ω|q " 0, and there is nothing to prove. Thus, let S 0 h ă |Ω|. Set β k " mint1, exppλ k h qu ď 1. As s is convex, we infer from Jensen's inequality and spβ k q " 0 for λ k h ą 0 that where in the last step we have used the monotonicity of k Þ Ñ S k h . With this preparation, we are able to verify the lower bound. As σ´is decreasing, we find that For the upper bound, we can assume that ş Ω exppλ k h qdx ě |Ω|, since otherwise, the inequality is trivially satisfied in view of σ`pvq ě 1. By the concavity of σ`, we can again apply the Jensen inequality: proving the claim.
Proposition 9 (Discrete entropy decay). Let ε ě 0 and let λ k h be a solution to (8). We assume that S 0 h ă |Ω|. Then there exists a constant C 1 ą 0, only depending on S 0 h , such that for all k P N, The proof is based on two properties: The diffusion drives the solution towards a constant, while the reaction term guarantees that there is only one (positive) steady state. In order to cope with the interplay of diffusion and reaction, we prove first the following lemma.
Lemma 10. Introduce for θ ą 0 the functions Then can be continuously extended to v " 0 (with value gp0q " 1{2) and it is decreasing with limits lim vÑ8 gpvq " 0 and lim vÑ´8 gpvq "`8. Therefore, gpvq ď gplog θq " M 1 pθq for all v ě log θ, showing the first inequality. For the second one, let v ď log θ. Then spe v q ă 1 for v ď 0 and the monotonicity of v Þ Ñ spe v q for v ě 0 implies that spe v q ď spθq. Thus, for any v P R, spe v q ď maxt1, spθqu " M 2 pθq, completing the proof.
Proof of Proposition 9. The idea of the proof is to split S k h into two integrals, for some suitably chosen α ą 0 and to estimate these integrals by the second and third terms on the left-hand side of the discrete entropy inequality (16).
We turn to the first integral on the right-hand side of (19). We claim that there exists a constant C ε 0 θ ą 0 such that To prove this inequality, we begin by showing that ş Ω exppλ k h {2qdx is bounded from below. Indeed, using the monotonicity of s˝exp in r0, 1s and of k Þ Ñ S k , This yields the lower bound Therefore, as long as is positive. Squaring this expression and integrating over tλ k h ď logpε 0 θqu thus does not change the inequality sign: Combining the estimate of Lemma 10 and the previous estimate, we arrive at This proves claim (20) with Next, we estimate the second integral on the right-hand side of (19). It follows from Lemma 10 that ż Therefore, (19) gives and solving this recursion shows the proposition.
Theorem 11 (Decay in the L 1 pΩq norm). Let the assumptions of Proposition 9 hold. Then there exists a constant C 2 ą 0, only depending on S 0 h and |Ω|, such that where η P p0, 1q and C 1 ą 0 are as in Proposition 9.
Proof. To simplify the notation, we set u " e λ k h andū " |Ω|´1 ş Ω e λ k h dx. Then the Csiszár-Kullback inequality (see, e.g., [4, (2.8)]) gives using the property spuq ě 0 for all u ě 0. We know from Lemma 8 thatū is bounded from above by σ`pS 0 h {|Ω|q. Hence, It remains to show that a similar estimate holds for |ū´1|. Since the entropy density s is convex, Jensen's inequality shows that (22) spūq " sˆ1 |Ω| It holds spvq ă 1 if and only if v ă e. Consequently, we haveū ă e. Applying the elementary inequality spuq ě pu´1q 2 pe´1q 2 for all 0 ď u ď e to u "ū and using (22) gives Thus, combining (21) and the previous inequality, we conclude that and the proof follows after applying Proposition 9.
Remark 12. We discuss the case S 0 h ě |Ω|. Fix △t P p0, 1q and L P N with L ą 1. Define λ k h " pL´kq`logp1´△tq, where z`" maxt0, zu denotes the positive part of z P R. Then e λ k h " p1´△tq L´k ă 1 for k ă L and e λ k h " 1 for k ě L. Consider the case L ą k " 1. Then, setting δ :" p1´△tq L´k , we estimate By the convexity of s, it follows that spuq´spvq ď pu´vqs 1 puq " pu´vq log u for all u, v ą 0. Since λ k h ď 0 for k ď L, (23) yields which directly implies the entropy inequality (16). This inequality is trivially satisfied for k ě L. However, it holds for L " 2k that This means that if S 0 h ě |Ω|, there exists no constant C ą 0 depending only on S 0 h such that (18) holds for all pλ k h q Ă L 2 pΩq satisfying the entropy inequality (16). Note that the constructed function e λ k h does not possess a uniform positive lower bound.

Analysis of the DG scheme: numerical convergence
We show first that the solutions to (8) are uniformly bounded in the DG norm (7) if the initial entropy S 0 h is bounded uniformly in h.
Lemma 13 (Uniform bound in DG norm). Let ε ě 0 and let λ k h be a solution to (8). Then there exists a constant C ą 0 such that Proof. We have shown in the proof of Lemma 7 that Then, by definition of the DG norm, Using the inequality u ď 2`spuq for u ě 0, applied to u " e λ k h , and the monotonicity of k Þ Ñ S k h , we find that If 2 mint1, C 2 inv u△t ď 1 then On the other hand, if 2 mint1, C 2 inv u△t ą 1, we have, again by the monotonicity of k Þ Ñ S k h , such that in either case, proving the lemma.
h Ñ u k´1 strongly in L 2 pΩq as pε, hq Ñ 0. Then there exists a unique strong solution u k P H 2 n pΩq to (24) 1 △t pu k´uk´1 q " ∆u k`uk p1´u k q in Ω, ∇u k¨n " 0 on BΩ such that e λ k h Ñ u k strongly in L 2 pΩq as pε, hq Ñ 0.
Proof. Let λ k h P V h be a solution to (8).
Step 1: We claim that there exists a subsequence pε i , h i q Ñ 0 such that e λ k h i Ñ u k strongly in L 2 pΩq as i Ñ 8.
Indeed, by assumption, the initial entropy pS 0 h i q iPN is bounded. Then Lemma 13 implies that e λ k h {2 is bounded in the DG norm uniformly in ε and h. By the compactness Lemma 4, there exists a subsequence pε i , h i q Ñ 0 and a function v k P H 1 pΩq satisfying Consequently, e λ k h i Ñ pv k q 2 ": u k strongly in L 1 pΩq. The discrete entropy inequality (16)  h i Ñ w k weakly in L 1 pΩq as i Ñ 8, for some function w k . We deduce from the strong L 1 convergence of e λ k h i , possibly for another subsequence, that e 2λ k h i Ñ pu k q 2 " w k a.e. in Ω. This implies that thus proving the desired L 2 convergence.
Step 2: We claim that for any φ P H 2 n pΩq X C 1 pΩq, it holds that 1 △t Since φ does not necessarily belong to V h , we cannot use it as a test function in the weak formulation (8). Therefore, let P h : C 0 pΩq Ñ C 0 pΩq X V h be the interpolation operator, defined, e.g., in [10, Section 2.3]. It possesses the following property [10, Section 3.1.6]: There exists a constant C I ą 0 such that for all K P T h and φ P H 2 pKq, for m ď 2 ď q such that m´d{q ď 2´d{2. In particular, for φ P H 2 pΩq and d ď 3, For given φ P H 2 n pΩq X C 1 pΩq, we choose the test function φ h i :" Here, we have used the fact that vφ h i w " 0 since φ h i is continuous. Note that (28) implies that φ h i Ñ φ strongly in L 8 pΩq as i Ñ 8. As e λ k´1 h i Ñ u k´1 strongly in L 2 pΩq, by assumption, we have for the left-hand side of (29): Similarly, as e λ k h i Ñ u k strongly in L 2 pΩq, we infer for the first and last integrals on the right-hand side of (29) that Inequality (15) shows that This implies that the second integral on the right-hand side of (29) converges to zero.
Next, we prove for the third integral on the right-hand side of (29) that Indeed, by the Hölder inequality, the interpolation property (27), and the discrete entropy inequality (16), we obtaiňˇˇˇÿ It remains to prove for the fourth integral on the right-hand side of (29) that To this end, we use the elementary inequality |tu∇vu| ď 2tuut|∇v|u for functions u, v with nonnegative u and the Cauchy-Schwarz inequality:ˇˇˇÿ We estimate both integrals separately. First, the multiplicative trace inequality in Lemma 2 shows that, for some constant C ą 0 and for faces or edges f " BK`X BK´, We deduce from (27), i.e.
Therefore, also using hpxq ď h i , we deduce from (30) thaťˇˇˇÿ where h i pxq " minth i,K`, h i,K´u for x P BK`X BK´. We claim that the sum on the right-hand side is bounded uniformly in h i . By Definition (10), such that we can estimate where we used (12) in the last step. The proof of Lemma 7 shows that Bpλ k We put together all the previous convergence results to infer that 1 △ Thus, inserting (29), all integrals involving φ h i cancel, and we end up with (26).
Step 3: We prove that the limit u k , derived in Step 1, is a solution to the very weak formulation for all φ P H 2 n pΩq X C 1 pΩq. For the proof, we pass to the limit h i Ñ 0 in each term of (26). Because of (25), we have ż The limit i Ñ 8 in the remaining expression is more delicate. Consider the first term in the definition of I i . Integrating by parts elementwise gives where we have used the fact that ∇φ has a continuous normal component across interelement boundaries. From the previous identity and the L 2 convergence of e We claim that since this implies that and thus shows (32). For the proof of (33), let x P BK`X BK´for two neighboring elements K`, K´P T h i and set λ˘:" λ k h i | K˘. We assume without loss of generality that λ`ě λ´since otherwise, we may exchange K`and K´. The definitions of the jump v¨w and average t¨u imply thaťˇˇv e λ´ˇˆp e λ`´λ´´1 qn`´pλ`´λ´qn`1 2 pe λ`´λ´`1 q˙¨∇φˇˇˇď e λ´ˇp e λ`´λ´´1 q´pλ`´λ´q 1 2 pe λ`´λ´`1 qˇˇˇˇ|∇φ| ": e λ´| gpλ`´λ´q||∇φ|, where gpsq " pe s´1 q`spe s`1 q{2 for s ě 0. A Taylor expansion shows that gpsq " g 2 pξqs 2 {2 "´ξe ξ s 2 {4 for some 0 ď ξ ď s. Therefore |gpsq| ď s 2 e 2s {4 for s ě 0, and we obtain e λ´| gpλ`´λ´q| ď e λ4 pλ`´λ´q 2 e 2pλ`´λ´q " 1 4 pλ`´λ´q 2 e 2λ`´λď The difference can be identified with the jump of λ k h i across f " BK`X BK´, while the sum corresponds to the average of e 2λ k h i in f . Thus, it follows from (34) thaťˇˇˇÿ By definition (10) of the stabilization factor, it holds that Using this estimate and the coercivity estimate (12) for the form B, we can writěˇˇˇÿ We know from the proof of Lemma 7 that Bpλ k h i ; λ k h i , λ k h i q ď C{△t is uniformly bounded. This proves our claim (33). Now, we can pass to the limit i Ñ 8 in (26), which yields (32).
Step 4: We claim that the solution u k P L 2 pΩq to (26) satisfies the regularity u k P H 1 pΩq and hence is a weak solution to (4)- (5). To this end, we use the duality method as in [6, p. 318]. Let T : L 2 pΩq Ñ H 2 n pΩq be defined by T v " u, where u solves the elliptic problem u´△t∆u " v in Ω, ∇u¨n " 0 on BΩ. By [25,Theorem 8.3.10], for v P C 8 0 pΩq, it holds that T v P H 2 n pΩq X C 1 pΩq. Then, introducing g :" u k´1`△ tu k p1´u k q, the very weak formulation (32) can be equivalently written as ż Ω u k pφ´△t∆φqdx " ż Ω gφdx for all φ P H 2 n pΩq X C 1 pΩq. Given v P C 8 0 pΩq, we set φ " T v, and the previous equation becomes (35) ż Ω u k vdx " ż Ω gT vdx.
As C 8 0 pΩq is dense in L 2 pΩq and T is continuous, (35) remains valid for all v P L 2 pΩq. Next, we denote by T 1 : H 2 n pΩq 1 Ñ L 2 pΩq the dual operator of T . According to [25,Theorem 8.3.10], the operator T can be extended to an operator T : L p pΩq X H 1 pΩq 1 Ñ W 2,p pΩq for 1 ă p ď 2. (This is basically a regularity statement for the elliptic problem.) We deduce from the Sobolev embedding theorem that W 2,p pΩq ãÑ C 0 pΩq for p ą 3{2 since d ď 3. Therefore, there exists an extension T˚: C 0 pΩq 1 Ñ L p 1 pΩq of T 1 , where p 1 " p{pp´1q ă 3. Now, g P L 1 pΩq Ă C 0 pΩq 1 . Then (35) implies that u k " T˚pgq P L p 1 pΩq for p 1 ă 3 and consequently, g P L q pΩq for q ă 3{2.
It remains to show that u k P H 1 pΩq. Since u k P L 2 pΩq, the elliptic problem Thus, pv m q is bounded in L 2 pΩq and it follows the existence of a subsequence which is not relabeled that v m á u k weakly in L 2 pΩq as m Ñ 8. Using v " v m´△ t∆v m in (35), it follows that We apply the Hölder inequality and use the Sobolev embedding H 1 pΩq ãÑ L q 1 pΩq for q 1 ď 6, knowing that g P L q pΩq for q ă 3{2: where 3 ă q 1 ď 6 and 1{q`1{q 1 " 1. The H 1 pΩq norm can be absorbed by the corresponding term on the right-hand side of (36), and we end up with This shows that pv m q is bounded in H 1 pΩq. Thus, there exists a subsequence (not relabeled) which converges weakly in H 1 pΩq to some function w P H 1 pΩq. Since v m á u k weakly in L 2 pΩq, we conclude that w " u k P H 1 pΩq.
Knowing that u k P H 1 pΩq solves (32), we can integrate by parts in the first term of the right-hand side of (32), leading to 1 △t ż Ω pu k´uk´1 qφdx "´ż Ω ∇u k¨∇ φdx`ż Ω u k p1´u k qφdx for all φ P H 2 n pΩq and, by density, for all φ P H 1 pΩq. Moreover, since u k P L 4 pΩq and consequently, u k p1´u k q P L 2 pΩq, elliptic regularity implies that u k P H 2 pΩq and ∇u k¨n " 0 on BΩ, i.e. u k P H 2 n pΩq. We conclude that u k solves (24).
Step 5: We prove the uniqueness of weak solutions to (4)-(5) to conclude the convergence of the whole sequence e λ k h to u k . Let u k , v k P H 1 pΩq be two weak solutions to (4)-(5).
Taking the difference of the corresponding weak formulations with the test function u k´vk , we obtain Thus, choosing △t ă 1, we infer that u k´vk " 0 in Ω.

Numerical results
We present some numerical results for the Fisher-KPP equation in one space dimension, 5.1. One group of species. Let D " 10´4 and u 0 pxq " 0.8 for 0 ă x ă 1{2, u 0 pxq " 0 else. Problem (37)-(38) models the evolution of one species initially concentrated in the domain p0, 1{2q. We solve this problem by using an implicit Euler scheme in time and a continuous P 1 finite-element discretization, both on a uniform mesh. The reaction term is treated implicitly. The Newton method with relaxation is used at each time step, up to convergence. The integrals are computed by using a Gauß-Legendre quadrature formula of order 8. Figure 1 shows the density upx, tq at various time instances. We observe that the finite-element solution u k h becomes negative even on the finer mesh, and it is pushed towards´8 in some region since u˚" 0 is a repulsive steady state. These results motivate the introduction of the exponential transformation u " exppλq. We are choosing the same initial datum as before but choosing u 0 pxq " 10´1 6 instead of u 0 pxq " 0 to allow for the exponential transformation. Figure 2 shows the solutions to the continuous P 1 finite-element approximation associated to problem (4)-(5) in the variable λ k h . The implicit nonlinear scheme is solved again by Newton's method with relaxation at each time step. The integrals are solved again by a Gauß-Legendre quadrature formula of order 8. Note that if p " 1, the integrals are of the type ş K e ax`b pcx`dqdx and thus can be integrated exactly. The discrete densities exppλ k h q are positive by construction. However, we need more relaxation in the Newton method when higher-order schemes p ą 1 are used, which slows down the algorithm.  Figure 2. Continuous P 1 finite-element discretization of problem (4)-(5) in the variable λ k , using N el " 20 elements (left) and N el " 40 elements (right). The time step size is in both cases △t " 1{6 and the end time is T " 20.
Therefore, we employ a discontinuous Galerkin method with polynomial order p ě 1 for problem (4)-(5) in the variable λ k h ; see scheme (8). The regularization term is not necessary for the numerics, i.e., we set ε " 0 in (8) for our simulations. Figure 3 illustrates the discrete solutions for polynomial orders p " 1, 2, 3, indicating that the method is stable with respect to the order. The jumps are due to the discontinuous Galerkin method. Figure 4 represents the discrete solutions with the same numerical parameters as in Figure 3 but with the initial datum u 0 pxq " 1 for 0 ă x ă 1{2 and u 0 pxq " 0 else. Also in this example, the lower and upper bounds 0 ď exppλ k h q ď 1 are always satisfied.

5.2.
Entropy decay. Proposition 9 shows that the discrete entropy S k h decays exponentially fast if S 0 h ă |Ω| " 1. To illustrate this behavior numerically, we consider the one-group model with initial condition u 0 pxq " 1 for 0 ă x ă 1{2 and u 0 pxq " 10´1 6 else. Figure 5 (left) shows that there are two different slopes. For small times, the entropy decay is rather slow. When the time step k is sufficiently large such that exppλ k h q ą ε for some ε ą 0, the reaction dominates, and the entropy decay becomes faster. This behavior becomes even  Figure 3. Reference solution computed from the P 1 finite-element scheme with N el " 300 elements and △t " 1{3 (left top) and solutions computed from the DG scheme with N el " 40 elements, △t " 1{3, and polynomial order p " 1 (right top), p " 2 (left bottom), and p " 3 (right bottom). The initial datum is u 0 pxq " 0.8 for 0 ă x ă 1{2 and u 0 pxq " 10´1 6 else.
more apparent in the case of pure diffusion (i.e. without reaction terms), illustrated in Figure 5 (right). We remark that in this situation, the total mass is conserved numerically. Figure 6 shows the entropy decay in semi-log scale for the initial data u 0 pxq " n for 0 ă x ă 1{n and u 0 pxq " 10´1 6 else for n " 3, 6, 12. Then 1dx " log n is larger than |Ω| " 1 if n ą e. We observe a region in which the decay rate is very small and it becomes smaller when n (and S 0 h ) increases. This indicates that the assumption S h 0 ă |Ω| is not just technical to derive exponential time decay.  Figure 4. Reference solution computed from the P 1 finite-element scheme with N el " 300 elements and △t " 1{3 (left top) and solutions computed from the DG scheme with N el " 40 elements, △t " 1{3, and polynomial order p " 1 (right top), p " 2 (left bottom), and p " 3 (right bottom). The initial datum is u 0 pxq " 1 for 0 ă x ă 1{2 and u 0 pxq " 10´1 6 else.

Traveling waves.
We are looking for traveling-wave solutions to (37) with D " 1. Setting upx, tq " φpsq with s " x´ct, the Fisher-KPP equation can be rewritten as a system of first-order differential equations: We choose the initial data φp0q " 1 and ψp0q "´10´1 0 . The (reference) traveling-wave solution φpsq, computed from (39) using the Matlab command ode45, is compared in Figure  7 with the DG solution exppλ k h q, computed from the DG scheme (8), and the continuous P 1 finite-element solution u h . The solutions are shown at the time instances t " 0, t " T {6, t " T {3, and t " T {2 (with T " 20). Both approximations are diffusive, i.e., the travelingwave speed is overestimated by the DG and finite-element solutions. On the finer mesh with N el " 80 elements, the DG solution is clearly less diffusive compared to the other  Figure 6. Entropy decay for the one-group model with the initial datum u 0 pxq " n for 0 ă x ă 1{n and u 0 pxq " 10´1 6 else. discrete solutions on the coarser mesh with N el " 50 elements. Better approximations are expected by using a higher-order time discretization. Structure preservation of higherorder temporal approximations is a delicate topic (see, e.g., [21]) and will be studied in a future work. Also an error analysis is postponed to a future work.
Appendix A. Linear elements: Existence of solutions for ε " 0 We show that the regularization term ε ş Ω λ k h φ h dx in the DG scheme (8) is not needed if we consider linear elements.  Figure 7. Comparison of the traveling-wave solution φpsq, the DG solution exppλ k h q (with N el " 50 or N el " 80 and △t " 1{3), and the finite-element solution u h .
The proof is based on the idea of the proof of [1, Lemma 3.10]. Let ε ą 0 and let λ k h P V h be a solution to (8) given by Proposition 6. In order to emphasize the dependency on ε, we write λ ε :" λ k h . Our goal is to derive an ε-uniform L 8 pΩq bound for λ ε .
Step 1: We derive first some estimates for e λε . Lemma 8 shows that where δ " σ´pS 0 h {|Ω|q and M " σ`pS 0 h {|Ω|q. Since we assumed that S 0 h ă |Ω|, we have δ ą 0. Using the coercivity estimate (12), the inequality (17), and e λε pe λε´1 qλ ε ě 0, we find that where we used in the last step the monotonicity of k Þ Ñ S k h , guaranteed by Lemma 7. Consequently, Step 2: We claim that for any ε ą 0, there exists an element K ε P T h and a constant µ ą 0, independent of ε, such that For the proof, let N P N be the number of elements in T h . The lower and upper bounds (40) imply the existence of an element K ε P T h such that By the mean-value theorem, there exists x ε P K ε such that Then (43) gives Since λ ε is a polynomial of degree one on K ε , by assumption, its gradient is constant on K ε , and we deduce from (41) and (43) that Combining the last two estimates, it follows for x P K ε that which shows the claim.
Step 3: We wish to prove a uniform L 8 bound for λ ε on the faces or edges of K ε . Let µ ą 0 and K ε P T h such that }λ ε } L 8 pKεq ď µ. Set K´:" K ε and consider neighboring elements K`P T h satisfying f P BK´X BK`‰ H. Furthermore, let λ˘" λ ε | K˘. We claim that there exists C µ ą 0, independent of ε, such that (44) }λ`} L 8 pf q ď C µ .
Step 4: Denote by pφ 0 , . . . , φ d q the basis of P 1 pK`q such that φ i pe j q " δ ij for i, j " 0, . . . , d, where the vertices e i of K`are ordered in such a way that e 0 R f . Then we can formulate λ`on K`as λ`pxq " where a ε i " λ`pe i q. Estimate (44) shows that λ`pe i q is uniformly bounded at the vertices a 1 , . . . , a d of K`, i.e. |a ε i | ď C µ for all i " 1, . . . , d.
Step 5: We wish to estimate the remaining vertex e ε 0 that is not an element of K´" K ε . We claim that there exist constants L µ ď U µ , being independent of ε, such that L µ ď a ε 0 ď U µ . We first prove the upper bound. Using the bound for a ε i for i " 1, . . . , d, we have λ`ě a ε 0 φ 0´d If a ε 0 ď 0, there is nothing to show. Otherwise, it follows from (43) that Then, setting c 0 " ş K`φ 0 dx, we infer that a ε 0 ď Me Cµd {c 0 ": U µ . The proof of the lower bound is more involved. Let f 0 be the face or edge that is opposite of the vertex e 0 , and let f 1 , . . . , f d be the remaining faces or edges. For later use, we note that the integrals Ipbq`Jpbq ď 1 for all b ď L 1 µ . We estimate |∇λ`¨n| ě |a ε 0 ||∇φ 0¨n |´d ÿ i"1 |a ε i ||∇φ 0¨n | ě |a ε 0 ||∇φ 0¨n |´C µ d max i"1,...,d |∇φ i |.
Step 6: Combining the previous steps, we infer that there exists a constant gpµq ą 0 such that }λ ε } L 8 pK`q ď gpµq. This estimate means that if λ ε is bounded in some element with constant µ, then λ ε is bounded in the neighboring elements with constant gpµq. Now, take an arbitrary element K P T h . Then there exists a finite sequence K 0 , K 1 , . . . , K m of elements with K 0 " K ε and K m " K such that K j´1 and K j are neighboring elements. Repeating the arguments of Steps 3-5, the bound }λ ε } L 8 pK 1 q ď µ 1 :" gpµq implies that }λ ε } L 8 pK 2 q ď gpµ 1 q " gpgpµqq. Thus, by iteration, }λ ε } L 8 pKq ď pg˝¨¨¨˝gq looooomooooon m times pµq.
The upper bound is independent of ε and holds for all elements K P T h . Consequently, pλ ε q is bounded in L 8 pΩq.
We deduce that there exists a subsequence (not relabeled) such that λ ε Ñ λ strongly in L 8 pΩq, recalling that V h is finite-dimensional. In fact, the convergence holds in any norm.
Thus, we can pass to the limit ε Ñ 0 in (8), and the limit equation is the same as (8) with ε " 0.