# Numerical Solution of a Class of Predator-Prey Systems with Complex Dynamics Characters Based on a Sinc Function Interpolation Collocation Method

Authors: Mingjing Du; Pengfei Ning; Yulan Wang
Journal: Complexity
 (2020)

---

## Abstract

Although many kinds of numerical methods have been announced for the predator-prey system, simple and efficient methods have always been the direction that scholars strive to pursue. Based on this problem, in this paper, a new interpolation collocation method is proposed for a class of predator-prey systems with complex dynamics characters. Some complex dynamics characters and pattern formations are shown by using this new approach, and the results have a good agreement with theoretical results. Simulation results show the effectiveness of the method.


---

## Body

## 1. Introduction

Mathematical models have played a significant role in understanding the dynamics of populations for over a century. These models have helped abstract key features of interactions between biological organisms that have given insight into a variety of observed phenomena in real populations.In [1], let ηt and ut be the population size of a prey and predator species, respectively, and suppose that these functions obey the Gause-type model:(1)∂η∂t=rη1−ηk−θηη+au,∂u∂t=I−ru−ψuu+cu.In [2], the authors give a predator-prey population dynamics. This model includes all three components (predator, prey, and subsidy) in a single spatial location:(2)∂η∂t=rη1−ηk−θηη+av,∂u∂t=I−ru−ψuu+cv,∂v∂t=εθηη+av+bψuu+cv−δv.Turing [3] proposed a dynamical mechanism, which has been extensively used to explain how Turing patterns are formed, and it is now known as Turing bifurcation or Turing instability. Reaction-diffusion systems have attracted increasing attention from the mathematical biologists in recent years to seek insights into the fascinating patterns that occur in living organisms and in ecological systems. As suggested in [2], a more realistic system could account for continuous spatial variation in the populations. In[4, 5], the authors assumed that these populations disperse randomly over a spatial region. This leads to the following set of predator-prey systems with diffusion:(3)∂η∂t=d1Δη+rη1−ηk−θηη+av,∂u∂t=d2Δu+I−ru−ψuu+cv,∂v∂t=d3Δv+εθηη+av+bψuu+cv−δv,where η=ηx,y,t and v=vx,y,t are the prey and predator populations, u=ux,y,t is the quantity of subsidy, t is time, r,k,θ,a,I,ψ,c,ε,η>0 are all nonnegative kinetic parameters, and d1, d2, and d3 are the positive diffusion coefficients for the three populations; ∇=∂2/∂x2+∂2/∂y2.The sufficient condition of locally asymptotically stability of the equilibrium has been obtained in [4–6]. Turing bifurcation occurs when the equilibrium state that is stable in absence of nondiffusion becomes unstable in presence of cross-diffusion.In the present paper, we consider the following reaction-diffusion system [7–9], which is described as(4)∂η∂t=d1Δη+a11η+a12u+a13v+f1η,u,v,∂u∂t=d2Δu+a21η+a22u+a23v+f2η,u,v,∂v∂t=d3Δv+a31η+a32u+a33v+f3η,u,v,where ηx,y,t,ux,y,t, and vx,y,t are the unknown functions. x,y∈Ω=a,b×c,d, t>0, the smooth boundary is ∂Ω, and homogeneous Neumann boundary condition, namely, ∂u/∂n∂Ω=∂v/∂n∂Ω=0. The parameters d1, d2, and d3 are the diffusion coefficients, and Δ=∂2/∂x2+∂2/∂y2 is the Laplacian operator, respectively. Parameter aij is the known parameter. fiη,u,v are the known nonlinear functions of η,u,andv.In recent years, reaction-diffusion systems [10, 11] have attracted increasing attention from the mathematical biologists to seek insights into the fascinating patterns that occur in living organisms and in ecological systems [5, 12, 13].In general, the exact solution of system (4) cannot be obtained, and the approximate solution must be obtained by numerical calculation. Although many kinds of numerical methods of the nonlinear reaction-diffusion system have been announced, such as finite difference method [14], B-spline method [15], finite element method [16, 17], spectral method [18–20], the perturbation method and variational iteration method [21, 22], barycentric interpolation collocation method [23–26], and reproducing kernel method [27, 28], this paper investigates some nonlinear diffusion predator-prey systems [5, 12, 13] based on a new interpolation collocation method, and the model (21) is adopted as an example to elucidate the solution process.


## 2. Bifurcation Analysis of the System

In this section, we give the Turing bifurcation conditions of system (4). We assume an equilibrium point of nondiffusive system (4) as E∗=η0,u0,v0, so(5)a11η0+a12u0+a13v0+f1η0,u0,v0=0,a21η0+a22u0+a23v0+f2η0,u0,v0=0,a31η0+a32u0+a33v0+f3η0,u0,v0=0.Now, we do linear stability analysis of the equilibrium pointE∗. We evaluate the Jacobian matrix A0 of the system at E∗ as(6)A0=a11+∂f1∂ηη0,u0,v0a12+∂f1∂uη0,u0,v0a13+∂f1∂vη0,u0,v0a21+∂f2∂ηη0,u0,v0a22+∂f2∂uη0,u0,v0a23+∂f2∂vη0,u0,v0a31+∂f3∂ηη0,u0,v0a32+∂f3∂uη0,u0,v0a33+∂f3∂vη0,u0,v0.In general, if real part of eigenvalues ofA0 is negative, then nondiffusive system (4) is stable. Next, we derive diffusion-driven instability conditions and show that these are a generalization of the classical diffusion-driven instability conditions in the absence of diffusion. The Jacobian matrix Ak of system (4) is(7)Ak=a11+∂f1∂ηη0,u0,v0a12+∂f1∂uη0,u0,v0a13+∂f1∂vη0,u0,v0a21+∂f2∂ηη0,u0,v0a22+∂f2∂uη0,u0,v0a23+∂f2∂vη0,u0,v0a31+∂f3∂ηη0,u0,v0a32+∂f3∂uη0,u0,v0a33+∂f3∂vη0,u0,v0−k2d1000d2000d3.Turing bifurcation occurs when the equilibrium state that is stable in absence of nondiffusion becomes unstable in presence of cross-diffusion. Thus, the only way for theE∗ to become an unstable point of cross-diffusion system (4) is that the real part of eigenvalues Ak is positive; then, diffusion system (4) is unstable.


## 3. Description of the Method

### 3.1. Periodic Sinc Function

Sinc functions are used in different areas of physics and mathematics. A periodic sinc functionSN is defined as(8)SNx=sinπx/h2π/htanx/2,where h=2π/N and SN is the interpolation function of periodic δ function.It can be proved thatSNxx⟶0=1, and the first-order derivative of function SNx at xj=jh,j=1,2,…,N is(9)SN′xj=0,j%N=0,−1jcotjh/22,j%N≠0.The two-order derivative of functionSNx at xj=jh,j=1,2,…,N is(10)SN″xj=−π23h2−16,j%N=0,−−1jcsc2jh/22,j%N≠0.Givenh>0, we define the following interpolation space:(11)SpanSNx−jh,j=1,2,…,N.LetIN be the interpolation operator for any function ux defined on 0,2π. The interpolation function INux of sequence u1,u2,…,uN can be written as(12)ux∼INux=∑m=1NumSNx−xm,where uj=ujh,j=1,…,N. It is not difficult to derive simple expressions of the nth derivatives of INux at xj=jh:(13)INunxj=∑m=1NumSNnxj−xm,where xj−xm=j−mh. INunxj can be written in the matrix form as follows:(14)INunx1INunx2INunx3⋮INunxN=SNnxNSNnxN−1SNnxN−2…SNnx1SNnx1SNnxNSNnxN−1…SNnx2SNnx2SNnx1SNnxN…SNnx3……………SNnxN−1SNnxN−2SNnxN−3…SNnxNu1u2u3⋮uN.Noting thenth differential matrix DNn:(15)DNn=SNnxNSNnxN−1SNnxN−2…SNnx1SNnx1SNnxNSNnxN−1…SNnx2SNnx2SNnx1SNnxN…SNnx3……………SNnxN−1SNnxN−2SNnxN−3…SNnxN,where the 1st differential matrix is as follows:(16)DN1=0⋱⋮⋱−cot1h/22−cot1h/22⋱cot2h/22⋱cot2h/22cot2h/22⋱−cot1h/22⋱−cot3h/22⋮⋱0⋱⋮−cot3h/22⋱−cot1h/22⋱cot2h/22cot2h/22⋱cot2h/22⋱−cot1h/22−cot1h/22⋱⋮⋱0,and the 2nd differential matrix is(17)DN2=−π23h2−16⋱⋮⋱12csc21h212csc21h2⋱−12csc22h2⋱12csc22h212csc22h2⋱12csc21h2⋱12csc23h2⋮⋱−π23h2−16⋱⋮12csc23h2⋱12csc21h2⋱12csc22h212csc22h2⋱−12csc22h2⋱12csc21h212csc21h2⋱⋮⋱−π23h2−16.


### 3.2. Interpolation Collocation Method

To solve system (4), we consider a regular region Ω=0,2π×0,2π, and the interval 0,2π is divided into N different nodes. h=2π/N, xj=jh,yj=jh,andj=1,2,…,N. Using periodic sinc function (12), we first approximate ηx,y,t, ux,y,t, and vx,y,t as(18)ηx,y,t∼INηx,y,t=∑i=1N∑j=1NSNx−xiSNy−yiηxi,yj,t,ux,y,t∼INux,y,t=∑i=1N∑j=1NSNx−xiSNy−yiuxi,yj,t,vx,y,t∼INvx,y,t=∑i=1N∑j=1NSNx−xiSNy−yivxi,yj,t.At collocation nodesxp,yq, the following relations hold:(19)ηl,kxp,yq,t∼INηl,kxp,yq,t=∂l+kηxp,yq,t∂xl∂yk=∑i=1N∑j=1NSNlxp−xiSNkyq−yiηxi,yj,t,ul,kxp,yq,t∼INul,kxp,yq,t=∂l+kuxp,yq,t∂xl∂yk=∑i=1N∑j=1NSNlxp−xiSNkyq−yiuxi,yj,t,vl,kxp,yq,t∼INvl,kxp,yq,t=∂l+kvxp,yq,t∂xl∂yk=∑i=1N∑j=1NSNlxp−xiSNkyq−yivxi,yj,t.Noting(20)η=η11,η21,…,ηN1,η12,η22,…,ηN2,η1N…,ηNNT,u=u11,u21,…,uN1,u12,u22,…,uN2,u1N…,uNNT,v=v11,v21,…,vN1,v12,v22,…,vN2,v1N…,vNNT,where ηi,j=ηxi,yj,t, ui,j=uxi,yj,t,andvi,j=vxi,yj,t,i,j=1,2,…,N.Therefore, formula (19) can be written in the matrix form as follows:(21)ηl,k=DNl,kη,ul,k=DNl,ku,vl,k=DNl,kv,where DNl,k=DNl⊗DNk is the Kronecker product of matrix DNl and DNk and DN0=IN, IN is the N order unit matrix, respectively.Employing equations (18)–(21), the discrete form of equation (4) can be written as(22)∂∂tηuv=d1D+a11D0,0a12D0,0a13D0,0a21D0,0d2D+a22D0,0a23D0,0a31D0,0a32D0,0d3D+a33D0,0ηuv=f1η,u,vf2η,u,vf3η,u,v.Here,(23)η,u,v=η11,…,ηN1,η12,…,ηNN,u11,…,uN1,u12…,uN2,u1N…,uNN,v11,…,vN1,v12,…,vN2,v1N…,vNN,D=DN2,0+DN0,2,D0,0=IN⊗IN,f1η,u,v,f2η,u,v,f3η,u,v=f1η11,u11,v11,…,f1ηNN,uNN,vNN,f2η11,u11,v11,…,f2ηNN,uNN,vNN,f3η11,u11,v11,…,f3ηNN,uNN,vNN.Thus, we have the following system of ODEs:(24)∂Wt∂t=SWt+FW,t,in which(25)W=ηuv,S=d1D+a11D0,0a12D0,0a13D0,0a21D0,0d2D+a22D0,0a23D0,0a31D0,0a32D0,0d3D+a33D0,0,FW,t=f1η,u,vf2η,u,vf3η,u,v.The RK-4 for ODEs (23) is given as(26)k1=RWn,k2=RWn+τ2k1,k3=RWn+τ2k2,k4=RWn+τk3,Wn+1=Wn+τ6k1+2k2+2k3+k4,where RW=SWt+FW,t. An initial condition is required for starting the integration process and at each time level, boundary conditions should be applied as follows:(27)Wn+1=Bxl,tn+1,in which, B denotes the boundary operator.


### 3.3. Stability of the Method

The stability of system (23) can be discussed according to the eigenvalues of spatial operator S. Let the eigenvalues of S be λi; then, one of the following conditions should be fulfilled for the stability of system (23):(i)
If allλi are real and −2.78<τλi<0(ii)
If allλi are pure imaginary and −22<τλi<22(iii)
Ifλi are complex, then τλi should be in the stability boundary of the RK-4 methodIt is stated that in [29], there may be some tolerance that real parts of eigenvalues may be small positive numbers if the eigenvalues are complex.


## 3.1. Periodic Sinc Function

Sinc functions are used in different areas of physics and mathematics. A periodic sinc functionSN is defined as(8)SNx=sinπx/h2π/htanx/2,where h=2π/N and SN is the interpolation function of periodic δ function.It can be proved thatSNxx⟶0=1, and the first-order derivative of function SNx at xj=jh,j=1,2,…,N is(9)SN′xj=0,j%N=0,−1jcotjh/22,j%N≠0.The two-order derivative of functionSNx at xj=jh,j=1,2,…,N is(10)SN″xj=−π23h2−16,j%N=0,−−1jcsc2jh/22,j%N≠0.Givenh>0, we define the following interpolation space:(11)SpanSNx−jh,j=1,2,…,N.LetIN be the interpolation operator for any function ux defined on 0,2π. The interpolation function INux of sequence u1,u2,…,uN can be written as(12)ux∼INux=∑m=1NumSNx−xm,where uj=ujh,j=1,…,N. It is not difficult to derive simple expressions of the nth derivatives of INux at xj=jh:(13)INunxj=∑m=1NumSNnxj−xm,where xj−xm=j−mh. INunxj can be written in the matrix form as follows:(14)INunx1INunx2INunx3⋮INunxN=SNnxNSNnxN−1SNnxN−2…SNnx1SNnx1SNnxNSNnxN−1…SNnx2SNnx2SNnx1SNnxN…SNnx3……………SNnxN−1SNnxN−2SNnxN−3…SNnxNu1u2u3⋮uN.Noting thenth differential matrix DNn:(15)DNn=SNnxNSNnxN−1SNnxN−2…SNnx1SNnx1SNnxNSNnxN−1…SNnx2SNnx2SNnx1SNnxN…SNnx3……………SNnxN−1SNnxN−2SNnxN−3…SNnxN,where the 1st differential matrix is as follows:(16)DN1=0⋱⋮⋱−cot1h/22−cot1h/22⋱cot2h/22⋱cot2h/22cot2h/22⋱−cot1h/22⋱−cot3h/22⋮⋱0⋱⋮−cot3h/22⋱−cot1h/22⋱cot2h/22cot2h/22⋱cot2h/22⋱−cot1h/22−cot1h/22⋱⋮⋱0,and the 2nd differential matrix is(17)DN2=−π23h2−16⋱⋮⋱12csc21h212csc21h2⋱−12csc22h2⋱12csc22h212csc22h2⋱12csc21h2⋱12csc23h2⋮⋱−π23h2−16⋱⋮12csc23h2⋱12csc21h2⋱12csc22h212csc22h2⋱−12csc22h2⋱12csc21h212csc21h2⋱⋮⋱−π23h2−16.


## 3.2. Interpolation Collocation Method

To solve system (4), we consider a regular region Ω=0,2π×0,2π, and the interval 0,2π is divided into N different nodes. h=2π/N, xj=jh,yj=jh,andj=1,2,…,N. Using periodic sinc function (12), we first approximate ηx,y,t, ux,y,t, and vx,y,t as(18)ηx,y,t∼INηx,y,t=∑i=1N∑j=1NSNx−xiSNy−yiηxi,yj,t,ux,y,t∼INux,y,t=∑i=1N∑j=1NSNx−xiSNy−yiuxi,yj,t,vx,y,t∼INvx,y,t=∑i=1N∑j=1NSNx−xiSNy−yivxi,yj,t.At collocation nodesxp,yq, the following relations hold:(19)ηl,kxp,yq,t∼INηl,kxp,yq,t=∂l+kηxp,yq,t∂xl∂yk=∑i=1N∑j=1NSNlxp−xiSNkyq−yiηxi,yj,t,ul,kxp,yq,t∼INul,kxp,yq,t=∂l+kuxp,yq,t∂xl∂yk=∑i=1N∑j=1NSNlxp−xiSNkyq−yiuxi,yj,t,vl,kxp,yq,t∼INvl,kxp,yq,t=∂l+kvxp,yq,t∂xl∂yk=∑i=1N∑j=1NSNlxp−xiSNkyq−yivxi,yj,t.Noting(20)η=η11,η21,…,ηN1,η12,η22,…,ηN2,η1N…,ηNNT,u=u11,u21,…,uN1,u12,u22,…,uN2,u1N…,uNNT,v=v11,v21,…,vN1,v12,v22,…,vN2,v1N…,vNNT,where ηi,j=ηxi,yj,t, ui,j=uxi,yj,t,andvi,j=vxi,yj,t,i,j=1,2,…,N.Therefore, formula (19) can be written in the matrix form as follows:(21)ηl,k=DNl,kη,ul,k=DNl,ku,vl,k=DNl,kv,where DNl,k=DNl⊗DNk is the Kronecker product of matrix DNl and DNk and DN0=IN, IN is the N order unit matrix, respectively.Employing equations (18)–(21), the discrete form of equation (4) can be written as(22)∂∂tηuv=d1D+a11D0,0a12D0,0a13D0,0a21D0,0d2D+a22D0,0a23D0,0a31D0,0a32D0,0d3D+a33D0,0ηuv=f1η,u,vf2η,u,vf3η,u,v.Here,(23)η,u,v=η11,…,ηN1,η12,…,ηNN,u11,…,uN1,u12…,uN2,u1N…,uNN,v11,…,vN1,v12,…,vN2,v1N…,vNN,D=DN2,0+DN0,2,D0,0=IN⊗IN,f1η,u,v,f2η,u,v,f3η,u,v=f1η11,u11,v11,…,f1ηNN,uNN,vNN,f2η11,u11,v11,…,f2ηNN,uNN,vNN,f3η11,u11,v11,…,f3ηNN,uNN,vNN.Thus, we have the following system of ODEs:(24)∂Wt∂t=SWt+FW,t,in which(25)W=ηuv,S=d1D+a11D0,0a12D0,0a13D0,0a21D0,0d2D+a22D0,0a23D0,0a31D0,0a32D0,0d3D+a33D0,0,FW,t=f1η,u,vf2η,u,vf3η,u,v.The RK-4 for ODEs (23) is given as(26)k1=RWn,k2=RWn+τ2k1,k3=RWn+τ2k2,k4=RWn+τk3,Wn+1=Wn+τ6k1+2k2+2k3+k4,where RW=SWt+FW,t. An initial condition is required for starting the integration process and at each time level, boundary conditions should be applied as follows:(27)Wn+1=Bxl,tn+1,in which, B denotes the boundary operator.


## 3.3. Stability of the Method

The stability of system (23) can be discussed according to the eigenvalues of spatial operator S. Let the eigenvalues of S be λi; then, one of the following conditions should be fulfilled for the stability of system (23):(i)
If allλi are real and −2.78<τλi<0(ii)
If allλi are pure imaginary and −22<τλi<22(iii)
Ifλi are complex, then τλi should be in the stability boundary of the RK-4 methodIt is stated that in [29], there may be some tolerance that real parts of eigenvalues may be small positive numbers if the eigenvalues are complex.


## 4. Numerical Experiment

In this section, some numerical experiments are studied to demonstrate the accuracy of the present method.Experiment 1.
We consider system (3). For, bifurcation analysis of Experiment 1, see references [4–6]. Using the present method, taking N=100, the parameter d1=0.001, d2=0.001, d3=0.001, r=1, θ=5, ψ=5, ε=0.1, b=0.1, δ=0.1, c=1, a=1, k=0.4, and I=0.3. Numerical solution and pattern of Experiment 1 with the different initial conditions ηx,y,0 and initial condition vx,y,0=1 at the different t are shown in Figures 1–7.
The present method could be extended for use in solving other nonlinear cross-diffusion systems.Figure 1
Numerical solution and pattern of Experiment1 with initial condition ηx,y,0=0.6sechx/0.2+y×randN.
(a)
(b)
(c)Figure 2
Numerical solution and pattern of Experiment1 with initial condition ηx,y,0=sinsechx/2+y2+0.5.
(a)
(b)
(c)Figure 3
Numerical solution and pattern of Experiment1 with initial condition ηx,y,0=onesN+sechx2/0.02−9+2⋅y/0.01.
(a)
(b)
(c)Figure 4
Numerical solution and pattern of Experiment1 with initial condition ηx,y,0=sincosx/2+y2+0.5.
(a)
(b)
(c)Figure 5
Numerical solution and pattern of Experiment1 with initial condition ηx,y,0=−sinx2+y/10.
(a)
(b)
(c)Figure 6
Numerical solution and pattern of Experiment1 with initial condition ηx,y,0=0.8sincosy2−x2+0.2.
(a)
(b)
(c)Figure 7
Numerical solution and pattern of Experiment1 with initial condition ηx,y,0=0.8sinsiny2−x2+0.2.
(a)
(b)
(c)Experiment 2.
Let us consider the following nonlinear diffusion system:(28)∂u∂t=au1−ub−buvbu+v+d11∇2u+d12∇2v,∂v∂t=bubu+v−cv+d21∇2u+d22∇2v.


### 4.1. Bifurcation Analysis of Experiment2

In this section, we derive non-cross-diffusion-driven stability conditions.If0<c<1anda>b1−c, then the nondiffusive Experiment 2 has three constant equilibria: the zero equilibrium E0=0,0; the boundary equilibrium E1=b,0; and the positive equilibrium E∗=u∗,v∗ with u∗=ba+c−1b/a,v∗=b21−ca+c−1b/ac.Now, we do stability analysis of the equilibrium pointE∗. The linearized system of Experiment 2 at E∗ is(29)∂u∂t∂v∂t=a11a12a21a22uv+d11d12d21d22∇2u∇2v,with a11=−a−b1−c2,a12=−c2,a21=b1−c2,anda22=−c1−c. Jacobian matrix Ak is(30)Ak=a11a12a21a22−k2d11d12d21d22,where k>0 is often called the wave number, and we have the following eigenvalue of Ak:(31)λ1,2k=12trAk±trAk2−4detAk,where(32)trAk=a11+a22−k2d11+d22=−a−b1−c2−c1−c−k2d11+d22,(33)detAk=d11d22−d12d21k4−d22a11+d11a22−d21a12−d12a21k2+a11a22−a12a21=d11d22−d12d21k4−−a−b1−c2d22−c1−cd11+c2d21−b1−c2d12+a−b1−c2c1−c+bc21−c2.From equations (31)–(33), it is clear that if 0<c<1,a>b1−c holds, then a11a22−a12a21>0. It is easy to see that a11+a22<0 if and only if the condition a>b1−c2−c1−c holds. Therefore, without diffusion, i.e., dij=0,i,j=1,2, when(34)0<c<1,a>b1−c,a>b1−c2−c1−chold; the positive equilibrium E∗ of Experiment 2 is locally asymptotically stable.In this section, we derive cross-diffusion-driven instability conditions and show that these are a generalization of the classical diffusion-driven instability conditions in the absence of cross-diffusion.We assume thatd11d22−d12d21>0. This implies that the product of the self-diffusion coefficients d11 and d22 should diffuse much faster than the product of the cross-diffusion coefficients d12 and d21.From equations (31)–(33), obviously, if d21=0, a>b1−c2 holds. In this case, the positive equilibrium E∗ is stable and the Turing bifurcation cannot occur. Thus, the Turing bifurcation may occur only when d21>0. Choose d21 as the bifurcation parameter.From equations (31)–(33), Turing bifurcation at d21=d21∗, and the parameter d21∗ is determined by the following equality:(35)a>b1−c2,a−b1−c2d22+c1−cd11−c2d21∗+b1−c2d12<0,a−b1−c2d22+c1−cd11−c2d21∗+b1−c2d12+2d11d22−d12d21∗a−b1−c2c1−c+bc21−c2=0.Based on the above analysis, we obtain ifd21>0,d11d22−d12d21>0 and the condition a>b1−c2 holds, then the positive equilibrium E∗ of Experiment 2 is locally asymptotically stable for d21<d21∗ and unstable for d21>d21∗.Based on the above analysis, the sufficient condition of locally asymptotically stable of the equilibrium is obtained.Theorem 1.
Suppose thata,b,c>0,dij≥0,i,j=1,2, and the condition 0<c<1,a>b1−c holds.(i)
Ifdij=0,i,j=1,2, and a>b1−c2−c1−c, then the positive equilibrium E∗ of Experiment 2 is locally asymptotically stable(ii)
Ifd21=0 and a>b1−c2, then the positive equilibrium E∗ of Experiment 2 is locally asymptotically stable for any d11,d22>0 and d12≥0(iii)
Ifd21>0,d11d22−d12d21>0 and a>b1−c2, then the positive equilibrium E∗ of Experiment 2 is locally asymptotically stable for d21<d21∗ and unstable for d21>d21∗Bifurcation diagram of Experiment2 with the different parameters are shown in Figures 8 and 9. In Figures 8 and 9, hopf curve 1 is a−b1−c2d22+c1−cd11−c2d21∗+b1−c2d12=0, Turing curve 1 is a−b1−c2d22+c1−cd11−c2d21∗+b1−c2d12+2d11d22−d12d21∗a−b1−c2c1−c+bc21−c2=0, and hopf curve 2 is a=maxb1−c,b1−c2−c1−c,0<c<1.Figure 8
Bifurcation diagram ind21−a space of Experiment 2 with c=0.78,b=0.78,d11=1.2,d12=0.2,andd22=0.9⊳.Figure 9
Bifurcation diagram ind21−b space of Experiment 2 with c=0.83,a=0.69,d11=2.71,d12=2.5,andd22=5.


### 4.2. Numerical Simulation of Experiment2

Using the present method, takingN=100, the parameter d11=1, d12=0.1, d21=10.9686, d22=15, a=0.8, b=1.06, and c=0.6. Numerical solution and pattern of Experiment 2 with the different initial conditions ux,y,0 at the different t are shown in Figures 10–16.Figure 10
Numerical solution and Turing pattern of Experiment2 with initial condition ux,y,0=sechπ−x2+y2 and vx,y,0=1.
(a)
(b)
(c)
(d)Figure 11
Numerical solution and Turing pattern of Experiment2 with initial condition ux,y,0=cos−πx−0.42+y+0.42−e−20x+0.42+y−0.42 and vx,y,0=1.
(a)
(b)
(c)
(d)Figure 12
Numerical solution and Turing pattern of Experiment2 with initial condition ux,y,0=sinπx−0.42+y+0.42+e−20x+0.42+y−0.42 and vx,y,0=1.
(a)
(b)
(c)
(d)Figure 13
Numerical solution and Turing pattern of Experiment2 with initial condition ux,y,0=−sinπx2+y2 and vx,y,0=1.
(a)
(b)
(c)
(d)Figure 14
Numerical solution and Turing pattern of Experiment2 with initial condition ux,y,0=cosπx2+y2 and vx,y,0=1.
(a)
(b)
(c)
(d)Figure 15
Numerical solution and Turing pattern of Experiment2 with initial condition ux,y,0=12/25−sin10x2−510−sin10y−220 and vx,y,0=1.
(a)
(b)
(c)
(d)Figure 16
Numerical solution and Turing pattern of Experiment2 with initial condition ux,y,0=12/25−sin10x2−510+sincosy−22 and vx,y,0=1.
(a)
(b)
(c)
(d)


## 4.1. Bifurcation Analysis of Experiment2

In this section, we derive non-cross-diffusion-driven stability conditions.If0<c<1anda>b1−c, then the nondiffusive Experiment 2 has three constant equilibria: the zero equilibrium E0=0,0; the boundary equilibrium E1=b,0; and the positive equilibrium E∗=u∗,v∗ with u∗=ba+c−1b/a,v∗=b21−ca+c−1b/ac.Now, we do stability analysis of the equilibrium pointE∗. The linearized system of Experiment 2 at E∗ is(29)∂u∂t∂v∂t=a11a12a21a22uv+d11d12d21d22∇2u∇2v,with a11=−a−b1−c2,a12=−c2,a21=b1−c2,anda22=−c1−c. Jacobian matrix Ak is(30)Ak=a11a12a21a22−k2d11d12d21d22,where k>0 is often called the wave number, and we have the following eigenvalue of Ak:(31)λ1,2k=12trAk±trAk2−4detAk,where(32)trAk=a11+a22−k2d11+d22=−a−b1−c2−c1−c−k2d11+d22,(33)detAk=d11d22−d12d21k4−d22a11+d11a22−d21a12−d12a21k2+a11a22−a12a21=d11d22−d12d21k4−−a−b1−c2d22−c1−cd11+c2d21−b1−c2d12+a−b1−c2c1−c+bc21−c2.From equations (31)–(33), it is clear that if 0<c<1,a>b1−c holds, then a11a22−a12a21>0. It is easy to see that a11+a22<0 if and only if the condition a>b1−c2−c1−c holds. Therefore, without diffusion, i.e., dij=0,i,j=1,2, when(34)0<c<1,a>b1−c,a>b1−c2−c1−chold; the positive equilibrium E∗ of Experiment 2 is locally asymptotically stable.In this section, we derive cross-diffusion-driven instability conditions and show that these are a generalization of the classical diffusion-driven instability conditions in the absence of cross-diffusion.We assume thatd11d22−d12d21>0. This implies that the product of the self-diffusion coefficients d11 and d22 should diffuse much faster than the product of the cross-diffusion coefficients d12 and d21.From equations (31)–(33), obviously, if d21=0, a>b1−c2 holds. In this case, the positive equilibrium E∗ is stable and the Turing bifurcation cannot occur. Thus, the Turing bifurcation may occur only when d21>0. Choose d21 as the bifurcation parameter.From equations (31)–(33), Turing bifurcation at d21=d21∗, and the parameter d21∗ is determined by the following equality:(35)a>b1−c2,a−b1−c2d22+c1−cd11−c2d21∗+b1−c2d12<0,a−b1−c2d22+c1−cd11−c2d21∗+b1−c2d12+2d11d22−d12d21∗a−b1−c2c1−c+bc21−c2=0.Based on the above analysis, we obtain ifd21>0,d11d22−d12d21>0 and the condition a>b1−c2 holds, then the positive equilibrium E∗ of Experiment 2 is locally asymptotically stable for d21<d21∗ and unstable for d21>d21∗.Based on the above analysis, the sufficient condition of locally asymptotically stable of the equilibrium is obtained.Theorem 1.
Suppose thata,b,c>0,dij≥0,i,j=1,2, and the condition 0<c<1,a>b1−c holds.(i)
Ifdij=0,i,j=1,2, and a>b1−c2−c1−c, then the positive equilibrium E∗ of Experiment 2 is locally asymptotically stable(ii)
Ifd21=0 and a>b1−c2, then the positive equilibrium E∗ of Experiment 2 is locally asymptotically stable for any d11,d22>0 and d12≥0(iii)
Ifd21>0,d11d22−d12d21>0 and a>b1−c2, then the positive equilibrium E∗ of Experiment 2 is locally asymptotically stable for d21<d21∗ and unstable for d21>d21∗Bifurcation diagram of Experiment2 with the different parameters are shown in Figures 8 and 9. In Figures 8 and 9, hopf curve 1 is a−b1−c2d22+c1−cd11−c2d21∗+b1−c2d12=0, Turing curve 1 is a−b1−c2d22+c1−cd11−c2d21∗+b1−c2d12+2d11d22−d12d21∗a−b1−c2c1−c+bc21−c2=0, and hopf curve 2 is a=maxb1−c,b1−c2−c1−c,0<c<1.Figure 8
Bifurcation diagram ind21−a space of Experiment 2 with c=0.78,b=0.78,d11=1.2,d12=0.2,andd22=0.9⊳.Figure 9
Bifurcation diagram ind21−b space of Experiment 2 with c=0.83,a=0.69,d11=2.71,d12=2.5,andd22=5.


## 4.2. Numerical Simulation of Experiment2

Using the present method, takingN=100, the parameter d11=1, d12=0.1, d21=10.9686, d22=15, a=0.8, b=1.06, and c=0.6. Numerical solution and pattern of Experiment 2 with the different initial conditions ux,y,0 at the different t are shown in Figures 10–16.Figure 10
Numerical solution and Turing pattern of Experiment2 with initial condition ux,y,0=sechπ−x2+y2 and vx,y,0=1.
(a)
(b)
(c)
(d)Figure 11
Numerical solution and Turing pattern of Experiment2 with initial condition ux,y,0=cos−πx−0.42+y+0.42−e−20x+0.42+y−0.42 and vx,y,0=1.
(a)
(b)
(c)
(d)Figure 12
Numerical solution and Turing pattern of Experiment2 with initial condition ux,y,0=sinπx−0.42+y+0.42+e−20x+0.42+y−0.42 and vx,y,0=1.
(a)
(b)
(c)
(d)Figure 13
Numerical solution and Turing pattern of Experiment2 with initial condition ux,y,0=−sinπx2+y2 and vx,y,0=1.
(a)
(b)
(c)
(d)Figure 14
Numerical solution and Turing pattern of Experiment2 with initial condition ux,y,0=cosπx2+y2 and vx,y,0=1.
(a)
(b)
(c)
(d)Figure 15
Numerical solution and Turing pattern of Experiment2 with initial condition ux,y,0=12/25−sin10x2−510−sin10y−220 and vx,y,0=1.
(a)
(b)
(c)
(d)Figure 16
Numerical solution and Turing pattern of Experiment2 with initial condition ux,y,0=12/25−sin10x2−510+sincosy−22 and vx,y,0=1.
(a)
(b)
(c)
(d)


## 5. Conclusions and Remarks

In this paper, a sinc function interpolation collocation method is proposed for a class of nonlinear diffusion predator-prey systems. Some complex dynamics characters and pattern formations of predator-prey systems with the same parameters and different initial conditions are shown by using this new approach. Simulation results show the effectiveness of the new method. The present method could be extended for use in solving other nonlinear reaction-diffusion systems. In further work, we will be devoted to studying some fractional-order predator-prey systems.



---

*Source: 5404851-2020-03-23.xml*