Throughout the history of mathematical modeling, traditional calculus has been the primary tool for explaining the behavior of various phenomena, ranging from the physics of kinetics and electromagnetism to the spread of diseases in biology and numerous other systems in between. For a significant period, it was widely accepted that classical calculus sufficed to model these systems comprehensively. However, contemporary research suggests that a more generalized form of calculus, known as fractional calculus, offers superior modeling capabilities for these phenomena [1,2,3,4,5]. The concept that fractional calculus might provide more accurate models than traditional calculus is still under investigation. The immediate implications of these enhanced models include a more detailed understanding of the specifics behind systematic evolutions and anomalous behaviors in many physical problems. Often, such nuances are overlooked when using traditional calculus models. By comprehending fractional calculus, we can elucidate patterns in physical phenomena that traditional methods have not yet captured. This deeper insight enhances our understanding of these problems and may drive innovative advancements in their modeling techniques [6,7].

A reaction-diffusion model is a mathematical framework used to describe the behavior of chemical or biological systems in which substances interact (reaction) and spread out in space (diffusion). This type of model is instrumental in understanding a wide range of phenomena, from chemical reactions and biological pattern formation to ecological dynamics and neural activity. A typical reaction-diffusion system is governed by a set of partial differential equations (PDEs) of the form where ${\vartheta }_i={\vartheta }_i(x,t)$ represents the concentration of the $i-$th substance at position $x$ and time $t$, $d_i$ is the diffusion coefficient of the $i-$th substance, indicating how it spreads through space. $\mathsf{\Delta }$ denotes the Laplacian operator, representing spatial diffusion and ${\mathcal{R}}_i$ is the reaction term, describing the interactions between the substances, often involving nonlinear functions of the concentrations.

Equation of the form (1) has been applied to model many physical phenomena. For example, in chemistry, reaction-diffusion models describe how chemical species react and diffuse in a medium. Classic examples include the Belousov-Zhabotinsky reaction, which exhibits oscillatory behavior, and the formation of patterns in chemical gardens [8]. Reaction-diffusion models are fundamental in developmental biology, explaining how patterns such as animal coat markings, spatial organization of cells, and morphogenesis (the development of the structure of an organism) arise [8,9]. Alan Turing’s seminal work introduced the concept that reaction-diffusion systems can generate stable patterns, known as Turing patterns [10]. In ecology, these models describe the spread of species, the interaction between predator and prey populations, and the dynamics of ecosystems [11,12,13,14,15,16]. They help understand phenomena like population waves, species invasion, and spatial distribution of organisms [17,18]. In neuroscience, reaction-diffusion models are used to simulate the propagation of electrical signals in neurons and the dynamics of neurotransmitter diffusion in the brain. They also model the activity patterns in neural networks, explaining wave-like behaviors and pattern formation in brain activity.

A fractional reaction-diffusion equation is an extension of the classical reaction-diffusion equation, incorporating the concept of fractional calculus. This approach allows for more accurate modeling of anomalous diffusion and non-local interactions, which are observed in many complex systems where the assumptions of standard diffusion are not valid. The benefits of fractional differential equations are particularly evident in modeling the mechanical and electrical characteristics of real materials, as well as in describing the rheological properties of rocks. Additionally, they are extensively used in various scientific and engineering disciplines, including fluid dynamics, practical applications, diffusive transport analogous to diffusion, electrical circuits, electromagnetic theory, and probability [1,2,5,19,20].

The fractional reaction-diffusion equation can be written as with conditions where ${\displaystyle \frac{{\partial }^{\tau }\vartheta (x,t)}{\partial t^{\tau }}}$ is the fractional derivative of order $0<\tau \le 1$, defined by the Caputo operator as

Clearly, the Caputo derivative ${\mathfrak{D}}_t^{\tau }$ is made up of derivatives ${\mathfrak{D}}_t^{\tau -1}$ and ${\mathcal{D}}_t$, that is,

Given this make-up, it follows that

where $\mathcal{R}(x,t)={\mathcal{J}}_t^{\tau -1}R(x,t),$ and ${\mathcal{J}}_t^{\eta }$ denotes the Riemann-Liouville fractional integral derivative of order $\eta >0$ formulated as

Readers are referred to [21] and (4), for details of the proof.

Numerical methods play a crucial role in approximating solutions for time-fractional reaction-diffusion equations, which model processes exhibiting anomalous diffusion and memory effects. Among the popular techniques, finite difference methods (FDM) are widely utilized due to their simplicity and straightforward implementation [22,23]. These methods discretize the time-fractional derivatives using schemes like the Grunwald-Letnikov or L1 approximation, enabling the transformation of the continuous problem into a system of algebraic equations. The stability and convergence of these methods are often analyzed to ensure accurate solutions, with recent studies focusing on improving their efficiency and accuracy for complex problems. In contrast, finite volume methods (FVM) and spectral methods offer alternative approaches that can provide higher accuracy, especially for problems with complex geometries or requiring high-resolution solutions [23,24]. Finite volume methods, which conserve fluxes across control volumes, are particularly effective for handling heterogeneous media and ensuring conservation laws. Spectral methods, leveraging orthogonal polynomials or Fourier series, provide exponential convergence for smooth problems and have been adapted for fractional derivatives using fractional calculus [25,26]. Recent advancements include hybrid methods that combine FDM, FVM, implicit-explicit (IMEX), spectral methods and a family of Adams-Bashforth methods [27,28], to leverage their respective strengths, offering enhanced stability and accuracy for simulating time-fractional reaction-diffusion processes in various applications [29,30].

The rest of the paper is organized as follows. In Section 2, we explore various numerical approximation methods and analyze their convergence properties. Section 3 introduces three significant dynamic systems and discusses their applications in various fields. Section 4 focuses on high-dimensional numerical experiments that examine pattern formation arising from these dynamic examples. Finally, the conclusion summarizes the findings and suggests future research directions.

Here, we provide details of numerical approximation of the Caputo reaction-diffusion equation. To discretize, we let ${\mathsf{\Omega }}_{\hslash }=\left\{x_j|0\le j\le P\right\}$ to be a equal mesh on integration interval $[0,L]$, where $x_j=j\hslash ,0\le j\le P$ with $\hslash =L/P$. In the same manner, we let ${\mathsf{\Omega }}_{\zeta }=\left\{t_s|0\le j\le Q\right\}$ be a uniform mesh of interval $[0,T]$, where $t_s=s\zeta ,0\le s\le Q$ with $\zeta =T/Q$. Assume $U=\left\{U_j|0\le j\le P,U_0=U_P=0\right\}$ and $Z=\left\{Z_j|0\le j\le P,Z_0=Z_P=0\right\}$ are the grid functions on ${\mathsf{\Omega }}_{\hslash }.$ We use a similar notation as suggested in [31]:

It is not difficult to verify that

According to Lubich [32], there exists a numerical method for approximating the Riemann-Liouville (RL) fractional integral.

where ${\rho }_0=1,{\rho }_i=(-1{)}^i\left(\begin{array}{l}1-\tau \\ i\end{array}\right)$ for $i\ge 1.$ By applying the above lemma, one gets a finite difference method for the approximation of (5) as where $0\le s\le Q-1,1\le i\le M-1,{\vartheta }_i^j$ stands for the numerical solution of $\vartheta (x_i,t_j),{\varphi }_i=\varphi \left(x_i\right),{\mathcal{R}}_i^{s+1}=\mathcal{R}(x_i,t_{s+1})$. The terms ${\vartheta }_{xx}(x,t)$ and ${\vartheta }_t(x,t)$ are respectively substituted with the central difference scheme and classical backward Euler method.

Next, we rearrange the method above to get where $\beta ={\displaystyle \frac{{\zeta }^{\tau }}{h^2}}.$ We shall refer to this as first method. Likewise, using the Crank-Nicolson method to discretize (6) results in another scheme

It should be mentioned that either of the first (9) or (10) schemes forms a tridiagonal system of linear equations, which can be solved at each time step using the Thomas algorithm.

In these sections, we present and briefly discuss some practical reaction-diffusion models that continue to be of significant interest due to their wide-ranging applications in various fields. These models are essential for understanding complex systems in both natural and engineered environments.

In this section, we carry out some numerical experiment to investigate dynamic behavior of the Allen-cahn equation (35), KPP-Fisher equation (36) and the Ginzburg-Landau equation (38) in one, two and three dimensional spaces on Alienware computer using the Matlab R2021a software.

At first, one requires to justify the performance of schemes (9) and (10) before solving the main problems. To achieve this, we consider the Caputo time-fractional diffusion problem (2) as

With choice we compute with exact solution $\vartheta (x,t)=t^2\mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi x\right)$. The maximum error is displayed in Table 1 for different instances of $t$ and fractional power $\tau$ with both schemes.

For all the experiments, we utilize the homogeneous (zero-flux) boundary conditions and the following initial conditions and

The dynamic behaviour of the fractional Allen-Cahn equation (35) in 1D, 2D and 3D are displayed in Figure 1, Figure 2 and Figure 3, respectively, with $\phi =10$. Figure 1 is obtained for varying $\tau$ and $\alpha =0.5$ depicts the exact behavior of the standard Allen-Cahn equation in [45]. In the experiments, the equation exhibits stable equilibria at $\vartheta =\pm 1$ and an unstable equilibrium at $\vartheta =0$. A notable aspect of this equation is its metastability. Solutions close to $\pm 1$ tend to have flat surface, with the interface between these regions remaining stable for extended periods before experiencing a sudden change. In 2D, we observed different spiral patterns as shown in Figure 2. The 3D dynamics for subdiffusive $(0<\tau <1)$ and superdiffusive $(1<\tau <2)$ is displayed in Figure 3 showing both stable and unstable evolution.