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

$$\begin{aligned} \frac{\partial u({\mathbf{x}} ,t)}{\partial t}=\frac{\partial ^2 u^2}{\partial x^2}+\frac{\partial ^2 u^2}{\partial y^2} +\sigma (u),\quad {\mathbf{x}} =(x,y)\in \Omega \subseteq {\mathbb{R}}^2,\quad t\in [0,T], \end{aligned}$$
(1)

with a given initial condition \(u({\mathbf{x}} ,0)\), where u denotes the population density, and \(\sigma\) represents the population supply due to births and deaths. We can consider a more general form for \(\sigma\) as \(hu^a(1-ru^b)\), in which \(h, a, r, b\in \mathbb {R}\). The field \(u({\mathbf{x}} ,t)\) gives the number of individuals’ per-unit volume at position \({\mathbf{x}}\) and time t, and it is integral over any subregion \(\Omega\) to give the total population of \(\Omega\) at time t. The field \(\sigma (u)\) gives the rate at which individuals are supplied (per-unit volume) exactly at \({\mathbf{x}}\) through the births and deaths. The flow of population from point to point is then depicted by the diffusion velocity \(v({\mathbf{x}} ,t)\), which provides the average velocity of those individuals who occupies \({\mathbf{x}}\) at the time t. The fields u, v, and \(\sigma\) should be consistent with the following law of population balance (for any regular subregion \(\Omega\) for the defined region during all the time t):

$$\begin{aligned} \frac{\text {d}}{\text {d} t}\int _\Omega u\text {d}V+\int _{\partial \Omega } u\times v\times n\text {d}V=\int _\Omega \sigma \text {d}V, \end{aligned}$$
(2)

where n is the outward unit normal vector to the boundary \(\partial \Omega\) of \(\Omega\). This equation means that the rate of change of population of \(\Omega\) addition to the rate at which individuals left \(\Omega\) across its boundary should be equal to the rate at which individuals are directly supplied to \(\Omega\) [4, 5].

Several papers have considered the existence, uniqueness, and regularity of weak solutions [4–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 general, the mesh-free methods can be grouped into two categories, first category uses the equation in weak form, for instance, element-free Galerkin method (EFG), improved element-free Galerkin method (IEFG), complex variable element-free Galerkin method (CVEFG), improved complex variable element-free Galerkin method (ICVEFG), meshless local Petrov–Galerkin method, and etc. [11–21], and the second category applies the strong form of the equation, for example, meshless collocation approach [22–28]. In addition, the meshless methods depend on how their shape or basis function constructed are different, and it might be based on interpolation or curve fitting and etc. [29–41].

In the current work, we testify spectral meshless radial point interpolation (SMRPI) method [42–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:

$$\begin{aligned} \frac{\partial u({\mathbf{x}} ,t)}{\partial t}\cong \frac{1}{\Delta t}(u^{k+1}({\mathbf{x}} )-u^{k}({\mathbf{x}} )). \end{aligned}$$
(3)

Moreover, we apply the following approximations using the Crank–Nicolson technique:

$$\begin{aligned}&\frac{\partial ^2 u^2}{\partial x^2}+\frac{\partial ^2 u^2}{\partial y^2}\cong \frac{1}{2} \left( \frac{\partial ^2 (u^{k+1})^2}{\partial x^2}+\frac{\partial ^2 (u^{k})^2}{\partial x^2}+ \frac{\partial ^2 (u^{k+1})^2}{\partial y^2}+\frac{\partial ^2 (u^{k})^2}{\partial y^2}\right) , \end{aligned}$$
(4)
$$\begin{aligned}&\sigma (u({\mathbf{x}} ,t))\cong \frac{1}{2}\left( \sigma (u^{k+1}({\mathbf{x}} ))+\sigma (u^k({\mathbf{x}} ))\right) , \end{aligned}$$
(5)

where \(u^{k}({\mathbf{x}} )=u({\mathbf{x}} ,k\Delta t)\) and \({\mathbf{x}} =(x,y)\). Applying the above approximation and impose them to the original Eq. (1), we are conducted to the following time discretization equation:

$$\begin{aligned} u^{k+1}=u^{k}+\frac{\Delta t}{2}\left\{ \frac{\partial ^2 (u^{k+1})^2}{\partial x^2}+\frac{\partial ^2 (u^{k})^2}{\partial x^2}+\frac{\partial ^2 (u^{k+1})^2}{\partial y^2}+\frac{\partial ^2 (u^{k})^2}{\partial y^2}+\sigma [u^{k+1}]+\sigma [u^k]\right\} . \end{aligned}$$
(6)

The basis functions in the frame of MLRPI

Consider a continuous function \(u({\mathbf{x}} )\) defined in a domain \(\Omega\), which is represented by a set of field nodes. The \(u({\mathbf{x}} )\) at a point of interest \({\mathbf{x}}\) is approximated in the form of

$$\begin{aligned} u({\mathbf{x}} )=\sum _{i=1}^nR_i({\mathbf{x}} )a_i+\sum _{j=1}^mp_j({\mathbf{x}} )b_j= {\mathbf{R}} ^T({\mathbf{x}} ){\mathbf{a}} +{\mathbf{P}} ^T({\mathbf{x}} ){\mathbf{b}} , \end{aligned}$$
(7)

where \(R_i({\mathbf{x}} )\) is a radial basis function (RBF), n is the number of RBFs, \(p_j({\mathbf{x}} )\) is monomial in the space coordinate \({\mathbf{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:

$$\begin{aligned} R({\mathbf{x}} )=r^{2m}\ln (r), \ \ \ m=1, 2, 3,{\ldots } \end{aligned}$$
(8)

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 \({\mathbf{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:

$$\begin{aligned} u({\mathbf{x}} )={\mathbf{{\Phi }}}^T({\mathbf{x}} ){\mathbf{U}} _s=\sum _{i=1}^n\phi _i({\mathbf{x}} )u_i, \end{aligned}$$
(9)

where \(\phi _i({\mathbf{x}} )\) are called the RPIM shape functions which have the Kronecker delta function property, that is

$$\begin{aligned} \phi _i({\mathbf{x}} _j)=\left\{ \begin{array}{ll} 1, &\quad i=j, \quad j=1, 2,{\ldots }, n, \\ 0, &\quad i\ne j, \quad i,j=1, 2,{\ldots }, n. \end{array}\right. \end{aligned}$$
(10)

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

$$\begin{aligned} \sum _{i=1}^n\phi _i({\mathbf{x}} )=1, \end{aligned}$$
(11)

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., \(\overline{\Omega }=\left( \Omega \cup \partial \Omega \right)\). On the other hand, as we noted in the previous section, n is depend on the point of interest \({\mathbf{x}}\) (so, we call it \(n_{\mathbf{x}}\)) in Eq. (9) which is the number of nodes included in support domain \(\Omega _{\mathbf{x}}\) corresponding to the point of interest \({\mathbf{x}}\). Therefore, we have \(n_{\mathbf{x}}\le N\) and Eq. (9) could be reformulated as

$$\begin{aligned} u({\mathbf{x}} )={\mathbf{\Phi}}^T({\mathbf{x}} ){\mathbf{U}}_s=\sum _{j=1}^N\phi _j({\mathbf{x}} )u_j. \end{aligned}$$
(12)

In fact, there is one shape (basis) function \(\phi _j({\mathbf{x}} )\), \(j=1,2,3,{\ldots },N\) corresponding to the node \({\mathbf{x}}\), we define \(\Omega _{\mathbf{x}}^c=\left\{ {\mathbf{x}}_j:{\mathbf{x}} _j\notin \Omega _{\mathbf{x}}\right\}\), and then, it is straightforward, from the previous section, to conclude

$$\begin{aligned} \forall {\mathbf{x}} _j\in \Omega _{\mathbf{x}}^c: \phi _j({\mathbf{x}})=0. \end{aligned}$$
(13)

The derivatives of \(u({\mathbf{x}} )\) are easily obtained as

$$\begin{aligned} \frac{\partial u({\mathbf{x}} )}{\partial x}=\sum _{j=1}^N\frac{\partial \phi _j({\mathbf{x}} )}{\partial x}u_j, \quad \frac{\partial u({\mathbf{x}} )}{\partial y}=\sum _{j=1}^N\frac{\partial \phi _j({\mathbf{x}} )}{\partial y}u_j, \end{aligned}$$
(14)

and also high derivatives of \(u({\mathbf{x}} )\) are clearly obtained as

$$\begin{aligned} \frac{\partial ^{s} u({\mathbf{x}} )}{\partial x^s}=\sum _{j=1}^N\frac{\partial ^{s}\phi _j({\mathbf{x}} )}{\partial x^s}u_j, \quad \frac{\partial ^{s} u({\mathbf{x}} )}{\partial y^s}=\sum _{j=1}^N\frac{\partial ^{s}\phi _j({\mathbf{x}} )}{\partial y^s}u_j, \end{aligned}$$
(15)

where \(\frac{\partial ^{s}(\cdot )}{\partial x^s}\) and \(\frac{\partial ^{s}(\cdot )}{\partial y^s}\) are the sth derivative with respect to x and y, respectively. By indicating \(u_x^{(s)}(\cdot )=\frac{\partial ^{s}(\cdot )}{\partial x^s}\) and \(u_y^{(s)}(\cdot )=\frac{\partial ^{s}(\cdot )}{\partial y^s}\), and substituting \({\mathbf{x}} ={\mathbf{x}} _i\) in Eq. (14), we can formulate the discrete differentiations process as a matrix-vector multiplications

$$\begin{aligned} \left( \begin{array}{l} u_{x_1}^\prime \\ u_{x_2}^\prime \\ \vdots \\ u_{x_N}^\prime \\ \end{array} \right) =\left( \begin{array}{llll} \frac{\partial \phi _1({\mathbf{x}} _1)}{\partial x} &{} \frac{\partial \phi _2({\mathbf{x}}_1)}{\partial x} &{} \cdots &{} \frac{\partial \phi _N({\mathbf{x}} _1)}{\partial x} \\ \frac{\partial \phi _1({\mathbf{x}} _2)}{\partial x} &{} \frac{\partial \phi _2({\mathbf{x}} _2)}{\partial x} &{} \cdots &{} \frac{\partial \phi _N({\mathbf{x}} _2)}{\partial x} \\ \vdots &{} \vdots &{} \ddots &{} \vdots \\ \frac{\partial \phi _1({\mathbf{x}} _N)}{\partial x} &{} \frac{\partial \phi _2({\mathbf{x}} _N)}{\partial x} &{} \cdots &{} \frac{\partial \phi _N({\mathbf{x}} _N)}{\partial x} \\ \end{array} \right) \left( \begin{array}{l} u_1 \\ u_2 \\ \vdots \\ u_N \\ \end{array} \right) , \end{aligned}$$
(16)
$$\begin{aligned} \left( \begin{array}{l} u_{y_1}^\prime \\ u_{y_2}^\prime \\ \vdots \\ u_{y_N}^\prime \\ \end{array} \right) =\left( \begin{array}{llll} \frac{\partial \phi _1({\mathbf{x}} _1)}{\partial y} &{} \frac{\partial \phi _2({\mathbf{x}} _1)}{\partial y} &{} \cdots &{} \frac{\partial \phi _N({\mathbf{x}} _1)}{\partial y} \\ \frac{\partial \phi _1({\mathbf{x}} _2)}{\partial y} &{} \frac{\partial \phi _2({\mathbf{x}} _2)}{\partial y} &{} \cdots &{} \frac{\partial \phi _N({\mathbf{x}} _2)}{\partial y} \\ \vdots &{} \vdots &{} \ddots &{} \vdots \\ \frac{\partial \phi _1({\mathbf{x}} _N)}{\partial y} &{} \frac{\partial \phi _2({\mathbf{x}} _N)}{\partial y} &{} \cdots &{} \frac{\partial \phi _N({\mathbf{x}} _N)}{\partial y} \\ \end{array} \right) \left( \begin{array}{l} u_1 \\ u_2 \\ \vdots \\ u_N \\ \end{array} \right) , \end{aligned}$$
(17)

we indicate above coefficients matrices as \({\mathbf{D}}x\) and \({\mathbf{D}}y\), respectively. In addition, clearly, we propose the following matrix-vector multiplications for high-order derivatives

$$\begin{aligned} U_x^{(s)}={\mathbf{D}}^{(s)}xU, \quad U_y^{(s)}={\mathbf{D}}^{(s)}yU, \end{aligned}$$
(18)

where

$$\begin{aligned}&U_x^{(s)}=\left( u_{x_1}^{(s)},u_{x_2}^{(s)},{\ldots },u_{x_N}^{(s)}\right) ^T, \end{aligned}$$
(19)
$$\begin{aligned}&U_y^{(s)}=\left( u_{y_1}^{(s)},u_{y_2}^{(s)},{\ldots },u_{y_N}^{(s)}\right) ^T, \end{aligned}$$
(20)
$$\begin{aligned}&{\mathbf{D}}^{(s)}x_{ij}=\frac{\partial ^s\phi _j({\mathbf{x}} _i)}{\partial x^s}, \end{aligned}$$
(21)
$$\begin{aligned}&{\mathbf{D}}^{(s)}y_{ij}=\frac{\partial ^s\phi _j({\mathbf{x}} _i)}{\partial y^s}, \end{aligned}$$
(22)
$$\begin{aligned}&U=(u_1,u_2,{\ldots },u_N)^T. \end{aligned}$$
(23)

Discretization and numerical implementation for the method

Equation (6) can be rewritten as

$$\begin{aligned} u^{k+1}&=u^{k}+\Delta t\left\{ u^{k+1}\frac{\partial ^2 u^{k+1}}{\partial x^2}+\left( \frac{\partial u^{k+1}}{\partial x}\right) ^2+u^{k}\frac{\partial ^2 u^{k}}{\partial x^2}+\left( \frac{\partial u^{k}}{\partial x}\right) ^2\right\} \nonumber \\&\quad +\Delta t\left\{ u^{k+1}\frac{\partial ^2 u^{k+1}}{\partial y^2}+\left( \frac{\partial u^{k+1}}{\partial y}\right) ^2+u^{k}\frac{\partial ^2 u^{k}}{\partial y^2}+\left( \frac{\partial u^{k}}{\partial y}\right) ^2\right\} \nonumber \\&\quad +\frac{\Delta t}{2}\left\{ \sigma [u^{k+1}]+\sigma [u^k]\right\} . \end{aligned}$$
(24)

To overcome nonlinearity, suppose u k + 1 ≊ u k in the right-hand side of the above equation. This is possible if the time step \(\Delta t\) be sufficiently small, therefore, Eq. (24) can be converted to

$$\begin{aligned} u^{k+1}=u^{k}+2\Delta t\left\{ u^{k}\frac{\partial ^2 u^{k}}{\partial x^2}+ \left( \frac{\partial u^{k}}{\partial x}\right) ^2+u^{k} \frac{\partial ^2 u^{k}}{\partial y^2}+\left( \frac{\partial u^{k}}{\partial y}\right) ^2\right\} + \Delta t\sigma [u^k]. \end{aligned}$$
(25)

Now, consider N scattered nodes on the boundary and domain of the problem (1), i.e., \(\overline{\Omega }=(\Omega \cup \partial \Omega )\). Assuming that \(u({\mathbf{x}} _i,k\Delta t),\ i= 1,2,{\ldots },N\) are known, our purpose is to compute \(u({\mathbf{x}} _i,(k+1)\Delta t),\ i= 1,2,{\ldots },N\). For nodes which are included in the inside of the domain, i.e., \({\mathbf{x}} _i \in \Omega\), to obtain the discrete system of algebraic equations, let us substitute approximate formulas (12) and (14), (15) into equation (25), then

$$\begin{aligned} u^{k+1}&=u^{k}+2\Delta t\left\{ \left( \sum _{j=1}^N\phi _j({\mathbf{x}} )u_j^{k}\right) \left( \sum _{j=1}^N\frac{\partial ^{2}\phi _j({\mathbf{x}} )}{\partial x^2}u_j^{k}\right) + \left( \sum _{j=1}^N\frac{\partial \phi _j({\mathbf{x}} )}{\partial x}u_j^{k}\right) ^2\right\} \nonumber \\&\quad +2\Delta t\left\{ \left( \sum _{j=1}^N\phi _j({\mathbf{x}} )u_j^{k}\right) \left( \sum _{j=1}^N\frac{\partial ^{2}\phi _j({\mathbf{x}} )}{\partial y^2}u_j^{k}\right) + \left( \sum _{j=1}^N\frac{\partial \phi _j({\mathbf{x}} )}{\partial y}u_j^{k}\right) ^2\right\} \nonumber \\&\quad +\Delta t\sigma \left[ \left( \sum _{j=1}^N\phi _j({\mathbf{x}} )u_j^{k}\right) \right] . \end{aligned}$$
(26)

Now, by setting \({\mathbf{x}} ={\mathbf{x}} _i\), \(i=1,2,3,{\ldots },N_\Omega\) (\(N_\Omega\) is the number of nodes in \(\Omega\)) in the above equation and then applying notations (19)–(23), we obtain

$$\begin{aligned} u_i^{k+1}&=u_i^{k}+2\Delta t\left\{ u_i\left( \sum _{j=1}^N{\mathbf{D}}^{(2)}x_{ij}u_j^{k}\right) + \left( \sum _{j=1}^N{\mathbf{D}}x_{ij}u_j^{k}\right) ^2\right\} \nonumber \\&\quad +2\Delta t\left\{ u_i\left( \sum _{j=1}^N{\mathbf{D}}^{(2)}y_{ij}u_j^{k}\right) + \left( \sum _{j=1}^N{\mathbf{D}}y_{ij}u_j^{k}\right) ^2\right\} \nonumber \\&\quad +\Delta t\sigma \left[ u_i\right] . \end{aligned}$$
(27)

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

$$\begin{aligned} u({\mathbf{x}} ,t)=g({\mathbf{x}} ,t) ,\quad x\in \partial \Omega ,\quad 0\le t\le T. \end{aligned}$$
(28)

For nodes which are located on the boundary \(\partial \Omega\), we set as

$$\begin{aligned} u^{k+1}({\mathbf{x}} _i)=g({\mathbf{x}} _i,(k+1)\Delta t),\quad i=1,2,{\ldots },N_{\partial \Omega }, \end{aligned}$$
(29)

where \(N_{\partial \Omega }\) is the number of nodes located on \(\partial \Omega\).

Clamped boundary conditions

In the case of clamped boundary conditions, it is usually included the following types of boundary conditions:

$$\begin{aligned} \frac{\partial u({\mathbf{x}} ,t)}{\partial n}=g({\mathbf{x}} ,t) ,\quad x\in \partial \Omega , \quad 0\le t\le T. \end{aligned}$$
(30)

Therefore, in this case, for nodes which are located on the boundary \(\partial \Omega\), we set as

$$\begin{aligned} {\mathbf{n}} _1({\mathbf{x}} _i)\frac{\partial u}{\partial x}\left( {\mathbf{x}} _i,k+1\triangle t\right) +{\mathbf{n}} _2({\mathbf{x}} _i) \frac{\partial u}{\partial y}\left( {\mathbf{x}} _i,(k+1)\Delta t\right) =g({\mathbf{x}} _i,(k+1)\Delta t), \end{aligned}$$
(31)

where \({\mathbf{n}} ={\mathbf{n}} _1({\mathbf{x}} _i)i+{\mathbf{n}} _2({\mathbf{x}} _i)j\) is the outward unit normal to the boundary \(\partial \Omega\) at \({\mathbf{x}} _i\in \partial \Omega\). Then, we have the following equations:

$$\begin{aligned} {\mathbf{n}} _1({\mathbf{x}} _i)\sum _{j=1}^N{\mathbf{D}}x_{ij}u_j^{k+1}+{\mathbf{n}} _2({\mathbf{x}} _i) \sum _{j=1}^N{\mathbf{D}}y_{ij}u_j^{k+1}=g({\mathbf{x}} _i,(k+1)\Delta t),\ \ \ i=1,2,{\ldots },N_{\partial \Omega }, \end{aligned}$$
(32)

where \(N_{\partial \Omega }\) is the number of nodes on \(\partial \Omega\).

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 \(\varepsilon _{\infty }\) defined by

$$\begin{aligned} \varepsilon _{\infty }(u)=\Vert u_\mathrm{{exact}}-u_\mathrm{{approx}}\Vert _{\infty }=\lbrace \vert u_\mathrm{{exact}} ({\mathbf{x}} _i,t)-u_\mathrm{{approx}}({\mathbf{x}} _i,t)\vert : i=1, 2,{\ldots }, N\rbrace , \end{aligned}$$
(33)

is used, where \(\lbrace u_\mathrm{{approx}}, p_\mathrm{{approx}}\rbrace\) 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:

$$\begin{aligned} \frac{\partial u({\mathbf{x}} ,t)}{\partial t}=\frac{\partial ^2 u^2}{\partial x^2}+\frac{\partial ^2 u^2}{\partial y^2} +hu^a(1-ru^b), \quad {\mathbf{x}} =(x,y)\in [0,1]^2, \quad t\in [0,T], \end{aligned}$$
(34)

with the initial condition

$$\begin{aligned} u({\mathbf{x}} ,0)=\sqrt{xy},\quad {\mathbf{x}} \in [0,1]^2, \end{aligned}$$
(35)

with \(h=\frac{1}{5}\), \(a=1\), and \(r=0\), and the exact solution is

$$\begin{aligned} u({\mathbf{x}} ,t)=e^{\frac{t}{5}}\sqrt{xy},\quad {\mathbf{x}} \in [0,1]^2. \end{aligned}$$
(36)

In this problem, we set \(\Delta t = 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.

Fig. 1
figure 1

SMRPI simulations and the absolute errors for Example 1

Full size image
Fig. 2
figure 2

Comparison of SMRPI and exact solutions for Example 1

Full size image

Example 2

Considering the following biological population equation:

$$\begin{aligned} \frac{\partial u({\mathbf{x}} ,t)}{\partial t}=\frac{\partial ^2 u^2}{\partial x^2}+ \frac{\partial ^2 u^2}{\partial y^2} +hu^a(1-ru^b),\quad {\mathbf{x}} =(x,y)\in [0,1]^2,\quad t\in [0,T], \end{aligned}$$
(37)

with the initial condition

$$\begin{aligned} u({\mathbf{x}} ,0)=e^{\frac{x+y}{3}},\quad {\mathbf{x}} \in [0,1]^2, \end{aligned}$$
(38)

with \(h=-1\), \(a=1\), \(r=-\frac{8}{9}\), and \(b=1\), and the exact solution is

$$\begin{aligned} u({\mathbf{x}} ,t)=e^{\frac{1}{3}(x+y)-t},\quad {\mathbf{x}} \in [0,1]^2. \end{aligned}$$
(39)

In this problem, we set \(\Delta t = 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".

Fig. 3
figure 3

SMRPI simulations and the absolute errors for Example 2

Full size image
Fig. 4
figure 4

Comparison of SMRPI and exact solutions for Example 2

Full size image

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.