A Vector Series Solution for a Class of Hyperbolic System of Caputo Time-Fractional Partial Differential Equations With Variable Coefficients

In this paper, we introduce a series solution to a class of hyperbolic system of time-fractional partial differential equations with variable coefficients. The fractional derivative has been considered by the concept of Caputo. Two expansions of matrix functions are proposed and used to create series solutions for the target problem. The first one is a fractional Laurent series, and the second is a fractional power series. A new approach, via the residual power series method and the Laplace transform, is also used to find the coefficients of the series solution. In order to test our proposed method, we discuss four interesting and important applications. Numerical results are given to authenticate the efficiency and accuracy of our method and to test the validity of our obtained results. Moreover, solution surface graphs are plotted to illustrate the effect of fractional derivative arrangement on the behavior of the solution.


INTRODUCTION
Many natural phenomena have been modeled through partial differential equations (PDEs), especially in physics, engineering, chemistry, and biology, as well as in humanities [1,2]. A wide range of PDEs can be classified under the name of hyperbolic PDEs that have the following general form [2][3][4][5][6]: subject to the following initial condition: The equations of compressible fluid flow and the Euler equations are examples of PDEs that can be reduced to hyperbolic PDEs when the effects of viscosity and heat conduction are neglected [6]. In addition, many mathematical models are appearing as hyperbolic systems of PDEs that have the following general form: U t (x, t) A(x, t)U x (x, t) + B(x, t)U(x, t) + F(x, t), x ∈ I, t ≥ 0, subject to 1 Γ(α) t t0 (t − τ) α− 1 u(x, τ)dτ, t > τ > t 0 ≥ 0. (6) For more details about the properties of the two previous definitions, readers can refer to the references [7][8][9][10][11][12]. The most useful properties that we need in this research are As mentioned, the definition of Caputo is one of the most important definitions of the FD, since it was and still is an appropriate and effective tool in the modeling of many natural phenomena in all sciences and fields. For example, but not limited to, the definition of Caputo has recently been used to construct a mathematical model to illustrate the impacts of deforestation on wildlife species [13], in a fractional investigation of bank data [14], to model the spread of hookworm infection [15], and newly to model and analyze the dynamics of novel coronavirus (COVID-19) [16].
It is difficult to find exact solutions (ESs) for the fractional differential and integral equations; for this reason, analytical and numerical methods are usually used to find approximate solutions (ASs) for those equations. In recent decades, many methods have been used to find analytical and numerical solutions for fractional differential and integral equations such as the variational iteration method [17], the Adomian decomposition method [18], the homotopy perturbation method [17], the homotopy analysis method, the fractional transform method [19], Green's function method [20], and other methods [21,22].
In the last five years, the residual power series method (RPSM) has achieved an advanced rank among the methods used to find ASs for many fractional differential and integral equations. It has been used in determining ESs and ASs for many equations such as homogeneous and non-homogeneous time-and space-fractional telegraph equation [23], timefractional Boussinesq-type and space-fractional Klein-Gordon-type equations [24], fractional multipantograph system [25], space-and time-fractional linear and nonlinear KdV-Burgers equation [26], multi-energy groups of neutron diffusion equations [27], and other equations. The RPSM is characterized by its ease and speed in finding solutions for equations in the form of a power series. In fact, the RPSM is a mechanism for finding the coefficients of the fractional power series (FPS) without having to find a recurrence relation that we normally obtain from equating the corresponding coefficients in the series. The RPSM is a good alternate proceeding for gaining analytic solutions for fractional PDEs.
Despite all the features we mentioned about the RPSM, we will present in this paper a major modification to the method. We use the concept of limit at infinity instead of the concept of FD in determining the coefficients of the power series solution (SS). As is well known, finding an FD manually is not easy and sometimes takes tens of minutes when it is calculated by software packages, while calculating the limit is much easier than calculating the FD manually and faster by compute. Indeed, the RPSM determines the coefficients of the power SS of the differential or integral equations, whereas the proposed technique determines the coefficients of the expansion that represents the Laplace transform (LT) of the solution. Therefore, we do not need FDs during the transaction-finding process. To be able to implement the new method, we need to provide two expansions of the matrix functions, one to represent the solution of the equation and the other to represent the LT of the solution. Moreover, the convergence of the introduced expansions is studied. In fact, the proposed method called the Laplace-RPSM (L-RPSM) was first introduced by the authors in [28] and used for introducing exact and approximate SSs to the linear and nonlinear neutral FDEs. El-Ajou [29] then adapted the new method in creating solitary solutions for the nonlinear time-fractional partial differential equations (T-FPDEs).
Several articles are interested in providing ASs to T-FPDEs of hyperbolic type, such as the Caputo time-fractional-order hyperbolic telegraph equation [30], hyperbolic T-FPDEs [31][32][33][34][35], the time-fractional diffusion equation [36], fractional advection-dispersion flow equations [37], and other hyperbolic equations. However, a limited number of research studies provided analytical and numerical solutions for hyperbolic systems of T-FPDEs. Kochubei [38] presented a numerical-analytical solution for homogeneous hyperbolic fractional systems, and Hendy et al. [39] introduced a solution for two-dimensional fractional systems that was provided by a separate contrast scheme. Therefore, more research is needed in providing analytical and numerical solutions for such systems due to their importance in many applications as mentioned above.
The present work aims to apply the L-RPSM to construct ASs of a hyperbolic system of T-FPDEs with variable coefficients in the sense of Caputo's FD, which are given in the form of the following model: subject to where x U(x, t) refers to Caputo's space-FD of order ß of the multivariable vector function U(x, t), and the definitions of all terms and variables in Eqs 11, 12 are the same as those in Eqs 3, 4.
This paper is organized as follows: In Section 2, the analysis of matrix FPS is prepared. In Section 3, the construction of FPS solution to a hyperbolic system of T-FPDEs with variable coefficients in the sense of Caputo's FD is presented using the L-RPSM. In Section 4, application models and graphical and numerical simulations are performed in order to illustrate the capability and the simplicity of the proposed method. Finally, the conclusion is presented in Section 5.

PRELIMINARIES OF MATRIX FPS
Here, we present some definitions and theories regarding matrix analysis and the matrix FPS, which are playing a central role in constructing the L-RPSM solution to a hyperbolic system of T-FPDEs with variable coefficients.
Lemma 2.1. If m − 1 < α ≤ m and m ∈ N, then Definition 2.4 For 0 < α ≤ 1, x ∈ I, and t ≥ t 0 , a matrix power series of the following form: is called a bivariate matrix FPS around t 0 , where t is an independent variable and H m (x) ∈ M r×k are matrix functions of the independent variable x called series coefficients.
Proof. Of the operator definition in Eqs 6, 13 we have based on the second mean value theorem for integral [4] Proof. From Theorem 2.1, it suffices to demonstrate that According to Lemma 2.1, it is easy to show that the formula is correct for n 0 and n 1. Thus, inductively, we prove the theorem as follows: Thus, the proof of Theorem 2.2 has been completed. Let us call the series (Eq. 17) the bivariate fractional matrix Taylor's formula (BFMTF) of the matrix function U(x, t). As any series, the tail of the series (Eq. 17), R n (x, t) is called the nth remainder for the Taylor series of U(x, t). The function P(x, t) U(x, t) − R n (x, t) is an approximate function for U(x, t), and the accuracy of the approximation increases as R n (x, t) decreases. Finding a bound for R n (x, t) gives an indication of the accuracy of the approximation P(x, t) ≈ U(x, t). The following theorem provides such a bound.
Proof. The definition of the remainder of the BFMTF of U(x, t) as in Eq. 17 is given by According to Theorem 2.2, the remainder can be expressed as Thus, the proof is completed. Note that when n → ∞, Taylor's formula (17) is of the form which can be applied throughout this work. Finally, it is worth to mention that if α 1, then the BFMTF (Eq. 21) becomes which is the bivariate classical matrix Taylor's formula of a matrix function.
are of exponential orders (EOs) λ ij and piecewise continuous functions (PCFs) on I × [t 0 , ∞), respectively. Then, and of EOs λ ij , respectively. Assume that U(x, t) can be represented as a BFMTF as in Eq. 21. Then, the inverse LT of U(x, t) has the following fractional matrix expansion (FME): where λ minλ ij , 1 ≤ i ≤ r, 1 ≤ j ≤ k, which can be applied directly throughout this work when t 0 0.
and at a fixed point x, then the norm of the remainder of the FME in Eq. 23 satisfies Proof. As it is assumed in the text of the theorem, suppose that As in Eq. 19, the remainder of the FME in Eq. 23 is Multiplying Eq. 26 by s 1+(n+1)α , we get Finally, for 0 ≤ λ < s ≤ c and fixed x, we have Thus, we reach the end of the proof.

APPLYING THE L-RPSM TO THE HYPERBOLIC SYSTEM OF T-FPDES
In this section, we construct an AS to the hyperbolic system of T-FPDEs with variable coefficients given in Eqs 11, 12 by using the L-RPSM. To achieve it, firstly, we apply the LT on both sides of Eq. 11, and use Lemma 2.2, and Eq. 12; then, we have where Let the solution of Eq. 30 have the following FME: where Of course, treating with a finite series is acceptable more than an infinite series. For this reason, the L-RPSM deals with a finite series while calculating coefficients of the SS. So, we express the kth truncated series (kth TS) of U(x, s) as follows: To apply the L-RPSM for determining the coefficients H m (x), m 1, 2, 3, . . . , k, in the kth TS in Eq. 32, we define the so-called residual matrix function (RMF) for Eq. 30 as and the kth residual matrix function (RMF k ) of the style form The main idea of the L-RPSM can be shown in the following clear facts related to the RMF and RMF k : 1. lim k → ∞ RMF k (x, s) RMF(x, s), x ∈ I, s > λ ≥ 0 2. RMF(x, s) 0 ∈ M r×1 , x ∈ I, s > λ ≥ 0 3. RMF(x, s) has an FME. So, we can express it as follows: where N m , m 1, 2, 3, . . ., are operators depending on the operators L and z β x .
6. Using the following fact determines the unknown coefficients H k (x) , k 1, 2, 3, . . ., in the FME (Eq. 31): Now, to find H 1 (x), in Eq. 32, substitute the 1st TS, U 1 (x, s) U0(x)) s + H1(x)) s 1+α , into the 1st RMF, RMF 1 (x, s), to get Multiply Eq. 38 by s 1+α to obtain Compute the limit to Eq. 39 as s → ∞, use the fact in Eq. 37, and solve the new obtained equation for H 1 (x) to have Similarly, to determine the second unknown coefficient in Eq. 32, H 2 (x), we substitute U 2 (x, s) U0(x)) s + H1(x)) s 1+α + H2(x)) s 1+2α into RMF 2 (x, s) to get Frontiers in Physics | www.frontiersin.org RMF 2 (x, s) Multiply Eq. 41 by s 1+2α and compute the limit at infinity for both sides of a new obtained equation, according to Eq. 37, to have In general, to determine the nth unknown coefficient in Eq. 32, H n (x), we substitute U n (x, s) U0(x)) s + H1(x)) s 1+α + . . . + Hn(x)) s 1+nα into RMF k (x, s) for k n, re-multiplying both sides of the new obtained formula by s 1+nα , and use the fact in Eq. 37 to obtain This procedure can be repeated for the required number of FME coefficients representing the solution of Eq. 30. Therefore, the kth approximation of the solution of Eq. 30 can be represented as the following finite series: If we act the inverse LT on both sides of Eq. 44, then we obtain the kth approximation of the solution of the initial value problem (IVP) (Eqs 11, 12), which takes the following expression: (45)

APPLICATIONS AND NUMERICAL SIMULATIONS
To test our proposed method, we present in this section four interesting and important applications. The first three applications are prepared so that the ES is already known, while the last application is prepared without knowing the solution in advance to test the predictability of the solution or obtain a suitable AS. Application 4.1. Consider the following homogeneous hyperbolic system of T-FPDEs with variable coefficients: subject to where To obtain an FME solution for this application using the L-RPSM, transform Eq. 46 to the Laplace space as follows: Let the solution of Eq. 48 have a form of the FME as in Eq. 31. According to the condition in Eq. 47, the first coefficient of the FME in Eq. 31, H 0 (x) U(x, 0) x 1 . Therefore, the kth TS of Eq. 31 takes the following expression: and the kth RMF of Eq. 48 is To find the first unknown coefficient Eq. 49, we put the 1st TS, U 1 (x, s) x 1 1 s + H 1 (x) s 1+α , into the 1st RMF to get the following abbreviated expression: Take the limit at infinity for Eq. 52 and use Eq. 37 to get Similarly, we can obtain the second unknown coefficient in Eq. s) to get the next summarized expression of the second RMF of Eq. 48: Multiply Eq. 54 by s 1+2α , apply the limit at infinity of the obtained equation, and use Eq. 37 to get If we repeat the previous procedure for k 3, 4, . . ., we can see that So, the ES for Eq. 48 will be as follows: If we apply the inverse LT on Eq. 57, then the SS of the IVP (Eqs 46, 47) will take the following form: which coincides with the ES. Application 4.2. Consider the following non-homogeneous hyperbolic system of T-FPDEs with variable coefficients: where According to the construction in Section 3, the LT of Eq. 59 can be represented by The kth TS of the FME of the solution of Eq. 61 has the following form: and the kth RMF of Eq. 61 is (63) So, to set the first unknown coefficient of Eq. 62, substitute Frontiers in Physics | www.frontiersin.org and use the result in Eq. 37 to obtain Again, substitute U 2 (x, s) into s 1+2α RMF 2 (x, s) to get According to the fact in Eq. 37, we have Similarly, anyone can check that H k (x) 0, for k 3, 4, . . .. So, the ES for Eq. 61 can be expressed as Thus, the FME solution of the IVP (Eqs 59, 60) would be as follows: which is identical to the ES U(x, t) e x t α x 2 + t 2α . Application 4.3. Consider the following non-homogeneous hyperbolic system of T-FPDEs with variable coefficients: subject to where and the ES is where E α (t) is the Mittag-Leffler function defined by the following expansion [40]: Mathematica 7 software has been used through a low-RAM PC for obtaining all numerical calculations and symbolism. Since the Mittag-Leffler function is an infinite expansion, it was difficult to perform the calculations using the Mittag-Leffler function as it is. For this, the fifth truncated series of the expansion in Eq. 73 was used throughout the calculations.
Like the previous applications, transform Eq. 70 to the Laplace space using the initial condition in Eq. 71 to read as follows: Let the solution of the algebraic equation (74) has an FME as in Eq. 31. Then, the kth TS of the FME of U(x, s) can be given by and the kth RMF of Eq. 74 is given by where H 1j (x; α, β) ∈ M 2×1 , j 1, 2, . . . , 7, are vector functions free from s. So, according to Eq. 37, we have Using the same previous approach, we find the following vector coefficients of Eq. 75: So, the fifth AS of Eq. 74 can be written as follows: Transforming the AS in Eq. 80 to the t-space by the inverse LT, we get the fifth approximation of the solution of the IVP (Eqs 70, 71) as follows: Obviously, there is a pattern between the terms of Eq. 81 that gives us the ES as in Eq. 72.
The mathematical behavior of the solution of the IVP (Eqs 70, 71) is illustrated next by plotting the three-dimensional space figures of the fifth approximation of the two components of the vector solution in Eq. 81 for different values of α and a fixed value of β 0.5. Figures 1A-C show the fifth AS, (U 1 ) 5 (x, t) and (U 2 ) 5 (x, t), when α 0.7, α 0.85, and α 1, respectively, on the square [0, 1] × [0, 1]. Figure 1D shows the ES expressed by Eq. 72 for α 1. Figures 1C,D show that the fifth AS of the IVP (Eq. 70, 71) is excellent compared to the ES, as well as in the previous cases, which have not been documented in order not to increase the numbers of graphs. It is known that, by increasing the number of terms in the series, the accuracy of the solution increases and, thus, the error of solution reduces; therefore, we can reduce the error of the solution by calculating more coefficients of the FME solution as in Eq. 31.
In the next application, the ES is unknown. Therefore, we are trying to find the ES or an appropriate approximation of the solution.
Application 4.4. Consider the following non-homogeneous hyperbolic system of T-FPDEs with variable coefficients: subject to where Similar to the previous applications, the LT of Eq. 82 is given by the kth TS of the expansion of the solution of Eq. 84 is given as and the kth RMF of Eq. 84 is given by According to the fact in Eq. 37, we can create, successively, the following first eight coefficients of the expansion in Eq. 31 for this application: (91) Tables 1, 2 show the values of ||RES 6 (x, t)|| and ||RES 7 (x, t)||, respectively, for different values of α. The data in the tables indicate that the norm of the residual error of the obtained AS decreases as (x, t) → (0, 0) as well as when α → 1. This indicates that the convergence of the BFMTF in Eq. 17 depends on t, x, and α as illustrated in Theorem 2.3. As we know, we can reduce the error in the FME solution as we increase the number of terms of the expansion. As we can see from the data in Tables 1, 2, the seventh approximation is more accurate than the sixth approximation. Anyway, it can be said that the L-RPSM is good at providing an accurate AS of a hyperbolic system of T-FPDEs with variable coefficients.

CONCLUSION
We have found that the ES for the hyperbolic system of T-FPDEs with variable coefficients is available if the solution is a linear combination of power functions or if it is a composite of an elementary function and a power function. In case the ES is not available, a good approximation of the solution can be obtained. The L-RPSM is an effective, accurate, easy, and speed technique in obtaining the values of coefficients for the SS. Through this work, we have presented a solution that may be missing for this kind of problem and we have opened the way for researchers to provide other ways to solve this class of equations. Moreover, the newly proposed technique can be used to construct many types of the ordinary or partial DEs of fractional order such as Lane-Emden, Boussinesq, KdV-Burgers, K(m, n), Klein-Gordon, and B(l, m, n) equations.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/ supplementary files.

AUTHOR CONTRIBUTIONS
The idea of this work, implementation, and output of this form was carried out by both authors.