A Numerical Investigation of Fractional-Order Nonlinear Models in Population Dynamics

ABSTRACT

This study presents a numerical approach for solving fractional-order formulations of classical population dynamics models (the Logistic, Richards, Gompertz, and Predator-Prey models). The aim is to verify the impact of incorporating non-integer orders in these models. The methodology is based on a discretization scheme tailored to fractional differential equations, combined with Newton’s iterative method, for solving the resulting nonlinear systems. Stability, convergence, and robustness of the proposed method were verified by means of the mesh refinement and sensitivity tests. Comparisons with classical analytical solutions and high-precision numerical results also demonstrated the method’s accuracy. The results highlight the significant role of the fractional order α in modulating system behavior: smaller values of α lead to memory effects that slow population growth and attenuate oscillations in predator-prey interactions. These results indicate the potential of fractional calculus to improve the modeling of complex population dynamics.

Keywords:
fractional calculus; population dynamics; differential equations

1 INTRODUCTION

Mathematical modeling plays an important role in understanding population dynamics, as it enables the formulation of equations that describe growth, interaction, and decline processes in biological systems. Classical models, such as the logistic and predator-prey equations, have long been employed to capture the essential features of such dynamics. However, it is well known that these models are limited in their ability to incorporate memory and hereditary effects, which are intrinsic to many biological and ecological phenomena. In this regard, fractional calculus has emerged in recent years as a powerful tool for extending traditional models by introducing derivatives of non-integer order, thereby offering a more flexible and accurate framework to capture the complex dynamical behaviors observed in real-world populations.

In this context, numerous studies have explored the application of fractional differential operators to classical population models. The stability, existence, uniqueness, and numerical solution of the fractional-order logistic equation were investigated using an Adams-type predictor-corrector method ⁹. Mathematical models for mosquito-borne diseases involving systems of fractional-order differential equations were studied in ²². The Caleman linearization technique was applied to seek an exact solution to the fractional logistic equation with a Caputo derivative ²⁹. However, it was later demonstrated that the real function proposed by ²⁹ is not, in fact, an exact solution of the fractional logistic equation ⁴. The modeling of hepatitis B using systems of fractional differential equations was investigated in ⁸. Two systems of fractional differential equations were employed to study the dynamics between tumor cells and the immune system ²³. An adaptation of the three-step Adams-Bashforth scheme was used to numerically solve fractional population dynamics models, including the logistic and predator-prey equations ¹². A version of the fractional logistic equation involving the Caputo-Fabrizio derivative was solved analytically in ¹⁶. The Laplace Residual Power Series Method (LRPSM) was applied to numerically solve a fuzzy population dynamics model with a Caputo fractional derivative ¹⁷. Solutions to a population model involving a fractional derivative were obtained through a scheme combining the Laplace transform and the generalized power series method ². An approach based on the Fibonacci wavelet method and the quasi-linearization technique was proposed to solve the fractional-order logistic equation ¹. A multilayer neural network for deep learning based on fractional differential equations was employed to study two fractional logistic models, and parallel computing was used to determine the optimal network architecture ²⁸. Power series expansions were used to solve a fractional population model incorporating an Allee effect ⁵. A fractional differential equation was employed to model global population growth, and comparisons with real data from 1910 to 2010 showed that the fractional model outperformed its classical integer-order counterpart ²⁷. A nonlinear system of fractional differential equations was employed to study hepatitis B, with model parameters estimated from real-world data ⁷. A technique based on the Elzaki transform was applied to solve a fractional population model ¹⁰. A fractional population model with carrying capacity was analyzed in ²⁴. A model involving a Caputo fractional derivative was employed to describe the transmission dynamics of HIV/AIDS ³. Finally, the dynamics of COVID-19 transmission were modeled through a system of equations with Caputo fractional derivatives ¹³.

Unlike the aforementioned studies, the present work focuses on the application of Grünwald-Letnikov fractional derivatives to four classical nonlinear models: three based on a single ordinary differential equation (the Logistic, Richards, and Gompertz models) and one based on a system of ordinary differential equations (the Predator-Prey model). It also addresses the numerical solution of the resulting nonlinear systems, with the objective of evaluating the potential of fractional models to describe real biological systems and of highlighting the precautions that must be taken when applying such models.

2 FRACTIONAL-ORDER DERIVATIVES

Several definitions of fractional-order derivatives are available in the literature ^15),(18),(21. In this paper, however, we exclusively employ the left-sided Grünwald-Letnikov derivative, which can be approximated as follows ¹⁸:

$\frac{d^{\alpha }y}{d{t}^{\alpha }}\approx \frac{1}{h^{\alpha }}\sum _{j=0}^{t/h}{w}_j^{\alpha }y(t-jh),$ (2.1)

where t denotes time, y = y(t) is the dependent variable, α is the order of the derivative, h is the time step, and $w_j^{\alpha }$ are the weights defined recursively by

$w_0^{\alpha }=1,w_j^{\alpha }=\left(1-\frac{\alpha +1}{j}\right)w_{j-1}^{\alpha },j=1,2,\dots .$ (2.2)

3 POPULATION DYNAMICS MODELS

3.1 Logistic Model

The classical logistic model, also known as the Verhulst model ²⁵, is given by

$\frac{dP}{dt}=r_{p}P\left(1-\frac{P}{M}\right),$ (3.1)

where P = P(t) is the population size at time t, r _p is the intrinsic growth rate, and M is the carrying capacity of the environment, representing the maximum sustainable population size. The fractional logistic model, in turn, is given by

$\frac{d^{\alpha }P}{d{t}^{\alpha }}=r_{p}P\left(1-\frac{P}{M}\right),$ (3.2)

which can be rewritten

$\frac{d^{\alpha }P}{d{t}^{\alpha }}=r_{p}P-\frac{r_p}{M}{P}^2.$ (3.3)

3.2 Richards Model

The Richards model ²⁰, also known as the generalized logistic model, is given by

$\frac{dP}{dt}=r_{p}P\left[1-{\left(\frac{P}{M}\right)}^m\right],$ (3.4)

where m is a shape parameter that controls the asymmetry of the growth curve. The fractional form of the Richards model is given by

$\frac{d^{\alpha }P}{d{t}^{\alpha }}=r_{p}P\left[1-{\left(\frac{P}{M}\right)}^m\right],$ (3.5)

which can be rewritten as

$\frac{d^{\alpha }P}{d{t}^{\alpha }}=r_{p}P-\frac{r_p}{M^m}{P}^{m+1}.$ (3.6)

3.3 Gompertz Model

The Gompertz model ¹¹ is given by

$\frac{dP}{dt}=r_{p}P\ln \left(\frac{M}{P}\right),$ (3.7)

and its fractional counterpart is

$\frac{d^{\alpha }P}{d{t}^{\alpha }}=r_{p}P\ln \left(\frac{M}{P}\right).$ (3.8)

3.4 Predator-Prey Model

The classical predator-prey model, based on the Lotka-Volterra equations ^{14), (26}, is defined as the system:

$\left\{\begin{array}{l}\frac{dP}{dt}={\beta }_{1p}P-{\beta }_{2p}PZ,\\ \frac{dZ}{dt}=-{\beta }_{1z}Z+{\beta }_{2z}PZ,\end{array}\right.$ (3.9)

where P(t) is the prey population and Z(t) is the predator population. The parameters are: β _1p (intrinsic growth rate of the prey), β _2p (predation rate coefficient), β _1z (natural death rate of the predator), and β _2z (growth rate of the predator per prey consumed).

The corresponding fractional model is given by:

$\left\{\begin{array}{l}\frac{d^{\alpha }P}{d{t}^{\alpha }}={\beta }_{1p}P-{\beta }_{2p}PZ\\ \frac{d^{\alpha }Z}{d{t}^{\alpha }}=-{\beta }_{1z}Z+{\beta }_{2z}PZ,\end{array}\right.$ (3.10)

4 DIMENSIONAL ANALYSIS AND DISCRETIZATION

The application of fractional derivatives to the Logistic, Richards, Gompertz, and Predator-Prey models requires special attention to preserve dimensional consistency. In the classical setting, the operator d/dt has units of time⁻¹ , whereas the fractional operator d ^α /dt ^α has units of time^−α . This dimensional inconsistency is a well-known issue in the fractional calculus literature and has been discussed extensively in, among others, ^{8), (18), (22), (23}.

Although this inconsistency could be corrected mathematically by redefining the model parameters, such as growth and decay rates, this approach would alter their physical interpretation. A common alternative is to introduce a dimensional parameter τ with units of time. Since (1/τ ^{1 - α} )(d ^α /dt ^α ) has units of time⁻¹ , dimensional consistency is restored ^{8), (22), (23}. Equivalently, this corresponds to multiplying the right-hand side of each fractional equation by the factor τ ^{1 - α} .

This approach is adopted throughout the present work. Consequently, the dimensionally consistent coefficients are written as r = τ ^{1 - α} r _p for the Logistic, Richards, and Gompertz models. For the Predator-Prey model, the coefficients become β ₁ = τ ^{1 - α} β _1p , β ₂ = τ ^{1 - α} β _2p , β ₃ = τ ^{1 - α} β _1z , and β ₄ = τ ^{1 - α} β _2z .

The resulting fractional-order models and their corresponding discretized systems are presented below.

4.1 Fractional Logistic Model

The resulting fractional-order logistic model is given by

$\frac{d^{\alpha }P}{d{t}^{\alpha }}=rP-\frac{r}{M}{P}^2.$ (4.1)

4.1.1 Discretization of the Logistic Model

Considering a temporal mesh defined by t _i = t ₀ + ih, with h = (t _N −t ₀)/N for i = 0, 1, . . . , N, and using Eq. (2.1), Eq. (4.1) leads to the following nonlinear system:

$\sum _{j=0}^i{w}_j^{\alpha }{P}_{i-j}-h^{\alpha }r{P}_i+\frac{h^{\alpha }r}{M}{P}_i^2=0,$ (4.2)

where i = 1, 2, . . . , N, and P ₀ is the initial population size.

4.2 Fractional Richards Model

The resulting fractional-order Richards model is given by

$\frac{d^{\alpha }P}{d{t}^{\alpha }}=rP-\frac{r}{M^m}{P}^{m+1}.$ (4.3)

4.2.1 Discretization of the Richards Model

Using the same temporal mesh as in the logistic model and applying Eq. (2.1) to Eq. (4.3), we obtain the following nonlinear system:

$\sum _{j=0}^i{w}_j^{\alpha }{P}_{i-j}-h^{\alpha }r{P}_i+\frac{h^{\alpha }r}{M^m}{P}_i^{m+1}=0,$ (4.4)

where i = 1, 2, . . . , N, and P ₀ is the initial population size.

4.3 Fractional Gompertz Model

The resulting fractional-order Gompertz model is given by

$\frac{d^{\alpha }P}{d{t}^{\alpha }}=rP\ln \left(\frac{M}{P}\right).$ (4.5)

4.3.1 Discretization of the Gompertz Model

Using the temporal mesh t _i = t ₀ + ih and Eq. (2.1), the discretized form of Eq. (4.5) results in:

$\sum _{j=0}^i{w}_j^{\alpha }{P}_{i-j}-h^{\alpha }r{P}_i\ln \left(\frac{M}{P_i}\right)=0,$ (4.6)

where i = 1, 2, . . . , N, and P ₀ is the initial population size.

4.4 Fractional Predator-Prey Model

The resulting fractional-order predator-prey system is given by

$\left\{\begin{array}{l}\frac{d^{\alpha }P}{d{t}^{\alpha }}={\beta }_{1}P-{\beta }_{2}PZ,\\ \frac{d^{\alpha }Z}{d{t}^{\alpha }}=-{\beta }_{3}Z+{\beta }_{4}PZ,\end{array}\right.$ (4.7)

4.4.1 Discretization of the Predator-Prey Model

Using the same temporal mesh and applying Eq. (2.1), the discretized version of the system (4.7) is:

$\sum _{j=0}^i{w}_j^{\alpha }{P}_{i-j}-h^{\alpha }{\beta }_1{P}_i+h^{\alpha }{\beta }_2{P}_i{Z}_i=0,$ (4.8)

$\sum _{j=0}^i{w}_j^{\alpha }{Z}_{i-j}+h^{\alpha }{\beta }_3{Z}_i-h^{\alpha }{\beta }_4{P}_i{Z}_i=0,$ (4.9)

where i = 1, 2, . . . , N, P ₀ is the initial prey population, and Z ₀ is the initial predator population.

5 SOLUTION OF NONLINEAR SYSTEMS

The numerical scheme resulting from the discretization of the fractional differential equations was solved using Newton’s method. This section provides a detailed description of the method and presents numerical tests conducted to assess its convergence, stability, and sensitivity, as well as the effects of time-step refinement.

5.1 Newton’s Method for Solving Nonlinear Systems

Consider a nonlinear system of equations of the form

$\mathbf{F}\left(\mathbf{X}\right)=\mathbf{0},$ (5.1)

where

$\mathbf{F}\left(\mathbf{X}\right)=\left[\begin{array}{c}{f}_1(x_1,\dots ,x_n)\\ f_2(x_1,\dots ,x_n)\\ ⋮\\ f_n(x_1,\dots ,x_n)\end{array}\right],\mathbf{X}=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right],$ (5.2)

and each f _i is a nonlinear function of n variables. Newton’s method is an iterative procedure used to approximate the solution of such systems. The method proceeds as follows ⁶:

• 1. Choose an initial approximation X ⁽⁰⁾ ∈ ℝⁿ .

• 2. Compute the Jacobian matrix J(X) of partial derivatives:

$J\left(\mathbf{X}\right)=\left[\begin{array}{ccc}\frac{\partial f_1}{\partial x_1}& \cdots & \frac{\partial f_1}{\partial x_n}\\ ⋮& \ddots & ⋮\\ \frac{\partial f_n}{\partial x_1}& \cdots & \frac{\partial f_n}{\partial x_n}\end{array}\right].$ (5.3)

• 3. For each iteration k = 0, 1, 2, . . . , solve the linear system

$J\left({\mathbf{X}}^{\left(k\right)}\right)\Delta {\mathbf{X}}^{\left(k\right)}=-\mathbf{F}\left({\mathbf{X}}^{\left(k\right)}\right),$ (5.4)

• and update the approximation:

${\mathbf{X}}^{(k+1)}={\mathbf{X}}^{\left(k\right)}+\Delta {\mathbf{X}}^{\left(k\right)}.$ (5.5)

• 4. Repeat the process until a convergence criterion is satisfied, such as:

$\Vert \Delta {\mathbf{X}}^{\left(k\right)}{\Vert }_{\infty }<\epsilon \text{or}\Vert \mathbf{F}\left({\mathbf{X}}^{\left(k\right)}\right){\Vert }_{\infty }<\epsilon ,$ (5.6)

• for a prescribed tolerance ε > 0.

5.2 Convergence and Stability Tests

To verify the mesh independence of the numerical solutions obtained for the fractional models, simulations were performed using different values of N. For the Logistic, Richards, and Gompertz models, N = 20, 40, 80, 160, 320, 640, and 1280, with t ₀ = 0, t _N = 10, α = 0.95, ε = 10⁻⁵ , P ₀ = 10, r = 0.8, m = 4, and M = 500. All numerical experiments and simulations presented in this work were conducted with τ = 1.

Figure 1 shows the results for the fractional Logistic model. As N increases, the solution curves become increasingly similar. The maximum relative error at t = 10 between the solutions for N = 640 and N = 1280 was approximately 0.1%, justifying the use of N = 640 in the subsequent simulations.

Figure 1:
Mesh independence test - Logistic model.

A similar pattern is observed for the Richards model (Figure 2), with a maximum relative error of 0.003% at t = 10, reinforcing the adequacy of N = 640.

Figure 2:
Mesh independence test - Richards model.

For the Gompertz model (Figure 3), the maximum relative error at t = 10 between the two finest meshes was 0.001%, and once again, N = 640 was adopted.

Figure 3:
Mesh independence test - Gompertz model.

For the Predator-Prey model, larger mesh sizes (N = 40 to 2560) were used due to the greater number of equations involved (2N). Simulations were performed with t ₀ = 0, t _N = 100, P ₀ = 40, Z ₀ = 9, α = 0.95, ε = 10⁻⁵ , β ₁ = 0.1, β ₂ = 0.02, β ₃ = 0.01, and β ₄ = 0.1. Figures 4 and 5 show the solution curves for P (prey) and Z (predator), respectively. The mesh with N = 1280 was selected, since the relative error between N = 1280 and N = 2560 was 0.8% for P and 3% for Z.

Figure 4:
Mesh independence test - prey population (P) in the Predator-Prey model.

Figure 5:
Mesh independence test - predator population (Z) in the Predator-Prey model.

As expected, reducing the time step h (i.e., increasing N) led to smaller differences in the solution curves, suggesting convergence of the method. To further assess this, sensitivity tests were conducted by introducing 1% perturbations in the initial conditions and, subsequently, in both the initial conditions and parameters. The objective was to quantify the impact of small variations on the solutions.

Table 1 summarizes the maximum percentage deviations observed in P or Z. With the exception of the Predator-Prey model, 1% perturbations in initial conditions led to deviations of the same magnitude. When perturbations were also applied to parameters, variations generally ranged from 1.4% to 4.4%, with the Predator-Prey model showing higher sensitivity (up to 7.6%).

Table 1:
Maximum pointwise percentage deviations due to perturbations in the initial conditions (ICs) or both ICs and model parameters.

Model	ICs only	ICs and parameters
Logistic	1.0%	3.3%
Richards	1.0%	4.4%
Gompertz	1.0%	2.5%
Predator-Prey (P)	2.5%	7.6%
Predator-Prey (Z)	2.9%	6.6%

The analytical solutions for the classical models (with integer-order derivatives) are given by:

Logistic model:

$P\left(t\right)=\frac{M}{1+{\displaystyle \frac{M-P_0}{P_0}}{e}^{-rt}},$ (5.7)

Richards model:

$P\left(t\right)=\frac{M}{{\left[1+\left({\left({\displaystyle \frac{M}{P_0}}\right)}^m-1\right)e^{-rmt}\right]}^{1/m}},$ (5.8)

Gompertz model:

$P\left(t\right)=M{e}^{-e^{-rt}\ln \left(\frac{M}{P_0}\right)},$ (5.9)

For α ≈ 1, the solutions of the fractional models are expected to approximate the classical solutions above. Figure 6 shows the relative errors over time between the classical analytical solutions and the numerical solutions of the fractional models with α = 0.9999. For the Predator-Prey model, the reference solution was computed using the classical Runge-Kutta method of order 4. In all cases, the maximum relative error remained below 2%, highlighting the high accuracy of the method.

Figure 6:
Relative error between the analytical (or high-precision numerical) solutions of the classical models and the fractional model solutions with α = 0.9999.

These tests demonstrate the robustness and reliability of the numerical method across different fractional population models. The consistent convergence with mesh refinement, along with sensitivity analyses, confirms that the chosen discretization ensures both stability and accuracy. Furthermore, the agreement between the solutions of the fractional models for α ≈ 1 and the corresponding classical solutions provides strong validation of the proposed approach.

6 RESULTS

This section presents the results obtained from the numerical solutions of the fractional models. In all cases, the same parameter values used in the previously described tests were employed. For the Logistic, Richards, and Gompertz models, results were obtained for t = 0 to t = 10 with α values of 1.00, 0.99, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, and 0.50. For the Predator-Prey model, results were obtained for t = 0 to t = 50, with α values of 1.00, 0.98, 0.96, 0.94, 0.92 and 0.90.

Figure 7 displays the results for the fractional logistic growth model. The solution curves exhibit characteristic sigmoidal growth behavior. The population initially grows slowly, then accelerates during an intermediate phase (approximately between t = 2 and t = 8, depending on the value of α), and finally the growth rate decreases as the population approaches its carrying capacity (M = 500).

Figure 7:
Results - Logistic model.

Higher values of α (e.g., 1.00, 0.99, and 0.95) result in significantly faster population growth. The population reaches the inflection point (where the growth rate is maximal and corresponds to half the carrying capacity) and approaches saturation near t = 10. For smaller values of α (e.g., 0.50, 0.55, and 0.60), the growth is substantially slower. The population reaches only between 200 and 300 at t = 10, remaining in a pronounced growth phase and still far from saturation.

In this context, the order of the derivative serves as a memory-related parameter that indirectly governs the system’s responsiveness to growth. Higher values of α indicate that the population grows more rapidly under ideal conditions (e.g., in scenarios with abundant initial resources).

Lower values of α, in contrast, result in slower initial growth, which can be particularly useful for modeling situations with limited early resources or representing populations that face initial adaptation challenges and progressively overcome them over time.

Figure 8 presents the results of the fractional Richards model. As observed in the logistic model, all curves exhibit a sigmoidal shape. The population initially grows slowly, accelerates during an intermediate phase, and then the growth rate decreases as it approaches M = 500.

Figure 8:
Results - Richards model.

The most notable feature here, compared to the previous plot, is that the main effect of varying α is a horizontal shift of the curves along the time axis. The overall shape of the curve and the saturation level remain virtually unchanged. However, the order of the derivative influences the timing of the accelerated growth phase and, consequently, the time required to reach saturation.

Although the timing is shifted, the slopes of the steepest central regions of each curve remain nearly identical. This suggests the existence of a threshold-like temporal activation, after which the population begins to grow rapidly, and that the value of α determines how early (or late) this threshold is reached.

From an application perspective, in situations where a specific event or condition initiates growth, the fractional order α could represent the efficiency with which this hypothetical trigger operates-leading to either a faster or slower onset of population growth. While this interpretation is intuitive, it warrants further investigation and may be explored more rigorously in future studies.

Figure 9 presents the results obtained for the fractional Gompertz model. For all values of α, the population initially grows at an increasing rate, and then the growth rate decreases as the population approaches the carrying capacity M. As in previous models, higher values of α (e.g., 1.00, 0.99, and 0.95) result in faster initial population growth and quicker convergence to saturation. For instance, for α = 1.00, the population reaches approximately 450 around t = 5. In contrast, for α = 0.50, growth is significantly slower, and the population remains in an accelerated expansion phase at t = 10.

Figure 9:
Results - Gompertz model.

Interestingly, during the early stage (t = 0 to t = 2), the curves for different α values are very close and even intersect around t = 2 and a population size of approximately 200. This suggests that the order of the derivative has a less pronounced effect at the beginning of the process. However, this behavior may be attributed to the attenuation of the fractional memory effect during the initial stages.

Another noteworthy observation is that the plateau reached by the curve for α = 0.50 appears significantly lower than M within the observed time interval. This may suggest that α influences not only the growth rate but also the system’s long-term equilibrium or apparent carrying capacity.

This is clarified in Figure 10, which shows the solution curves extended to t = 200. For larger values of t, P continues to grow and gradually tends toward M, albeit more slowly, with P(50) = 447, P(100) = 463, P(150) = 470, and P(200) = 474.

Figure 10:
Results for t = 0 to t = 200 (α = 0.5) - Gompertz model.

Figure 11 presents the results of the Predator-Prey model, where the black dashed line represents the prey population and the red dotted line represents the predator population. The graphs exhibit cyclic oscillations in both populations, a characteristic behavior of such systems. Peaks in the prey population (P) generally precede those in the predator population (Z). This pattern aligns with biological reasoning: an increase in prey provides more food for predators, leading to a subsequent rise in the predator population. This, in turn, causes a decline in the prey population due to intensified predation, which is followed by a decline in the predator population due to food scarcity.

Figure 11:
Comparison between P (black dashed line) and Z (red dotted line) in the Predator-Prey model for different values of α.

For α = 1, the system exhibits high-amplitude oscillations in both populations. The prey population reaches peaks of approximately 40 and troughs near zero, while the predator population peaks around 20, also reaching values close to zero at the minima. For α = 0.98, the oscillations remain of high amplitude, though a slight reduction is noticeable. For smaller values of α, this decreasing trend continues, and the oscillations become progressively damped, leading the populations to fluctuate within a narrower range.

The plots suggest that the order of the derivative may play a significant role in the stability and coexistence dynamics of the Predator-Prey system. Higher values of α promote large-amplitude oscillations, which may pose a potential risk to the long-term persistence of both species due to population troughs approaching zero. In contrast, lower values of α (e.g., around 0.90) tend to guide the system toward a more stable state, in which oscillations are damped and the populations stabilize at a coexistence equilibrium.

7 CONCLUSION

This study presented a numerical framework for solving the Logistic, Richards, Gompertz, and Predator-Prey population models in their fractional-order formulations. The discretization scheme, coupled with Newton’s iterative method, proved to be effective in addressing the non-linearities inherent to the models under consideration. The mesh refinement and sensitivity analyses confirmed both the numerical stability and the convergence of the solutions. Furthermore, comparisons with classical analytical solutions and high-precision numerical results validated the method’s accuracy and reliability.

The fractional order α proved to be an important parameter influencing system dynamics. It introduces memory effects that modulate both growth rates and oscillatory behavior. Notably, lower values of α resulted in slower growth in single-species models and contributed to the attenuation of oscillations in the predator-prey dynamics. These effects highlight the enhanced flexibility of fractional-order models in capturing complex population behaviors that may not be adequately represented by classical integer-order formulations.

This study contributes theoretically by highlighting the feasibility and advantages of fractional modeling in representing the complex temporal dynamics of population systems. In real-world contexts, several studies have demonstrated the superiority of fractional models over classical ones. The study by ¹⁹, for instance, which modeled the spread of COVID-19, showed that the fractional model produced predictions that were closer to the observed data, a result similar to that reported by ¹³. Nevertheless, certain limitations are acknowledged, particularly regarding the selection of appropriate fractional orders in practical applications. Future research could address these challenges by incorporating data-driven parameter estimation and applying the proposed approach to real cases. Such efforts would enhance its practical relevance and predictive capability.

Practically, the fractional population models and numerical techniques presented herein offer promising tools for ecological studies in which memory and hereditary effects are significant. Possible applications include the modeling of species with delayed responses to environmental changes, epidemics characterized by incubation periods, and systems strongly impacted by the past.

In summary, the findings reinforce the potential of fractional calculus in population dynamics, enabling the representation of more realistic biological behaviors. The proposed numerical methodology offers a robust foundation for future studies seeking to explore fractional models in broader ecological or epidemiological contexts.

Acknowledgments

I acknowledge the Santa Catarina State University for providing the infrastructure required for the development of this research.

Data availability

Datasets related to this article are available upon request to the corresponding author.

REFERENCES

• ¹ S. Ahmed, S. Jahan & K.S. Nisar. Hybrid Fibonacci wavelet method to solve fractional-order logistic growth model. Mathematical Methods in the Applied Sciences, 46 (2023), 16218-16231.
• ² M. Alaroud, A.K. Alomari, N. Tahat & A. Ishak. Analytical Computational Scheme for Multivariate Nonlinear Time-Fractional Generalized Biological Population Model. Fractal and Fractional, 7 (2023).
• ³ S.A. Alharbi & N.A. Almuallem. Computational Study of a Fractional-Order HIV Epidemic Model with Latent Phase and Treatment. Fractal and Fractional , 9 (2025).
• ⁴ I. Area, J. Losada & J.J. Nieto. A note on the fractional logistic equation. Physica A, 444 (2016), 182-187.
• ⁵ I. Area & J.J. Nieto. On the fractional Allee logistic equation in the Caputo sense. Examples and Counterexamples, 4 (2023).
• ⁶ R.L. Burden & J.D. Faires. “Numerical Analysis”. Brooks/Cole, Cengage Learning, 9 ed. (2010).
• ⁷ L.C. Cardoso, , P.F. Toro, R.F. Camargo & F.L.P. Santos. A computational and fractional mathematical model for assessing hepatitis B with optimal therapy. Journal of Biological Systems, 32 (2024), 1203-1223.
• ⁸ L.C. Cardoso, F.L.P. Santos & R.F. Camargo. Analysis of fractional-order models for hepatitis B. Computational and Applied Mathematics, 37 (2018), 4570-4586.
• ⁹ A.M.A. El-Sayed, A.E.M. El-Mesiry & H.A.A. El-Saka. On the fractional-order logistic equation. Applied Mathematics Letters, 20 (2007), 817-823.
• ¹⁰ T.M. Elzaki & M.Z. Mohamed. A Novel Analytical Method for the Exact Solution of the Fractional-Order Biological Population Model. Acta Mechanica et Automatica, 18 (2024), 564-570.
• ¹¹ B. Gompertz. On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. Philosophical Transactions of the Royal Society of London, 115 (1825), 513-583.
• ¹² H. Jafari, R.M. Gangi, N.S. Nkomo & Y.P. Lv. A numerical study of fractional order population dynamics model. Results in Physics, 27 (2021).
• ¹³ M. Khalighi, L. Lahti, F. Ndäırou, P. Rashkov & D.F.M. Torres. Fractional modelling of COVID-19 transmission incorporating asymptomatic and super-spreader individuals. Mathematical Biosciences, 380 (2025).
• ¹⁴ A.J. Lotka. “Elements of Physical Biology”. Williams and Wilkins, Baltimore (1925).
• ¹⁵ K.S. Miller & B. Ross. “An Introduction to the Fractional Calculus and Fractional Differential Equations”. John Wiley and Sons, New York (1993).
• ¹⁶ J.J. Nieto. Solution of a fractional logistic ordinary differential equation. Applied Mathematics Letters , 123 (2022).
• ¹⁷ M.N. Oqielat, T. Eriqat, O. Ogilat, Z. Odibat, Z. Al-Zhour & I. Hashim. Approximate solutions of fuzzy fractional population dynamics model. European Physical Journal Plus, 137 (2022).
• ¹⁸ I. Podlubny. “Fractional Differential Equations”. Academic Press, San Diego (1999).
• ¹⁹ K. Rajagopal, N. Hasanzadeh, F. Parastesh, I.I. Hamarash, S. Jafari & I. Hussain. A fractional-order model for the novel coronavirus (COVID-19) outbreak. Nonlinear Dyn, 101 (2020).
• ²⁰ F.J. Richards. A flexible growth function for empirical use. Journal of Experimental Botany, 10(2) (1959), 290-301.
• ²¹ S.G. Samko, A.A. Kilbas & O.I. Marichev. “Fractional Integrals and Derivatives: Theory and Applications”. Gordon and Breach Science Publishers, Amsterdam (1993).
• ²² T. Sardar, S. Rana, S. Bhattacharya, K. Al-Khaled & J. Chattopadhyay. A generic model for a single strain mosquito-transmitted disease with memory onthe host and the vector. Mathematical Biosciences , 263 (2015), 18-36.
• ²³ J.G. Silva, A.C.O. Ribeiro, R.F. Camargo, P.F.A. Mancera & F.L.P. Santos. Stability analysis and numerical simulations via fractional calculus for tumor dormancy models. Commun Nonlinear Sci Numer Simulat, 72 (2019), 528-543.
• ²⁴ J. Singh, R. Agrawal & D. Baleanu. Dynamical analysis of fractional order biological population model with carrying capacity under Caputo-Katugampola memory. Alexandria Engineering Journal, 91 (2024), 394-402.
• ²⁵ P.F. Verhulst. Notice on the law that a population follows in its growth. Correspondance Mathématique et Physique, 10 (1838), 113-121.
• ²⁶ V. Volterra. Fluctuations in the abundance of a species considered mathematically. Nature, 118 (1926), 558-560.
• ²⁷ O.K. Wanassi & D.F.M. Torres. An integral boundary fractional model to the world population growth. Chaos, Solitons and Fractals, 168 (2023).
• ²⁸ J.L. Wei, G.C. Wu, B.Q. Liu & J.J. Nieto. An optimal neural network design for fractional deep learning of logistic growth. Neural Computing and Applications, 35 (2023), 10837-10846.
• ²⁹ B.J. West. Exact solution to fractional logistic equation. Physica A , 429 (2015), 103-108.

Associate editor:

Rosana Sueli Motta Jafelice

Publication Dates

• Publication in this collection
  10 Aug 2026
• Date of issue
  2026

History

• Received
  25 Aug 2025
• Accepted
  09 June 2026