Higher-order convergence analysis for interior and boundary layers in a semi-linear reaction-diffusion system networked by a $k$-star graph with non-smooth source terms

1.   Introduction

Taxis-diffusion and aggregation equations are widely studied in the context of biological populations (see [10,15,16,21] for instance). They describe cell communities which react to external stimuli and form aggregates of organisms (pattern formation), such as bacterial colonies, slime mold or cancer cells. The Patlak-Keller-Segel model [18] is the most famous system and we are interested in the following generalization

Here, $u(x, t) \geq 0$ represents the density of a given quantity (e.g. cells or bacteria population) and the initial data $u^0(x)$ is a given nonnegative smooth function. As for the function $v$, which models a molecular concentration, we choose either the case of the Fokker-Planck (FP in short) equation, where $v(x)$ is known

or the case of the generalized Keller-Segel (GKS in short) equation, where

Depending on the modeling choice for $\varphi(u)$, solutions to 1 can blow up in finite time depending upon a critical mass (see [6,22]) or reach stationary profiles in the form of peaks or plateaus [24] (pattern formation by Turing instability). The high nonlinearities due to the advection term make problem 1 mainly untractable through analytical methods. Thus, it is important to provide reliable numerical methods avoiding non-physical oscillations and numerical instabilities even when dealing with non-smooth solutions. The main properties that one wishes to preserve in a numerical method are

$\rm (P1)$ positivity property, since we are dealing with densities or concentrations,

$\rm (P2)$ mass conservation, because no-flux boundary conditions are imposed,

$\rm (P3)$ preservation of discretized steady states of the form

where $\mu$ is a parameter related to the mass of $u$, and

$\rm (P4)$ energy dissipation

where $G(u)$ is a primitive of $g(u)$ and the value of $\kappa$ differs for the two cases we study here, namely

The aims of our work are first to recall two points of view for the derivation of the above energy inequality, second to use them for the construction of conservative, finite volume numerical schemes preserving energy dissipation to solve equation 1, third to make numerical comparisons in the case of complex patterns in order to distinguish physical instabilities from numerical artifacts. The two different derivations of the energy dissipation use two symmetrization strategies: the gradient flow or the Scharfetter-Gummel approach. It turns out that they lead to two strategies for discretization of problem 1. We prove that the proposed schemes statisfy properties 4-7 and because we build implicit schemes, there is no limitation on the time step in the fully discrete case.

There exist other works which propose schemes for the resolution of problems in the form 1. For instance, finite elements methods are used, see [25] and references therein. Optimal transportation schemes for Keller-Segel systems are introduced in [5]. The papers [8] and [9] propose a finite-volume method able to preserve the above properties, including energy dissipation, at the semi-discrete level or with an explicit in time discretization, using the gradient flow approach, see also [4]. The symmetrization using the Scharfetter-Gummel approach is used in [20] where properties similar to ours are proved. However, the results do not include sensitivity saturation. To the best of our knowledge, our work is the first to propose implicit in time methods, without time step limitation (CFL condition), for which we are able to prove that, under generic conditions, the energy decreases at both semi-discrete and discrete level. Moreover, we build an alternative to the gradient flow approach applying the Scharfetter-Gummel strategy [26] for the discretization of drift-diffusion equations 1 with a general saturation function $\varphi$.

The paper is organized as follows. In Section 2, we present in more details our assumptions for the equation 1. We also explain some modeling choices in particular for the nonlinearity $\varphi (u)$ and on the choice of the kernel $K$. Section 3 is devoted to the introduction of the two approaches, gradient flow or Scharfetter-Gummel, and to how we use the continuous version of energy dissipation to derive the schemes. In Section 4 and Section 5, we show how the aforementioned two approaches lead to two different numerical methods, developed from the semi-discrete (only in space discretization) level to the fully discretized scheme. In particular, using a general result recalled in Appendix A about monotone schemes, we prove that the proposed schemes are well-posed and satisfy the fundamental properties 4-7. These theoretical results are illustrated in Section 6 by numerical simulations: we compare the gradient flow and the Scharfetter-Gummel schemes with the upwind approach, typically used to solve this kind of models.

2.   Modeling assumptions

The standard biological interpretation of 1 ([14,21,23]) provides us with some further properties of the nonlinearities which we describe now.

Chemotactic sensitivity. The function $\varphi (u)$ is called chemotactic sensitivity. It determines how the random movement of particles of density $u$ is biased in the direction of the gradient of $v$. In order to include the different choices of $\varphi$, as $\varphi (u) = u$ as in the Keller-Segel or drift-diffusion model, or the logistic case $\varphi (u) = u(1-u)$, or the generalized case $\varphi (u) = u e^{-u}$, we use the formalism

More precisely, we consider two cases for the smooth function $\psi$,

or

In the case 11 we only consider solutions which satisfy $u\in [0, M]$.

It is convenient to introduce the notations

where $a$ is a constant chosen so that $g$ is well-defined (depending on $\varphi$). For instance, for the standard case $\varphi (u) = u$, $a = 1$ and one obtains $g(u) = \ln (u)$ and $G(u) = u\ln(u) -u$. For functions $\varphi$ satisfying 10, a natural hypothesis which is related to blow-up is the following

an assumption which, as we see it later, appears naturally when it comes to the well-posedness of numerical schemes.

Note that under assumption 10 and if $\psi$ is bounded, solutions exist globally and are uniformly bounded [11]. Under assumption 11, if $0\leq u^0 \leq M$, the solution is also globally defined and satisfies $0\leq u(t, \cdot) \leq M$ for all times [16].

Expression of the drift $v$. The convolution expression for $v$ as a function of $u$ has been widely used in recent studies [1,2,3]. It also comes from the Keller-Segel model [15,18,23], where the equation for the cells density in 1 is complemented with a parabolic equation for the chemoattractant concentration $v$. Since the chemoattractant is supposed to diffuse much quicker than the cells density, we can consider a simplified form of the Keller-Segel system and couple 1 with the elliptic equation for $v$

This equation leads to 3 using the Green function given by the positive and symmetric kernel $K(x, y)$ defined as

3.   Energy dissipation

Energy dissipation is the most difficult property to preserve in a discretization and methods might require corrections [17]. Therefore, it is useful to recall how it can be derived simply for the continuous equation. We focus on two different strategies, that lead to two different discretization approaches, the gradient flow approach and the Scharfetter-Gummel approach.

3.1.   The gradient flow approach to energy.

Using the notations 12, the equation for $u$ can be rewritten as

so that

Consequently, we find, in the Fokker-Planck case

and in the generalized Keller-Segel case

because, thanks to the symmetry assumption on $K$ and by using the definition 3 of $v$, we have

3.2.   The Scharfetter-Gummel approach to energy.

Inspired from the case of electric forces in semi-conductors, the equation for $u$ can be rewritten as

so that

It is immediate to see that the last term has the negative sign while the time derivative term is exactly the same as in the gradient flow approach.

At the continuous level, these two calculations are very close to each other. However, they lead to the construction of different discretizations. The gradient flow point of view is used for numerical schemes by [10], the Scharfetter-Gummel approach is used in [20].

4.   Semi-discretization

We give here our notations for the semi-discretization. We consider a (small) space discretization $\Delta x = \frac{1}{I}$, $I \in \mathbb{N}$. The mesh is centered at $x_i = i \Delta x$ for $i = 1, \ldots, I$, with endpoints $x_{i+1/2} = (i+1/2) \Delta x$ for $i = 1, \ldots, I$. Therefore, our computational domain is always shifted and takes the form $\big( \frac { \Delta x }{2}, (I+ \frac 12) \Delta x \big)$. Finally, the mesh is formed by the intervals

The semi-discrete approximation of $u(x, t)$ at a given time $t$, interpreted in the finite volume sense ([7,13,19]), is denoted by

As for the discretization of $v$ for $i = 1, \ldots, I$, $v_i(t)$ stands for $v(x_i)$ in the FP case, while it is defined by $\sum_{j = 1}^I K_{ij} u_j(t)$ in the GKS case, where

Integration on the interval $I_i$ yields fluxes $F_{i+ \frac 12}(t)$ which will approximate the quantity $\left( \frac { \partial u }{ \partial x} - \varphi (u) \frac { \partial v}{ \partial x}\right)$ at $x_{i+ \frac 12}$ for $i = 0, \ldots, I$ through the interval interfaces. The boundary conditions require $F_{ \frac 12}(t) = F_{I+ \frac 12}(t) = 0$, and the whole problem will be to properly define $F_{i+ \frac 12}(t)$ for $i = 1, \ldots, I-1$. This will be chosen depending on the strategy, either through a gradient flow or a Scharfetter-Gummel approach, as well as the equation considered (FP or GKS).

The mass conservative form of 1 leads to a finite volume semi-discrete scheme

We use the definition 8 for $\kappa$ and set for $i = 1, \dots, I$

The semi-discrete energy is then

4.1.   The gradient flow approach

Using the form 15 of equation 1, we define the semi-discrete flux as

The precise expression of $\varphi _{i+1/2}$ is not relevant for our present purpose which is to preserve the energy dissipation property. However, for stability considerations it is useful to upwind, an issue which we shall tackle when we consider the full discretization.

Then, the semi-discrete energy form is obtained after multiplication by $(g (u_i) -v_i)$ and yields

Therefore, we find the semi-discrete form of energy dissipation

4.2.   The Scharfetter-Gummel approach

We choose to discretize the form 16, defining the semi-discrete flux as

where, again, the specific form of the interpolant $\big( e^{v-g(u)} \varphi (u) \big)_{i+1/2}$ is not relevant here.

As above, the semi-discrete energy form follows upon multiplication by $g (u_i)-v_i$ and reads

Summing up, the semi-discrete form of energy dissipation here writes

4.3.   Discrete steady states

Steady states make the energy derivative vanish which imposes both in the gradient flow and the Scharfetter-Gummel approaches that $\big( g(u_{i+1}) - v_{i+1}\big) = \big( g(u_{i}) - v_{i} \big)$. In other words they are given, up to a constant $\mu$, as the discrete version of 6,

We recall from [24] that in the GKS case, there are several steady states and the constant ones can be unstable.

5.   Fully discrete schemes

For the time discretization, we consider (small) time steps $\Delta t >0$, and set $t^n = n \Delta t$ for $n \in \mathbb{N}$. The discrete approximation of $u(x, t)$ is now

Integration on the interval $I_i$ yields fluxes $F_{i+ \frac 12}^n$ which will approximate the quantity $\left( \frac { \partial u }{ \partial x} - \varphi (u) \frac { \partial v}{ \partial x}\right)$ at $x_{i+ \frac 12}$ for $i = 1, \ldots, I-1$ through the interval interfaces, and at time $t^n$ for $n\in \mathbb{N}$. The boundary conditions require $F_{ \frac 12}^n = F_{I+ \frac 12}^n = 0$ for all $n \in \mathbb{N}$.

To achieve the time discretization, and restricting our analysis to the Euler scheme, we write the time discretization $\frac {du_i(t)}{dt}$ as $\frac {u_i^{n+1} - u_i^{n}}{ \Delta t }$. Therefore, the full discretization of 17 reads for $n \in \mathbb{N}$ as

The issue here is to decide which terms (in $u$ and $v$) should be discretized with implicit or explicit schemes based on fully discrete energy dissipation. We claim that, apart from the interpolant, we need to make the terms in $u_i$ implicit and, for the GKS case, the terms in $v_i$ explicit, a fact on which we now elaborate.

We define the energy at the discrete level through

and

The computation made in the semi-discrete case, $\frac {dE_i (t) }{dt} = \frac {du_i (t) }{dt} (g(u_i(t)) - v_i(t))$, extends to the fully discrete setting and leads to the following constraint on the energy

Here, $u_i^{\alpha_n} : = \alpha u_i^{n} + (1-\alpha) u_i^{n+1}$ and $v_i^{\beta_n} : = \beta v_i^{n} + (1-\beta) v_i^{n+1}$. The convexity of $G(\cdot)$ motivates the choice of an implicit discretization for $u$, i.e. $\alpha = 0$, because for $i = 1, \ldots, I$, there holds

Regarding the term in $uv$, only the case of GKS needs to be fixed and we thus require

It is natural to try and balance the terms by choosing a semi-explicit discretization with $\beta = \frac{1}{2}$, which yields

with the last term vanishing due to the symmetry of $K$.

However, implicit and explicit time discretizations for $v$ can also be considered at the expense of adding hypotheses on the kernel $K$. Indeed, for a given $0 \leq \beta \leq 1$, we find

As a consequence, an explicit (resp. implicit) scheme is suitable for the time discretization of $v$ provided that $K$ is a non-negative (resp. non-positive) symmetric kernel. Since $K$ is a non-negative symmetric kernel for the Generalized Keller-Segel equation 3, for simplicity we choose an explicit discretization for $v$.

Finally, we note that the interpolant does not play any role for energy discretization and we can use the simplest explicit or implicit discretization (both in $u$ and $v$), so as to make the analysis of the scheme as simple as possible.

5.1.   The gradient flow approach

We consider the full discretization of 18 and define the fully discrete flux in 21 as

At this level, we need to define the form of the interpolant $\varphi (u)_{i+1/2}^{n+1}$. From the theorem in Appendix A, we use an upwind technique in order to ensure well-posedness and monotonicity properties of the scheme. Thus, for $i = 1, \dots, I-1$, we define

Proposition 1 (Fully discrete gradient flow scheme). We assume either 10 and 13, or 11 and give the $u_i^0 \geq 0$. Then, the scheme 21-22-23 has the following properties:

(ⅰ) the solution $u_i^{n}$ exists and is unique, for all $i = 1, \dots, I$, and $n\geq 1$;

(ⅱ) it satisfies $u_i^{n} \geq 0$, and $u^n_i \leq M$ for the case 11, if it is initially true;

(ⅲ) the steady states $g(u_i) -v_i = \mu$ are preserved;

(ⅳ) the discrete energy dissipation inequality is satisfied

Notice that this theorem does not state a uniform bound in the case 10 and 13.

Proof. (ⅰ) We prove that the scheme satisfies the hypotheses of the theorem in Appendix A. We set

Then, the simplest case is when $\varphi$ satisfies 11, since clearly $u_i^{n+1} \equiv 0$ and $u_i^{n+1} \equiv M$ are respectively a sub- and supersolution. When $\varphi$ satisfies 10 and 13, $u_i^{n+1} \equiv 0$ is again a subsolution, while for the supersolution we choose $\bar{U}^{n+1}_i = g^{-1}(C +v_i^n)$. Such a choice indeed makes the flux terms vanish:

Thus $\bar{U}^{n+1}_i$ is a supersolution as soon as $g^{-1}(C + v_i^n) \geq u_i^n$, which holds when $C$ is taken to be large enough because we recall that assumption 13 ensures that $g(u)$ tends to $+\infty$ as $u$ tends to $+\infty$.

Moreover, the scheme is monotone since

where

so that $[x]_{+} \geq 0, \; [x]_{-} \leq 0$ for all $x$.

(ⅱ) Positivity of discrete solutions and the upper bound in the logistic case follow from the subsolution and supersolution built in step (ⅰ).

(ⅲ) Preservation of steady states at the discrete level follows immediately from the form we have chosen for the fully discrete fluxes.

(ⅳ) For the energy inequality, we remark that the contribution regarding time discretization is treated in the introduction of the present section. The space term is exactly treated as in the corresponding subsection of Section 4.

5.2.   The Scharfetter-Gummel approach

In 21, the fully discrete Scharfetter-Gummel flux reads

As for the gradient flow approach, we need the upwind technique to get a scheme which satisfies the hypotheses in Appendix A. So, we set for $i = 1, \dots, I-1$

Proposition 2 (Fully discrete Scharfetter-Gummel scheme). We assume either 10 and 13, or 11 and give the $u_i^0 \geq 0$. Then, the scheme 21-22-23 has the following properties:

(ⅰ) the solution $u_i^{n}$ exists and is unique, for all $i = 1, \dots, I$, and $n\geq 1$;

(ⅱ) it satisfies $u_i^{n} \geq 0$, and $u^n_i \leq M$ for the case 11, if it is initially true;

(ⅲ) the steady states $g(u_i) -v_i = \mu$ are preserved;

(ⅳ) the discrete energy dissipation inequality is satisfied

Proof. We argue exactly as for the gradient flow approach.

5.3.   The upwinding approach

The upwind scheme is driven by simplicity and, in 21, the fluxes are defined by

with

as in 23, but this time depending on the sign of $v_{i+1}^{n} - v_{i}^{n}$.

Proposition 3 (Fully discrete upwind scheme). We assume either 10 and 13, or 11 and give the $u_i^0 \geq 0$. Then, the scheme 21-22-23 has the following properties:

(i) the solution $u_i^{n}$ exists and is unique, for all $i = 1, \dots, I$, and $n\geq 1$;

(ii) it satisfies $u_i^{n} \geq 0$, and $u^n_i \leq M$ for the case 11, if it is initially true.

Proof. As for the gradient flow approach, the above choice makes the scheme monotone, because

Thus, arguing as for the gradient flow approach and relying on the results in Appendix A, existence and uniqueness of the discrete solution as well as preservation of the initial bounds follow immediately.

Thus, choice 24 enables to prove that the scheme is well-defined, satisfies $u_i^n \geq 0$ and preserves the bound $u_i^n \leq M$ for the case 11, but the energy dissipation inequality is lost. Also the steady states, in this case, are defined by the nonlinear relation $u_{i+1} - u_{i} = \varphi (u)_{i+1/2} (v_{i+1} - v_{i})$ which are usually not in the form 20.

6.   Numerical simulations

6.1.   The Fokker-Planck equation, $\varphi (u)=u$

We first present the numerical implementation of the Fokker-Planck equation with $\varphi (u) = u$. We here compare the Scharfetter-Gummel and upwind approaches (the gradient flow gives the same solution as the Scharfetter-Gummel and is thus not presented here). Both these schemes have error convergence of order $1$ in space, as it can be easily checked.

We consider a first case with $\chi/ D = 24$, with $I = 100$ and an initial density $u^0 = 1$. We take the velocity field as

In Figure 1, we compare the approximate stationary solutions obtained with the upwind scheme (blue, dashed line) and the Scharfetter-Gummel scheme (red line) with the exact stationary solution (black line), which in this case has the form $u(x) = C e^{\chi v(x)/D}$, with $C = \left( \int_{0}^{1} e^{\chi v(x)/D} dx \right)^{-1}$ so that the mass of the stationary solution is normalized. In this first case, the two schemes have no significant differences; this is a major difference with the Keller-Segel equation, as we show it in the next subsection.

6.2.   The nonlinear Keller-Segel equation

We turn to the equation 1 coupled with 3 for two nonlinear forms of the chemotactic sensitivity function: the logistic form $\varphi (u) = u (1-u)$ and the exponential form $\varphi (u) = u e^{-u}$. The goal is to compare the discrete solutions obtained with the three numerical approaches presented above when patterns arise, namely when Turing instabilities drive the formation of spatially inhomogeneous solutions (we refer to [21] for an introduction to this topic). To this end, we slightly modify the original equation 1 to

in order to emphasize the coefficients driving the instabilities: $D>0$, the constant diffusion coefficient and $\chi>0$, the chemosensitivity. The concentration of the chemoattractant $v$ remains driven by (3).

We first consider the logistic case with $\chi / D = 40$. We take as initial condition a random spatial perturbation of the constant steady state $u^0 = 0.5$ and solve the equation with 100 equidistant points in $[0, 1]$.

Figure 2 shows the evolution in time of the density $u_i^n$ obtained with the Scharfetter-Gummel (red line) and the gradient flow schemes (black, dashed line). After a rather short time, the initial spatial perturbation evolves, as expected, in spatially inhomogeneous patterns: a series of "steps" arise in the regions where the concentration of the chemoattractant is greater. After some time, a structure with a smaller number of steps forms when the two central plateaus merge. It is worth noticing that, even if the transitions from one structure to another happen very quickly, the time period during which these structures remain unchanged grows with the number of transitions that occurred. In [24], these intermediate patterns are called metastable, and this peculiar phenomenon is explained in details.

As for the schemes, Figure 2 shows that the Scharfetter-Gummel and the gradient flow approaches give the same solution; no difference can be spotted. This is not true for the upwind approach. In Figure 3, we compare the solutions to the Scharfetter-Gummel (red line) and the upwind (blue, dashed line) schemes. The upwind solution transitions faster from one metastable structure to the following than the Scharfetter-Gummel one. In fact, as proved above, the latter preserves discrete stationary profiles which, using the no-flux boundary conditions, solve the equation

From 26, it is clear that, in the logistic case, the expected stationary solutions are 0-1 plateaus (or "steps") connected by a sigmoid curve which is increasing or decreasing when $v$ is. We refer again to [24] and also to [12] for a detailed study of the stationary solutions and their properties for the logistic Keller-Segel system. In Figures 4a and 4b, we compare two stationary solutions to the Scharfetter-Gummel and upwind schemes, at time $t = 50$ and $t = 9000$ respectively. The Scharfetter-Gummel scheme approximates the 0-1 plateaus of metastable solutions better than the upwind scheme, whose solutions have smoother edges. We hypothesize that, since the Scharfetter-Gummel scheme preserves the metastable profiles better, it will also better preserve the time during which the solution will remain very close to a metastable state, up until the next transition. Consequently, this would mean that the upwind scheme accelerates the real dynamics.

Moreover, in Figure 4c we compare the $L^{\infty}$ dynamics of the three schemes, computing the quantity $||u^n- u^{n-1}||_{\infty}/||u^{n-1}||_{\infty}$ for each $n$. The peaks shown by this figure correspond to the transitions from one profile to another. Observe that, for both solutions of the Scharfetter-Gummel and the gradient flow scheme, the two peaks are further away in time than the ones from the upwind scheme: this confirms that the upwind solution is in advance when it comes to transitioning. Nevertheless, from $t \approx 6000$, the relative errors of the upwind solution are consistently greater that the ones from the two other approaches, thus confirming that only the Scharfetter-Gummel and the gradient flow schemes better preserve the exact discrete stationary profiles as well as the metastable ones. Also notice that none of the schemes produce overshoot, due to our upwinding of the term in $\psi(u)$.

Next, we consider an exponentially decreasing form of the chemotactic sensitivity function with $\chi / D = 24$. Again, we take as initial condition a random spatial perturbation of the constant steady state $u^0 = 0.7$ and solve the equation on 100 equidistant points. The evolution in time of discrete solutions obtained with the three numerical approaches are compared in Figures 5 and 7. In this model, no initial upper bound for the solution is imposed, so that the cells aggregate "naturally" where the chemoattractant has the greatest concentration, resulting in profiles without the plateaus observed in the logistic model. However, solutions face the same kind of transitions observed before, evolving from one stationary profile to another. As before, the Scharfetter-Gummel and the gradient flow approaches give the same solutions (Figure 5), while the solution of the upwind scheme evolves faster. In Figures Figures 6a and 6b, we compare stationary profiles obtained with the different approaches while in Figure 6c we compare dynamics of the solutions. This last figure shows that, as for the logistic model, smaller errors can be expected for the Scharfetter-Gummel and gradient flow approaches when steady states are reached.

7.   Conclusion

In the context of the Generalized Keller-Segel system, we have presented constructions of numerical schemes which extend previous works [10,20], built on two different views of energy dissipation. Our construction unifies these two views, the gradient flow and Scharfetter-Gummel symmetrizations. Our schemes preserve desirable continuous properties: positivity, mass conservation, exact energy dissipation, discrete steady states. Being correctly tuned between implicit and explicit discretization, they can handle large time steps without CFL condition.

The present work is motivated by experiments of breast cancer cells put in a 3D structure mimicking the conditions they meet in vivo, namely in the extracellular matrix. After a few days, images of 2D sections show that cells have organized as spheroids, a phenomenon believed to be driven by chemotaxis. The spheroids can then be interpreted as Turing patterns for Keller-Segel type models and it is crucial to use appropriate schemes for them to be distinguishable from actual steady states or numerical artifacts. Comparing 2D simulations of such models with these experimental images will be the subject of future work.

In fact, it is important to remark that the schemes we presented here in 1D could be easily extended to rectangular domains, without loss of properties 4-7. However, it remains a perspective to treat more general geometries in a multi-dimensional setting with our approach.

Acknowledgments

The authors acknowledge partial funding from the "ANR blanche" project Kibord [ANR-13-BS01-0004] funded by the French Ministry of Research.

B.P. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 740623).

Appendix A. Existence and uniqueness for monotone schemes

We recall sufficient conditions for which an implicit Euler discretization in time can be solved, independently of the step-size. This is the case for a monotone scheme. The proof relies on the existence of sub- and supersolutions, and thus also yields the preservation of positivity and other pertinent bounds as we have used in Section 5.

We consider the problem of finding a unique solution $\left(u_i^{n+1}\right)$ to the nonlinear equation arising in Section 5 which reads

We write a general proof for a scheme of the form

where we consider the problem of finding a solution $(u_i)$ (which stands for $u_i^{n+1}$).

Here we assume that the $f_i$ are given (it stands for $u^n_i$) and that the $A_{i+ \frac 12}$ are Lipschitz continuous and, a.e.,

and there are a supersolution $(\bar U_i)_{ i = 1\; \ldots , I}$ and a subsolution $(\underline U_i)_{ i = 1\; \ldots , I}$ such that for all $i = 1, \ldots, I,$

We build a solution of 28 using an evolution equation

Theorem A.1. Assume 29 and the existence of a subsolution and of a supersolution. Then,

(ⅰ) For a supersolution (resp. subsolution) initial data, the dynamics 32 satisfies $\frac {d \bar u_i(t)}{dt} \leq 0$ (resp. $\frac {d \underline u_i(t)}{dt} \geq 0$) for all times $t\geq 0$, and thus $\bar u_i(t)$ is a supersolution (resp. subsolution) for all times.

(ⅱ) A subsolution is smaller than a supersolution.

(ⅲ) $\bar u_i(t)$ and $\underline u_i(t)$ converge to the same solution of 28.

Proof. (ⅰ) We prove the statement with the supersolution. We set

Since the $A_{i+ \frac 12}$ are Lipschitz continuous, from equation 32 we deduce that the quantities $\frac {d \bar u_i(t)}{dt}$ are also Lipschitz continuous. From Rademacher's Theorem, the $z_i$ are differentiable a.e. and their a.e. derivatives are also their distributional derivatives.

Differentiating equation 32, we obtain for $i = 1, \ldots I$, and for a.e. $t>0$

The solution cannot change sign and thus for $i = 1, \ldots , I$, $z_i(t) \leq 0$ for all times.

(ⅱ) Consider $\underline u, \; \bar u$ sub (super) solutions. Set $w = \underline u- \bar u$ and we want to prove that $w\leq 0$.

We write for $i = 1, \ldots , I$

For $i = 1, \ldots , I$, we multiply by ${\rm sgn}_+(w_i): = 1_{w_i >0}$ and add the relations to conclude that

with

For each of the these terms, we show that $J_{i+ \frac 12} \geq 0$, $K_{i+ \frac 12} \geq 0$, as follows. Only the case when the + signs in the right brackets are different has to be considered. Assume for instance that

Then, we have by assumption 29,

Therefore $\sum_{i = 1}^{I} (w_i)_+ \leq 0$ and this implies $w_i \leq 0$ for all $i$. From (ⅰ) and (ⅱ) we infer that the subsolution increases and is bounded from above by the supersolution, so that it exists for all times. Similarly, the supersolution exists for all times and thus, so does the solution $u_i(t)$ to 32, so that we can speak of its convergence.

(ⅲ) This is clear since the limits are solutions.