Dynamical Analysis of Two-Dimensional Fractional-Order-in-Time Biological Population Model Using Chebyshev Spectral Method

Abstract

In this study, we investigate the application of fractional calculus to the mathematical modeling of biological systems, focusing on fractional-order-in-time partial differential equations (FTPDEs). Fractional derivatives, especially those defined in the Caputo sense, provide a useful tool for modeling memory and hereditary characteristics, which are problems that are frequently faced with integer-order models. We use the Chebyshev spectral approach for spatial derivatives, which is known for its faster convergence rate, in conjunction with the $L1$ scheme for time-fractional derivatives because of its high accuracy and robustness in handling nonlocal effects. A detailed theoretical analysis, followed by a number of numerical experiments, is performed to confirmed the theoretical justification. Our simulation results show that our numerical technique significantly improves the convergence rates, effectively tackles computing difficulties, and provides a realistic simulation of biological population dynamics.

1. Introduction

Fractional calculus is an extension of conventional calculus that explores the possibility of taking real-numbered powers of the differentiation and integration operators. Unlike classical calculus, which only works with integer orders, fractional calculus accepts integrals and derivatives of any fractional order, positive or negative. The fractional integral of order $\alpha$ is a generalization of the n-fold integral applied n times, where $\alpha$ can be any real number. It is typically defined using a convolution with a power function, leading to an integral operator known as the Riemann–Liouville integral. Fractional calculus has a wide range of applications in disciplines including biology, engineering, physics, economics, and more. It is employed in physics to simulate systems that have memory and inherited characteristics, such as viscoelastic materials or processes in which the system’s future state is dependent on both its past behavior and its current state. Fractional calculus is helpful to engineers in control theory, especially for designing controllers that provide reliable performance under a variety of circumstances [1,2,3]. Fractional derivatives are incorporated into standard partial differential equations in order to generate fractional partial differential equations (FPDEs), which enable the modeling of memory and hereditary features in a variety of materials and systems. Time-fractional differential equations, in which the order of the time derivative is a fraction instead of an integer, are the subject of this research. When modeling systems where the rate of change is not constant but rather depends on the process’s whole history, these equations are especially helpful. Most recently, the application of fractional-order equation-based vegetation–water in an arid flat environment has been studied in [4].

The most common form of a time-fractional partial differential equation involves the Caputo or Riemann–Liouville definition of fractional derivatives. Time-fractional PDEs are extensively used to model physical phenomena exhibiting anomalous diffusion, which is diffusion where the flux at any point in space depends on the historical states of the system, rather than being directly proportional to the gradient of concentration as in classical diffusion. These models have been used to simulate subdiffusion or superdiffusion processes in complex systems where particle trajectories diverge from those of traditional Brownian motion in domains like physics. When molecular movement is impeded by barriers or binding effects that lead to anomalous diffusion, biology uses it to describe the transport dynamics within cells. In finance, on the other hand, it is employed to capture memory and hereditary effects in financial markets, enhancing asset price and risk management modeling [5,6,7,8].

The goal of the latest developments in fractional calculus has been to enhance the numerical techniques for solving fractional differential equations, which are crucial for real-world applications across a range of industries. Several novel techniques have been presented by researchers in an effort to more effectively handle the nonlocal characteristics included in fractional derivatives. Adaptive methods are among them; they improve accuracy and save computing costs by modifying the time steps based on the behavior of the answer. The ability to numerically solve differential equations with fractional derivatives, which are more difficult by nature because of their nonlocal nature, has been made possible by advances in computational mathematics. Modern numerical techniques and algorithms are still developing, which makes fractional calculus more useful in both industry and scientific study [9,10,11,12,13,14]. As fractional calculus advances, it will offer a more comprehensive mathematical framework that goes beyond the constraints of traditional differential and integral calculus, providing up new possibilities for modeling and understanding the natural world.

Biological population models are crucial and have contributed significantly to our understanding of ecology, public health, and the behavior of the environment. Because they aid in forecasting and controlling species interactions, population growth, and the impact of changing environmental circumstances on ecosystems, these models are crucial to ecology. They serve as crucial tools in conservation biology as well, helping scientists identify a species’ vulnerability and develop effective conservation strategies. In order to evaluate programs, predict the beginning of epidemics, and comprehend how infectious diseases propagate, biological population models are essential to public health research [15,16]. A more comprehensive study on the role of fractional calculus in modeling a biological phenomena can be found in [17]. Fractional derivatives are used in biology to simulate anomalous diffusion, which deviates from the standard equations of diffusion in cells and tissues due to its ability to capture memory and hereditary properties inherent in such systems as compared to the integer-order differential equation. Due to the fact that these mathematical models, which consist of fractional-order derivatives, have a nonlocal nature, involve complicated boundary conditions and memory effects, and therefore lack closed-form solutions, the analytical solution of fractional differential equations is more complicated, and sometimes it is even not possible to find it. These models are especially helpful in ecosystems where organisms have lengthy lifespans or where environmental changes occur slowly, which reduces the efficacy of standard models.

A sophisticated method of applying fractional calculus to study biological population dynamics is introduced by the idea of a two-dimensional fractional-order-in-time biological population model. By utilizing fractional derivatives, which are expansions of ordinary derivatives, this mathematical framework expands upon the traditional population models by including memory and hereditary features. This model takes into account both the history and present population sizes when determining each population’s growth rate. This is accomplished through the use of fractional differential equations, in which the degree to which the past influences the dynamics of the present is indicated by the derivative’s order, typically between 0 and 1. These models are better able to capture the intricacies of biological processes that, as is common in many ecological systems, display long-term memory or power-law waiting times. Predators and prey are two examples of interacting animals that could be included in a two-dimensional model. Fractional derivatives are useful in modeling scenarios in which the dynamics of both populations in the future are influenced throughout time by the residual effects of past interactions (e.g., resource competition or predation rates). When contrasted to traditional models, this may result in oscillations and stability conditions that are more realistic. In general, biologists and ecologists can benefit greatly from the use of two-dimensional fractional-order models in biological populations since they offer a greater understanding of the dynamics of intricate ecosystems where past conditions have a substantial influence on future states. Analytically solving population models is often infeasible for complex boundary and initial conditions; thus, numerical solutions such as finite difference methods, finite element methods, and spectral methods are employed. These methods involve discretizing the time and space variables and approximating the fractional derivatives using techniques adapted for their nonlocal nature. The integration of fractional calculus into partial differential equations offers a powerful tool for exploring and describing dynamic systems with memory and hereditary characteristics, which are not adequately modeled by classical PDEs for these models. The development of robust numerical methods continues to be a vital area of research, enabling the practical application of FPDEs across various scientific disciplines [18,19].

Fractional-order diffusion equations, which are generalizations of classical diffusion equations and address super-diffusive flow processes, are among the most basic instances of the former. A large portion of published work on FPDEs has been focused on fractional diffusion equations. In [20], Agrawal uses a finite sine transform method to find the general solution of fractional-order diffusion equations. A theoretical framework using the least-squares finite element technique has been investigated in [21], while a spectral collocation scheme for two-dimensional nonlinear fractional diffusion equations and a radial basis function approximation method is used in [22,23], respectively. A mathematical model based on nonlinear fractional-order equations for the description of the behavior of viscoplastic materials was developed in [24], while a fractional advection–dispersion equation has been numerically investigated in [25,26]. In order to find the scale-invariant solution of the TFDE in terms of the Wright function, Gorenflo et al. employed the similarity approach and the Laplace transform method [27,28]. Numerous authors have examined these models in analytical and numerical frameworks. A few of these researchers have attempted to find analytical solutions for differential equations using temporal fractions. For instance, time-fractional diffusion-wave equations were taken into consideration by Schneider and Wyss and Wyss [29,30]. The temporal fractional PDE using finite difference in fractional time and a higher numerical scheme was investigated in detail in [31]. Capturing memory and hereditary properties is crucial in biological systems due to their inherent dependence on past states. Fractional differential equations are particularly adept at modeling these characteristics because they incorporate nonlocal properties, meaning the rate of change at any point in time depends on all previous states. However, these equations are complex and often lack closed-form solutions, making analytical solutions challenging and necessitating robust numerical methods. The complexity of fractional derivatives, especially their nonlocal nature, poses significant computational challenges that need to be addressed with efficient numerical techniques.

In this work, we consider solving numerically a fractional-order-in-time biological population model in two dimensions. We use the $L1$ scheme for fractional-order-in-time derivatives. This scheme is very useful in handling the nonlocal properties of fractional-order-in-time derivatives. For spatial derivatives, we use an efficient numerical scheme based on the Chebyshev spectral method. Spectral methods are a class of techniques used in numerical analysis to solve differential equations whose solutions are essentially represented by the sum of globally defined basis functions. These functions are typically orthogonal or trigonometric polynomials, depending on the boundary conditions and nature of the problem. When the solution has smooth properties, spectral methods are well known for their high accuracy and exponential rates of convergence, which greatly exceed the performance provided by conventional finite difference or finite element methods. The most notable method among these is the Chebyshev spectral method due to its particular effectiveness in handling problems with complex geometries or boundary conditions. This method utilizes Chebyshev polynomials, which are a set of orthogonal polynomials defined on the interval $[-1,1]$. The main advantage of Chebyshev polynomials is their near-minimal property for polynomial approximation. This property significantly minimizes the maximum error between the numerical solution and the true solution. The efficiency of this approach is further increased by the employment of Chebyshev nodes, or the roots of these polynomials, which assist prevent the Runge phenomenon and guarantee stability and dependability in numerical approximations [32,33,34,35,36,37,38,39,40]. Several studies have applied fractional calculus to model biological systems. For instance, fractional derivatives have been used to simulate anomalous diffusion in cells and tissues, capturing more realistic dynamics than traditional models. Hattaf and Yousfi explored the global stability of fractional diffusion equations in biological systems, highlighting the advantages of fractional models in capturing complex biological behaviors [41]. However, many of these studies focus on one-dimensional models or rely on analytical solutions, which are not feasible for complex systems. This paper addresses this gap by developing a numerical scheme for a two-dimensional fractional-order-in-time biological population model [42,43].

The rest of the paper is organized as follows: In Section 2, we present some preliminaries and some basic definitions, which play a key role in the analysis of the scheme, followed by the mathematical description of the model and the discretization scheme. Section 4 consists of a detailed error analysis, followed by numerical examples in Section 5. Section 6 consists of concluding remarks.

2. Preliminaries and Some Basic Definitions

In this section, we present some basic definitions related to fractional calculus theory and orthogonal polynomials that are useful in the error analysis of our proposed numerical scheme [44,45,46,47].

Definition 1.

A real function $g\left(\tau \right)$, $\tau >0$ is defined to belong to the space $C_{\nu }$, $\nu \in \mathbb{R}$, if there exists a real number $q>\nu$, such that $g\left(\tau \right)={\tau }^q{g}_1\left(\tau \right)$, where $g_1\left(\tau \right)\in C(0,\infty )$, and g is in the space $C_{n\nu }$ if and only if $g^{\left(n\right)}\in C_{\nu }$, $n\in \mathbb{N}$.

Definition 2.

The Riemann–Liouville fractional integral operator $I^{\beta }$ ($\beta \ge 0$) of a function $g\in C_{\nu }$, $\nu \ge -1$ is given by

$$I^{\beta }g\left(\tau \right)=\frac{1}{\Gamma \left(\beta \right)}{\int }_0^{\tau }{(\tau -s)}^{\beta -1}g\left(s\right)ds,(\beta \ge 0)$$

(1)

where $\Gamma (·)$ is the gamma function. The properties of the operator $I^{\beta }$ include

$$I^{\beta }{I}^{\gamma }g\left(\tau \right)=I^{\beta +\gamma }g\left(\tau \right),(\beta \ge 0,\gamma \ge 0)$$

(2)

$$I^{\beta }{\tau }^{\delta }=\frac{\Gamma (1+\delta )}{\Gamma (1+\delta +\beta )}{\tau }^{\beta +\delta },(\delta \ge -1)$$

(3)

Definition 3.

The Caputo fractional derivative $C^{\alpha }$ of a function $g\left(\tau \right)$ is defined as

$${ }^C{D}_{\tau }^{\alpha }g\left(\tau \right)=\frac{1}{\Gamma (n-\alpha )}{\int }_0^{\tau }{g}^{\left(n\right)}\left(s\right)\frac{ds}{{(\tau -s)}^{\alpha +1-n}},(n-1<\mathit{Re}\left(\alpha \right)\le n,n\in \mathbb{N})$$

(4)

The properties of the Caputo fractional derivative include

$${ }^C{D}_{\tau }^{\alpha }{\tau }^{\beta }=\frac{\Gamma (1+\beta )}{\Gamma (1+\beta -\alpha )}{\tau }^{\beta -\alpha },$$

(5)

$$\left({I^{\alpha }}^C{D}^{\alpha }\right)g\left(\tau \right)=g\left(\tau \right)-\sum _{k=0}^{n-1}\frac{g^{\left(k\right)}(0+)}{k!}{\tau }^k,$$

(6)

The Caputo fractional derivative is chosen for its ability to incorporate conventional initial and boundary conditions in problem formulations.

Definition 4.

The Mittag–Leffler function $E_{\alpha ,\beta }\left(z\right)$ for complex numbers z and parameters α and β is defined by the series

$$E_{\alpha ,\beta }\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\Gamma (\alpha k+\beta )},$$

(7)

where $\alpha >0$ and β are real numbers, and $\Gamma (·)$ denotes the Gamma function. This function generalizes the exponential function, which is recovered as a special case when $\alpha =1$ and $\beta =1$.

$E_{\alpha ,\beta }\left(z\right)$ represents a two-parameter family of functions, where $\alpha$ and $\beta$ allow for various forms of convergence and divergence characteristics depending on their values. When $\alpha =1$ and $\beta =1$, the Mittag–Leffler function simplifies to the exponential function $e^z$, that is $E_{1,1}\left(z\right)=e^z$.

Definition 5.

Chebyshev polynomials, represented as $C_n\left(x\right)$, are established over the domain $[-1,1]$ and can be formulated by the expression

$$C_n\left(x\right)=cos(narccos\left(x\right)),n=0,1,2,\dots$$

(8)

Definition 6.

The collection of Chebyshev polynomials $\left\{C_n\left(x\right)\right\}$ set within the interval $[-1,1]$ adhere to this orthogonality condition with respect to the weighted scalar product

$${\langle C_i,C_j\rangle }_{\omega }:={\int }_{-1}^1{C}_i\left(x\right)C_j\left(x\right)\omega \left(x\right)dx=\left\{\begin{array}{cc}0,\hfill & \mathit{if}i\ne j,\hfill \\ \pi ,\hfill & \mathit{if}i=j=0,\hfill \\ \frac{\pi }{2},\hfill & \mathit{if}i=j\ne 0.\hfill \end{array}\right.$$

(9)

where the weighting function $\omega \left(x\right)$ is defined by

$$\omega \left(x\right)=\frac{1}{\sqrt{1-x^2}}.$$

(10)

Definition 7.

Let ${\delta }_n$ denote the scaling coefficients for orthogonality:

$${\delta }_n=\left\{\begin{array}{cc}2,\hfill & \mathit{if}n=0,\hfill \\ 1,\hfill & \mathit{if}n>0.\hfill \end{array}\right.$$

(11)

Considering practical approaches in mathematical representations, when examining polynomials of a degree at most N, the matrix representation of the weighted scalar product can be succinctly illustrated as

$$H=\mathrm{diag}\left\{h_{ii}\right\},h_{ii}:=\frac{\sigma }{2}{\delta }_n$$

(12)

Here, the matrix H is diagonal with size $(N+1)\times (N+1)$, and each diagonal element $h_{ii}$ is influenced by the coefficient ${\delta }_n$ from the previous equation.

Definition 8.

A smooth, continuous function $\sigma \left(\theta \right)$, defined on $[-1,1]$, can be approximated using Chebyshev polynomials $C_n\left(\theta \right)$ as follows:

$$\sigma \left(\theta \right)\approx \sum _{n=0}^N{\tilde{\sigma }}_n{C}_n\left(\theta \right),$$

(13)

where N represents the truncation level in the spectral method, and ${\tilde{\sigma }}_n$ is the Chebyshev expansion coefficient for σ. This coefficient is calculated by

$${\tilde{\sigma }}_n=\frac{1}{{\delta }_n}{\int }_{-1}^1\sigma \left(\theta \right)C_n\left(\theta \right)\omega \left(\theta \right)d\theta ,{\delta }_n={\int }_{-1}^1{C}_n^2\left(\theta \right)\sqrt{1-{\theta }^2}d\theta =\left\{\begin{array}{cc}\pi ,\hfill & \mathit{if}n=0,\hfill \\ \frac{\pi }{2},\hfill & \mathit{if}n>0.\hfill \end{array}\right.$$

(14)

Definition 9.

In the computational process known as the Chebyshev forward transformation, the Gauss–Lobatto integration method is typically employed for computing the weighted integral:

$${\tilde{\sigma }}_n\approx \frac{1}{{\delta }_n}\sum _{i=0}^N\sigma \left({\theta }_i\right)C_n\left({\theta }_i\right)w_i$$

(15)

where N represents the number of integration nodes, matching the truncation level. These nodes and weights are specified as

$${\theta }_0=1,{\theta }_N=-1,{\theta }_i=cos\left(\frac{\pi i}{N}\right),w_0=w_N=\frac{\pi }{2N},w_i=\frac{\pi }{N}$$

(16)

Definition 10.

For any continuous and smooth function $\sigma \left(\theta \right)$ defined over the interval $[-1,1]$, if the derivative ${\sigma }^{\prime }\left(\theta \right)$ maintains smoothness, it can be expanded using $C_n\left(\theta \right)$ as

$${\sigma }^{\prime }\left(\theta \right)\approx \sum _{n=0}^N{\tilde{\sigma }}_n^{\prime }{C}_n\left(\theta \right),$$

(17)

It is established that the expansion coefficient ${\tilde{\sigma }}_n^{\prime }$ for the derivative ${\sigma }^{\prime }\left(\theta \right)$ and the coefficient ${\tilde{\sigma }}_n$ for the original function satisfy the relation

$${\tilde{\sigma }}_n^{\prime }\approx \frac{2}{{\delta }_n}\sum _{q=n+1,q+n\mathit{odd}}^{N}q{\tilde{\sigma }}_q,n\ge 0$$

(18)

This relation can then be reformulated in matrix terminology:

$${\tilde{\sigma }}^{\prime }=D\tilde{\sigma },$$

(19)

where ${\tilde{\sigma }}^{\prime }$ and $\tilde{\sigma }$ are arrays with dimensions $1\times (N+1)$, defined as ${\tilde{\sigma }}^{\prime }=[{\tilde{\sigma }}_0^{\prime },{\tilde{\sigma }}_1^{\prime },\dots ,{\tilde{\sigma }}_N^{\prime }]$ and $\tilde{\sigma }=[{\tilde{\sigma }}_0,{\tilde{\sigma }}_1,\dots ,{\tilde{\sigma }}_N]$. The matrix D is an upper triangular square matrix, sized $(N+1)\times (N+1)$.

Definition 11.

The natural Sobolev norms, appropriate for gauging approximation errors within the Chebyshev framework, integrate the Chebyshev weight into the quadratic means of the error and its derivatives over the span $(-1,1)$. Thus, we establish the weighted Sobolev norm as

$${\parallel g\parallel }_{{\mathcal{H}}_m^w(-1,1)}={\left(\sum _{n=0}^m{\parallel g^{\left(n\right)}\parallel }_{L_w^2(-1,1)}^2\right)}^{1/2}.$$

(20)

The associated Hilbert space is denoted ${\mathcal{H}}_m^w(-1,1)$, where:

• $\parallel g^{\left(n\right)}{\parallel }_{L_w^2(-1,1)}$ signifies the $L^2$-norm of the n-th derivative of g, weighted over the interval $(-1,1)$.
• ${\mathcal{H}}_m^w(-1,1)$ is the weighted Sobolev space capturing the behavior of functions and their derivatives up to order m under the weighted norm.

3. Fractional-Order-in-Time Dispersal in Population Dynamics

According to biologists, migration or dispersal has a significant impact on the regulation of species populations. The diffusion of a biological species in a given region C is described by three functions of position $\mathbf{x}=(\eta ,\xi )$ and time t: the population supply $s(\mathbf{x},t)$, the diffusion velocity $\mathbf{u}(\mathbf{x},t)$, and the population density $\rho (\mathbf{x},t)$. The population density $\rho (\mathbf{x},t)$ quantifies the number of individuals per unit volume at position $\mathbf{x}$ and at time t. By integrating $\rho$ over any subregion D of C, the total population in D at time t is obtained. The term $s(\mathbf{x},t)$ indicates the rate at which individuals are added or removed per unit volume at position $\mathbf{x}$ due to births and deaths. The diffusion velocity $\mathbf{u}(\mathbf{x},t)$ represents the mean velocity of individuals at position $\mathbf{x}$ at time t, facilitating the description of population movement from one location to another.

The following equation governs the dynamics:

$$\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}{\int }_D\rho dV+{\int }_{\partial D}\rho \mathbf{u}·\mathbf{n}dA={\int }_{D}sdV,$$

(21)

where $\mathbf{n}$ is the outward unit normal on the boundary $\partial D$. The Caputo fractional derivative is used to interpret the derivative. This basic equation means that the rate at which new individuals are delivered directly to D must equal the rate at which the population within D changes plus the rate at which people leave D over its boundary. With $s=s\left(\rho \right)$ and $\mathbf{u}=\lambda \left(\rho \right)\nabla \rho$ as assumptions, one can derive the following two-dimensional nonlinear degenerate parabolic partial differential equation for $\rho$, where $\lambda \left(\rho \right)>0$ for $\rho >0$ and ∇ is the gradient operator:

$$\frac{{\partial }^{\alpha }\rho }{\partial t^{\alpha }}=\Delta \left(s\left(\rho \right)\right)+s\left(\rho \right),t>0,(\eta ,\xi )\in {\mathbb{R}}^2,$$

(22)

where the order of the fractional derivative with respect to time t is denoted by $\alpha$. The temporal fractional-order biological population model (TFBPM) is the name given to this equation. The analysis presents $\varphi \left(\rho \right)={\rho }^2$ as a specific instance for modeling animal populations, as explored by Gurney and Nisbet [48]. The migrations usually happen as a result of individuals traveling down the population density gradient, which moves more quickly at higher densities, in search of less congested areas to breed in a model that took into account an animal walking through a rectangular grid was created to mimic this behavior. With each step, the animal might either remain in the same spot or migrate toward the area with the lowest density. The size of the population density gradient at the relevant grid boundary determines the probability distribution for these movements. This model leads to

$$\frac{{\partial }^{\alpha }\rho }{\partial t^{\alpha }}=\Delta \left({\rho }^2\right)+s\left(\rho \right),t>0,(\eta ,\xi )\in {\mathbb{R}}^2,$$

(23)

with the initial condition $\rho (\eta ,\xi ,0)$ provided. When $\alpha =1$, this equation simplifies to the normal biological population model (NBPM):

$$\frac{{\partial }^{\alpha }\rho }{\partial t^{\alpha }}=\Delta \left({\rho }^2\right)+s\left(\rho \right),t>0,(\eta ,\xi )\in {\mathbb{R}}^2.$$

(24)

Additionally, various properties such as Hölder estimates and solutions of this model have been explored.

Constitutive equations for $s\left(\rho \right)$ may include the following:

• Malthusian Law: $s\left(\rho \right)=c\rho$, where c is a constant.
• Verhulst Law: $s\left(\rho \right)=c_1\rho -c_2{\rho }^2$, where $c_1,c_2$ are positive constants.
• Porous Media: $s\left(\rho \right)=c{\rho }^q$, where $c>0$ and $0<q<1$.

For a generalized form, consider

$$\frac{{\partial }^{\alpha }\rho }{\partial t^{\alpha }}=\frac{{\partial }^2{\rho }^2}{\partial {\eta }^2}+\frac{{\partial }^2{\rho }^2}{\partial {\xi }^2}+h{\rho }^a(1-\mathsf{\ell }{\rho }^b),t>0,(\eta ,\xi )\in {\mathbb{R}}^2,$$

(25)

subject to some appropriate initial conditions and where $h,a,\mathsf{\ell },b$ are real numbers. Under some specific parameter conditions, both the Verhulst and the Malthusian laws are covered by Equation (25).

Discretization Methodology

We consider the time-fractional partial differential equation given by Equation (25):

The $L1$ scheme is used to discretize the fractional time derivative $\frac{{\partial }^{\alpha }\rho }{\partial t^{\alpha }}$. This scheme is particularly suitable for fractional derivatives due to its capability to handle the memory effect inherent in such derivatives. The $L1$ approximation at a time $t_n$ is defined as: The $L1$ scheme approximates the fractional time derivative as follows:

$$\frac{{\partial }^{\alpha }{\rho }_{n+1}}{\partial t^{\alpha }}=\frac{{\left(\delta t\right)}^{-\alpha }}{\Gamma (2-\alpha )}\sum _{k=0}^n\left({\rho }_{n+1-k}-{\rho }_{n-k}\right)\left({(k+1)}^{1-\alpha }-k^{1-\alpha }\right)+O\left(\delta t^{2-\alpha }\right),$$

(26)

where $\delta t$ is the time step, and $\alpha$ is the fractional order of the derivative.

The correction term $B_n$, defined as

$$B_n=\frac{{\left(\delta t\right)}^{-\alpha }}{\Gamma (2-\alpha )}\sum _{k=1}^n\left({(k+1)}^{1-\alpha }-k^{1-\alpha }\right)({\rho }_{n+1-k}-{\rho }_{n-k}),$$

(27)

becomes zero when $n=0$. Therefore, the discretization formula simplifies as

$$\frac{{\partial }^{\alpha }{\rho }_{n+1}}{\partial t^{\alpha }}=\left\{\begin{array}{cc}\frac{{\left(\delta t\right)}^{-\alpha }}{\Gamma (2-\alpha )}({\rho }_{n+1}-{\rho }_n)+B_n+O\left(\delta t^{2-\alpha }\right),\hfill & \alpha \in (0,1),\hfill \\ \frac{{\rho }_{n+1}-{\rho }_n}{\delta t}+O\left(\delta t\right),\hfill & \alpha =1.\hfill \end{array}\right.$$

(28)

The Chebyshev spectral method involves representing the solution $\rho$ as a series expansion in terms of Chebyshev polynomials. These polynomials are particularly effective for approximating functions on bounded intervals due to their excellent approximation properties and the clustering of nodes at the endpoints, which can help in resolving boundary layer effects. The approximation of $\rho$ at a fixed time $t_n$ is given by

$$\rho (\eta ,\xi ,t_n)\approx \sum _{i=0}^N\sum _{j=0}^N{c}_{ij}^n{T}_i\left(\eta \right)T_j\left(\xi \right),$$

(29)

where $T_i$ are Chebyshev polynomials and $c_{ij}^n$ are the spectral coefficients.

The derivatives in the spatial domain are computed using the derivative properties of Chebyshev polynomials:

$$\frac{{\partial }^2{\rho }^2}{\partial {\eta }^2}\approx \sum _{i=0}^N\sum _{j=0}^N{\left(c_{ij}^n\right)}^2\frac{d^2}{d{\eta }^2}{T}_i\left(\eta \right)T_j\left(\xi \right),$$

(30)

$$\frac{{\partial }^2{\rho }^2}{\partial {\xi }^2}\approx \sum _{i=0}^N\sum _{j=0}^N{\left(c_{ij}^n\right)}^2{T}_i\left(\eta \right)\frac{d^2}{d{\xi }^2}{T}_j\left(\xi \right).$$

(31)

Combining the discretized forms of the time and space derivatives, we obtain the fully discretized version of the the model Equation (25), given by

$$\begin{array}{cc}\hfill \frac{1}{\Gamma (2-\alpha )}\sum _{k=0}^n\frac{{\rho }^{k+1}-{\rho }^k}{\Delta t^{\alpha }}& {((k+1)\Delta t)}^{1-\alpha }-{(k\Delta t)}^{1-\alpha }\hfill \\ \hfill & =\sum _{i=0}^N\sum _{j=0}^N\left[{\left(c_{ij}^n\right)}^2\left(\frac{d^2}{d{\eta }^2}{T}_i\left(\eta \right)+\frac{d^2}{d{\xi }^2}{T}_j\left(\xi \right)\right)\right]\hfill \\ \hfill & +h{\left(\sum _{i=0}^N\sum _{j=0}^N{c}_{ij}^n{T}_i\left(\eta \right)T_j\left(\xi \right)\right)}^a\left(1-\mathsf{\ell }{\left(\sum _{i=0}^N\sum _{j=0}^N{c}_{ij}^n{T}_i\left(\eta \right)T_j\left(\xi \right)\right)}^b\right).\hfill \end{array}$$

(32)

4. Error Analysis

This section deals with the error analysis of our proposed numerical scheme. Before the main results, we state some useful results in the form of lemmas [49].

Lemma 1

(Estimate for the Truncation Error). Consider the truncation error of a function u when approximated by its truncated Chebyshev series, denoted as $P_{N}u$, where

$$P_{N}u=\sum _{k=0}^N{\widehat{u}}_k{T}_k,$$

and ${\widehat{u}}_k$ represents the Chebyshev-series coefficients of u. The truncation error $u-P_{N}u$, measured in the weighted $L^2$ norm over the interval $(-1,1)$, satisfies the inequality

$$\parallel u-P_N{u\parallel }_{L_w^2(-1,1)}\le C{N}^{-m}{\parallel u\parallel }_{{\mathcal{H}}_{m;N}^w(-1,1)},$$

(33)

for all functions u belonging to the weighted Sobolev space ${\mathcal{H}}_m^w(-1,1)$ with $m\ge 0$. Here, C is a constant that depends on N and m, and ${\mathcal{H}}_{m;N}^w(-1,1)$ represents the Sobolev space characterized by integrating the function and its derivatives up to order m, each weighted by the Chebyshev weight.

Lemma 2

(Interpolation Error Estimate). The interpolation error associated with approximating a function u by its interpolant $I_{N}u$, which is defined at Chebyshev Gauss points across three different families. For this interpolant, which belongs to the polynomial space $P_N$, an important error estimates. The error between the function u and its interpolant, measured in the weighted $L^2$ norm over the interval $[-1,1]$, adheres to the following bound:

$$\parallel u-I_N{u\parallel }_{L_w^2(-1,1)}\le C{N}^{-m}{\parallel u\parallel }_{H_{m;N}^w(-1,1)},$$

(34)

$$\parallel u-I_N{u\parallel }_{L_w^{\infty }(-1,1)}\le C{N}^{\frac{1}{2}-m}{\parallel u\parallel }_{H_{m;N}^w(-1,1)},$$

(35)

where C represents a constant and $m\ge 1$ reflects the smoothness level of the function. This estimate is valid under the condition that u is an element of the weighted Sobolev space $H_m^w(-1,1)$. This space is characterized by considering functions that maintain their derivatives up to the m-th order, each weighted appropriately over the interval.

Lemma 3

(Integration Error Estimate). Given a function u from the weighted Sobolev space $H_m^w(-1,1)$ with $m\ge 1$ and a polynomial φ from the space $P_N$, the error produced by applying a Gauss-type quadrature formula for integration relative to the Chebyshev weight can be estimated as follows. Consider the integral of the product of u and φ, weighted by the Chebyshev weight $w\left(x\right)$, and its approximation by the quadrature formula:

$$\left|{\int }_{-1}^{1}u\left(x\right)\phi \left(x\right)w\left(x\right)dx-{(u,\phi )}_N\right|\le C{N}^{-m}{\parallel u\parallel }_{H_{m;N}^w(-1,1)}{\parallel \phi \parallel }_{L_w^2(-1,1)},$$

(36)

where

• ${(u,\phi )}_N$ denotes the approximation of the integral by the Gauss-type quadrature.
• C is a constant dependent on N and m, indicating the rate at which the error diminishes as the polynomial degree N or the smoothness m of the function u increases.
• ${\parallel u\parallel }_{H_{m;N}^w(-1,1)}$ and ${\parallel \phi \parallel }_{L_w^2(-1,1)}$ are the norms measuring the magnitude of u and φ in their respective function spaces.

One can state and prove the following error estimates for the fractional-order time discretization.

Theorem 1

([31]). Let ρ be the exact solution of the model Equation (25), and let ${\left\{{\rho }_k\right\}}_{k=0}^K$ be the time-discrete solution of with the initial condition ${\rho }_0\left(x\right)=\rho (x,0)$; then, the following error estimates hold:

1.

For $0\le a<1$,

$$\parallel \rho \left(t_k\right)-{\rho }_k{\parallel }_1\le c_{u,a}{T}^a\Delta t^{2-a},k=1,2,\dots ,K,$$

(37)

where $c_{u,a}:=\frac{c_u}{1-a}$ with $c_u$ a constant defined in (3.3).

2.

As $a\to 1$,

$$\parallel \rho \left(t_k\right)-{\rho }_k{\parallel }_1\le c_{u}T\Delta t,k=1,2,\dots ,K.$$

(38)

Similarly for the space discretization, using Lemmas 1–3, one can state and prove the following results of exponential order of convergence.

Theorem 2.

Let ${\left\{{\rho }_k^N\right\}}_{k=0}^K$ be the solution of the problem given in Equation (25) with the initial condition ${\rho }_0^N$ taken to be $I_N\rho \left(0\right)$, and ${\left\{{\rho }_k\right\}}_{k=0}^K$ be the numerical solution. Suppose ${\rho }_k\in H^m\left(K\right)\cap H_1^0\left(K\right)$, $m>1$; then,

• For $0\le \alpha <1$,

  $$\parallel {\rho }_k-{\rho }_k^N{\parallel }_{L_2}\le c_{\alpha }\Delta t^{\alpha }{N}^{-m}\underset{0\le j\le k}{max}{\parallel {\rho }_j\parallel }_m,k=1,2,\dots ,K,$$

  (39)
  where $c_{\alpha }=\frac{c}{1-\alpha }$ with c a constant.
• As $\alpha \to 1$,

  $$\parallel {\rho }_k-{\rho }_k^N{\parallel }_{L_2}\le c\Delta t^{-1}{N}^{-m}\underset{0\le j\le k}{max}{\parallel {\rho }_j\parallel }_m,k=1,2,\dots ,K,$$

  (40)
  where c only depends on T.

Proof. 

Assume $r^N={\rho }_{k+1}^N-w^N$ where $w^N$ is any function in $P_0^N\left(K\right)$. By direct calculation,

$$A^N(r^N,r^N)=A({\rho }_{k+1}-w^N,r^N)+A(w^N,r^N)-A^N(w^N,r^N)+F^N\left(r^N\right)-F\left(r^N\right),$$

(41)

which leads to

$$\begin{array}{cc}\hfill \parallel r^N{\parallel }_{1,N}^2& \le \parallel {\rho }_{k+1}-w^N{\parallel }_1{\parallel r^N\parallel }_1+\left|A(w^N,r^N)-A^N(w^N,r^N)\right|\hfill \\ \hfill & +\left|F\left(r^N\right)-F^N\left(r^N\right)\right|\forall w^N\in P_0^{N-1}\left(K\right).\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill \left|F\left(r^N\right)-F^N\left(r^N\right)\right|& =\left|(1-b_1)\left[({\rho }_k,r^N)-{({\rho }_k^N,r^N)}^N\right]\right.\hfill \\ \hfill & +\sum _{j=1}^{k-1}(b_j-b_{j+1})\left[({\rho }_{k-j},r^N)-{({\rho }_{k-j}^N,r^N)}^N\right]\hfill \\ \hfill & +\left.b_k\left[({\rho }_0,r^N)-{({\rho }_0^N,r^N)}^N\right]\right|.\hfill \end{array}$$

(43)

$$(g,r^N)-{(g^N,r^N)}^N\le (c{N}^{-m}{\parallel g\parallel }_m+\parallel g-g^N{\parallel }_{L_2})\parallel r^N{\parallel }_{L_2}.$$

(44)

Applying this to $g={\rho }_j$; $g^N={\rho }_j^N$ for all $j=1,2,\dots ,k$, we derive

$$\begin{array}{cc}\hfill \left|F\left(r^N\right)-F^N\left(r^N\right)\right|& \le (1-b_1)\parallel e_k^N{\parallel }_{L_2}+\sum _{j=1}^{k-1}(b_j-b_{j+1}){\parallel e_{k-j}^N\parallel }_{L_2}\hfill \\ \hfill & +b_k\parallel e_0^N{\parallel }_{L_2}+c{N}^{-m}\underset{0\le j\le k}{max}\parallel {\rho }_j{\parallel }_m{\parallel r^N\parallel }_{L_2}.\hfill \end{array}$$

(45)

This shows that the order of convergence of our combined numerical scheme is $O(\delta t^{2-\alpha }+N^{-m})$ in the $L_2$ norm in the case of spatial coordinates. Similarly, if we use the $L_{\infty }$ norm in the space coordinate, then the error estimate is of $O(\delta t^{2-\alpha }+N^{\frac{1}{2}-m})$.    □

5. Numerical Results

We present in this section a number of numerical examples to confirm our theoretical justification, as given in Theorems 1 and 2. In our simulations, we use a fixed time $t=0.1$ in comparison between the exact solution and the solution obtained by using the proposed scheme. In addition to this, in all of our simulations, we use a PC with processor 12th Gen Intel(R) Core(TM) i7-1255U and 16 GB RAM. The average CPU time in all these simulations is 0.03 s for collocation points from 15–50. The algorithm for these simulations is given in Appendix A.

Example 1.

Consider the model Equation (25) with $a=0,b=1$, that is, we have the following linear fractional-order-in-time biological population model of the form

$$\frac{{\partial }^{\alpha }\rho }{\partial t^{\alpha }}=\frac{{\partial }^2{\rho }^2}{\partial {\eta }^2}+\frac{{\partial }^2{\rho }^2}{\partial {\xi }^2}+h\rho ,$$

(46)

subject to the initial condition:

$${\rho }_0=\sqrt{\eta \xi }.$$

(47)

Equation (46), subject to the initial condition given in Equation (47), corresponds to Malthusian law. The exact solution is given by

$$\rho (\eta ,\xi ,t)=\sqrt{\eta \xi }{e}^{ht}{\mathcal{E}}_{\alpha }h{t}^{\alpha },$$

(48)

where ${\mathcal{E}}_{\alpha }h{t}^{\alpha }$ is the Mittag–Leffler function. When $\alpha \to 1$, then the exact solution becomes

$$\rho (\eta ,\xi ,t)=\sqrt{\eta \xi }\sum _{n=0}^{\infty }\frac{{\left(ht\right)}^n}{\Gamma (1+n)}=\sqrt{\eta \xi }{e}^{ht}.$$

(49)

In our simulations, we use the fractional-order $\alpha =0.5$. We compare our results for different collocation points with the exact solution as shown in Figure 1, Figure 2 and Figure 3. The error behavior between the exact and numerical solution is shown in Figure 4. One can see that the error between the exact and approximate solution decays exponentially by increasing the number of collocation points. These results emphasize the reliability of our numerical scheme, which is crucial for accurately predicting population dynamics in more detailed biological studies, making it practical for real-world biological population models where computational resources may be limited.

Example 2.

Consider the model Equation (25) with $a=1,b=1$. This gives us Verhulst law, and the model Equation (25) becomes

$$\frac{{\partial }^{\alpha }\rho }{\partial t^{\alpha }}=\frac{{\partial }^2{\rho }^2}{\partial {\eta }^2}+\frac{{\partial }^2{\rho }^2}{\partial {\xi }^2}+h\rho (1-\mathsf{\ell }\rho ),$$

(50)

subject to the initial condition

$${\rho }_0=e^{\sqrt{\frac{h\mathsf{\ell }}{8}}(\eta +\xi )},$$

(51)

The exact solution is given by

$${\rho }_0=e^{\sqrt{\frac{h\mathsf{\ell }}{8}}(\eta +\xi )}{\mathcal{E}}_{\alpha }h{t}^{\alpha }$$

(52)

For $\alpha \to 1$, the exact solution has the form

$${\rho }_0=e^{\sqrt{\frac{h\mathsf{\ell }}{8}}(\eta +\xi )+ht}.$$

(53)

The numerical experiments are performed between exact and approximate solution for $\alpha =0.5$ and different collocations points, as shown in Figure 5, Figure 6 and Figure 7. The error behavior between the exact and numerical solution is given in

Table 1. Example 2: Error behavior of exact vs. numerical solution.

N	${\mathit{L}}_{\mathit{\infty }}$	${\mathit{L}}_2$	N	${\mathit{L}}_{\mathit{\infty }}$	${\mathit{L}}_2$
6	$9.541\times {10}^{-1}$	$4.086\times {10}^{-1}$	8	$4.842\times {10}^{-2}$	$2.485\times {10}^{-2}$
10	$8.145\times {10}^{-3}$	$3.139\times {10}^{-3}$	12	$2.670\times {10}^{-4}$	$1.210\times {10}^{-4}$
14	$1.650\times {10}^{-5}$	$7.114\times {10}^{-6}$	16	$4.706\times {10}^{-7}$	$2.099\times {10}^{-7}$
18	$2.132\times {10}^{-8}$	$6.028\times {10}^{-9}$	20	$5.213\times {10}^{-10}$	$1.906\times {10}^{-10}$

. Again, one can see how fast the error is decreasing while increasing the number of collocation points. The numerical solution accurately captures the behavior of the exact solution, which represents the population growth with a carrying capacity, highlighting how the population stabilizes over time, which is a key aspect in ecological studies. The high-level accuracy is essential for making reliable predictions, especially in scenarios involving critical decisions about species conservation or resource management.

Example 3.

In this example, we choose $a=-1,b=1$. This gives us the model Equation (25) of the form

$$\frac{{\partial }^{\alpha }\rho }{\partial t^{\alpha }}=\frac{{\partial }^2{\rho }^2}{\partial {\eta }^2}+\frac{{\partial }^2{\rho }^2}{\partial {\xi }^2}+h({\rho }^{-1}-\mathsf{\ell }),$$

(54)

subject to the initial condition

$${\rho }_0=\sqrt{\frac{h\mathsf{\ell }}{4}{\eta }^2+\frac{h\mathsf{\ell }}{4}{\xi }^2+\xi +5},$$

(55)

The exact solution is given by

$$\rho (\eta ,\xi ,t)={\rho }_0+\sum _{n=0}^{\infty }\frac{\left(h{t}^{\alpha }\right){\rho }_0}{\Gamma (1+(n+1\left)\alpha \right)}+{\left(\frac{h{t}^{\alpha }}{{\rho }_0^2}\right)}^n$$

(56)

For $\alpha \to 1$, the transform exact solution is

$$\rho (\eta ,\xi ,t)={\rho }_0+\frac{ht}{{\rho }_0}{e}^{\frac{ht}{{\rho }_0^2}.}$$

(57)

Figure 8, Figure 9 and Figure 10 show the comparison between the actual and numerical solution for $\alpha =0.5$ with different collocations points. Our simulations confirmed that our numerical method performed very well against the exact solution. Further validation of exponential error decays is shown in

Table 2. Example 3: Error behavior of exact vs. numerical solution.

N	${\mathit{L}}_{\mathit{\infty }}$	${\mathit{L}}_2$	N	${\mathit{L}}_{\mathit{\infty }}$	${\mathit{L}}_2$
6	$1.901\times {10}^0$	$7.400\times {10}^{-1}$	8	$1.951\times {10}^{-1}$	$1.009\times {10}^{-1}$
10	$2.218\times {10}^{-2}$	$6.456\times {10}^{-3}$	12	$1.093\times {10}^{-3}$	$5.585\times {10}^{-4}$
14	$7.824\times {10}^{-5}$	$2.584\times {10}^{-5}$	16	$3.480\times {10}^{-6}$	$1.295\times {10}^{-6}$
18	$1.039\times {10}^{-7}$	$3.863\times {10}^{-8}$	20	$6.763\times {10}^{-9}$	$1.602\times {10}^{-9}$

. These results illustrate the model’s capability to handle nonlinear dynamics, which are common in biological systems where interactions between species or within populations can lead to complex behaviors.

6. Conclusions

This study has successfully shown how fractional calculus may be used to describe fractional-order-in-time partial differential equations, which are used to simulate the dynamic behavior of biological populations. It has been shown that using the $L1$ scheme for time-fractional derivatives is a reliable and accurate way to deal with the nonlocal properties included in these kinds of equations. This method is especially crucial for capturing the genetic and memory components that are essential to biological systems. Additionally, it has been demonstrated that the Chebyshev spectral method for spatial derivatives is a highly effective method that offers an exponential rate of convergence, making it perfect for managing intricate boundary conditions and geometries that are common in biological modeling. Our algorithms provide a powerful tool for approximation solutions in situations where analytical methods are not sufficient, as confirmed by our theoretical and numerical investigation. Our approach is accurate and reliable, as evidenced by the well-aligned numerical experiments with the theoretical predictions.