A Comparative Study of the Fractional-Order Belousov–Zhabotinsky System

Abstract

In this article, we present a modified strategy that combines the residual power series method with the Laplace transformation and a novel iterative technique for generating a series solution to the fractional nonlinear Belousov–Zhabotinsky (BZ) system. The proposed techniques use the Laurent series in their development. The new procedures’ advantages include the accuracy and speed in obtaining exact/approximate solutions. The suggested approach examines the fractional nonlinear BZ system that describes flow motion in a pipe.

1. Introduction

Fractional calculus (FC) is considered a branch of mathematical analysis that generalizes traditional calculus by allowing non integer order integral and derivative. The concept of FC dates back to the 17th century when Leibniz and L’Hopital studied the possibility of defining a derivative of non-integer order. However, the modern development of FC began in the 19th century with the work of Liouville, Riemann, and Grunwald, who independently introduced fractional derivatives and integrals [1,2,3,4]. FC has many application in various field of engineering and science, including chemistry, physics, economics, biology, and finance. For instance, it has been applied to models complex scheme such as viscoelastic materials, electrical circuit, and signal processing. FC also plays a critical role in the analysis of anomalous diffusion and stochastic processes, where the conventional tools of calculus are not sufficient. FC has attract increases attentions in recent decades, and researchers have explored new applications and properties of fractional derivatives and integrals. The development of numerical methods for FC has also opened up new possibilities for solving fractional differential equations (FDEs), making it a valuable tools in several scientific fields and engineering problems [5,6,7,8,9].

Partial differential equations (PDEs) are mathematical equations that described physical processes involving multiple independent variables, such as time and space. They play a crucial role in modeling many phenomena in science and engineering, including fluid dynamics, electromagnetism, quantum mechanics, and plasma physics [10,11,12,13,14,15,16,17,18,19]. The important area of partial differential equations is the fractional system of PDEs, which extends the traditional integer-order calculus to non-integer orders [20,21,22]. The fractional system of PDEs has achieved important attentions in recent decades due to its effect to model complex phenomena that cannot be described by the classical integer-order models. The fractional system of PDEs has a broad ranges of implementation in many field, including geophysics, biology, finance, physics, and engineering. It provides a powerful tool for modeling and analyzing complex systems that exhibit long-range dependence, memory, and fractal behavior [23,24,25,26,27]. In this essay, we provide an overview of the fractional system of PDEs, including its definition, properties, and applications. We also discuss the numerical methods used to solve these equations and some open problems in this field [28,29,30].

Many application of FDEs in applied sciences such as electro-dynamics, accounting, chaos ideas, biological populations design and fluid mechanic digital signals, FDEs are more many areas of sciences [31,32,33,34,35]. In this study, we aim to employ an efficient analytical approach to address nonlinear differential equations of arbitrary order (ODEs). By doing so, it is possible to improve the accuracy of analysis in related fields through the use of FDEs. Numerous approaches have been devised to address this issue; One of the methods employed is the Adomian decomposition method (DM) [36], the reduced differential transform method [37], variational iteration method [38], the Elzaki DM [39,40], the iterative transformation technique [41], the Natural DT [42], the homotopy perturbation technique (HPT) [43], and so on [44,45,46,47,48].

The BZ reaction is a group of chemical responses that exhibits oscillatory behavior. These reactions involve the catalytic oxidative stress of numerous reductants, typically natural compounds, by bromic acetone in an acidic aqueous solvent, facilitated by transition-metal ions. Most BZ reactions occur in a homogeneous phase. One significant advantage of the BZ reaction is that it enables the observation of the formation of intricate patterns over space and time, visible to the naked eye, on a convenient sentient timeframe of tens of secs and a spatial extent of a few millimeters. The BZ reaction can produce several thousand oscillatory cycles in a closed system, allowing the study of chemical waves and patterns without requiring a constant supply of reactants [49]. The literature reports various mathematical methods develop to obtain numeric result over a specifics ranges or to approximate solutions using a limited number of terms in an iterative computational series. These available methods include the Laplace iterative method [50], the variational iteration technique [51], homotopy analysis perturbation methods [52], Adomian’s decomposition method [53], and the residual power series method (RPSM) [54].

The RPSM is a commonly applied technique to solve integral and differential equations of both integer and fractional orders. The RPSM is a method introduced by Omar et al. mathematicians, in 2013. It is designed to be a fast and simple way to calculate the coefficient of power series solution for differential fuzzy equations. The approach involves suppose that the result to the equations can be expressed as a power series and then finding the coefficients of this series [55]. Unlike other methods that require perturbation, linearization, or discretization, the RPSM provides a straightforward solution for highly linear and nonlinear equations without these requirements. It has been used to analysis a varieties of non-linear ODEs and PDEs with different orders and classes. For instance, the RPSM was used to solve the generalized Lane–Emden equation, to approximate solutions to the nonlinear fractional Korteweg-De Vries Equation-Burger equation, and to predicts the solitary pattern results of nonlinear fractional dispersive PDEs (Abu Arqub, 2013; Al-Khaled & Abu Arqub, 2017; El-Kalla et al., 2021) [56,57,58,59]. The RPSM has several advantages over other analytical and numerical approaches. Firstly, it does not require a recursive connection or comparison of coefficients of related terms. Secondly, it provides a simplified technique to ensures the convergences of the series result by reducing the associated residual errors. Thirdly, the residual power series method does not suffers from mathematical rounds error and does not consume significant time or memory. Fourthly, it can be immediately used to the proposed model by selecting an appropriates starting condition approximations, without requiring any conversion when transitioning from lower to higher orders (Al-Khaled & Abu Arqub, 2017; El-Kalla et al., 2021) [60,61,62]. In this paper, we utilize the Laplace RPSM (LRPSM) to obtain a precise solutions for nonlinear fractional PDEs. By integrating the RPSM with the LT, we present a renewabile algorithmic method that generates insightful results through a convergent series. The fractional Caputo derivative enables us to categorize the PDEs quantitatively [63,64,65]. The exact analytical results obtained through this methodology provide a valuable tool for analyzing complex system dynamics, especially for computational fractional PDEs (FPDEs) [66,67,68].

2. Preliminaries

Definition 1.

The Caputo fractional derivative of a feature $\mu (\varsigma ,\tau )$ of order α can be expressed as follows [69]

$${}^C{D}_{\tau }^{\alpha }\mu (\varsigma \varsigma ,\tau )=J_{\tau }^{n-\alpha }{\mu }^n(\varsigma ,\tau ),n-1<\alpha \le n,t>0,$$

(1)

where m is a natural number, and $J_{\tau }^{\alpha }$ represents the Riemann–Liouville fractional integral of $\mu (\varsigma ,t)$ of order α expressed as

$$J_{\tau }^{\sigma }\mu (\varsigma ,\tau )=\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^{\tau }{(\tau -\rho )}^{\alpha -1}\mu (\varsigma ,\rho )d\rho .$$

(2)

Definition 2.

The Laplace transform (LT) of $\mu (\varsigma ,\tau )$ is defined by [69]

$$\mu (\varsigma ,s)={\mathcal{L}}_{\tau }\left[\mu (\varsigma ,\tau )\right]={\int }_0^{\infty }{e}^{-s\tau }\mu (\varsigma ,\tau )d\tau ,s>\alpha ,$$

(3)

whereas the inverse of the LT reads

$$\mu (\varsigma ,\tau )={\mathcal{L}}_{\tau }^{-1}\left[\mu (\varsigma ,s)\right]={\int }_{l-i\infty }^{l+i\infty }{e}^{s\tau }\mu (\varsigma ,s)ds,l=Re\left(s\right)>l_0.$$

(4)

Lemma 1.

Assume that $u(\varsigma ,\tau )$ is a piecewises continuous terms, and $U(\varsigma ,s)=\mathcal{L}\left[u\right(\varsigma ,\tau \left)\right]$, we obtain

1.

$\mathcal{L}\left[J_{\tau }^{\alpha }u(\varsigma ,\tau )\right]=\frac{U(\varsigma ,s)}{s^{\alpha }},\alpha >0.$

2.

$\mathcal{L}\left[D_{\tau }^{\alpha }u(\varsigma ,\tau )\right]=s^{\sigma }U(\varsigma ,s)-{\sum }_{k=0}^{m-1}{s}^{\alpha -k-1}{u}^k(\varsigma ,0),n-1<\alpha \le n.$

3.

$\mathcal{L}\left[D_{\tau }^{n\alpha }u(\varsigma ,\tau )\right]=s^{n\alpha }U(\varsigma ,s)-{\sum }_{k=0}^{n-1}{s}^{(n-k)\alpha -1}{D}_{\tau }^{k\alpha }u(\varsigma ,0),0<\alpha \le 1.$

Proof. 

For the proof, see Ref. [69]. □

Theorem 1.

Consider a piecewise continuous function $u(\varsigma ,\tau )$ defined on the interval $I\times [0,\infty )$ and possessing exponential order ζ. The Laplace transform of $u(\varsigma ,\tau )$, $U(\varsigma ,s)$, has a fractional representation as follows:

$$U(\varsigma ,s)=\sum _{n=0}^{\infty }\frac{f_n\left(\varsigma \right)}{s^{1+n\alpha }},0<\alpha \le 1,\varsigma \in I,s>\varsigma .$$

(5)

Then, $f_n\left(\zeta \right)=D_t^{n\sigma }u(\zeta ,0)$.

Proof. 

For the proof, see Ref. [69]. □

Remark 1.

The inverse of the LT of Equation (5) reads [69]:

$$u(\varsigma ,\tau )=\sum _{i=0}^{\infty }\frac{D_{\tau }^{\alpha }u(\varsigma ,0)}{\mathsf{\Gamma }(1+i\alpha )}{\tau }^{i\left(\zeta \right)},0<\varsigma \le 1,t\ge 0.$$

(6)

3. General Implementation of the Suggested Methods

3.1. Laplace Residual Power Series Method (LRPSM)

Consider the following FPDE

$$\begin{array}{cc}\hfill & D_{\tau }^{\rho }\mu (\varsigma ,\tau )+N\left[\mu (\varsigma ,\tau )\right]+R\left[\mu (\varsigma ,\tau )\right]=0,\mathrm{where}0<\rho \le 1,\hfill \end{array}$$

(7)

which is subjected to the initial condition (IC):

$$\begin{array}{c}\hfill \mu (\varsigma ,\tau )=f_0\left(\varsigma \right).\end{array}$$

(8)

Applying the Laplace transform to Equation (7) and use Equation (8), we get

$$\begin{array}{c}\hfill \mu (\varsigma ,s)-\frac{f_0(\varsigma ,s)}{s}+\frac{1}{s^{\rho }}{\mathcal{L}}_{\tau }\left[N\left[{\mathcal{L}}_{\tau }^{-1}\left[\mu (\varsigma ,s)\right]\right]+A\left[\mu (\varsigma ,\tau )\right]\right]=0.\end{array}$$

(9)

It is assumed that the solution to Equation (9) can be expressed using the following expansion

$$\begin{array}{cc}\hfill & \mu (\varsigma ,s)=\sum _{n=0}^{\infty }\frac{f_n(\varsigma ,s)}{s^{n\rho +1}},\hfill \end{array}$$

(10)

and the kth-truncated term series reads

$$\begin{array}{c}\hfill \mu (\varsigma ,s)=\frac{f_0(\varsigma ,s)}{s}+\sum _{n=1}^k\frac{f_n(\varsigma ,s)}{s^{n\rho +1}},n=1,2,3,4\cdots \end{array}$$

(11)

$$\begin{array}{c}\hfill {\mathcal{L}}_{\tau }Res(\varsigma ,s)=\mu (\varsigma ,s)-\frac{f_0(\varsigma ,s)}{s}+\frac{1}{s^{\rho }}{\mathcal{L}}_{\tau }\left[N\left[{\mathcal{L}}_{\tau }^{-1}\left[\mu (\varsigma ,s)\right]\right]+A\left[\mu (\varsigma ,\tau )\right]\right].\end{array}$$

(12)

In addition, the kth LRF is:

$$\begin{array}{c}\hfill {\mathcal{L}}_{\tau }Re{s}_k(\varsigma ,s)={\mu }_k(\varsigma ,s)-\frac{f_0(\varsigma ,s)}{s}+\frac{1}{s^{\rho }}{\mathcal{L}}_{\tau }\left[N\left[{\mathcal{L}}_{\tau }^{-1}\left[{\mu }_k(\varsigma ,s)\right]\right]+A\left[{\mu }_k(\varsigma ,\tau )\right]\right].\end{array}$$

(13)

The above coefficients can be calculated by recursively solving the following system using $f_n(\varsigma ,s)$.

$$\begin{array}{cc}\hfill & \underset{s\to \infty }{lim}{s}^{k\alpha +1}{\mathcal{L}}_{\tau }Re{s}_{\mu ,k}(\alpha ,s)=0,k=1,2,\cdots .\hfill \end{array}$$

(14)

Finally, the inverse of the LT to Equation (10) is considered to obtain the kth analytical result of ${\mu }_k(\varsigma ,\tau )$.

Theorem 2.

Consider the following FPDE in $D\subset {\mathbb{R}}^n$ with $n\ge 1$ and $s\in (0,1]$:

$${\mathcal{D}}_{\tau }^{s}u(\varsigma ,\tau )=f(\varsigma ,\tau ),\varsigma \in D,\tau >0,$$

(15)

where $\mathcal{D}{\tau }^s$ denotes the fractional Caputo operator of order s with regard to τ, and f is a given component. Suppose that $u(\varsigma ,\tau )$ is sufficiently smooth and satisfies suitable initial and/or boundary conditions.

Let $u_n(\varsigma ,\tau )$ be the Laplace Residual power series approximation to $u(\varsigma ,\tau )$, which can be obtained by solving the following iterated problem

$$\begin{array}{ccc}\hfill s{\int }_0^{\infty }{e}^{-s\tau }{\tau }^{s-1}\Delta u_{n+1}(\varsigma ,\tau )d\tau & =\Delta u_n(\varsigma ,0)-f(\varsigma ,0),s{\int }_0^{\infty }{e}^{-s\tau }{\tau }^{s-1}\Delta u_{n+1}(\varsigma ,\tau )d\tau \hfill & \hfill ={\mathcal{D}}_{\tau }^s\Delta u_n(\varsigma ,\tau )-{\mathcal{D}}_{\tau }^{s}f(\varsigma ,\tau ),\tau >0,\end{array}$$

where Δ denotes the Laplacian operator with respect to ς.

Then, under suitable conditions on the initial/boundary data and the function f, the sequence $u_n$ converges to the unique solution u of the FPDE (15) in a suitable norm, as $n\to \infty$.

3.2. General Application of NIM

To discuss the fundamental concept of the new iterative approach, we examine the functional equation in a broad sense:

$$\begin{array}{c}\hfill \mu \left(\varsigma \right)=f\left(\varsigma \right)+N\left(\mu \right(\varsigma \left)\right),\end{array}$$

(16)

Let N be a nonlinear operator that maps from a Banach space B to itself, and let f be an unknown function.

$$\begin{array}{cc}\hfill & \mu \left(\varsigma \right)=\sum _{i=0}^{\infty }{\mu }_i\left(\varsigma \right).\hfill \end{array}$$

(17)

The non-linear terms can be expressed as

$$\begin{array}{cc}\hfill & N(\sum _{i=0}^{\infty }{\mu }_i\left(\varsigma \right))=N\left({\mu }_0\right)+\sum _{i=0}^{\infty }\left[N\left(\sum _{j=0}^i{\mu }_j\left(\varsigma \right)\right)-N\left(\sum _{j=0}^{i-1}{\mu }_j\left(\varsigma \right)\right)\right].\hfill \end{array}$$

(18)

Inserting Equations (17) and (18) into (16), we obtain

$$\begin{array}{c}\hfill \sum _{i=0}^{\infty }{\mu }_i\left(\varsigma \right)=f+N\left({\mu }_0\right)+\sum _{i=0}^{\infty }\left[N(\sum _{j=0}^i{\mu }_j\left(\varsigma \right))-N(\sum _{j=0}^{i-1}{\mu }_j\left(\varsigma \right))\right].\end{array}$$

(19)

The following recurrence relations are introduced

$$\begin{array}{cc}\hfill & {\mu }_0=f,\hfill \\ \hfill & {\mu }_1=N\left({\mu }_0\right),\hfill \\ \hfill & {\mu }_2=N({\mu }_0+{\mu }_1)-N\left({\mu }_0\right),\hfill \\ \hfill & {\mu }_{n+1}=N({\mu }_0+{\mu }_1+\cdots {\mu }_n)-N({\mu }_0+{\mu }_1+\cdots {\mu }_{n-1}),n=1,2,3\cdots .\hfill \end{array}$$

(20)

Then,

$$\begin{array}{cc}\hfill & ({\mu }_0+{\mu }_1+\cdots {\mu }_n)=N({\mu }_0+{\mu }_1+\cdots {\mu }_n),n=1,2,3\cdots ,\hfill \\ \hfill & \mu =\sum _{i=0}^{\infty }{\mu }_i\left(\varsigma \right)=f+N\left(\sum _{i=0}^{\infty }{\mu }_i\left(\varsigma \right)\right).\hfill \end{array}$$

(21)

4. Appropriate Algorithmic Approach

In this section, we represent a viable technique for investigating nonlinear fractional PDEs, utilizing a novel iterative approach.

$$\begin{array}{c}\hfill D_{\tau }^{\alpha }\mu (\varsigma ,\tau )=A(\mu ,\partial \mu )+B(\varsigma ,\tau ),m-1<\alpha \le m,mϵN,\end{array}$$

(22)

with the IC

$$\begin{array}{c}\hfill \frac{{\partial }^k}{\partial {\tau }^k}\mu (\varsigma ,0)=h_k\left(\varsigma \right),k=0,1,2,3\cdots m-1,\end{array}$$

(23)

The nonlinear function A is dependent on $\mu$ and the partial derivative of $\mu$ with respect to both $\varsigma$ and t. B represents the source function. With the implementation of the new iterative method, the initial value problem described in Equations (22) and (23) can be expressed as an equivalent integral equation

$$\begin{array}{c}\hfill \mu (\varsigma ,\tau )=\sum _{k=0}^{m-1}{h}_k\left(\varsigma \right)\frac{t^k}{k!}+I_{\tau }^{\mu }\left(A\right)+I_{\tau }^{\mu }\left(B\right)=f+N\left(\mu \right),\end{array}$$

(24)

with

$$\begin{array}{c}\hfill f=\sum _{k=0}^{m-1}{h}_k\left(\varsigma \right)\frac{t^k}{k!}+I_t^{\alpha }\left(B\right),\end{array}$$

(25)

$$\begin{array}{c}\hfill N\left(\omega \right)=I_{\tau }^{\alpha }\left(A\right).\end{array}$$

(26)

Theorem 3

(Convergence of New Iterative Method). Let $u^{\left(k\right)}$ be the sequence generated by the new iterative method for solving the following FPDE

$$\begin{array}{cc}\hfill {\mathcal{D}}_{\tau }^{\alpha }u(\varsigma ,\tau )& =\mathcal{L}u(\varsigma ,\tau )+f(\varsigma ,\tau ),0<\tau \le T,\varsigma \in \Omega \hfill \\ \hfill u(\varsigma ,0)& =u_0\left(\varsigma \right),\varsigma \in \Omega ,\hfill \end{array}$$

(27)

where ${\mathcal{D}}_{\tau }^{\alpha }$ is the Caputo fractional derivative of order $\alpha \in (0,1]$, $\mathcal{L}$ is a linear differential operator, $f(\varsigma ,\tau )$ is a given function, $u_0\left(\varsigma \right)$ is the initial condition, Ω is a bounded domain in ${\mathbb{R}}^n$ with smooth boundary $\partial \mathsf{\Omega }$, and $T>0$ is the final time.

The following conditions should be hold:

1.

The operator $\mathcal{L}$ is uniformly elliptic and has smooth coefficients.

2.

The initial condition $u_0\left(\varsigma \right)$ is bounded and measurable.

3.

The function $f(\varsigma ,\tau )$ is continuous and locally bounded in $\mathsf{\Omega }\times (0,T]$.

4.

The solution $u(\varsigma ,\tau )$ of (27) satisfies the following estimate:

$${\left|u(\varsigma ,\tau )\right|}_{L^{\infty }(\Omega )}\le M\left(\tau \right),0<\tau \le T,$$

(28)

where $M\left(\tau \right)$ is a nonnegative and increasing function.

Then, the sequence $u^{\left(k\right)}$ converges uniformly to the unique solution $u(\varsigma ,\tau )$ of (27) as $k\to \infty$. Moreover, the convergence rate is at least of order $O\left({\lambda }^{2k}\right)$, where $\lambda \in (0,1)$ is the relaxation parameter used in the iterative method.

5. Applications

5.1. Implementation of the LRPSM

The following fractional BZ model is considered,

$$\begin{array}{cc}\hfill & D_{\tau }^{\alpha }\mu (\varsigma ,\tau )=\mu (\varsigma ,\tau )(1-\mu (\varsigma ,\tau )-r\nu (\varsigma ,\tau ))+{\mu }_{\varsigma \varsigma }(\varsigma ,\tau ),\hfill \\ \hfill & D_{\tau }^{\alpha }\nu (\varsigma ,\tau )=-a\mu (\varsigma ,\tau )\nu (\varsigma ,\tau )+{\nu }_{\varsigma \varsigma }(\varsigma ,\tau ),0<\alpha \le 1.\hfill \end{array}$$

(29)

Example 1.

Consider the fractional BZ with $r=2$ and $a=3$:

$$\begin{array}{cc}\hfill & D_{\tau }^{\alpha }\mu (\varsigma ,\tau )-\mu (\varsigma ,\tau )-\frac{{\partial }^2\mu (\varsigma ,\tau )}{\partial {\varsigma }^2}+{\mu }^2(\varsigma ,\tau )+2\mu (\varsigma ,\tau )\nu (\varsigma ,\tau )=0,\hfill \\ \hfill & D_{\tau }^{\alpha }\nu (\varsigma ,\tau )-\frac{{\partial }^2\nu (\varsigma ,\tau )}{\partial {\varsigma }^2}+3\mu (\varsigma ,\tau )\nu (\varsigma ,\tau )=0,\mathit{where}0<\alpha \le 1,\hfill \end{array}$$

(30)

subject to the following ICs:

$$\begin{array}{cc}\hfill & \mu (\varsigma ,0)=-\frac{1}{2}\left(1-tan{h}^2(\frac{\varsigma }{2})\right),\hfill \\ \hfill & \nu (\varsigma ,0)=-\frac{1}{2}+tanh(\frac{\varsigma }{2})+\frac{1}{2}tan{h}^2(\frac{\varsigma }{2}).\hfill \end{array}$$

(31)

Using the LT of Equation (30), we obtain

$$\begin{array}{cc}\hfill & \mu (\varsigma ,s)-\frac{\frac{1}{2}\left(1-tan{h}^2\left(\frac{\varsigma }{2}\right)\right)}{s}-\frac{1}{s^{\alpha }}[\mu (\varsigma ,s)+\frac{{\partial }^2\mu (\varsigma ,s)}{\partial {\varsigma }^2}-2\mathcal{L}\left({\mathcal{L}}_{\tau }^{-1}\left[\mu (\varsigma ,\tau )\right]{\mathcal{L}}_{\tau }^{-1}\left[\nu (\varsigma ,\tau )\right]\right)\hfill \\ & -\mathcal{L}{\left({\mathcal{L}}_{\tau }^{-1}\left[\mu (\varsigma ,s)\right]\right)}^2]=0,\hfill \\ \\ \hfill & \nu (\varsigma ,s)+\frac{-\frac{1}{2}+tanh\left(\frac{\varsigma }{2}\right)+\frac{1}{2}tan{h}^2\left(\frac{\varsigma }{2}\right)}{s}-\frac{1}{s^{\alpha }}\left[\frac{{\partial }^2\mu (\varsigma ,s)}{\partial {\varsigma }^2}-3\mathcal{L}\left({\mathcal{L}}_{\tau }^{-1}\left[\mu (\varsigma ,s)\right]{\mathcal{L}}_{\tau }^{-1}\left[\nu (\varsigma ,\tau )\right]\right)\right]=0,\hfill \end{array}$$

(32)

and so the $kth$-truncated term series are

$$\begin{array}{cc}\hfill & \mu (\varsigma ,s)=-\frac{\frac{1}{2}\left(1-tan{h}^2\left(\frac{\varsigma }{2}\right)\right)}{s}+\sum _{n=1}^k\frac{f_n(\varsigma ,s)}{s^{n\alpha +1}},\hfill \\ \hfill & \nu (\varsigma ,s)=\frac{-\frac{1}{2}+tanh\left(\frac{\varsigma }{2}\right)+\frac{1}{2}tan{h}^2\left(\frac{\varsigma }{2}\right)}{s}+\sum _{n=1}^k\frac{g_n(\varsigma ,s)}{s^{n\alpha +1}},\hfill \\ & k=1,2,3,4\cdots .\hfill \end{array}$$

(33)

The residual Laplace function is given by

$$\begin{array}{cc}\hfill & {\mathcal{L}}_{t}Re{s}_u(\zeta ,s)=\mu (\varsigma ,s)-\frac{\frac{1}{2}\left(1-tan{h}^2\left(\frac{\varsigma }{2}\right)\right)}{s}-\frac{1}{s^{\alpha }}[\mu (\varsigma ,s)+\frac{{\partial }^2\mu (\varsigma ,s)}{\partial {\varsigma }^2}-2\mathcal{L}\left({\mathcal{L}}_{\tau }^{-1}\left[\mu (\varsigma ,\tau )\right]{\mathcal{L}}_{\tau }^{-1}\left[\nu (\varsigma ,\tau )\right]\right)\hfill \\ & -\mathcal{L}{\left({\mathcal{L}}_{\tau }^{-1}\left[\mu (\varsigma ,s)\right]\right)}^2],\hfill \\ \\ \hfill & \mathcal{L}Re{s}_u(\zeta ,s)=\nu (\varsigma ,s)+\frac{-\frac{1}{2}+tanh\left(\frac{\varsigma }{2}\right)+\frac{1}{2}tan{h}^2\left(\frac{\varsigma }{2}\right)}{s}-\frac{1}{s^{\alpha }}\left[\frac{{\partial }^2\mu (\varsigma ,s)}{\partial {\varsigma }^2}-3\mathcal{L}\left({\mathcal{L}}_{\tau }^{-1}\left[\mu (\varsigma ,s)\right]{\mathcal{L}}_{\tau }^{-1}\left[\nu (\varsigma ,\tau )\right]\right)\right],\hfill \end{array}$$

(34)

and the kth-LRFs are:

$$\begin{array}{cc}\hfill & \mathcal{L}Re{s}_u(\zeta ,s)={\mu }_k(\varsigma ,s)-\frac{\frac{1}{2}\left(1-tan{h}^2\left(\frac{\varsigma }{2}\right)\right)}{s}-\frac{1}{s^{\alpha }}[{\mu }_k(\varsigma ,s)+\frac{{\partial }^2{\mu }_k(\varsigma ,s)}{\partial {\varsigma }^2}-2\mathcal{L}\left({\mathcal{L}}_{\tau }^{-1}\left[{\mu }_k(\varsigma ,\tau )\right]{\mathcal{L}}_{\tau }^{-1}\left[{\nu }_k(\varsigma ,\tau )\right]\right)\hfill \\ & -\mathcal{L}{\left({\mathcal{L}}_{\tau }^{-1}\left[{\mu }_k(\varsigma ,s)\right]\right)}^2],\hfill \\ \\ \hfill & \mathcal{L}Re{s}_u(\zeta ,s)={\nu }_k(\varsigma ,s)+\frac{-\frac{1}{2}+tanh\left(\frac{\varsigma }{2}\right)+\frac{1}{2}tan{h}^2\left(\frac{\varsigma }{2}\right)}{s}-\frac{1}{s^{\alpha }}\left[\frac{{\partial }^2{\mu }_k(\varsigma ,s)}{\partial {\varsigma }^2}-3\mathcal{L}\left({\mathcal{L}}_{\tau }^{-1}\left[{\mu }_k(\varsigma ,s)\right]{\mathcal{L}}_{\tau }^{-1}\left[{\nu }_k(\varsigma ,\tau )\right]\right)\right].\hfill \end{array}$$

(35)

We calculate $fk(\varsigma ,s)$, with k ranging from 1 to infinity, by putting the $kth$ truncated series from Equation (33) into the $kth$ Laplace residual term in Equation (35), multiplying the solution equations by $s^{k\alpha +1}$, and then solving respectively the limit $lims\to \infty \left(s^{k\alpha +1}\mathcal{L}Resu,k(\zeta ,s)\right)=0$, for each $k=1,2,3,\cdots$. The first few terms in this calculation are as follows:

$$\begin{array}{cc}\hfill & f_0(\varsigma ,s)=-\frac{1}{2}\left(1-tan{h}^2\left(\frac{\varsigma }{2}\right)\right),\hfill \\ \hfill & g_0(\varsigma ,s)=-\frac{1}{2}+tanh\left(\frac{\varsigma }{2}\right)+\frac{1}{2}tan{h}^2\left(\frac{\varsigma }{2}\right),\hfill \\ \hfill & f_1(\varsigma ,s)=csc{h}^3\left(\varsigma \right)sin{h}^4\left(\frac{\varsigma }{2}\right),\hfill \\ \hfill & g_1(\varsigma ,s)=\frac{-1+tanh\left(\frac{\varsigma }{2}\right)}{1+cosh\left(\varsigma \right)},\hfill \\ \hfill & f_2(\varsigma ,s)=\frac{8e^{\varsigma }\left(2+e^{\varsigma }(-5+e^{\varsigma })\right)}{{(1+e^{\varsigma })}^5},\hfill \\ \hfill & g_2(\varsigma ,s)=\frac{2e^{\varsigma }\left(3+e^{\varsigma }(13+e^{\varsigma }(-31+7e^{\varsigma })))\right)}{{(1+e^{\varsigma })}^5}.\hfill \end{array}$$

(36)

Put the values of $f_k(\zeta ,s)$, $k=1,2,3,\cdots$ into Equation (33), we obtain

$$\begin{array}{cc}\hfill & \mu (\varsigma ,s)=-\frac{\frac{1}{2}\left(1-tan{h}^2\left(\frac{\varsigma }{2}\right)\right)}{s}+\frac{csc{h}^3\left(\varsigma \right)sin{h}^4\left(\frac{\varsigma }{2}\right)}{s^{\alpha +1}}+\frac{8e^{\varsigma }\left(2+e^{\varsigma }(-5+e^{\varsigma })\right)}{{(1+e^{\varsigma })}^5{s}^{2\alpha +1}}+\cdots .\hfill \\ \hfill & \nu (\varsigma ,s)=\frac{-\frac{1}{2}+tanh\left(\frac{\varsigma }{2}\right)+\frac{1}{2}tan{h}^2\left(\frac{\varsigma }{2}\right)}{s}+\frac{-1+tanh\left(\frac{\varsigma }{2}\right)}{1+cosh\left(\varsigma \right)}\frac{1}{s^{\alpha +1}}+\hfill \\ & \frac{2e^{\varsigma }\left(3+e^{\varsigma }(13+e^{\varsigma }(-31+7e^{\varsigma })))\right)}{{(1+e^{\varsigma })}^5{s}^{2\alpha +1}}+\cdots .\hfill \end{array}$$

(37)

Using the inverse of the LT, we obtain

$$\begin{array}{cc}\hfill & \mu (\varsigma ,\tau )=-\frac{1}{2}\left(1-tan{h}^2\left(\frac{\varsigma }{2}\right)\right)+\frac{csc{h}^3\left(\varsigma \right)sin{h}^4\left(\frac{\varsigma }{2}\right)}{\mathsf{\Gamma }(\alpha +1)}{\tau }^{\alpha }+\frac{8e^{\varsigma }\left(2+e^{\varsigma }(-5+e^{\varsigma })\right)}{{(1+e^{\varsigma })}^5\mathsf{\Gamma }(2\alpha +1)}{\tau }^{2\alpha }+\cdots .\hfill \\ \hfill & \nu (\varsigma ,\tau )=-\frac{1}{2}+tanh\left(\frac{\varsigma }{2}\right)+\frac{1}{2}tan{h}^2\left(\frac{\varsigma }{2}\right)+\frac{-1+tanh\left(\frac{\varsigma }{2}\right)}{\left(1+cosh\left(\varsigma \right)\right)\mathsf{\Gamma }(\alpha +1)}{\tau }^{\alpha }+\hfill \\ & \frac{2e^{\varsigma }\left(3+e^{\varsigma }(13+e^{\varsigma }(-31+7e^{\varsigma })))\right)}{{(1+e^{\varsigma })}^5\mathsf{\Gamma }(2\alpha +1)}{\tau }^{2\alpha }+\cdots .\hfill \end{array}$$

(38)

5.2. Implementation of NIM

By applying the RL integral $I_{\tau }^{\alpha }$ to both sides of Equation (30) and using Equation (31), we obtain the equivalent integral form:

$$\begin{array}{cc}\hfill & \mu (\varsigma ,\tau )=-\frac{1}{2}\left(1-tan{h}^2\left(\frac{\varsigma }{2}\right)\right)+I_{\tau }^{\alpha }\left[\mu (\varsigma ,\tau )+\frac{{\partial }^2\mu (\varsigma ,\tau )}{\partial {\varsigma }^2}-{\mu }^2(\varsigma ,\tau )-2\mu (\varsigma ,\tau )\nu (\varsigma ,\tau )\right],\hfill \\ \hfill & \nu (\varsigma ,\tau )=-\frac{1}{2}+tanh\left(\frac{\varsigma }{2}\right)+\frac{1}{2}tan{h}^2\left(\frac{\varsigma }{2}\right)+I_{\tau }^{\alpha }\left[\frac{{\partial }^2\mu (\varsigma ,\tau )}{\partial {\varsigma }^2}-3\mu (\varsigma ,\tau )\nu (\varsigma ,\tau )\right].\hfill \end{array}$$

(39)

Using the NIM formulation that is discussed in Section 3, we obtain

$$\begin{array}{cc}\hfill & f_0(\varsigma ,s)=-\frac{1}{2}\left(1-tan{h}^2\left(\frac{\varsigma }{2}\right)\right),g_0(\varsigma ,s)=-\frac{1}{2}+tanh\left(\frac{\varsigma }{2}\right)+\frac{1}{2}tan{h}^2\left(\frac{\varsigma }{2}\right),\hfill \\ \hfill & f_1(\varsigma ,s)=\frac{8e^{\varsigma }{\tau }^{\alpha }}{(1+e^{\varsigma }{)}^3\mathsf{\Gamma }(\alpha +1)},g_1(\varsigma ,s)=\frac{2e^{\varsigma }(-3+e^{\varsigma }){\tau }^{\alpha }}{(1+e^{\varsigma }{)}^3\mathsf{\Gamma }(\alpha +1)},\hfill \\ \hfill & f_2(\varsigma ,s)=-\frac{8e^{\varsigma }{\tau }^{2\alpha }}{{(1+e^{\varsigma })}^6}\left(\frac{{(1+e^{\varsigma })}^2(4+3e^{\varsigma })}{\mathsf{\Gamma }(2\alpha +1)}+\frac{24{\left(2\right)}^{2\alpha }{e}^{\varsigma }{\tau }^{\alpha }\mathsf{\Gamma }(\frac{1}{2}+\alpha )}{\sqrt{\pi }\mathsf{\Gamma }(1+\alpha )\mathsf{\Gamma }(1+3\alpha )}\right),\hfill \\ \hfill & g_2(\varsigma ,s)=\frac{2e^{\varsigma }{\tau }^{2\alpha }}{{(1+e^{\varsigma })}^6}(-\frac{96e^{\varsigma }{\tau }^{\alpha }}{{\mathsf{\Gamma }}^3(1+\alpha )}+\frac{(1+e^{\varsigma }(-3+e^{\varsigma }(25+e^{\varsigma }(-19+e^{\varsigma }))))}{\mathsf{\Gamma }(1+2\alpha )}\hfill \\ \hfill +& \frac{24e^{\varsigma }(1+e^{\varsigma })(-2+sinh\left(\varsigma \right))}{{\mathsf{\Gamma }}^2(1+\alpha )}).\hfill \end{array}$$

(40)

The expressions for the result of $\mu (\varsigma ,\tau )$ and $\nu (\varsigma ,\tau )$ read

$$\begin{array}{cc}\hfill & \mu \left(\xi \right)=\sum _{i=0}^{\infty }{\mu }_i(\varsigma ,\tau ),\hfill \\ \hfill & \nu \left(\xi \right)=\sum _{i=0}^{\infty }{\nu }_i(\varsigma ,\tau ),\hfill \end{array}$$

(41)

which for some limited terms, we obtain

$$\begin{array}{cc}\hfill & \mu (\varsigma ,\tau )=-\frac{1}{2}\left(1-tan{h}^2\left(\frac{\varsigma }{2}\right)\right)+\frac{8e^{\varsigma }{\tau }^{\alpha }}{(1+e^{\varsigma }{)}^3\mathsf{\Gamma }(\alpha +1)}\hfill \\ & -\frac{8e^{\varsigma }{\tau }^{2\alpha }}{{(1+e^{\varsigma })}^6}\left(\frac{{(1+e^{\varsigma })}^2(4+3e^{\varsigma })}{\mathsf{\Gamma }(2\alpha +1)}+\frac{24{\left(2\right)}^{2\alpha }{e}^{\varsigma }{\tau }^{\alpha }\mathsf{\Gamma }(\frac{1}{2}+\alpha )}{\sqrt{\pi }\mathsf{\Gamma }(1+\alpha )\mathsf{\Gamma }(1+3\alpha )}\right)+\cdots ,\hfill \\ \hfill & \nu (\varsigma ,\tau )=-\frac{1}{2}+tanh\left(\frac{\varsigma }{2}\right)+\frac{1}{2}tan{h}^2\left(\frac{\varsigma }{2}\right)+\frac{2e^{\varsigma }(-3+e^{\varsigma }){\tau }^{\alpha }}{(1+e^{\varsigma }{)}^3\mathsf{\Gamma }(\alpha +1)}+\hfill \\ & \frac{2e^{\varsigma }{\tau }^{2\alpha }}{{(1+e^{\varsigma })}^6}(-\frac{96e^{\varsigma }{\tau }^{\alpha }}{{\mathsf{\Gamma }}^3(1+\alpha )}+\frac{(1+e^{\varsigma }(-3+e^{\varsigma }(25+e^{\varsigma }(-19+e^{\varsigma }))))}{\mathsf{\Gamma }(1+2\alpha )}\hfill \\ \hfill +& \frac{24e^{\varsigma }(1+e^{\varsigma })(-2+sinh\left(\varsigma \right))}{{\mathsf{\Gamma }}^2(1+\alpha )})+\cdots .\hfill \end{array}$$

(42)

In Figure 1, two-dimensional (2D) representations of the LRPSM and NIM solutions for $\mu (\varsigma ,\tau )$ and $\nu (\varsigma ,\tau )$ are displayed at different fractional-orders, with $\tau =0.1$. Figure 2 displays 3D profile of the LRPSM and NIM solutions for $\mu (\varsigma ,\tau )$ and $\nu (\varsigma ,\tau )$ with varying fractional-orders. In Figure 3, 3D representations of the LRPSM and NIM solutions for $\mu (\varsigma ,\tau )$ are shown at different fractional-orders. The comparison between the NIM and LRPSM solutions, with the absolute error, for both $\mu (\varsigma ,\tau )$ and $\nu (\varsigma ,\tau )$ can be found in

Table 1. The numerical results of $\mu (\varsigma ,\tau )$.

Comparison of the NIM and LRPSM Solutions with Absolute Error for $\mathit{\mu }(\mathit{\varsigma },\mathit{\tau })$
$t=0.095$	$eta$	NIM Solution	LRPSM Solution	Exact Solution	NIM Error	LRPSM Error	Generalized Taylor series
	1.00	−0.39842	−0.38887	−0.38884	0.009557	0.0000290375	0.0000280395
	1.20	−0.35982	−0.35121	−0.3512	0.008622	0.0000109524	0.0000108523
	1.40	−0.32044	−0.31277	−0.31277	0.00767	0.0000034773	0.0000034883
	1.60	−0.28182	−0.27507	−0.27509	0.006738	0.0000139555	0.0000139654
	1.80	−0.24515	−0.23927	−0.23929	0.005855	0.0000207583	0.0000207789
	2.00	−0.21121	−0.20615	−0.20617	0.00504	0.0000244851	0.0000244978
$t=0.010$	1.00	−0.39843	−0.38887	−0.38884	0.009557	0.0000290375	0.0000290375
	1.20	−0.35982	−0.35121	−0.3512	0.008622	0.0000109524	0.0000109524
	1.40	−0.32044	−0.31277	−0.31277	0.00767	0.0000034779	0.0000034779
	1.60	−0.28182	−0.27507	−0.27509	0.006738	0.0000139555	0.0000139555
	1.80	−0.24515	−0.23927	−0.23929	0.005855	0.0000207583	0.0000207583
	2.00	−0.21121	−0.20615	−0.20617	0.00504	0.0000244851	0.0000244851

and

Table 2. The numerical results of $\nu (\varsigma ,\tau )$.

Comparison of NIM and LRPSM Solutions with Absolute Error for $\mathit{\nu }(\mathit{\varsigma },\mathit{\tau })$
$t=0.095$	$eta$	NIM Solution	LRPSM Solution	Exact Solution	NIM Error	LRPSM Error	Generalized Taylor series
	2.0	0.551882	0.551883	0.553454	0.001571810	0.001570570	0.001570570
	2.2	0.621168	0.621170	0.622512	0.001343914	0.001342650	0.001342650
	2.4	0.681399	0.681401	0.682546	0.001140890	0.001139630	0.001139630
	2.6	0.733241	0.733241	0.734202	0.000962688	0.000961455	0.000961455
	2.8	0.777484	0.777485	0.778292	0.000808172	0.000807004	0.000807004
	3.0	0.814978	0.814979	0.815654	0.000675545	0.000674467	0.000674467
$t=0.010$	1.00	−0.39843	−0.38887	−0.38884	0.009557	0.0000290375	0.0000290375
	2.0	0.552432	0.552443	0.557129	0.00469676	0.004685690	0.004685690
	2.2	0.621710	0.621722	0.625723	0.00401319	0.004001820	0.004001820
	2.4	0.681911	0.681922	0.685316	0.00340511	0.003393710	0.003393710
	2.6	0.733708	0.733719	0.736580	0.00287199	0.002860890	0.002860890
	2.8	0.777903	0.777913	0.780313	0.00241016	0.002399650	0.002399650
	3.0	0.815346	0.815356	0.817360	0.00201405	0.002004340	0.002004340

, respectively.

6. Conclusions

The fractional-order nonlinear Belousov–Zhabotinsky system has been analyzed via two analytical approaches known as the hybrid residual power series method with the Laplace transformation and a novel iterative technique. The analytical solution of the given problem was calculated and compared with obtained solutions using the proposed techniques. It was observed from the numerical examples that the obtained results were completely identical to the exact solutions. Actual examples demonstrated the accuracy of the suggested methods. Moreover, the suggested methods are characterized by being highly efficient with fewer calculations. Furthermore, the suggested approaches can easily be widely used for resolving various fractional-order partial differential equation nonlinear systems. Finally, the proposed methods can be used to interpret and analyze many non-linear phenomena that arises in plasma physics, such as soliton waves, rogue waves, shock waves, etc. [10,11,12,13,14,15,16,17,18,19].