Modified Chebyshev collocation method for delayed predator–prey system

In this study, the approximate solutions of the predator–prey system with delay have been obtained by using the modified Chebyshev collocation method. The main technique is that this method transforms the original problem into a system of nonlinear algebraic equations. By using the residual function of the operator equations, error differential equations are constructed and thus the approximate solutions are corrected. A numerical example is given to confirm the reliability and applicability of the method, and comparisons with existing results are given. The numerical results show that the obtained solutions are in good agreement with earlier studies.


Introduction
The predator-prey model is the essential model in studying the population dynamics of many species, which was initially proposed by Lotka and Volterra (see [1,2]). It has wide application in various research areas, such as chemical processes (see [3,4]), bioparticles granulation (see [5]), the interaction of microorganisms and ecosystems (see [6,7]). In recent years, many researchers have worked on the Lotka-Volterra type predator-prey system. In particular, Zhu et al. [8] focused on competitive Lotka-Volterra model in random environments. Li et al. [7] studied the canard phenomenon for predator-prey systems with response functions of Holling types. Badri et al. [9] dealt with the stabilization of the feasible equilibrium point of a special class of nonlinear quadratic systems known as Lotka-Volterra systems.
In order to get a more realistic model, the time delay has been taken into account in the predator-prey system. The analytical and dynamical aspects of such time delay models have been studied extensively by many researchers (see [10][11][12]). Also, there have been many studies interested in obtaining numerical solutions to the predator-prey system. For example, Capobianco [13] solved the numerical solution of Lotka-Volterra by using high performance parallel numerical methods. Susmita Paul [14] explained how to solve the Lotka-Volterra predator-prey model by using the Runge-Kutta-Fehlberg (RKF) method. Gokmen [15] used Taylor's collocation method to find the numerical solution of the predator-prey system with delay.
Chebyshev polynomials have become very important in numerical analysis. They are widely used because of their advantages, such as the roots of the first kind of Chebyshev polynomials (Gauss-Lobatto nodes) being used in polynomial interpolation for minimizing the Runge phenomena, providing the best uniform approximation of polynomials in continuous functions (see [35][36][37]). Most commonly used techniques with Chebyshev polynomials have been examined in [38][39][40] and the references therein.
Motivated by the above discussion, we are mainly interested in applying a modified Chebyshev collocation method for the time-delay predator-prey model in [10] as follows: with initial conditions y 1 (0) = α, where y 1 (t) and y 2 (t) are interpreted as the densities of prey and predator respectively, r 1 > 0 is the growth rate of prey in the absence of predators, a 11 > 0 denotes the selfregulation constant of prey, a 12 > 0 describes the predation of prey by predators, r 2 > 0 is the death rate of predators in the absence of prey, a 21 > 0 is the conversion rate for predators, a 22 > 0 describes the intraspecific competition among predators and τ is the generation time of the prey species, α, β are constant. The objective of this paper is to obtain the approximation solutions of system (1) in the form of truncated Chebyshev series. The primary benefit of this method is the nonlinear term that can easily be dealt with without any extra efforts. Other advantages include this method being nondifferential, nonintegral, and easily implemented on a computer.
This paper is organized as follows: In Sect. 2, a brief review of the shifted Chebyshev polynomial and its properties is provided. In Sect. 3, we apply the collocation method for system (1) using the shifted Chebyshev polynomial. In Sect. 4, we construct a fundamental matrix equation for system (1). In Sect. 5, we introduce the technique of residual error correction in order to check the accuracy of the method. Finally, a numerical example is presented to verify the efficiency and accuracy of the proposed method.

Shifted Chebyshev polynomials and their properties
In this section, we introduce Chebyshev polynomials. The Chebyshev polynomials are the sets of orthogonal polynomials and they are simply related to the trigonometric functions (see [41,42]) by the formula T n (cos θ ) = cos(nθ ) with θ ∈ [0, π]. The Chebyshev polynomial T n (x) of the first kind is a polynomial in x of degree n, defined by the following relation [43]: Since we use polynomial on t ∈ [0, L] for any real number L > 0, we can obtain the shifted Chebyshev polynomials T * n (t) = T n ( 2t L -1) by introducing the change of variable x = 2t/L -1, t ∈ [0, L]. The shifted Chebyshev polynomial T * n (t) satisfies the recurrence relation as follows: with the boundary condition T * n (0) = (-1) n , T * n (L) = 1.
And T * n (t) satisfies the discrete orthogonality condition where the interpolation points t k are chosen to be the Chebyshev-Gauss-Lobatto associated with the interval [0, L] and t k = L 2 (1 -cos(k π N )), k = 0, 1, 2, . . . , N . The summation symbol with double primes denotes a sum with both the first and last term halved [43].

Method of solution
Continuous and bounded functions y s (t) (s = 1, 2) can be approximated in terms of shifted Chebyshev polynomials in the interval [0, L] as follows: Using the discrete orthogonality relation (3), coefficient c sk in (4) is given by Our aim is to obtain the unknown coefficients y s (t i ) for i = 0, 1, 2, . . . , N , and the method of solution being considered should be programmable in a computer. From equations (4) and (5) we can obtain the function y sN (t) as follows: where We know that where m1, m2, and m3 are respectively N , 0, 2N for odd N and 0, 2N , 0 for even N . Then from the above equation we can write y sN (t) as follows:

Fundamental matrix equation for system (1)
To obtain the fundamental matrix equations of system (1), we substitute equations (6) and (7) into system (1). We get the fundamental matrix system Let We write equations (8) as follows: Then we can rewrite equations (9) as follows: where and By substituting the interpolation points t k into equations (10), we have two nonlinear and Let and In addition, by the initial value, we have Thus, replacing first rows of the argument matrix W 1 , W 2 by T(t 0 ), we have and Then we can rewrite equations (11) and (12) as follows: where

Error estimation and residual correction
This section is devoted to checking the accuracy of our method. Since the exact solution of system (1) cannot be obtained, we will use the residual correction to obtain better approximate solutions. Residual correction is a process when the obtained approximate solution is substituted into the original equation, and a system whose solution is the error corresponding to the approximate solution is obtained. In the sequence, substituting the approximate solution y sN (t) (s = 1, 2) into system (1), we obtain where E sN (t) (s = 1, 2) denotes the residual functions. We define the error corresponding to y 1N (t) and y 2N (t) as follows: and e 2N (t) = y 2 (t)y 2N (t).
Substituting y 1 (t) and y 2 (t) into system (1), we have Similar to system (1), this system is also nonlinear delay differential system with initial values e 1N (0) = 0 and e 2N (0) = 0. The unknown functions are e 1N (t) and e 2N (t). We will apply the method of Sect. 3 to Eq. (19) in order to obtain the approximate solutions. Let e 1N,M (t) and e 2N,M (t) be the estimation solutions of errors e 1N (t) and e 2N (t). We can obtain the new approximate solutions as follows: Then y 1N,M (t) and y 2N,M (t) are the correction solutions which are more accurate than y 1N (t) and y 2N (t). We will use the residual functions to measure the accuracy of numerical solutions by using y 1N,M (t) and y 2N,M (t) instead of y 1N (t) and y 2N (t). In the next section we use an example to demonstrate the above idea.

Numerical application
In this section, we demonstrate our method by a detailed example. We give the values of approximate solutions y sN (t), (s = 1, 2) at selected points of the given interval for different N values.
Example 1 ( [44]) We consider the following system: with α = 1 and β = 0.2. In order to obtain y 1N (t) and y 2N (t) with N = 5, 6, and 7, we apply the method of Sect. 3 for Eq. (20). Then we have y 15 (t) = 0.9141 -9.0596 × 10 -3 T * 1 (t) + 9.1229 × 10 -2 T * 2 (t)  Approximate solutions for the prey population with N = 5, 6, and 7 by the present method and the Taylor collocation method [44] The approximate solutions of prey and predator and comparison with the results of Ref. [44] are presented in Fig. 1 and Fig. 2. Figures show that the proposed method preserves the positivity of the solutions, which is the part of the solutions of Eq. (20). To examine their accuracy, we considered the absolute residual errors of these approximate solutions. Figure 3 plots the absolute residual errors of Example 1. In Table 1, we list the absolution residual errors of the present method and Ref. [44]. It is seen from the table and figures that the absolute residual error values are decreasing as we increase the parameter N , which are in good agreement with the results given in Ref. [44].
For implementation residual error correction in Sect. 5, we apply again the method of Sect. 3 for Eq. (20) with choosing N = 4 and M = 5, 6. The approximate solutions of y 14 (t) and y 24 (t) are found as follows: y 14 (t) = 0.91279 -1.1906 × 10 -2 T * 1 (t) + 6.5714 × 10 -2 T * 2 (t) -1.4116 × 10 -2 T * 3 (t) -4.5221 × 10 -3 T * 4 (t), Figure 2 Approximate solutions for the predator population with N = 5, 6, and 7 by the present method and the Taylor collocation method [44] y 24 (t) = 0.19698 -6.3365 × 10 -3 T * 1 (t) + 7.2034 × 10 -3 T * 2 (t) To realize the error approximate concept in Sect. 5 with N = 4 and M = 5, 6, the estimated errors are obtained, namely Then we can obtain our improved approximate solutions: and y 14,6 (t) = y 14 (t) + e 14,6 (t), y 24,6 (t) = y 24 (t) + e 24,6 (t).    The improvement values of y 14 (t) and y 24 (t) with M = 5, 6 are given in Fig. 4. It revealed that the proposed technique preserved the positive solutions of the given delayed prey-predator system. In order to understand how much improvement is provided by this scheme, the absolute residual errors of the original approximate solutions y 14 (t) and y 24 (t) are shown together with those of corrected solutions in Fig. 5. In Table 2, we list the residual errors of the improvement solutions at some point on our interval. In view of Fig. 5 and Table 2, absolute residual errors of y 14,5 (t) and y 24,5 (t) are smaller than those of y 14 (t) and y 24 (t) respectively, and absolute residual errors of y 14,6 (t) and y 24,6 (t) are smaller. Hence, we can comment that residual correction in general provides a certain improvement in the approximate solutions for Eq. (20).

Conclusion
In this paper, a modified Chebyshev collocation method based on the residual correction technique is presented to solve the Lotka-Volterra model with delay. An efficient error estimation can be made by using this technique. The key advantages of this approach are its low-cost computing and simplicity of implementation. Also the present method has the ability to convert the given problem into a system of mathematical equations, which can be solved easily using MATLAB or MAPLE software. Our numerical results are compared with the numerical results of [44]. The results show that they are in good correspondence with the results obtained in [44]. Based on the above facts, the modified Chebyshev collocation method is a powerful mathematical tool to obtain the numerical solutions of a nonlinear system. In future the proposed method will be applied to the fractional Lotka-Volterra biological model with and without a time delay. It is hoped that the biological relevance of the numerical results, such as stability and chaotic behavior, can be obtained. Similarly, numerical techniques may be designed for fractional reaction diffusion systems.