A numerical approach for a general class of reaction-diffusion systems

# A numerical approach for a general class of the spatial segregation of reaction-diffusion systems arising in population dynamics

Avetik Arakelyan Institute of Mathematics, National Academy of Sciences of Armenia, 0019 Yerevan, Armenia  and  Rafayel Barkhudaryan Institute of Mathematics, National Academy of Sciences of Armenia and American University of Armenia, 0019 Yerevan, Armenia
###### Abstract.

In the current work we consider the numerical solutions of equations of stationary states for a general class of the spatial segregation of reaction-diffusion systems with population densities. We introduce a discrete multi-phase minimization problem related to the segregation problem, which allows to prove the existence and uniqueness of the corresponding finite difference scheme. Based on that scheme, we suggest an iterative algorithm and show its consistency and stability. For the special case we show that the problem gives rise to the generalized version of the so-called two-phase obstacle problem. In this particular case we introduce the notion of viscosity solutions and prove convergence of the difference scheme to the unique viscosity solution. At the end of the paper we present computational tests, for different internal dynamics, and discuss numerical results.

###### Key words and phrases:
Free boundary, Two-phase obstacle problem, Reaction-diffusion systems, Finite difference, Viscosity solution
###### 2000 Mathematics Subject Classification:
35R35, 65N06, 65N22, 92D25

## 1. Introduction and known results

### 1.1. The setting of the problem

In recent years there have been intense studies of spatial segregation for reaction-diffusion systems. The existence of spatially inhomogeneous solutions for competition models of Lotka-Volterra type in the case of two and more competing densities have been considered in [11, 12, 13, 15, 21, 20, 25]. The aim of this paper is to study the numerical solutions for a certain class of the spatial segregation of reaction-diffusion system with population densities.

Let be a connected and bounded domain with smooth boundary and be a fixed integer. We consider the steady-states of competing species coexisting in the same area . Let denotes the population density of the component with the internal dynamic prescribed by .

We call the -tuple a segregated state if

 ui(x)⋅uj(x)=0, a.e.  for i≠j, x∈Ω.

The problem amounts to

 (1) Minimize E(u1,⋯,um)=∫Ωm∑i=1(12|∇ui(x)|2+Fi(x,ui(x)))dx

over the set

 S={(u1,…,um)∈(H1(Ω))m:ui≥0,ui⋅uj=0,ui=ϕion∂Ω},

where for and on the boundary .

We assume that

 Fi(x,s)=∫s0fi(x,v)dv,

where is Lipschitz continuous in uniformly continuous in and .

###### Remark 1.

Functions ’s are defined only for non negative values of s (recall that our densities ’s are assumed non negative); thus we can arbitrarily define such functions on the negative semiaxis. For the sake of convenience, when , we will let . This extension preserves the continuity due to the conditions on defined above. In the same way, each is extended as an even function.

###### Remark 2.

We emphasize that for the case the assumption is that for all the functions are nonnegative and uniformly continuous in . Also for simplicity, throughout the paper we shall call both and as internal dynamics.

###### Remark 3.

We would like to point out that the only difference between our minimization problem (1) and the problem discussed by Conti,Terrachini and Verzini [12], is the sign in front of the internal dynamics . In our case, the plus sign of allows to get rid of some additional conditions, which are imposed in [12, Section ]. Those conditions are important to provide coercivity of a minimizing functional in [12]. But in our case the above given conditions together with convexity assumption on with respect to the variable are enough to conclude which in turn implies coercivity of a functional (1).

In order to speak on the local properties of the population densities, let us introduce the notion of multiplicity of a point in .

###### Definition 1.

The multiplicity of the point is defined by:

 m(x)=card{i:measure(Ωi∩B(x,r))>0,∀r>0},

where .

For the local properties of the same results as in [12] with the opposite sign in front of the internal dynamics hold. Below, for the sake of clarity, we write down those results from [12] with appropriate changes.

###### Lemma 1.

(Proposition in [12]) Assume that then the following holds:

1. If then there exists such that for every ;

 ui≡0onB(x0,r).
2. If then there are and such that in

 Δui=fi(x,ui), uj≡0for j≠i.
3. If then there are and such that for every and , we have and in

 Δ(ui−uj)=fi(x,(ui−uj))χ{ui>uj}−fj(x,−(ui−uj))χ{ui

Next, we state the following uniqueness Theorem due to Conti, Terrachini and Verzini.

###### Theorem 1 (Theorem 4.2 in [12]).

Let the functional in minimization problem (1) be coercive and moreover each be convex in the variable for all . Then, the problem (1) has a unique minimizer.

This theorem will play a crucial role in studying the difference scheme, especially for the case where we will reformulate it as a generalized two-phase obstacle problem. Note that in this case, the problem will be reduced to:

 (2) Minimize E(u1,u2)=∫Ω2∑i=1(12|∇ui(x)|2+Fi(x,ui(x)))dx,

over the set

 S={(u1,u2)∈(H1(Ω))2:ui≥0,u1⋅u2=0,ui=ϕion∂Ω}.

Here with property on the boundary .

### 1.2. Known results

In last years there has been much interest given to study the numerical approximation of reaction-diffusion type equations. For instance, the equations arising in the study of population ecology when high competitive interactions between different species occurs.

We refer the reader to [12, 16, 17, 18, 19, 20, 21] for models involving Dirichlet boundary data. A complete analysis of the stationary case has been studied in [12]. Also numerical simulation for the spatial segregation limit of two diffusive Lotka-Volterra models in presence of strong competition and inhomogeneous Dirichlet boundary conditions is provided in [26].

In the work [6] Bozorgnia proposed two numerical algorithms for the problem (1) with the internal dynamics . The finite element approximation is based on the local properties of the solution. In this case the author was able to provide the convergence of the method. The second approach is a finite difference method, but lack of its analysis in [6]. This finite difference method has been generalized in [9] for the case of non-negative . In [9] the authors present a numerical consistent variational system with strong interaction, and provide disjointness condition of populations during the iteration of the scheme. In this case the proposed algorithm is lack of deep analysis, especially for the case of three and more competing populations.

The present work concerns to close that gap and provides theoretical results for finite difference scheme, with competing populations and general internal dynamics satisfying certain suitable conditions. We introduce the discrete analogue of minimization problem and prove the existence and uniqueness of the difference scheme. Moreover for the special case we introduce viscosity solution and prove the convergence of corresponding difference scheme.

### 1.3. Notations

We will make the notations for the one-dimensional and two-dimensional cases parallely. For the sake of simplicity, we will assume that in one-dimensional case and in two-dimensional case in the rest of the paper, keeping in mind that the method works also for more complicated domains.

Let be a positive integer, and

 xi=−1+ih,yi=−1+ih,i=0,1,...,N.

We use the notation and (or simply , where is one- or two-dimensional index) for the finite difference scheme approximation to and , respectively. Concerning the boundary functions we assume they are extended to be zero everywhere outside the boundary for all . The discrete approximation for these functions will be and respectively (or simply , where is one- or two-dimensional index). Note that for the case with two population densities, we additionally will use the following notations:

 gi=ϕ1i−ϕ2i=ϕ1(xi)−ϕ2(xi)

and

 gi,j=ϕ1i,j−ϕ2i,j=ϕ1(xi,yj)−ϕ2(xi,yj),

in one- and two-dimensional cases, respectively.

In this paper we will use also notations , (not to be confused with functions ).

Denote

 N={i: 0≤i≤N}orN={(i,j): 0≤i,j≤N},
 No={i: 1≤i≤N−1}orNo={(i,j): 1≤i,j≤N−1},

in one- and two- dimensional cases, respectively, and

 ∂N=N∖No.

In one-dimensional case, we consider the following approximation for Laplace operator: for any ,

 Δhvi≡Lhvi=vi−1−2vi+vi+1h2,

and for two-dimensional case we introduce the following 5-point stencil approximation for Laplacian:

 Δhvi,j≡Lhvi,j=vi−1,j+vi+1,j−4vi,j+vi,j−1+vi,j+1h2

for any .

## 2. Segregation problem with m≥2 population densities

### 2.1. The minimizing functional and existence of difference scheme

In this section, we introduce the discrete counterpart of the spatial segregation problem for the general case. As in the previous section we assume that are convex in the variable for all and satisfy the properties stated in introduction.

In the rest of the paper the following notation

 ^zk=zk−∑j≠kzj,

for elements , will play a crucial role. We focus on the following functional

 (3) Jh(v1,v2,…,vm)=−12m∑l=1(Lh^vl,vl)+m∑l=1(∑α∈NFl(xα,vlα))−m∑l=1(Lh^ϕl,vl),

defined over the set

 (4) S={(v1,v2,…,vm)∈(H)m:vpα≥0,vpα⋅vqα=0,p≠q,vpα=0on∂N},

where

 H={v=(vα):vα∈R, α∈N}.

Here for and , , the inner product is defined by

 (w,v)=∑α∈N0wα⋅vα.
###### Theorem 2.

In view of definitions (3) and (4), the following minimization problem

 (5) infSJh(v1,v2,…,vm)

has a solution.

###### Proof.

We are going to prove that functional (3) is coercive in the set . Due to the standard arguments of calculus of variations coercivity and lower semi-continuity will imply the existence of minimizers over the closed set . To this end, we observe that

 −12m∑l=1(Lh^vl,vl)=−12m∑l=1(Lhvl,vl)+∑1≤i

For every fixed and such that we have . Thus,

 (Lhvi,vj)=∑{vjα>0}Lhviα⋅vjα≥0

due to the simple fact that implies , which in turn yields . Therefore

 −12m∑l=1(Lh^vl,vl)≥−12m∑l=1(Lhvl,vl)≥C⋅m∑l=1(vl)2

for some constant . Recalling that we finally obtain that the functional (3) is coercive. ∎

###### Proposition 1.

If an element solves the following minimization problem:

 infSJh(v1,v2,…,vm),

then for every and

 (6) Lh(^ulα+^ϕlα)=fl(xα,ulα),wheneverulα>0.
###### Proof.

First of all, it is easy to verify that, if two vectors and belong to the set then for arbitrary we have

 ((^u1+ε^v1)+,(^u2+ε^v2)+,…,(^um+ε^vm)+)∈S,

and

 ((^u1−ε^v1)+,(^u2−ε^v2)+,…,(^um−ε^vm)+)∈S.

Let a vector be a minimizer to our discrete minimization problem (5). If we obtain:

 Jh(w1,w2,…,wm)−Jh(u1,u2,…,um)≥0,

for every . Thus,

 −12(m∑l=1(Lh^wl,wl)−m∑l=1(Lh^ul,ul))++m∑l=1(∑α∈N(Fl(xα,wlα)−Fl(xα,ulα)))−m∑l=1(Lh^ϕl,wl−ul)≥0.

We choose the vector such that for every and the following condition holds

 ^ulα⋅^vlα≥0.

This implies the following identity:

 wlα=(^ulα+ε^vlα)+=(^ulα)++ε(^vlα)+=ulα+εvlα.

Hence,

 (7) −12(m∑l=1∑α∈NLh(^ulα+ε^vlα)⋅(ulα+εvlα)−m∑l=1(Lh^ul,ul))++m∑l=1(∑α∈N(Fl(xα,ulα+εvlα)−Fl(xα,ulα)))−εm∑l=1(Lh^ϕl,vl)=
 (8) =−12(m∑l=1∑α∈NLh^ulα⋅ulα+ε2m∑l=1∑α∈NLh^vlα⋅vlα++εm∑l=1(∑α∈N(Lh^ulα⋅vlα+Lh^vlα⋅ulα+2Lh^ϕlα⋅vlα))−m∑l=1(Lh^ul,ul))++m∑l=1(∑α∈N(Fl(xα,ulα+εvlα)−Fl(xα,ulα)))≥0.

Next, dividing both sides in (8) by and letting we arrive at:

 (9) −12m∑l=1∑α∈N(Lh^ulα⋅vlα+Lh^vlα⋅ulα+2Lh^ϕlα⋅vlα)+m∑l=1(∑α∈Nfl(xα,ulα)⋅vlα)≥0.

Let . We set as follows:

 (10) ⎧⎪ ⎪⎨⎪ ⎪⎩vlα=0,for alll≠l0,α∈N,vl0α0=ul0α0,vl0α=0,for allα∈N∖α0.

It is easy to see, that the chosen vector satisfies . Therefore, the inequality (9) holds for this vector, which means it can be substituted into (9). We clearly obtain

 −12m∑l=1∑α∈N(Lh^ulα⋅vlα+Lh^vlα⋅ulα+2Lh^ϕlα⋅vlα)=−Lh(^ul0α0+^ϕl0α0)⋅ul0α0

and

 m∑l=1(∑α∈Nfl(xα,ulα)⋅vlα)=fl0(xα0,ul0α0)⋅ul0α0.

Hence,

 (11) −Lh(^ul0α0+^ϕl0α0)+fl0(xα0,ul0α0)≥0.

In the same way, for every and we apparently have:

 Jh((^u1−ε^v1)+,(^u2−ε^v2)+,…,(^um−ε^vm)+)−Jh(u1,u2,…,um)≥0.

In this case we choose the vector such that for every and the following condition holds

 ^vlα⋅(^ulα−^vlα)≥0.

This implies the following identity:

 (^ulα−ε^vlα)+=(^ulα)+−ε(^vlα)+=ulα−εvlα.

After proceeding the same steps as above, we obtain

 (12) 12m∑l=1∑α∈N(Lh^ulα⋅vlα+Lh^vlα⋅ulα+2Lh^ϕlα⋅vlα)−m∑l=1(∑α∈Nfl(xα,ulα)⋅vlα)≥0.

Here again we choose as in (10), which apparently satisfies . Therefore, we can substitute the vector into (12), which will lead to

 (13) Lh(^ul0α0+^ϕl0α0)−fl0(xα0,ul0α0)≥0.

Thus, in light of (11) and (13) we obtain

 (14) Lh(^ul0α0+^ϕl0α0)=fl0(xα0,ul0α0),wheneverul0α0>0.

###### Lemma 2.

Let an element solves the minimization problem (5), and moreover assume that for some and . Then

 Lh(^ulα+^ϕlα)≤fl(xα,ulα).
###### Proof.

Assume that for some fixed and . There are two possibilities: either or . Let then we take a vector as follows:

 (15) ⎧⎪ ⎪⎨⎪ ⎪⎩vlα=0,for alll≠l0,α∈N,vl0α0=1,vl0α=0,for allα∈N∖α0.

It is clear that is satisfied for every and . Hence, one can substitute the vector (15) into (9). This implies the inequality (11), namely

 Lh(^ul0α0+^ϕl0α0)≤fl0(xα0,ul0α0).

Now, if we assume that then there exists some such that . In this case according to Proposition 1, we have

 Lh(^uq0α0+^ϕq0α0)=fq0(xα0,uq0α0)≥0.

Thus,

 Lh(^ul0α0+^ϕl0α0)=−Lh(^uq0α0+^ϕq0α0)−2∑r≠l0,q0Lh(urα0+ϕrα0)≤0≤fl0(xα0,ul0α0).

Next theorem shows the existence of the corresponding finite difference scheme for two or more population densities. Note that for the particular case we obtain the difference scheme for the so-called multi-phase obstacle problem.

###### Theorem 3.

If an element is a minimizer to (5), then the vector where solves the discrete system:

 (16) ⎧⎪ ⎪⎨⎪ ⎪⎩wlα=max(−fl(xα,wlα)h24+¯¯¯¯wlα−∑p≠l¯¯¯¯wpα,0),α∈Nowlα=ϕlα,α∈∂N

for every and . Here for a given uniform mesh on we define

 ¯¯¯¯wlα=14[wl(xi−1,yj)+wl(xi+1,yj)+wl(xi,yj−1)+wl(xi,yj+1)]

as the average of for all neighbour points of .

###### Proof.

If solves the discrete system (16), then one can easily see that for every and the following property holds:

 zpα⋅zqα=0,wheneverp≠q,

which in turn implies that belongs to the set . Thus the solution to the discrete system (16) itself implies the disjointness property.

Now, in view of Proposition 1 and Lemma 2, we clearly infer that for every minimizer to (5), the vector satisfies the system (16). ∎

###### Remark 4.

We would like to emphasize that the theorem statement holds also for more general uniformly discretized schemes without restriction on the number of stencil (neighbour) points corresponding to a discrete Laplacian.

Observe that we have only used to prove the existence of the solution to the discrete system (16). The convexity assumption on will be needed to prove the uniqueness of the discrete scheme which is done in Section 2.2.

### 2.2. Uniqueness of difference scheme

In this section our goal is to show the uniqueness of the difference scheme, which solves the discrete system (16). To this aim, we need two auxiliary lemmas.

For the sake of convenience we denote by the set of all closest neighbour points corresponding to a mesh point .

###### Lemma 3.

Let the functions be nondecreasing with respect to the variable . We take any two elements and in . If and are satisfying the discrete system (16), then the following equation holds:

 maxN(^ulα−^vlα)=max{ulα≤vlα}(^ulα−^vlα),

for all .

###### Proof.

We argue by contradiction. Suppose for some we have

 (17) ^ul0α0−^vl0α0=maxN(^ul0α−^vl0α)=max{ul0α>vl0α}(^ul0α−^vl0α)>max{ul0α≤vl0α}(^ul0α−^vl0α).

Then taking into account the following simple chain of inclusions

 (18) {ulα>vlα}⊂{^ulα>^vlα}⊂{ulα≥vlα},

we obviously see that implies . On the other hand, the discrete system (16) gives us

 Lh(^ul0α0+^ϕl0α0)=fl0(xα0,ul0α0)andLh(^vl0α0+^ϕl0α0)≤fl0(xα0,vl0α0).

Therefore

 Lh(^ul0α0−^vl0α0)≥fl0(xα0,ul0α0)−fl0(xα0,vl0α0)≥0.

Thus,

 ^ul0α0−^vl0α0≤14∑{δ∈nbr(α0)}(^ul0δ−^vl0δ),

which implies that