On the Approximation of a Nonlinear Biological Population Model Using Localized Radial Basis Function Method

Abstract: A localized radial basis function meshless method is applied to approximate a nonlinear biological population model with highly satisfactory results. The method approximates the derivatives at every point corresponding to their local support domain. The method is well suited for arbitrary domains. Compared to the finite element and element free Galerkin methods, no integration tool is required. Four examples are demonstrated to check the efficiency and accuracy of the method. The results are compared with an exact solution and other methods available in literature.


Introduction
Biological population models have attracted many researchers in the last few years. Consider the nonlinear biological population model [1] ρ t (x, t) = ρ 2 xx (x, t) + ρ 2 yy (x, t) + σ(ρ), x = (x, y) ∈ Ω ⊆ R d , d = 1, 2, where σ is a function of ρ(x, t) subject to the initial condition ρ(x, 0), the solution ρ(x, t) denotes the density of population in a region Γ, and σ represent the population supply due to births and deaths. The fields ρ(x, t), V (diffusion velocity) and σ must obey the following population balance law for a regular subregion of region Γ for all time t.
where n denotes the outward unit normal to the boundary ∂R of R. Physically, the equation describes that rate of population directly supplied to R is equal to the sum of the rate of change of population and the rate at which individuals leave R across the boundary [1][2][3]. Several methods have been applied to consider the uniqueness, existence and regularity of these models. In recent work, the biological population model (BPM) has been numerically approximated by various methods [4][5][6]. Meshless procedures have attracted researchers because they can be applied to a complex shaped domain with high dimensional problems. These methods include the radial basis function (RBF) [7][8][9], smooth particle hydrodynamics methods (SPH) [10], reproducing kernel particle method (RKPM) [11], element free Galerkin method (EFG) [11], and meshless local Petrov Galerkin method (MLPG) [12].
Since from discovery of mutiquadrics by Hardy [13], the RBFs (radial basis function) attracted researchers. Thin plate splines were proposed by Douchen [14]. Edwards Kansa for the first time used RBFs to solve PDEs [15], which greatly improved the numerical solution of PDEs. Presently RBFs are applied to solve a lot of mathematical models arising in engineering and science. During the past decades, robust kernel based methods have been developed by Atluri [12], Buhumann [16], Wendland [17], Fasshaur [18,19] and Levesely [20] and many others [21][22][23][24][25]. The convergence theory of the RBFs method was proved by Schaback [26]. The importance of RBFs for solving PDES is easy to implement, robustness and applicability to complicated and multi dimensional problems.
To overcome the problems occurring during global RBFs, partition of unity, domain decomposition method and greedy algorithms have been discovered [26][27][28]. Another alternative is the implementation of local RBFs (LRBFs) [29][30][31]. In LRBFs, small domains centred at each points are used to give rise to differentiation weights. In the present scenario, we have applied local RBFs to solve the nonlinear biological model.

Spatial Approximation
The multivariate data interpolation problem involves recovering an unknown function ρ : where the centres x 1 , ..., x M ∈ Ω and Ω ⊂ R 2 is bounded domain and the centres are taken anywhere in the domain Ω. In local RBF approximation techniques, at every centre x j ∈ Ω, the local interpolant where η j = [η 1 , ..., η m ] is expansion coefficients and x j − x k is the Euclidean distance between x j and x k . Here, φ(r) is a kernel function (Multiquadrics) defined for r ≥ 0 and Ω j ⊂ Ω is the local support domain for each node x j and have m nearest centres with centre x j . The M number of m × m linear systems is where matrix entries of N j are a j kj = φ( x k − x j ), where x k and x j ∈ Ω j . The matrix N j is interpolation matrix corresponding to centre x j . Every system is solved for unknown coefficients corresponding to each local domain Ω j . Applying the differential operator L to Equation (3) Lρ(x j , t) = ∑ The expression in (5) can be written as element wise product of two vectors, where η j is the m × 1 vector of unknown coefficients, and Z j is 1 × m vector with entries To eliminate the unknown coefficients, Equation (4) on substituting the values of η j from (8) in (6) we have , where, gives weight to centre x j . For all centres, it gives where W is M × M sparse differentiation matrix, each row has and M − m zeros elements, where m is the stencil size in each local domain Ω i .

Temporal Approximation
After spatial local RBF approximation, we obtained the following system of ODEs In our case F(ρ) = WP. For discretization any built-in Matlab ODEs solver can be used, such as ode113 or ode45, and ρ 0 has been taken as initial solution vector. Note that ode45 uses an explicit Runge-Kutta (4,5) formula [32]. For computing ρ(t n ), needs solution at the preceding time ρ(t n−1 ). For ODEs computation, the fourth-order Runge-Kutta method has been taken and the time step δt is taken manually [33]. The method of lines approach is the solution of system of ODEs using finite difference formula (FD) in t. The rule of thumb says if all the eigenvalues of the spatial discretization operator scaled using δt, lies in the stability region of the time discretization operator, the method of lines will be stable.
In addition we may use some other efficient methods for time integrations such as exponential integrator schemes [34,35].

Choosing Optimal Shape Parameter
A large number of kernel functions such as Multiquadrics, Thin Plate Splines, Quadratics and Inverse Mutiquadrics can be used. We used Multiquadrics φ(r) = 1 + (rc) 2 with a shape parameter c. The numerical solution of the problem also depends on the parameter c. Finding most suitable value of shape parameter analytically for various RBFs is still an open problem. A large number of papers are available for finding good value of shape parameter (see for example [36] and the references therein). The matrix condition number (MCN) is used to check the sensitivity of a linear system and accuracy of the solution. Good conditioning results require small values of shape parameter and large separation distance between centres. Both results cannot occur at the same time, which is known as the uncertainty principle [37]. The smallest error occurs when system condition number is kept in range of Algorithm 1.

Numerical Experiments
This section is devoted to demonstrate the accuracy and performance of LRBF method by solving the three test problems of the population model (1). The accuracy is tested in terms of maximum absolute error (MAE) and root mean squared error (RMSE) defined by the following equations where ρ(x j , t) and ρ(x j , t) are the exact and the approximate solutions, respectively.
The present problem is solved by the proposed method in the square domain Ω = [0, 1] 2 for various number of nodes in the global domain Ω and local sub-domains Ω i . For time integration, the RK4 method with time step δt = 0.0001 is used. The sparsity of differentiation matrix W and distribution of nodes in the global and local sub-domains corresponding to an interior and boundary points are shown in Figure 1. Figure 2 shows the point wise error and Figure 3 shows the problem solution ρ(x, 0.5) at different times. The results for M = 400 total number of nodes in the global domain, different stencil size ns of local sub-domain, matrix condition number (MCN) and shape parameter are shown in Table 1. The comparison of our method is made with the results in [5] in Table 2, which shows that our results are more accurate and stable.

Example 2
The biological population model defined by (1) with the initial condition with the parameters values h = 1 96 , α = −1, h = 48 and β = 1. The exact solution of this problem given in [6] is The problem is solved in the square domain Ω = [0, 1] 2 , the numerical scheme is advanced in time with step size δt = 0.0001. Figure 4 compares the exact and numerical solutions and Figure 5 shows the problem solution ρ(x, 0.5) at different times.

Example 3
The biological population model defined by (1) with initial condition with the values of parameters h = −1, α = 1, h = −8 9 and β = 1. The exact solution [6] is given by The results are obtained by the present local kernel based method in the domain [0.1] 2 . Figure 6 compares the exact and numerical solution and Figure 7 shows the point wise error and the model solution ρ(x, 0.5) at different times. The results for different M and ns are obtained and the matrix condition number MCN and shape parameter c are shown in Table 3. The comparison of the present method is made with [5] in Table 4.

Example 4
In this problem the present method is applied over irregular domains for example Ameoba and Cassini like shaped domains for solving same model discussed in Example 1 above. The parametric equations of such types irregular boundaries ∂Ω are defined as where for Ameoba-shaped boundary r(θ) is defined by r(θ) = e sin(θ) sin 2 (2θ) + e cos θ cos 2 (2θ), while for Cassini-shaped boundary r(θ) is defined by The irregular domains with inner and boundary stencils are shown in Figure 8. The exact and numerical solutions corresponding to Ameoba-like domain with 100 boundary points and Cassini-shape domain with M = 220 are shown in Figures 9 and 10 and Table 5, respectively.

Concluding Remarks
In the present paper, the nonlinear biological population model has been numerically simulated using localized RBF method. The multiquadric RBF has been used as the radial kernel. The shape parameter has been chosen using the uncertainty principle. The nonlinear equation is converted to a system of ODEs, which is then advanced in time by the RK4 method. The numerical results are compared with the benchmark problems. The present numerical scheme is more accurate and stable as compared to other numerical methods for solving the biological model considered in this work.