A novel Jacobi operational matrix for numerical solution of multi-term variable-order fractional differential equations

In this article, we introduce a numerical technique for solving a class of multi-term variable-order fractional differential equation. The method depends on establishing a shifted Jacobi operational matrix (SJOM) of fractional variable-order derivatives. By using the constructed (SJOM) in combination with the collocation technique, the main problem is reduced to an algebraic system of equations that can be solved numerically. The bound of the error estimate for the suggested method is investigated. Numerical examples are introduced to illustrate the applicability, generality, and accuracy of the proposed technique. Moreover, many physical applications problems that have the multi-term variable-order fractional differential equation formulae such as the damped mechanical oscillator problem and Bagley-Torvik equation can be solved via the presented method. Furthermore, the proposed method will be considered as a generalization of many numerical techniques.

In all modelled problems via the fractional differential equations (FDEs), the main concern of the mathematician researchers is to obtain the analytical or the numerical solutions for these equations. But, in many problems, obtaining the exact solution is more complicated. Consequently, there are great efforts of researchers to acquire accurate techniques for gaining a numerical solution for FDEs. Various techniques have been dragged special interest, such as in [40][41][42]. Also, one of the most familiar methods that is used is the spectral methods which lead to accurate approximate solutions because their basis is considered as a linear combination of the orthogonal polynomials (see, [43][44][45]).
Moreover, the spectral methods for solving the fractional-order differential equations basically depends on a set of orthogonal polynomials. One of these polynomials is the classical Jacobi polynomials that indicated by P (α,β) n (x) (n ≥ 0, α > −1 & β > −1), [46]. Jacobi polynomials have been used extensively in practical applications and mathematical analysis because it has the advantages of gaining the numerical solutions in α and β parameters. Thus, it would be beneficial to carry out a systematic study via Jacobi polynomials with general indexes α and β and this clearly considered one of the aims and the novelty of this manuscript in addition to the extension of the time interval t ∈ [0, l]. Furthermore, in recent decades, numerical study for VFDEs is widely used [47,48]. Therefore, it appeared a large number of techniques for obtaining the numerical solution of these equations (see for instance [49][50][51].
Here, the basic idea of this article generalization of the orthogonal polynomials in the base of solution by introducing a new shifted Jacobi operational matrix of the fractional derivative for solving numerically the multi-term variable-order FDEs in the following form: where D μ(t) u(t) and D γ i (t) u(t) (i = 1, 2, . . . , k) are the variable-order fractional derivatives defined in the Caputo sense. Note that: if μ(t) and γ i (t) (i = 1, 2, . . . , k) are constants, then Equation (1) have the form: Also note that: our new suggestion method considered as a generalization for every numerical solution of the multi-term fractional-order problem that can be solved by many polynomials such as all Chebyshev polynomials, Legendre polynomials, Gegenbauer polynomials, Lucas polynomials, Vieta-Lucas polynomials, and Fibonacci polynomials.
In what follows: begin by reviewing certain definitions and properties of the variable-order fractional derivatives in Section 2. In Section 3, the Jacobi polynomials and their shifted ones are introduced in addition to the matrix formula of the approximate solution via the shifted Jacobi polynomials. In Section 4, the shifted Jacobi operational matrix (SJOP) of the variable-order derivatives is derived. In Section 5, the bound error of the investigated method is obtained. The numerical results of the suggested method are obtained and compared with other methods 6. Conclusions are obvious in the last Section.

Preliminaries
In this preparative section, we introduce briefly main required definitions of the Capout's fractional derivative in variable-order functions which help us in what follows.

Definition 2.1 ([19,48]):
The Caputo variable-order fractional derivatives definition for the function u(t) ∈ C m [0, b] can be defined as: At the beginning time and 0 < μ(t) < 1, we have: where a and b constant amounts. According to Equation (3), as in [19,48] we have:

Jacobi polynomials and their properties
The beginning of this section by introducing essential basic facts of the Jacobi polynomials and their shifted ones, in addition to deriving some important tools that help us for developing the suggested method.

Jacobi polynomials and their shifted ones
The well-known Jacobi polynomials, P (α,β) n (x), are orthogonal polynomials of degree n in x defined on [−1, 1] as the following analytical form [46]: these polynomials given in Equation (8) can be generated by: where a α,β With the starting values For using the polynomials of Equation (8) on the interval t ∈ [0, l], modification of the variable x = (2t/l − 1) will be implemented. Hence, we have so-called the shifted Jacobi polynomials P (α,β) n (2t/l − 1) which indicated by P * (α,β) n (t). Thus, the analytical formula for these polynomials is given by: where P * (α,β) n and Also, P * (α,β) n (t) can be obtained via the recurrence relation: Where b α,β With the starting values P * (α,β) 0 (t) = 1 and Moreover, P * (α,β) n (t) satisfy the orthogonality relation on [0, l] as: where ω (α,β) = x β (l − x) α is the weight function of P * (α,β) n (t). Furthermore, since the shifted Jacobi polynomials that are given in Equation (12) with general indexes (α > −1 & β > −1). Then can claim of these polynomials; the most usually used orthogonal polynomials in numerical methods; the shifted Legendre polynomials L * n (t); the shifted Gegenbauer polynomials G * (α,β) n (t). In addition to the shifted Chebyshev polynomials of the first, second, third and fourth kinds which denoted by T * n (t), U * n (t), V * n (t) and W * n (t), respectively. All these orthogonal polynomials are contacted with P * (α,β) n (t) through the relations:

Function approximation via shifted Jacobi polynomials
can be expanded as the following expression: [46] where c i are the coefficients of the series. Hence, by taking n−terms of the series in Equation (19) we obtain: where in Equation (20)

is a vector and it's enters elements have values given by
(21) Now, if we suppose Then, according to Equation (22), the vector φ(t) can be expressed as: where A (α,β) is (n + 1) × (n + 1) of Equation (23) is a square matrix that specified by: For instance, if n = 4, α = β = 0 then A (matrix of shifted Legendre polynomials) is given by: Therefore, using Equation (23), we obtain Remark 3.1: For obtaining the square matrix A for all other orthogonal polynomials that related to the shifted Jacobi polynomials such as shifted Legendre polynomials, shifted Gegenbauer polynomials and shifted Chebyshev polynomials of all kinds, the coefficient of these polynomials according to the relation given in Equation (18) must be taken into consideration. For example, if n = 4, α = 1 2 , β = −1 2 then the square matrix A that derived in Equation (24) became (matrix of the fourth kind shifted Chebyshev polynomials) as follows:

Shifted Jacobi polynomials operational matrix (SJOM)
The shifted Jacobi operational matrix of fractional variable-order for supporting the numerical solution of the main problem will be investigated in this section. Hence, the problem will be transformed into the algebraic system of equations which solved numerically at the collocation points.
Firstly, D μ(t) φ(t) can be deduced as the following: Using Equation (7) in Equation (27), we obtain: Where Using Equation (25), we have . Now, we can obtain the fractional variable order the approximated function that given in Equation (20) as the following: Secondly, D γ i (t) φ(t), i = 1, 2, . . . , k can be obtained by follow the same way that proposed for obtaining the (SJOP) of D μ(t) φ(t), as follows: Where The square matrix Via, Equations (31) and (34), then Equation (1) converted to: Now, using t i = l(2i + 1)/(2n + 2), i = 0, 1, . . . , n. Then Equation (35) can be turned into the following system: Here, the system in Equation (36) in addition to the initial and boundary conditions can be solved numerically to determine the unknown vector C. Hence, the numerical solution that defined in Equation (20) can be calculated.

Error estimate
Theorem 5.1: Consider u(t) ∈ [0, l] be n-th times continuously differentiable and u n (t) be the best square approaches of u(t) given in Equation (20) then, we claim where Proof: Using Taylor expansion for u(t) as follows: where t 0 ∈ [0, l] and ξ ∈]t 0 , t[. Assumẽ then According to, u n (t) that given in Equation (20), we obtain (42) can rewritten as: Since, ω (α,β) = t β (l − t) α then, Hence, the proof is completed.
Corollary 5.1: Using Equation (37) in case of α = β = 0, then I = l. This implies that we can obtain the error bound in case of using the shifted Legendre polynomials [50]. (37) if α = β = 1 2 , then I = l 2 π/8 which is in complete agreement with the error bound that obtained via the second kind shifted Chebyshev polynomials which introduced by Liu and et al. [48]. (37)

Numerical experiments
To illustrate the accuracy, applicability, and generality of the presented method, several test examples are carried out in our hand section. Also, these examples support our theoretical discussion and give more general results where the indexes of the shifted Jacobi polynomials α, β carried out the results at different values of them. For comparing our computational results, the maximum absolute error E max between the analytical and numerical solutions will be used. Where E max defined as the following:

Example 6.1 ([48,49]):
Consider the multi-term fractional variable-order differential equation Where The analytical solution is u(t) = (2 − t 2 /2). This example will be solved in different cases as follows: Through the investigated technology in Section 4 with (n + 1) finite terms of the approximate solution that given in Equation (20) then we obtain Here, using the collection points t i that defined in Equation (45) for solving numerically the system of equations that given in Equation (49) via the collocation method. Then, we claim Hence, the system of Equations (50) solved numerically via Newton's iteration method and the vector C = [c 0 , c 1 , . . . , c n ] T is obtained. Therefore, the numerical solution in Equation (20) is calculated. Now, solving the algebraic system of Equations (50) with n = 2 and l = 1. Thus, we can write Table of the numerical results will not be presented in Case 1 because the exact solution is obtained. The three unknown coefficients with several choices of α and β are presented in Table 1 ( Figure 1). Moreover, Table 2 shows E max that obtained via the suggested method in the case of α = β = 0 and those given in [48] with several values of l and n for the fractional variable-order derivatives μ(t), γ 1 (t), γ 2 (t) and γ 3 (t) that given in Case 1. From these results, we can deduce that these results identically the analytical   [48] and proposed technique at α = β = 0 for Example 6.1 Case 1 with disjoint of l and n.
Method in [48] Our proposed method solution and more accurate than these results that obtained by Liu and et al. [48]. Also, the method in [48] depends on the second kind shifted Chebyshev polynomials which consider as a particular case of our suggested technique that depends on the shifted Jacobi polynomials. Furthermore, as seen from Table 3, the vector C T = [c 0 , c 1 , . . . , c n ] obtained is mainly composed of three terms, namely, [c 0 , c 1 , c 2 ]. Also, E max at l = 1 and distinct choices of n, α and β for the fractional variable-order derivatives μ(t), γ 1 (t), γ 2 (t) and γ 3 (t) is shown in the same Table. Table 3, shows achieving a full agreement with the analytical solution of the given problem via a few terms of P * (α,β) n (t).

Remark 6.1:
If α = 1 2 , β = −1 2 and l = 1, 2&4, then the numerical results that obtained by our technique for Example 6.1 Case 1 is in complete agreement with the results that obtained by authors in [49]. This case considered as a particular case of Case 1 because the multi-term fractional variable-order differential equation that given in Equation (46) becomes multi-term constant-order fractional differential equation. Table 4 gives E max that gained through the introduced technique for Example 6.1 Case 2 at distinct values of n, α = β = 0 and l = 1. From the first look at Table 4, the suggested gives highly accurate results comparing with others [40,44,48]. Table 5, presented Table 3. Computational results of Example 6.1 Case 1 for disjoint choices of α and β at l = 1..  Table 4. Comparison the E max for the methods that were given in [40,44,48], and the proposed technique for Example 6.1 Case 2 where l = 1 α = −0.5 and β = 0.5.  Table 5. Comparison E max between results in [48] and our results at α = −β = −1 2 for Example 6.1 Case 2 with disjoint choices of l and n.
Method in [48] Proposed method the comparison of E max between numerical results given in [48] and our method at α = −β = −1 2 for Example 6.1 Case 2 with different values of l and n. The values of E max are better than that given in [48], in addition to the extension the interval from [0, 1] to [0, 2] and [0, 4] which is not solved in [44]. Moreover, Figure 2 shows the analytical and the approximate solutions (left figure), the absolute error (right figure) for Example 6.1 Case 2 at n = 5, l = 6 and α = β = 0. These results in a whole compact with the numerical results in [50].

Remark 6.2:
At α = 1 2 , β = −1 2 and l = 1, 2 & 4, numerical results that obtained by our technique coincidence with the results that obtained by the authors in [49]. Moreover, we can obtain numerical results in Case 2 for several choices of the indexes (α, β) at disjoint values of n and the length of the interval [0, l]. These results approximately give the analytic solution. Moreover when the values of the fractional parametes are integeres or more than one, we must have this into our accounts in the calculting of the square A (became matrix of the ordinary differentiate in the integer case). Example 6.2: Assume the following equation: This example will take several forms depends on the known coefficients a, b and c, the initial conditions u 0 , u 1 , the known function g(t), and the variable of fractional order derivatives μ(t), γ (t) as the following cases: Case 1: [49] Suppose that The analytical solution of Example 6.2 Case 1 is u(t) = t 2 . Now, take only three terms of the numerical series of the solution and follow the same steps that followed    Table 6. c 0 , c 1 and c 2 of Example 6.2 Case 1 for several choices of α and β at n = 2, l = 1.
Which is in full accordance with the exact solution. Also, Figures 3 and 4 support the results (Table 6).

Remark 6.3:
If α = −β = 1 2 and l = 1, 2 & 4, the numerical results that obtained by our technique for Example 6.2 Case 1 is congruent to the results that obtained by the authors in [49].
Case 2: [20,49] Suppose that a = 1, b = −10, c = 1, u 0 = 5, The analytical solution of Example 6.2 Case 2 is u(t) = 5(1 + t) 2 . Again, we apply the proposed technique that proposed Section 4 and follows the same steps that followed in Example 6.1, we have the numerical results of this Example, Example 6.2, in this case (Case 2).
The absolute error that obtained via the suggested method for Example 6.2 Case 2 at α = β = 1 and n = 2 are given Table 7. By looking at Table 7, it transparent that the claimed results through the proposed technique are more accurate than others [20,49]. Moreover, the numerical results can be obtained for various values of α and β. Table 7. Comparison E max of the methods were given in [20,49] and the suggested method for Example 6.2 Case 2 at 1] In [20], n = 4 I n [ 49], n = 2 The The exact solution of Example 6.2 Case 3 is u(t) = t 2 .
In this case, the physical system whose behaviour is governed by an (ODE) given by: where, a symbolize the mass of the particle linked to the spiral, b is a measurement of the power of the damper, c refers to the spring stiffness, and g(t) is the applied outer force. t time and u(t) is the displacement of the mass from its rest position. Now, using technique given in Section 4 at α = β = 0 and n = 2. Then the exact solution of the damped mechanical oscillator Equation (58) that given in Case 3 of Example 6.2 is obtained. Also, we can claim the numerical results of Example 6.2 Case 3 for several values of the indexes α, β such as α = β = 1, α = β = 1 2 and α = −β = ∓ 1 2 at n = 2. These results are comparable with the exact solution.

Remark 6.4:
If α = −β = 1 2 and n = 2, l = 1. Then, the numerical results that obtained by our technique for Example 6.2 Case 3 is coincident with the results obtained by the authors in [49]. Therefore, the method that introduced in paper [49] considered as a special case of our suggested method in this article.
The exact solution of Example 6.2 case 4 is u(t) = t 2 .
The three unknown coefficients will be obtained in Table 8. These unknowns are calculated by applying the suggested method that proposed in Section 4 in the case of n = 2.
Consequently, we can rewrite the approximate solution as Which is the exact solution.
It is clear that with the aid of Table 8, the exact solution for Example 6.2 Case 4 is obtained at n = 2, l = 1 via the steps of a solution by the introduced method. Although, the best results of the same problem which obtained by the authors given in [21] are 6 × 10 −10 and 2 × 10 −10 . Moreover, the authors of the same refereed article obtained the exact solution when n → ∞. It is obvious that the suggested technique is more accurate than other methods given in [21]. Remark 6.5: If α = β = 0 and l = 1. Then, the approximate results that gained through the presented technique for Example 6.2 Case 4 is exactly like the results that obtained by the authors in [43] which considered as a special case of our suggested technique.

Conclusions
In this article, the operational matrix of shifted Jacobi polynomials is investigated to solve the multi-term FDEs and multi-term variable order FDEs numerically. The basic idea behind this method is converted the main problem into a system of algebraic equations that can be solved numerically. Also, the presented technique is applied for several choices of the indexes (α, β) and disjoint fractional variable-order fractional derivatives μ(t), γ (t) which implies us with the flexibility to solve many fractional order problems. Therefore, we enabled to solve many different applications in physics such as the damped mechanical oscillator and Bagley-Torvic equation. Moreover, the error estimate Table 8. c 0 , c 1 and c 2 of Example 6.2 Case 4 for several choices of α and β at n = 2, l = 1. has been discussed. Furthermore, numerical examples are introduced in many cases for several choices of the parameters of the original FDEs to demonstrate the accuracy, efficiency, generality, and applicability of the introduced technique. From the obtained results, only a few terms of the shifted Jacobi polynomials are required to acquire satisfactory results. Finally, in future work, we can use the proposed method for solving non-linear problems in one and two dimensional in addition to the linear problems in two dimensional. All computed numerical results obtained via the MATLAB program and +e data used to support the findings of this study are included within the article.