Numerical study of multi-order fractional differential equations with constant and variable coefficients

In this manuscript, a numerical method based on the conjunction of Paraskevopoulos's algorithm and operational matrices is developed to solve numerically the multi-order linear and nonlinear fractional differential equations. By means of this conjunction, the multi-order problem is decomposed into a system of differential equations of fractional order which are then solved by employing operational matrices approach. The accuracy and efficiency of the method is examined by taking some examples. In addition, the numerical results presented in Pak et al. [Analytical solutions of linear inhomogeneous fractional differential equation with continuous variable coefficients. Adv Differ Equ. 2019;2019(256):1–22] are improved in our work.


Introduction
Fractional calculus (FC) is a study of derivatives and integration of arbitrary order. At the initial phases of its development, it was considered as an abstract mathematical idea with nearly no applications. But, now the situation is entirely different with FC and it considered to be the most useful and important topic among the scientific community. The various important phenomena of science and engineering have been well described with the use of fractional derivatives, including, partial bed-load transport, diffusion model, dynamics of earthquakes, viscoelastic systems, biological systems, chaos, and wave propagation, (see [1][2][3][4][5][6][7]) and references therein. Furthermore, FDEs have been used in the modelling of important physical phenomena appearing in the study of control theory, electrostatic, study of polymers, visco-elasticity, signal and image processing phenomena, computer networking, and mathematical biology (see [8][9][10][11][12][13]) and references therein.
The exact solution of most of the fractional differential equations (FDEs) is difficult to find due to the involvement of the fractional differential and integral operators in the problem. So, it is a natural way to develop the numerical methods to solve them numerically. Operational matrices approach coupled with orthogonal polynomials is one of the efficient and stable numerical tool used to find the approximate solutions of fractional ordinary and partial differential equations (see [14][15][16]). This approach is easy in use which has the ability to transform the FDEs into a system of algebraic equations. Finally, this approach approximates the solution as basis vectors of orthogonal polynomials (see [17][18][19][20][21]).
Motivated by aforementioned studies, we propose a numerical method based on the extension of Paraskevopoulos's method in conjunction with shifted Legendre polynomials (SLPs) to solve numerically the multiorder linear and nonlinear FDEs with constant and variable coefficients. Under this approach, the multiorder fractional problem is decomposed into a system of FDEs which are then solved by using operational matrices approach. Finally, the solution is approximated as a basis vectors of orthogonal SLPs. We also compare the numerical results obtained using our method with the results obtained in [22]. We observe that our proposed method improves the numerical results presented in [22]. For instance, (see Example 7.5, Figure 1 and Table 1).
The manuscript is organized as follows: Some essential definitions, properties and results of FC are given in Section 2. In Section 3, the properties of SLPs are recalled. In Section 4, a generalized operational matrix in the sense of Caputo fractional differential operator is studied. In Section 5, the multi-order fractional problems are reduced into a system of FDEs. In Section 6, a numerical algorithm based on Paraskevopoulos's method in conjunction with operational matrices approach is developed to solve linear and nonlinear multi-order FDEs with constant and  variable coefficients. The numerical accuracy of the proposed method is examined by taking some examples in Section 7. Finally, in Section 9, the manuscript ends with a brief conclusion and some remarks.

Preliminary remarks on FC
Some important definitions and results are recalled in this section which are indispensable to construct the proposed numerical algorithm.
For the Caputo derivative, we have the following observations: Note the following basic results: Proof: By definition of Caputo differential operator, we have where m is an integer defined as m − 1 < α < m. Now, according to our assumption, n ≥ m which implies that u (m) (t) is a continuous function, and u (m) (t) ∈ C n−m [0, 1]. Consequently, for t → 0, the integral in (6) vanishes. Hence, C D α u(0) = 0.

Properties of SLPs
The following recurrence formulae are used to evaluate the orthogonal Legendre polynomials (LPs) on the interval of orthogonality [−1, 1] (see [23]): Now by setting, x = 2t−1 in (7), the SLPs on the interval of orthogonality [0, 1] can be expressed with the aid of the following analytical form, (see [23]) The orthogonality condition of SLPs can be expressed as

Functions approximation
Suppose u(t) ∈ L 2 [0, 1], then it can be demonstrated as a basis vectors of SLPs, (see [23]) Using (9), the coefficients e i can be determined as For the practicality, considering the first n + 1-terms of (10), we have where and
and for m even, we have

Lemma 4.2:
Let P j (t) be the SLPs as defined in (8), then Proof: Using the definition of SLPs, we have Using (17), the SLPs of degree α − 1 is demonstrated as Applying Caputo differential operator, and its linearity to (18), we have Since, α > α − 1, then using the results of (4)-(5), we have Consequently, Theorem 4.1: Let (t) be a function vector of SLPs and suppose α > 0, then where D (α) is a generalized derivative operational matrix in Caputo sense, demonstrated as and ϒ j,l,k is given as Proof: Using (8), (5), and Lemma 4.2, we have We can approximate t k−α by first n + 1-terms of SLPs, as where By inserting, (26) into (25), we have e l,k P l (t).
(28) By putting the value of e l,k in (28), and after some simplification, we have For simplicity, we have Now (29) can be expressed as (31) can be written in vector notation as Moreover, using Lemma 4.2, we have (33) Equations (32) and (33) together prove the required result.

Decomposition of multi-order FDEs into a system of FDEs
Consider the multi-order FDEs as can be a nonlinear in general, and C D α is a Caputo derivative of order α defined in (3). Problem (34) can be decomposed into a system of FDEs, as Set u 1 = u and suppose, Hence the claim is verified. If β 1 / ∈ N, then by Lemma 2.1, C D β 1 u 1 (0) = 0, and as Hence u 3 = u (m) . Further, we define: Claim: The process will be continued until the decomposition of the problem (34) into a system of FDEs.

Implementation to the problem
In this section, a multi-order fractional differential equation is decomposed into a system of FDEs using the algorithm developed in Section 5.
Consider the following multi-order fractional differential equation subject to the following initial conditions Set, u(t) = u 1 , we have

Application of operational matrices method of SLPs
In this section, we apply the Legendre operational matrix method to find approximate solution of linear and nonlinear multi-order FDEs by converting them into a system of FDEs.

Algorithm for the linear multi-order FDEs
Consider the linear system of FDEs established after implementation of the method studied in Section 5: with initial conditions u (j) where γ ij , k i , and d ij are real constants. The differential operator C D α represents the Caputo fractional derivative of order α. We can approximate u i (t) and g i (t) by SLPs as, Here , e i1 , . . . , e im , ] is unknown vector. Using (22) and (44), we have Employing Equation (44)-(45) and calculating the residuals R i m (t) for the system (42), it yields Now using typical Tau method [25], we can generate (m − p i + 1) linear algebraic equations by applying the following: Equations (48), and (49) together generate n(m + 1) linear algebraic equations, corresponds to (42). By solving these equations for unknown coefficients of the vector T i , we consequently obtained the solutions u i (t) for (i = 1, 2, . . . , n) of (42) under initial conditions (43).

Algorithm for the nonlinear multi-order FDEs
In this section, we develop the numerical algorithm to solve nonlinear multi-order FDEs with the support of the theory developed in Section 5: Consider, together with initial conditions defined in (43). Here F i generally a nonlinear function for i = 1, 2, . . . , n. Substituting the approximation of u i (t) and C D α i u i (t) into above Equation (50), we have Equations (51) are exactly satisfied at the first m − p i + 1 roots of P m+1 (t). So, we collocate these equations at root points. Consequently, these equations together with (49) produce n(m + 1) nonlinear algebraic equations with unknown coefficients of T i . Finally, approximated solution u i (t) for (i = 1, 2, . . . , n) can be derived by solving these nonlinear algebraic equations using Newton's method.

Illustrative examples
This section is devoted to demonstrate the applicability of our developed numerical method by taking some test examples. We solve numerically the multi-order linear and nonlinear FDEs in this section. The numerical results are demonstrated by using tables and plots. The software MATLAB is used for numerical calculations and simulations.

Example 7.1: Consider the following multi-order fractional Bagley-Torvik equation
subject to the integer-order initial conditions The source term G 1 is given as The exact solution at β 1 = 2, β 2 = 3 2 is given below Using the algorithm studied in Section 5, the problem (52)-(53) can be converted into a system of FDEs, as Then the functions u 1 (t), and u 2 (t), and the source term G 1 (t) can be approximated using the first three terms of SLPs as (56) The residuals for the system (54) can be evaluated using Equation (55), and the fractional operational matrix as Finally, the solution can be expressed as = e 10 P 0 (t) + e 11 P 1 (t) + e 12 P 2 (t) Remark 7.1: By increasing the number of terms of SLPs, the exact solution and the approximate solution of the problem, (52)-(53) coincide. For example, by considering the first four terms of SLPs, we have = e 10 P 0 (t) + e 11 P 1 (t) + e 12 P 2 (t) + e 13 P 3 (t) The maximum absolute error is zero because the exact solution and the approximate solution are same. Similarly by considering the first five terms of SLPs, we have the same result.

Example 7.2:
Consider the following inhomogeneous nonlinear fractional differential equation subject to the integer-order initial conditions The source term G 1 is given as where β 1 = 3, and β 2 = 2.5. The exact solution is given below Using the algorithm studied in Section 5, the problem (62)-(63) can be converted into a system of FDEs, as Table 2. Legendre series coefficients (e 10 , e 11 , e 12 , e 13 ) for Example 7.2 obtained by using our proposed method (PM) and the Jacobi series coefficients (e 10 = c 0 , e 11 = c 1 , e 12 = c 2 , e 13 = c 3 ) obtained by using the method ( [19], Example 4), when α = β = 0, and n = 3 = scale level = the number of terms of SLPs.
Coefficients PM Method proposed in [19] Tables 2-4 show that the results obtained by using our PM are in a good agreement with the results obtained by using the method in ( [19], Example 4).

Example 7.3: Consider the following inhomogeneous nonlinear FDE
subject to the integer-order initial conditions The source term G 1 is given as where β 1 = 1.5, and β 2 = 2. The exact solution is given below Using the algorithm studied in Section 5, the problem (65)-(66) can be converted into a system of FDEs, as

Example 7.4: Consider the following inhomogeneous FDE
subject to the integer-order initial conditions The source term G 1 is given as where 1 < β 1 ≤ 2, and β 2 = 1. The exact solution at β 1 = 2 can be easily determined, given as Using the algorithm studied in Section 5, the problem (68)-(69) for a fixed value of β 1 = 1.75, can be converted into a system of FDEs, as Then the functions u 1 (t), and u 2 (t) can be approximated by using the first three terms of shifted Lagendre polynomials as Now using Equation (51), the system (70) can also be written as Considering the initial conditions and using Equations (48) and (73) Finally, the solution is approximated as = e 10 P 0 (t) + e 11 P 1 (t) + e 12 P 2 (t) = 0.73662t 2 + 0.27244t + 0.00001 Example 7.5: Consider the following inhomogeneous nonlinear FDE subject to the integer-order initial conditions The source term G 1 is given as where β 1 = 1.5, and β 2 = 0.8. The exact solution is given below Using the algorithm studied in Section 5, the problem (79)-(80) can be converted into a system of FDEs, as

Results and discussion
We check the accuracy and stability of our developed numerical algorithm by solving different examples. The exact and numerical solutions are compared at different scale levels and for different choices of fractional  We examine that as fractional order approaches to the exact order, the approximate solution approaches to the exact solution, see Figure 4. It is also observed that with the increase in scale levels, the approximate solutions are in a good agreement with the exact solutions, see Figures 2, 5, and Table 5. We also compute the absolute errors and observe that with the increase in scale levels, the amount of absolute errors decreases significantly, see Figure 3, Tables 7, and 1. It is worth to mention that even at low scale level, our approximate solution is in a good agreement with the exact solution, see Figure 6. We determine the analytical solution for each example using our proposed numerical algorithm. For example, in Example 7.1, the analytical solution is obtained for fractional order derivatives using the proposed numerical algorithm. The obtained numerical results are     then compared with the exact solution to test the accuracy of the proposed method. Similarly, the analytical solutions for Examples 7.2 and 7.3 are computed for fractional order derivatives using the same procedure as used in Example 7.1. However, in Examples 7.4 and 7.5, the analytical solutions are determined for integer order derivatives to test the accuracy and applicability of the proposed numerical algorithm for different choices of the fractional parameter, β 1 . We observe that the results computed using the proposed algorithm are in a good agreement with the exact solution at various values of β 1 .

Conclusion
We proposed an approximate method to solve linear and nonlinear multi-order FDEs by extending the Paraskevopoulos's method for orthogonal SLPs. Our proposed method is based on decomposing the multiorder FDEs into a system of FDEs and then the resultant system is solved using the operational matrix approach. We checked the accuracy and efficiency of the method by solving various linear and nonlinear FDEs with constant and variable coefficients. We examined that with the increase of scale level, the approximate solutions were in a good agreement with the exact solutions. We also demonstrated the high efficiency of the method by determining the amount of absolute error and observed that as we increased the scale level, this amount was decreased significantly. It is important to note that only considering the few terms of SLPs, we obtained the satisfactory results. In addition, our proposed method has advantage to other methods, like Homotopy perturbation method, because in our case, the perturbation, linearization, or discritization are not necessary to be implemented. We also improved the numerically obtained results presented in [22]. Finally, our proposed method is fit to solve both linear and nonlinear problems of fractional order with constant and variable coefficients.