An efficient computational scheme for nonlinear time fractional systems of partial differential equations arising in physical sciences

In this paper, we broaden the utilization of a beautiful computational scheme, residual power series method (RPSM), to attain the fractional power series solutions of nonhomogeneous and homogeneous nonlinear time-fractional systems of partial differential equations. This paper considers the fractional derivatives of Caputo-type. The approximate solutions of given systems of equations are calculated through the utilization of the provided initial conditions. This iterative scheme generates the fast convergent series solutions with conveniently determinable components. The implementation of this numerical scheme clearly exhibits its effectiveness, reliability and easiness regarding the procedure of the solution, as well as its better approximation. The repercussions of the fractional order of Caputo derivatives on solutions are depicted through graphical presentations for various particular cases.


Introduction
Physical phenomena of the miscellaneous fields of engineering and science can be modelled very appropriately by utilizing fractional partial differential equations (FPDEs). Presently the theory of fractional calculus is equipped with fantastic tools to describe the dynamical behaviour, and memory related characteristics of scientific systems and processes. Various authors have used fractional differential equations (FDEs) in the modelling and analysis of scientific phenomena in different fields of knowledge [1][2][3]. Presently the theory of fractional calculus has been widely utilized in various fields and it is growing very fast in developing models due to its relation with memory and fractals which are abundant in real physical systems. Fractional order modelling minimizes the inaccuracy that arises from the ignorance of significant real parameters. It permits a greater degree of freedom in the model compared to an integer-order system. FDEs are equipped with magnificent techniques for the characterization of hereditary and memory characteristics which are simply ignored by the integer order system. In addition, they are also suitable in modelling the behaviour of real systems and also relevant in the investigation of dynamical systems. FDEs are also appropriate in case of modelling of systems with longer-range interactivity both in time and space. Fractional order systems are normally associated to the systems of memory which allows the incorporation of several types of information. The stability region increases in case of a fractional order system as compared to its integer order framework. Fractional calculus also provides non-local fractional derivative operators and numerical results with high accuracy. In addition, fractional order systems ultimately converge to the integer-order systems.
The non-local property of the fractional operator is the most advantageous feature in this context. The theory of fractional calculus develops many generalizations with respect to non-local characteristics of fractional operators, enhanced degree of freedom, maximum utilization of information, and these characteristics only occur in the case of fractional order systems and not of integer-order systems. The numerical schemes provided by fractional calculus originate the deeper understanding of complex systems and reduce the computational work regarding the solution procedure. Accurate analytical solutions are not found easily in the case of FDEs. Thus in the last two decades, many iterative schemes such as Adomian decomposition method (ADM) [4][5][6], Variational iteration method (VIM) [7][8][9][10], Homotopy perturbation method (HPM) [11,12], Homotopy perturbation transform method (HPTM), residual power series method (RPSM), etc., have been developed to determine the numerical solutions of several classes of fractional ordinary differential equations (ODEs) and partial differential equations (PDEs).
During the last decade, FDEs have gained popularity due to their extensive usefulness and relevance regarding the investigation of behaviour of real physical models. In the recent years, many authors have tried to solve various types of PDEs in different fields and obtained numerical solutions by using different approximation techniques. Some authors have acquired exact solutions of time-fractional PDEs by employing latest iterative technique [13]. Ceser et al. [14] have used a differential transform scheme to handle the system of FPDEs. In 2013, Neamaty et al. [15] utilized an iterative Laplace transform scheme and wavelet operational method, respectively, to handle a system of FPDEs. Babolian et al. [16] have used the combination of ADM and Spectral scheme to acquire the solution of FPDEs. In addition, the exact solutions of FPDEs by using the sub-equation method have been derived by Bekir et al. [17]. Wang et al. [18] have applied a fractional modified sub-equation scheme, and Gupta et al. [19] have employed Laplace transform to solve a system of FPDEs. Very recently a variety of FPDEs have been analysed by using Laplace VIM and Laplace ADM [20], modified HPM [21], HPM [22]. Neamaty et al. [23] have tried variational homotopy perturbation iteration method (VHPIM) and Zayed et al. [24] have utilized complex transformation for nonlinear FPDEs arising in mathematical physics.
Additionally, Singh et al. [25] have presented the analysis of coupled fractional Burger's equations by employing a homotopy technique. In 2017, FPDEs have been handled by using Muntz-Legendre polynomial approach [26], fractional complex differential transform scheme [27], and two-dimensional Laplace transform [28]. Recently the FPDEs with Caputo-Fabrizio fractional operator have been investigated through HPTM by Gomez-Aguilar et al. [29]. Singla et al. [30] have derived solutions of Nizhnik-Novikov-Veselov (NNV) system, Burger system and Navier-Stokes equations by using a generalized Lie symmetry approach. Generalized solitary wave solutions for fractional Klein-Gordon equation have been derived by an exp-function scheme [31]. In 2018, FPDEs were analysed by a B-spline polynomial approach [32] and hybrid Laplace transform technique [33]. Wang et al. [34] have gained solutions for FPDEs with proportional delay with iterative RPSM and HATM schemes.
Recently a bunch of new mathematical models have been studied by some authors with various singular and non-singular fractional derivative operators. Baleanu et al. [35] have handled the fractional Lagrangian and fractional Hamilton's equations of the movement of a particle in a round-shaped cavity with fractional derivative in Caputo sense and with Mittag-Leffler kernel. Recently some researchers have studied new mathematical models with fractional derivative operators related to the field of medical science and diseases. Baleanu et al. [36] have studied a newly generated fractional model for a tumor-immune surveillance with singular and non-singular derivative operators. In 2019, Jajarmi et al. [37] investigated a mathematical model of dengue fever outbreak based on FDEs. In this sequence, a new fractional mathematical model of diabetes and tuberculosis co-existence with a Mittag-Leffler non-singular derivative operator has also been analysed by Jajarmi et al. [38].
Over a decade ago, Baleanu et al. [39] derived the Lagrangians with linear velocities with RL fractional derivatives. Recently Baleanu and his co-workers have contributed to the classical and fractional study of physical systems with different types of fractional derivative operators. Baleanu et al. [40] have explored unusual characteristics of a physical model in the form of fractional Euler-Lagrange equations elucidating the motion of a capacitor microphone within non-singular derivative. Investigation of new features of a fractional model of spring pendulum and two coupled pendulums have been also carried out by Baleanu et al. [41,42]. Some time ago a group of authors studied the fractional optimal control problems by using certain approximation schemes. Mohammadi et al. [43] presented a numerical scheme based on hybrid Chelyshkov functions (HCFs) to handle a class of fractional optimal control problems (FOCPs). In this sequence, a new approximation approach has been utilized to investigate the nonlinear FOCPs by Jajarmi et al. [44].
This paper presents the execution of RPSM to compute the fractional power series (FPS) solutions of nonhomogeneous and homogeneous nonlinear systems of FPDEs arising in physical sciences as given below: Problem 1 Solve the nonhomogeneous nonlinear fractional system Problem 2 Solve the homogeneous nonlinear fractional system subject to the initial conditions where 0 < α ≤ 1 signifies the order of fractional time derivative.
Various definitions have been established for a fractional derivative. But this paper deals with a fractional derivative of Caputo type because it handles initial value problems (IVPs) in a very efficient way. The RPSM beautifully handles the systems of nonlinear ODEs and PDEs and computes the analytical solutions in a Taylor series form. It was first proposed and utilized by O.A. Arqub [45] to investigate the fuzzy differential equations so as to evaluate the coefficients of a power series solution. In addition, RPSM has been also applied to Lane-Emden equation, composite and non-composite DEs, regular IVP and fractional order boundary value problems (BVPs) [46][47][48][49]. Many authors have used RPSM to solve the problems of diverse streams of science and engineering [50][51][52][53][54][55].
Recently, Kumar et al. [56][57][58] have presented the analysis of the fractional exothermic reactions model, the fractional vibration equation and new fractional SIRS-SI malaria disease model with a variety of fractional operators with different kinds of kernel. Very recently a new investigation of Drinfeld-Sokolov-Wilson (DSW) equation was carried out by Bhatter et al. [59]. Singh et al. [60] have explored the new features of fractional Biswas-Milovic model with Mittag-Leffler kernel. In 2019, fractional Black-Scholes equations have been studied via RPSM by Dubey et al. [61].
The suggested nonhomogeneous and homogeneous systems of nonlinear fractional PDEs with Caputo derivatives, to the best of our knowledge, have not been handled utilizing RPSM in the available literature, yet. The main reason to choose the RPSM for application is that it effectively produces the approximate analytical solution of FDEs in the power series form with high accuracy and fast convergence. The primary purpose of the present work is to analyse the given systems of fractional PDEs with Caputo time derivatives through RPSM. The present work derives the fractional power series (FPS) solutions of the above-mentioned systems of FPDEs and further investigates the variations of numerical results obtained in a series form regarding time and several values of fractional parameter α through the two-dimensional and three-dimensional graphs.
There are various methods in the available literature that have been used to handle the above mentioned systems of FPDEs. Many researchers like Jafari et al. [62,63]; Wazwaz [64]; Aminikhah et al. [65] and Koçak et al. [66] have successfully utilized a variety of schemes to handle these systems of nonlinear fractional PDEs. But the RPSM has never been used to handle these systems of FDEs. The novelty of the present paper lies in the smooth implementation of RPSM to the systems of nonlinear PDEs arising in physical sciences, along with depiction of variation of numerical results regarding various values of fractional parameter α. The other novel factor in the proposed work is the handling of the set of nonlinear PDEs with Caputo fractional derivative (CFD). However, Riemann-Liouville (RL) fractional integral and derivative operators have been used significantly in the derivation of other fractional integral and derivative operators. But various problems of applied nature require conveniently suitable format of fractional derivatives, for instance, the Caputo derivative which delivers easy to deal with initial conditions having clear elucidation for the FDEs. This fact promotes the CFD as a more appropriate operator in comparison to the Riemann-Liouville definition regarding the application. However, it is noteworthy that the Caputo derivative demands the possibility of evaluation of the nth derivative of a function which makes its scope broader than its RL substitute. But a large class of mathematical functions that appear in applications satisfy this requirement.
This paper significantly presents the implementation of the RPSM to the systems of PDEs of nonlinear nature with Caputo fractional derivatives. This recently generated iterative scheme provides the solutions in analytical Taylor series form for the systems of ODEs and PDEs of both linear and nonlinear type. This semi-analytic technique actually utilizes the generalized Taylor series expansion accompanying the residual error function to build up the solution in a form of FPS expansion for a class of nonlinear FDEs without incorporating any additional restrictions. Utilizing the idea of residual error, a numerical solution is achieved straightforwardly in the form of truncated series. The primary benefit of this technique over the others is that it is directly applicable to the stated system with a proper pick of an initial guess approximation and it also lessens the complications originating in calculation of intricate terms. The remaining part of the manuscript is structured as follows: Sect. 2 describes the elemental concepts and formulae connected to the field of fractional calculus and theorems of fractional power series. Section 3 presents the fundamental procedure of RPSM. In Sects. 4 and 5, we present the application of RPSM to Problems 1 and 2, respectively, and related graphs are also presented there. Section 6 provides the numerical results and discussion. Finally, Sect. 7 records the concluding remarks.

Preliminaries and notations
This part formally presents some essential definitions related to the stream of the fractional calculus and some fundamental theorems related to the FPS expansion as follows: The fractional integral operator of order μ (μ ≥ 0) of g(y) of Riemann-Liouville (RL) type is defined as

Definition 2.2 ([67])
The fractional derivative of g(y), y > 0, in the RL sense is defined as

Definition 2.3 ([67])
The β-order fractional derivative operator of Caputo type for the function g(y), y > 0, is formulated as where D κ signifies the κ-order Newtonian derivative.
The following expressions hold true in the context of Caputo derivative: The fractional derivative of order α > 0 in the Caputo sense is stated as .
Here ζ signifies a variable and ℘ (z) are the coefficients of the expansion series.

Theorem 2.2 ([70, 71]) Let the multiple FPS expansion at
which is actually the formula of a generalized Taylor series. The classical Taylor series for γ = 1 is produced as

RPSM: fundamental procedure
In this section, the basic procedure of residual power series (RPS) algorithm is provided.
For the purpose of illustration of fundamental procedure of RPSM for fractional order PDEs, the ensuing fractional PDE is taken which is about to be mentioned here: subject to the conditions: where D λγ ζ = ∂ λγ ∂ζ λγ specifies the fractional derivative in Caputo sense; F[ϑ] and E[ϑ] characterize the nonlinear and linear operators in ϑ and H(ϑ, ζ ) denotes a continuous function.
As stated by RPS algorithm, the solution of equations (3.1)-(3.2) can be expressed as an FPS expansion about the point ζ = 0. Now it is assumed that the solution of the above mentioned PDE acquires the series expansion form as Therefore, the approximate analytical solution for equations (3.1) and (3.2) in the shape of an infinite FPS is provided by the RPS algorithm. In this sequence, the th truncated series of w(ϑ, ζ ), i.e. w (ϑ, ζ ) is defined as Clearly, the initial condition (3.2) is satisfied by w(ϑ, ζ ) thus from equation (3.3), the following algebraic expression is obtained: (3.5) From equation (3.4), the initial guess approximation of w(ϑ, ζ ) may be expressed as Now the expansion series of equation (3.4) can be rewritten as Now to compute the values of multipliers h λ (ϑ), λ = 1, 2, 3, . . . , in the series expansion of w (ϑ, ζ ) in equation (3.7), the residual function is expressed as (3.8) and the th residual function, Re s w, is expressed in the following form: Re s (ϑ, 0) = 0, 0 < γ ≤ 1, ϑ ∈ I, 0 ≤ ζ < , = 1, ∞. (3.10) Thus all the desired coefficients of multiple FPS of equations (3.1) and (3.2) can be easily calculated by endorsing the above-mentioned procedure.

Application of RPSM for Problem 2
The given homogeneous nonlinear system of FPDEs is subject to the conditions: The RPS algorithm assumes the following set of solutions for the above system of equations (1.2) as an FPS expansion about ζ = 0 as Clearly, u(ϑ, η, ζ ), v(ϑ, η, ζ ) and w(ϑ, η, ζ ) satisfy the initial conditions which can be writ- Hence the initial guess approximations of u(ϑ, η, ζ ), v(ϑ, η, ζ ) and w(ϑ, η, ζ ) are obtained Hence .

Numerical results and discussion
In this section, the numerical values of the functions u(ϑ, ζ ), v(ϑ, ζ ), u(ϑ, η, ζ ), v(ϑ, η, ζ ), and w(ϑ, η, ζ ) are computed for various values of fractional parameter α = 1 3 , 2 3 and classical motion α = 1. Figures 1-8 illustrate Problem 1, and Figs. 9-20 represent Problem 2. Figures 1 and 2 graphically show the two-dimensional variations of u(ϑ, ζ ), v(ϑ, ζ ) with  Figure 9 depicts that u(ϑ, η, ζ ) decreases with increasing time ζ and decreasing α for ϑ = 1. It is also observed from Figs. 10 and 11 that v(ϑ, η, ζ ) and w(ϑ, η, ζ ) increase with increasing ζ and decreasing α for ϑ = 1. Figures 12-20 graphically portray the three-dimensional representation of variations of u(ϑ, η, ζ ), v(ϑ, η, ζ ) and w(ϑ, η, ζ ) with respect to ϑ and ζ for the value η = 0.4. This study undoubtedly deduces that the achieved approximate analytical solutions continuously depend on the fractional parameter α. Two-and three-dimensional graphical demonstrations ensure the high accuracy of the generated results using RPSM by involving four iterations only. Nevertheless, the numerical results can be generated more precisely by computing additional iterative terms. The approximate analytical solutions and numerical results of the given sets of PDEs obtained via RPSM are in good correspondence with the results obtained via HAM [63], VIM [64], Laplace transform new homotopy perturbation method (LTNHPM) [65], and three-dimensional transform method [72]. The other methods and RPSM yield accurate solutions, but RPSM is simpler compared to other methods. The most significant feature of the RPS algorithm is that it solves a chain of single variable linear equations to compute the required coefficients of a power series solution which makes the solution procedure very easy and straightforward in comparison to other semi-analytical techniques.

Concluding remarks
This article extends the implementation of an iterative scheme based on the RPS algorithm to acquire approximate solutions for a nonlinear time-fractional nonhomogeneous and homogeneous systems of PDEs. The repercussions of the fractional order of Caputo time derivatives on the numerical solutions are well depicted through graphical representations for various specific cases. The study authenticates that RPSM supplies accurate numerical solutions without any restrictive assumptions on nonlinear PDEs. This technique has a clear advantage over other techniques in the sense that it yields approximate solutions of the problems by utilizing the contextual choice of an initial guess approximation. This technique delivers rapidly convergent series and also minimizes the complications that arise in the evaluation procedure of certain intricate terms. This paper beautifully checks the applicability strength and potency of RPSM to generate numerical solutions of a system of FPDEs. The stated systems of nonlinear FPDEs in this manuscript can also be investigated by the accurate discretization method which was earlier applied by Hajipour et al. [73] to one-, two-, and three-dimensional highly nonlinear Bratu-type problems. In this discretization method, nonlinear equations are discretized via a fourth-order nonstandard computational finite difference formula, and finally the given problem is reduced to the solution of a highly nonlinear algebraic system. In addition, the semi-implicit finite difference weighted essentially non-oscillatory (WENO) scheme of the sixth order utilized by Hajipour et al. [74] to solve the nonlinear heat equation can also be tried to investigate the described systems of FPDEs in this manuscript. These numerical schemes probably can provide new insights regarding application to FPDEs and further new conclusions in the future.