Comprehensive Numerical Analysis of Time-Fractional Reaction-Diffusion Models with Applications to Chemical and Biological Phenomena

1. Introduction

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

$${\displaystyle \frac{\partial {\vartheta }_i}{\partial t}}=d_i\mathsf{\Delta }{\vartheta }_i+{\mathcal{R}}_i({\vartheta }_1,{\vartheta }_2,\dots ,{\vartheta }_n)$$

(1)

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

$${\displaystyle \frac{{\partial }^{\tau }\vartheta (x,t)}{\partial t^{\tau }}}=d{\mathsf{\Delta }}^2\vartheta (x,t)+\mathcal{R}\left(\vartheta \right(x,t\left)\right),x\in [0,L],0<\tau \le 1,\mathrm{o}\mathrm{r}1,\tau \le \mathrm{2,0}<t\le T$$

(2)

with conditions

$$\vartheta (0,t)=0,\vartheta (L,t)=0,\vartheta (x,0)=g_0\left(x\right),{\displaystyle \frac{\partial \vartheta (x,0)}{\partial t}}=g_L\left(x\right),t>0,x\in [0,L],L\gg 0,$$

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

$${\displaystyle \frac{{\partial }^{\tau }\vartheta (x,t)}{\partial t^{\tau }}}={\mathfrak{D}}_t^{\tau }\vartheta (x,t)={\displaystyle \frac{1}{\beta (2-\tau )}}{\int }_0^t\mathrm{‍}(t-\xi {)}^{1-\alpha }{\displaystyle \frac{{\partial }^2\vartheta (x,\xi )}{\partial {\xi }^2}}d\xi .$$

(3)

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

$$\begin{array}{ll}{\mathfrak{D}}_t^{\tau }\vartheta (x,t)& ={\displaystyle \frac{1}{\beta (1-(\tau -1)}}{\int }_0^t\mathrm{‍}(t-\xi {)}^{-(\tau -1)}{\displaystyle \frac{\partial }{\partial \xi }}{\displaystyle \frac{\partial \vartheta (x,\xi )}{\partial \xi }}d\xi \\ & ={\mathfrak{D}}_t^{\tau -1}{\mathcal{D}}_t\vartheta (x,t).\end{array}$$

(4)

Given this make-up, it follows that

Theorem 1.1

The Caputo reaction-diffusion equation (2) is equivalent to integro- differential equation

$$\vartheta (x,t)=g_L\left(x\right)+{\displaystyle \frac{d}{\beta (\tau -1)}}{\int }_0^t\mathrm{‍}(t-\xi {)}^{(\tau -2)}{\displaystyle \frac{\partial \vartheta (x,\xi )}{\partial x^2}}d\xi +\mathcal{R}(x,t),$$

(5)

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

$${\mathcal{J}}_t^{\eta }R(x,t)={\displaystyle \frac{1}{\beta \left(\eta \right)}}{\int }_0^t\mathrm{‍}(t-\xi {)}^{\eta -1}R(x,\xi )d\xi .$$

(6)

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.

2. Numerical Methods

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]:

$$\begin{array}{l}{\psi }_{+}{U}_j=U_{j+1}-U_j,{‖U‖}_{\infty }=\underset{1\le j\le P-1}{\mathrm{m}\mathrm{a}\mathrm{x}}\left|U_1\right|{‖U‖}_1={\displaystyle \sum _{j=1}^{P-1}}\mathrm{‍}\hslash \left|U_j\right|,\end{array}$$

$$<U,Z>={\displaystyle \sum _{j=1}^{P-1}}\mathrm{‍}\mathrm{\hslash }{U}_j{Z}_j,{‖U‖}^2=<U,{\psi }^2{U}_j=<U,U>,{\psi }^2{U}_j=U_{j+1}-2U_j+U_{j-1}.$$

It is not difficult to verify that

$$<{\psi }^{2}U,Z>=-<{\psi }_{+}U,{\psi }_{+}Z>.$$

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

Lemma 2.1

Suppose $R\left(t\right)=t^{\eta -1}\vartheta \left(t\right)$ with  $\eta >0$ and $g$ continuous, then

$${\mathcal{J}}_t^{\tau -1}R\left(t_{s+1}\right)={\zeta }^{\tau -1}{\sum }_{i=0}^{s+1}\mathrm{‍}{\rho }_{i}R\left(t_{s+1-i}\right)+\mathcal{O}\left(\zeta \right),$$

(7)

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

$${\displaystyle \frac{{\vartheta }_i^{s+1}-{\vartheta }_i^s}{\zeta }}={\varphi }_i+d{\zeta }^{\tau -1}{\sum }_{i=0}^{s+1}\mathrm{‍}{\mathsf{\Omega }}_i{\displaystyle \frac{{\vartheta }_{i+1}^{s+1-i}-2{\vartheta }_i^{s+1-i}+{\vartheta }_{i-1}^{s+1-i}}{h^2}}+{\mathcal{R}}_i^{s+1},$$

(8)

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

$${\vartheta }_i^{s+1}-{\vartheta }_i^s=c\beta {\sum }_{i=0}^{s+1}\mathrm{‍}{\rho }_i{\psi }^2{\vartheta }_i^{s+1-i}+\zeta {\varphi }_i+\zeta {\mathcal{R}}_i^{s+1},$$

(9)

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

$${\vartheta }_i^{s+1}-{\vartheta }_i^s={\displaystyle \frac{d\beta }{2}}{\sum }_{i=0}^{s+1}\mathrm{‍}{\rho }_i{\psi }^2{\vartheta }_i^{s+1-i}+{\displaystyle \frac{d\beta }{2}}{\sum }_{i=0}^s\mathrm{‍}{\rho }_i{\psi }^2{\vartheta }_i^{s-i}+\zeta {\varphi }_i+{\displaystyle \frac{\zeta }{2}}\left({\mathcal{R}}_i^{s+1}+{\mathcal{R}}_i^s\right).$$

(10)

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.

2.1. Convergence and Solvability Properties of the Difference Schemes

Hereafter, let $D_1$ denote a positive constant that does not depend on $\zeta ,h,i,$ or $n$, though its value may vary in different contexts. We begin by applying the Gershgorin circle theorem (GCT) to demonstrate the solvability of numerical formulae (9) and (10).

Theorem 2.1

The schemes (9) and (10) are uniquely-solvable.

Proof. 

Since the two methods are similar, we present the solvability of only scheme (9) here in what follows

$$\begin{array}{ll}& -d\beta {\vartheta }_{i-1}^{s+1}+(1+2d\beta ){\vartheta }_i^{s+1}-d\beta {\vartheta }_{i-1}^{s+1}\\ & ={\vartheta }_i^s+d\beta {\sum }_{i=1}^{s+1}\mathrm{‍}{\rho }_i\left({\vartheta }_{i+1}^{s+1-i}-2{\vartheta }_i^{s+1-i}+{\vartheta }_{i-1}^{s+1-i}\right)+\zeta {\varphi }_i+\zeta F_i^{s+1}.\end{array}$$

(11)

Assume

$$\mathcal{A}=\left(\begin{array}{llll}1+2d\beta & -d\beta & \dots & \dots \\ -d\beta & 1+2d\beta & -d\beta & \dots \\ .& .& .& .\\ .& .& .& .\\ .& -d\beta & 1+2d\beta & -d\beta \\ \dots & \dots & -d\beta & 1+2d\beta \\ & & & \end{array}\right),\mathcal{B}=\left(\begin{array}{llll}-2& 1& \dots & \dots \\ 1& -2& 1& \dots \\ .& .& .& .\\ .& .& .& .\\ .& 1& -2& 1\\ \dots & \dots & 1& -2\\ & & & \end{array}\right)$$

(12)

$$\mathcal{C}=\left(\begin{array}{lllll}1-2d\beta {\rho }_1& d\beta {\rho }_1& \dots & \dots & \\ d\beta & 1-2d\beta {\rho }_1& d\beta & \dots & \dots \\ .& .& .& .& \\ .& .& .& .& \\ .& -d\beta {\rho }_1& 1-2d\beta {\rho }_1& -d\beta {\rho }_1& \\ \dots & \dots & -d\beta {\rho }_1& 1-2d\beta {\rho }_1& \\ & & & & \end{array}\right),$$

(13)

$${\vartheta }^n={\left({\vartheta }_1^s,{\vartheta }_2^s,\dots ,{\vartheta }_{P-2}^s,{\vartheta }_{P-1^s}\right)}^T,\varphi ={\left({\varphi }_1,{\varphi }_2,\dots ,{\varphi }_{1}P-2,{\varphi }_{1}P-1\right)}^T,$$

$${\mathcal{R}}^s={\left({\mathcal{R}}_1^s,{\mathcal{R}}_2^s,\dots ,{\mathcal{R}}_{P-2,{\mathcal{R}}_{P-1}^s}^s\right)}^T.$$

Bear in mind that we adopt the homogeneous boundary conditions for our Caputo time-fractional reaction-diffusion problem (2), which allows equation (11) to take the matrix form

$$\mathcal{A}{\vartheta }^{s+1}=\mathcal{C}{\vartheta }^s+d\beta {\sum }_{i=2}^{s+1}\mathrm{‍}{\rho }_i\mathcal{B}{U}^{s+1-i}+\zeta \left(\varphi +F^{s+1}\right).$$

(14)

By the GCT, the matrice $A$ is invertible, which proves that scheme (9) is solvable.

The subsequent lemma, which addresses the nonnegative nature of specific real quadratic forms with a convolution structure, is attributed to [31]. It is crucial here for the convergence analysis study.

Lemma 2.2

Suppose ${\left\{a_n\right\}}_{s=0}^{\infty }$ is monotonically decreasing function of positive real numbers with condition $a_{s+1}+a_{s-1}\le 2a_s(s\le 1),$ then for any integer $i>0$ and real vector $(U_1,U_2,\dots ,U_N{)}^T\in {\mathbb{R}}^N,$ satisfies that

$${\sum }_{s=0}^{N-1}\mathrm{‍}\left({\sum }_{p=0}^s\mathrm{‍}{a}_p{U}_{s+1-p}\right)U_{s+1}\ge 0.$$

Lemma 2.3

The sequence ${\left\{{\rho }_s\right\}}_{s=0}^{\infty }$ defined in (7) holds for ${\rho }_n>0,{\rho }_s<{\rho }_{s-1},{\rho }_{s+1}+{\rho }_{s-1}>2{\rho }_s(s\le 1).$

Proof. 

We observe

$$\begin{array}{ll}{\rho }_s& =(-1{)}^s\left(\begin{array}{l}1-\tau \\ i\end{array}\right)=(-1{)}^s{\displaystyle \frac{\beta (2-\tau )}{\beta (s+1)\beta (2-\tau -s)}}\\ & =(-1)(-1{)}^{s-1}{\displaystyle \frac{\beta (2-\tau )\beta (2-\tau -s)}{s\beta \left(s\right)\beta (3-\tau -s)}}={\displaystyle \frac{s-(2-\tau )}{s}}{\rho }_{s-1}\\ & =\left(1-{\displaystyle \frac{2-\tau }{s}}\right){\rho }_{s-1}.\end{array}$$

(15)

Since $1<\tau <2$ and ${\rho }_0=1,0<1-{\displaystyle \frac{2-\tau }{s+1}}<1,$ then ${\rho }_s>0,{\rho }_s<{\rho }_{s-1}.$

In addition,

$$\begin{array}{ll}{\rho }_{s+1}-2{\rho }_s+{\rho }_{s-1}& =\left(1-{\displaystyle \frac{2-\tau }{s+1}}\right)\left(1-{\displaystyle \frac{2-\tau }{s}}\right){\rho }_{s-1}\\ & -2\left(1-{\displaystyle \frac{2-\tau }{s}}\right){\rho }_{s-1}+{\rho }_{s-1}\\ & ={\displaystyle \frac{(2-\tau )(3-\tau )}{s(s+1)}}{\rho }_{s-1}>0.\end{array}$$

(16)

The proof is completed.

2.2. Analysis of Convergence for Scheme (9)

Assume

$$\left|{\vartheta }_{tt}\right|\le D_1,\left|{\vartheta }_{xxxx}\right|\le D_1\mathrm{f}\mathrm{o}\mathrm{r}(x,t)\in [0,L]\times [0,T],$$

(17)

then for (9), we obtain

$$\vartheta (x_i,t_{s+1})-\vartheta (x_i,t_s)=d\beta {\sum }_{i=0}^{s+1}\mathrm{‍}{\rho }_i{\psi }^2\vartheta (x_i,t_{s+1-i})+\zeta {\varphi }_i+\zeta F_i^{s+1}+O({\zeta }^2+\zeta h^2).$$

(18)

Theorem 2.2

Suppose $\vartheta (x,t)$ and $\left\{{\vartheta }_i^s\right\}$ are the numerical solutions to time-fractional reaction-diffusion problems (2) and (9) respectively, and that $\vartheta (x,t)$ meets the smoothness criteria of (17). Then, for sufficiently small $\zeta$ and $h$,

$$‖e^s‖=O(\zeta +h^2),s\le 1,$$

(19)

where $e_i^s=\vartheta (x_i,t_s)-{\vartheta }_i^s$.

Proof. 

By taking (9) from (18), one gets the error equation

$$e_i^{s+1}-e_i^s=d\beta {\sum }_{i=0}^s\mathrm{‍}{\rho }_i{\psi }^2{e}_i^{s+1-i}+r_i^{s+1},\mathrm{f}\mathrm{o}\mathrm{r}0\le s\le Q-1,1\le i\le P-1$$

(20)

where $r_i^{s+1}=O({\zeta }^2+\zeta h^2).$

Next, we multiply on both sides of (20) by $h{e}_i^{s+1}$ and add up for $i$ from 1 to $M-1$, to get

$${‖e^{s+1}‖}^2-<e^s,e^{s+1}>=d\beta {\sum }_{i=0}^s\mathrm{‍}{\rho }_i<{\psi }^2{e}^{s+1-i},e^{s+1}>+<r^{s+1},e^{s+1}>.$$

(21)

In particular,

$${‖e^{s+1}‖}^2-<e^s,e^{s+1}>=d\beta {\sum }_{i=0}^s\mathrm{‍}{\rho }_i<{\psi }^2{e}^{s+1-i},e^{s+1}>+<r^{s+1},e^{s+1}>,$$

$${‖e^s‖}^2-<e^{s-1},e^s>=d\beta {\sum }_{i=0}^{s-1}\mathrm{‍}{\rho }_i<{\psi }^2{e}^{s-i},e^s>+<r^s,e^s>,$$

$${‖e^2‖}^2-<e^1,e^2>=d\beta {\sum }_{i=0}^{s-1}\mathrm{‍}{\rho }_i<{\psi }^2{e}^{2-i},e^2>+<r^2,e^2>,$$

$${‖e^1‖}^2=d\beta {\sum }_{i=0}^{s-1}\mathrm{‍}{\rho }_i<{\psi }^2{e}^{1-i},e^1>+<r^1,e^1>,$$

(22)

Applying the inequality $<V,W>\le ‖V‖\cdot ‖W‖\le {\displaystyle \frac{1}{2}}\parallel V{|}^2+{\displaystyle \frac{1}{2}}{‖W‖}^2,$ we have

$${‖e^{s+1}‖}^2-{\displaystyle \frac{1}{2}}{‖e^{s+1}‖}^2-{\displaystyle \frac{1}{2}}{‖e^s‖}^2=d\beta {\sum }_{i=0}^s\mathrm{‍}{\rho }_i<{\psi }^2{e}^{s+1-i},e^{s+1}>+<r^{s+1},e^{s+1}>,$$

$${‖e^s‖}^2-{\displaystyle \frac{1}{2}}{‖e^s‖}^2-{\displaystyle \frac{1}{2}}{‖e^{s-1}‖}^2=d\beta {\sum }_{i=0}^s\mathrm{‍}{\rho }_i<{\psi }^2{e}^{s-i},e^s>+<r^s,e^s>,$$

$${‖e^2‖}^2-{\displaystyle \frac{1}{2}}{‖e^2‖}^2-{\displaystyle \frac{1}{2}}{‖e^{s1}‖}^2=d\beta {\sum }_{i=0}^s\mathrm{‍}{\rho }_i<{\psi }^2{e}^{2-i},e^2>+<r^2,e^2>,$$

$${‖e^1‖}^2=d\beta {\sum }_{i=0}^0\mathrm{‍}{\rho }_i<{\psi }^2{e}^{1-i},e^1>+<r^1,e^1>,$$

By summing up the inequities above results to

$$\begin{array}{ll}{‖e^{s+1}‖}^2& \le 2d\beta {\sum }_{p=0}^s\mathrm{‍}{\sum }_{i=0}^s\mathrm{‍}{\rho }_i<{\psi }^2{e}^{p+1-i},e^{p+1}>+2{\sum }_{p=0}^{s+1}\mathrm{‍}<r^p,e^p>\\ & =-2d\beta {\sum }_{p=0}^s\mathrm{‍}{\sum }_{i=0}^s\mathrm{‍}{\rho }_i<\psi e^{p+1-i},{\psi }_{+}{e}^{p+1}>+2{\sum }_{p=0}^{s+1}\mathrm{‍}<r^p,e^p>\\ & \le 2\sqrt{L}\mathcal{M}(s+1)D_1({\zeta }^2+\zeta h^2)\le 2\sqrt{L}\mathcal{M}T{D}_1\left(\zeta +h^2\right)\end{array}$$

where $\mathcal{M}={\mathrm{m}\mathrm{a}\mathrm{x}}_{1\le p\le Q}‖e^p‖,$ and have applied the inequality

$$\parallel e^p{|}_1\le \sqrt{Ph}‖e^p‖\le \sqrt{L}‖e^p‖,$$

and Lemma 11 result. Thus the condition (19) holds.

Using a method analogous to that employed in the convergence proof, we can demonstrate the stability of the scheme (9), that is,

$$‖{\vartheta }^{s+1}‖\le D_1\left(‖\varphi ‖+‖\varphi ‖+\zeta {\sum }_{i=1}^{s+1}\mathrm{‍}‖{\mathcal{R}}^i‖\right),0\le s\le Q-1.$$

(23)

2.3. Analysis of Convergence for Scheme (10)

Assume ${\vartheta }_{xxxx}$ and ${\vartheta }_t$ are continuous and differentiable on interval $[0,L]\times [0,T]$, with $\left|{\vartheta }_{tt}\right|\le D_1{t}^{\tau -2}$. This assumption is likely the minimal condition required for the existence of the Caputo-in-time derivative of order $0<\tau <0$ (for subdiffusion process) or $1<\tau <2$ (in the case of superdiffusion scenarios).

For any $R\left(t\right)\in D^1\left(\right[0,T\left]\right)\bigcap \mathrm{‍}{D}^3\left(\right(0,T\left]\right)$ satisfying $\left|R\text{'}\text{'}\left(t\right)\right|\le D_1{t}^{\tau -2}$ and $\left|f\text{'}\text{'}\left(t\right)\right|\le D_1{t}^{\tau -3}$ as $t\to 0$. We will approximate ${\mathcal{J}}_t^{\tau }-1R\left(t\right)$ numerically. We can generalize the result in [33] to the formulas as follows

For any $R\left(t\right)\in D^1\left(\right[0,T\left]\right)\bigcap \mathrm{‍}{D}^3\left(\right(0,T\left]\right)$ that satisfies $\left|R\text{'}\text{'}\left(t\right)\right|\le D_1{t}^{\tau -2}$ and $\left|R\text{'}\text{'}\left(t\right)\right|\le D_1{t}^{\tau -3}$ as $t\to 0$, we approximate ${\mathcal{J}}_t^{\tau -1}R\left(t\right)$ numerically. The result from [33] can be generalized to the following formulas.

$${\mathcal{J}}_t^{\tau -1}R\left(t_{s+1/2}\right)={\displaystyle \frac{1}{2}}\left[{\mathcal{J}}_t^{\tau -1}f\left(t_s\right)+{\mathcal{J}}_t^{\tau -1}R\left(t_{s+1}\right)\right]+\mathcal{O}\left({\zeta }^2{t}_s^{\tau -3}\right),s\ge 1$$

(24)

and

$${\mathcal{J}}_t^{\tau -1}R\left(t_{1/2}\right)={\displaystyle \frac{1}{2}}{\mathcal{J}}_t^{\tau -1}R\left(t_1\right)+\mathcal{O}\left({\zeta }^{\tau -1}\right).$$

(25)

Since $t_s^{\tau -3}\le 2^{3-\tau }{t}_{s+1}^{\tau -3}$ for $s\le 1,$ then results to

$${\mathcal{J}}_t^{\tau -1}R\left(t_{s+1/2}\right)={\displaystyle \frac{1}{2}}\left[{\mathcal{J}}_t^{\tau -1}R\left(t_s\right)+{\mathcal{J}}_t^{\tau -1}R\left(t_{s+1}\right)\right]+\mathcal{O}\left({\zeta }^2{t}_{s+1}^{\tau -3}\right),\mathrm{f}\mathrm{o}\mathrm{r}s\ge 0.$$

(26)

Furthermore, we can verify that

$$R\text{'}\left(t_{s+1/2}\right)={\displaystyle \frac{R\left(t_{s+1}\right)-R\left(t_s\right)}{\zeta }}+O\left({\zeta }^2{t}_{s+1}^{\tau -3}\right),s\ge 0.$$

(27)

When equations (7), (26) and (27) are combined, we have

$$\begin{array}{ll}\vartheta (x_i,t_{s+1})-\vartheta (x_i,t_s)& ={\displaystyle \frac{d\beta }{2}}{\sum }_{i=0}^{s+1}\mathrm{‍}{\rho }_i{\psi }^2\vartheta (x_i,t_{s+1-i})+{\displaystyle \frac{d\beta }{2}}{\sum }_{i=0}^s\mathrm{‍}{\rho }_i{\psi }^2\vartheta (x_i,t_{s-i})\\ & +\zeta {\varphi }_i+{\displaystyle \frac{\zeta }{2}}\left({\mathcal{R}}_i^{s+1}+{\mathcal{R}}_i^s\right)+\mathcal{O}\left({\zeta }^2+\zeta h^2\right).\end{array}$$

(28)

Next, we report the error estimate for the second scheme (10) as follows.

Theorem 2.3

Suppose that the solution to Caputo reaction-diffusion equation (2) meets the given smoothness requirements. Denote $u(x,t)$ as the solution of (2) and $\left\{u_i^n\right\}$ be the solution of numerical scheme (10). Then, as $\zeta$ and $h$ independently approach zero

$$\left|e^s\right|=\mathcal{O}\left(\zeta +h^2\right),s\ge 1.$$

(29)

Proof. 

Taking (7) from (28), we get

$$\begin{array}{ll}{e}_i^{s+1}-e_i^s& ={\displaystyle \frac{d\beta }{2}}{\sum }_{i=0}^s\mathrm{‍}{\rho }_i\left({\psi }^2{e}_i^{s+1-i}+{\psi }^2{e}_i^{s-i}\right)\\ & +R_i^{s+1},0\le s\le Q-1,1\le i\le P-1.\end{array}$$

(30)

where $R_i^{s+1}=O\left({\zeta }^3{t}_{s+1}^{\tau -3}+{\zeta }^2+\zeta h^2\right).$ We multiply $h(e_i^{s+1}+e_i^s)$ on both sides of equation (30), and adding up in terms of $i$, to have

$$\begin{array}{ll}{‖e^{s+1}‖}^2-{‖e^s‖}^2& ={\displaystyle \frac{d\beta }{2}}{\sum }_{i=0}^s\mathrm{‍}{\rho }_i<{\psi }^2\left(e^{s+1-i}+e^{s-i}\right),e^{s+1}+e^s>\\ & <R^{s+1},e^{s+1}+e^s>.\end{array}$$

(31)

Adding over $n$ results

$$\begin{array}{ll}{‖e^{s+1}‖}^2-{‖e^s‖}^2& ={\displaystyle \frac{d\beta }{2}}{\sum }_{p=0}^s\mathrm{‍}{\sum }_{i=0}^s\mathrm{‍}{\rho }_i<{\psi }^2\left(e^{p+1-i}+e^{p-i}\right),e^{p+1}+e^p>\\ & +{\sum }_{p=0}^s\mathrm{‍}<R^{p+1},e^{p+1}+e^p>\\ & \le -{\displaystyle \frac{D\beta }{2}}{\sum }_{p=0}^s\mathrm{‍}{\sum }_{i=0}^s\mathrm{‍}{\rho }_i<{\psi }_{+}\left(e^{p+1-i}+e^{p-i}\right),{\psi }_{+}\left(e^{p+1}+e^p\right)>\\ & +{\sum }_{p=0}^s\mathrm{‍}{\left|R^{p+1}\right|}_{\infty }\left(‖e^{p+1}+e^p‖+{‖e^p‖}_1\right)\\ & \le D_1\mathcal{M}\left(\zeta +h^2+{\zeta }^{\tau }{\sum }_{p=0}^s\mathrm{‍}(p+1{)}^{\tau }-3\right).\end{array}$$

(32)

For superdiffusive process ($1<\tau <2$), the series ${\sum }_{p=0}^s\mathrm{‍}(p+1{)}^{\tau -3}$ in (32) is convergent, which satisfies

$$‖e^s‖=\mathcal{O}(\zeta +h^2),\mathrm{f}\mathrm{o}\mathrm{r}s\ge 1.$$

(33)

In the same way, it can be verified that the numerical scheme (10) is stable and satisfies the inequality condition

$$‖{\vartheta }^{s+1}‖\le D_1\left(‖\varphi ‖+‖\varphi ‖=\zeta {\sum }_{i=0}^s\mathrm{‍}‖{\mathcal{R}}^{i+1}+{\mathcal{R}}^i‖\right),0\le s\le Q-1.$$

(34)

Based on the error estimates (19) and (29), both schemes exhibit first-order accuracy in time and second-order accuracy in space. However, Scheme 2 has less stringent smoothness requirements compared to Scheme 1, allowing it to be applied across a wider range of scenarios.

In the works [31,33,34], numerical discussions are presented for a partial integro-differential equation similar to (6) with a fractional order of $\tau =1.5$. Sanz-Serna introduces a temporal semi-discrete first-order algorithm. Lopez proposes a difference scheme based on the backward Euler method, with an error estimation of $\mathcal{O}\left(\zeta \right|\mathrm{l}\mathrm{n}\zeta {|}^{1/2}+h^2)$. Tang employs the Crank-Nicolson method and the product trapezoidal method to develop a difference scheme with an accuracy of $\mathcal{O}({\zeta }^{3/2}+h^2)$. However, it remains uncertain whether this approach can be extended to any $\tau \in \left(\mathrm{1,2}\right)$.

3. Model Equations

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.

3.1. Allen-Cahn Equation

The Allen-Cahn equation, named after John W. Cahn and Sam Allen, is a pivotal reaction-diffusion equation in mathematical physics. This equation is crucial for modeling the process of phase separation in multi-component alloy systems, which includes phenomena such as order-disorder transitions. Phase separation is a fundamental process where a homogeneous mixture of components evolves into distinct regions or phases, each enriched in different components. This process is driven by the system’s tendency to minimize its free energy. The Allen-Cahn (AC) equation provides a mathematical framework to describe the temporal evolution of this separation, capturing the complex interplay between reaction kinetics and diffusion [35].

Mathematically, the Allen-Cahn equation is typically written as:

$${\mathfrak{D}}_t^{\tau }\vartheta =\alpha {\mathsf{\Delta }}^2\vartheta -\mathcal{N}\left(\vartheta \right)$$

(35)

Here, $\vartheta =\vartheta (x,t)$ represents the order parameter, which indicates the local state of the system (e.g., the concentration of a particular component). The term $\alpha \mathsf{\Delta }\vartheta$ represents the diffusion term, where $\mathsf{\Delta }$ is the Laplacian operator, reflecting how the order parameter spreads out over time. The function $\mathcal{N}\left(\vartheta \right)$ is a nonlinear term derived from the free energy of the system, often taking the form of the derivative of a double-well potential, such as $\mathcal{N}\left(\vartheta \right)=\vartheta -{\vartheta }^3$.

In reaction-diffusion dynamics, the AC equation combines the effects of local reactions ()described by $\mathcal{N}\left(\vartheta \right))$ and diffusion (described by $\alpha \mathsf{\Delta }\vartheta$). This combination allows it to model the competition between different phases in a material. In multi-component alloy systems, the equation is particularly useful for studying order-disorder transitions. These transitions occur when a system changes from a disordered state, where components are randomly distributed, to an ordered state, where components are arranged in a regular pattern. The equation also describes the dynamics of interfaces between different phases. As time progresses, these interfaces move and evolve according to the balance between diffusive spreading and the local reaction kinetics [36]. In materials science, the Allen-Cahn equation is widely used to simulate and understand the microstructural evolution of alloys. For instance, it can predict how different phases nucleate, grow, and coarsen over time, providing insights into the material’s properties.

The Allen-Cahn equation has been extensively studied both analytically and numerically. Analytical approaches often involve studying the stability and bifurcation of solutions, which provide insights into the conditions under which different phases form and evolve. Numerical simulations, on the other hand, allow for the detailed visualization of phase separation processes and the exploration of complex scenarios that are analytically intractable. The AC equation is closely related to the Cahn-Hilliard equation, which also describes phase separation but focuses on conserved quantities. While the Allen-Cahn equation models non-conserved order parameters (like the local state of order), the Cahn-Hilliard equation deals with conserved quantities (like concentration). Both equations are fundamental in the study of pattern formation and are used to describe different aspects of phase separation dynamics.

3.2. The KPP-Fisher Equation

In mathematical physics, the KPP-Fisher equation, named after Andrey Kolmogorov, Ivan Petrovsky, Nikolai Piskunov, and Ronald Fisher, is a prominent partial differential equation. It is also commonly referred to as the KPP equation, the Fisher equation [37], or the Fisher-KPP equation [38,39]. This equation is significant in various fields, including population genetics, ecology, and the study of reaction-diffusion systems.

The general KPP-Fisher equation is typically expressed as:

$${\mathfrak{D}}_t^{\tau }\vartheta =\delta {\mathsf{\Delta }}^2\vartheta +\kappa \vartheta (1-{\vartheta }^q)$$

(36)

Here, $\vartheta =\vartheta (x,t)$ represents the density of the population or the concentration of a chemical substance at position $x$ and time $t$. The term $\delta {\mathsf{\Delta }}^2\vartheta ={\displaystyle \frac{{\partial }^2\vartheta }{\partial x^2}}$ is the Laplacian operator that represents the diffusion component, where $\delta$ is the diffusion coefficient, indicating how the substance spreads out in space. The term $\kappa \vartheta (1-{\vartheta }^q),q>0$ is the reaction component, where $\kappa$ is the growth rate of the population or reaction rate of the chemical substance. It should be noted that whenever $\kappa =1$ and $q=2$, we recover the Allen-Cahn equation (35). The KPP-Fisher equation is related to other reaction-diffusion equations like the Allen-Cahn equation and the Cahn-Hilliard equation. While the KPP-Fisher equation primarily focuses on traveling wave solutions and logistic growth, the Allen-Cahn equation models phase separation with non-conserved order parameters and the Cahn-Hilliard equation deals with conserved quantities in phase separation processes.

The Fisher equation was originally derived to describe the spread of an advantageous gene in a population. The term $\kappa \vartheta (1-{\vartheta }^q)$ models logistic growth, where $\vartheta$ grows quickly when it is small, but slows down as it approaches the carrying capacity (normalized to 1). One of the notable features of the KPP-Fisher equation is its ability to describe traveling wave solutions. These are wave-like fronts that move with a constant speed, representing the spread of a gene, population, or chemical concentration through space over time. The minimum speed of these waves is given by $v=2\sqrt{\kappa \delta }$. In an ecological context, the equation models the spread of a species through a habitat. It accounts for both the natural diffusion of the species and their growth and competition dynamics. Similarly in the context of chemical reactions, the equation describes the spread of reactants in a medium, combining the effects of chemical kinetics and molecular diffusion.

Analytically, the KPP-Fisher equation has been studied for its traveling wave solutions, stability of these waves, and the long-term behavior of solutions. Understanding the minimum wave speed and the conditions for wave formation are key analytical challenges. Numerical simulations are extensively used to explore the dynamics of the equation in more complex scenarios, such as heterogeneous environments or higher-dimensional spaces.

3.3. Ginzburg-Landau Equation

The Ginzburg-Landau equation, named after Vitaly Ginzburg and Lev Landau, is a fundamental equation in the field of nonlinear dynamics and pattern formation. It describes the evolution of small disturbances near a finite wavelength bifurcation, where a system transitions from a stable to an unstable state. This equation plays a critical role in understanding how complex patterns emerge in various physical, chemical, and biological systems.

The Ginzburg-Landau equation is typically given as [40,41]

$${\mathfrak{D}}_t^{\tau }\vartheta =\sigma \vartheta +\delta {\mathsf{\Delta }}^2\vartheta -\beta |\vartheta {|}^2\vartheta$$

(37)

where $\vartheta =\vartheta (x,t)$ is the complex amplitude of the disturbance, $\sigma$ is a small parameter representing the distance from the bifurcation point, $\delta$ is a coefficient related to the spatial diffusion of the amplitude, and $\beta$ is nonlinear coefficient. The real part of $\vartheta$ is particularly important for describing the physical aspects of the system’s state.

At the onset of a finite wavelength bifurcation, the system becomes unstable for a critical wavenumber ${\omega }_c$ which is non-zero. This instability leads to the growth of disturbances characterized by this wavenumber. Near the bifurcation, the disturbances evolve with a Fourier mode corresponding to ${\omega }_c$ and a slowly varying amplitude $\vartheta$. The Ginzburg-Landau equation governs the dynamics of this amplitude $\vartheta$, capturing how the disturbances grow, saturate, and interact over time. It should be mentioned that for oscillatory behavior, $\vartheta$ satisfies the novel complex Ginzburg-Landau equation [42] The Ginzburg-Landau equation is typically given as

$${\mathfrak{D}}_t^{\tau }\vartheta =(1+i\delta ){\mathsf{\Delta }}^2\vartheta +\vartheta -(1+i\beta )\vartheta |\vartheta {|}^2\vartheta$$

(38)

where $\delta$ and $\beta$ are treated as real constants. The solution of (38) often results to two important modes. The non-oscillatory modes which represent steady-state patterns that do not change in time once they have fully developed. The disturbances grow until they reach a stable amplitude and form stationary spatial structures, and the oscillatory modes which lead to time-dependent patterns that oscillate as they evolve. The amplitude $\vartheta$ varies not only in space but also in time, leading to complex temporal behavior.

The Ginzburg-Landau equation is crucial in the study of pattern formation. It describes how simple initial disturbances can grow and form intricate spatial and temporal patterns [43]. This is observed in a wide range of systems, from chemical reactions (like the Belousov-Zhabotinsky reaction) to fluid dynamics (such as Rayleigh-Benard convection). In the context of superconductivity, the equation helps describe the behavior of the order parameter near the critical temperature. It provides insights into the formation of vortices and other phenomena in superconducting materials. For nonlinear optics, the GL equation models the propagation of light in nonlinear media, explaining the formation of optical patterns and solitons in lasers and other optical systems, and in biology, it helps understand the development of patterns in animal skins, the dynamics of populations, and other phenomena where spatial and temporal variations are critical.

Analytically, the Ginzburg-Landau equation is studied to understand the stability and bifurcation behavior of solutions. Researchers investigate how solutions change as parameters like $\sigma$ and $\delta$ vary, and how nonlinear interactions lead to saturation and pattern formation. Numerically, simulations of the Ginzburg-Landau equation reveal detailed behavior of the amplitude $\vartheta$ over time and space [44]. These simulations help visualize complex phenomena that are analytically intractable, providing deeper insights into the dynamics of the system.

4. Numerical Experiments and Results

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

$${\mathfrak{D}}_t^{\tau }\vartheta =d{\mathsf{\Delta }}^2\vartheta (x,t)+\mathcal{R}(x,t),x\in (0,L),t\in (0,T)\mathrm{a}\mathrm{n}\mathrm{d}\tau \in \left(\mathrm{0,1}\right).$$

(39)

With choice

$$\mathcal{R}(x,t)=\left({\displaystyle \frac{2t^{2-\tau }}{\mathsf{\Gamma }(3-\tau )}}+4{\pi }^2{t}^2\right),d=1,\vartheta (x,0)=0,\vartheta (0,t)=0,$$

we compute with exact solution $\vartheta (x,t)=t^2\mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi x\right)$. The maximum error

$${\mathcal{E}}_{max}=\underset{0\le j\le P}{\mathrm{m}\mathrm{a}\mathrm{x}}|{\overline{\vartheta }}_j-{\vartheta }_j|,{\overline{\vartheta }}_j\left(\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\right),{\vartheta }_j\left(\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\right)$$

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

$$\begin{array}{ll}\vartheta (x,0)& =\mathrm{e}\mathrm{x}\mathrm{p}(-20((x-\phi /3{)}^2)/\phi )-\mathrm{e}\mathrm{x}\mathrm{p}(-20\left(\right(x-\phi /2{)}^2)/\phi )\\ & +\mathrm{e}\mathrm{x}\mathrm{p}(-20(x-\phi {)}^2/\phi ),\end{array}$$

(40)

$$\begin{array}{ll}\vartheta (x,y,0)& =\mathrm{e}\mathrm{x}\mathrm{p}(-20((x-\phi /3{)}^2+(y-\phi /3{)}^2)/\phi )\\ & -\mathrm{e}\mathrm{x}\mathrm{p}(-20((x-\phi /2{)}^2+(y-\phi /2{)}^2)/\phi )\\ & +\mathrm{e}\mathrm{x}\mathrm{p}(-20((x-\phi {)}^2+(y-\phi {)}^2)/\phi ),\end{array}$$

(41)

$$\vartheta (x,y,0)=2\mathrm{e}\mathrm{x}\mathrm{p}(-20(x-1{)}^2)+2.05\mathrm{e}\mathrm{x}\mathrm{p}(-10(y-1{)}^2)+\mathrm{e}\mathrm{x}\mathrm{p}(-20(x+y{)}^2),$$

(42)

and

$$\begin{array}{ll}\phi (x,y,z,0)& =\mathrm{e}\mathrm{x}\mathrm{p}(-20((x-\phi /3{)}^2+(y-\phi /3{)}^2+(z-\phi /3{)}^2)/\phi )\\ & -\mathrm{e}\mathrm{x}\mathrm{p}(-20((x-\phi /2{)}^2+(y-\phi /2{)}^2+(z-\phi /2{)}^2)/\phi )\\ & +\mathrm{e}\mathrm{x}\mathrm{p}(-20((x-\phi {)}^2+(y-\phi {)}^2+(z-\phi {)}^2)/\phi ).\end{array}$$

(43)

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.

The Allen-Cahn equation is a fundamental tool in mathematical physics and materials science, providing deep insights into the process of phase separation and order-disorder transitions in multi-component systems. Its ability to model the interplay between diffusion and reaction kinetics makes it indispensable for understanding and predicting the behavior of complex materials.

Figure 4, Figure 5, Figure 6 and Figure 7 illustrate the numerical solutions of the KPP-Fisher equation (36) in one, two, and three dimensions. The 1D dynamics show a stable spatiotemporal evolution, as depicted in Figure 4 (see caption for details). The 2D dynamics resemble the Allen-Cahn distribution due to similarities between the two equations. Simulations with the initial conditions (41) and (42) produce the results shown in Figure 5 and Figure 6, respectively. The complex dynamics of the 3D evolution, starting from the initial condition in (43), are presented in Figure 7.

The KPP-Fisher equation is a fundamental partial differential equation in mathematics and applied sciences. Its ability to model the spread of populations, genes, or chemical substances makes it a versatile tool in population genetics, ecology, and chemical kinetics. The equation’s traveling wave solutions provide deep insights into how advantageous genes spread through a population, how species expand their habitats, and how chemical reactions propagate through space. Understanding and solving the KPP-Fisher equation continues to be a significant area of research, with applications in various scientific and engineering disciplines.

Finally, we present numerical results for fractional complex Ginzburg-Landau (CGL) equation (38) as displayed in Figure 8, Figure 9, Figure 10 and Figure 11. The CGL equation is known for its oscillatory and spiral distribution behavior, this ie evident in the 2D and 3D dynamics.

The Ginzburg-Landau equation is a cornerstone in the study of nonlinear dynamics and pattern formation as evident Figure 8, Figure 9, Figure 10 and Figure 11. By describing the evolution of disturbances near a finite wavelength bifurcation, it provides a powerful framework for understanding how complex structures emerge in a variety of physical, chemical, and biological systems. Its applications span numerous fields, making it a versatile and essential tool in both theoretical and applied research.

5. Conclusions

This study focused on developing efficient numerical methods utilizing fractional difference schemes in the Caputo sense. The stability and convergence of these schemes were thoroughly analyzed. They were applied to solve the Allen-Cahn equation, the KPP-Fisher equation, and the complex Ginzburg-Landau equation across one, two, and three dimensions. These examples highlight the versatility and importance of reaction-diffusion models in various scientific and engineering disciplines. The ability to model and predict the behavior of complex systems makes these models invaluable tools for research and practical applications. The numerical experiments conducted revealed a variety of spatiotemporal patterns, showcasing their potential applications in physics, biology, and engineering processes. Future research will extend the presented methodology to address coupled real-life problems, enhancing its practical applicability.