Fractional partial differential equations (FPDEs) are a generalization of classical partial differential equations (PDEs) by incorporating fractional derivatives of arbitrary order, providing a versatile framework for describing a broad range of phenomena in physical, chemical, biological, and financial processes with anomalous transport mechanisms [1]. These equations capture non-local and non-linear aspects that are not captured by classical PDEs [2]. Real-world physical models often entail substantial uncertainty arising from numerous variables [3].

Solving FPDEs numerically poses challenges due to their inherently non-local and non-linear nature [2, 4, 5].

A fractional diffusion equation can be expressed in the form [1]

where t>0, L ≤ x ≤ R, and 0 < α ≤ 1, with initial condition given by Equation 2:

and boundary conditions: u(L, t) = f(t), and u(R, t) = g(t), where  0Dtαu(x,t) in Equation 1 is the Caputo fractional derivative of order α, and is given by Equation 3 as follows:

where m−1 < α < m, m∈ℕ [6].

In the realm of FPDEs, diffusion equation via the Caputo operator can be used in modeling transport of diffusion processes, by offering a profound mechanism for describing memory and hereditary properties of materials and processes. This aspect is crucial for our work with a fractional diffusion experiment, where the traditional integer-order models fall short in capturing the complex dynamics of anomalous diffusion prevalent in many natural and engineering systems. Incorporating the Caputo operator allows for a more accurate representation of the non-local and history-dependent behavior of diffusion processes, enhancing the model's ability to predict and analyze real-world scenarios effectively [7, 8]. Furthermore, the insights from Awadalla et al. [9] and Shaikh and Qureshi [10] shed some light on contributions of the Caputo operator toward modeling complex dynamical systems, including its application in fractional optimal control models and bifurcation analysis of human syncytial respiratory virus transmission dynamics. Awadalla et al. [9] showed the operator's efficacy in applications ranging from fractional optimal control models to the bifurcation analysis of human respiratory virus transmission dynamics, reinforcing the Caputo operator's indispensability in capturing some dynamics of such systems.

Mathematical modeling of pattern formation in coral reefs using fractional differential equations shows an elegant approach to capturing the complex dynamics of ecological systems [11]. Fractional models are particularly adept at describing phenomena with memory effects and spatial heterogeneity, characteristics inherent to ecological patterns observed in coral reefs. A fractional partial differential equation modeling coral reef growth and interaction can be expressed as [12]:

where N(x, t) represents the coral population density at location x and time t, D_N is the diffusion coefficient capturing the spread of coral larvae, α and 2β denote the fractional orders of time and space derivatives reflecting memory and spatial dispersal effects, R(N, C) is the growth rate dependent on coral and nutrient concentration C, and M(N) represents natural mortality. Such models allow for understanding of how local interactions and environmental conditions influence large-scale pattern formation, offering insights into conservation strategies and reef resilience [13, 14].

Recent advancements in numerical methods for solving non-linear partial differential equations (PDEs) encompass notable techniques such as the optimal homotopy continuation method [15], variational homotopy method [16], and a hybrid iterative approach [17]. Qureshi et al. [18] used a three-step numerical scheme showing ninth-order convergence in solving non-linear equations. Comparative analyses revealed its superior efficiency when compared to established schemes in computational science. In addition, the authors, in [19], investigated a novel optimal iterative algorithm with fourth-order accuracy tailored for root-finding in real functions.

Analytic solutions of most FPDEs cannot be obtained explicitly. In fact, the closed form solutions for many time-fractional differential equations are not available [1, 20]. Some nonstandard finite difference (NSFD) schemes were first devised by Mickens [21].

The primary merit of the NSFD schemes lies in their ability to overcome the inherent numerical instabilities associated with classical finite difference methods [21]. The development of NSFD methods adheres to some fundamental rules for its practical implementation [22–26]:

where ϕ(k,λ)=k+O(k2).

Recently, NSFD methods have been constructed by many researchers for solving differential equations. Agarwal and El-Sayed proposed a nonstandard finite difference method and a Chebyshev collocation method to solve a fractional order diffusion equation [27]. Momani et al. [25] constructed a nonstandard implicit Euler method for solving a class of fractional partial differential equations. Tijani and Appadu [28] constructed an unconditionally positive nonstandard finite difference scheme for a mathematical model of biofilm formation on a medical implant. The model employed uses the bistable Allen-Cahn partial differential equation, which is a generalization of Fisher's equation. Appadu and Tijani [29] obtained the numerical solution of a 1-D generalized Burgers-Huxley equation under specified initial and boundary conditions, and they used Forward Time Central Space (FTCS) and a nonstandard finite difference scheme. Agbavon and Appadu [26] constructed a nonstandard finite difference schemes to solve the FitzHugh-Nagumo equation with specified initial and boundary conditions under three different regimes. Kehinde et al. [30] obtained solution to a two-dimensional semilinear singularly perturbed semilinear convection-diffusion problem by constructing a nonstandard finite difference method. Jejeniwa et al. [31] used three methods: the Kowalic-Murty scheme, Lax-Wendroff scheme, and nonstandard finite difference (NSFD) scheme, to solve 1D and 2D convective diffusion equations. They looked at the cases when the advection velocity was much greater than of the diffusion coefficient and when the coefficient of diffusion was much greater than the advection velocity. They also analyzed the dispersion properties of the three methods.

Fractional subdiffusion and superdiffusion are described by the parameter α, which determines the anomalous nature of diffusion. In subdiffusion, α ranges from 0 to 1 (0 < α < 1), indicating that the mean squared displacement (MSD) increases slower than in classical diffusion. This behavior is modeled by the time-fractional diffusion equation:

where D is the diffusion coefficient and ∂^α/∂t^α represents a Caputo derivative, which is used to incorporate memory effects into the system [5]. This model is particularly relevant in biological and environmental contexts where diffusion is hindered by physical barriers and complex internal structures [32].

Conversely, for superdiffusion, the α values lie between 1 and 2 (1 < α < 2), where the MSD grows faster than it does in normal diffusion, suggestive of long jumps and persistent directional movement [33]. This dynamic is modeled by

where ∇^α represents the fractional Laplacian of order α. This equation is commonly applied in modeling dynamics driven by long-range interactions, such as in environmental science for the rapid spread of pollutants and in financial markets for capturing the large fluctuations typical of stock prices and market indices [34]. These applications are detailed in Anomalous Diffusion by Metzler and Klafter, which explores the relevance of these models in turbulent flows and financial markets [34]. See also [35] for fractional Kinetic equations and [36] for some physical applications of fractal operators.

This study is novel in many ways, and the study we carried out is quite different from that by Stynes et al. [37]. Stynes et al. [37] investigated the regularity of solutions for a general equation of the form

and obtained bounds for the L_∞ error. They then constructed an implicit classical finite difference scheme to solve three problems described by the following equations:

There is only one figure shown in the numerical profile for problem 3 which they considered. They calculated the numerical order of convergence using the implicit classical finite difference scheme for α = 0.4, 0.6, and 0.8 and showed that r=2-αα gives the optimal rate of convergence and that other values of α cause some deviation between the theoretical and numerical rates of convergence.

In this study, we have constructed four explicit finite difference methods to solve one of the three problems considered by Stynes et al. [37]. The first novelty is that such methods were not used before to solve that problem. Second, we compared standard and nonstandard finite difference methods using conformable and Caputo operators with regard to the study of stability and conditions for the following:

There are many concluding remarks on this study, as described below. We find that FTCSCO and NSFDCO give issues with the numerical rate of convergence. FTCSCA and NSFDCA are quite reliable methods to solve the problem we considered.

We show that the range of values of the time step size (at a given spatial step size of π64) for preservation of positivity is almost similar to that for stability for the NSFDCA scheme. The stability analyses of FTCSCA and NSFDCA are not very straightforward as the expressions for the amplification factors are quite complicated.

The organization of this study is briefly explained. In Section 2, the problem chosen is elaborated upon. Section 3 provides background information on conformable and Caputo fractional order derivatives. In Section 4, we derive the FTCSCO scheme, a finite difference scheme employing the conformable derivative. A comprehensive investigation into the stability and consistency analyses of the FTCSCO scheme is conducted. Section 5 introduces the NSFDCO scheme; a nonstandard finite difference scheme utilizing the conformable approximation. A thorough study on the positivity condition and consistency of the NSFDCO scheme is presented. Sections 6, 7 focus on construction of two novel finite difference schemes, incorporating FTCS and Caputo approximations (referred to as FTCSCA), and NSFD combined with Caputo approximations to derive a scheme known to be NSFDCA. In Section 8, we present results and tabulate the numerical rate of convergence in time. Conclusion is provided in Section 9.

The homogeneous time-fractional diffusion equation [5] is given by Equation 4 below:

where 0 < α ≤ 1, 0 ≤ x ≤ L, with boundary conditions given by Equation 5:

We note that K is the diffusion coefficient, L is the length of the domain, and α is the fractional order derivative.

In this study, we solve the non-homogeneous time-fractional diffusion equation [37] given by Equation 6:

where x∈[0, π], t∈(0, 1], with the initial condition and boundary conditions given by Equations 7 and 8 respectively:

We chose the spatial step size as Δx=π64. The spatial and temporal step sizes are denoted by h and k, respectively.

We note there is no known exact solution for this problem.

The numerical rate of convergence R_T is calculated as [28, 38, 39]

and the discrete maximum norm errors given by Equation 9:

Luchko [40] has proved the existence and uniqueness of a classical solution to the PDE

Stynes et al. [37] considered the problem in Equation 6, where the Caputo derivative is approximated by L₁ scheme and a classical finite difference operator is used to discretize ∂2u∂x2. They derived bounds on the L_∞ error and obtained the L_∞ errors when α = 0.2, 0.4, 0.6, and 0.8.

We should point out here that our numerical simulations were conducted using a Dell computer equipped with a Windows 11 operating system. The system specifications include an Intel Core i5 processor and 256 GB of RAM, with 8 GB of memory.

There exist several definitions of fractional derivatives of order α>0, with the most widely used being the Riemann-Liouville (RL), Caputo, and conformable fractional derivatives [cf. [5, 41–44]].

Definition 1. The Riemann-Liouville fractional integral is defined as follows by Equation 10:

where α>0 and Γ(α)=∫0∞xα-1e-x dx is the Euler Gamma function [cf. [1]].

Definition 2. The Caputo time-fractional derivative operator of order α>0 (m−1 < α ≤ m, m∈ℕ) for a real-valued function u(x, t) is defined as follows [6] by Equation 11:

Similarly, the Caputo space-fractional derivative operator  aDxαu(x,t) of order α>0 (m−1 < α ≤ m), m∈ℕ can be defined [43, 45]. It is worth noting that if u is sufficiently smooth [46], the fractional derivative Dtαu recovers the typical first-order derivative u′(t) as α → 1⁻ [45].

Definition 3. [41] For a function g:[0, ∞] → ℝ, the conformable fractional derivative of g of order α is defined by

Fundamental concepts and properties of conformable calculus are detailed in [41], and it is noteworthy that the conformable derivative is chosen to preserve some classical properties of standard calculus [42]. Given that other popular fractional derivatives such as Caputo and Riemann-Liouville lack certain natural properties of derivatives, including product rule, quotient rule, and chain rule, some authors, such as Khalil et al. [41] and Abdelwajad [42], motivate the study of the conformable derivative to fill these gaps and maintain some natural properties of derivatives [cf. [1, 45]].

Abdelwajad [42] approximated the time-fractional derivative using conformable approximation:

We next prove Equation 13. Using Equation 12, we have Equation 14.

Let k = η t^1−α. This gives Equation 15 and the steps involved are shown:

Hence, we can approximate ∂αu(x,t)∂tα by k1-α[Uin+1-Uink], if we choose to use a forward difference approximation for ∂u∂t in Equation 13.

For further study on conformable fractional derivatives, interested readers are referred to the following studies: [41, 42, 47, 48].

Mickens is the architect of nonstandard finite difference methods, and the rules for the construction of these methods are given in [21]. We construct a new scheme denoted as NSFDCO scheme. We first approximate the fractional derivative using conformable operators, and then, nonstandard finite difference approximations are used to approximate the derivatives. The NSFDCO scheme is given by

where ϕ(k) = e^k−1 and ψ(h) = 1−e^−h.

By multiplying Equation 24 by ϕ(k)·k^α−1, we obtain Equation 25 and Equation 26