A novel approach for numeric study of 2D biological population model

: In the present paper, the modified cubic B-spline differential quadrature method (MCB-DQM) has been implemented for the numerical computation of two-dimensional biological population model (BPM). The method is based on differential quadrature in which the weighting coefficients are computed by using MCB as a set of basis functions. We present three test problems to confirm the efficiency and accuracy of the method for BPM, which shows that the MCB-DQM solutions are in good agreement with the results obtained by the recent schemes: improved element-free Galerkin method by Zhang et al. and element-free kp-Ritz method by Cheng et al. The order of convergence of MCB-DQM for the solutions of BPM is shown to be quadratic.


Introduction
The dispersal or emigration plays a significant role in the regulation of population of some species. In the recent years, the biological problems have been received a much attention among the researchers due to its key role in the regulation of population of some species. The diffusion of a biological species in a region ℝ 2 is described by three functions: (i) population density ( , t), (ii) diffusion

PUBLIC INTEREST STATEMENT
The biological problems received much attention among the researchers due to its key role in the regulation of population of some species. Indeed, the computation of the exact solution each differential equation (DE) is too tough task. Therefore, the development of new schemes for numerical solution of PDEs with accurate, stable and efficient, is of vital importance. Modified cubic B-spline differential quadrature method (MCB-DQM) is implemented for solving two-dimensional biological population model (BPM). This is a new DQM with modified cubic B-splines as new set of basis functions. Three problems are considered to confirm the efficiency and accuracy of the method for BPM. MCB-DQM solutions are in good agreement with the recent results obtained using improved element-free Galerkin method by Zhang et al. and element-free kp-Ritz method by Cheng et al. The order of convergence of the proposed scheme is quadratic. velocity ( , t) and (iii) population supply ( , t) (Gurtin & Maccamy, 1977), where = (x, y) ∈ ℝ 2 and t is time.
The population density represents the number of individuals per unit volume at ( , t), and it's integral over any sub-region of region ℝ 2 gives the total population of at time t whereas ( , t) denotes the rate at which individuals are supplied in per unit volume at position at time t by births and deaths. The diffusion velocity ( , t) represents the average velocity of the individuals occupying the position at the time t, and it represents the flow of population from one point to another point.
The entities , and must be consistent with law of population balance (Gurney & Nisbet, 1975;Lu, 2000), given below in Equation (1): for every regular sub-region of ℝ 2 at time t where is the outward unit normal to the boundary of the region . The law (1) shows that "the sum of the rate of change of population of and the rate at which the individuals leave the region across its boundary is equal to the rate at which the individuals are supplied directly to the region ". For = ( ) and = − ( )∇ , where ( ) > 0 for > 0 and ∇ is the Laplace operator, the following nonlinear degenerate parabolic partial differential equation (DE) for is reduced to where ∇ 2 ≡ 2 x 2 + 2 y 2 and ( ) be function of (Gurtin & Maccamy, 1977). In Gurney and Nisbet (1975), ( ) is used as a special case in the modeling of the population of animals. The movements were made generally either by mature animals driven by mature invaders or by the young animals just reaching maturity moving out of their parental territory to establish breeding territory of their own. In each case, it is much more probable to consider that they will be directed toward nearby vacant territory. Therefore, in this model, the movement will take place almost exclusively down the population density gradient and will be more rapid at high population densities than at low ones. To model this situation, the authors further assumed a walk through a rectangular grid, in which at each step, an animal may either stay at its present location or may move toward the lowest population density (Gurney & Nisbet, 1975). The probability distribution between these two possibilities being determined by the magnitude of the population density gradient at the grid side is concerned. This model is correspond to Equation (2) with ( ) = 2 as Some properties of Equation (3) such as Holder-estimates and its solutions were studied in Lu (2000). Equation (3) is the generalization of the following lows: (a) For ( ) = c , c is a constant (Malthusian Law Gurtin and Maccamy, 1977). Gurtin & Maccamy, 1977).
In the recent years, a lot of numerical techniques have been developed for the numerical computation of time-dependent partial DEs (Arora, Mittal, & Singh, 2014;Arora & Singh, 2013;Cheng, Zhang, & Liew, 2014;Shivanian, 2013;Singh & Kumar, 2016c;Zhang et al., 2014). The existence, regularity and uniqueness of the weak solutions for degenerate parabolic equations (see, Aronson, 1986;Dibnedetto, 1993;Giuggioli & Kenkre, 2003;Gurney & Nisbet, 1975;Gurtin & Maccamy, 1977;Jager & Lu, 1997;Lu, 2000). The analytical solutions for time-fractional order population problems were obtained using fractional reduced transform method (Singh & Kumar, 2016a, 2016bSingh & Srivastava, 2015;Srivastava, Mishra, Kumar, Singh, & Awasthi, 2014) by . The computation of various type of population problems is done by using homotopy perturbation method (HPM) (Liu, Li, & Zhang, 2011;Roul, 2010), homotopy analysis method (HAM) (Arafa, Rida, & Mohamed, 2011), and variational iteration method (Shakeri & Dehghan, 2006). A meshless local radial point interpolation numerical method to simulate a nonlinear partial integro-DE arising in population dynamics is given by Shivanian (2013). Recently, the numerical computation of BPM model (4) is done by Cheng et al. (2014) using element-free kp-Ritz method, and by Zhang et al. (2014) using improved element-free Galerkin method (IEFGM).    The main aim of this paper is to implement "modified cubic B-spline differential quadrature method (MCB-DQM)" (Arora & Singh, 2013;Singh & Bianca, 2016) for the numerical computation of 2D BPM. The scheme is based on differential quadrature where the modified cubic-B-spline functions are used as set of basis functions to determine the weighting coefficients. The MCB-DQM is used to convert the given system of PDEs into a system of first order ODEs, in time. The resulting system of ODEs is solved using the SSP-RK54 scheme. Three test problems are considered to demonstrate the accuracy and utility of the method. The maximum absolute errors L ∞ and L 2 errors in the MCB-DQM solutions have been compared with the errors due to some recent exiting schemes.
The article is organized as follows. In Section 2, the description of the MCB-DQM is given. In Sections 3, procedure for implementation of method is described. Three numerical examples are given to establish the applicability and accuracy of the proposed method in Section . Section 5 concludes our study.  Bellman, Kashef, and Casti (1972) have introduced the DQM. DQM is an approximation to the derivatives of a function that is obtained by means of the weighted sum of the functional values at certain discrete points. In DQM, the weighting coefficients are evaluated using various test functions such as spline functions, sinc function, Lagrange interpolation polynomials, modified (extended) cubic B-spline, etc. (Arora & Singh, 2013;Korkmaz & Dag, 2011, 2013Quan & Chang, 1989a, 1989bShu, Chen, Xue, & Du, 2001;Shu & Richards, 1992;Singh, Arora, & Singh, 2016;Singh & Kumar, 2016c;Zhong, 2004). The weighting coefficients are dependent on the spatial grid spacing only, so one can assume NM grid

Description of MCB-DQM
The approximate value of nth-order partial derivative of the function (x, y, t) with respect to x at (x i , y j ) is given by and the approximate value of mth-order partial derivative of the function (x, y, t) with respect to y at (x i , y j ) is given by ij denote the weighting coefficients of the partial derivatives n ij m ij

Computation of the weighting coefficients a (r)
ij and b (r) ij (r = 1, 2) To find the weighting coefficients a (1) ij , the modified cubic B-spline k (x), k = 1, 2, … , N is used in Equation (6) due to fixed y axis. The first-order derivative approximation at the grid point (x i , y j ) is By Lemma 1 and Equations (9) and (10) is reduced into a tridiagonal system of linear equations where Φ = [ ij ] is the coefficient matrix of order N whose ith row is given by , and the coefficient vector �⃗ Now, we apply the well-known "Thomas algorithm" to solve the resulting tridiagonal system of equations which provides the vector � ⃗ a [i], that is, the weighting coefficients a (1) Similarly, the weighting coefficients b (1) ij can be computed by means of the modified cubic B-spline k (y), k = 1, 2, … , M in Equation (7).
Using the coefficients a (1) ij and b (1) ij , the weighting coefficients of higher order spatial derivatives can be obtained using the Shu's (2000) recursive formulae where r represents the order of the spatial derivative, a (r) ij and b (r) ij are the weighting coefficients of rth-order partial derivatives of u(x, y, t) at the point (x i , y j ) with respect to x and y, respectively. In particular, the weighting coefficients a (2) ij , b (2) ij of order 2 can be obtained by taking r = 2 in (12).

Implementation of MCB-DQM for 2D BPM
After putting the values of the spatial derivatives approximated using MCB-DQM, Equation (4) can be rewritten as On implementing the conditions (5), Equation (13) reduced to where L is nonlinear differential operator defined by The resulting system of initial values problem (14) can be solved by various schemes available, in literature. Strong stability preserving time discretization methods were introduced to address the need for nonlinear stability properties in the time discretization, as well as the spatial discretization, of hyperbolic Partial DEs. Strong stability preserving Runge-Kutta (SSP-RK) schemes have some excellent properties such as large regions of absolute stability and low storage, SSP-RK54 is strongly stable scheme with Euler time discretizations for nonlinear hyperbolic PDEs (see Gottlieb, 2005;Spiteri & Ruuth, 2002, 2004 and the papers therein). This is why, we prefer SSP-RK54 scheme as defined below, to solve the system of first-order ODEs (14): Lemma 2 (Prenter, 1975) If the function u ∈ C 4 [a, b] such that where j (x) is the cubic B-spline function. Then for any partition of [a, b] with the grids distributed uniformly , and if ∈ C 4 ( , ), then from Lemma 2, we get (1) = m + 0.391752226571890 △ tL ( Therefore, we get This shows that the order of convergence of the MCB-DQM for BPM is two, which is confirmed numerically from Tables 4 and 8.

Numerical studies
This section deals with the main goal, the numerical study of three test problems of BPM approximated by MCB-DQM with SSP-RK54 scheme. The accuracy and the efficiency of the method is measured by evaluating the computational order, the L 2 -error norm and the maximum error norm (L ∞ ), and their computational time.

The errors in Example 3 for large time levels 2 ≤ t ≤ 10 with h x = h y = 1∕12
Errors t = 2 t = 5 t = 8 t = 10  The solutions are obtained with △t = 0.0001. In Table 1, MCB-DQM solutions at time 0.1 ≤ t ≤ 0.5 are compared with the solutions in Zhang et al. (2014) and the exact solutions. The L 2 , L ∞ errors and CPU time at t = 0.1, 0.4 is reported in Table 2 while the L 2 and L ∞ error norms for large time 5 ≤ t ≤ 20 is reported in Table 3. Table 4 shows that the convergence order for the BPM is quadratic. The computational time in seconds, CPU(s), for different time intervals t ≤ 20 and various grid sizes with △t = 0.0001 is also reported in the above Tables 1-4. The absolute errors in the BPM at different time levels t ≤ 1 is depicted in Figure 1, and the surface solution is depicted in Figure Table 5, and L 2 and L∞ error norms for △t = 0.0001 and different values of h x = h y are reported in Table 6, while the L 2 and L ∞ error norms for large time levels 5 ≤ t ≤ 20 are reported in Table 7. The computational time in seconds, CPU(s), for different time intervals t ≤ 20 and various grid sizes with △t = 0.0001 is also reported in the above Tables 5-8. The contour and surface plots of absolute errors in the BPM at t = 0.1, 0.5, 1.0 are depicted in Figure 3, and the surface solutions at these time intervals is depicted in Figure 4. It is evident that our results are in good agreement with the results obtained in Zhang et al. (2014) and exact solutions. Table 8 shows that the convergence quadratic. The behavior of ( , t) shows the similar characteristics as depicted in Zhang et al. (2014).
Example 3 We consider BPM (4) with the parameters h = 1∕5, = 1 and r = 0 in the region = [0, 1] 2 , subject to the initial condition:  The solutions are computed with △t = 0.0001. In Table 9, the computed results, for h x = h y = 1∕13, are compared with the recent results in Zhang et al. (2014) and the exact solutions. The L 2 and L ∞ error norms for large time levels 2 ≤ t ≤ 10 are reported in Table 10  seconds, CPU(s), for different time intervals t ≤ 10, and h x = h y = 1∕12 with △t = 0.0001 is also reported in the Tables 9 and 10. The contour plots of absolute errors in the BPM at t = 0.1, 0.2, 0.3, 0.5 are depicted in Figure 5, and the surface solutions at these time intervals is depicted in Figure 6. It is evident that our results are in good agreement with the results obtained in Zhang et al. (2014) and exact solutions. The behavior of ( , t) shows the similar characteristics as depicted in Cheng et al. (2014), Zhang et al. (2014).

Conclusions
In the present paper, the MCB-DQM has been implemented for the numerical computation of 2D BPM. We present three test problems to confirm the efficiency and accuracy of the method for BPM, which shows that the the MCB-DQM solutions are in good agreement with that of obtained by the recent results using IEFGM by Zhang et al. (2014) and element free kp-Ritz method by Cheng et al. (2014).
Further, the scheme is tested for large time for the given problems which shows that the results obtained for large time levels are in good agreement with the exact solutions. It is found that the order of convergence of MCB-DQM for the solutions of BPM is quadratic.