A convergent finite-volume scheme for nonlocal cross-diffusion systems for multi-species populations

An implicit Euler finite-volume scheme for a nonlocal cross-diffusion system on the one-dimensional torus, arising in population dynamics, is proposed and analyzed. The kernels are assumed to be in detailed balance and satisfy a weak cross-diffusion condition. The latter condition allows for negative off-diagonal coefficients and for kernels defined by an indicator function. The scheme preserves the nonnegativity of the densities, conservation of mass, and production of the Boltzmann and Rao entropies. The key idea is to ``translate'' the entropy calculations for the continuous equations to the finite-volume scheme, in particular to design discretizations of the mobilities, which guarantee a discrete chain rule even in the presence of nonlocal terms. Based on this idea, the existence of finite-volume solutions and the convergence of the scheme are proven. As a by-product, we deduce the existence of weak solutions to the continuous cross-diffusion system. Finally, we present some numerical experiments illustrating the behavior of the solutions to the nonlocal and associated local models.

where   are some constants.The kernel functions   : T → R are periodically extended to R, and  = ( 1 , . . .,   ) is the solution vector.If we define   =  0 , where  ∈ {1, . . ., } and  0 is the Dirac measure, we can rewrite   as   () = -We prove the existence of solutions to the finite-volume scheme, which are nonnegative componentwise, conserve the discrete mass, and satisfy discrete versions of the entropy inequalities ( 6) and (7).-We show that the discrete solutions converge to a weak solution to (1)-(3) when the mesh size tends to zero.
As a by-product, this proves the existence of a weak solution to (1)-( 2).-We illustrate numerically the rate of convergence (in space) in the   -norm as well as the rate of convergence in different metrics of the solution to the nonlocal system towards the solution of the local one (localization limit).Moreover, we illustrate the segregation phenomenon exhibited by the solutions to (1)- (3); see [4].
The paper is organized as follows.The numerical scheme and our main results are introduced in Section 2. We prove the existence of discrete solutions in Section 3, while the proof of the convergence of the scheme is presented in Section 4. In Section 5, numerical experiments are given, Appendix A contains some auxiliary results, and we show in Appendix B that indicator kernels generally do not fulfill inequality (5).
We recall the definition of the space BV(T) of functions of bounded variation.A function  ∈  1 (T) belongs to BV(T) if its total variation TV(), given by TV() = sup
We consider the one-dimensional equations mainly for notational simplicity.In several space dimensions  > 1, we infer uniform estimates in spaces with weaker integrability than in one space dimension, because of Sobolev embeddings.Thanks to the positive definiteness condition on   ℓ−ℓ ′ , we obtain a bound for   in the discrete  2 (0,  ;  1 (T)) norm, which allows us to conclude, together with the Rao entropy estimate, by the discrete Gagliardo-Nirenberg inequality, a bound for   in  2 (  ), which is sufficient to estimate the product       ().In the one-dimensional situation, this procedure simplifies; see Lemma 4.5.We discuss the multidimensional case in Remark 4.7.
Our results also hold if  = 0, since the condition  > 0 provides an estimate for   in the discrete norm of  2 (0,  ;  1,1 (T)), while the positive definiteness condition on   ℓ−ℓ ′ allows us to conclude a stronger bound in the discrete norm of  2 (0,  ;  1 (T)).Notice that kernels of the type   = 1  satisfy Hypothesis (H3) (for suitable   and   ).
Condition  0 ∈  2 (T; R  ) in Hypothesis (H2) is needed to obtain a finite initial Rao entropy   ( 0 ).For the existence result, the assumption on the kernels can be weakened to   ∈  1 (T).The boundedness condition on   in Hypothesis (H3) is needed in the proof of the convergence of the scheme.
This theorem is proved by solving a fixed-point problem based on a topological degree argument, similar as in [18].For this, we formulate (11) in terms of the entropy variable   =   log   and regularize the equations by adding the discrete analog of −∆  +   .The regularization ensures the coercivity in the variable   .After transforming back to the original variable   = exp(  /  ), we obtain automatically the positivity of   (and nonnegativity after passing to the limit  → 0).Like on the continuous level, the derivation of the discrete entropy inequalities (20) and (21) relies on the detailed-balance condition     =     for all ,  = 1, . . ., .
For our second main result, we need to introduce some notation.We define the "diamond" cell of the dual mesh  ℓ+1/2 = ( ℓ ,  ℓ+1 ) with center  ℓ+1/2 .These cells define another partition of T. The gradient of  ∈   is then defined by We also introduce a sequence of space-time discretizations (  ) ∈N indexed by the mesh size The proof of Theorem 2.6 is based on suitable estimates uniform with respect to ∆  and ∆  , derived from the discrete entropy inequalites.A discrete version of the Aubin-Lions lemma from [15] yields the strong convergence of a subsequence of solutions (  ) to ( 11)- (13).The most technical part is the identification of the limit function as a weak solution to (1)-(2).

Proof of Theorem 2.5
Theorem 2.5 is proved by induction over  = 1, . . .,   .We first regularize the problem and prove the existence of an approximate solution by using a topological degree argument for the fixed-point problem.The discrete entropy inequalities yield a priori estimates independent of the approximation parameter.The de-regularization limit is performed thanks to the Bolzano-Weierstraß theorem.

Continuity of 𝐹
We fix  ∈ {1, . . ., }, multiply (22) by   ,ℓ , and sum over ℓ ∈ : The left-hand side can be rewritten by using discrete integration by parts (or summation by parts): The first term on the right-hand side of ( 23) is estimated by the Cauchy-Schwarz inequality, taking into account that  ∈   , which implies a finite discrete  2 (T) norm for  ,ℓ = exp( ,ℓ /  ): where here and in the following  > 0, (∆, ) > 0, etc. are generic constants with values changing from line to line.We split the second term on the right-hand side of (23) into two parts: where For  1 , we use discrete integration by parts, the Cauchy-Schwarz inequality, and the fact that  ∈   : Using discrete integration by parts, and definition (13) of  ,ℓ , we obtain where For  21 , because of the bound in   , we can estimate  ,ℓ+1/2 ≤ max{ ,ℓ+1 ,  ,ℓ } ≤ ().Then, thanks to the Cauchy-Schwarz inequality, we obtain For  22 , applying the discrete analog ( 17) of the rule where we used the notation of Section 2.1.Similarly to  21 , we infer that Then, by the Cauchy-Schwarz inequality and the discrete convolution inequality from Lemma A.3 in Appendix A, Combining these estimates, we deduce from ( 23) that ‖   ‖ 1,2, ≤ (∆, ).

Existence of a fixed point
We show that  :   → R  admits a fixed point by using a topological degree argument.We recall that the Brouwer topological degree is a mapping deg :  → Z, where see Chapter 1, Theorem 3.1 from [11] for details and properties.If we show that any solution (  , ) for sufficiently large values of  > 0, then we deduce from the invariance by homotopy that deg This implies that there exists   ∈   such that ( −  )(  ) = 0, which is the desired fixed point.Let (  , ) be a fixed point of   =  (  ).If  = 0, there is nothing to show.Therefore, let  > 0. Then for all ℓ ∈  and  = 1, . . ., , where   ,ℓ = exp(  ,ℓ /  ), and the fluxes ℱ  ,ℓ±1/2 are defined as in (12) with   ,ℓ replaced by   ,ℓ .We multiply the previous equation by ∆  ,ℓ , sum over ℓ ∈ ,  = 1, . . ., , and use discrete integration by parts as in (24): For the first term on the right-hand side, we use   ,ℓ =   log   ,ℓ and the convexity of ℎ() = (log  − 1): )︀ .
Recalling definition (18) of ℋ  , this shows that Like in Section 3.2, we split the second term in (26) into two parts: We use discrete integration by parts, the definition   ,ℓ =   log   ,ℓ , and the elementary inequality for ,  > 0 to estimate the first term: For the second term  4 , we use discrete integration by parts and   ,ℓ =   log   ,ℓ again as well as property ( 14) (discrete chain rule): Then, inserting definition (3) of   ,ℓ and using the discrete analog (17) of where We insert ( − 1) −1 ∑︀ ̸ = 1 = 1 and ∑︀ ℓ ′ ∈ ∆ = 1 (note that m(T) = 1) in  41 and split the resulting sum into two parts: We exchange  and  as well as ℓ and ℓ ′ in the second term, which leads to Similarly, we distinguish between  <  and  >  in  42 and exchange  and  as well as ℓ and ℓ ′ in the sum over  > , leading to By Remark 2.2, we have The sum of  41 and  42 can be written as a quadratic form in D ℓ    and D ℓ ′    with the matrix   ℓ−ℓ ′ , defined in (19).This shows that Collecting the estimates for  3 and  4 in (27), we deduce from (26) the following regularized discrete entropy inequality: We proceed with the topological degree argument.We set  = 1 + (ℋ  ( −1 )/(∆)) 1/2 .Then (28) implies that and hence   ̸ ∈   .We infer that deg( − ,   , 0) = 1 and consequently,  admits a fixed point.Note that we did not use the estimate for    in the seminorm | • | 1,2, at this point, such that  = 0 is admissible here (and also in the following two subsections).

Discrete Rao entropy inequality
We prove inequality (21).To this end, we multiply (11) by ∆    ,ℓ and sum over ℓ ∈ ,  = 1, . . ., : For the first term in (29), we use the definition of   ,ℓ : where We rewrite  5 and  6 according to Combining the second terms in  5 and  6 , using similar computations as for  4 in Section 3.3, and applying Hypothesis (H3) shows that the second term of  5 +  6 is nonnegative leading to Then it holds that Now, we split the second term in (29) again into two parts: where We reformulate  7 by using discrete integration by parts: Then, with similar computations as for  4 in Section 3.3, we obtain Finally, the term  8 can be rewritten as Hence, we infer from (29) that which proves (21).Finally, conservation of mass follows from summing (11) over ℓ ∈  and observing that the sum over the numerical fluxes vanishes.This ends the proof of Theorem 2.5.

Proof of Theorem 2.6
To prove the convergence of the scheme, we derive first some uniform estimates and then apply a discrete Aubin-Lions compactness lemma.

Uniform estimates
Let (  ) ∈N be a sequence of finite-volume solutions to (11)-( 13) associated to the mesh   and constructed in Theorem 2.5.The conservation of mass and the discrete entropy inequalities (20) and (21) show that, after summing over  = 1, . . .,    , max where  > 0 denotes here and in the following a constant independent of the mesh size   = max{∆  , ∆  }, but possibly depending on  0 and  .Because of the positive definiteness of   ℓ−ℓ ′ , we conclude a bound for    in the norm ‖ • ‖ 1,2, .
Lemma 4.1.Let the assumptions of Theorem 2.6 hold.Then there exists  > 0 independent of   (but depending on the positive definiteness constant   ) such that for all  ∈ N,  = 1, . . ., , Proof.We infer from (20) that Together with the first bound in (30), this finishes the proof.
Lemma 4.2.Let the assumptions of Theorem 2.6 hold.Then there exists a constant  > 0 independent of   (but depending on ) such that for all  ∈ N,  = 1, . . ., , Moreover, there exists another constant, still denoted by  > 0 and independent of   , such that Proof.As m(T) = 1, thanks to the Cauchy-Schwarz inequality, Using (31), this shows that To show the discrete  ∞ (T) bound, we apply the continuity of the embedding BV(T) ˓→  ∞ (T) (in one space dimension).We conclude that, for  = 1, . . ., , For the last part, we estimate as follows: Then we deduce from the elementary inequality ( Summing over , we infer that where we used Lemma 4.2 for the last inequality.At this point, we need the discrete  2 (0,  ;  1 (T)) bound of ( , ).This ends the proof.
Next, we show a uniform bound for the discrete time derivative.

Compactness
We claim that the estimates from Lemmas 4.2 and 4.3 are sufficient to conclude the relative compactness of (  ) ∈N .In fact, the result follows from the discrete Aubin-Lions lemma Theorem 3.4 from [15] if the following two properties are satisfied: -For any (  ) ∈N ⊂   such that sup ∈N ‖  ‖ 1,2, ≤  for some  > 0, there exists  ∈  2 (T) satisfing, up to a subsequence,   →  strongly in  2 (T).This property follows from Theorem 14.1 from [13].-If   →  strongly in  2 (T) and ‖  ‖ −1,2, → 0 as  → ∞, then  = 0.This property can be replaced by the condition that ‖ • ‖ 1,2, and ‖ • ‖ −1,2, are dual norms with respect to the  2 (T) norm, which is the case Remark 6 from [15].A more detailed proof can be found in Proposition 10 from [18].
Let us now adapt the Gagliardo-Nirenberg inequality to our situation.Let  = 1, . . .,    be fixed.We first apply Lemma A.4 with  =  = 2: Proof.We follow the strategy of Corollary 14 from [17].First, we rewrite   ,ℓ defined in (13).By a change of variables, we have Since we know that   −  , → 0 strongly in  2 (  ), it is sufficient to prove that   *   −    → 0 strongly in  2 (  ).For this, we write By Young's convolution inequality, we have Setting (, ) =   ( − ) −   ( ℓ − ) for  ∈  ℓ and  ∈ T, we estimate Since ( , ) is bounded in  2 (  ), it remains to verify that the first factor converges to zero as ∆  → 0. This follows from the density of continuous functions in  2 (T).Indeed, let  > 0 and    be continuous such that The last term is smaller than  if we choose ∆  sufficiently small.We have shown that sup ||≤Δ ‖  ( + •) −   ‖ 2  2 (T) → 0 as  → ∞ and   *   −    → 0 strongly in  2 (  ).This proves the first part of the lemma.
Thanks to (32), we have shown that (    , ) ∈N is bounded in  2 (  ).Hence, up to a subsequence,     , ⇀  weakly in  2 (  ).The first part of the proof shows that  =     (), finishing the proof.
It remains to prove that |  30 −   3 | → 0. First, using discrete integration by parts, we rewrite   3 as well as   30 as Then we find that Thanks to the regularity of   , there exists a constant  independent of   such that We obtain a similar expression if we integrate     over ( ℓ+1/2 ,  ℓ+1 ).Thus, since we have It follows for  ∈ {1, . . ., } with  ̸ = , using the discrete analog (17) of At this point, we need the regularity condition   ∈  ∞ (T) from Hypothesis (H3).Hence, it holds that It remains to apply the Cauchy-Schwarz inequality to conclude that Finally, we infer from Lemma 4.2 that |  3 −   30 | → 0 as  → ∞.Here, we need the discrete  2 (0,  ;  1 (T)) bound for   , which follows if   > 0. This concludes the proof of Lemma 4.6.

Numerical experiments
In this section, we present several numerical experiments to illustrate the behavior of the scheme.The scheme was implemented in one space dimension using Matlab.In all the subsequent numerical tests, we choose the upwind mobility (15).The code is available at https://gitlab.tuwien.ac.at/asc/nonlocal-crossdiff.Our code is an adaptation of that one developed in [17] for the approximation of the nonlocal SKT system.We refer the reader to Section 6.1 of [17] for a complete presentation of the different methods used to implement the scheme.
5.1.Test case 1. Rate of convergence in space for various   -norms, convolution kernels, and initial data We investigate the rate of convergence in space of the scheme at final time  = 1.In all test cases of this section, we consider  = 2 species,  = 10 −4 , the coefficient matrix  = (  ) 1≤,≤2 given by  = (︂ 0.1251 0.25 1 2 )︂ , and  1 = 4,  2 = 1.We consider various initial data and kernels.More precisely, we choose 0 1 () = cos (2) + 1,  0 2 () = sin (2 − /2) + 1, (41) and the kernels First, we consider a mesh of  init = 32 cells and the time step size ∆ init = 1/64.Then, starting from this initial mesh, we refine the mesh in space by doubling the number of cells and halving the time step size, i.e.  new = 2 old and ∆ new = ∆ old /2.This refinement of the meshes is in agreement with the first-order convergence rate of the Euler discretization in time and the expected first-order convergence rate in space of the scheme, due to the choice of the upwind mobility in the numerical fluxes.As exact solutions to system (1)-( 3) are not explicitly known, we refine the mesh in space and time until  end = 2048 and ∆ end = 1/4096, and we consider the solutions of the scheme obtained for  end and ∆ end as reference solutions.The error is computed between the reference solutions and the solutions obtained for  = 1024 cells and ∆ = 1/2048 at final time  = 1.Finally, using linear regression in logarithmic scale, we present in Table 1 the experimental order of convergence in the  1 and  ∞ -norms.As expected, we observe a rate of convergence around one.In Table 1, the numbers in bold letters denote the number of the test case available in our code (see the file loadTestcase.m).

Test case 2. Rate of convergence of the localization limit in various metrics
In the second test case, following [17], we evaluate numerically the rate of convergence of the localization limit.More precisely, for some sequences of kernels converging towards the Dirac measure  0 , we compute the rate of convergence in different metrics of the solutions to scheme (10)-( 13) towards its local version, i.e.   =  0 for all ,  = 1, . . ., .At the continuous level, one can show, by adapting the approach of [19], that the localization limit holds thanks to a compactness method; see also [12] for the SKT system.However, so far no explicit rate of convergence is available.The goal of this numerical test is to obtain a better insight into this rate of convergence.Besides, it also illustrates Remark 2.4.
We consider the following parameters (for all 6 test cases of this section):  = 3 species, diffusion parameter  = 10 −4 , coefficient matrix and  1 = 4,  2 = 2,  3 = 2.We choose the final time  = 1, a mesh of  = 512 cells, and the time step size ∆ = 10 −3 .Furthermore, we take the nonsmooth initial data and the smooth initial data The kernels are chosen according to In our experiments, starting from  init = 2 7 ∆, we successively halve  until we reach  = ∆.For each value of , we compute the solutions to the nonlocal scheme (10)-( 13) at final time.We evaluate the  1 ,  ∞ , and Wasserstein distance  1 between the solution to the nonlocal scheme and the solution to the local one (for this, it is enough to set  = 0 in our code).Since we work in one space dimension, we can explicitly compute the Wasserstein distance  1 ; see Chapter 2 from [24].The rates of convergence are estimated by linear regression (in log scale) and the results are presented in Table 2. Surprisingly, we observe a slightly better rate of convergence in the case of nonsmooth initial data.As before, the names in bold letters in Table 2 denote the name of the test case available in our code (see the file loadTestcase.m).

Test case 3. Segregation phenomenon
In this numerical experiment, we set  = 0.Under the assumptions  = 2 species,   = 1, and   =  0 for ,  = 1, 2, it has been shown in [4] that if the initial data are segregated (initial data with disjoint supports) then the solutions remain segregated (i.e., they have disjoint supports) for all time.The main goal of this subsection is to illustrate the segregation pattern due to the nonlocal terms, i.e.   ̸ =  0 .We expect that the solutions to the nonlocal model, given segregated initial data, are completely segregated, and that there exists a small region, i.e. a "gap" between the supports of the species, with a size that is related to the radius of the interaction kernels.Let us notice that in the subsequent test cases, Hypothesis (H3) is never satisfied.However, we did not encounter any numerical issues with our code.
We launched the code for a mesh of 512 cells and the time step size ∆ = 10 −4 .In the case of  = 2 species, we considered the initial data
In Figures 1 and 2, we present the segregation pattern at time  = 0.02 and  = 0.2 obtained for the local model,   =  0 , and the nonlocal model with For small times, the support of the species extends until reaching the support of another species.In the local model, the species slightly mix (due to numerical diffusion), while we observe a "gap" between the supports of the solutions in the nonlocal model.This "gap" is of order 0.1 which is the size of the radius of the kernels   .Similar numerical results have been observed in Section 6 of [8] but using different kernel functions and two species only.

Test case 4. Dissipation of entropy
In the last numerical experiment, we plot the two entropies ℋ  (()) and ℋ  (()) over time in semilogarithmic scale to illustrate the entropy production as proved in Theorem 2.5.We set  = 1.5, ∆ = 10 −4 , use a mesh of  = 512 cells, and choose  = 2 species.The remaining parameters are taken as in Section 5.1; see Table 1 and the test cases therein.As expected, the entropies are decreasing functions of time.The Rao entropy decays first quickly but then stabilizes slowly, while the Boltzmann entropy takes more time to stabilize.

Figure 1 .
Figure 1.Comparison of the segregation pattern for two species at times  = 0.02 (top) and  = 0.2 (bottom) obtained from the local model (left) and nonlocal model (right).The solutions are almost in the steady state at  = 0.2.

Figure 2 .
Figure 2. Comparison of the segregation patterns for three species at times  = 0.02 (top) and  = 0.2 (bottom) obtained from the local model (left) and nonlocal model (right).The solutions are almost in the steady state at  = 0.2.

Figure 3 .
Figure 3. Temporal decay of the Boltzmann and Rao entropies for test cases 15 (left) and 16 (right) in semi-logarithmic scale.
, →   holds in   (  ) for every  < 2/( − 1) (instead of  < 6 in the one-dimensional case) and in particular in  2 (  ).Thus, the statement of Lemma 4.4 holds, and we have ∇   , ⇀ ∇  () weakly in  2 (  ), where ∇  denotes the discrete gradient.In particular,  , ∇   , ⇀   ∇  () weakly in  4/3 (  ) as in the one-dimensional case.From this point on, the convergence of the scheme follows the lines of Section 4.3.⊓ ⊔

Table 1 .
Orders of convergence in the  1 and  ∞ norms in space at final time  = 1 for different kernels and initial data.

Table 2 .
Rates of convergence of the localization limit in the  1 ,  ∞ and  1 metric for different initial data and kernels.