Numerical simulation for a class of predator–prey system with homogeneous Neumann boundary condition based on a sinc function interpolation method

For the nonlinear predator–prey system (PPS), although a variety of numerical methods have been proposed, such as the difference method, the finite element method, and so on, but the efficient numerical method has always been the direction that scholars strive to pursue. Based on this question, a sinc function interpolation method is proposed for a class of PPS. Numerical simulations of a class of PPS with complex dynamical behaviors are performed. Time series plots and phase diagrams of a class of PPS without self-diffusion are shown. The pattern is obtained by setting up different initial conditions and the parameters in the system according to Turing bifurcation condition. The numerical simulation results have a good agreement with theoretical results. Simulation results show the effectiveness of the method.


Introduction
The PPS is a basic ecological system that exists widely in nature and is an essential component of ecosystems such as oceans, lakes, wetlands, forests, and grasslands. The predatory process plays an important role in promoting life evolution, maintaining ecological balance, and maintaining biodiversity. Therefore, research on PPS is crucial to the exploration of the fundamental nature of ecosystems. In [1,2], a PPS with Beddington-DeAngelis-type functional response is proposed and analyzed. In [3], a PPS with general Holling type functional response is given. In [4], a modified Leslie-Gower-type PPS with Holling's type II functional response is studied. In [5], Paul and Ghosh gave preypredator-generalist predator system of the following form: © The Author(s) 2020. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
with the initial conditions where x, y, and z are respectively the prey biomass, predator biomass, and top predator biomass at any time t, r and k are respectively the intrinsic growth rate and the environmental carrying capacity of prey species. α 0 , β 0 , γ 0 , α 1 , β 1 , γ 1 , m 1 , m 2 are all parameters. When studying the spatial distribution structure and maintaining biodiversity [6,7] of predators and prey populations, the reaction diffusion system [8][9][10][11][12][13] can more accurately describe the interaction between predators and preys. In the PPS, spatial diffusion is reflected in the predator's efforts to catch up with the prey, and the prey's efforts to escape the predator's pursuit. In [5], if diffusion behavior of predator-prey biomass is considered, the following system can be obtained [14]: where η = η(x, y, t), u = u(x, y, t), and v = v(x, y, t) are respectively the prey biomass, predator biomass, and top predator biomass at any time t, r and k are respectively the intrinsic growth rate and environmental carrying capacity of the prey species. α 0 , β 0 are respectively the predation rate of the predator and top predator on prey species. α 1 and β 1 are respectively measure of the conversion rate of prey species to its predator species and γ 1 is the conversion rate of predator species to the top predator species. γ 0 is the predation rate of top predator on predator species. m 1 and m 2 are respectively the natural death rate of the predator and top predator; and d 1 , d 2 , and d 3 are positive diffusion coefficients, , the smooth boundary is ∂Ω, homogeneous Neumann boundary condition, namely ∂u ∂n | ∂Ω = ∂v ∂n | ∂Ω = 0. For the nonlinear PPS, a variety of numerical methods have been proposed, such as finite difference method [15], B-spline method [16], finite element method [17,18], spectral method [19][20][21], the perturbation method and variational iteration method (VIM) [22,23], barycentric interpolation collocation method (BICM) [24][25][26][27], reproducing kernel method (RKM) [28][29][30][31][32], etc. Nevertheless, the efficient numerical method has always been the direction that scholars strived to pursue. Based on this question, a sinc function interpolation method is proposed for a class of PPS. PPS (2) has at most five equilibrium points as follows:
The Jacobian matrix of nondiffusive system (2) at arbitrary equilibrium point (η, u, v) is given as follows: The Jacobian matrix A λ of system (2) is Turing bifurcation occurs when the equilibrium state is stable in absence of nondiffusion, but it becomes unstable in presence of cross-diffusion. Thus, if there exists λ, the the equilibrium state becomes an unstable point of the cross-diffusion system (2), and if the real part of eigenvalues A k is positive, then diffusion system (2) is unstable.

Description of the sinc function interpolation method
To solve system (2), we consider a regular region N , x j = jh, y j = jh, j = 1, 2 . . . , N . Sinc functions are used in different areas of physics and mathematics. A periodic sinc function where h = 2π N , S N is the interpolation function of periodic δ function. It can be proved that S N (x)| x→0 = 1. S N (x ix j ) is an N order unit matrix, respectively.
Let I N be the interpolation operator such that for functions η(x, y, t), u(x, y, t), and v(x, y, t) defined on [0, 2π] with homogeneous Neumann boundary condition, the interpolation functions I N η(x, y, t), I N u(x, y, t), and I N v(x, y, t) y j , t) can be written as follows [33]: At collocation nodes (x p , y q ), the following relations hold: Therefore, formula (8) can be written as the following matrix form: 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,0) = I N ⊗ I N , D (0) N = I N , I N is an N order unit matrix, respectively. Employing Eqs. (7), (9), and (10), the discrete form of Eq. (2) can be written as follows: Using ode45 in MATLAB to solve Eq. (11) with different initial conditions, we can get the numerical solution of system (2).

Numerical experiments
In this section, we give some numerical illustrations for better explanation of the above analytical results using different initial conditions and parameters.
Taking the parameters β 0 = β 1 = 0, c 1 = 1, c 2 = 1, c 3 = 1, and using the present method, time series plots for Experiment 1 with different parameters are given in Fig. 2 Table 2 Figure 2 Time series plots for Experiment 1 with different parameters, for parameters see Table 2 Figure 3 Phase diagram of Experiment 1 with the value of the parameter, for parameters see Table 2 Figure 4 Phase diagram of Experiment 1 with the value of the parameter, for parameters see Table 2 diagram of Experiment 1 with the value of the parameter is shown in Figs. 3-4. Figure 2 shows that the coexistence equilibrium P 4 exists.   Table 3 Figure 6 Numerical solution and pattern of Experiment 2, for initial conditions see Table 3 Figure 7 Numerical solution and pattern of Experiment 2, for initial conditions see Table 3 Figure 8 Numerical solution and pattern of Experiment 3, for initial conditions see Table 3 Figure 9 Numerical solution and pattern of Experiment 3, for initial conditions see Table 3 Figure 10 Numerical solution and pattern of Experiment 3, for initial conditions see Table 3 5 Conclusions