Exploring new routes in fractional modeling: analytical solutions of Burgers-type systems via Caputo-Hadamard and ϕ-Caputo derivatives

Abstract

The study of fractional-order partial differential equations has gained significant attention due to its ability to model complex physical systems with memory effects and hereditary properties. Among these systems is the Burgers equation, which serves as a fundamental model for describing nonlinear waves, turbulence, and the behavior of viscous fluids. Recent advancements in fractional calculus have led to the development of generalized fractional operators, such as the Caputo-Hadamard and ϕ-Caputo fractional derivatives. These operators offer greater flexibility and precision in capturing temporal and spatial nonlocal effects. This paper provides a comprehensive review of analytical and semi-analytical methods for solving these systems, with a particular emphasis on the Residual Power Series Method (RPSM) and the New Iterative Method (NIM). Both methods demonstrate superior convergence rates and accuracy. By combining generalized fractional operators with iterative algorithms, new approaches can be developed to model real-world phenomena in fields like physics, engineering, and applied mathematics. This work not only summarizes recent progress but also establishes a strong foundation for future research into complex fractional-order systems.

Introduction

Fractional calculus has become an effective mathematical method in the modeling of real-world situations with memory and invariances. The significance of the fractional calculus theory is indeed noteworthy especially for nonlinear partial differential equation (PDEs), being the traditional integer-order models unable to grasp the complex behavior that is typical of many real world phenomena [1, 2]. Among these nonlinear PDEs, system of Burgers equation has attracted much attention because it is a well-known PDE model which describes various motion and phenomena in the fields of fluid dynamics, gas dynamics, traffic flow, turbulence, and acoustics [3–5]. The classical Burgers equation, which was initially proposed in 20 century, is a basic simplification of the Navier-Stokes equations, and is given in its standard form as:

1

$\begin{array}{rl}& {\mathcal{D}}_{\mathfrak{Q}}^{\varsigma ;\varphi }\mathrm{\Omega }+\mathrm{\Omega }\frac{\partial \mathrm{\Omega }}{\partial \mathfrak{P}}-\frac{{\partial }^2\mathrm{\Omega }}{\partial {\mathfrak{P}}^2}=0,\text{ where }0<\varsigma \le 1.\end{array}$ This equation includes both nonlinear convection and linear diffusion, and therefore serves as the example of a “first-order model” for the formation and dissipation of shock waves [6]. The Burgers system deals with this extension when considering coupled systems of multiple interacting equations, which is widely used kinetic models that describe more complicated interactions in multi-physical aspects [7, 8].

The increasing popularity for fractional-order modeling has prompted scientists to extend the Burgers equation using different fractional derivatives definitions. Among them, the Caputo-Hadamard fractional derivative has attracted much attention, since it can handle singular kernels and remains physical relevant by using it as a nonlocal operator [9–11]. With the help of the Caputo sense of initial condition treatment and Hadamard framework for logarithmic-based scaling, the derivative enables model to efficiently and effectively model complex temporal behaviors [12, 13].

Recent developments in fractional calculus have proposed the use of the generalized operator called the ψ-Hilfer derivative, which is a more general operator that offers a generalized framework that brings together and extends a variety of other existing fractional derivatives [14]. It has been developed further theoretically, including the formulation of a type of Leibniz rule involving the ψ-Hilfer operator, to give it a broader range of applications to nonlinear problems [15]. Based on them, the analysis of fractional differential equations, such as the fractional p-logistic equation, and its derivatives has provided important understanding of the existence and nonexistence of positive solutions, which illustrates the more general applicability of generalized fractional operators in nonlinear analysis [16]. Most recently, the introduction of the ϕ-Caputo fractional derivative has increased the set of tools available for the description of nonlinear dynamics [17, 18]. This operator contains a function $\varphi (\mathfrak{Q})$, which serves as a kernel modifying function, and can be chosen to provide a more flexible and accurate description of memory effects. The utility of this operator in tackling different types of differential equations has also been investigated by some other researchers which further shows its highly adaptive feature and high effectiveness for complex modeling problems [19, 20].

For analytical or semi-analytical solutions for such fractional-order systems, the Residual Power Series Method (RPSM) is one of the promising methods. Jointly developed as an alternative to perturbative solutions, RPSM seeks the solution iteratively and refines it at each iteration by reducing the residual error in each approximation [21–23]. It has been widely used for various linear and nonlinear problems, and can provide high accuracy and rapid convergence under some conditions (e.g., see [24, 25]).

New iterative methods (NIMs) with RPSM have recently come into popularity as well. These schemes are based on the scheme of successive approximations, expanding on trial functions and correcting them using residual functions or correction functionals [26–28]. They offer strong tools to handle the nonlinearities and memory nature of fractional PDEs. The application of these techniques to the fractional Burgers type systems has been illustrated in several recent studies and shown to provide a promising answer in terms of solution behaviour, convergence and computational costs [29–31].

The embedding of the Burgers system with modern type of fractional operators i.e. Caputo- Hadamard and ϕ-Caputo derivative brings with itself new mathematical challenges especially in connection with available analytics. Thus, the importance of the modern iterative and semi-analytical techniques is inevitable to construct accurate solutions. This field has expanded rapidly in recent years, underlying the theoretical aspects of fractional calculus and the application of numerical or semi-numerical algorithms [32–34]. Consider the nonlinear system of Burger’s equation:

2

$\begin{array}{rl}& {\mathcal{D}}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\mathrm{\Omega }}_1-\frac{{\partial }^2{\mathrm{\Omega }}_1}{\partial {\mathfrak{P}}^2}-2{\mathrm{\Omega }}_1\frac{\partial {\mathrm{\Omega }}_1}{\partial \mathfrak{P}}+{\mathrm{\Omega }}_2\frac{\partial {\mathrm{\Omega }}_1}{\partial \mathfrak{P}}+{\mathrm{\Omega }}_1\frac{\partial {\mathrm{\Omega }}_2}{\partial \mathfrak{P}}=0,\\ & {\mathcal{D}}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\mathrm{\Omega }}_2-\frac{{\partial }^2{\mathrm{\Omega }}_2}{\partial {\mathfrak{P}}^2}-2{\mathrm{\Omega }}_2\frac{\partial {\mathrm{\Omega }}_2}{\partial \mathfrak{P}}+{\mathrm{\Omega }}_2\frac{\partial {\mathrm{\Omega }}_1}{\partial \mathfrak{P}}+{\mathrm{\Omega }}_1\frac{\partial {\mathrm{\Omega }}_2}{\partial \mathfrak{P}}=0,\text{ where }0<\varsigma \le 1.\end{array}$ Several researchers, for example, Odibat, Abdeljawad, Atangana based much of their research on the foundational theories of generalized fractional operators and applications to nonlinear problems [35–37]. They lay the foundations for modeling, stability, and solution structure in fractional PDEs. The implications of these methods are not just theoretical, they are instead expanded into real-world aspects such as viscoelastic materials and anomalous diffusions and bioengineering systems [38–40]. Moreover, the complementary characteristics of RPSM with other techniques, such as Variational Iteration Method (VIM), Adomian Decomposition Method (ADM) and Homotopy Analysis Method (HAM), have been demonstrated to be promising [41–43]. Hybrid approaches generally accelerate the convergence and improve the accuracy, and are useful in multi-term ones or coupled systems. Furthermore, comparisons of the RPSM with finite element, finite difference, and spectral methods have demonstrated the efficiency and simplicity of RPSM [44, 45].

The current paper contributes to the analysis of the fractional-order Burgers equations by using the RPSM, and novel iterative algorithms that are especially useful in solving analytical problems and numerical issues [46–49]. The given methodology, based on recent advances in adaptive minimization of residuals schemes and symbolic computation software, makes it possible to conduct detailed, parametric studies and sensitivity analyses of any fractional model. Our primary contribution is that we derive the analytical solutions of the Burgers system with ϕ-caputo fractional derivative, thus generalizing the classical integer-order theory to more complicated nonlocal theories. Two interesting benefits of this extension are the memory and hereditary effects which are absent at the classical level are naturally captured by the fractional operators, and the additional degrees of freedom available through the ϕ-Caputo kernels give the models more freedom to match complex dynamical behaviour. In addition to confirming the superiority of the suggested fractional framework over classical methods, these contributions give a clear methodology to describe and study the nonlinear fractional systems, both analytically and numerically.

Preliminaries

In this section, we review key properties and definitions from the theory of fractional calculus, which are utilized throughout this paper.

Definition 2.1

[50, 51] Assume the function φ defined on $I=[a,b]$, in relation to function ϕ of an order ς, such that $\varphi \in C^1(I)$ is a growing function, ${\varphi }^{\prime }(\mathfrak{Q})\ne 0$ for each $\mathfrak{Q}\in I$, and $\varsigma >0$. The following is the definition of the left fractional integral of φ:

3

$\begin{array}{r}{\mathcal{J}}_a^{\varsigma ;\varphi }[\phi (\mathfrak{P},\mathfrak{Q})]=\frac{1}{\mathrm{\Gamma }(\varsigma )}{\int }_a^{\mathfrak{Q}}{\varphi }^{\prime }(\sigma ){(\varphi (\mathfrak{Q})-\varphi (\sigma ))}^{\varsigma -1}\phi (\mathfrak{P},\sigma )\mathrm{d}\sigma .\end{array}$

4

$\begin{array}{r}{\mathcal{J}}_a^{0;\varphi }[\phi (\mathfrak{P},\mathfrak{Q})]=\phi (\mathfrak{P},\mathfrak{Q}).\end{array}$

Theorem 2.1

([51])

Let$\varsigma >0,m\in \mathbb{N},I$be the interval$-\mathrm{\infty }\le a<\mathrm{\infty }$, and$\phi ,\varphi \in C^m(I)$two functions, whereϕis increasing and${\varphi }^{\prime }(\mathfrak{Q})\ne 0$, for all$\mathfrak{Q}\in I$. The leftϕ-Caputo fractional derivative ofφof an orderςis given by:

5

$\begin{array}{r}{}_a{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }[\phi (\mathfrak{P},\mathfrak{Q})]={\mathcal{J}}_a^{m-\varsigma ;\varphi }{\left(\frac{1}{{\varphi }^{\prime }(\mathfrak{Q})}\frac{\partial }{\partial \mathfrak{Q}}\right)}^m\phi (\mathfrak{P},\mathfrak{Q}),\end{array}$such that$m=[\varsigma ]+1$for$\varsigma \notin \mathbb{N},m=\varsigma$for$\varsigma \in \mathbb{N}$.

In order to simplify the notation, we employ the abbreviated form.

6

$\begin{array}{r}{\phi }^{[m];\varphi }(\mathfrak{P},\mathfrak{Q})={\left(\frac{1}{{\varphi }^{\prime }(\mathfrak{Q})}\frac{\partial }{\partial \mathfrak{Q}}\right)}^m\phi (\mathfrak{P},\mathfrak{Q}).\end{array}$$\varphi (\mathfrak{Q})=\mathfrak{Q}$ yields the Caputo fractional derivative in [52], whereas $\varphi (\mathfrak{Q})=ln(\mathfrak{Q})$ yields the Caputo-Hadamard fractional derivative in [53].

Theorem 2.2

([53])

Let$\mathrm{\Re }(\varsigma )\ge 0,m=[\mathrm{\Re }(\varsigma )]+1$. If$\phi \in C_{\delta }^m([u,v])$, with$(0<u<v<\mathrm{\infty })$, and

7

$\begin{array}{r}A{C}_{\delta }^m([u,v])=\{g:[u,v]\to \mathbb{C}:{\delta }^{m-1}g(\mathfrak{P})\in AC[u,v],\delta =\mathfrak{P}\frac{\mathrm{d}}{\mathrm{d}\mathfrak{P}}\},\end{array}$then${}_u{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;ln(\mathfrak{Q})}\phi (\mathfrak{Q})$exist everywhere on$[u,v]$:

8

$\begin{array}{r}{}_u{}^{\mathcal{CH}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma }\phi (\mathfrak{Q})=\frac{1}{\mathrm{\Gamma }(m-\varsigma )}{\int }_a^{\mathfrak{Q}}{(ln\frac{\mathfrak{Q}}{\xi })}^{m-\varsigma -1}{\left(\mathfrak{Q}\frac{\mathrm{d}}{\mathrm{d}\mathfrak{Q}}\right)}^m\phi (\xi )\frac{\mathrm{d}\xi }{\xi }.\end{array}$

Lemma 2.3

([53])

Let$\phi \in A{C}_{\delta }^m([u,v])$, and$\varsigma \in \mathbb{C}$, then

9

$\begin{array}{r}{\mathcal{J}}_u^{\varsigma ;ln(\mathfrak{Q})}({}_u{}^{\mathcal{CH}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma }\phi (\mathfrak{Q}))=\phi (\mathfrak{Q})-\sum _{k=0}^{n-1}\frac{{\delta }^k\phi (a)}{k!}{(ln\frac{\mathfrak{Q}}{u})}^k.\end{array}$

Theorem 2.4

([53])

Let$\phi \in C^m([u,v])$, and$m-1<\varsigma \le m,m\in \mathbb{N},\tau >0$, then

1. If$\phi (\mathfrak{Q})={(\varphi (\mathfrak{Q})-\varphi (u))}^{\tau -1}$, then

   10

   $\begin{array}{r}{\mathcal{J}}_u^{\varsigma ;\varphi }[\phi (\mathfrak{Q})]=\frac{\mathrm{\Gamma }(\varsigma )}{\mathrm{\Gamma }(\varsigma +\tau )}{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma +\tau -1}.\end{array}$

2. 11

   $\begin{array}{r}{\mathcal{J}}_u^{\varsigma ;\varphi }({}_u{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }[\phi (\mathfrak{P},\mathfrak{Q})])=\phi (\mathfrak{P},\mathfrak{Q})-\sum _{k=0}^{n-1}\frac{{\phi }^{[k];\varphi }(\mathfrak{P},u)}{k!}{(\varphi (\mathfrak{Q})-\varphi (u))}^k.\end{array}$

Residual power series (RPS) method

In this section, we discuss the RPS method through a new approach.

Definition 3.1

A power series representation of the form

12

$\begin{array}{r}\sum _{m=0}^{\mathrm{\infty }}{H}_m{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }=C_0+C_1{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }+C_2{(\varphi (\mathfrak{Q})-\varphi (u))}^{2\varsigma }+\cdots ,\end{array}$ is called a fractional power series around a, such that $\mathfrak{Q}$ is a variable, $H_m$s are constants called the coefficients of the series, where $0\le n-1<\varsigma \le n,n\in \mathbb{N}$, and $\varphi (\mathfrak{Q})\ge \varphi (u)$.

Theorem 3.1

Assume that the fractional power series (FPS) notation ofφat$\varphi (u)$has the following form.

13

$\begin{array}{r}\phi (\mathfrak{Q})=\sum _{m=0}^{\mathrm{\infty }}{H}_m{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma },\end{array}$where$0\le n-1<\varsigma \le n,\varphi (u)\le \varphi (\mathfrak{Q})\le \varphi (u)+\mathcal{R}$.

If$\phi (\mathfrak{Q}),{}_u{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{m\varsigma ;\varphi }\phi (\mathfrak{Q})\in C[\varphi (u),\varphi (u)+\mathcal{R}]$where${}_u{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{m\varsigma ;\varphi }\phi \in C[\varphi (u),\varphi (u)+\mathcal{R}]$for$m=1,2,3,\dots$ , then the coefficients$H_m$are given as:

14

$\begin{array}{r}{H}_m=\frac{{}_u{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{m\varsigma ;\varphi }\phi (u)}{\mathrm{\Gamma }(m\varsigma +1)},\end{array}$where${}_u{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{m\varsigma ;\varphi }={\underbrace{{}_u{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }\circ \cdots {\circ }_u^{\mathcal{C}}{\mathcal{D}}_{\mathfrak{Q}}^{\varsigma ;\varphi }}}_{m\mathit{\text{ times}}}$, with$\mathcal{R}$is the radius of convergence.

Theorem 3.2

A power series of the form${\sum }_{m=0}^{\mathrm{\infty }}{\phi }_m(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }$is called a multiple FPS about$\varphi (\mathfrak{Q})=\varphi (u)$of the form

15

$\begin{array}{r}\mathrm{\Omega }(\mathfrak{P},\mathfrak{Q})=\sum _{m=0}^{\mathrm{\infty }}{\phi }_m(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma },\mathfrak{P}\in I,\varphi (u)\le \varphi (\mathfrak{Q})\le \varphi (u)+\mathcal{R}.\end{array}$If${}_u{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{m\varsigma ;\varphi }\mathrm{\Omega }(\mathfrak{P},\mathfrak{Q}),m=0,1,2,3,\dots$are continuous on$I\times (\varphi (u),\varphi (u)+\mathcal{R})$, then

16

$\begin{array}{r}{\phi }_m(\mathfrak{P})=\frac{{}_u{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{m\varsigma ;\varphi }\mathrm{\Omega }(\mathfrak{P},{\varphi }^{-1}(u))}{\mathrm{\Gamma }(m\varsigma +1)},\end{array}$where${}_u{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{m\varsigma ;\varphi }={\underbrace{{}_u{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }\circ \cdots {\circ }_u^{\mathcal{C}}{\mathcal{D}}_{\mathfrak{Q}}^{\varsigma ;\varphi }}}_{m\mathit{\text{ times}}}$, and$\mathcal{R}={min}_{\mathcal{C}\in I}{\mathcal{R}}_H$with${\mathcal{R}}_H$is the radius of convergence of the FPS${\sum }_{m=0}^{\mathrm{\infty }}{\phi }_m(H){(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }$.

It is evident from the last theorem that $m+1$ dimensional function can be obtained in the same way as the following corollary $\varphi (\mathfrak{Q})=\varphi (u)$.

Corollary 3.3

Assume that the multiple FPS notation of$\mathrm{\Omega }(\mathfrak{P},\mathfrak{Q})$at$\varphi (\mathfrak{Q})=\varphi (u)$has the following form:

17

$\begin{array}{r}\mathrm{\Omega }(\mathfrak{P},\mathfrak{Q})=\sum _{m=0}^{\mathrm{\infty }}{\phi }_m(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma },\end{array}$

18

$\begin{array}{r}\mathfrak{P}\in I_1\times I_2,\varphi (u)\le \varphi (\mathfrak{Q})\le \varphi (u)+\mathcal{R}.\end{array}$If${}_u{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{m\varsigma ;\varphi }\mathrm{\Omega }(\mathfrak{P},\mathfrak{Q}),m=0,1,2,3,\dots$are continuous on$I_1\times I_2\times (\varphi (u),\varphi (u)+\mathcal{R})$, then

19

$\begin{array}{r}{\phi }_m(\mathfrak{P})=\frac{{}_u{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{m\varsigma ;\varphi }\mathrm{\Omega }(\mathfrak{P},{\varphi }^{-1}(u))}{\mathrm{\Gamma }(m\varsigma +1)},\end{array}$where${}_u{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{m\varsigma ;\varphi }={\underbrace{{}_u{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }\circ \cdots {\circ }_u^{\mathcal{C}}{\mathcal{D}}_{\mathfrak{Q}}^{\varsigma ;\varphi }}}_{m\mathit{\text{ times}}}$, and$\mathcal{R}={min}_{\mathcal{C}\in I_1\times I_2}{\mathcal{R}}_{H,K}$in which${\mathcal{R}}_{H,K}$is the radius of convergence of the FPS${\sum }_{m=0}^{\mathrm{\infty }}{\phi }_m(H,K){(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }$.

Generalization of the solutions of the fractional system of PDEs using the ϕ-Caputo fractional derivative

In this section, we apply the RPS method for solving the system of nonlinear PDEs of time fractional ϕ-Caputo fractional derivative in the following form:

20

$\begin{array}{rl}& {}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\mathrm{\Omega }}_1+\mathcal{R}{\mathrm{\Omega }}_1+\mathcal{N}{\mathrm{\Omega }}_1+{\mathbb{H}}_1=0,\\ & {}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\mathrm{\Omega }}_2+\mathcal{R}{\mathrm{\Omega }}_2+\mathcal{N}{\mathrm{\Omega }}_2+{\mathbb{H}}_2=0.\end{array}$ The starting criteria are satisfied by ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$, and they can be rewritten as:

21

$\begin{array}{rl}{\mathrm{\Omega }}_1(\mathfrak{P},0)& =f(\mathfrak{P}),\\ {\mathrm{\Omega }}_2(\mathfrak{P},0)& =g(\mathfrak{P}).\end{array}$ As a result, we may get the first estimate of ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$ as:

22

$\begin{array}{rl}{{\mathrm{\Omega }}_1}_0(\mathfrak{P},0)& =f_0(\mathfrak{P})=f(\mathfrak{P}),\\ {{\mathrm{\Omega }}_2}_0(\mathfrak{P},0)& =g_0(\mathfrak{P})=g(\mathfrak{P}).\end{array}$ Therefore, the approximate series solution of Eq. (20) could be rewritten as:

23

$\begin{array}{rl}{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})& =f(\mathfrak{P})+\sum _{m=1}^{\mathrm{\infty }}{f}_m(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)},\\ {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})& =g(\mathfrak{P})+\sum _{m=1}^{\mathrm{\infty }}{g}_m(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)}.\end{array}$ We use ${{\mathrm{\Omega }}_1}_i$ and ${{\mathrm{\Omega }}_2}_i$ to represent the i-th truncated series of ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$, respectively, in the next step:

24

$\begin{array}{rl}{{\mathrm{\Omega }}_1}_i(\mathfrak{P},\mathfrak{Q})& =f(\mathfrak{P})+\sum _{m=1}^i{f}_m(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)},\\ {{\mathrm{\Omega }}_2}_i(\mathfrak{P},\mathfrak{Q})& =g(\mathfrak{P})+\sum _{m=1}^i{g}_m(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)},\end{array}$ for $i=1,2,3,\dots$ For Eq. (20), we define the residual functions ${\text{Res}}_{{\mathrm{\Omega }}_1}$ and ${\text{Res}}_{{\mathrm{\Omega }}_2}$ as follows:

25

$\begin{array}{rl}{\text{Res}}_{{\mathrm{\Omega }}_1}& ={}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})+\mathcal{R}{\mathrm{\Omega }}_1+\mathcal{N}{\mathrm{\Omega }}_1+{\mathbb{H}}_1,\\ {\text{Res}}_{{\mathrm{\Omega }}_2}& ={}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})+\mathcal{R}{\mathrm{\Omega }}_2+\mathcal{N}{\mathrm{\Omega }}_2+{\mathbb{H}}_2.\end{array}$ The i-th truncated residual functions are thus

26

$\begin{array}{rl}{\text{Res}}_{{{\mathrm{\Omega }}_1}_i}& ={}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{{\mathrm{\Omega }}_1}_i(\mathfrak{P},\mathfrak{Q})+\mathcal{R}{{\mathrm{\Omega }}_1}_i+\mathcal{N}{{\mathrm{\Omega }}_1}_i+{{\mathbb{H}}_1}_i,\\ {\text{Res}}_{{{\mathrm{\Omega }}_2}_i}& ={}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{{\mathrm{\Omega }}_2}_i(\mathfrak{P},\mathfrak{Q})+\mathcal{R}{{\mathrm{\Omega }}_2}_i+\mathcal{N}{{\mathrm{\Omega }}_2}_i+{{\mathbb{H}}_2}_i.\end{array}$ Based on [50, 54, 55], ${lim}_{i\to \mathrm{\infty }}{\text{Res}}_i=\text{Res}(\mathfrak{P},\mathfrak{Q})$, and $\text{Res}(\mathfrak{P},\mathfrak{Q})=0$ for each $\varphi (\mathfrak{Q})\in [\varphi (u),\varphi (u)+\mathcal{R}]$ and $\mathfrak{P}\in \mathbb{R}$, with $\mathcal{R}$ is a non-negative real number representing the radius of convergence. Hence, ${}_a{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{k\varsigma ;\varphi }\text{Res}(\mathfrak{P},\mathfrak{Q})=0$. Given that a constant function’s fractional derivative in the Caputo sense is zero, the fractional derivative ${}_a{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{k\varsigma ;\varphi }$ of $\text{Res}(\mathfrak{P},\mathfrak{Q})$ and ${\text{Res}}_i(\mathfrak{P},\mathfrak{Q})$ correspond to $\varphi (\mathfrak{Q})=\varphi (u)$ for each $i=0,1,2,\dots$ .

If we set $\varphi (u)=0$, and $r=i-1$, we get

27

$\begin{array}{rl}{}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{(i-1)\varsigma ;\varphi }{\text{Res}}_{{{\mathrm{\Omega }}_1}_i}(\mathfrak{P},0)& =0,\\ {}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{(i-1)\varsigma ;\varphi }{\text{Res}}_{{{\mathrm{\Omega }}_2}_i}(\mathfrak{P},0)& =0.\end{array}$ Now, we use the RPS technique to obtain the form of the coefficients $f_m(\mathfrak{P})$, and $g_m(\mathfrak{P})$, where $m=1,2,3,\dots ,i$.

First, we enter the i-th shortened ${\mathrm{\Omega }}_1,{\mathrm{\Omega }}_2$, series into Eq. (26). Second, we determine the formula for the fractional derivative of ${}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{(i-1)\varsigma ;\varphi }$ for both ${\text{Res}}_{{{\mathrm{\Omega }}_1}_i}(\mathfrak{P},\mathfrak{Q})$ and ${\text{Res}}_{{{\mathrm{\Omega }}_2}_i}(\mathfrak{P},\mathfrak{Q})$, where $i=1,2,3,\dots$ Last, we solve the algebraic system (26) that was acquired.

Solution of Burger equation through ϕ-Caputo RPS method

28

$\begin{array}{rl}& {\mathcal{D}}_{\mathfrak{Q}}^{\varsigma ;\varphi }\mathrm{\Omega }+\mathrm{\Omega }\frac{\partial \mathrm{\Omega }}{\partial \mathfrak{P}}-\frac{{\partial }^2\mathrm{\Omega }}{\partial {\mathfrak{P}}^2}=0,\\ & \text{ where }0<\varsigma \le 1.\end{array}$ Starting constrant

29

$\begin{array}{rl}\mathrm{\Omega }(\mathfrak{P},0)& =\frac{1}{2}-\frac{1}{2}tanh\left(\frac{\mathfrak{P}}{4}\right).\end{array}$ By using RPS technique, we suppose the following series solution to the problem.

30

$\begin{array}{rl}\mathrm{\Omega }(\mathfrak{P},\mathfrak{Q})& =\sum _{m=0}^{\mathrm{\infty }}{f}_m(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)},\end{array}$ The first guess for $m=0$ may be obtained with starting condition Eq. (29), we get

31

$\begin{array}{rl}\mathrm{\Omega }(\mathfrak{P},\mathfrak{Q})& =\frac{1}{2}-\frac{1}{2}tanh\left(\frac{\mathfrak{P}}{4}\right)+\sum _{m=1}^{\mathrm{\infty }}{f}_m(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)}.\end{array}$ The following is the construction of the series of recommended solutions:

32

$\begin{array}{rl}{\mathrm{\Omega }}_i(\mathfrak{P},\mathfrak{Q})& =\frac{1}{2}-\frac{1}{2}tanh\left(\frac{\mathfrak{P}}{4}\right)+\sum _{m=1}^i{f}_m(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)}.\end{array}$ Defining the residual function as:

33

$\begin{array}{rl}{\text{Res}}_{\mathrm{\Omega }}& ={}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }\mathrm{\Omega }(\mathfrak{P},\mathfrak{Q})+\mathrm{\Omega }\frac{\partial \mathrm{\Omega }}{\partial \mathfrak{P}}-\frac{{\partial }^2\mathrm{\Omega }}{\partial {\mathfrak{P}}^2}.\end{array}$ The truncated residual functions for the i-th term are given as:

34

$\begin{array}{rl}{\text{Res}}_{{\mathrm{\Omega }}_i}& ={}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\mathrm{\Omega }}_i(\mathfrak{P},\mathfrak{Q})+{\mathrm{\Omega }}_i\frac{\partial {\mathrm{\Omega }}_i}{\partial \mathfrak{P}}-\frac{{\partial }^2{\mathrm{\Omega }}_i}{\partial {\mathfrak{P}}^2}.\end{array}$ Put Eq. (32) in Eq. (34), we obtain

35

$\begin{array}{rl}{\text{Res}}_{{\mathrm{\Omega }}_i}=& {}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }\left(\sum _{m=1}^i{f}_m\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)}\right)\\ & +(\sum _{m=1}^i\frac{f_m(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}-\frac{1}{2}tanh\left(\frac{\mathfrak{P}}{4}\right)+\frac{1}{2})\\ & \times (\sum _{m=1}^i\frac{\frac{\partial }{\partial \mathfrak{P}}{f}_m(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}-\frac{1}{8}{\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right))\\ & -\sum _{m=1}^i\frac{\frac{{\partial }^2}{\partial {\mathfrak{P}}^2}{f}_m(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}-\frac{1}{16}tanh\left(\frac{\mathfrak{P}}{4}\right){\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right).\end{array}$ Put $i=1$ in Eq. (32) to obtain the first approximation.

36

$\begin{array}{rl}{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})& =\frac{1}{2}-\frac{1}{2}tanh\left(\frac{\mathfrak{P}}{4}\right)+f_1(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}.\end{array}$ Consequently, these equations may be substituted into Eq. (35) to get the residual truncated functions as follows:

37

$\begin{array}{rl}{\text{Res}}_{{\mathrm{\Omega }}_1}=& f_1+(\frac{f_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}-\frac{1}{2}tanh\left(\frac{\mathfrak{P}}{4}\right)+\frac{1}{2})\\ & \times (\frac{\frac{\partial }{\partial \mathfrak{P}}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}-\frac{1}{8}{\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right))\\ & -\frac{\frac{{\partial }^2}{\partial {\mathfrak{P}}^2}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}-\frac{1}{16}tanh\left(\frac{\mathfrak{P}}{4}\right){\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right).\end{array}$ Next, we compute the truncated residual functions at $\varphi (\mathfrak{Q})=0$ to get:

38

$\begin{array}{rl}{\text{Res}}_{{\mathrm{\Omega }}_1}& =f_1-\frac{1}{16}{\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right),\end{array}$ and by Eq. (27), we know that

39

$\begin{array}{r}{\text{Res}}_{{\mathrm{\Omega }}_1}(\mathfrak{P},0)=0.\end{array}$ Solving for $f_1$ to obtain the following result.

40

$\begin{array}{rl}{f}_1& =\frac{1}{16}{\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right).\end{array}$ The two terms estimated RPS solution is as follows:

41

$\begin{array}{rl}{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})& =\frac{1}{2}-\frac{1}{2}tanh\left(\frac{\mathfrak{P}}{4}\right)+\frac{1}{16}{\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right)\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}.\end{array}$ Put $i=2$ in Eq. (32) to get the next approximation.

42

$\begin{array}{rl}{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})& =\frac{1}{2}-\frac{1}{2}tanh\left(\frac{\mathfrak{P}}{4}\right)+f_1\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}+f_2\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}.\end{array}$ Consequently, by putting these equations into Eq. (35), the residual truncated functions that follow may be derived as:

43

$\begin{array}{rl}{\text{Res}}_{{\mathrm{\Omega }}_2}=& f_1+f_2\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & +\frac{f_1^{\prime }(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{2\mathrm{\Gamma }(\varsigma +1)}+\frac{f_1(\mathfrak{P})f_1^{\prime }(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }{(\varsigma +1)}^2}\\ & +\frac{f_2^{\prime }(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{2\mathrm{\Gamma }(2\varsigma +1)}\\ & +\frac{f_2(\mathfrak{P})f_1^{\prime }(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{3\varsigma }}{\mathrm{\Gamma }(\varsigma +1)\mathrm{\Gamma }(2\varsigma +1)}+\frac{f_1(\mathfrak{P})f_2^{\prime }(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{3\varsigma }}{\mathrm{\Gamma }(\varsigma +1)\mathrm{\Gamma }(2\varsigma +1)}\\ & +\frac{f_2(\mathfrak{P})f_2^{\prime }(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{4\varsigma }}{\mathrm{\Gamma }{(2\varsigma +1)}^2}\\ & -\frac{tanh\left(\frac{\mathfrak{P}}{4}\right)f_1^{\prime }(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{2\mathrm{\Gamma }(\varsigma +1)}-\frac{tanh\left(\frac{\mathfrak{P}}{4}\right)f_2^{\prime }(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{2\mathrm{\Gamma }(2\varsigma +1)}\\ & -\frac{f_1^{″}(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & -\frac{f_2^{″}(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}-\frac{f_1(\mathfrak{P}){\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{8\mathrm{\Gamma }(\varsigma +1)}\\ & -\frac{f_2(\mathfrak{P}){\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{8\mathrm{\Gamma }(2\varsigma +1)}\\ & -\frac{1}{16}{\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right).\end{array}$ Next, applying operator ${}_0{}^{\mathcal{C}}\mathcal{D}_t^{\varsigma ;\varphi }$ into Eq. (43), and putting $f_1$, $f_2$ and $\varphi (\mathfrak{Q})=0$, we get

44

$\begin{array}{rl}{}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\text{Res}}_{{\mathrm{\Omega }}_2}(\mathfrak{P},0)& =f_2-\frac{1}{64}tanh\left(\frac{\mathfrak{P}}{4}\right){\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right).\end{array}$ From Eq. (27), we can deduce the subsequent result.

45

$\begin{array}{rl}{}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\text{Res}}_{{\mathrm{\Omega }}_2}(\mathfrak{P},0)& =0.\end{array}$ Consequently, $f_2$ is written as:

46

$\begin{array}{rl}{f}_2& =\frac{1}{64}tanh\left(\frac{\mathfrak{P}}{4}\right){\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right).\end{array}$ The three terms approximated RPS solution is given as:

47

$\begin{array}{rl}{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})=& \frac{1}{2}-\frac{1}{2}tanh\left(\frac{\mathfrak{P}}{4}\right)+\frac{1}{16}{\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right)\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & +\frac{1}{64}tanh\left(\frac{\mathfrak{P}}{4}\right){\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right)\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}.\end{array}$ By using the same procedure for $i=3,4,5,\dots$ , we obtain the subsequent result.

The exact solution for $\varphi (\mathfrak{Q})=\mathfrak{Q}$, and $a\to 1$ is given as:

48

$\begin{array}{rl}\mathrm{\Omega }(\mathfrak{P},\mathfrak{Q})& =\frac{1}{2}-\frac{1}{2}tanh(\frac{\mathfrak{P}}{4}-\frac{\mathfrak{Q}}{8}).\end{array}$

Solution of system of Burger equation through ϕ-Caputo RPS method

49

$\begin{array}{rl}& {\mathcal{D}}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\mathrm{\Omega }}_1-\frac{{\partial }^2{\mathrm{\Omega }}_1}{\partial {\mathfrak{P}}^2}-2{\mathrm{\Omega }}_1\frac{\partial {\mathrm{\Omega }}_1}{\partial \mathfrak{P}}+{\mathrm{\Omega }}_2\frac{\partial {\mathrm{\Omega }}_1}{\partial \mathfrak{P}}+{\mathrm{\Omega }}_1\frac{\partial {\mathrm{\Omega }}_2}{\partial \mathfrak{P}}=0,\\ & {\mathcal{D}}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\mathrm{\Omega }}_2-\frac{{\partial }^2{\mathrm{\Omega }}_2}{\partial {\mathfrak{P}}^2}-2{\mathrm{\Omega }}_2\frac{\partial {\mathrm{\Omega }}_2}{\partial \mathfrak{P}}+{\mathrm{\Omega }}_2\frac{\partial {\mathrm{\Omega }}_1}{\partial \mathfrak{P}}+{\mathrm{\Omega }}_1\frac{\partial {\mathrm{\Omega }}_2}{\partial \mathfrak{P}}=0,\\ & \text{ where }0<\varsigma \le 1.\end{array}$ Starting constraints

50

$\begin{array}{rl}{\mathrm{\Omega }}_1(\mathfrak{P},0)& =cos(\mathfrak{P}),\\ {\mathrm{\Omega }}_2(\mathfrak{P},0)& =cos(\mathfrak{P}).\end{array}$ By using RPS technique, we suppose the following series solution to the problem.

51

$\begin{array}{rl}{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})& =\sum _{m=0}^{\mathrm{\infty }}{f}_m(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)},\\ {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})& =\sum _{m=0}^{\mathrm{\infty }}{g}_m(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)}.\end{array}$ The first guess for $m=0$ may be obtained with starting conditions Eq. (50), we get

52

$\begin{array}{rl}{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})& =cos(\mathfrak{P})+\sum _{m=1}^{\mathrm{\infty }}{f}_m(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)},\\ {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})& =cos(\mathfrak{P})+\sum _{m=1}^{\mathrm{\infty }}{g}_m(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)}.\end{array}$ The following is the construction of the series of recommended solutions:

53

$\begin{array}{rl}{{\mathrm{\Omega }}_1}_i(\mathfrak{P},\mathfrak{Q})& =cos(\mathfrak{P})+\sum _{m=1}^i{f}_m(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)},\\ {{\mathrm{\Omega }}_2}_i(\mathfrak{P},\mathfrak{Q})& =cos(\mathfrak{P})+\sum _{m=1}^i{g}_m(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)}.\end{array}$ Defining the residual functions as:

54

$\begin{array}{rl}{\text{Res}}_{{\mathrm{\Omega }}_1}& ={}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})-\frac{{\partial }^2{\mathrm{\Omega }}_1}{\partial {\mathfrak{P}}^2}-2{\mathrm{\Omega }}_1\frac{\partial {\mathrm{\Omega }}_1}{\partial \mathfrak{P}}+{\mathrm{\Omega }}_2\frac{\partial {\mathrm{\Omega }}_1}{\partial \mathfrak{P}}+{\mathrm{\Omega }}_1\frac{\partial {\mathrm{\Omega }}_2}{\partial \mathfrak{P}},\\ {\text{Res}}_{{\mathrm{\Omega }}_2}& ={}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})-\frac{{\partial }^2{\mathrm{\Omega }}_2}{\partial {\mathfrak{P}}^2}-2{\mathrm{\Omega }}_2\frac{\partial {\mathrm{\Omega }}_2}{\partial \mathfrak{P}}+{\mathrm{\Omega }}_2\frac{\partial {\mathrm{\Omega }}_1}{\partial \mathfrak{P}}+{\mathrm{\Omega }}_1\frac{\partial {\mathrm{\Omega }}_2}{\partial \mathfrak{P}}.\end{array}$ The truncated residual functions for the i-th term are given as:

55

$\begin{array}{rl}{\text{Res}}_{{{\mathrm{\Omega }}_1}_i}& ={}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{{\mathrm{\Omega }}_1}_i(\mathfrak{P},\mathfrak{Q})-\frac{{\partial }^2{{\mathrm{\Omega }}_1}_i}{\partial {\mathfrak{P}}^2}-2{{\mathrm{\Omega }}_1}_i\frac{\partial {{\mathrm{\Omega }}_1}_i}{\partial \mathfrak{P}}+{{\mathrm{\Omega }}_2}_i\frac{\partial {{\mathrm{\Omega }}_1}_i}{\partial \mathfrak{P}}+{{\mathrm{\Omega }}_1}_i\frac{\partial {{\mathrm{\Omega }}_2}_i}{\partial \mathfrak{P}},\\ {\text{Res}}_{{{\mathrm{\Omega }}_2}_i}& ={}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{{\mathrm{\Omega }}_2}_i(\mathfrak{P},\mathfrak{Q})-\frac{{\partial }^2{{\mathrm{\Omega }}_2}_i}{\partial {\mathfrak{P}}^2}-2{{\mathrm{\Omega }}_2}_i\frac{\partial {{\mathrm{\Omega }}_2}_i}{\partial \mathfrak{P}}+{{\mathrm{\Omega }}_2}_i\frac{\partial {{\mathrm{\Omega }}_1}_i}{\partial \mathfrak{P}}+{{\mathrm{\Omega }}_1}_i\frac{\partial {{\mathrm{\Omega }}_2}_i}{\partial \mathfrak{P}}.\end{array}$ Put Eq. (53) in Eq. (55), we obtain

56

$\begin{array}{rl}{\text{Res}}_{{{\mathrm{\Omega }}_1}_i}=& {}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }\left(\sum _{m=1}^i{f}_m\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)}\right)\\ & +(\sum _{m=1}^i\frac{g_m(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}+cos(\mathfrak{P}))\\ & \times (\sum _{m=1}^i\frac{\frac{\partial f_m(\mathfrak{P})}{\partial \mathfrak{P}}{(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}-sin(\mathfrak{P}))\\ & +(\sum _{m=1}^i\frac{f_m(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}+cos(\mathfrak{P}))\\ & \times (\sum _{m=1}^i\frac{\frac{\partial g_m(\mathfrak{P})}{\partial \mathfrak{P}}{(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}-sin(\mathfrak{P}))\\ & -2(\sum _{m=1}^i\frac{f_m(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}+cos(\mathfrak{P}))\\ & \times (\sum _{m=1}^i\frac{\frac{\partial f_m(\mathfrak{P})}{\partial \mathfrak{P}}{(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}-sin(\mathfrak{P}))\\ & -\sum _{m=1}^i\frac{\frac{{\partial }^2{f}_m(\mathfrak{P})}{\partial {\mathfrak{P}}^2}{(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}+cos(\mathfrak{P}),\end{array}$

57

$\begin{array}{rl}{\text{Res}}_{{{\mathrm{\Omega }}_2}_i}=& {}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }\left(\sum _{m=1}^i{g}_m\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{m\varsigma }}{\mathrm{\Gamma }(m\varsigma +1)}\right)\\ & +(\sum _{m=1}^i\frac{g_m(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}+cos(\mathfrak{P}))\\ & \times (\sum _{m=1}^i\frac{\frac{\partial f_m(\mathfrak{P}}{\partial \mathfrak{P}}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}-sin(\mathfrak{P}))\\ & +(\sum _{m=1}^i\frac{f_m(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}+cos(\mathfrak{P}))\\ & \times (\sum _{m=1}^i\frac{\frac{\partial g_m(\mathfrak{P})}{\partial \mathfrak{P}}{(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}-sin(\mathfrak{P}))\\ & -2(\sum _{m=1}^i\frac{g_m(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}+cos(\mathfrak{P}))\\ & \times (\sum _{m=1}^i\frac{\frac{\partial g_m(\mathfrak{P})}{\partial \mathfrak{P}}{(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}-sin(\mathfrak{P}))\\ & -\sum _{m=1}^i\frac{\frac{{\partial }^2{g}_m(\mathfrak{P})}{\partial {\mathfrak{P}}^2}{(\varphi (\varphi )-\varphi (a))}^{\varsigma m}}{\mathrm{\Gamma }(m\varsigma +1)}+cos(\mathfrak{P}).\end{array}$

Put $i=1$ in Eq. (53) to obtain the first approximation.

58

$\begin{array}{rl}{{\mathrm{\Omega }}_1}_1(\mathfrak{P},\mathfrak{Q})& =cos(\mathfrak{P})+f_1(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)},\\ {{\mathrm{\Omega }}_2}_1(\mathfrak{P},\mathfrak{Q})& =cos(\mathfrak{P})+g_1(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}.\end{array}$ The next step is to find the first residual functions by putting Eq. (58) in Eq. (56) and Eq. (57) in the manner described below:

59

$\begin{array}{rl}{\text{Res}}_{{{\mathrm{\Omega }}_1}_1}=& f_1+\frac{g_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }{(\varsigma +1)}^2}+\frac{f_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }{(\varsigma +1)}^2}\\ & -\frac{2f_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }{(\varsigma +1)}^2}\\ & -\frac{cos(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}-\frac{\frac{{\partial }^2}{\partial {\mathfrak{P}}^2}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & +\frac{f_1(\mathfrak{P})sin(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & +\frac{cos(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}-\frac{g_1(\mathfrak{P})sin(x){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & +cos(\mathfrak{P}),\end{array}$

60

$\begin{array}{rl}{\text{Res}}_{{{\mathrm{\Omega }}_2}_1}=& g_1+\frac{g_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }{(\varsigma +1)}^2}+\frac{f_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }{(\varsigma +1)}^2}\\ & +\frac{cos(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & -\frac{f_1(\mathfrak{P})sin(x){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}-\frac{2g_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }{(\varsigma +1)}^2}\\ & -\frac{cos(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & -\frac{\frac{{\partial }^2}{\partial {\mathfrak{P}}^2}{g}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}+\frac{g_1(\mathfrak{P})sin(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & +cos(\mathfrak{P}).\end{array}$ Next, we compute the truncated residual functions at $\varphi (\mathfrak{Q})=0$ to get:

61

$\begin{array}{rl}{\text{Res}}_{{{\mathrm{\Omega }}_1}_1}& =f_1+cos(\mathfrak{P}),\\ {\text{Res}}_{{{\mathrm{\Omega }}_2}_1}& =g_1+cos(\mathfrak{P}).\end{array}$ From Eq. (27), we can deduce the subsequent result.

62

$\begin{array}{r}{\text{Res}}_{{{\mathrm{\Omega }}_1}_1}(\mathfrak{P},0)=0,{\text{Res}}_{{{\mathrm{\Omega }}_2}_1}(\mathfrak{P},0)=0.\end{array}$ Solving for $f_1$ and $g_1$ to obtain the following result.

63

$\begin{array}{rl}{f}_1& =-cos(\mathfrak{P}),\\ g_1& =-cos(\mathfrak{P}).\end{array}$ The two terms estimated RPS solution is as follows:

64

$\begin{array}{rl}{{\mathrm{\Omega }}_1}_1(\mathfrak{P},\mathfrak{Q})& =cos(\mathfrak{P})-cos(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)},\\ {{\mathrm{\Omega }}_2}_1(\mathfrak{P},\mathfrak{Q})& =cos(\mathfrak{P})-cos(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}.\end{array}$

Put $i=2$ in Eq. (53) to get the next approximation.

65

$\begin{array}{rl}{{\mathrm{\Omega }}_1}_2(\mathfrak{P},\mathfrak{Q})& =cos(\mathfrak{P})+f_1\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}+f_2\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)},\\ {{\mathrm{\Omega }}_2}_2(\mathfrak{P},\mathfrak{Q})& =cos(\mathfrak{P})+g_1\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}+g_2\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}.\end{array}$ Consequently, these equations may be substituted into Eq. (56) and Eq. (57) to get the residual truncated functions as follows:

66

$\begin{array}{rl}{\text{Res}}_{{{\mathrm{\Omega }}_1}_2}=& f_1+f_2\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & +\frac{g_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }{(\varsigma +1)}^2}+\frac{f_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }{(\varsigma +1)}^2}\\ & +\frac{g_2(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{3\varsigma }}{\mathrm{\Gamma }(\varsigma +1)\mathrm{\Gamma }(2\varsigma +1)}\\ & +\frac{g_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{3\varsigma }}{\mathrm{\Gamma }(\varsigma +1)\mathrm{\Gamma }(2\varsigma +1)}+\frac{f_2(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{3\varsigma }}{\mathrm{\Gamma }(\varsigma +1)\mathrm{\Gamma }(2\varsigma +1)}\\ & +\frac{f_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{3\varsigma }}{\mathrm{\Gamma }(\varsigma +1)\mathrm{\Gamma }(2\varsigma +1)}\\ & +\frac{g_2(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{4\varsigma }}{\mathrm{\Gamma }{(2\varsigma +1)}^2}+\frac{f_2(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{4\varsigma }}{\mathrm{\Gamma }{(2\varsigma +1)}^2}\\ & -\frac{2f_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }{(\varsigma +1)}^2}\\ & -\frac{2f_2(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{3\varsigma }}{\mathrm{\Gamma }(\varsigma +1)\mathrm{\Gamma }(2\varsigma +1)}-\frac{2f_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{3\varsigma }}{\mathrm{\Gamma }(\varsigma +1)\mathrm{\Gamma }(2\varsigma +1)}\\ & -\frac{2f_2(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{4\varsigma }}{\mathrm{\Gamma }{(2\varsigma +1)}^2}\\ & -\frac{cos(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}-\frac{cos(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}\\ & -\frac{\frac{{\partial }^2}{\partial {\mathfrak{P}}^2}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & -\frac{\frac{{\partial }^2}{\partial {\mathfrak{P}}^2}{f}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}+\frac{f_2(\mathfrak{P})sin(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}\\ & +\frac{f_1(\mathfrak{P})sin(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & +\frac{cos(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}+\frac{cos(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}\\ & -\frac{g_1(\mathfrak{P})sin(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & -\frac{g_2(\mathfrak{P})sin(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}+cos(\mathfrak{P}),\end{array}$

67

$\begin{array}{rl}{\text{Res}}_{{{\mathrm{\Omega }}_2}_2}=& g_1+g_2\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & +\frac{g_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }{(\varsigma +1)}^2}+\frac{f_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }{(\varsigma +1)}^2}\\ & +\frac{g_2(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{3\varsigma }}{\mathrm{\Gamma }(\varsigma +1)\mathrm{\Gamma }(2\varsigma +1)}\\ & +\frac{g_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{3\varsigma }}{\mathrm{\Gamma }(\varsigma +1)\mathrm{\Gamma }(2\varsigma +1)}+\frac{f_2(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{3\varsigma }}{\mathrm{\Gamma }(\varsigma +1)\mathrm{\Gamma }(2\varsigma +1)}\\ & +\frac{f_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{3\varsigma }}{\mathrm{\Gamma }(\varsigma +1)\mathrm{\Gamma }(2\varsigma +1)}\\ & +\frac{g_2(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{4\varsigma }}{\mathrm{\Gamma }{(2\varsigma +1)}^2}+\frac{f_2(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{4\varsigma }}{\mathrm{\Gamma }{(2\varsigma +1)}^2}\\ & +\frac{cos(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ \begin{array}{rl}& +\frac{cos(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{f}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}-\frac{f_1(\mathfrak{P})sin(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & -\frac{f_2(\mathfrak{P})sin(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}\end{array}\\ & -\frac{2g_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }{(\varsigma +1)}^2}-\frac{2g_2(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{3\varsigma }}{\mathrm{\Gamma }(\varsigma +1)\mathrm{\Gamma }(2\varsigma +1)}\\ & -\frac{2g_1(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{3\varsigma }}{\mathrm{\Gamma }(\varsigma +1)\mathrm{\Gamma }(2\varsigma +1)}\\ & -\frac{2g_2(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_2(\mathfrak{P}){(\varphi (\mathfrak{P})-\varphi (a))}^{4\varsigma }}{\mathrm{\Gamma }{(2\varsigma +1)}^2}-\frac{cos(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & -\frac{cos(\mathfrak{P})\frac{\partial }{\partial \mathfrak{P}}{g}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}\\ & -\frac{\frac{{\partial }^2}{\partial {\mathfrak{P}}^2}{g}_1(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}-\frac{\frac{{\partial }^2}{\partial {\mathfrak{P}}^2}{g}_2(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}\\ & +\frac{g_1(\mathfrak{P})sin(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & +\frac{g_2(\mathfrak{P})sin(\mathfrak{P}){(\varphi (\mathfrak{Q})-\varphi (a))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}+cos(\mathfrak{P}).\end{array}$ Next, applying operator ${}_0{}^{\mathcal{C}}\mathcal{D}_t^{\varsigma ;\varphi }$ into Eq. (66) and Eq. (67), and put $f_1$, $g_1$ and $\varphi (\mathfrak{Q})=0$, we obtain

68

$\begin{array}{rl}{}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\text{Res}}_{{{\mathrm{\Omega }}_1}_2}(\mathfrak{P},0)& =f_2-cos(\mathfrak{P}),\\ {}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\text{Res}}_{{{\mathrm{\Omega }}_2}_2}(\mathfrak{P},0)& =g_2-cos(\mathfrak{P}).\end{array}$ Equation (27) as a foundation, we obtain

69

$\begin{array}{rl}{}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\text{Res}}_{{{\mathrm{\Omega }}_1}_2}(\mathfrak{P},0)& =0,\\ {}_0{}^{\mathcal{C}}\mathcal{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\text{Res}}_{{{\mathrm{\Omega }}_2}_2}(\mathfrak{P},0)& =0.\end{array}$ Solving Eq. (68) for $f_2$ and $g_2$, we get

70

$\begin{array}{rl}{f}_2& =cos(\mathfrak{P}),\\ g_2& =cos(\mathfrak{P}).\end{array}$ The three terms approximate series solution becomes

71

$\begin{array}{rl}{{\mathrm{\Omega }}_1}_2(\mathfrak{P},\mathfrak{Q})& =cos(\mathfrak{P})-cos(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}+cos(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)},\\ {{\mathrm{\Omega }}_2}_2(x,t)& =cos(\mathfrak{P})-cos(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}+cos(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)},\end{array}$ By using the same procedure for $i=3,4,5,\dots$ , we obtain the subsequent result. $\begin{array}{rl}{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})=& cos(\mathfrak{P})[1-\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}+\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}\\ & -\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{3\varsigma }}{\mathrm{\Gamma }(3\varsigma +1)}+\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{4\varsigma }}{\mathrm{\Gamma }(4\varsigma +1)}+\cdots ]\\ =& cos(\mathfrak{P})\sum _{k=0}^{\mathrm{\infty }}\frac{{(-1)}^k{(\varphi (\mathfrak{Q})-\varphi (u))}^{k\varsigma }}{\mathrm{\Gamma }(k\varsigma +1)},\\ {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})=& cos(\mathfrak{P})[1-\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}+\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}\\ & -\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{3\varsigma }}{\mathrm{\Gamma }(3\varsigma +1)}+\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{4\varsigma }}{\mathrm{\Gamma }(4\varsigma +1)}+\cdots ]\\ =& cos(\mathfrak{P})\sum _{k=0}^{\mathrm{\infty }}\frac{{(-1)}^k{(\varphi (\mathfrak{Q})-\varphi (u))}^{k\varsigma }}{\mathrm{\Gamma }(k\varsigma +1)}.\end{array}$ The exact solution for $\varphi (\mathfrak{Q})=\mathfrak{Q}$, and $u\to 1$ is given as:

72

$\begin{array}{rl}{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})& =cos(\mathfrak{P})e^{-\mathfrak{Q}},\\ {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})& =cos(\mathfrak{P})e^{-\mathfrak{Q}}.\end{array}$

General procedure of ϕ-Caputo New Iterative Method (NIM)

New iterative method

Assume the general PDE system of time fractional-order:

73

$\begin{array}{rl}& {\mathrm{\Omega }}_1(\mathfrak{P})=f(\mathfrak{P})+N_1({\mathrm{\Omega }}_1(\mathfrak{P}),\\ & {\mathrm{\Omega }}_2(\mathfrak{P})=g(\mathfrak{P})+N_2({\mathrm{\Omega }}_2(\mathfrak{P}),\end{array}$ where f and g are known functions and $N_1$ and $N_2$ are nonlinear operators from a Banach space $B\to B$. $\mathfrak{P}=({\mathfrak{P}}_1,{\mathfrak{P}}_2,\dots ,{\mathfrak{P}}_n)$ in this case. In the series form, we look for the solutions ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$ of Eq. (73):

74

$\begin{array}{rl}& {\mathrm{\Omega }}_1(\mathfrak{P})=\sum _{i=0}^{\mathrm{\infty }}{{\mathrm{\Omega }}_1}_i(\mathfrak{P}),\\ & {\mathrm{\Omega }}_2(\mathfrak{P})=\sum _{i=0}^{\mathrm{\infty }}{{\mathrm{\Omega }}_2}_i(\mathfrak{P}).\end{array}$ The nonlinear operator $N_1$ and $N_2$ can be decomposed as:

75

$\begin{array}{rl}& N_1\left(\sum _{i=0}^{\mathrm{\infty }}{{\mathrm{\Omega }}_1}_i\right)=N_1({{\mathrm{\Omega }}_1}_0)+\sum _{i=1}^{\mathrm{\infty }}\{N_1\left(\sum _{j=0}^i{{\mathrm{\Omega }}_1}_j\right)-N_1\left(\sum _{j=0}^{i-1}{{\mathrm{\Omega }}_1}_j\right)\},\\ & N_2\left(\sum _{i=0}^{\mathrm{\infty }}{{\mathrm{\Omega }}_2}_i\right)=N_2({{\mathrm{\Omega }}_2}_0)+\sum _{i=1}^{\mathrm{\infty }}\{N_2\left(\sum _{j=0}^i{{\mathrm{\Omega }}_2}_j\right)-N_2\left(\sum _{j=0}^{i-1}{{\mathrm{\Omega }}_2}_j\right)\}.\end{array}$

Putting Eq. (74) and Eq. (75) into Eq. (73) to get the subsequent result:

76

$\begin{array}{rl}& \sum _{i=0}^{\mathrm{\infty }}{{\mathrm{\Omega }}_1}_i=f+N_1({{\mathrm{\Omega }}_1}_0)+\sum _{i=1}^{\mathrm{\infty }}\{N_1\left(\sum _{j=0}^i{{\mathrm{\Omega }}_1}_j\right)-N_1\left(\sum _{j=0}^{i-1}{{\mathrm{\Omega }}_1}_j\right)\},\\ & \sum _{i=0}^{\mathrm{\infty }}{{\mathrm{\Omega }}_2}_i=g+N_2({{\mathrm{\Omega }}_2}_0)+\sum _{i=1}^{\mathrm{\infty }}\{N_2\left(\sum _{j=0}^i{{\mathrm{\Omega }}_2}_j\right)-N_2\left(\sum _{j=0}^{i-1}v{{\mathrm{\Omega }}_2}_j\right)\}.\end{array}$ Defining the recurrence relations as:

77

$\begin{array}{rl}& {{\mathrm{\Omega }}_1}_0=f,\\ & {{\mathrm{\Omega }}_1}_1=N_1({{\mathrm{\Omega }}_1}_0),\\ & {{\mathrm{\Omega }}_1}_{m+1}=N_1\left(\sum _{j=0}^m{{\mathrm{\Omega }}_1}_j\right)-N_1\left(\sum _{j=0}^{m-1}{{\mathrm{\Omega }}_1}_j\right),m=1,2,\dots \end{array}$

78

$\begin{array}{rl}& {{\mathrm{\Omega }}_2}_0=g,\\ & {{\mathrm{\Omega }}_2}_1=N_2({{\mathrm{\Omega }}_2}_0),\\ & {{\mathrm{\Omega }}_2}_{m+1}=N_2\left(\sum _{j=0}^m{{\mathrm{\Omega }}_2}_j\right)-N_2\left(\sum _{j=0}^{m-1}{{\mathrm{\Omega }}_2}_j\right),m=1,2,\dots \end{array}$ This satisfies:

79

$\begin{array}{rl}& \sum _{j=1}^{m+1}{{\mathrm{\Omega }}_1}_j=N_1\left(\sum _{j=0}^m{{\mathrm{\Omega }}_1}_j\right),m=1,2,\dots ,\\ & \sum _{j=1}^{m+1}{{\mathrm{\Omega }}_2}_j=N_2\left(\sum _{j=0}^m{{\mathrm{\Omega }}_2}_j\right),m=1,2,\dots \end{array}$ and

80

$\begin{array}{rl}& \sum _{i=0}^{\mathrm{\infty }}{{\mathrm{\Omega }}_1}_i=f+N_1\left(\sum _{i=0}^{\mathrm{\infty }}{{\mathrm{\Omega }}_1}_i\right),\\ & \sum _{i=0}^{\mathrm{\infty }}{{\mathrm{\Omega }}_2}_i=f+N_2\left(\sum _{i=0}^{\mathrm{\infty }}{{\mathrm{\Omega }}_2}_i\right).\end{array}$ The k-term approximate solutions are given by ${\mathrm{\Omega }}_1\approx {\sum }_{i=0}^{k-1}{{\mathrm{\Omega }}_1}_i$ and ${\mathrm{\Omega }}_2\approx {\sum }_{i=0}^{k-1}{{\mathrm{\Omega }}_2}_i$.

Solving PDEs using ϕ-Caputo NIM

Assume the system of PDEs of order ς:

81

$\begin{array}{rl}{D}_{\mathfrak{Q}}^{\varsigma ;\varphi }{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})& =A({\mathrm{\Omega }}_1,\partial {\mathrm{\Omega }}_1)+B(\mathfrak{P},\mathfrak{Q}),\\ D_{\mathfrak{Q}}^{\varsigma ;\varphi }{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})& =C({\mathrm{\Omega }}_2,\partial {\mathrm{\Omega }}_2)+D(\mathfrak{P},\mathfrak{Q}),m-1<\varsigma \le m,m\in \mathbb{N},\\ \frac{{\partial }^k{\mathrm{\Omega }}_1}{\partial {\mathfrak{Q}}^k}(\mathfrak{P},0)& =h_k(\mathfrak{P}),\\ \frac{{\partial }^k{\mathrm{\Omega }}_2}{\partial {\mathfrak{Q}}^k}(\mathfrak{P},0)& =g_k(\mathfrak{P}),k=0,1,\dots ,m-1,\end{array}$ where B and D are source functions and A and C are nonlinear functions of ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$ and their partial derivatives ($\partial {\mathrm{\Omega }}_1$) and ($\partial {\mathrm{\Omega }}_2$). The integral equation is equal to the initial value problem using the inverse property of the ϕ-Caputo derivative:

82

${\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})=\sum _{k=0}^{m-1}{h}_k(\mathfrak{P})\frac{{\mathfrak{Q}}^k}{k!}+{\mathcal{J}}_t^{\varsigma ;\varphi }B(\mathfrak{P},\mathfrak{Q})+{\mathcal{J}}_t^{\varsigma ;\varphi }A({\mathrm{\Omega }}_1,\partial {\mathrm{\Omega }}_1)\equiv f+N_1({\mathrm{\Omega }}_1),$

83

${\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})=\sum _{k=0}^{m-1}{g}_k(\mathfrak{P})\frac{{\mathfrak{Q}}^k}{k!}+{\mathcal{J}}_{\mathfrak{Q}}^{\varsigma ;\varphi }D(\mathfrak{P},\mathfrak{Q})+{\mathcal{J}}_{\mathfrak{Q}}^{\varsigma ;\varphi }C({\mathrm{\Omega }}_2,\partial {\mathrm{\Omega }}_2)\equiv w+N_2({\mathrm{\Omega }}_2),$ where $\begin{array}{rl}f& :=\sum _{k=0}^{m-1}{h}_k(\mathfrak{P})\frac{{\mathfrak{Q}}^k}{k!}+{\mathcal{J}}_{\mathfrak{Q}}^{\varsigma ;\varphi }B(\mathfrak{P},\mathfrak{Q}),\\ N_1({\mathrm{\Omega }}_1)& :={\mathcal{J}}_{\mathfrak{Q}}^{\varsigma ;\varphi }A({\mathrm{\Omega }}_1,\partial {\mathrm{\Omega }}_1).\end{array}$ and $\begin{array}{rl}w& :=\sum _{k=0}^{m-1}{g}_k(\mathfrak{P})\frac{{\mathfrak{Q}}^k}{k!}+{\mathcal{J}}_{\mathfrak{Q}}^{\varsigma ;\varphi }D(\mathfrak{P},\mathfrak{Q}),\\ N_2({\mathrm{\Omega }}_2)& :={\mathcal{J}}_{\mathfrak{Q}}^{\varsigma ;\varphi }C({\mathrm{\Omega }}_2,\partial {\mathrm{\Omega }}_2).\end{array}$ Applying the recurrence relations Eq. (77) and Eq. (78) yields the solution to Eq. (82).

Solution of Burger equation through ϕ-Caputo NIM

84

$\begin{array}{rl}& {\mathcal{D}}_{\mathfrak{Q}}^{\varsigma ;\varphi }\mathrm{\Omega }=\frac{{\partial }^2\mathrm{\Omega }}{\partial {\mathfrak{P}}^2}-\mathrm{\Omega }\frac{\partial \mathrm{\Omega }}{\partial \mathfrak{P}},\\ & \text{ where }0<\varsigma \le 1,\end{array}$ with the initial conditions listed below:

85

$\begin{array}{r}\mathrm{\Omega }(\mathfrak{P},0)=\frac{1}{2}-\frac{1}{2}tanh\left(\frac{\mathfrak{P}}{4}\right),\end{array}$ The integral from of Eq. (84) and Eq. (85) are given as:

86

$\begin{array}{rl}& \mathrm{\Omega }(\mathfrak{P},\mathfrak{Q})=\frac{1}{2}-\frac{1}{2}tanh\left(\frac{\mathfrak{P}}{4}\right)+{\mathcal{J}}_t^{\varsigma ;\varphi }(\frac{{\partial }^2\mathrm{\Omega }}{\partial {\mathfrak{P}}^2}-\mathrm{\Omega }\frac{\partial \mathrm{\Omega }}{\partial \mathfrak{P}})\\ & =f+N(\mathrm{\Omega }),\end{array}$ where $f=w=\frac{1}{2}-\frac{1}{2}tanh\left(\frac{\mathfrak{P}}{4}\right)$ and $N(\mathrm{\Omega })={\mathcal{J}}_t^{\varsigma ;\varphi }(\frac{{\partial }^2\mathrm{\Omega }}{\partial {\mathfrak{P}}^2}-\mathrm{\Omega }\frac{\partial \mathrm{\Omega }}{\partial \mathfrak{P}})$. Using Eq. (77), we obtain:

87

$\begin{array}{rl}& {\mathrm{\Omega }}_0(\mathfrak{P},\mathfrak{Q})=\frac{1}{2}-\frac{1}{2}tanh\left(\frac{\mathfrak{P}}{4}\right),\\ & {{\mathrm{\Omega }}_1}_1(\mathfrak{P},\mathfrak{Q})=\frac{1}{8\mathrm{\Gamma }(\varsigma +1)(cosh\left(\frac{\mathfrak{P}}{2}\right)+1)}\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}.\end{array}$ The four-terms approximate solutions are given as:

88

$\begin{array}{rl}& \mathrm{\Omega }(\mathfrak{P},\mathfrak{Q})\\ & =\frac{1}{2}-\frac{1}{2}tanh\left(\frac{\mathfrak{P}}{4}\right)+\frac{1}{8\mathrm{\Gamma }(\varsigma +1)(cosh\left(\frac{\mathfrak{P}}{2}\right)+1)}\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}\\ & +\frac{1}{512}{(\varphi (\mathfrak{Q})-\varphi (u))}^{2\varsigma }tanh\left(\frac{\mathfrak{P}}{4}\right){\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right)\\ & \times (\frac{4^{\varsigma }{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }\mathrm{\Gamma }(\varsigma +\frac{1}{2}){\text{sech}}^2\left(\frac{\mathfrak{P}}{4}\right)}{\sqrt{\pi }\mathrm{\Gamma }(\varsigma +1)\mathrm{\Gamma }(3\varsigma +1)}+\frac{8}{\mathrm{\Gamma }(2\varsigma +1)}).\end{array}$

Solution of system of Burger equation through ϕ-Caputo NIM

89

$\begin{array}{rl}{}_0{}^{\mathcal{C}}\mathcal{D}_t^{\varsigma ;\varphi }{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})=& \frac{{\partial }^2{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})}{\partial {\mathfrak{P}}^2}+2{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}}\\ & -{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}}-{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}},\end{array}$

90

$\begin{array}{rl}{}_0{}^{\mathcal{C}}\mathcal{D}_t^{\varsigma ;\varphi }{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})=& \frac{{\partial }^2{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})}{\partial {\mathfrak{P}}^2}+2{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}}\\ & -{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}}-{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}},\\ & \text{ where }0<\varsigma \le 1.\end{array}$ Initial conditions:

91

$\begin{array}{r}{\mathrm{\Omega }}_1(\mathfrak{P},0)=cos(\mathfrak{P}),\end{array}$

92

$\begin{array}{r}{\mathrm{\Omega }}_2(\mathfrak{P},0)=cos(\mathfrak{P}).\end{array}$ The integral form of Eq. (89) and Eq. (90) are given as;

93

$\begin{array}{rl}{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})=& cos(\mathfrak{P})+{\mathcal{J}}_t^{\varsigma ;\varphi }(\frac{{\partial }^2{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})}{\partial {\mathfrak{P}}^2}+2{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}}\\ & -{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}}-{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}})\\ =& f+N_1({\mathrm{\Omega }}_1),\end{array}$

94

$\begin{array}{rl}{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})=& cos(\mathfrak{P})+{\mathcal{J}}_t^{\varsigma ;\varphi }(\frac{{\partial }^2{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})}{\partial {\mathfrak{P}}^2}+2{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}}\\ & -{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}}-{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}})\\ =& w+N_2({\mathrm{\Omega }}_2),\end{array}$ where $f=w=cos(\mathfrak{P})$ and $N_1({\mathrm{\Omega }}_1)={\mathcal{J}}_t^{\varsigma ;\varphi }(\frac{{\partial }^2{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})}{\partial {\mathfrak{P}}^2}+2{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}}-{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}}-{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}})$ and $N_2({\mathrm{\Omega }}_2)={\mathcal{J}}_t^{\varsigma ;\varphi }(\frac{{\partial }^2{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})}{\partial {\mathfrak{P}}^2}+2{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})\times \frac{\partial {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}}-{\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}}-{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})\frac{\partial {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})}{\partial \mathfrak{P}})$. Using Eq. (77) and Eq. (78), we obtain:

95

$\begin{array}{rl}& {{\mathrm{\Omega }}_1}_0(\mathfrak{P},\mathfrak{Q})=cos(\mathfrak{P}),{{\mathrm{\Omega }}_1}_1(\mathfrak{P},\mathfrak{Q})=-cos(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)},\end{array}$

96

$\begin{array}{rl}& {{\mathrm{\Omega }}_2}_0(\mathfrak{P},\mathfrak{Q})=cos(\mathfrak{P}),{{\mathrm{\Omega }}_2}_1(\mathfrak{P},\mathfrak{Q})=-cos(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}.\end{array}$ The four-terms approximate solutions are given as:

97

$\begin{array}{rl}{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})=& cos(\mathfrak{P})-cos(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}+cos(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}\\ & -cos(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{3\varsigma }}{\mathrm{\Gamma }(3\varsigma +1)},\\ {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})=& cos(\mathfrak{P})-cos(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}+cos(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}\\ & -cos(\mathfrak{P})\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{3\varsigma }}{\mathrm{\Gamma }(3\varsigma +1)}.\end{array}$ The final solution through $\varphi -Caputo$ NIM algorithm can also be written as: $\begin{array}{rl}{\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})=& cos(\mathfrak{P})[1-\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}+\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}\\ & -\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{3\varsigma }}{\mathrm{\Gamma }(3\varsigma +1)}+\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{4\varsigma }}{\mathrm{\Gamma }(4\varsigma +1)}+\cdots ]\\ =& cos(\mathfrak{P})\sum _{k=0}^{\mathrm{\infty }}\frac{{(-1)}^k{(\varphi (\mathfrak{Q})-\varphi (u))}^{k\varsigma }}{\mathrm{\Gamma }(k\varsigma +1)},\\ {\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})& =cos(\mathfrak{P})[1-\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{\varsigma }}{\mathrm{\Gamma }(\varsigma +1)}+\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{2\varsigma }}{\mathrm{\Gamma }(2\varsigma +1)}\\ & -\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{3\varsigma }}{\mathrm{\Gamma }(3\varsigma +1)}+\frac{{(\varphi (\mathfrak{Q})-\varphi (u))}^{4\varsigma }}{\mathrm{\Gamma }(4\varsigma +1)}+\cdots ]\\ =& cos(\mathfrak{P})\sum _{k=0}^{\mathrm{\infty }}\frac{{(-1)}^k{(\varphi (\mathfrak{Q})-\varphi (u))}^{k\varsigma }}{\mathrm{\Gamma }(k\varsigma +1)}.\end{array}$

Numerical result and discussion

This section deals with numerical and graphical discussion of the coupled fractional Burgers equations by the use of two semi-analytical methods: The ϕ-Caputo RPS and the ϕ-Caputo NIM. The main purpose is to determine the validity, convergence, and properties of the solutions for different fractional orders $0<\varsigma \le 1$. To visualise the dynamic response of the solution $\mathrm{\Omega }(\mathfrak{P},\mathfrak{Q})$ for different fractional orders via both RPSM and NIM for the first problem, the Figs. 1 and 2 have been plotted at time $\mathfrak{Q}=0.1$. The 3D surface plot of Fig. 1 illustrates the solution in terms of fractional order, which exposes the continuous and smooth profiles for variable fractional order, and it reflects the nonlocal memory effect of fractional calculus. In Fig. 2, subfigure (a) shows 2D profiles for RPSM and subfigure (b) the results by making use of NIM, which are favourably in agreement with the expected analysis trend. The comparison Table 1 also supports the accuracy of both methods by recording the numerical results with the exact solution for $p=1$. The absolute errors make it clear that both RPSM and NIM tend to result in a very accurate approximation, although RPSM has a slight tendency of yielding smaller error magnitudes in the initial condition but still rather good approximation is achieved by NIMs at the beginning and expansion of the domain. These results demonstrate the credibility and efficiency of the methods of the present work for solving the problems of fractional-order involving the ϕ-Caputo derivative. For the second problem, 3D surface plots of ${\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})$ and ${\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})$ as functions of $\mathfrak{Q}$ and $\mathfrak{P}$ are shown in Figs. 3 and 5 for $\mathfrak{Q}=0.1$ and for different values of ς. These surfaces show the evolution of the solutions in space $\mathfrak{P}$ as the fractional parameter varies. Importantly, when ς is decreased from 1.0 to lower values, the diffusive behavior dominates. This is a direct consequence of the memory factor present in fractional derivatives, smaller ς will induce faster decrease of the evolutionary amplitude, that is, the former one has a stronger memory factor than the latter one. Figure 4 shows the effect of the fractional-order parameter ς on the solution ${\mathrm{\Omega }}_1$. It was clearly observed that the oscillatory behavior is very different as compared to the classical case $\varsigma =1.0$. Reduced fractional orders $\varsigma =0.4,0.6$, and 0.8 change the amplitude and phase of oscillations, with damping-like effects and with the capability of fractional models to reveal memory and hereditary behavior that is not found in integer-order models. Similarly, Fig. 6 demonstrates the dynamics of the ${\mathrm{\Omega }}_2$ with fractional orders providing smoother curves with lower amplitude and slower decay than the classical case. The solutions are converged to the integer-order model as the value of the parameter ς approach to 1, and the framework is therefore shown to be consistent. In general, Figs. 4 and 6 show that the classical models are generalized by the ϕ-Caputo fractional-order formulations to nonlocality and memory as well, and therefore provide nonlinear transport and wave propagation problems with more rich dynamics. Tables 2, 3 show numerical computation with the exact solution, for $\varsigma =1.0$, and $\mathfrak{Q}=0.1$. For both of the tables, we observe that the absolute errors for ${\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})$ and ${\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})$ are very small, indicating a very good fit and a high numerical precision. This is an indication to the fact that both RPS, and NIM could strongly approximate the true solution of the fractional Burgers system.

[See PDF for image]

Figure 1

Different fractional-order surface analysis of the solution $\mathrm{\Omega }(\mathfrak{P},\mathfrak{Q})$ for $\mathfrak{Q}=0.1$ using the proposed methods, and its comparison with the classical integer-order case ($\varsigma =1$)

[See PDF for image]

Figure 2

Different fractional-order 2D analysis of the solution $\mathrm{\Omega }(\mathfrak{P},\mathfrak{Q})$ for $\mathfrak{Q}=0.1$, (a) represents the fractional-order comparison of RPSM, and (b) represents the fractional-order comparison of NIM with the classical integer-order case ($\varsigma =1$)

[See PDF for image]

Figure 3

Different fractional-order surface analysis of the solution ${\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})$ for $\mathfrak{Q}=0.1$ using the proposed methods, and its comparison with the classical integer-order case ($\varsigma =1$)

[See PDF for image]

Figure 4

Different fractional-order 2D analysis of the solution ${\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})$ for $\mathfrak{Q}=0.1$ using the proposed methods, and its comparison with the classical integer-order case ($\varsigma =1$)

[See PDF for image]

Figure 5

Different fractional-order surface analysis of the solution ${\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})$ for $\mathfrak{Q}=0.1$ using the proposed methods, and its comparison with the classical integer-order case ($\varsigma =1$)

[See PDF for image]

Figure 6

Different fractional-order 2D analysis of the solution ${\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})$ for $\mathfrak{Q}=0.1$ using the proposed methods, and its comparison with the classical integer-order case ($\varsigma =1$)

Table 1. Comparison of $\varphi -Caputo$ RPSM and NIM solutions with Exact Values and Absolute Errors for $\varsigma =1$

$\mathfrak{P}$	RPSM Solution	NIM Solution	Exact	Absolute Error (ς = 1)
				RPSM	NIM
0	0.5062500	0.5062500	0.5062500	3.2550 × 10−7	3.2550 × 10−7
1	0.3834340	0.3834340	0.3834330	2.5262 × 10−7	3.9352 × 10−7
2	0.2738850	0.2738850	0.2738850	9.4008 × 10−8	2.8009 × 10−7
3	0.1861840	0.1861840	0.1861840	3.9602 × 10−8	1.0757 × 10−7
4	0.1218530	0.1218530	0.1218530	1.0083 × 10−7	1.3377 × 10−8
5	0.0776294	0.0776294	0.0776295	1.0592 × 10−7	6.2497 × 10−8
6	0.0485681	0.0485681	0.0485682	8.6062 × 10−8	6.6819 × 10−8
7	0.0300319	0.0300319	0.0300320	6.1734 × 10−8	5.3796 × 10−8
8	0.0184331	0.0184331	0.0184331	4.1341 × 10−8	3.8208 × 10−8
9	0.0112619	0.0112619	0.0112619	2.6604 × 10−8	2.5401 × 10−8
10	0.0068611	0.0068611	0.0068611	1.6721 × 10−8	1.6267 × 10−8

Table 2. Comparison of Approximate solution through $\varphi -Caputo$ RPSM, $\varphi -Caputo$ NIM and Exact Solutions with Absolute Errors for $\mathfrak{Q}=0.1$

$\mathfrak{P}$	${\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})$	${\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})$	Exact	Absolute Error (ς = 1.0)
				${\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})$	${\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})$
0.0	0.990050	0.990050	0.990050	4.1583481102946 × 10−10	4.1583481102946 × 10−10
0.1	0.985104	0.985104	0.985104	4.1375736170579 × 10−10	4.1375736170579 × 10−10
0.2	0.970315	0.970315	0.970315	4.0754588592761 × 10−10	4.0754588592761 × 10−10
0.3	0.945831	0.945831	0.945831	3.9726211209512 × 10−10	3.9726211209512 × 10−10
0.4	0.911896	0.911896	0.911896	3.8300917992728 × 10−10	3.8300917992728 × 10−10
0.5	0.868850	0.868850	0.868850	3.6492930899356 × 10−10	3.6492930899356 × 10−10
0.6	0.817123	0.817123	0.817123	3.4320324360237 × 10−10	3.4320324360237 × 10−10
0.7	0.757232	0.757232	0.757232	3.1804792133272 × 10−10	3.1804792133272 × 10−10
0.8	0.689774	0.689774	0.689774	2.8971480769968 × 10−10	2.8971480769968 × 10−10
0.9	0.615425	0.615425	0.615425	2.5848700957453 × 10−10	2.5848700957453 × 10−10
1.0	0.534926	0.534926	0.534926	2.2467638860490 × 10−10	2.2467638860490 × 10−10

Table 3. Comparison of Approximate solution through $\varphi -Caputo$ RPSM, $\varphi -Caputo$ NIM and Exact Solutions with Absolute Errors for $\mathfrak{Q}=0.1$

$\mathfrak{P}$	${\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})$	${\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})$	Exact	Absolute Error (ς = 1.0)
				${\mathrm{\Omega }}_1(\mathfrak{P},\mathfrak{Q})$	${\mathrm{\Omega }}_2(\mathfrak{P},\mathfrak{Q})$
0.0	0.904833	0.904833	0.904837	4.0847026261392 × 10−6	4.0847026261392 × 10−6
0.1	0.900313	0.900313	0.900317	4.0642961269821 × 10−6	4.0642961269821 × 10−6
0.2	0.886797	0.886797	0.886801	4.0032805243006 × 10−6	4.0032805243006 × 10−6
0.3	0.864420	0.864420	0.864424	3.9022654660936 × 10−6	3.9022654660936 × 10−6
0.4	0.833407	0.833407	0.833410	3.7622602609976 × 10−6	3.7622602609976 × 10−6
0.5	0.794066	0.794066	0.794070	3.5846637952419 × 10−6	3.5846637952419 × 10−6
0.6	0.746791	0.746791	0.746795	3.3712505536076 × 10−6	3.3712505536076 × 10−6
0.7	0.692055	0.692055	0.692058	3.1241528910541 × 10−6	3.1241528910541 × 10−6
0.8	0.630403	0.630403	0.630406	2.8458397254294 × 10−6	2.8458397254294 × 10−6
0.9	0.562453	0.562453	0.562456	2.5390918698686 × 10−6	2.5390918698686 × 10−6
1.0	0.488884	0.488884	0.488886	2.2069742477426 × 10−6	2.2069742477426 × 10−6

The minimal error stability in spatial position $\mathfrak{P}$ reflects the perfectly uniform convergence of the solution techniques. In addition, both of the solutions from both methods are consistently dependent on each other proving the good consequence of the algorithms. The obtain results are evidence that both the ϕ-Caputo RPS and ϕ-Caputo NIM are highly accurate and whose solutions are quantitatively consistent. The fractional order ς is an important parameter that fundamentally influences the diffusivity and wave attenuation. For $\varsigma =1$, the model converges to its classical limits but for lower ς, the system is subdiffusive implying that it represents systems that exhibit long concepts in the memory or transience in transport phenomena. The graphical and tabulated results are proof that the fractional models can be executed with the presented methods.

Conclusions

In this work, we have explored the analytical and semi-analytical strategies for solving generalized fractional-order Burgers systems, emphasizing their critical role in modeling complex physical processes characterized by memory and hereditary effects. By integrating generalized fractional operators such as the Caputo-Hadamard and ϕ-Caputo derivatives into modern solution frameworks like the RPSM and NIM, we demonstrate enhanced accuracy, stability, and convergence in addressing nonlinear and ill-posed fractional models. Our findings underscore the significance of generalized fractional calculus not merely as a theoretical advancement but as a practical tool for capturing temporal and spatial nonlocalities inherent in real-world phenomena. The synergy between these advanced operators and iterative algorithms offers promising new directions for mathematical modeling in physics, engineering, and applied sciences. This study provides a solid foundation for future investigations, paving the way for more robust numerical and analytical approaches to tackle a broader class of fractional differential equations. Further research may focus on extending these techniques to multi-dimensional systems, exploring hybrid methods, and validating theoretical models against experimental data from physical systems exhibiting anomalous diffusion and nonlinear dynamics.

Acknowledgements

This research has been funded by Scientific Research Deanship at University of Ha’il-Saudi Arabia through project number < RG-24-050>.

Author contributions

Conceptualization, Z.A. R.S and T.S.A.; methodology, R.S.; software, S.A.; validation, M.A.; formal analysis, M.A.A.; investigation, Z.A.; resources, T.S.A.; data curation, S.J.; writing-original draft preparation, M.A.; writing-review and editing, M.A.A.; visualization, S.J.; supervision, S.A.; project administration, T.S.A. All authors have read and agreed to the published version of the manuscript.

Funding information

There was no financial support funding for this study.

Data availability

No data sets were generated or analysed during the current study.

Declarations

Consent for publication

“Not applicable” in this section.

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.

References

1. Podlubny, I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications; 1998; Amsterdam, Elsevier: 198

2. Kilbas, A.A.; Srivastava, H.M.; Trujillo, J.J. Theory and Applications of Fractional Differential Equations; 2006; Amsterdam, Elsevier: [DOI: https://dx.doi.org/10.1016/S0304-0208(06)80001-0] 204

3. Burgers, J.M. Solutions of the linear diffusion equation with a boundary condition referring to a parabola. The Nonlinear Diffusion Equation: Asymptotic Solutions and Statistical Problems; 1974; pp. 46-71. [DOI: https://dx.doi.org/10.1007/978-94-010-1745-9_5]

4. Kumar, N.K. A review on Burgers’ equations and it’s applications. J. Inst. Sci. Technol.; 2023; 28, 2 pp. 49-52. [DOI: https://dx.doi.org/10.3126/jist.v28i2.61073]

5. Boritchev, A.: Turbulence for the generalised Burgers equation (2013). arXiv preprint. arXiv:1304.6814

6. Whitham, G.B. Linear and Nonlinear Waves; 2011; New York, Wiley:

7. Wazwaz, A.M., Alhejaili, W., El-Tantawy, S.A.: Physical multiple shock solutions to the integrability of linear structures of Burgers hierarchy. Phys. Fluids 35(12) (2023)

8. Hussain, M.; Haq, S.; Ghafoor, A.; Ali, I. Numerical solutions of time-fractional coupled viscous Burgers’ equations using meshfree spectral method. Comput. Appl. Math.; 2020; 39, 14036538 [DOI: https://dx.doi.org/10.1007/s40314-019-0985-3] 6.

9. Jarad, F.; Abdeljawad, T.; Baleanu, D. Caputo-type modification of the Hadamard fractional derivatives. Adv. Difference Equ.; 2012; 2012, 12992066 [DOI: https://dx.doi.org/10.1186/1687-1847-2012-142] 142.

10. Green, C.W.H.; Liu, Y.; Yan, Y. Numerical methods for Caputo-Hadamard fractional differential equations with graded and non-uniform meshes. Mathematics; 2021; 9, 21 [DOI: https://dx.doi.org/10.3390/math9212728] 2728.

11. Baleanu, D.; Diethelm, K.; Scalas, E.; Trujillo, J.J. Fractional Calculus: Models and Numerical Methods; 2012; Singapore, World Scientific: 3

12. Atangana, A., Baleanu, D.: New fractional derivatives with nonlocal and non-singular kernel: theory and application to heat transfer model (2016). arXiv preprint. arXiv:1602.03408

13. Abdeljawad, T. On conformable fractional calculus. J. Comput. Appl. Math.; 2015; 279, pp. 57-66.3293309 [DOI: https://dx.doi.org/10.1016/j.cam.2014.10.016]

14. Sousa, J.V.D.C.; De Oliveira, E.C. On the ?-Hilfer fractional derivative. Commun. Nonlinear Sci. Numer. Simul.; 2018; 60, pp. 72-91.3762997 [DOI: https://dx.doi.org/10.1016/j.cnsns.2018.01.005]

15. Sousa, J.V.D.C.; De Oliveira, E.C. Leibniz type rule: ?-Hilfer fractional operator. Commun. Nonlinear Sci. Numer. Simul.; 2019; 77, pp. 305-311.3948892 [DOI: https://dx.doi.org/10.1016/j.cnsns.2019.05.003]

16. Cen, J., Sousa, J.V.D.C., Costa, F.S.: Existence and nonexistence of positive solutions for fractional p-logistic equation. Complex Var. Elliptic Equ., 1–11 (2025)

17. Wahash, H.A.; Abdo, M.S.; Saeed, A.M.; Panchal, S.K. Singular fractional differential equations with ψ-Caputo operator and modified Picard’s iterative method. Appl. Math. E-Notes; 2020; 20, pp. 215-229.4117066

18. Patnaik, S.; Semperlotti, F. Application of variable-and distributed-order fractional operators to the dynamic analysis of nonlinear oscillators. Nonlinear Dynam.; 2020; 100, 1 pp. 561-580. [DOI: https://dx.doi.org/10.1007/s11071-020-05488-8]

19. Derbazi, C., Al-Mdallal, Q.M., Jarad, F., Baitiche, Z.: Some qualitative properties of solutions for nonlinear fractional differential equation involving two Φ–Caputo fractional derivatives (2021). arXiv preprint. arXiv:2108.13758

20. Hattaf, K. A new class of generalized fractal and fractal-fractional derivatives with non-singular kernels. Fractal Fract.; 2023; 7, 5 [DOI: https://dx.doi.org/10.3390/fractalfract7050395] 395.

21. Inc, M.; Korpinar, Z.S.; Al Qurashi, M.M.; Baleanu, D. A new method for approximate solutions of some nonlinear equations: residual power series method. Adv. Mech. Eng.; 2016; 8, 4 [DOI: https://dx.doi.org/10.1177/1687814016644580] 1687814016644580.

22. Liaqat, M.I.; Etemad, S.; Rezapour, S. A novel analytical Aboodh residual power series method for solving linear and nonlinear time-fractional partial differential equations with variable coefficients. AIMS Math.; 2022; 7, 9 pp. 16917-16948.4468454 [DOI: https://dx.doi.org/10.3934/math.2022929]

23. Mirzaei, H.R.; Bayat, F.; Miralikhani, K. A semi-analytic approach for determining Marx generator optimum setup during power transformers factory test. IEEE Trans. Power Deliv.; 2020; 36, 1 pp. 10-18. [DOI: https://dx.doi.org/10.1109/TPWRD.2020.2982270]

24. Saratha, S.R.; Yildirim, A. A novel approach to fractional differential equations via JSN transform coupled fractional residual power series method. Phys. Scr.; 2024; 99, 12 [DOI: https://dx.doi.org/10.1088/1402-4896/ad92bb] 125285.

25. Biswas, P.; Andallah, L.S.; Eusha-Bin-Hafiz, K.M. Estimation of thermal pollution using numerical simulation of energy equation coupled with viscous Burgers’ equation. Amer. J. Comput. Math.; 2022; 12, 3 pp. 306-313. [DOI: https://dx.doi.org/10.4236/ajcm.2022.123020]

26. He, J.H. Variational iteration method a kind of non-linear analytical technique: some examples. Int. J. Non-Linear Mech.; 1999; 34, 4 pp. 699-708. [DOI: https://dx.doi.org/10.1016/S0020-7462(98)00048-1]

27. Noor, M.A.; Waseem, M.; Noor, K.I.; Ali, M.A. New iterative technique for solving nonlinear equations. Appl. Math. Comput.; 2015; 265, pp. 1115-1125.3373551

28. Rashidi, M.; Tashakori, S.; Kalhori, H.; Bahmanpour, M.; Li, B. Iterative-based impact force identification on a bridge concrete deck. Sensors; 2023; 23, 22 [DOI: https://dx.doi.org/10.3390/s23229257] 9257.

29. Kumar, N.; Majumdar, R.; Singh, S. Physics-based preconditioning of Jacobian free Newton Krylov for Burgers’ equation using modified nodal integral method. Prog. Nucl. Energy; 2019; 117, [DOI: https://dx.doi.org/10.1016/j.pnucene.2019.103104] 103104.

30. Odibat, Z.; Momani, S. Numerical methods for nonlinear partial differential equations of fractional order. Appl. Math. Model.; 2008; 32, 1 pp. 28-39. [DOI: https://dx.doi.org/10.1016/j.apm.2006.10.025]

31. Gupta, M.; Prakash Shukla, J.; Singh, B.K. Study of dynamical behavior of Burgers’ equation in quantum modeling. Math. Methods Appl. Sci.; 2025; 48, 7 pp. 7890-7905.4889498 [DOI: https://dx.doi.org/10.1002/mma.9378]

32. Tarasov, V.E. Fractional Dynamics: Applications of Fractional Calculus to Dynamics of Particles, Fields and Media; 2011; Berlin, Springer:

33. Magin, R.L. Fractional calculus in bioengineering: a tool to model complex dynamics. Proceedings of the 13th International Carpathian Control Conference (ICCC); 2012; New York, IEEE Press: pp. 464-469. [DOI: https://dx.doi.org/10.1109/CarpathianCC.2012.6228688]

34. Diethelm, K. An efficient parallel algorithm for the numerical solution of fractional differential equations. Fract. Calc. Appl. Anal.; 2011; 14, 3 pp. 475-490.2837642 [DOI: https://dx.doi.org/10.2478/s13540-011-0029-1]

35. Odibat, Z. A finite difference-based Adams-type approach for numerical solution of nonlinear fractional differential equations: a fractional Lotka-Volterra model as a case study. J. Comput. Nonlinear Dyn.; 2025; 20, 14794477 [DOI: https://dx.doi.org/10.1115/1.4066885] 011004.

36. Abdeljawad, T.: Fractional operators with generalized Mittag-Leffler kernels and their iterated differintegrals. Chaos 29(2) (2019)

37. Bas, E.; Ozarslan, R. Real world applications of fractional models by Atangana-Baleanu fractional derivative. Chaos Solitons Fractals; 2018; 116, pp. 121-125.3887016 [DOI: https://dx.doi.org/10.1016/j.chaos.2018.09.019]

38. Islam, M.: Reconstruct the Two Temperature Generalized Thermo-elasticity with Memory Dependent Derivative. Calcutta Institute of Theoretical Physics, p. 55 (2023)

39. Mainardi, F. Fractional calculus: some basic problems in continuum and statistical mechanics. Fractals and Fractional Calculus in Continuum Mechanics; 1997; pp. 291-348. [DOI: https://dx.doi.org/10.1007/978-3-7091-2664-6_7]

40. Moghaddam, B.P.; Machado, J.T.; Babaei, A. A computationally efficient method for tempered fractional differential equations with application. Comput. Appl. Math.; 2018; 37, 3 pp. 3657-3671.3826050 [DOI: https://dx.doi.org/10.1007/s40314-017-0522-1]

41. Khirsariya, S.; Rao, S.; Chauhan, J. Semi-analytic solution of time-fractional Korteweg-de Vries equation using fractional residual power series method. Results Nonlinear Anal.; 2022; 5, 3 pp. 222-234. [DOI: https://dx.doi.org/10.53006/rna.1024308]

42. Dehghan, M.; Manafian, J.; Saadatmandi, A. Solving nonlinear fractional partial differential equations using the homotopy analysis method. Numer. Methods Partial Differential Equations; 2010; 26, 2 pp. 448-479.2605472 [DOI: https://dx.doi.org/10.1002/num.20460]

43. Odibat, Z. On Legendre polynomial approximation with the VIM or HAM for numerical treatment of nonlinear fractional differential equations. J. Comput. Appl. Math.; 2011; 235, 9 pp. 2956-2968.2771279 [DOI: https://dx.doi.org/10.1016/j.cam.2010.12.013]

44. Darassi, M.H., Yasin, O., Ahmed, M.:. A semi-analytical method to solve the Fitzhugh-Nagumo equation

45. Mao, Z.; Shen, J. Hermite spectral methods for fractional PDEs in unbounded domains. SIAM J. Sci. Comput.; 2017; 39, 5 pp. A1928-A1950.3696054 [DOI: https://dx.doi.org/10.1137/16M1097109]

46. Sigl, J. Nonlinear residual minimization by iteratively reweighted least squares. Comput. Optim. Appl.; 2016; 64, 3 pp. 755-792.3506232 [DOI: https://dx.doi.org/10.1007/s10589-016-9829-x]

47. Khan, N.A.; Ahmad, S. Framework for treating non-linear multi-term fractional differential equations with reasonable spectrum of two-point boundary conditions. AIMS Math.; 2019; 4, 4 pp. 1181-1202.4137028 [DOI: https://dx.doi.org/10.3934/math.2019.4.1181]

48. Sweilam, N.H.; Khader, M.M.; Al-Bar, R.F. Numerical studies for a multi-order fractional differential equation. Phys. Lett. A; 2007; 371, 1–2 pp. 26-33.2419274 [DOI: https://dx.doi.org/10.1016/j.physleta.2007.06.016]

49. Almoneef, A.A.; Barakat, M.A.; Hyder, A.A. Analysis of the fractional HIV model under proportional Hadamard-Caputo operators. Fractal Fract.; 2023; 7, 3 [DOI: https://dx.doi.org/10.3390/fractalfract7030220] 220.

50. Alqahtani, A.M.; Mihoubi, H.; Arioua, Y.; Bouderah, B. Analytical solutions of time-fractional Navier-Stokes equations employing homotopy perturbation-Laplace transform method. Fractal Fract.; 2024; 9, 1 [DOI: https://dx.doi.org/10.3390/fractalfract9010023] 23.

51. Almeida, R. A Caputo fractional derivative of a function with respect to another function. Commun. Nonlinear Sci. Numer. Simul.; 2017; 44, pp. 460-481.3554830 [DOI: https://dx.doi.org/10.1016/j.cnsns.2016.09.006]

52. Samko, S.G. Fractional Integrals and Derivatives. Theory and Applications; 1993;

53. Gambo, Y.Y.; Jarad, F.; Baleanu, D.; Abdeljawad, T. On Caputo modification of the Hadamard fractional derivatives. Adv. Difference Equ.; 2014; 2014, 13213915 [DOI: https://dx.doi.org/10.1186/1687-1847-2014-10] 10.

54. Company, R.; Jodar, L.; Pintos, J.R. A numerical method for European option pricing with transaction costs nonlinear equation. Math. Comput. Modelling; 2009; 50, 5–6 pp. 910-920.2569252 [DOI: https://dx.doi.org/10.1016/j.mcm.2009.05.019]

55. Gazizov, R.K.; Ibragimov, N.H. Lie symmetry analysis of differential equations in finance. Nonlinear Dynam.; 1998; 17, 4 pp. 387-407.1671842 [DOI: https://dx.doi.org/10.1023/A:1008304132308]