Convergence of a Finite Volume Scheme for a System of Interacting Species with Cross-Diffusion
aa r X i v : . [ m a t h . NA ] F e b CONVERGENCE OF A FINITE VOLUME SCHEME FOR A SYSTEM OFINTERACTING SPECIES WITH CROSS-DIFFUSION
Jos´e A. Carrillo
Department of Mathematics, Imperial College LondonLondon SW7 2AZ, United Kingdom
Francis Filbet
Institut de Math´ematiques de Toulouse , Universit´e Paul SabatierToulouse, France
Markus Schmidtchen
Department of Mathematics, Imperial College LondonLondon SW7 2AZ, United Kingdom
Abstract.
In this work we present the convergence of a positivity preserving semi-discrete finitevolume scheme for a coupled system of two non-local partial differential equations with cross-diffusion. The key to proving the convergence result is to establish positivity in order to obtaina discrete energy estimate to obtain compactness. We numerically observe the convergence toreference solutions with a first order accuracy in space. Moreover we recover segregated stationarystates in spite of the regularising effect of the self-diffusion. However, if the self-diffusion or thecross-diffusion is strong enough, mixing occurs while both densities remain continuous. Introduction
In this paper we develop and analyse a numerical scheme for the following non-local interactionsystem with cross-diffusion and self-diffusion ∂ρ∂t = ∂∂x (cid:18) ρ ∂∂x ( W ⋆ ρ + W ⋆ η + ν ( ρ + η )) + ǫ ∂ρ ∂x (cid:19) ,∂η∂t = ∂∂x (cid:18) η ∂∂x ( W ⋆ η + W ⋆ ρ + ν ( ρ + η )) + ǫ ∂η ∂x (cid:19) , (1)governing the evolution of two species ρ and η on an interval ( a, b ) ⊂ R for t ∈ [0 , T ). The systemis equipped with nonnegative initial data ρ , η ∈ L ( a, b ) ∩ L ∞ + ( a, b ). We denote by m the massof ρ and by m the mass of η , respectively, m = Z ba ρ ( x ) d x, and m = Z ba η ( x ) d x. Mathematics Subject Classification.
Primary: 74S10; 65M12; 92C15; Secondary: 45K05; 92D25; 47N60.
Key words and phrases.
Finite volume methods; Integro-partial differential equations; Population dynamics (gen-eral); Developmental biology, pattern formation;
On the boundary x = a and b , we prescribe no-flux boundary conditions ρ ∂∂x ( W ⋆ ρ + W ⋆ η + ν ( ρ + η ) + ǫρ ) = 0 ,η ∂∂x ( W ⋆ η + W ⋆ ρ + ν ( ρ + η ) + ǫη ) = 0 , such that the total mass of each species is conserved with respect to time t ≥
0. While the self-interaction potentials W , W ∈ C b ( R ) model the interactions among individuals of the samespecies (also referred to as intraspecific interactions), the cross-interaction potentials W , W ∈ C b ( R ) encode the interactions between individuals belonging to different species, i.e. interspecific interactions. Here C b ( R ) denotes the set of twice continuously differentiable functions on R withbounded derivatives. Notice that the convolutions W ij ⋆ ψ , with ψ a density function defined on[ a, b ], are defined by extending the density ψ by zero outside the interval [ a, b ]. The two positiveparameters ǫ, ν > th century biologists found that thedispersal rate of certain insects depend on the density itself, leading to the nonlinear diffusion termswe incorporated in the model, cf. [36, 37, 34, 13, 32]. At the same time we would like to stress thatthe self-diffusion terms are relevant for the convergence analysis below.It is the interplay between the non-local interactions of both species and their individual and jointsize-exclusion, modelled by the non-linear diffusion [5, 4, 6, 40, 12, 10], that leads to a large varietyof behaviours including complete phase separation or mixing of both densities in both stationaryconfigurations and travelling pulses [11, 19].While their single species counterparts have been studied quite intensively [35, 21, 41, 15] andreferences therein, related two-species models like the system of our interest, Eq. (1), have onlyrecently gained considerable attention [25, 11, 19, 24, 17]. One of the most striking phenomenaof these interaction models with cross-diffusion is the possibility of phase separation. Since theseminal papers [33, 5] established segregation effects for the first time for the purely diffusive systemcorresponding to (1) for W ij ≡ i, j ∈ { , } and ǫ = 0, many generalisations were presented. Theyinclude reaction-(cross-)diffusion systems [3, 7, 17] and references therein, and by adding non-localinteractions [11, 19, 24, 2] and references therein. Ref. [24] have established the existence of weaksolutions to a class of non-local systems under a strong coercivity assumption on the cross-diffusionalso satisfied by system (1).Typical applications of these non-local models comprise many biological contexts such cell-celladhesion [39, 38, 16], for instance, as well as tumour models [31, 26], but also the formation of thecharacteristic stripe patterns of zebrafish can be modelled by these non-local models [42]. Systemsof this kind are truly ubiquitous in nature and we remark that ‘species’ may not only refer tobiological species but also to a much wider class of (possibly inanimate) agents such as planets,physical or chemical particles, just to name a few.Since system (1) is in conservative form a finite volume scheme is a natural choice as a numericalmethod. This is owing to the fact that, by construction, finite volume schemes are locally conser-vative: due the divergence theorem the change in density on a test cell has to equal the sum ofthe in-flux and the out-flux of the same cell. There is a huge literature on finite volume schemes,first and foremost [28]. Therein, the authors give a detailed description of the construction of suchmethods and address convergence issues. Schemes similar to the one proposed in Section 2 have ONVERGENCE OF A FINITE VOLUME SCHEME 3 been studied in [9] in the case of nonlinear degenerate diffusion equations in any dimension. Asimilar scheme for a system of two coupled PDEs was proposed in [22]. Later, the authors in [14]generalised the scheme proposed in [9] including both local and non-local drifts. The scheme wasthen extended to two species in [19]. All the aforementioned schemes have in common that theypreserve nonnegativity – a property that is also crucial for our analysis.Before we define the finite volume scheme we shall present a formal energy estimate for the con-tinuous system. The main difficulty in this paper is to establish positivity and reproducing thecontinuous energy estimate at the discrete level. The remainder of the introduction is dedicated topresenting the aforementioned energy estimate. Let us considerdd t Z ba ρ log ρ d x = Z ba log ρ ∂ρ∂t d x = Z ba log ρ ∂∂x (cid:18) ρ ∂∂x ( W ⋆ ρ + W ⋆ η + ν ( ρ + η ) + ǫρ ) (cid:19) d x = − Z ba ρ ∂∂x ( W ⋆ ρ + W ⋆ η + ν ( ρ + η ) + ǫρ ) ∂∂x (log ρ ) d x. Upon rearranging we getdd t Z ba ρ log ρ d x + ν Z ba ∂∂x ( ρ + η ) ∂ρ∂x d x + ǫ Z ba (cid:12)(cid:12)(cid:12)(cid:12) ∂ρ∂x (cid:12)(cid:12)(cid:12)(cid:12) d x = Z ba ( W ′ ⋆ ρ + W ′ ⋆ η ) ∂ρ∂x d x. A similar computation for η yieldsdd t Z ba η log η d x + ν Z ba ∂∂x ( ρ + η ) ∂η∂x d x + ǫ Z ba (cid:12)(cid:12)(cid:12)(cid:12) ∂η∂x (cid:12)(cid:12)(cid:12)(cid:12) d x = Z ba ( W ′ ⋆ η + W ′ ⋆ ρ ) ∂η∂x d x, whence, upon adding both, we obtaindd t Z ba [ ρ log ρ + η log η ] d x + ν Z ba (cid:12)(cid:12)(cid:12)(cid:12) ∂σ∂x (cid:12)(cid:12)(cid:12)(cid:12) d x + ǫ Z ba (cid:12)(cid:12)(cid:12)(cid:12) ∂ρ∂x (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) ∂η∂x (cid:12)(cid:12)(cid:12)(cid:12) ! d x = D ρ + D η , where σ = ρ + η and D ρ := Z ba ( W ′ ⋆ ρ + W ′ ⋆ η ) ∂ρ∂x d x, D η := Z ba ( W ′ ⋆ η + W ′ ⋆ ρ ) ∂η∂x d x, denote the advective parts associated to ρ and η , respectively. The advective parts can be controlledby using the weighted Young inequality to get |D ρ | = (cid:12)(cid:12)(cid:12)(cid:12)Z ba ( W ′ ⋆ ρ + W ′ ⋆ η ) ∂ρ∂x d x (cid:12)(cid:12)(cid:12)(cid:12) ≤ α Z ba | W ′ ⋆ ρ + W ′ ⋆ η | d x + α Z ba (cid:12)(cid:12)(cid:12)(cid:12) ∂ρ∂x (cid:12)(cid:12)(cid:12)(cid:12) d x, JOS´E A. CARRILLO, FRANCIS FILBET, AND MARKUS SCHMIDTCHEN for some α >
0. In choosing 0 < α < ǫ we obtaindd t Z ba [ ρ log ρ + η log η ] d x + ν Z ba (cid:12)(cid:12)(cid:12)(cid:12) ∂σ∂x (cid:12)(cid:12)(cid:12)(cid:12) d x + (cid:16) ǫ − α (cid:17) Z ba (cid:12)(cid:12)(cid:12)(cid:12) ∂ρ∂x (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) ∂η∂x (cid:12)(cid:12)(cid:12)(cid:12) ! d x ≤ C ρ + C η α , (2)where C ρ = k W ′ ⋆ ρ + W ′ ⋆ η k L and C η = k W ′ ⋆ η + W ′ ⋆ ρ k L . From the last line, Eq. (2), wemay deduce bounds on the gradient of each species as well as on their sum. As mentioned abovethe crucial ingredient for this estimate is the positivity of solutions.The rest of this paper is organised as follows. In the subsequent section we present a semi-discretefinite volume approximation of system (1) and we present the main result, Theorem 2.4. Section3 is dedicated to establishing positivity and to the derivation of a priori estimates. In Section 4we obtain compactness, pass to the limit, and identify the limiting functions as weak solutionsto system (1). We conclude the paper with a numerical exploration in Section 6. We study thenumerical order of accuracy and discuss stationary states and phase segregation phenomena.2. Numerical scheme and main result
In this section we introduce the semi-discrete finite volume scheme for system (1). To begin with,let us introduce our notion of weak solutions.
Definition 2.1 (Weak solutions.) . A couple of functions ( ρ, η ) ∈ L (0 , T ; H ( a, b )) is a weaksolution to system (1) if it satisfies − Z ba ρ ϕ (0 , · ) d x = Z T Z ba (cid:20) ρ (cid:18) ∂ϕ∂t + (cid:18) − ν ∂σ∂x + ∂V ∂x (cid:19) ∂ϕ∂x (cid:19) + ǫ ρ ∂ ϕ∂x (cid:21) d x d t, (3a) and − Z ba η ϕ (0 , · ) d x = Z T Z ba (cid:20) η (cid:18) ∂ϕ∂t + (cid:18) − ν ∂σ∂x + ∂V ∂x (cid:19) ∂ϕ∂x (cid:19) + ǫ η ∂ ϕ∂x (cid:21) d x d t, (3b) respectively, for any ϕ ∈ C ∞ c ([0 , T ) × ( a, b ); R ) . Here we have set V k = − W k ⋆ ρ − W k ⋆ η , for k ∈ { , } , and σ = ρ + η , as above. Notice that the existence of weak solutions to system (1) will follow directly from the convergenceof the numerical solution. Indeed, our analysis relies on a compactness argument which does notsuppose a priori existence of solution to system (1).To this end we first define the following space discretisation of the domain.
Definition 2.2 (Space discretisation) . To discretise space, we introduce the mesh T := [ i ∈ I C i , where the control volumes are given by C i = [ x i − / , x i +1 / ) for all i ∈ I := { , . . . , N } . We assumethat the measure of the control volumes are given by | C i | = ∆ x i = x i +1 / − x i − / > , for all i ∈ I .Note that x / = a , and x N +1 / = b . ONVERGENCE OF A FINITE VOLUME SCHEME 5 x i − / x i − / x i +1 / x i +3 / x i − x i x i +1 C i C i − C i +1 x / = a x N +1 / = b Figure 1.
Space discretisation according to Definition 2.2.We also define x i = ( x i +1 / + x i − / ) / C i and set ∆ x i +1 / = x i +1 − x i for i = 1 , . . . , N −
1. We assume that the mesh is regular in the sense that there exists ξ ∈ (0 ,
1) suchthat for h := max ≤ i ≤ N { ∆ x i } (4) ξ h ≤ ∆ x i ≤ h, and, as a consequence, ξ h ≤ ∆ x i +1 / ≤ h , as well.On this mesh we shall now define the semi-discrete finite volume approximation of system (1). Thediscretised initial data are given by the cell averages of the continuous initial data, i.e. ρ i := 1∆ x i Z C i ρ ( x ) d x, and η i := 1∆ x i Z C i η ( x ) d x, (5)for all i ∈ I . Next, we introduce the discrete versions of the cross-diffusion and the interactionterms. We set ( V ) i := − N X j =1 ∆ x j (cid:0) W i − j ρ j + W i − j η j (cid:1) , ( V ) i := − N X j =1 ∆ x j (cid:0) W i − j η j + W i − j ρ j (cid:1) , (6)where W i − jkl = 1∆ x j Z C j W kl ( | x i − s | )d s, (7)for k, l = 1 ,
2, and U i := − (cid:0) ρ i + η i ) , for the cross-diffusion term, respectively. Then the scheme reads d ρ i d t ( t ) = − F i +1 / − F i − / ∆ x i , d η i d t ( t ) = − G i +1 / − G i − / ∆ x i , (8a) JOS´E A. CARRILLO, FRANCIS FILBET, AND MARKUS SCHMIDTCHEN for i ∈ I . Here the numerical fluxes are given by F i +1 / = h ν (d U ) + i +1 / + (d V ) + i +1 / i ρ i + h ν (d U ) − i +1 / + (d V ) − i +1 / i ρ i +1 − ǫ ρ i +1 − ρ i ∆ x i +1 / , G i +1 / = h ν (d U ) + i +1 / + (d V ) + i +1 / i η i + h ν (d U ) − i +1 / + (d V ) − i +1 / i η i +1 − ǫ η i +1 − η i ∆ x i +1 / , (8b)for i = 1 , . . . , N −
1, with the numerical no-flux boundary condition F / = F N +1 / = 0 , and G / = G N +1 / = 0 , (8c)where we introduced the discrete gradient d u i +1 / asd u i +1 / := u i +1 − u i ∆ x i +1 / . As usual, we use ( z ) ± to denote the positive (resp. negative) part of z , i.e. ( z ) + := max( z, , and ( z ) − := min( z, . At this stage, the numerical flux (8b) may look strange since • the cross-diffusion term is approximated as a convective term using that ∂∂x (cid:18) ρ ∂∂x ( ρ + η ) (cid:19) = ∂∂x (cid:18) ρ ∂σ∂x (cid:19) with σ = ρ + η and ∂σ∂x is considered as a velocity field. This treatment has already beenused in [9] and allows to preserve the positivity of both discrete densities ( ρ, η ) (see Lemma3.1), which is crucial for the convergence analysis. • In this new formulation, the velocity field is split in two parts both treated by an upwindscheme. One part comes from the cross-diffusion part, and the second one comes from thenon-local interaction fields. This splitting is crucial to recovering a consistent dissipativeterm for the discrete energy estimate corresponding to Eq. (2).
Definition 2.3 (Piecewise constant approximation) . For a given mesh T h we define the approximatesolution to system (1) by ρ h ( t, x ) := ρ i ( t ) , and η h ( t, x ) := η i ( t ) , for all ( t, x ) ∈ [0 , T ] × C i , with i = 1 , . . . , N . Moreover, we define the following approximations ofthe gradients d ρ h ( t, x ) = ρ i +1 − ρ i ∆ x i +1 / , and d η h ( t, x ) = η i +1 − η i ∆ x i +1 / for ( t, x ) ∈ [0 , T ) × [ x i , x i +1 ) , for i = 1 , . . . , N − . Furthermore, in order to define d ρ h and d η h onthe whole interval ( a, b ) we set them to zero on ( a, x ) and ( x N , b ). ONVERGENCE OF A FINITE VOLUME SCHEME 7
Notice that the discrete gradients (d ρ h , d η h ) are piecewise constant just like ( ρ h , η h ) however noton the same partition of the interval ( a, b ). We have set out all definitions necessary to formulatethe convergence of the numerical scheme (8). Theorem 2.4 (Convergence to a weak solution.) . Let ρ , η ∈ L ( a, b ) ∩ L ∞ + ( a, b ) be some initialdata and Q T := (0 , T ) × ( a, b ) . Then, ( i ) there exists a nonnegative approximate solution ( ρ h , η h ) in the sense of Definition 2.3; ( ii ) up to a subsequence, this approximate solution converges strongly in L ( Q T ) to ( ρ, η ) ∈ L ( Q T ) , where ( ρ, η ) is a weak solution as in Definition 2.1. Furthermore we have ρ , η ∈ L (0 , T ; H ( a, b )) ; ( iii ) as a consequence system (1) has a weak solution. A priori estimates
This section is dedicated to deriving a priori estimates for our system. In order to do so we requirethe positivity of approximate solutions and their conservation of mass, respectively. The followinglemma guarantees these properties.
Lemma 3.1 (Existence of nonnegative solutions and conservation of mass) . Assume that the ini-tial data ( ρ , η ) in non-negative. Then there exists a unique nonnegative approximate solution ( ρ h , η h ) h> to the scheme (8a)-(8b). Furthermore, the finite volume scheme conserves the initialmass of both densities.Proof. On the one hand we notice that the right-hand side of (8a)-(8b) is locally Lipschitz withrespect to ( ρ i , η i ) ≤ i ≤ N . Hence, we may apply the Cauchy-Lipschitz theorem to obtain a uniquecontinuously differentiable local-in-time solution.On the other hand to prove that this solution is global in time, we show the nonnegativity of thesolution together with the conservation of mass and argue by contradiction.On a given mesh, let some initial data, ρ i (0) , η i (0) ≥
0, be given for i = 1 , . . . N . We rewrite thescheme in the following way.d ρ i d t ( t ) = − F i +1 / − F i − / ∆ x i = 1∆ x i (cid:0) A i ρ i + B i ρ i +1 + C i ρ i − (cid:1) , (9)where A i = ν (d U ) − i − / + (d V ) − i − / − ν (d U ) + i +1 / − (d V ) + i +1 / − ǫ (cid:18) ρ i ∆ x i +1 / + ρ i ∆ x i − / (cid:19) ,B i = − ν (d U ) − i +1 / − (d V ) − i +1 / + ǫ ρ i +1 x i +1 / ,C i = ν (d U ) + i − / + (d V ) + i − / + ǫ ρ i − x i − / . Then let t ⋆ ≥ i.e. t ⋆ = sup { t ≥ | , ρ i ( s ) ≥ , for all s ∈ [0 , t ] , and i = 1 , . . . , N } . If t ⋆ < ∞ , then there exists a nonincreasing sequence ( t k ) k ∈ N such that t k > t ⋆ , t k → t ⋆ as k → ∞ and there exists i k ∈ { , . . . , N } verifying ρ i k ( t k ) < , ∀ k ∈ N . JOS´E A. CARRILLO, FRANCIS FILBET, AND MARKUS SCHMIDTCHEN
Since the index i k takes a finite number of integer values, we can extract a nonincreasing subsequenceof ( t k ) k ∈ N still labeled in the same manner such that there exists an index j ∈ { , . . . , N } and ρ j ( t k ) < , ∀ k ∈ N , where t k → t ⋆ , as k goes to infinity.Also note by continuity of ( ρ i ) ≤ i ≤ N , we have that ρ i ( t ⋆ ) ≥ i ∈ { , . . . , N } .By the above computation, Eq. (9), we see that, if ρ j +1 ( t ⋆ ) > ρ j − ( t ⋆ ) >
0, theneither B j ( t ⋆ ) > C j ( t ⋆ ) > ρ j d t ( t ⋆ ) = 1∆ x j ( A j ρ j ( t ⋆ ) + B j ρ j +1 ( t ⋆ ) + C j ρ j − ( t ⋆ ))= 1∆ x j ( B j ρ j +1 ( t ⋆ ) + C j ρ j − ( t ⋆ )) > , hence there exists τ > t ∈ [ t ⋆ , t ⋆ + τ ), we have ρ j ( t ) > ρ j ( t ⋆ ) = 0, whichcannot occur since ρ j is continuous and for t k > t ⋆ , ρ j ( t k ) < k ∈ N with t k → t ⋆ when k goes to infinity.If ρ j − ( t ⋆ ) = ρ j ( t ⋆ ) = ρ j +1 ( t ⋆ ) = 0 then by uniqueness of the solution, we have that ρ j ≡ t ≥ t ⋆ , which contradicts again that ρ j ( t k ) < k ∈ N large enough.Finally we get the conservation of mass,dd t Z ba ρ h ( t, x )d x = N X i =1 ∆ x i dd t ρ i = N X i =1 ∆ x i F i +1 / − F i − / ∆ x i = F N +1 / − F / = 0 , by the no-flux condition. Analogously, the second species remains nonnegative and its mass isconserved as well. As a consequence of the control of the L -norm of ( ρ h , η h ) we can extend thelocal solution to a global, nonnegative solution. (cid:3) Now, we are ready to study the evolution of the energy of the system on the semi-discrete level.The remaining part of this section is dedicated to proving the following lemma – an estimate similarto (2) for the semi-discrete scheme (8).
Lemma 3.2 (Energy control) . Consider a solution of the semi-discrete scheme (8a) - (8b) . Thenwe have dd t N X i =1 ∆ x i [ ρ i log ρ i + η i log η i ] + N − X i =1 ∆ x i +1 / h ν | d U i +1 / | + ǫ (cid:0) | d ρ i +1 / | + | d η i +1 / | (cid:1)i ≤ C ǫ , where the constant C ǫ > is given by (10) C ǫ = ( b − a ) ǫ (cid:0) (cid:0) k W ′ k L ∞ + k W ′ k L ∞ (cid:1) m + (cid:0) k W ′ k L ∞ + k W ′ k L ∞ (cid:1) m (cid:1) . Proof.
Upon using the scheme, Eq. (8a), we getdd t N X i =1 ∆ x i ρ i log ρ i = − N X i =1 ( F i +1 / − F i − / ) log ρ i . ONVERGENCE OF A FINITE VOLUME SCHEME 9
By discrete integration by parts and the no-flux condition, Eq. (8c), we obtaindd t N X i =1 ∆ x i ρ i log ρ i = N − X i =1 ∆ x i +1 / F i +1 / dlog ρ i +1 / = ν N − X i =1 ∆ x i +1 / (cid:18) (d U ) + i +1 / ρ i + (d U ) − i +1 / ρ i +1 (cid:19) dlog ρ i+1 / + N − X i =1 ∆ x i +1 / (cid:18) (d V ) + i +1 / ρ i + (d V ) − i +1 / ρ i +1 (cid:19) dlog ρ i+1 / , − ǫ N − X i =1 (cid:0) ρ i +1 − ρ i (cid:1) dlog ρ i+1 / , where, in the last equality, we substituted the definition of the numerical flux, Eq. (8b). Let usdefine ˜ ρ i +1 / := ρ i +1 − ρ i log ρ i +1 − log ρ i , if ρ i = ρ i +1 ,ρ i + ρ i +1 , else , (11)for i ∈ { , . . . , N − } , and note that then ˜ ρ i +1 / ∈ [ ρ i , ρ i +1 ] by concavity of the log. Reorderingthe terms, we obtaindd t N X i =1 ∆ x i ρ i log ρ i − N − X i =1 ∆ x i +1 / h ν d U i +1 / ˜ ρ i +1 / − ǫ ρ i +1 / i dlog ρ i +1 / = ν N − X i =1 ∆ x i +1 / (cid:16) (d U ) + i +1 / ( ρ i − ˜ ρ i +1 / ) + (d U ) − i +1 / ( ρ i +1 − ˜ ρ i +1 / ) (cid:17) dlog ρ i +1 / + N − X i =1 ∆ x i +1 / (cid:16) (d V ) + i +1 / ( ρ i − ˜ ρ i +1 / ) + (d V ) − i +1 / ( ρ i +1 − ˜ ρ i +1 / ) (cid:17) dlog ρ i +1 / + N − X i =1 ∆ x i +1 / ˜ ρ i +1 / d V ,i +1 / dlog ρ i +1 / . (12)Thus, using ˜ ρ i +1 / ∈ [ ρ i , ρ i +1 ] and the monotonicity of log, we note that ( ρ i − ˜ ρ i +1 / )dlog ρ i +1 / (cid:0) ν (d U ) + i +1 / + (d V ) + i +1 / (cid:1) ≤ , ( ρ i +1 − ˜ ρ i +1 / )dlog ρ i +1 / (cid:0) ν (d U ) − i +1 / + (d V ) − i +1 / (cid:1) ≤ . (13)This is easy to see, for if ρ i = ρ i +1 we observe dlog ρ i +1 / = 0 and Eqs. (13) hold with equality. Inthe case of ρ i < ρ i +1 we observe( ρ i − ˜ ρ i +1 / ) | {z } ≤ dlog ρ i +1 / | {z } ≥ (cid:0) ν (d U ) + i +1 / + (d V ) + i +1 / (cid:1)| {z } ≥ ≤ , while, for ρ i > ρ i +1 there also holds( ρ i − ˜ ρ i +1 / ) | {z } ≥ dlog ρ i +1 / | {z } ≤ (cid:0) ν (d U ) + i +1 / + (d V ) + i +1 / (cid:1)| {z } ≥ ≤ , whence we infer the inequality. The same argument can be applied in order to obtain the secondline of Eq. (13). Thus we may infer from Eq. (12) thatdd t N X i =1 ∆ x i ρ i log ρ i − N − X i =1 ∆ x i +1 / (cid:16) ν ˜ ρ i +1 / d U i +1 / − ǫ ρ i +1 / (cid:17) dlog ρ i +1 / ≤ N − X i =1 ∆ x i +1 / ˜ ρ i +1 / dlog ρ i +1 / d V ,i +1 / . Note that the definition of ˜ ρ i +1 / in Eq. (11), is consistent with the case ρ i = ρ i +1 and there holds˜ ρ i +1 / dlog ρ i +1 / = d ρ i +1 / , whence we getdd t N X i =1 ∆ x i ρ i log ρ i − ν N − X i =1 ∆ x i +1 / d ρ i +1 / d U i +1 / + ǫ N − X i =1 ∆ x i +1 / d ρ i +1 / dlog ρ i +1 / ≤ N − X i =1 ∆ x i +1 / d ρ i +1 / d V ,i +1 / . Furthermore, we notice that12 d ρ i +1 / dlog ρ i +1 / = ρ i +1 + ρ i ρ i +1 / | d ρ i +1 / | ≥ | d ρ i +1 / | , where we employed Eq. (11). Hence we havedd t N X i =1 ∆ x i ρ i log ρ i − ν N − X i =1 ∆ x i +1 / d ρ i +1 / d U i +1 / + ǫ N − X i =1 ∆ x i +1 / | d ρ i +1 / | ≤ N − X i =1 ∆ x i +1 / d ρ i +1 / d V ,i +1 / . (14)A similar computation can be applied to the second species, which yieldsdd t N X i =1 ∆ x i η i log η i − ν N − X i =1 ∆ x i +1 / d η i +1 / d U i +1 / + ǫ N − X i =1 ∆ x i +1 / | d η i +1 / | ≤ N − X i =1 ∆ x i +1 / d η i +1 / d V ,i +1 / . (15) ONVERGENCE OF A FINITE VOLUME SCHEME 11
Upon adding up equations (14) and (15), we obtaindd t N X i =1 ∆ x i [ ρ i log ρ i + η i log η i ] − ν N − X i =1 ∆ x i +1 / d U i +1 / d ( ρ + η ) i +1 / + ǫ N − X i =1 ∆ x i +1 / (cid:18) | d ρ i +1 / | + | d η i +1 / | (cid:19) ≤ R h , (16)where R h is given by R h = N − X i =1 ∆ x i +1 / (cid:2) d V ,i +1 / d ρ i +1 / + d V ,i +1 / d η i +1 / (cid:3) . Finally we notice that(17) R h ≤ N − X i =1 ∆ x i +1 / " | d V ,i +1 / | α + α | d ρ i +1 / | + | d V ,i +1 / | α + α | d η i +1 / | , for any α >
0, by Young’s inequality. Observing, that for k = 1 , | d V k,i +1 / | ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N − X j =1 ρ j Z C j W k ( x i +1 − y ) − W k ( x i − y )∆ x i +1 / d y (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N − X j =1 η j Z C j W k ( x i +1 − y ) − W k ( x i − y )∆ x i +1 / d y (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ k W ′ k k L ∞ m + k W ′ k k L ∞ m and by conservation of positivity and mass, it gives for k = 1 , α N − X i =1 ∆ x i +1 / | d V k,i +1 / | ≤ ( b − a )2 α (cid:0) m k W ′ k k L ∞ + m k W ′ k k L ∞ (cid:1) . Thus Eq. (17) becomes R h ≤ α N − X i =1 ∆ x i +1 / (cid:0) | d ρ i +1 / | + | d η i +1 / | (cid:1) + C α , where C α is given in Eq. (10). Finally, substituting the latter estimate into Eq. (16), we obtaindd t N X i =1 ∆ x i [ ρ i log ρ i + η i log η i ] + ν N − X i =1 ∆ x i +1 / | d U i +1 / | + ǫ − α N − X i =1 ∆ x i +1 / (cid:0) | d ρ i +1 / | + | d η i +1 / | (cid:1) ≤ C α , for any solution ( ρ i ) i ∈ I , ( η i ) i ∈ I of the semi-discrete scheme (8). Hence choosing α = ǫ/ (cid:3) Corollary 3.3 (A priori bounds) . Let ( ρ i ) i ∈ I , ( η i ) i ∈ I be solutions of the semi-discrete scheme (8) .Then there exists a constant C > such that Z T N − X i =1 ∆ x i +1 / (cid:0) | d ρ i +1 / | + | d η i +1 / | + | d U i +1 / | (cid:1) d t ≤ C. Proof.
Using the fact that x log x ≥ − log( e ) /e , i.e. x log x is bounded from below, yields X i ∈ I ∆ x i [ ρ i log ρ i + η i log η i ]( t ) ≥ − e ) e ( b − a ) =: − C . Hence the discrete version of the classical entropy functional is bounded from below. Therefore, weintegrate the inequality of Lemma 3.2 in time and get Z T N − X i =1 ∆ x i +1 / h ν | d U i +1 / ( t ) | + ǫ (cid:0) | d ρ i +1 / ( t ) | + | d η i +1 / ( t ) | (cid:1)i d t ≤ C ǫ T + C + Z ba ρ | log ρ | + η | log η | d x, which proves the statement. (cid:3) Thanks to a classical discrete Poincar´e inequality [28, Lemma 3.7] and [8], we get uniform L -estimates on the discrete approximation ( ρ h , η h ) h> . Lemma 3.4.
Let ( ρ i ) i ∈ I , ( η i ) i ∈ I be the numerical solutions obtained from scheme (8) . Then thereholds k ρ h k L ( Q T ) + k η h k L ( Q T ) ≤ C, for some constant C > independent of h > . Proof of Theorem 2.4
This section is dedicated to proving compactness of both species, the fluxes, and the regularisingporous-medium type diffusion. Upon establishing the compactness result we identify the limits asweak solutions in the sense of Definition 2.1.First by application of Lemma 3.1, we get existence and uniqueness of a nonnegative approximatesolution ( ρ h , η h ) to (8a)-(8b). Hence the first item of Theorem 2.4 is proven. Now let us investigatethe asymptotic h → Strong compactness of approximate solutions.
We shall now make use of the aboveestimates in order to obtain strong compactness of both species, ( ρ h , η h ) in L ( Q T ). Lemma 4.1 (Strong compactness in L ( Q T )) . Let ( ρ h , η h ) h> be the approximation to system (1) obtained by the semi-discrete scheme (8) . Then there exist functions ρ, η ∈ L ( Q T ) such that ρ h → ρ, and η h → η, strongly in L ( Q T ) , up to a subsequence. ONVERGENCE OF A FINITE VOLUME SCHEME 13
Proof.
We invoke the compactness criterion by Aubin and Lions [23]. Accordingly, a set P ⊂ L (0 , T ; B ) is relatively compact if P is bounded in L (0 , T ; X ) and the set of derivatives { ∂ t ρ (cid:12)(cid:12) ρ ∈ P } is bounded in a third space L (0 , T ; Y ), whenever the involved Banach spaces satisfy X ֒ → ֒ → B ֒ → Y , i.e. the first embedding is compact and the second one continuous. For our purpose wechoose X := BV ( a, b ), B := L ( a, b ), and Y := H − ( a, b ). The first embedding is indeed compact, e.g. Ref. [1, Theorem 10.1.4] and the second one is continuous.In the second step we show the time derivatives are bounded in L (0 , T ; H − ( a, b )). To this end,let ϕ ∈ C ∞ c (( a, b )). Throughout, we write h· , ·i for h· , ·i H − ,H for the dual pairing. Making use ofthe scheme, there holds (cid:28) d ρ h d t , ϕ (cid:29) = N X i =1 Z C i d ρ i d t ϕ d x = − N X i =1 F i +1 / − F i − / ∆ x i Z C i ϕ d x, having used the scheme, Eq. (8a). Next we set ϕ i := 1∆ x i Z C i ϕ d x, perform a discrete integration by parts and use the no-flux boundary conditions, Eq. (8c), to obtain (cid:28) d ρ h d t , ϕ (cid:29) = N − X i =1 F i +1 / ( ϕ i +1 − ϕ i ) . Using the definition of the numerical flux, Eq. (8b), we get (cid:28) d ρ h d t , ϕ (cid:29) = N − X i =1 h(cid:16) ν (d U ) + i +1 / + (d V ) + i +1 / (cid:17) ρ i + (cid:16) ν (d U ) − i +1 / + (d V ) − i +1 / (cid:17) ρ i +1 i ( ϕ i +1 − ϕ i ) − ǫ N − X i =1 ρ i +1 − ρ i ∆ x i +1 / ( ϕ i +1 − ϕ i )Let us begin with the self-diffusion part. Using the Cauchy-Schwarz inequality, we estimate thediscrete gradient and ρ itself by Corollary 3.3 and Lemma 3.4 ǫ Z T N − X i =1 ρ i +1 − ρ i ∆ x i +1 / ( ϕ i +1 − ϕ i ) d t ≤ ǫ (cid:13)(cid:13)(cid:13)(cid:13) ∂ϕ∂x (cid:13)(cid:13)(cid:13)(cid:13) L ∞ Z T N − X i =1 ∆ x i +1 / | d ρ i +1 / | ( ρ i +1 + ρ i ) d t ≤ ǫ (cid:13)(cid:13)(cid:13)(cid:13) ∂ϕ∂x (cid:13)(cid:13)(cid:13)(cid:13) L ∞ Z T N − X i =1 ∆ x i +1 / | d ρ i +1 / | d t ! / Z T N X i =1 ξ − ∆ x i | ρ i | d t ! / ≤ C k ϕ k H ( a,b ) , (18)where we used that ϕ ′ ∈ H ⊂ L ∞ . Next, we address the cross-diffusion and non-local interactions terms using the same argument. Forinstance for the cross-diffusive part, we have ν Z T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N − X i =1 h (d U ) + i +1 / ρ i + (d U ) − i +1 / ρ i +1 i(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) | ϕ i +1 − ϕ i | d t ≤ ν √ ξ (cid:13)(cid:13)(cid:13)(cid:13) ∂ϕ∂x (cid:13)(cid:13)(cid:13)(cid:13) L ∞ Z T N − X i =1 ∆ x i +1 / | d U i +1 / | d t ! / Z T N X i =1 ∆ x i ρ i d t ! / ≤ C k ϕ k H ( a,b ) , where we used Corollary 3.3 and Lemma 3.4 again. The non-local interaction term is estimated inthe same way, thus there holds Z T (cid:12)(cid:12)(cid:12)(cid:12)(cid:28) d ρ h d t , ϕ (cid:29)(cid:12)(cid:12)(cid:12)(cid:12) d t ≤ C k ϕ k H ( a,b ) . By density of C ∞ c (( a, b )) in H ( a, b ) we may infer the boundedness of ( d ρ h d t ) h> in L (0 , T ; H − ( a, b )),which concludes the proof. (cid:3) From the latter result we can prove the convergence of the discrete advection field d V ,h and d V ,h defined as in Definition 2.3. Lemma 4.2.
For any ≤ p ≤ ∞ and k ∈ { , } , the piecewise constant approximation d V k,h converges strongly in L (0 , T ; L ( a, b )) to − ( W ′ k ⋆ ρ + W ′ k ⋆ η ) , where ( ρ, η ) corresponds to thelimit obtained in Lemma 4.1.Proof. Let k ∈ { , } . For each i = 0 , . . . , N −
1, and x ∈ [ x i , x i +1 ) we haved V k, h ( x ) = (d V k ) i +1 / = − N X j =1 Z C j W k ( x i +1 − y ) − W k ( x i − y )∆ x i +1 / ρ j d y − N X j =1 Z C j W k ( x i +1 − y ) − W k ( x i − y )∆ x i +1 / η j d y. We define V ′ k,h and V ′ k as V ′ k,h ( x ) := − W ′ k ⋆ ρ h − W ′ k ⋆ η h ,V ′ k ( x ) := − W ′ k ⋆ ρ − W ′ k ⋆ η. On the one hand from the strong convergence of ( ρ h , η h ) to ( ρ, η ) in L (0 , T ; L ( a, b )) and theconvolution product’s properties, we obtain(19) k V ′ k,h − V ′ k k L (0 ,T ; L ( a,b )) → , when h → . ONVERGENCE OF A FINITE VOLUME SCHEME 15
On the other hand, we have for any x ∈ [ x i , x i +1 ) | d V k, h ( x ) − V ′ k, h ( x ) | ≤ N X j =1 Z C j (cid:12)(cid:12)(cid:12)(cid:12) W k ( x i +1 − y ) − W k ( x i − y )∆ x i +1 / − W ′ k ( x − y ) (cid:12)(cid:12)(cid:12)(cid:12) ρ j d y, + N X j =1 Z C j (cid:12)(cid:12)(cid:12)(cid:12) W k ( x i +1 − y ) − W k ( x i − y )∆ x i +1 / − W ′ k ( x − y ) (cid:12)(cid:12)(cid:12)(cid:12) η j d y, ≤ (cid:0) k W ′′ k k L ∞ m + k W ′′ k k L ∞ m (cid:1) h, hence there exists a constant C > | d V k, h ( x ) − V ′ k, h ( x ) | ≤ C p h . Integrating over x ∈ [ x i , x i +1 ) and summing over i ∈ { , . . . , N − } , we get that(20) k d V k,h − V ′ k k L (0 ,T ; L ( a,b )) → , when h → . Notice that ( x , x N ) ⊂ ( a, b ) where x → a and x N → b as h →
0. From Eqs. (19) and (20) we getthat k d V k,h − V ′ k, k L (0 ,T ; L ( a,b )) goes to zero as h tends to zero. (cid:3) Weak compactness for the discrete gradients.
In the previous section we have estab-lished the strong L -convergence of both species, ( ρ h ) h> and ( η h ) h> . However, in order to be ableto pass to the limit in the cross-diffusion term ρ h (d ρ h + d η h ) we need to establish weak convergencein the discrete gradients in L . This is done in the following proposition. Proposition 1 (Weak convergence of the derivatives) . The discrete spatial derivatives, defined inDefinition 2.3, satisfy d β h converges weakly to ∂β∂x in L ( Q T ) and β ∈ L (0 , T ; H ( a, b )) , where β ∈ { ρ, η, U } Proof.
Take β ∈ { ρ, η, U } , hence from Lemma 4.1, we know that β h → β strongly in L ( Q T ).Furthermore, from Corollary 3.3 we also deduce that d β h weakly converges to some function r ∈ L ( Q T ).Let us show that β ∈ L (0 , T, H ( a, b )) and r = ∂β∂x . First, we have for any t ∈ [0 , T ] and any ϕ ∈ C ∞ c ((0 , T ) × ( a, b )), Z Q T β h ( t ) ∂ϕ∂x d x = Z T N X i =1 β i ( t ) (cid:2) ϕ ( t, x i +1 / ) − ϕ ( t, x i − / ) (cid:3) d t = − Z T N − X i =1 ∆ x i +1 / d β i +1 / ( t ) ϕ ( t, x i +1 / ) d t, having used discrete integration by parts and the fact that ϕ is compactly supported, i.e. ϕ ( t, x N +1 / ) = ϕ ( t, x / ) = 0. Then, by Definition 2.3 on the discrete gradient, we may consider (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)Z T N − X i =1 Z x i +1 x i d β i +1 / ϕ ( t, x )d x d t + Z T Z ba β h ∂ϕ∂x d x d t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ Z T N − X i =1 Z x i +1 x i (cid:12)(cid:12) d β i +1 / (cid:12)(cid:12) (cid:12)(cid:12) ϕ ( t, x ) − ϕ ( t, x i +1 / ) (cid:12)(cid:12) d x d t ≤ (cid:13)(cid:13)(cid:13)(cid:13) ∂ϕ∂x (cid:13)(cid:13)(cid:13)(cid:13) ∞ Z T N − X i =1 ∆ x i +1 / (cid:12)(cid:12) d β i +1 / (cid:12)(cid:12) d t ! / Z T N − X i =1 ∆ x i +1 / d t ! / ≤ (cid:13)(cid:13)(cid:13)(cid:13) ∂ϕ∂x (cid:13)(cid:13)(cid:13)(cid:13) ∞ C / T / √ b − a h, This yields the statement, when h →
0, for we have Z T N − X i =1 Z x i +1 x i d β i +1 / ϕ ( t, x ) d x d t + Z T Z ba β h ∂ϕ∂x d x d t → , (21)which proves that d β h converges weakly to ∂β∂x , as h → β ∈ L (0 , T ; H ( a, b )). (cid:3) Passing to the limit.
We have now garnered all information necessary to prove Theorem2.4. For brevity we shall only show the convergence result for ρ , as it follows for η similarly, usingthe same arguments. Let ϕ ∈ C ∞ c ([0 , T ) × ( a, b )) be a test function. We introduce the followingnotations: E h := Z T Z ba ρ h ∂ϕ∂t d x d t + Z ba ρ h (0) ϕ (0) d x, A h := Z T Z ba d V ,h ρ h ∂ϕ∂x d x d t, C h := ν Z T Z ba d U h ρ h ∂ϕ∂x d x d t, D h := ǫ Z T Z ba ρ h ∂ ϕ∂x d x d t. and ε ( h ) := E h + A h + C h + D h . On the other hand, we set ϕ i ( t ) = 1∆ x i Z C i ϕ ( t, x ) d x, and multiply the scheme, Eq. (8a), by the test function and integrate in time and space to get(22) E h + A ,h + C ,h + D ,h = 0 , ONVERGENCE OF A FINITE VOLUME SCHEME 17 where A ,h := N − X i =1 Z T ∆ x i +1 / h (d V ) + i +1 / ρ i + (d V ) − i +1 / ρ i +1 i d ϕ i +1 / ( t ) d t, C ,h := ν N − X i =1 Z T ∆ x i +1 / h (d U ) + i +1 / ρ i + (d U ) − i +1 / ρ i +1 i d ϕ i +1 / ( t ) d t, D ,h := − ǫ N − X i =1 Z T (cid:2) ρ i +1 − ρ i (cid:3) d ϕ i +1 / ( t ) d t. When h tends to zero and from the strong convergence of ( ρ h , η h ) h> to ( ρ, η ) in L ( Q T ), thestrong convergence of (d V k,h ) h> to V ′ k in L ( Q T ) and the weak convergence of the discrete gradient(d U h ) h> to − ∂σ∂x in L ( Q T ), it is easy to see that ε ( h ) → Z T Z ba (cid:26) ρ (cid:20) ∂ϕ∂t + (cid:18) ∂V ∂x − ν ∂∂x ( ρ + η ) (cid:19) ∂ϕ∂x (cid:21) + ǫ ρ ∂ ϕ∂x (cid:27) d x d t + Z ba ρ (0) ϕ (0) d x, when h →
0. Therefore it suffices to prove that ε ( h ) →
0, as h goes to zero, which will be achievedby proving that A h − A ,h , C h − C ,h and D h − D ,h vanish in the limit h → The self-diffusion part D h − D ,h . On the one hand, after a simple integration we get D h = ǫ N X i =1 Z T ρ i ( t ) (cid:20) ∂ϕ∂x ( t, x i +1 / ) − ∂ϕ∂x ( t, x i − / ) (cid:21) d t = − ǫ N − X i =1 Z T (cid:2) ρ i +1 ( t ) − ρ i ( t ) (cid:3) ∂ϕ∂x ( t, x i +1 / ) d t. Hence, we have D h − D ,h = − ǫ N − X i =1 Z T (cid:2) ρ i +1 ( t ) − ρ i ( t ) (cid:3) (cid:20) ∂ϕ∂x ( t, x i +1 / ) − d ϕ i +1 / ( t ) (cid:21) d t and observing that (cid:12)(cid:12)(cid:12)(cid:12) ∂ϕ∂x ( t, x i +1 / ) − d ϕ i +1 / ( t ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:13)(cid:13)(cid:13)(cid:13) ∂ ϕ∂x (cid:13)(cid:13)(cid:13)(cid:13) L ∞ h, we obtain, in conjunction with the Cauchy-Schwarz inequality and the a priori bounds establishedin Corollary 3.3, that and Lemma 3.4, |D h − D ,h | ≤ ǫ (cid:13)(cid:13)(cid:13)(cid:13) ∂ ϕ∂x (cid:13)(cid:13)(cid:13)(cid:13) L ∞ N − X i =1 Z T ∆ x i +1 / | d ρ i +1 / | d t ! / k ρ h k L ( Q T ) ξ / h (23) ≤ C h, in the virtue of the estimate Eq. (18).
The cross-diffusion part.
Let us now treat the cross-diffusion part. This term is more complicatedsince it involves the piecewise constant functions ρ h and d U h , which are not defined on the samemesh. Thus, on the one hand we reformulate the discrete cross-diffusion term C ,h as C ,h = C ,h + C ,h with C ,h = ν N − X i =1 Z T ∆ x i +1 / (d U ) − i +1 / [ ρ i +1 − ρ i ] d ϕ i +1 / ( t ) d t and C ,h = ν N − X i =1 Z T ∆ x i +1 / ρ i d U i +1 / d ϕ i +1 / ( t ) d t, where a direct computation and the application of Corollary 3.3 and Lemma 3.4 yield |C ,h | ≤ ν (cid:13)(cid:13)(cid:13)(cid:13) ∂ϕ∂x (cid:13)(cid:13)(cid:13)(cid:13) L ∞ k d U h k L ( Q T ) k d ρ h k L ( Q T ) h (24) ≤ C h.
On the other hand, the term C h can be rewritten as C h = ν Z T N − X i =1 d U i +1 / ( t ) Z x i +1 x i ρ h ∂ϕ∂x d x d t. Since Z x i +1 x i ρ h ∂ϕ∂x d x = ρ i (cid:2) ϕ ( t, x i +1 / ) − ϕ ( t, x i ) (cid:3) + ρ i +1 (cid:2) ϕ ( t, x i +1 ) − ϕ ( t, x i +1 / ) (cid:3) , = [ ρ i − ρ i +1 ] (cid:2) ϕ ( t, x i +1 / ) − ϕ ( t, x i ) (cid:3) + ρ i +1 [ ϕ ( t, x i +1 ) − ϕ ( t, x i )] , the term C h can be decomposed as C h = C ,h + C ,h with C ,h = − ν Z T N − X i =1 d U i +1 / [ ρ i +1 − ρ i ] (cid:2) ϕ ( t, x i +1 / ) − ϕ ( t, x i ) (cid:3) d t and C ,h = ν Z T N − X i =1 d U i +1 / ρ i [ ϕ ( t, x i +1 ) − ϕ ( t, x i )] d t. Similarly to (24), the first term C ,h can be estimated as(25) |C ,h | ≤ C h, whereas the second term C ,h is compared to C ,h |C ,h − C ,h | ≤ ν Z T N − X i =1 ∆ x i +1 / | d U i +1 / | ρ i (cid:12)(cid:12)(cid:12)(cid:12) ϕ ( t, x i +1 ) − ϕ ( t, x i )∆ x i +1 / − dϕ i +1 / ( t ) (cid:12)(cid:12)(cid:12)(cid:12) d t. Using a second order Taylor expansion of ϕ at x i and x i +1 , it yields that (cid:12)(cid:12)(cid:12)(cid:12) ϕ ( t, x i +1 ) − ϕ ( t, x i )∆ x i +1 / − dϕ i +1 / ( t ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ C h,
ONVERGENCE OF A FINITE VOLUME SCHEME 19 hence we get from Corollary 3.3 and Lemma 3.4 that(26) |C ,h − C ,h | ≤ C h. from Corollary 3.3 and Lemma 3.4. Gathering Eqs. (24), (25), and (26), we finally obtain that(27) |C h − C ,h | ≤ C h.
The advective part.
The evaluation of A h − A ,h is along the same lines of the cross-diffusion terms C h − C ,h since the latter is treated as an advective term. Hence, thanks to Lemma 4.2, we get that(28) |A h − A ,h | ≤ C h.
Finally by definition of ε ( h ) and using Eq. (22) together with Eqs. (23), (27), and (28), we obtain | ε ( h ) | = | − ( A ,h + C ,h + D ,h ) + A h + C h + D h |≤ |A h − A ,h | + |C h − C ,h | + |D h − D ,h |≤ C h, that is, ε ( h ) →
0, when h →
0, which proves that ( ρ, η ) is a weak solution to Eq. (1). This provesthe second item of Theorem 2.4.Finally the last item concerning the existence of solutions to (1) is a direct consequence of theconvergence. 5.
A Fully Discrete Implicit Scheme
In this section we shall comment on a discrete-in-time version of the semi-discrete scheme (8). Tothis end we replace the time derivative in Eq. (8a) by simple forward differences and obtain thefollowing implicit and fully-discrete scheme ρ n +1 i − ρ ni ∆ t = − F n +1 i +1 / − F n +1 i − / ∆ x i ,η n +1 i − η ni ∆ t = − G n +1 i +1 / − G n +1 i − / ∆ x i , (29)where ∆ t >
0. System (29) gives rise to two approximating sequences ( ρ ni ) ≤ i ≤ N and ( η ni ) ≤ i ≤ N ,for 0 ≤ n ≤ M where M := ⌈ T / ∆ t ⌉ and the discrete time instances t n := n ∆ t ; cf. Theorem 5.1.Here the numerical fluxes are given by F n +1 i +1 / = h ν (d U n +1 ) + i +1 / + (d V n ) + i +1 / i ρ n +1 i + h ν (d U n +1 ) − i +1 / + (d V n ) − i +1 / i ρ n +1 i +1 − ǫ ρ n +1 i +1 ) − ( ρ n +1 i ) ∆ x i +1 / , G n +1 i +1 / = h ν (d U n +1 ) + i +1 / + (d V n ) + i +1 / i η n +1 i + h ν (d U n +1 ) − i +1 / + (d V n ) − i +1 / i η n +1 i +1 − ǫ η n +1 i +1 ) − ( η n +1 i ) ∆ x i +1 / , (29a) for i = 1 , . . . , N −
1, with the numerical no-flux boundary condition F n +11 / = F n +1 N +1 / = 0 , and G n +11 / = G n +1 N +1 / = 0 , (29b)for n = 0 , . . . , M . Recall that(d U n +1 ) ± i +1 / = ( ρ n +1 i +1 + η n +1 i +1 ) − ( ρ n +1 i + η n +1 i )∆ x i +1 / ! ± , and (d V nk ) ± i +1 / = (cid:18) − N X j =1 Z C j W k ( x i +1 − y ) − W k ( x i − y )∆ x i +1 / ρ nj d y − N X j =1 Z C j W k ( x i +1 − y ) − W k ( x i − y )∆ x i +1 / η nj d y (cid:19) ± , for k = 1 , ρ h ( t, x ) := ρ ni , and η h ( t, x ) := η ni , for all ( t, x ) ∈ [ t n , t n +1 ) × C i , with i = 1 , . . . , N , and n = 0 , . . . , M . Moreover, we define the discreteapproximation of the spatial gradients asd x ρ h ( t, x ) = ρ ni +1 − ρ ni ∆ x i +1 / , and d x η h ( t, x ) = η ni +1 − η ni ∆ x i +1 / , for ( t, x ) ∈ [ t n , t n +1 ) × [ x i , x i +1 ), for i = 1 , . . . , N − n = 0 , . . . , M . As above, we set the discretegradients to zero on ( a, x ) and ( x N , b ). Furthermore, we define the discrete time derivative asd t ρ h ( t, x ) = ρ n +1 i − ρ ni ∆ t , and d t η h ( t, x ) = η ni − η ni ∆ t , for ( t, x ) ∈ [ t n , t n +1 ) × C i , for i = 1 , . . . , N and n = 0 , . . . , M − Theorem 5.1 (Existence and uniqueness result) . Let ρ i , η i be nonnegative initial data with mass m and m , respectively, and assume the following CFL condition
16 ( m + m ) ∆ t ( ξ h ) < . (30) Then there exists a unique nonnegative solution ( ρ ni , η ni ) to scheme (29) , (29a) and (29b) .Proof. We show existence first and prove uniqueness later. Suppose we are given ( ρ ni ) ≤ i ≤ N and( η ni ) ≤ i ≤ N from some previous iteration. In order to construct the next iteration we shall employBrouwer’s fixed point theorem. It is easy to verify that the set X := ( ( ρ, η ) ∈ R N | ∀ ≤ i ≤ N : ρ i , η i ≥ , N X i =1 ∆ x i ρ i ≤ m , and N X i =1 ∆ x i η i ≤ m ) , ONVERGENCE OF A FINITE VOLUME SCHEME 21 is a convex and compact subset of R N × R N . Hence, we define the fixed point operator S : X → R N × R N by setting ( ρ ⋆ , η ⋆ ) = S ( ρ, η ) where ( ρ ⋆ , η ⋆ ) are implicitly given as ρ ⋆i = ρ ni − ∆ t ∆ x i (cid:16) F ⋆i +1 / − F ⋆i − / (cid:17) , η ⋆i = η ni − ∆ t ∆ x i (cid:16) G ⋆i +1 / − G ⋆i − / (cid:17) for any ( ρ, η ) ∈ X , where F ⋆ and G ⋆ denote the numerical fluxes F ⋆i +1 / = h ν (d U ) + i +1 / + (d V n ) + i +1 / i ρ ⋆i + h ν (d U ) − i +1 / + (d V n ) − i +1 / i ρ ⋆i +1 − ǫ ( ρ i +1 + ρ i )2 ρ ⋆i +1 − ρ ⋆i ∆ x i +1 / , G ⋆i +1 / = h ν (d U ) + i +1 / + (d V n ) + i +1 / i η ⋆i + h ν (d U ) − i +1 / + (d V n ) − i +1 / i η ⋆i +1 − ǫ ( η i +1 + η i )2 η ⋆i +1 − η ⋆i ∆ x i +1 / , (31)for i = 1 , . . . , N − U is computed from ( ρ, η ), with the numerical no-flux boundarycondition F ⋆ / = F ⋆N +1 / = 0 , and G ⋆ / = G ⋆N +1 / = 0 . (32)Notice that for any given ( ρ, η ) ∈ X , the viscosity terms in front of the discrete gradients involved inthe definition of the fluxes F ⋆i +1 / and G ⋆i +1 / are indeed nonnegative, hence the couple ( ρ ⋆ , η ⋆ ) is welldefined since it corresponds to the unique solution of a classical fully implicit scheme in time withan upwind discretisation for the convective terms and a centred approximation for diffusive terms[29]. Moreover, since ρ n and η n are nonnegative and using the monotonicity of the numerical fluxwith respect to ( ρ ⋆ , η ⋆ ), we prove that both densities ρ ⋆ and η ⋆ are also nonnegative. Furthermore,using the nonnegativity and the no-flux conditions, we get k ρ ⋆ k L = N X i =1 ∆ x i ρ ⋆i = m , and k η ⋆ k L = N X i =1 ∆ x i η ⋆i = m , which yields that S ( X ) ⊂ X . Finally, S is continuous as the composition of continuous func-tions. Thus, we may apply Brouwer’s fixed point theorem to infer the existence of a fixed point,( ρ n +1 , η n +1 ). It now remains to show uniqueness of the fixed point.To treat in a systematic way the boundary conditions and simplify the presentation, we defineghost values for ( ρ, η ) by setting for α ∈ { ρ, η } , and k ∈ { , } , α N +1 = α N α = α and (d V nk ) N +1 / = (d V nk ) / = 0 . Then we consider two solutions (˜ ρ, ˜ η ) and ( ρ, η ) to (29). Setting h ( x ) := x / s := ρ − ˜ ρ and r := η − ˜ η , we get after substituting the two solutions to (29), for i = 1 , . . . , N , s i = − ∆ t ∆ x i (cid:16) (d V n ) + i +1 / s i + (d V n ) − i +1 / s i +1 − (d V n ) + i − / s i − − (d V n ) − i − / s i (cid:17) + ∆ t ∆ x i (cid:18) [ h ( ρ i +1 ) − h (˜ ρ i +1 )] − [ h ( ρ i ) − h (˜ ρ i )]∆ x i +1 / − [ h ( ρ i ) − h (˜ ρ i )] − [ h ( ρ i − ) − h (˜ ρ i − )]∆ x i − / (cid:19) − ∆ t ∆ x i (cid:16) (d U ) + i +1 / s i + (d U ) − i +1 / s i +1 − (d U ) + i − / s i − − (d U ) − i − / s i (cid:17) − ∆ t ∆ x i (cid:16) [(d U ) + i +1 / − (d ˜ U ) + i +1 / ] ˜ ρ i + [(d U ) − i +1 / − (d ˜ U ) − i +1 / ] ˜ ρ i +1 (cid:17) + ∆ t ∆ x i (cid:16) [(d U ) + i − / − (d ˜ U ) + i − / ] ˜ ρ i − + [(d U ) − i − / − (d ˜ U ) − i − / ] ˜ ρ i (cid:17) and a similar relation for ( r i ) ≤ i ≤ N . Applying a Taylor expansion on h ( ρ ) = h (˜ ρ ) + h ′ (ˆ ρ ) s , with ˆ ρ a convex combination of ρ and ˜ ρ , we may write(33) (∆ x i + ∆ t A i ) s i = ∆ t ( B i − s i − + C i +1 s i +1 ) − ∆ t (cid:16) [(d U ) + i +1 / − (d ˜ U ) + i +1 / ] ˜ ρ i + [(d U ) − i +1 / − (d ˜ U ) − i +1 / ] ˜ ρ i +1 (cid:17) +∆ t (cid:16) [(d U ) + i − / − (d ˜ U ) + i − / ] ˜ ρ i − + [(d U ) − i − / − (d ˜ U ) − i − / ] ˜ ρ i (cid:17) , where A i , B i − , C i +1 are nonnegative coefficients given by A i = +(d V n ) + i +1 / − (d V n ) − i − / + (d U ) + i +1 / − (d U ) − i − / + h ′ (ˆ ρ i )∆ x i +1 / + h ′ (ˆ ρ i )∆ x i − / ,B i − = +(d V n ) + i − / + (d U ) + i − / + h ′ (ˆ ρ i − )∆ x i − / ,C i +1 = − (d V n ) − i +1 / − (d U ) − i +1 / + h ′ (ˆ ρ i +1 )∆ x i +1 / . Now, we multiply equation (33) by sign( s i ) and sum over i = 1 , . . . , N , hence using that x → x ± isLipschitz continuous and observing that A i = B i + C i , with A i , B i , C i ≥
0, it yields N X i =1 ∆ x i | s i | ≤ t k ˜ ρ k ∞ N X i =1 (cid:18) | s i +1 | + | r i +1 | + | s i | + | r i | ∆ x i +1 / + | s i − | + | r i − | + | s i | + | r i | ∆ x i − / (cid:19) and in a similar way, N X i =1 ∆ x i | r i | ≤ t k ˜ η k ∞ N X i =1 (cid:18) | r i +1 | + | s i +1 | + | r i | + | s i | ∆ x i +1 / + | r i − | + | s i − | + | r i | + | s i | ∆ x i − / (cid:19) . Gathering these latter inequalities and from (4), it gives that N X i =1 ∆ x i ( | s i | + | r i | ) ≤ k (˜ ρ, ˜ η ) k ∞ ∆ t ( ξ h ) N X i =1 ∆ x i ( | s i | + | r i | ) . ONVERGENCE OF A FINITE VOLUME SCHEME 23
Finally, from the nonnegativity and the preservation of mass, we have k (˜ ρ, ˜ η ) k ∞ ≤ m + m ξh , hence under the condition (30), we conclude that s = r = 0 and the uniqueness follows. (cid:3) It is worth to mention here that the condition (30) is not optimal since we only use the discrete L -estimate on ρ and η to control the discrete gradient and the L ∞ -norm.We are now in position to state for (29), (29a) and (29b) an analogous results to the semi-discretecase Theorem 5.2 (Convergence to a weak solution of the implicit Euler discretisation) . Under theassumptions of Theorem 5.1, let ρ , η ∈ L ( a, b ) ∩ L ∞ + ( a, b ) be some initial data and Q T :=(0 , T ) × ( a, b ) as above. Then, given two nonnegative sequences ( ρ ni ) ≤ i ≤ N and ( η ni ) ≤ i ≤ N satisfying (29) , (29a) and (29b) , for any n ∈ { , . . . , M } , then ( i ) up to a subsequence, the piecewise constant approximations converge strongly in L ( Q T ) to ( ρ, η ) ∈ L ( Q T ) , where ( ρ, η ) is a weak solution as in Definition 2.1. Furthermore we have ρ , η ∈ L (0 , T ; H ( a, b )) ; ( ii ) in particular, system (1) has a weak solution.Sketch of the proof of Theorem 5.2. It is easily observed that the total mass is conserved due to thediscrete no-flux boundary conditions, cf. (29b). Together with the nonnegativity we were able toprove a semi-discrete version of the energy estimate Eq. (2) which is at the heart of the convergenceresult. Similarly as above, we are able to prove a fully discrete version of the energy estimate whichthen reads M X n =0 ∆ t N X i =1 ∆ x i ρ n +1 i log ρ n +1 i − ρ ni log ρ ni ∆ t + η n +1 i log η n +1 i − η ni log η ni ∆ t + M X n =0 ∆ t N − X i =1 ∆ x i +1 / h ν | d U i +1 / | + ǫ (cid:16) | d ρ n +1 i +1 / | + | d η n +1 i +1 / | (cid:17)i ≤ C ǫ , with C ǫ > x (log x −
1) since ρ n +1 i (log ρ n +1 i − − ρ ni (log ρ ni − ≤ log( ρ n +1 i )( ρ n +1 i − ρ ni ) . Multiplying this expression by ∆ x i / ∆ t and summing over i = 1 , . . . N and n = 0 , . . . , M we obtain M X n =0 ∆ t N X i =1 ∆ x i ρ n +1 i (log ρ n +1 i − − ρ ni (log ρ ni − t ≤ M X n =0 ∆ t N X i =1 ∆ x i log( ρ n +1 i ) ρ n +1 i − ρ ni ∆ t . The right-hand side is then substituted by the scheme and simplified along the lines of the proof ofLemma 3.2. Since the computations are exactly the same we omit them here for brevity and onlynote that it is important to set˜ ρ n +1 i +1 / := ρ n +1 i +1 − ρ n +1 i log ρ n +1 i +1 − log ρ n +1 i , if ρ n +1 i = ρ n +1 i +1 ,ρ n +1 i + ρ n +1 i +1 , else , (34) to obtain the right sign in the numerical artefacts in Eq. (13), which now read ( ρ n +1 i − ˜ ρ n +1 i +1 / )dlog ρ n +1 i +1 / (cid:0) ν (d U ) + i +1 / + (d V ) + i +1 / (cid:1) ≤ , ( ρ n +1 i +1 − ˜ ρ n +1 i +1 / )dlog ρ n +1 i +1 / (cid:0) ν (d U ) − i +1 / + (d V ) − i +1 / (cid:1) ≤ . Following the lines of the proof above we obtain the fully discrete a priori bounds N X n =0 ∆ t N − X i =1 ∆ x i +1 / (cid:16) | d x ρ n +1 i +1 / | + | d x η n +1 i +1 / | + | d x U i +1 / | (cid:17) ≤ C, for some constant C >
0, where we used ‘d x ’ to denote the discrete spatial gradient as before. Again,an application of Aubin-Lions Lemma provides relative compactness in the space L ((0 , T ) × ( a, b ))of the piecewise constant interpolations ρ h , η h . As above, the discrete gradients are uniformlybounded in L ((0 , T ) × ( a, b )) and their weak convergence is a consequence of the Banach-AlaogluTheorem. Identifying the limits as a weak solution to system (1) is shown in the same way as aboveand we leave it as an exercise for the reader. (cid:3) Remark 1 (Explicit Scheme) . We do not consider an explicit scheme here since its analysis ismuch more complicated due to the lack of uniform estimates. Indeed, an explicit scheme requiresa CFL condition on the time step which appears in the stability analysis or the energy estimateprovided in Lemma 3.2. In our case, it leads to ρ n +1 log ρ n +1 − ρ n log ρ n ∆ t = ρ n +1 − ρ n ∆ t (log ρ n −
1) + ( ρ n +1 − ρ n ) (cid:12)(cid:12) ρ n +1 / (cid:12)(cid:12) ∆ t , where ρ n +1 / belongs to the interval [ ρ n , ρ n +1 ] or [ ρ n +1 , ρ n ] . The control of this reminder termwould require some lower bounds estimates on the density ρ (see for instance [30] ). Numerical examples and validation
In this section we perform some numerical simulations of system (1) using our scheme, Eqs. (8).In Section 6.1 we test our scheme by computing the error between the numerical simulation and abenchmark solution on a finer grid. Furthermore we determine the numerical convergence order. InSection 6.2 we compute the numerical stationary states of system (1) and discuss the implicationof different cross-diffusivities and the self-diffusivities, respectively.Let us note here, that in the case of no regularising porous-medium diffusion, i.e. ǫ = 0, and certainsingular potentials the stationary states of system (1) are even known explicitly [19]. This allowsus to compare the numerical solution directly to the analytical stationary state in Section 6.2.1.Throughout the remainder of this section we apply the scheme (8) to system (1) using differentself-diffusions, ǫ , and cross-diffusions, ν .6.1. Error and numerical order of convergence.
This section is dedicated to validating ourmain convergence result, Theorem 2.4. Due to the lack of explicit solutions we compute the nu-merical solution on a fine grid and consider it a benchmark solution. We then compute numericalapproximations on coarser grids and study the error in order to obtain the numerical convergenceorder.
ONVERGENCE OF A FINITE VOLUME SCHEME 25
In all our simulations we use a fourth order Runge-Kutta scheme to solve the ordinary differentialequations Eqs. (8) with initial data Eqs. (5). The discrepancy between the benchmark solutionand numerical solutions on coarser grids is measured by the error e := ∆ t M X k =1 ∆ x N X i =1 | ρ ex ( t k , x i ) − ρ ( t k , x i ) | + | η ex ( t k , x i ) − η ( t k , x i ) | !! / . Here ρ ex , η ex denote the benchmark solutions. We use this quantity to study the convergence ofour scheme as the grid size decreases.6.1.1. No non-local interactions.
Let us begin with the purely diffusive system. We consider system(1) without any interactions, i.e. W ij ≡
0, for i, j = 1 ,
2, and we choose ν = 0 . ǫ = 0 . x = 2 − on the time interval [0 ,
10] with ∆ t = 0 .
05. Figure 2(a)shows the convergence for symmetric initial data whereas Figure 2(b) shows the convergence of thesame system in case of asymmetric initial data. In both cases we overlay a line of slope one and weconclude that the numerical convergence is of order one. −9 −8 −7 −6 −5 −4 −3 −2 Δx −6 −5 −4 −3 −2 −1 e rr o r errorlineΔofΔslopeΔ1 (a) Symmetric initial data. ρ ( x ) = η ( x ) = [7 , . −9 −8 −7 −6 −5 −4 −3 −2 Δx −5 −4 −3 −2 −1 e rr o r errorlineΔofΔslopeΔ1 (b) Asymmetric initial data. ρ ( x ) = [5 , and η ( x ) = [10 , . Figure 2.
In the purely diffusive system all interaction kernels are set to zero. Bothgraphs show the convergence to the benchmark solution. The triangular markersdenote the discrepancy between the numerical solution and the benchmark solution.A line of slope one is superimposed for the ease of comparison. On the left we startthe system with symmetric initial conditions, on the right we start with asymmetricinitial data. In both cases the scheme has numerical convergence order 1.6.1.2.
Gaussian cross-interactions.
Next, we add non-local self-interaction and cross-interactions.We choose smooth Gaussians with different variances. These potentials, like the related, moresingular Morse potentials, are classical in mathematical biology since oftentimes the availabilityof sensory information such as sight, smell or hearing is spatially limited [27, 20, 18]. For the intraspecific interaction we use W ( x ) = W ( x ) = 1 − exp (cid:18) − | x | × . (cid:19) , while we choose W ( x ) = − W ( x ) = 1 − exp (cid:18) − | x | × . (cid:19) , for the interspecific interactions.We consider system (1) with the diffusive coefficients ν = 0 . ǫ ∈ { . , . } , and we initialisethe system with ρ ( x ) = η ( x ) = c (cid:0) ( s − . . − s ) (cid:1) + , on the domain [0 , c is such that ρ and η have unit mass. Figure 3 depicts thesimulation with Gaussian kernels as interaction potentials. In Figures 3(a) & 3(b) we present theerror plots corresponding to ǫ ∈ { . , . } . Again we observe convergence to the reference solutionwith a first order accuracy in space. −9 −8 −7 −6 −5 −4 −3 −2 Δx −5 −4 −3 −2 −1 e rr o r errorlineΔofΔslopeΔ1 (a) Convergence to benchmark solution for individualdiffusion constant ǫ = 0 . −9 −8 −7 −6 −5 −4 −3 −2 Δx −5 −4 −3 −2 −1 e rr o r errorlineΔofΔslopeΔ1 (b) Convergence to benchmark solution for individualdiffusion constant ǫ = 0 . Figure 3.
We choose Gaussian interaction kernels of different strengths and rangesfor the self-interaction and the cross-interaction, respectively. The graphs show thenumerical convergence order in the cases of ǫ = 0 . ǫ = 0 . General behaviour of solutions and stationary states.
In this section we aim to studythe asymptotic behaviour of system (1) numerically. Let us begin by going back to the set upof Section 6.1.2. We note that the potentials were chosen in such a way that there is an attrac-tive intraspecific force, and the cross-interactions are chosen as attractive-repulsive explaining thesegregation observed in Figure 4(a). For larger self-diffusivity, ǫ = 0 .
5, we see that some mixingoccurs. In the absence of any individual diffusion we would have expected adjacent species with ajump discontinuity at their shared boundary [19]. However, this phenomenon is no longer possible
ONVERGENCE OF A FINITE VOLUME SCHEME 27 as we have a control on the gradients of each individual species, by Lemma 3.3, rendering jumpsimpossible thus explaining the continuous transition. ρη (a) Stationary state ( ǫ = 0 . ρη (b) Stationary state ( ǫ = 0 . Figure 4.
We choose Gaussian interaction kernels of different strengths and rangesfor the self-interaction and the cross-interaction, respectively. The graphs corre-spond to the simulationed stationary states for self-diffusivities ǫ = 0 .
1, and ǫ = 0 . Case of singular potentials.
In this section we go beyond the limit of what we could provein this paper. On the one hand we consider more singular potentials and on the other hand weconsider vanishing individual diffusion. We study system (1) for ǫ ∈ { , . , . , . , . } , and ν ∈ { . , . } . Here the potentials are given by W ( x ) = W ( x ) = x / . for the self-interaction terms and W ( x ) = | x | = ± W ( x ) . for the cross-interactions. The system is posed on the domain [0 ,
5] with a grid size of ∆ x = 2 − .Note that the case of ǫ = 0 corresponds to the absence of individual diffusion, see Figures 5(a) &5(c). By virtue of Corollary 3.3, it is the individual diffusion that regularises the stationary states, inthe sense that we will not observe any discontinuities in either ρ or η . As we add individual diffusionwe can see the immediate regularisation. While stationary states may still remain segregated, as itis shown in Figure 5(b), adjacent solutions are not possible anymore (see Figure 5(d)). In the caseof attractive-repulsive interspecific interactions, i.e. W ( x ) = | x | = − W ( x ) . analyticalρη (a) In the case ν = 0 . ǫ = 0, we obtain a great agree-ment of the numerically computed stationary states andthe analytical stationary states described in [19]. analyticalε = 0.00ε = 0.04ε = 0.09 (b) Adding individual diffusion may still lead to segre-gated stationary states. However both species remaincontinuous as they mix. ρηanalytical (c) The case ν = 0 . ǫ = 0 leads to adjacent station-ary states. Again we see an excellent agreement of thenumerical stationary states and the analytical ones [19]. analyticalε = 0.00ε = 0.02ε = 0.04ε = 0.06 (d) The regularising effect of the individual diffusion, byCorollary 3.3, becomes apparent immediately. Instanta-neously both species become continuous as they start tointermingle in a small region. This region grows as wekeep increasing the individual diffusion. Figure 5.
We pushed our numerical scheme to see how it performs in regimes inwhich we are unable to prove convergence. We chose Newtonian cross-interactionsin the attractive-repulsive case. The red curves denote the symmetric stationarystates of η while the blue curves are the stationary distributions of ρ . The differentline widths and styles correspond to varying values of ǫ .we expect both species to segregate [19]. We initialise the system with the following symmetricinitial data ρ ( x ) = η ( x ) = c (cid:0) ( x − − x ) (cid:1) + , ONVERGENCE OF A FINITE VOLUME SCHEME 29 as symmetric initial data are known to approach stationary states [19]. Here the constant c nor-malises the mass of ρ and η to one.Figures 5(a) & 5(c) show our scheme performs well even in regimes we are unable to show con-vergence due to the lack of regularity in the potentials as well as the lack of regularity due tothe absence of the porous medium type self-diffusion. While the schemes developed in [14, 19]are asymptotic preserving, their convergence to weak solutions of the respective equations couldnot be established. We reproduce the steady states of [19] that exhibit phase separation phenom-ena. Figure 6 displays the stationary state in the case ν = 0 .
09 and ǫ = 0. Even in the case ofno regularising individual diffusion and with Newtonian cross-interactions we observe a numericalconvergence order of one. ρη (a) Strictly segregated stationary state in the absence ofindividual diffusion, ǫ = 0, and Newtonian potentials. −8 −7 −6 −5 −4 −3 −2 Δx −3 −2 −1 e rr o r errorlineΔofΔslopeΔ1 (b) Convergence to benchmark solution in the case of ǫ = 0 and Newtonian attractive-repulsive potentials. Figure 6.
Stationary state and numerical convergence order in the case of New-tonian attractive-repulsive cross-interactions. Even though our estimates fail in theanalysis above we observe a numerical convergence order of one.In the case of attractive-attractive self-interactions, i.e. W ( x ) = | x | = W ( x ), we observe aninteresting phenomenon. Even in the absence of the individual diffusion, i.e. ǫ = 0, some additionalmixing occurs even though we expect sharp boundaries, see Figure 7, due to numerical diffusion.This is in contrast to the finite volume schemes proposed in [14, 19].7. Conclusion
In this paper we presented a finite volume scheme for a system of non-local partial differentialequations with cross-diffusion. We were able to reproduce a continuous energy estimate on thediscrete level for our scheme. These discrete estimates for the approximate solution are enoughto get compactness results and we are able to identify the limit of the approximate solutions asa weak solution of the equation. We complement the analytical part with numerical simulations.These back up our convergence result and we are also able to apply the scheme in cases in whichwe cannot show convergence. To this end we pushed the scheme to regimes of singular potentialsalso lacking the regularising porous medium type self-diffusion terms. Comparing them with the analyticalρη
Figure 7.
We choose m = 0 . m = 0 . Acknowledgements
JAC was partially supported by EPSRC grant number EP/P031587/1.
References [1] H. Attouch, G. Buttazzo, and G. Michaille.
Variational analysis in Sobolev and BV spaces: applications to PDEsand optimization . SIAM, 2014.[2] J. Berendsen, M. Burger, and J.-F. Pietschmann. On a cross-diffusion model for multiple species with nonlocalinteraction and size exclusion.
Nonlinear Analysis , 159:10–39, 2017.[3] M. Bertsch, R. Dal Passo, and M. Mimura. A free boundary problem arising in a simplified tumour growthmodel of contact inhibition.
Interfaces and Free Boundaries , 12(2):235–250, 2010.[4] M. Bertsch, M. Gurtin, and H. D. On a degenerate diffusion equation of the form c ( z ) t = ϕ ( z x ) x with applicationto population dynamics. Journal of differential equations , 67(1):56–89, 1987.[5] M. Bertsch, M. Gurtin, H. D., and L. Peletier. On interacting populations that disperse to avoid crowding:preservation of segregation.
Journal of mathematical biology , 23(1):1–13, 1985.[6] M. Bertsch, M. Gurtin, and D. Hilhorst. On interacting populations that disperse to avoid crowding: the caseof equal dispersal velocities.
Nonlinear Analysis: Theory, Methods & Applications , 11(4):493–499, 1987.[7] M. Bertsch, D. Hilhorst, H. Izuhara, and M. Mimura. A nonlinear parabolic-hyperbolic system for contactinhibition of cell-growth.
Differ. Equ. Appl , 4(1):137–157, 2012.[8] M. Bessemoulin-Chatard, C. Chainais-Hillairet, and F. Filbet. On discrete functional inequalities for some finitevolume schemes.
IMA J. Num. Analysis , 35(3):1125–1149, 2015.[9] M. Bessemoulin-Chatard and F. Filbet. A finite volume scheme for nonlinear degenerate parabolic equations.
SIAM Journal on Scientific Computing , 34(5):B559–B583, 2012.[10] M. Bruna, M. Burger, H. Ranetbauer, and M.-T. Wolfram. Cross-diffusion systems with excluded-volume effectsand asymptotic gradient flow structures.
Journal of Nonlinear Science , 27(2):687–719, 2017.[11] M. Burger, M. Di Francesco, S. Fagioli, and A. Stevens. Sorting phenomena in a mathematical model for twomutually attracting/repelling species. preprint arXiv:1704.04179 , 2017.
ONVERGENCE OF A FINITE VOLUME SCHEME 31 [12] V. Calvez and J. A. Carrillo. Volume effects in the Keller-Segel model: energy estimates preventing blow-up.
J.Math. Pures Appl. (9) , 86(2):155–175, 2006.[13] E. A. Carl. Population control in arctic ground squirrels.
Ecology , 52(3):395–413, 1971.[14] J. A. Carrillo, A. Chertock, and Y. Huang. A finite-volume method for nonlinear nonlocal equations with agradient flow structure.
Communications in Computational Physics , 17(01):233–258, 2015.[15] J. A. Carrillo, Y.-P. Choi, and M. Hauray. The derivation of swarming models: mean-field limit and wassersteindistances. In
Collective dynamics from bacteria to crowds , pages 1–46. Springer, 2014.[16] J. A. Carrillo, A. Colombi, and M. Scianna. Adhesion and volume constraints via nonlocal interactions lead tocell sorting. preprint arXiv:1706.08969 , 2017.[17] J. A. Carrillo, S. Fagioli, F. Santambrogio, and M. Schmidtchen. Splitting schemes & segregation in reaction-(cross-)diffusion systems. arXiv preprint arXiv:1711.05434 , 2017.[18] J. A. Carrillo, Y. Huang, and S. Martin. Explicit flock solutions for Quasi-Morse potentials.
European J. Appl.Math. , 25(5):553–578, 2014.[19] J. A. Carrillo, Y. Huang, and M. Schmidtchen. Zoology of a non-local cross-diffusion model for two species. arXiv preprint arXiv:1705.03320 , 2017.[20] J. A. Carrillo, S. Martin, and V. Panferov. A new interaction potential for swarming models.
Phys. D , 260:112–126, 2013.[21] J. A. Carrillo, R. J. McCann, C. Villani, et al. Kinetic equilibration rates for granular media and related equations:entropy dissipation and mass transportation estimates.
Revista Matematica Iberoamericana , 19(3):971–1018,2003.[22] C. Chainais-Hillairet and F. Filbet. Asymptotic behaviour of a finite-volume scheme for the transient drift-diffusion model.
IMA journal of numerical analysis , 27(4):689–716, 2007.[23] X. Chen, A. J¨ungel, and J.-G. Liu. A note on Aubin-Lions-Dubinski˘ı lemmas.
Acta Appl. Math. , 133:33–43,2014.[24] M. Di Francesco, A. Esposito, and S. Fagioli. Nonlinear degenerate cross-diffusion systems with nonlocal inter-action.
Nonlinear Analysis , 169:94–117, 2018.[25] M. Di Francesco and S. Fagioli. Measure solutions for non-local interaction pdes with two species.
Nonlinearity ,26(10):2777, 2013.[26] P. Domschke, D. Trucu, A. Gerisch, and M. A. J. Chaplain. Mathematical modelling of cancer invasion: impli-cations of cell adhesion variability for tumour infiltrative growth patterns.
J. Theoret. Biol. , 361:41–60, 2014.[27] M. R. D’Orsogna, Y.-L. Chuang, A. L. Bertozzi, and L. S. Chayes. Self-propelled particles with soft-core inter-actions: patterns, stability, and collapse.
Physical review letters , 96(10):104302, 2006.[28] R. Eymard, T. Gallou¨et, and R. Herbin. Finite volume methods.
Handbook of numerical analysis , 7:713–1018,2000.[29] R. Eymard, T. Gallou¨et, and R. Herbin. Finite volume methods. In
Handbook of numerical analysis, Vol. VII ,Handb. Numer. Anal., VII, pages 713–1020. North-Holland, Amsterdam, 2000.[30] F. Filbet and M. Herda. A finite volume scheme for boundary-driven convection-diffusion equations with relativeentropy structure.
Numer. Math. , 137(3):535–577, 2017.[31] A. Gerisch and M. A. J. Chaplain. Mathematical modelling of cancer cell invasion of tissue: Local and non-localmodels and the effect of adhesion.
Journal of Theoretical Biology , 250(4):684–704, 2008.[32] W. Gurney and R. Nisbet. The regulation of inhomogeneous populations.
Journal of Theoretical Biology ,52(2):441–457, 1975.[33] M. E. Gurtin and A. Pipkin. A note on interacting populations that disperse to avoid crowding.
Quarterly ofApplied Mathematics , 42(1):87–94, 1984.[34] Y. Itˆo. The growth form of populations in some aphids, with special reference to the relation between populationdensity and the movements.
Researches on Population Ecology , 1(1):36–48, 1952.[35] A. Mogilner and L. Edelstein-Keshet. A non-local model for a swarm.
Journal of Mathematical Biology , 38(6):534–570, 1999.[36] M. Morisita. Population density and dispersal of a water strider. gerris lacustris: Observations and considerationson animal aggregations.
Contributions on Physiology and Ecology, Kyoto University , 65:1–149, 1950.[37] M. Morisita. Dispersion and population pressure: experimental studies on the population density of an ant-lion,glenuroides japonicus m’l (2).
Jap. J. Ecol , 4(71):9, 1954. [38] H. Murakawa and H. Togashi. Continuous models for cell–cell adhesion.
Journal of theoretical biology , 374:1–12,2015.[39] K. J. Painter, J. M. Bloomfield, J. A. Sherratt, and A. Gerisch. A nonlocal model for contact attraction andrepulsion in heterogeneous cell populations.
Bull. Math. Biol. , 77(6):1132–1165, 2015.[40] K. J. Painter and T. Hillen. Volume-filling and quorum-sensing in models for chemosensitive movement.
Can.Appl. Math. Q. , 10(4):501–543, 2002.[41] C. M. Topaz, A. L. Bertozzi, and M. A. Lewis. A nonlocal continuum model for biological aggregation.
Bulletinof mathematical biology , 68(7):1601–1623, 2006.[42] A. Volkening and B. Sandstede. Modelling stripe formation in zebrafish: an agent-based approach.
Journal ofThe Royal Society Interface , 12(112), 2015.
E-mail address : [email protected] E-mail address : [email protected] E-mail address : [email protected] .0 0.5 1.0 1.5 2.0 2.5 3.00.000.250.500.751.001.251.501.75.0 0.5 1.0 1.5 2.0 2.5 3.00.000.250.500.751.001.251.501.75