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

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.


Introduction
Mathematical models have played a significant role in understanding the dynamics of populations for over a century. ese 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 u(t) be the population size of a prey and predator species, respectively, and suppose that these functions obey the Gause-type model: In [2], the authors give a predator-prey population dynamics.
is model includes all three components (predator, prey, and subsidy) in a single spatial location: 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.
is leads to the following set of predator-prey systems with diffusion: where η � η(x, y, t) and v � v(x, y, t) are the prey and predator populations, u � u(x, y, t) is the quantity of subsidy, t is time, r, k,θ, a, I, ψ, c, ε, η > 0 are all nonnegative kinetic parameters, and d 1 , d 2 , and d 3 are the positive diffusion coefficients for the three populations; ∇ � (z 2 /zx 2 ) + (z 2 /zy 2 ). e sufficient condition of locally asymptotically stability of the equilibrium has been obtained in [4][5][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 reactiondiffusion system [7][8][9], which is described as zη zt � d 1 Δη + a 11 η + a 12 u + a 13  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][19][20], the perturbation method and variational iteration method [21,22], barycentric interpolation collocation method [23][24][25][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.

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 , u 0 , v 0 ), so Now, we do linear stability analysis of the equilibrium point E * . We evaluate the Jacobian matrix A 0 of the system at E * as In general, if real part of eigenvalues of A 0 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. e Jacobian matrix A k of system (4) is 2 Complexity Turing bifurcation occurs when the equilibrium state that is stable in absence of nondiffusion becomes unstable in presence of cross-diffusion. us, the only way for the E * to become an unstable point of cross-diffusion system (4) is that the real part of eigenvalues A k is positive; then, diffusion system (4) is unstable.

Description of the Method
3.1. Periodic Sinc Function. Sinc functions are used in different areas of physics and mathematics. A periodic sinc function S N is defined as where h � (2π/N) and S N is the interpolation function of periodic δ function. It can be proved that S N (x)| x⟶0 � 1, and the first-order derivative of function S N (x) at x j � jh, j � 1, 2, . . . , N is Given h > 0, we define the following interpolation space: Let I N be the interpolation operator for any function u(x) defined on [0, 2π]. e interpolation function I N u(x) of sequence u 1 , u 2 , . . . , u N can be written as where u j � u(jh), j � 1, . . . , N. It is not difficult to derive simple expressions of the nth derivatives of I N u(x) at x j � jh: where can be written in the matrix form as follows: Noting the nth differential matrix D (n) N : Complexity 3 where the 1st differential matrix is as follows: and the 2nd differential matrix is h � (2π/N), x j � jh, y j � jh, and j � 1, 2, . . . , N. Using periodic sinc function (12), we first approximate η(x, y, t), u(x, y, t), and v(x, y, t) as At collocation nodes (x p , y q ), the following relations hold: where η i,j � η(x i , y j , t), u i,j � u(x i , y j , t), and v i,j � v(x i , y j , t), i, j � 1, 2, . . . , N. erefore, formula (19) can be written in the matrix form as follows: where D (l,k) N � D (l) N ⊗ D (k) N is the Kronecker product of matrix D (l) N and D (k) N and D (0) N � I N , I N is the N order unit matrix, respectively.
Employing equations (18)- (21), the discrete form of equation (4) can be written as Here, Complexity 5 us, we have the following system of ODEs: in which e RK-4 for ODEs (23) is given as where R(W) � SW(t) + F(W, t). An initial condition is required for starting the integration process and at each time level, boundary conditions should be applied as follows: in which, B denotes the boundary operator.

Stability of the
Method. e 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 − 2

then τλ i should be in the stability boundary of the RK-4 method
It 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.

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][5][6]. Using the present method, taking N � 100, the parameter Numerical solution and pattern of Experiment 1 with the different initial conditions η(x, y, 0) and initial condition v(x, y, 0) � 1 at the different t are shown in Figures 1-7. e present method could be extended for use in solving other nonlinear cross-diffusion systems.

Experiment 2.
Let us consider the following nonlinear diffusion system:

Bifurcation Analysis of Experiment 2.
In this section, we derive non-cross-diffusion-driven stability conditions. If 0 < c < 1 and a > b(1 − c), then the nondiffusive Experiment 2 has three constant equilibria: the zero equilibrium E 0 � (0, 0); the boundary equilibrium E 1 � (b, 0); and the positive equilibrium Now, we do stability analysis of the equilibrium point E * . e linearized system of Experiment 2 at E * is zu zt zv zt where k > 0 is often called the wave number, and we have the following eigenvalue of A k :               Complexity      Complexity Figure 13: Numerical solution and Turing pattern of Experiment 2 with initial condition u(x, y, 0) � − sin(π(x 2 + y 2 )) and v(x, y, 0) � 1.  Complexity   Figure 14: Numerical solution and Turing pattern of Experiment 2 with initial condition u(x, y, 0) � cos(π(x 2 + y 2 )) and v(x, y, 0) � 1.      hold; 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 that d 11 d 22 − d 12 d 21 > 0. is implies that the product of the self-diffusion coefficients d 11 and d 22 should diffuse much faster than the product of the cross-diffusion coefficients d 12 and d 21 .
From equations (31)-(33), obviously, if d 21 � 0, a > b(1 − c 2 ) holds. In this case, the positive equilibrium E * is stable and the Turing bifurcation cannot occur. us, the Turing bifurcation may occur only when d 21 > 0. Choose d 21 as the bifurcation parameter.
From equations (31)-(33), Turing bifurcation at d 21 � d * 21 , and the parameter d * 21 is determined by the following equality: Based on the above analysis, we obtain if d 21 > 0, d 11 d 22 − d 12 d 21 > 0 and the condition a > b(1 − c 2 ) holds, then the positive equilibrium E * of Experiment 2 is locally asymptotically stable for d 21