Numerical simulation based on meshless technique to study the biological population model

A kind of spectral meshless radial point interpolation method is proposed to degenerate parabolic equations arising from the spatial diffusion of biological populations and satisfactory agreements is archived. This method is based on collocation methods with mesh-free techniques as a background. In contrast to the ﬁnite-ele-ment method and those meshless methods based on Galerkin weak form, such as element-free Galerkin, there is no integration tools in this approach. Furthermore, some numerical experiments are given to validate the accuracy of the results and stability of the present method.


Introduction
The biological population problems have attracted much attention and research recently [1,2]. Biologists believe that dispersal or emigration play a key role in the regulation of population of some species. A persuasive example of this suggestion occurs in a paper by Carl [3], whose observations on a population of arctic ground squirrels showed that this species achieves population control by the dispersal of animals from densely populated areas of favorable habitat into unfavorable areas, where burrow sites are not available, and where they are subjected to intensive predation.
Consider the following nonlinear biological population model with a given initial condition uðx; 0Þ, where u denotes the population density, and r represents the population supply due to births and deaths. We can consider a more general form for r as hu a ð1 À ru b Þ, in which h; a; r; b 2 R. The field uðx; tÞ gives the number of individuals' per-unit volume at position x and time t, and it is integral over any subregion X to give the total population of X at time t. The field rðuÞ gives the rate at which individuals are supplied (per-unit volume) exactly at x through the births and deaths. The flow of population from point to point is then depicted by the diffusion velocity vðx; tÞ, which provides the average velocity of those individuals who occupies x at the time t. The fields u, v, and r should be consistent with the following law of population balance (for any regular subregion X for the defined region during all the time t): where n is the outward unit normal vector to the boundary oX of X. This equation means that the rate of change of population of X addition to the rate at which individuals left X across its boundary should be equal to the rate at which individuals are directly supplied to X [4,5]. Several papers have considered the existence, uniqueness, and regularity of weak solutions [4][5][6][7][8][9] for general degenerate parabolic equations of Eq. (1). Numerical solutions of the biological population equation have seldom been explored and investigated, though some numerical struggles have been done in this field. Modeling of the biological population problem (1) has been explored using the element-free Kp-Ritz method in [1]. An improved element-free Galerkin method for numerical modeling of the biological population problem (1) has also been applied by Zhang et. al. [2]. Shakeri and Dehghan obtained numerical solution of the model using He's variational iteration method [10].
In the current work, we testify spectral meshless radial point interpolation (SMRPI) method [42][43][44] on the problem (1) and then make simulations on two numerical experiments which leads to satisfactory results. A technique based on radial point interpolation is adopted to construct shape functions, also called basis functions, using radial basis functions. These shape functions have delta function property in the frame work of interpolation; therefore, they convince us to impose boundary conditions directly. The time derivatives are approximated by the finite-difference time-stepping method. This method is based on collocation methods with mesh-free techniques as a background. In contrast to those meshless methods based on weak form, there is no integration tools in this approach. Therefore, the computational complexity of SMRPI method seems to be of low order.

The time discretization approximation
In the current work, let us to employ a time-stepping scheme to evaluate the time derivative. To this end, the following first-order finite-difference approximation for the time derivative operator is adopted: Moreover, we apply the following approximations using the Crank-Nicolson technique: rðuðx; tÞÞ ffi where u k ðxÞ ¼ uðx; kDtÞ and x ¼ ðx; yÞ. Applying the above approximation and impose them to the original Eq. (1), we are conducted to the following time discretization equation: The basis functions in the frame of MLRPI Consider a continuous function uðxÞ defined in a domain X, which is represented by a set of field nodes. The uðxÞ at a point of interest x is approximated in the form of where R i ðxÞ is a radial basis function (RBF), n is the number of RBFs, p j ðxÞ is monomial in the space coordinate x, and m is the number of polynomial basis functions. Coefficients a i and b j are unknown which should be determined. In the current work, we use the thin plate spline (TPS) as radial basis functions in Eq. (7) which is defined as follows: To determine a i and b j in Eq. (7), we consider a support domain, such as a disk or square, surrounding the point of interest x and use all nodes included in the support domain to form a system of equations based on Eq. (7). In this way, coefficients of a i and b j are obtained. Therefore, by the idea of interpolation, Eq. (7) is converted to the following form: where / i ðxÞ are called the RPIM shape functions which have the Kronecker delta function property, that is This is because the RPIM shape functions are created to pass thorough nodal values. Moreover, the shape functions are the partitions of unity, that is for more details about RPIM shape functions and the way they are constructed, the readers are referred to see [42].

Operational matrices of high-order derivatives
The essential tools of the current approach is operational matrices which is constructed and described briefly in this section. In fact, the operational matrices make the method easy to apply the high-order differential equations which are really difficult to handle by the majority of methods, especially those techniques based on weak forms. Consider the total N nodes for covering the domain of the problem, i.e., On the other hand, as we noted in the previous section, n is depend on the point of interest x (so, we call it n x ) in Eq. (9) which is the number of nodes included in support domain X x corresponding to the point of interest x. Therefore, we have n x N and Eq. (9) could be reformulated as In fact, there is one shape (basis) function / j ðxÞ, j ¼ 1; 2; 3; . . .; N corresponding to the node x, we define , and then, it is straightforward, from the previous section, to conclude The derivatives of uðxÞ are easily obtained as we indicate above coefficients matrices as Dx and Dy, respectively. In addition, clearly, we propose the following matrix-vector multiplications for high-order derivatives where U ðsÞ y ¼ u ðsÞ y 1 ; u ðsÞ y 2 ; . . .; u ðsÞ Discretization and numerical implementation for the method Equation (6) can be rewritten as To overcome nonlinearity, suppose u k + 1 u u k in the right-hand side of the above equation. This is possible if the time step Dt be sufficiently small, therefore, Eq. (24) can be converted to Now, consider N scattered nodes on the boundary and domain of the problem (1), i.e., X ¼ ðX [ oXÞ. Assuming that uðx i ; kDtÞ; i ¼ 1; 2; . . .; N are known, our purpose is to compute uðx i ; ðk þ 1ÞDtÞ; i ¼ 1; 2; . . .; N. For nodes which are included in the inside of the domain, i.e., x i 2 X, to obtain the discrete system of algebraic equations, let us substitute approximate formulas (12) and (14), (15) into equation (25), then Now, by setting x ¼ x i , i ¼ 1; 2; 3; . . .; N X (N X is the number of nodes in X) in the above equation and then applying notations (19)-(23), we obtain

Enforcement of boundary conditions
There are several techniques to enforcing essential boundary conditions in meshless methods, such as the use of penalty methods, Lagrange multipliers and modified variational principles, etc. In the current work, the essential boundary conditions are imposed directly.

Simply-supported boundary conditions
In the case of simply-supported boundary conditions, we have uðx; tÞ ¼ gðx; tÞ; x 2 oX; 0 t T: For nodes which are located on the boundary oX, we set as where N oX is the number of nodes located on oX.

Clamped boundary conditions
In the case of clamped boundary conditions, it is usually included the following types of boundary conditions: ouðx; tÞ on ¼ gðx; tÞ; x 2 oX; 0 t T: Therefore, in this case, for nodes which are located on the boundary oX, we set as where n ¼ n 1 ðx i Þi þ n 2 ðx i Þj is the outward unit normal to the boundary oX at x i 2 oX. Then, we have the following equations: where N oX is the number of nodes on oX.

Numerical simulations
In this section, we aim to demonstrate that the SMRPI approach has a wider applications by testifying two examples of the type in Eq. (1). To show the accuracy and convergence of the method, maximum absolute error e 1 defined by e 1 ðuÞ ¼ ku exact À u approx k 1 ¼ fju exact ðx i ; tÞ À u approx ðx i ; tÞj : i ¼ 1; 2; . . .; Ng; ð33Þ is used, where fu approx ; p approx g are exact and numerical SMRPI solutions, respectively. In implementing the SMRPI method, we consider support domain as a disk with the radius r s ¼ 4:2h, where h ¼ 0:1 is the distance length between two nodes in x-or y-directions.
Example 1 Considering the following biological population equation: with the initial condition with h ¼ 1 5 , a ¼ 1, and r ¼ 0, and the exact solution is In this problem, we set Dt ¼ 0:00001. Figure 1 shows SMRPI simulations and the absolute errors in some levels of specific values of the time. Figure 2 compares the exact and approximate SMRPI solutions.
Example 2 Considering the following biological population equation: with the initial condition with h ¼ À1, a ¼ 1, r ¼ À 8 9 , and b ¼ 1, and the exact solution is uðx; tÞ ¼ e In this problem, we set Dt ¼ 0:00001 as well. Figure 3 shows SMRPI simulations and the absolute errors in some levels of specific values of the time. Figure 4 compares the exact and approximate SMRPI solutions. As it is seen, the SMRPI and the exact solutions are not distinguishable, while we have adopted a very simple idea to overcome the nonlinearity, as it has been pointed out in Sect. ''Discretization and numerical implementation for the method''.

Conclusions
In this paper, the biological population equation has been investigated using the spectral meshless radial point interpolation method. The shape (basis) functions constructed by radial point interpolation augmented to monomials have been employed to approximate the 2D displacement field. A system of nonlinear discrete equations is obtained through application of the SMRPI. The nonlinear equation system is solved by iteration with a very simple scheme. A mesh-free method does not require a mesh to discretize the domain of the problem under consideration, and the approximate solution is constructed entirely based on a set of scattered nodes. In SMRPI technique, in contrast to those meshless methods based on weak form, there is no integration tools in this approach. Therefore, the computational complexity of SMRPI method seems to be of low order. The numerical results are in excellent agreement with exact solutions.