
==== Front
Sci Rep
Sci Rep
Scientific Reports
2045-2322
Nature Publishing Group UK London

38555293
58073
10.1038/s41598-024-58073-z
Article
Dynamics of the time-fractional reaction–diffusion coupled equations in biological and chemical processes
Ghafoor Abdul abdulghafoor@kust.edu.pk

1
Fiaz Muhammad 1
Hussain Manzoor 2
Ullah Asad asad@ujs.edu.cn

34
Ismail Emad A. A. 5
Awwad Fuad A. 5
1 https://ror.org/057d2v504 grid.411112.6 0000 0000 8755 7717 Institute of Numerical Sciences, Kohat University of Science and Technology, Kohat, 26000 KP Pakistan
2 https://ror.org/00pb7yd26 Department of Mathematics, Faculty of Sciences and Technology, Women University of Azad Jammu and Kashmir, Bagh, Azad Kashmir Pakistan
3 https://ror.org/03jc41j30 grid.440785.a 0000 0001 0743 511X School of Finance and Economics, Jiangsu University, 301, Xuefu Road, Jingkou District, Zhenjiang, 212013 Jiangsu China
4 grid.513214.0 Department of Mathematical Sciences, University of Lakki Marwat, Lakki Marwat, 28420 Khyber Pakhtunkhwa Pakistan
5 grid.56302.32 0000 0004 1773 5396 Department of Quantitative Analysis, College of Business Administration, King Saud University, P.O. Box 71115, Riyadh, 11587 Saudi Arabia
30 3 2024
30 3 2024
2024
14 754931 12 2023
25 3 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
This paper aims to demonstrate a numerical strategy via finite difference formulations for time fractional reaction–diffusion models which are ubiquitous in chemical and biological phenomena. The time-fractional derivative is considered in the Caputo sense for both linear and nonlinear problems. First, the Caputo derivative is replaced with a quadrature formula, then an implicit method is used for the remaining part. In the linear case, the proposed strategy reduces the time fractional models into linear simultaneous equations. In nonlinear cases, Quasilinearization is utilized to tackle the nonlinear parts. With this strategy, solutions of the fractional system transform into linear algebraic systems which are easy to solve. Next, the Von Neumann method is implemented to examine the stability of the scheme which discloses that the scheme is unconditionally stable. Further, the applicability of the presented scheme is tested with different linear and nonlinear models which include the one dimensional Schnakenberg and Gray–Scott models, and one and two dimensional Brusselator models. To analyze the accuracy of the present technique two norms namely, L∞ and L2, and relative error are addressed. Moreover, the obtained outcomes are shown tabulated and graphically which identifies that the scheme properly works for the time fractional reaction–diffusion systems.

Keywords

Fractional calculus
Implicit scheme
Caputo fractional derivative
Brusselator model
Schnakenberg model
Gray–Scott model
Stability analysis
Subject terms

Energy science and technology
Mathematics and computing
issue-copyright-statement© Springer Nature Limited 2024
==== Body
pmcIntroduction

Reaction–diffusion models (RDMs) play a vital role in describing various spatial patterns like mazes, stripes, and spots through chemical operations in cells. RDMs theory has been started from the pioneer work of Turing1 which explored the importance of pattern formation via RDMs and biological processes. Particularly, this theory describes that uniform stability of the system remains in the absence of diffusion parameters while different spatial pattern formations can be realized in the presence of reaction and diffusion. Many authors identified the usage of RDMs models in various scientific and engineering disciplines. For example pattern formation in hydra2, shell pigmentation3, animal coat markings4 and many other for which the readers may refer to see5. The aforementioned applications show that RDMs are ubiquitous in different areas of science.

RDMs are highly non-linear and its closed form solution is a challenging task. Therefore, numerical techniques are the alternative remedies to capture the dynamics of non-linear models. Several computational strategies have been advised in the literature to determine the numerical solutions of non-linear RDMs related to pattern formation. For example, Ersoy6 established a computational algorithm for the study of RDMs using an exponential cubic B-spline. Onarcan et al.7 proposed a numerical based on trigonometric cubic B-spline to solve RDMs. Similarly, finite difference-based techniques8,9 and finite element method10 have been used to solve the RDMs. Mittal and his co-author developed solved RDMs by modified cubic B-spline coupled with differential quadrature. Korkmaz et al.11 investigated the motion of different patterns modeled by a special case of RDMs.

Different RDMs are presented in available studies. For example, Brusselator model which is also known as tri-molecular chemical reaction system which demonstrates the existence of Turing instability in the autocatalytic reactions12. Another famous system is the Gray-Scott model which delineates spot-like patterns that remain self-reproducing structures in the whole domain11. Besides, RDMs consist of the well-known Schnakenberg model also known as activator-depleted model. Physically, the Schnakenberg model elucidates a chemical reaction between the chemical substances namely: activator and an inhibitor. The fast disappearance of inhibitor than activator leads to Turing instability.

The study of the aforementioned models was limited to integer calculus which means that all RDMs involve integer order derivatives. The fractional version of these models is still an open area of research. In this case, the order of the involved derivatives is arbitrary which lies in well known branch of mathematics known as factional calculus(FC)13. The theory of FC was commenced when Leibnitz asked L’Hôpital about the half derivative and then the idea was extended in different directions by renowned mathematicians like Riemann, Abel, Laplace, and Euler. Recently many applications of fractional partial differential equations (FPDEs) have been observed in different fields of science and engineering such as bioengineering14, solid mechanics, nonlinear oscillations of an earthquake, fluid-dynamic traffic model15, economics, and anomalous transport16. Fractional models dominate over classical models in the sense of memory effect which means that the next stat of the dynamical system will equally defend on the preceding and present stats. FPDEs are subdivided into three categories, those having space fractional derivatives are known as space fractional models. Similarly, those having time fractional derivatives and combined space and time fractional derivatives are called time-fractional and space-time fractional models. Further, the FPDEs are divided on the basis of order which is either constant or a function of space and time variables which are called consultant and variable order FPDEs, respectively. In the last two decades, the time fractional partial differential equations (TFPDEs) attracted a lot of researchers due to the time evolution in fractions.

Numerous approaches have been discussed in existing work to solve FPDEs and TFPDEs. Such as, Chen et al.17 presented fractional percolation problem via implicit finite difference method. Huang and Zhao18 implemented two distinct alternating direction implicit numerical stratagems for the numerical estimation of linear non-linear super diffusion equations. In19 finite element methods (FEMs) have been suggested to solve TFPDEs like fractional advection–dispersion systems, and Fokker–Planck model20 having fractional time and space derivatives. Like finite difference techniques and FEMs, Meshfree numerical procedures have been implemented for FPDEs and TFPDEs. For example, Uddin21 explored radial basis functions (RBFs) meshfree strategy for the numerical simulations of TFPDEs. Hussain et al.22 solved the time-fractional coupled KdV equations via meshless spectral numerical technique. Aghdam et al.23 utilized shifted Legendre polynomials and proposed a computational strategy for the numerical estimation of fractional advection–diffusion equations. Similarly, in24 the authors developed an efficient method to simulate the fractional Black-Scholes model for European options via Gegenbauer polynomials and Caputo derivatives. Moreover, a higher order numerical scheme based on quadratic interpolation for time and Chebyshev collocation for space was addressed for the solutions of space-time fractional advection–diffusion equation25. Mesgarani et al.26 proposed a numerical scheme for the space fractional advection–diffusion equation b coupling second-order accurate difference approach and second kind shifted Chebyshev polynomials.

Besides, we address some more recent work related to pattern formation. Such as Owolabi et al.27 proposed Fourier spectral method for the study of unveil complex Turing patterns arising in autocatalytic reactions via the Caputo fractional derivative. Alqhtani et al.28 investigated spatiotemporal and chaotic patterns with fractional order prey-predator models. Alqhtani et al.29 studied the Caputo fractional derivatives in predator-prey models to explore, the formation of spatiotemporal patterns through local and global stability analysis. The authors in30 presented different pattern formation scenarios of various superdiffusive system involving Caputo operator. Owolabi et al.31 developed a mathematical model using fractional-order super-diffusive processes for some emergent pattern formation in predator-prey system. Owolabi and Baleanu32 described several emergent patterns for diffusive turing-like models with fractional operator. Similarly in33 the authors studied the spatial patterns of modified prey-predator via diffusion-driven instability with associated chaotic behaviors. Moreover, Owolabi and Patidar34 discussed the higher order solutions strategies for the stiff time-dependent PDEs and its spatiotemporal dynamics of reaction–diffusion systems. Owolabi35 explored the pattern formation of different fractional reaction–diffusion systems. Pindza and Owolabi36 implemented a numerical scheme for the space fractional reaction–diffusion equations using Fourier spectral strategy for space and exponential integrator for time.

The current paper addresses the numerical solutions of time-fractional reactional models TFRDMs. There is sufficient space to study the dynamics of the TFRDMs. Our target is to consider these models in fractional form and explain the solutions numerically. The next goal is to elaborate on the stability analysis of the system. For verification of the proposed scheme, some test problems will be addressed. In those problems having an exact solution will be treated as artificial solutions for fractional cases. Moreover, three types of boundary conditions will be encountered.

The leftover paper is prescribed the following way. In “Basic definitions” some preliminaries related to FC, are reported. The proposed methodology and the boundary conditions are presented in “Description of the method for one-dimensional TFRDMs”. The stability of the scheme is given in “Stability analysis”. Algorithm for the proposed method is presented in “Algorithm”. Finally, the test paradigms and the conclusions are drawn in “Numerical experiments” and “Description of the method for two-dimensional TFRDMs”, respectively.

Basic definitions

In this section of the manuscript, some basic definitions related to FC are discussed.

Definition

The fractional derivative of χ in Caputo sense37 is defined as follow:1 cDtαχ(x,t)=∂αχ(x,t)∂tα=1Γ(1-α)∫0t∂χ(x,s)∂s(t-s)-αds,0<α<1∂χ∂t,α=1.

Quadrature rule

The Caputo fractional derivative of χ can be approximated by the quadrature formula given by21:2 ∂αχ(x,tn+1)∂tα=Aα∑k=0nβkα(χn-k+1-χn-k)+O(τ2-α),if0<α<1,χn+1-χnτ+O(τ)ifα=1,

Where χ(x,tn)=χn, Aα=τ-αΓ(2-α); τ is time step and βkα=(k+1)1-α-(k)1-α.

Definition

The one and two parameters Mittag-Leffler function is defined as follows37:3 Eη(z)=∑k=0∞zkΓ(ηk+1)η>0,

4 Eη,ζ(z)=∑k=0∞zkΓ(ηk+ζ)η,ζ>0.

The Caputo derivative of exponential and trigonometric functions can be expressed in the form of the Mittag-Leffler function which is described below:0CDxα(eλx)=λnxn-αE1,n-α+1(λx), where n=⌈α⌉

0CDxαsin(λx)=12i((iλ)nxn-αE1,n-α+1(iλx)-(iλ)nxn-αE1,n-α+1(iλx)), where n=⌈α⌉

0CDxαcos(λx)=12((iλ)nxn-αE1,n-α+1(iλx)+(iλ)nxn-αE1,n-α+1(iλx)) , where n=⌈α⌉.

Description of the method for one-dimensional TFRDMs

Consider the following TFRDM38–41:5 ∂αU(x,t)∂tα=a1Uxx+b1U+c1V+d1U2V+e1UV+m1UV2+κ1+F1(x,t),x∈Ω,t>0∂αV(x,t)∂tα=a2Vxx+b2U+c2V+d2U2V+e2UV+m2UV2+κ2+F2(x,t),x∈Ω,t>0,

In Eq. (5) U(x,t) and V(x,t) are two dependent variables describes the dynamics where 0<α≤1, aj,bj,cj,dj,ej,mj and κj, are real constants for each j=1,2 which described some physical interpretation discussed in problem 5.1. F1 and F2 are the source terms and Ω=[x0,xN] is the spacial domain is divided into M subintervals where width of each subinterval is dx=xN-x0M. The associated initial conditions (ICs) are:6 U(x,0)=U0(x),V(x,0)=V0(x)x∈Ω.

The corresponding boundary conditions (BCs) are categorized in the following way:Type 1:7 U(x0,t)=ε0(t),U(xN,t)=δ0(t),t>0V(x0,t)=ε1(t),V(xN,t)=δ1(t),t>0.

Type 2: 8 Ux(x0,t)=ε0(t),U(xN,t)=δ0(t),t>0Vx(x0,t)=ε1(t),V(xN,t)=δ1(t),t>0.

Type 3: 9 Ux(x0,t)=ε0(t),Ux(xN,t)=δ0(t),t>0Vx(x0,t)=ε1(t),Vx(xN,t)=δ1(t),t>0.

Now, based on the variation of the coefficient, Eq. (5) reduces to a variety of different linear and non-linear models. The models which will be under investigation, are listed in Table 138–41:Table 1 Coefficients of different RDMs.

Problems	a1	a2	b1	b2	c1	c2	d1	d2	e1	e2	m1	m2	κ1	κ2	
Linear	d	d	-a	0	1	-b	0	0	0	0	0	0	0	0	
Brusselator	ε1	ε2	-B-1	B	0	0	1	− 1	0	0	0	0	A	0	
Schnakenberg	1	d	-γ	0	0	0	γ	-γ	0	0	0	0	γa	γb	
Gray-Scott	ε1	ε2	-f	0	0	-f-k	0	0	0	0	− 1	1	f	0	

Using Eq. (2) and an implicit scheme to Eq. (5) the resultant is:10 Aα∑k=0nβkαUn-k+1-Un-k=a1∂2U∂x2n+1+b1Un+1+c1Vn+1+d1U2Vn+1+e1UVn+1+m1UV2n+1+κ1+F1n+1,

11 Aα∑k=0nβkαVn-k+1-Vn-k=a2∂2V∂x2n+1+b2Un+1+c2Vn+1+d2U2Vn+1+e2UVn+1+m2UV2n+1+κ2+F2n+1,

where Un=U(x,tn), Vn=V(x,tn), F1n=F1(x,tn) and F2n=F2(x,tn)

Further simplification of the above Eqs. (10–11) leads to:12 a1∂2U∂x2n+1+b1Un+1+c1Vn+1+d1U2Vn+1+e1UVn+1+m1UV2n+1-AαUn+1=-AαUn-F1n+1-κ1+Aα∑k=1nβkαUn-k+1-Un-k,

13 a2∂2V∂x2n+1+b2Un+1+c2Vn+1+d2U2Vn+1+e2UVn+1+m2UV2n+1-AαVn+1=-AαVn-κ2-F2n+1+Aα∑k=1nβkαVn-k+1-Vn-k.

The nonlinear terms U2Vn+1, UVn+1, and UV2n+1 are linearized using Qausilinearization42:14 U2Vn+1=2UnVnUn+1-2U2nVn+U2nVn+1

15 UVn+1=VnUn+1-UVn+UnVn+1

16 UV2n+1=2UnVnVn+1-2UnV2n+V2nUn+1

Inserting the central difference approximation of Uxxn+1 and Vxxn+1 together with non-linear terms Eqs. (14–16), in Eqs. (12) and (13) eventually gives linear system of equations with the compact form given by:17 AWn+1=Wn,

where the dimensions of the matrices of A and W are described below:For type 1 boundary conditions, the order of A and W will be (2M-2)×(2M-2) and (2M-2)×1, respectively.

For type 2 boundary conditions, the order of A and W will be 2M×2M and 2M×1, respectively.

For type 3 boundary conditions, the order of A and W will be (2M+2)×(2M+2) and (2M+2)×1, respectively.

Stability analysis

In this part of the manuscript, the stability analysis of the numerical method for the TFRDMs model is discussed. Assume the coefficients for this model from the first row of Table 1, then the system reduces to:18 dλUj-1n+1-(2dλ+a+Aα)Ujn+1+dλUj+1n+1+Vjn+1=-AαUjn+Aα∑k=1nβkαUjn-k+1-Ujn-k

19 dλVj-1n+1-(2dλ+b+Aα)Vjn+1+dλVj+1n+1=-AαVjn+Aα∑k=1nβkαVjn-k+1-Vjn-k.

Theorem

The implicit scheme (18-19), for the system (5), with α∈(0,1) on x∈[x0,xN] is unconditionally stable.

Proof

Following the procedure43 we assume:20 Ujn=ξneiωjh,Vjn=ηneiωjh,

where i=-1. Using Eq. (20) in Eqs. (18–19) and some algebraic manipulation leads to:21 ξn+12dλAα(cos(ωh)-1)-aAα-1+ηn+1Aα=-ξn+∑k=1nβkα(ξn-k+1-ξn-k),

22 ηn+12dλAα(cos(ωh)-1)-bAα-1=-ηn+∑k=1nβkα(ηn-k+1-ηn-k),

From Eq. (22) we get:23 ηn+1=ηn+∑k=1nβkα(ηn-k-ηn-k+1)2dλAα(1-cos(ωh))+bAα+1,

Since 2dλAα(1-cos(ωh)+bAα+1≥1, hence it follows:24 ηn+1≤ηn+∑k=1nβkα(ηn-k-ηn-k+1).

For n=025 η1≤η0.

In a similar fashion the following hold:η2≤η1+β1α(η0-η1).

If two consecutive approximations are closed then their difference approaches zero hence:η2≤η1.

Continuing in this way one can write26 ηn+1≤ηn≤⋯≤η0.

Now putting Eq. (23) in Eq. (21) we obtain:27 ξn+12dλAα(cos(ωh)-1)-aAα-1+ηn+∑k=1nβkα(ηn-k-ηn-k+1)2dλAα(1-cos(ωh))+bAα+1Aα=-ξn+∑k=1nβkα(ξn-k+1-ξn-k),

which gives:28 ξn+1≤ξn+∑k=1nβkα(ξn-k-ξn-k+1)2dλAα(1-cos(ωh))+aAα+1.

From earlier discussion, the following result can be deduced:29 ξn+1≤ξn≤⋯≤ξ0.

Thus, ξn+1=|Ujn+1|≤ξ0=|Uj0|=|fj| and ηn+1=|Vjn+1|≤ξ0=|Vj0|=|fj|. These results imply that ||Un||l2≤||f||l2, ||Vn||l2≤||f||l2 which is the stability condition.

Algorithm

Input: 0<α≤1, aj,bj,cj,dj,ej,mj and κj for j=1,2 and time step size τ.

Output:Solve the system of PDEs numerically using FDM to evaluate approximate solutions Un and Vn.

Step 1: Generate computational grid over the domain Ω, discretizing the spatial dimensions into a finite set of points.

Step 2: Apply central differences formula to discretize the spatial derivatives in the PDEs and quadrature formula for time derivative.

Step 3: Set n=0.

Step 4: Calculate matrices A and Wn.

Step 5:Calculate Wn+1 by using inversion method from AWn+1=Wn,.

Step 6:Start time loop n=1:N

Step 7:Repeat the step 4-5.

Numerical experiments

In this section, the proposed scheme is applied to some linear and non-linear RDMs. To assure the accuracy, global relative error (RE), L2 and L∞ have been used which are defined below:30 RE=∑j=1N|χjn+1-χjn|2|χjn+1|12,

where χjn+1 and χjn are the approximate solutions at two consecutive time levels.31 L2=∑j=1Nχjext-χjapp212,

32 L∞=max1≤j≤Nχjext-χjapp,

where χjext and χjapp are the exact and approximate solutions, respectively.

Problem 5.1 (linear model)

Consider the coefficients from the first row of Table 1 the resultant equation is given by38–41:33 ∂αU(x,t)∂tα=d∂2U(x,t)∂x2-aU(x,t)+V(x,t)+F1∂αV(x,t)∂tα=d∂2V(x,t)∂x2-bV(x,t)+F2.

The corresponding ICs and BCs are:34 U(x,0)=2cos(x),0≤x≤π/2,V(x,0)=(a-b)cos(x),0≤x≤π/2,

35 Ux(0,t)=0,U(π/2,t)=0,0<t≤1,Vx(0,t)=0,V(π/2,t)=0,0<t≤1.

The closed-form solutions of this system are:U(x,t)=e-(a+d)t+e-(b+d)tcosx,V(x,t)=(a-b)e-(b+d)tcosx.

Using the fact (0CDxα(eλx)=λnxn-αE1,n-α+1(λx)) the associated source terms are extracted as follows:F1=-cos(x)t1-α[(a+d)E1,2-α(-(a+d)t)+(b+d)E1,2-α(-(b+d)t)]+cos(x)[(a+d)e-(a+d)t+(b+d)e-(b+d)t]F2=-(a-b)(b+d)cos(x)t1-α(a+d)E1,2-α(-(b+d)t)+(a-b)(b+d)cos(x)e-(b+d)t.

For numerical simulations three cases are addressed here.

Diffusion-dominated

For diffusion-dominated case the parameters considered are a = 0.1, b = 0.01 and d = 1.

Reaction-dominated

For the reaction-dominated case the parameters considered are a = 2, b = 1 and d = 0.001

Reaction-dominated with stiff reaction

For this case the selection of parameters are a = 100, b = 1 and d = 0.001.

This problem has been solved numerically, with the suggested technique, and the obtained results are noted in the form of tabulated and graphical forms. In Tables 2, 3, 4, 5, 6 and 7, the L2 and L∞ error norms are recorded for different times and α. The consecutive Tables 2 and 3 show the results for diffusion-dominated case, Tables 4 and 5 for reaction-dominated case and Tables 6 and 7 for reaction-dominated with stiff reaction case, respectively. Tabulated simulations reveal the good performance of the present technique. Similarly, solutions profile of exact versus numerical solutions are displayed in Figs. 1, 2 and 3 for diffusion dominated, reaction dominated and reaction dominated with stiff reaction using α=0.5. From the graphical results it is plainly visible, that both solutions are in good agreement.Table 2 L∞ norms for distinct values of α at different time levels in the diffusion-dominated when a = 0.1, b = 0.01, d = 1.

t	α=0.25	α=0.5	α=0.75	α=1	
L∞ of U	L∞ of V	L∞ of U	L∞ of V	L∞ of U	L∞ of V	L∞ of U	L∞ of V	
0.001	6.0820e−03	2.7371e−04	6.0245e−03	2.7112e−04	5.9195e−03	2.6639e−04	5.7664e−03	2.5950e−04	
0.01	6.0349e−03	2.7170e−04	6.0126e−03	2.7069e−04	5.9832e−03	2.6937e−04	5.9586e−03	2.6826e−04	
0.05	5.7908e−03	2.6118e−04	5.7825e−03	2.6080e−04	5.7755e−03	2.6048e−04	5.7753e−03	2.6047e−04	
0.075	5.6416e−03	2.5473e−04	5.6360e−03	2.5448e−04	5.6328e−03	2.5433e−04	5.6321e−03	2.5447e−04	
0.1	5.4960e−03	2.4844e−04	5.4924e−03	2.4827e−04	5.4916e−03	2.4823e−04	5.4914e−03	2.4818e−04	

Table 3 L2 norms for distinct values of α at different time levels in the diffusion-dominated when a = 0.1, b = 0.01, d = 1.

t	α=0.25	α=0.5	α=0.75	α=1	
L2 of U	L2 of V	L2 of U	L2 of V	L2 of U	L2 of V	L2 of U	L2 of V	
0.001	5.0461e−02	2.2745e−03	3.3379e−02	1.5027e−03	2.2302e−02	1.0037e−03	1.4967e−02	6.7357e−04	
0.01	5.5703e−02	2.5131e−03	4.4871e−02	2.0224e−03	3.4863e−02	1.5700e−03	2.7260e−02	1.2273e−03	
0.05	5.6607e−02	2.5591e−03	5.1516e−02	2.3267e−03	4.6751e−02	2.1082e−03	4.1236e−02	1.8538e−03	
0.075	5.5951e−02	2.5315e−03	5.2535e−02	2.3720e−03	5.0255e−02	2.2612e−03	4.7814e−02	2.1335e−03	
0.1	5.5298e−02	2.5025e−03	5.3686e−02	2.4178e−03	5.3054e−02	2.4174e−03	5.2420e−02	2.2854e−03	

Table 4 L∞ norms for distinct values of α at different time levels in the reaction-dominated when when a = 2, b = 1, d = 0.001.

t	α=0.25	α=0.5	α=0.75	α=1	
L∞ of U	L∞ of V	L∞ of U	L∞ of V	L∞ of U	L∞ of V	L∞ of U	L∞ of V	
0.001	4.6895e−03	2.3665e−03	3.5473e−03	1.7817e−03	1.9870e−03	9.9530e−04	5.7697e−04	2.8865e−04	
0.01	4.8733e−03	2.4754e−03	4.3649e−03	2.2087e−03	3.7403e−03	1.8860e−03	2.8101e−03	1.4127e−03	
0.05	4.7159e−03	2.4461e−03	4.5045e−03	2.3300e−03	4.3338e−03	2.2334e−03	4.2020e−03	2.1545e−03	
0.075	4.5710e−03	2.4002e−03	4.4164e−03	2.3132e−03	4.3221e−03	2.2553e−03	4.3050e−03	2.2335e−03	
0.1	4.4232e−03	2.3505e−03	4.3049e−03	2.2822e−03	4.2558e−03	2.2476e−03	4.3001e−03	2.2567e−03	

Table 5 L2 norms for distinct values of α at different time levels in the reaction-dominated when when a = 2, b = 1, d = 0.001.

t	α=0.25	α=0.5	α=0.75	α=1	
L2 of U	L2 of V	L2 of U	L2 of V	L2 of U	L2 of V	L2 of U	L2 of V	
0.001	7.3047e−03	3.7330e−03	4.3563e−03	2.1923e−03	2.0693e−03	1.0366e−03	5.7805e−04	2.8918e−04	
0.01	8.2597e−03	4.2811e−03	6.3431e−03	3.2291e−03	4.6695e−03	2.3454e−03	3.1145e−03	1.5238e−03	
0.05	9.2742e−03	4.5919e−03	1.0909e−02	4.2832e−03	1.3506e−02	4.2313e−03	1.8929e−02	4.7582e−03	
0.075	1.1636e−02	4.9199e−03	1.8456e−02	5.5718e−03	2.5948e−02	6.5884e−03	3.9181e−02	8.8837e−03	
0.1	1.6020e−02	5.6583e−03	2.9267e−02	7.7966e−03	4.2602e−02	1.0189e−02	6.5769e−02	1.4717e−02	

Table 6 L∞ norms for distinct values of α at different time levels in the reaction-dominated with stiff when when a = 100, b = 1, d = 0.001.

t	α=0.25	α=0.5	α=0.75	α=1	
L∞ of U	L∞ of V	L∞ of U	L∞ of V	L∞ of U	L∞ of V	L∞ of U	L∞ of V	
0.001	3.2176e−03	2.9633e−01	8.6635e−04	6.7405e−02	6.0743e−04	1.1838e−02	1.4373e−04	2.0600e−03	
0.01	4.3451e−03	4.1886e−01	2.6476e−03	1.8405e−01	1.5799e−03	6.2873e−02	8.0842e−04	2.0158e−02	
0.05	5.0192e−03	4.9547e−01	4.5914e−03	3.4635e−01	3.4836e−03	2.0612e−01	2.8933e−03	1.0922e−01	
0.075	5.1232e−03	5.0579e−01	4.6426e−03	3.8353e−01	4.1662e−03	2.5651e−01	2.5435e−03	1.5425e−01	
0.1	5.2140e−03	5.1566e−01	5.0360e−03	3.8781e−01	4.2630e−03	2.6075e−01	4.1022e−03	1.6462e−01	

Table 7 L2 norms for distinct values of α at different time levels in the reaction-dominated with stiff when when a = 100, b = 1, d = 0.001.

t	α=0.25	α=0.5	α=0.75	α=1	
L2 of U	L2 of V	L2 of U	L2 of V	L2 of U	L2 of V	L2 of U	L2 of V	
0.001	3.2467e−03	2.9931e−01	1.0402e−03	6.7622e−02	3.0531e−03	1.1873e−02	7.2142e−04	2.0664e−03	
0.01	4.6633e−03	4.2629e−01	4.3442e−03	1.8500e−01	4.0301e−03	6.3150e−02	8.506e−04	2.1582e−02	
0.05	1.1079e−02	5.0874e−01	1.0922e−02	3.5317e−01	6.8222e−03	2.2716e−01	8.4562e−03	2.0929e−01	
0.075	1.2536e−02	5.2131e−01	1.1626e−02	4.0113e−01	1.0644e−02	3.2946e−01	9.2769e−03	2.1055e−01	
0.1	6.9750e−02	5.3264e−01	1.6586e−02	4.2573e−01	1.3483e−02	4.5355e−01	1.0593e−02	7.3185e−01	

Figure 1 Solution profile with absolute error of U and V at t=0.001 for diffusion dominant case.

Figure 2 Solutions profile of U and V at t=0.001 for reaction dominant case.

Figure 3 Solutions profile of U and V at t=0.001 for reaction dominant with stiff reaction.

Problem 5.2 (Brusselator model)

The exact solution this model and the next three nonlinear models is not available. For these cases ignored the corresponding source term. The numerical results in non-fractional form are given in literature. We tried the proposed technique for different values of α and observed its graphical behavior which approaches towards the integer value of α=1. In these cases the relative error is measured between two consecutive time levels. On behalf of these arguments we claim that proposed method works for such nonlinear fractional models. Using the coefficients from the second row of Table 1, which gives the following Brusselator model of kinetics for two chemical components38–41:36 ∂αU(x,t)∂tα=ε1Uxx-B+1U+U2V+A,∂αV(x,t)∂tα=ε2Vxx+BU-U2V,

with ICs:37 U(x,0)=0.5,V(x,0)=1+5x+14tanh(20x)-14tanh(20(x-1)),

and the BCs are:38 Ux(0,t)=0,U(1,t)=0,Vx(0,t)=0,V(1,t)=0.

Simulations are performed for the parameters: ε1=ε2=0.00002,A=1,B=3.4, used in41 over the problem domain [0,1] . The relative error values for U and V at various values of α and different time levels are presented in Table 8. Similarly, the density values for periodical motion are given in Table 9. Tabulated data discloses that the proposed technique gives good results in terms of error values which are comparable with the previous results in literature. Graphically, the numerical solutions for integer derivative and times t=3,6,10.7,13.7 are shown in Fig. 4, which predict the oscillatory behavior in chemical reactions. In Fig. 5 the results are calculated for various values of fractional derivative, which retain the structure of oscillations, only the peaks are different due to fractional values of α. From a graphical view two things are clear, the results agree with classical solutions and the scheme does not alter the physical meaning of the model with fractional derivative.Table 8 Relative error for various α values at various time levels when dt=0.01 and dx=0.05 of the Brusselator model.

t	α=0.25	α=0.5	α=0.75	α=1	
RE of U	RE of V	RE of U	RE of V	RE of U	RE of V	RE of U	RE of V	
3	2.0354e−10	2.2025e−10	3.7863e−04	3.5430e−03	1.5726e−02	7.0784e−03	1.7728e−02	7.0979e−03	
6	9.9154e−15	4.4068e−15	2.8500e−04	1.6795e−03	2.2394e−02	3.7622e−03	3.4562e−04	8.2614e−04	
10	7.5003e−16	2.0311e−15	2.1909e−02	6.8632e−03	2.1484e−02	1.9623e−03	1.9945e−02	7.7037e−03	
13	6.0902e−15	6.8191e−16	4.1097e−04	3.7381e−03	6.2843e−03	1.6586e−03	1.0856e−02	4.6599e−03	

Table 9 Density value for periodic motion of Brusselator model.

Density	x	0	0.2	0.4	0.6	0.8	1	
U	3	0.294851285042341	0.342874476982005	0.439887206799160	4.084153286206512	0.868264657607814	0.631669227084894	
6	0.455133646517029	4.637328613106386	1.362737460708937	0.316810767726646	0.341809435938578	0.353024226227808	
10.7	0.315507295348274	0.340444936787278	0.426904326384855	4.422181773179734	0.975057358624334	0.703916805055259	
13.7	0.453599334757512	3.647001834026506	1.500682343948168	0.319195589894421	0.338228804864134	0.348914444829436	
V	3	3.687773918047040	4.682788687309119	5.410800581895748	0.844858452761222	2.299626799142092	2.604506569694734	
6	5.499253473873485	1.150072687778154	1.840424829740125	3.662292719056844	4.649894995801780	4.809935414837560	
10.7	3.680063713497682	4.627752556692187	5.392449423645793	0.770900658236533	2.184852828836223	2.501560571021531	
13.7	5.491790039896567	2.426212381897124	1.738222712948466	3.605397600513641	4.591809581911599	4.754792117638447	

Figure 4 Numerical illustration of Brusselator model for U and V when α=1..

Figure 5 Numerical illustration of Brusselator model for U and V when α=0.75 and α=0.9.

Problem 5.3 (Gray–Scott model)

Here, we take the values of coefficients from the third row of Table 1 which gives the Gray-Scot model38–4139 ∂αU(x,t)∂tα=ε1Uxx-UV2+f1-U,∂αV(x,t)∂tα=ε2Vxx+UV2-f+kV.

The ICs are:40 U(x,0)=1-12sin100πx-L2LV(x,0)=14sin100πx-L2L,

and the Dirichlet type of BCs are given below:41 U(-L,t)=1,U(L,t)=1V(-L,t)=1,V(L,t)=1.

Simulations are carried out for the spatial domain [-L,L], with parameter values L=50, ε1=1,ε2=0.01, f=0.02, k=0.066, which are taken from41. The achieved error values for differena longues of α and time are presented in Table 10. Tabulated simulation reveals that the scheme works well for large time. Numerical solutions are plotted in Fig. 6 for integer cases which show that the outcome is in good agreement with available results in previous work. Likewise, in Fig. 7 the simulations are noted for fractional cases which indicate that the graphical behavior is nearly similar to integer when α=0.9. In Fig. 8 three dimensions numerical solutions are shown which show the wave type pattern.Figure 6 Numerical illustration for Gray–Scott model for U and V when α=1..

Table 10 Relative error for various α values at various time levels when dt=100 and dx=1 of the Gray-Scott model.

t	α=0.25	α=0.5	α=0.75	α=1	
RE of U	RE of V	RE of U	RE of V	RE of U	RE of V	RE of U	RE of V	
2000	1.5586e−03	2.4029e−03	2.5015e−03	1.7205e−03	5.3181e−07	1.3595e−06	1.7787e−09	1.3074e−09	
2400	2.2024e−04	1.5784e−04	6.6420e−04	3.6827e−04	7.6539e−08	1.2872e−07	2.6013e−11	1.9011e−11	
2800	3.8928e−05	2.7964e−05	1.9440e−04	7.4461e−05	1.1049e−08	1.2020e−08	3.3840e−13	4.1099e−13	
4000	8.5794e−07	1.7319e−06	6.2215e−06	1.2775e−06	3.3900e−11	7.6449e−12	1.0078e−13	1.0261e−13	

Figure 7 Numerical illustration for Gray–Scott model in 2D Graph when dt=0.5, dx=0.25 α=0.9.

Figure 8 Numerical illustration of Gray–Scott model in 3D Graph when dt=0.5, dx=0.25 α=0.9..

Problem 5.4 (Schnakenberg model)

Choosing the coefficient values from the fourth row of Table 1, we obtain the following Schnakenberg model38–41:42 ∂αU(x,t)∂tα=d1Uxx+γa1-U+U2V,∂αV(x,t)∂tα=d2Vxx+γb1-U2V,

where U and V denote the activator and inhibitor concentrations, respectively, and d is diffusion coefficients, γ, a, and b are biological reaction rate constants. The following are the associated ICs:43 U(x,0)=0.919145+11000∑j=125cos(2πjx)jV(x,0)=0.937903+11000∑j=125cos(2πjx)j,

and the BCs are:44 Ux(0,t)=0,Ux(1,t)=0Vx(0,t)=0,Vx(1,t)=0.

To solve this model, we use a=0.126779,b=0.792366,d=10, and γ=5000,10000. The computed simulations in the form of relative error are noted in Table 11. From tabulated values, it is obvious that computed solutions are pretty much accurate. The one-dimensional solution profiles are plotted in Fig. 9 for integer and fractional values of α which provide a clear picture of oscillatory motion. Moreover, the solutions in the three-dimensional form are presented in Fig. 10.Table 11 Relative error for various α values at various time levels when t=2.5 and dx=0.01 of the Schnakenberg model when γ=10,000.

dt	α=0.25	α=0.5	α=0.75	α=1	
RE of U	RE of V	RE of U	RE of V	RE of U	RE of V	RE of U	RE of V	
0.1	5.1715e−15	4.7112e−15	1.4203e−15	1.1778e−15	2.3672e−15	2.3556e−15	2.5858e−15	6.5671e−15	
0.01	4.0789e−15	3.4977e−15	2.9864e−15	3.7118e−15	3.0592e−15	2.1771e−15	4.9530e−15	7.2095e−15	
0.005	2.4401e−15	5.8533e−15	4.0061e−15	2.6768e−15	4.2246e−15	1.5347e−15	3.7512e−15	2.1414e−15	
0.001	3.1685e−15	1.6061e−15	1.9302e−15	1.3661e−15	4.5524e−15	6.7813e−15	3.1685e−15	5.1038e−15	

Figure 9 Numerical illustration for Schnakenberg model in one dimensional graph for U and V for γ=5000 and γ=10,000 at t=2.5 when dt=0.001, dx=0.005 and α=0.75,0.9 and 1.

Figure 10 Numerical illustration of the Schnakenberg model in 3D graph for U and V for γ=5000 and γ=10,000 at t=2.5 when dt=0.001, dx=0.005 and α=0.9 and 1.

Ethics approval and consent to participate

In this work, no materials of any other person are used. We compared the data for which the relevant reference is given.

Description of the method for two-dimensional TFRDMs

Here, we extend the proposed strategy for two-dimensional TFRDMs44,45:45 ∂αU(x,y,t)∂tα=ε1Uxx+Uyy-B+1U+U2V+A+F1(x,y,t),∂αV(x,y,t)∂tα=ε2Vxx+Vyy+BU-U2V+F2(x,y,t),

where ε1,ε2,A and B are taken as given in44,45, and F1 and F2 are the source functions to be determined via exact solutions. The spatial domain for this problem is Ω=[x0,xN]×[y0,yN], along with the following ICs and BCs:46 U(x,y,0)=U0(x,y),V(x,y,0)=V0(x,y)x,y∈Ω.

47 U(x0,y,t)=ε0(y,t),U(xN,y,t)=δ0(y,t),U(x,y0,t)=ω0(x,t),U(x,yN,t)=ξ0(x,t),t>0V(x0,y,t)=ε1(y,t),V(xN,y,t)=δ1(y,t),V(x,y0,t)=ω1(x,t),V(x,yN,t)=ξ1(x,t).

Using the stratagem discussed earlier, Eq. (45) reduces to:48 ∂αU∂tα=ε1∂2U∂x2+∂2U∂y2n+1-β+1Un+1-U2Vn+1+A+F1n+1,

49 ∂αV∂tα=ε2∂2V∂x2+∂2V∂y2n+1+βUn+1-U2Vn+1+F2n+1,

where Un=U(x,y,tn), Vn=V(x,y,tn), F1n=F1(x,y,tn) and F2n=F2(x,y,tn).

Plugging the quadrature rule for fractional derivative in Eq. (48) and Eq. (49) the new equations are given below:50 Aα∑k=0nβkαUn-k+1-Un-k=ε1∂2U∂x2+∂2U∂y2n+1-β+1Un+1-U2Vn+1+A+F1n+1,

51 Aα∑k=0nβkαVn-k+1-Vn-k=ε2∂2V∂x2+∂2V∂y2n+1+βUn+1-U2Vn+1+F2n+1.

After one term expansion Eqs. (6.6–6.7) transform to:52 AαUn+1-Un+Aα∑k=1nβkαUn-k+1-Un-k=ε1∂2U∂x2+∂2U∂y2n+1-β+1Un+1-U2Vn+1+A+F1n+1,

53 AαVn+1-Vn+Aα∑k=1nβkαVn-k+1-Vn-k=ε2∂2V∂x2+∂2V∂y2n+1+βUn+1-U2Vn+1+F2n+1.

Rearranging the terms in the above system the results are:54 ε1∂2U∂x2+∂2U∂y2n+1-β+1Un+1-U2Vn+1-AαUn+1=-Un-F1n+1-A+Aα∑k=1nβkαUn-k+1-Un-k,

55 ε2∂2V∂x2+∂2V∂y2n+1+βUn+1-U2Vn+1-AαVn+1=-Vn-F2n+1+Aα∑k=1nβkαVn-k+1-Vn-k.

Linearization of the nonlinear term U2Vn+1 is tackled by the following formula42:56 U2Vn+1=2UnVnUn+1-2U2nVn+U2nVn+1

Inserting the approximation of the nonlinear term, and the central difference approximation for the involved derivative, the linear system of equations can be obtained which are given in compact form:AWn+1=Wn,

where A=(2M-2)2×(2M-2)2 and Wn+1=(2M-2)2×1 matrices.Figure 11 Plotting of U and V at t = 5 when α=1 and dt=0.001.

Figure 12 Solution profiles of U and V at t=0.1 when α=0.5 and dt=0.001.

Problem 5.4 (2D Brusselator model)

Here, the following two-dimensional Brusselator model of kinetics for two chemical components is considered for validation of the scheme44,45:57 ∂αU(x,y,t)∂tα=ε1Uxx+Uyy-B+1U+U2V+A+F1(x,y,t),∂αV(x,y,t)∂tα=ε2Vxx+Vyy+BU-U2V+F2(x,y,t),

along with the following conditions:58 U(x,y,0)=e(-x-y)V(x,y,0)=e(x+y),x,y∈Ω.

59 U(0,y,t)=e(-y-0.5t),U(1,y,t)=e(-1-y-0.5t),U(x,0,t)=e(-x-0.5t)),U(x,1,t)=e(-x-1-0.5t),t>0V(0,y,t)=e(y+0.5t),V(1,y,t)=e(1+y+0.5t),V(x,0,t)=e(x+0.5t),V(x,1,t)=e(x+1+0.5t).

Table 12 L∞ norms for U and V at t = 2 and dt = 0.001 of Brusselator model in 2D for α=1..

dx	proposed method L∞	proposed method L2	proposed method RMS	45 L∞	46 L∞	
U	V	U	V	U	V	U	V	U	V	
1/10	2.7029e−05	9.6987e−04	1.6094e−04	5.9584e−03	1.6094e−05	5.9584e−04	2.7094e−05	1.7571e−03	7.6449e−05	3.6792e−03	
1/20	1.3277e−05	4.7984e−04	1.4357e−04	5.3227e−03	7.5564e−06	2.8014e−04	2.3714e−05	7.5561e−04	2.3469e−05	1.2749e−03	
1/25	1.0147e−05	3.4828e−04	1.3718e−04	4.8381e−03	5.7160e−06	2.0159e−04	1.3115e−05	1.4635e−03	1.5346e−05	8.3510e−04	

The associated source terms can be adjusted using the exact solution and the formula:0CDxα(eλx)=λnxn-αE1,n-α+1(λx),

where E1,n-α+1(:) is the Mittag-Leffler function defined earlier. The closed-form solution to the above problem is:U(x,y,t)=e(-x-y-0.5t),V(x,y,t)=e(x+y+0.5t).

We solve this model for the parameters ε1=ε2=0.25,A=0,B=1,dx=dt. Also, the spatial and temporal domains are [0,1]×[0,1], and [0,5], respectively. In Table 12 the obtained error norms are matched with the previous norms given in the papers for integer case44,45. It is noticed that computed outcomes are matchable with available solutions. Further, the scheme is tested for fractional values of α and the achieved results are reported in the form of L∞ and L2 norms in Tables 13 and 14, respectively. From these tables one can see, that the proposed scheme works for fractional cases as well. For further clarification, the solutions are sketched for integer and fractional cases in Figs. 11 and 12, which show that exact and numerical solutions are promised well.Table 13 L∞ norm of Brusselator model in 2D for different α and time.

t	α=0.25	α=0.5	α=0.75	α=1	
L∞ of U	L∞ of V	L∞ of U	L∞ of V	L∞ of U	L∞ of V	L∞ of U	L∞ of Vs	
0.001	9.6134e−05	3.6464e−02	1.6767e−05	6.6861e−03	6.5949e−06	2.6495e−03	1.7943e−06	7.2387e−04	
0.01	1.4283e−04	5.2265e−02	4.7277e−05	1.8750e−02	3.4592e−05	1.3754e−02	1.7389e−05	7.0160e−03	
0.05	2.1876e−04	7.7603e−02	1.2651e−04	4.9876e−02	1.3946e−04	5.5138e−02	7.6450e−05	3.0925e−02	
0.075	2.5953e−04	9.1730e−02	1.9894e−04	7.4616e−02	2.2068e−04	8.7199e−02	1.0665e−04	4.3279e−02	
0.1	3.0504e−04	1.0727e−01	2.9193e−04	1.0733e−01	3.2034e−04	1.2444e−01	1.3295e−04	5.4176e−02	

Table 14 L2 norm of Brusselator model in 2D for different α and time.

t	α=0.25	α=0.5	α=0.75	α=1	
L2 of U	L2 of V	L2 of U	L2 of V	L2 of U	L2 of V	L2 of U	L2 of V	
0.001	3.0527e−04	1.1432e−01	4.2678e−05	1.6998e−02	1.6061e−05	6.4519e−03	4.2605e−06	1.7188e−03	
0.01	5.1760e−04	1.8633e−01	1.3539e−04	5.3386e−02	9.3857e−05	3.7239e−02	4.2242e−05	1.7042e−02	
0.05	8.1083e−04	2.8673e−01	4.5437e−04	1.7421e−01	4.5566e−04	1.7771e−01	2.0426e−04	8.2508e−02	
0.075	9.7574e−04	3.4300e−01	7.3575e−04	2.7506e−01	7.6843e−04	2.9778e−01	3.0099e−04	1.2176e−01	
0.1	1.1632e−03	4.0580e−01	1.0933e−03	3.9948e−01	1.1542e−03	4.4491e−01	3.9502e−04	1.6010e−01	

Conclusions and future plan

In this work, an implicit scheme has been addressed to solve RDCMs in fractional form. The involved fraction derivative and spatial derivatives were approximated with a well-known L1 formula (Quadrature rule) and finite differences, respectively. Next, the stability of the scheme was investigated via Von Neumann analysis. Moreover, the scheme has been tested with different linear and nonlinear problems and the outcomes were compared with exact and existing results in literature. From tabulated simulations and graphical solutions, it has been observed that the proposed scheme works well for RDMs and can be used for such complicated problems having no exact solutions. In future the proposed methodology can be extended to the three dimensional problems coupling with different variable order fractional derivatives like, Caputo Fabrizio, and Atangana–Baleanu–Caputo etc. Moreover, the strategy can be tested for variable order local fractional derivative problems as well.

Acknowledgements

Researchers supporting project number (RSPD2024R1060), King Saud University, Riyadh, Saudi Arabia.

Author contributions

Conceptualization, A.G., A.U.; software, A.G., M.F., E.A.A.I.; writing—original draft preparation, F.A.A, M.F, A.G.; formal analysis, A.G., A.U, M.H.; validation, A.G, A.U., E.A.A.I.; methodology, A.G., A.U., M.F.; investigation, A.G., M.H., F.A.A. A.U.; project administration, A.G; M.F.; funding acquisition, E.A.A.I and F.A.A. All authors have read and agreed to the published version of the manuscript. The author confirms that the work described has not been published before.

Data availability

All data generated or analysed during this study are included in this published article.

Competing interests

The authors declare no competing interests.

Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Turing AM The chemical basis of morphogenesis Bull. Math. Biol. 1990 52 153 197 10.1016/S0092-8240(05)80008-4 2185858
Turing, A. M. The chemical basis of morphogenesis. Bull. Math. Biol. 52, 153–197 (1990).2185858 10.1016/S0092-8240(05)80008-4
2. Meinhardt, H.Vol. 118 (New York, 1982).
3. Meinhardt H The algorithmic beauty of sea shells springer-verlag 1995 Heidelberg New-York
Meinhardt, H. The algorithmic beauty of sea shells springer-verlag (New-York, Heidelberg, 1995).
4. Murray JD A pre-pattern formation mechanism for animal coat markings J. Theor. Biol. 1981 88 161 199 10.1016/0022-5193(81)90334-9
Murray, J. D. A pre-pattern formation mechanism for animal coat markings. J. Theor. Biol. 88, 161–199 (1981).10.1016/0022-5193(81)90334-9
5. Murray JD Mathematical biology II: spatial models and biomedical applications 2001 New York Springer
Murray, J. D. Mathematical biology II: spatial models and biomedical applications Vol. 3 (Springer, New York, 2001).
6. Ersoy, O. & Dag, I. Numerical solutions of the reaction diffusion system by using exponential cubic b-spline collocation algorithms. Open Phys. 13 (2015).
7. Onarcan AT Adar N Dag I Trigonometric cubic b-spline collocation algorithm for numerical solutions of reaction-diffusion equation systems Comput. Appl. Math. 2018 37 6848 6869 10.1007/s40314-018-0713-4
Onarcan, A. T., Adar, N. & Dag, I. Trigonometric cubic b-spline collocation algorithm for numerical solutions of reaction-diffusion equation systems. Comput. Appl. Math. 37, 6848–6869 (2018).10.1007/s40314-018-0713-4
8. Chou C-S Zhang Y-T Zhao R Nie Q Numerical methods for stiff reaction-diffusion systems Discr. Contin. Dyn. Syst. B 2007 7 515
Chou, C.-S., Zhang, Y.-T., Zhao, R. & Nie, Q. Numerical methods for stiff reaction-diffusion systems. Discr. Contin. Dyn. Syst. B 7, 515 (2007).
9. Özuğurlu E A note on the numerical approach for the reaction-diffusion problem to model the density of the tumor growth dynamics Comput. Math. Appl. 2015 69 1504 1517 10.1016/j.camwa.2015.04.018
Özuğurlu, E. A note on the numerical approach for the reaction-diffusion problem to model the density of the tumor growth dynamics. Comput. Math. Appl. 69, 1504–1517 (2015).10.1016/j.camwa.2015.04.018
10. Madzvamuse A Chung AH The bulk-surface finite element method for reaction-diffusion systems on stationary volumes Finite Elem. Anal. Des. 2016 108 9 21 10.1016/j.finel.2015.09.002
Madzvamuse, A. & Chung, A. H. The bulk-surface finite element method for reaction-diffusion systems on stationary volumes. Finite Elem. Anal. Des. 108, 9–21 (2016).10.1016/j.finel.2015.09.002
11. Korkmaz, A., Ersoy, O. & Dag, I. Motion of patterns modeled by the gray-scott autocatalysis system in one dimension. arXiv preprint arXiv:1605.09712 (2016).
12. Sahin, A. Numerical solutions of the reaction-diffusion equations with B-spline finite element method. Ph.D. thesis, Ph. D. Thesis. Turkey: Doctoral dissertation. Department of Mathematics ... (2009).
13. Podlubny, I. Fractional differential equations: An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications (Elsevier, 1998).
14. Magin, R. Fractional calculus in bioengineering, part 1. Crit. Rev. Biomed. Eng. 32 (2004).
15. He J Some applications of nonlinear fractional differential equations and their approximations Bull. Sci. Technol. 1999 15 86 90
He, J. Some applications of nonlinear fractional differential equations and their approximations. Bull. Sci. Technol. 15, 86–90 (1999).
16. Miller, K. S. & Ross, B. An introduction to the fractional calculus and fractional differential equations (Wiley, 1993).
17. Chen S Liu F Anh V A novel implicit finite difference method for the one-dimensional fractional percolation equation Num. Algorithms 2011 56 517 535 10.1007/s11075-010-9402-0
Chen, S., Liu, F. & Anh, V. A novel implicit finite difference method for the one-dimensional fractional percolation equation. Num. Algorithms 56, 517–535 (2011).10.1007/s11075-010-9402-0
18. Huang, J., Zhao, Y., Arshad, S., Li, K. & Tang, Y. Alternating direction implicit schemes for the two-dimensional time fractional nonlinear super-diffusion equations. J. Comput. Math. 37 (2019).
19. Jiang Y Ma J High-order finite element methods for time-fractional partial differential equations J. Comput. Appl. Math. 2011 235 3285 3290 10.1016/j.cam.2011.01.011
Jiang, Y. & Ma, J. High-order finite element methods for time-fractional partial differential equations. J. Comput. Appl. Math. 235, 3285–3290 (2011).10.1016/j.cam.2011.01.011
20. Deng W Finite element method for the space and time fractional fokker-planck equation SIAM J. Numer. Anal. 2009 47 204 226 10.1137/080714130
Deng, W. Finite element method for the space and time fractional fokker-planck equation. SIAM J. Numer. Anal. 47, 204–226 (2009).10.1137/080714130
21. Uddin M Haq S Rbfs approximation method for time fractional partial differential equations Commun. Nonlinear Sci. Numer. Simul. 2011 16 4208 4214 10.1016/j.cnsns.2011.03.021
Uddin, M. & Haq, S. Rbfs approximation method for time fractional partial differential equations. Commun. Nonlinear Sci. Numer. Simul. 16, 4208–4214 (2011).10.1016/j.cnsns.2011.03.021
22. Hussain M Haq S Ghafoor A Meshless spectral method for solution of time-fractional coupled kdv equations Appl. Math. Comput. 2019 341 321 334
Hussain, M., Haq, S. & Ghafoor, A. Meshless spectral method for solution of time-fractional coupled kdv equations. Appl. Math. Comput. 341, 321–334 (2019).
23. Esmaeelzade Aghdam, Y., Mesgarani, H. & Asadi, Z. Estimate of the fractional advection-diffusion equation with a time-fractional term based on the shifted legendre polynomials. J. Math. Model. 731–744 (2023).
24. Aghdam, Y. E., Mesgarani, H., Amin, A. & Gómez-Aguilar, J. An efficient numerical scheme to approach the time fractional black–scholes model using orthogonal gegenbauer polynomials. Comput. Econ. 1–14 (2023).
25. Aghdam YE Mesgarani H Moremedi G Khoshkhahtinat M High-accuracy numerical scheme for solving the space-time fractional advection-diffusion equation with convergence analysis Alex. Eng. J. 2022 61 217 225 10.1016/j.aej.2021.04.092
Aghdam, Y. E., Mesgarani, H., Moremedi, G. & Khoshkhahtinat, M. High-accuracy numerical scheme for solving the space-time fractional advection-diffusion equation with convergence analysis. Alex. Eng. J. 61, 217–225 (2022).10.1016/j.aej.2021.04.092
26. Mesgarani H Rashidnina J Esmaeelzade Aghdam Y Nikan O The impact of chebyshev collocation method on solutions of fractional advection-diffusion equation Int. J. Appl. Comput. Math. 2020 6 149 10.1007/s40819-020-00903-5
Mesgarani, H., Rashidnina, J., Esmaeelzade Aghdam, Y. & Nikan, O. The impact of chebyshev collocation method on solutions of fractional advection-diffusion equation. Int. J. Appl. Comput. Math. 6, 149 (2020).10.1007/s40819-020-00903-5
27. Owolabi, K. M., Agarwal, R. P., Pindza, E., Bernstein, S. & Osman, M. S. Complex turing patterns in chaotic dynamics of autocatalytic reactions with the caputo fractional derivative. Neural Comput. Appl. 1–27 (2023).
28. Alqhtani M Owolabi KM Saad KM Pindza E Spatiotemporal chaos in spatially extended fractional dynamical systems Commun. Nonlinear Sci. Numer. Simul. 2023 119 107118 10.1016/j.cnsns.2023.107118
Alqhtani, M., Owolabi, K. M., Saad, K. M. & Pindza, E. Spatiotemporal chaos in spatially extended fractional dynamical systems. Commun. Nonlinear Sci. Numer. Simul. 119, 107118 (2023).10.1016/j.cnsns.2023.107118
29. Alqhtani M Owolabi KM Saad KM Spatiotemporal (target) patterns in sub-diffusive predator-prey system with the caputo operator Chaos Solitons Fract. 2022 160 112267 10.1016/j.chaos.2022.112267
Alqhtani, M., Owolabi, K. M. & Saad, K. M. Spatiotemporal (target) patterns in sub-diffusive predator-prey system with the caputo operator. Chaos Solitons Fract. 160, 112267 (2022).10.1016/j.chaos.2022.112267
30. Owolabi KM Pindza E Atangana A Analysis and pattern formation scenarios in the superdiffusive system of predation described with caputo operator Chaos Solitons Fract. 2021 152 111468 10.1016/j.chaos.2021.111468
Owolabi, K. M., Pindza, E. & Atangana, A. Analysis and pattern formation scenarios in the superdiffusive system of predation described with caputo operator. Chaos Solitons Fract. 152, 111468 (2021).10.1016/j.chaos.2021.111468
31. Owolabi KM Karaagac B Baleanu D Dynamics of pattern formation process in fractional-order super-diffusive processes: A computational approach Soft. Comput. 2021 25 11191 11208 10.1007/s00500-021-05885-0
Owolabi, K. M., Karaagac, B. & Baleanu, D. Dynamics of pattern formation process in fractional-order super-diffusive processes: A computational approach. Soft. Comput. 25, 11191–11208 (2021).10.1007/s00500-021-05885-0
32. Owolabi KM Baleanu D Emergent patterns in diffusive turing-like systems with fractional-order operator Neural Comput. Appl. 2021 33 12703 12720 10.1007/s00521-021-05917-8
Owolabi, K. M. & Baleanu, D. Emergent patterns in diffusive turing-like systems with fractional-order operator. Neural Comput. Appl. 33, 12703–12720 (2021).10.1007/s00521-021-05917-8
33. Owolabi KM Jain S Spatial patterns through diffusion-driven instability in modified predator-prey models with chaotic behaviors Chaos Solitons Fract. 2023 174 113839 10.1016/j.chaos.2023.113839
Owolabi, K. M. & Jain, S. Spatial patterns through diffusion-driven instability in modified predator-prey models with chaotic behaviors. Chaos Solitons Fract. 174, 113839 (2023).10.1016/j.chaos.2023.113839
34. Owolabi KM Patidar KC Higher-order time-stepping methods for time-dependent reaction–diffusion equations arising in biology Appl. Math. Comput. 2014 240 30 50
Owolabi, K. M. & Patidar, K. C. Higher-order time-stepping methods for time-dependent reaction–diffusion equations arising in biology. Appl. Math. Comput. 240, 30–50 (2014).
35. Owolabi KM Mathematical analysis and numerical simulation of patterns in fractional and classical reaction-diffusion systems Chaos Solitons Fract. 2016 93 89 98 10.1016/j.chaos.2016.10.005
Owolabi, K. M. Mathematical analysis and numerical simulation of patterns in fractional and classical reaction-diffusion systems. Chaos Solitons Fract. 93, 89–98 (2016).10.1016/j.chaos.2016.10.005
36. Pindza E Owolabi KM Fourier spectral method for higher order space fractional reaction-diffusion equations Commun. Nonlinear Sci. Numer. Simul. 2016 40 112 128 10.1016/j.cnsns.2016.04.020
Pindza, E. & Owolabi, K. M. Fourier spectral method for higher order space fractional reaction-diffusion equations. Commun. Nonlinear Sci. Numer. Simul. 40, 112–128 (2016).10.1016/j.cnsns.2016.04.020
37. Das, S. Initialized differintegrals and generalized calculus. In Functional Fractional Calculus, 271–322 (Springer, 2011).
38. Sahin, A. Numerical solutions of the reaction–diffusion equations with B-spline finite element method Ph. D. Ph.D. thesis, dissertation. Department of Mathematics. Eskişehir Osmangazi University ... (2009).
39. Ersoy, O. & Dag, I. Numerical solutions of the reaction diffusion system by using exponential cubic b-spline collocation algorithms. Open Phys. 13 (2015).
40. Onarcan AT Adar N Dag I Trigonometric cubic b-spline collocation algorithm for numerical solutions of reaction-diffusion equation systems Comput. Appl. Math. 2018 37 6848 6869 10.1007/s40314-018-0713-4
Onarcan, A. T., Adar, N. & Dag, I. Trigonometric cubic b-spline collocation algorithm for numerical solutions of reaction-diffusion equation systems. Comput. Appl. Math. 37, 6848–6869 (2018).10.1007/s40314-018-0713-4
41. Hepson OE Numerical simulations of kuramoto-sivashinsky equation in reaction-diffusion via galerkin method Math. Sci. 2021 15 199 206 10.1007/s40096-021-00402-8
Hepson, O. E. Numerical simulations of kuramoto-sivashinsky equation in reaction-diffusion via galerkin method. Math. Sci. 15, 199–206 (2021).10.1007/s40096-021-00402-8
42. Haq S Uddin M A meshfree interpolation method for the numerical solution of the coupled nonlinear partial differential equations Eng. Anal. Boundary Elem. 2009 33 399 409 10.1016/j.enganabound.2008.06.005
Haq, S. et al. A meshfree interpolation method for the numerical solution of the coupled nonlinear partial differential equations. Eng. Anal. Boundary Elem. 33, 399–409 (2009).10.1016/j.enganabound.2008.06.005
43. Murio DA Implicit finite difference approximation for time fractional diffusion equations Comput. Math. Appl. 2008 56 1138 1145 10.1016/j.camwa.2008.02.015
Murio, D. A. Implicit finite difference approximation for time fractional diffusion equations. Comput. Math. Appl. 56, 1138–1145 (2008).10.1016/j.camwa.2008.02.015
44. Mittal R Jiwari R Numerical study of two-dimensional reaction-diffusion brusselator system by differential quadrature method Int. J. Comput. Methods Eng. Sci. Mech. 2011 12 14 25 10.1080/15502287.2010.540300
Mittal, R. & Jiwari, R. Numerical study of two-dimensional reaction-diffusion brusselator system by differential quadrature method. Int. J. Comput. Methods Eng. Sci. Mech. 12, 14–25 (2011).10.1080/15502287.2010.540300
45. Haq S Ali I Nisar KS A computational study of two-dimensional reaction-diffusion brusselator system with applications in chemical processes Alex. Eng. J. 2021 60 4381 4392 10.1016/j.aej.2021.02.064
Haq, S., Ali, I. & Nisar, K. S. A computational study of two-dimensional reaction-diffusion brusselator system with applications in chemical processes. Alex. Eng. J. 60, 4381–4392 (2021).10.1016/j.aej.2021.02.064
46. Ali A Haq S A computational modeling of the behavior of the two-dimensional reaction-diffusion brusselator system Appl. Math. Model. 2010 34 3896 3909 10.1016/j.apm.2010.03.028
Ali, A. et al. A computational modeling of the behavior of the two-dimensional reaction-diffusion brusselator system. Appl. Math. Model. 34, 3896–3909 (2010).10.1016/j.apm.2010.03.028
