Numerical solution of the Bagley–Torvik equation using shifted Chebyshev operational matrix

In this study, an efficient numerical scheme based on shifted Chebyshev polynomials is established to obtain numerical solutions of the Bagley–Torvik equation. We first derive the shifted Chebyshev operational matrix of fractional derivative. Then, by the use of these operational matrices, we reduce the corresponding fractional order differential equation to a system of algebraic equations, which can be solved numerically by Newton’s method. Furthermore, the maximum absolute error is obtained through error analysis. Finally, numerical examples are presented to validate our theoretical analysis.


Introduction
Over the past few decades, many natural phenomena have been successfully modeled using the fractional differential equation [1][2][3][4][5][6]. As an example, the authors in [7] constructed the following equation that describes the motion of a rigid plate immersed in a Newtonian fluid: and where the operator D 3 2 is a Liouville-Caputo derivative, A = 0, B are constants, and the function f (t) is known. The existence and uniqueness of the problem have been established in [8,9]. Many methods have been developed to deal with the problem in the literature [10][11][12][13][14]. Podlubny also investigated this equation and introduced an approximate analytical solution by Green's function in his monograph [15]. Ray and Bera [16] adopted © The Author(s) 2020. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. a semi-analytical method for solving Bagley-Torvik equation and obtained the same solution as Podlubny's solution. Rajarama and Chakraverty [17] adopted the Sumudu transformation method to obtain the analytical solution of the problem. Cenesiz et al. [18][19][20] suggested a Taylor polynomial along with the collocation method for dealing with a class of fractional differential equations including the Bagley-Torvik equation. In [21][22][23][24], the wavelet method was used to deal with the problems. Diethelm and Ford [25] solved the problem by using Adams predictor and corrector methods. In [26][27][28] the spectral collocation method based on a hybrid function and Chebyshev polynomial were employed to handle the equation. Moreover, shifted Legendre polynomial based Galerkin and collocation methods were utilized in delay Bagley-Torvik equations in [29]. Most recently, Hou and Ji [30,31] introduced Jacobi polynomials and Laplace transform together with Laguerre polynomials to solve the equation.
Papers [26,27] and [32] focused on the Chebyshev polynomial method for the Bagley-Torvik equation. In these studies, the operational matrix of fractional integration and Tau method were applied to tackle the problem. The objective of the current study is to develop a modified Chebyshev spectral collocation method to handle the Bagley-Torvik equation. We generate the operational matrices of derivative for shifted Chebyshev polynomials in the physical space. Thereafter, we obtain a discrete numerical scheme, in which the nonhomogeneous terms are not approximated. A rigorous error analysis in L ∞ -norm is provided.

The fractional integration and differentiation
In this section, we mainly introduce the widely used Riemann-Liouville fractional integral and Liouville-Caputo fractional derivative.

Shifted Chebyshev polynomials and their properties
The well-known Chebyshev polynomials are defined on the interval [-1, 1] and are obtained by expanding the following formulae: To use these polynomials on the interval t ∈ [0, L] for any real L > 0, we introduce the change of variable x = 2t/L -1, 0 ≤ t ≤ L, and obtain the shifted Chebyshev polynomials T * Ln (t) = T n (2t/L -1).
The shifted Chebyshev polynomials T * Ln (t) satisfy the recurrence relation where T * L0 (t) = 1, T * L1 (t) = 2t/L-1. The analytic form of the shifted Chebyshev polynomials T * Li (t) of degree i is given by where T * where the interpolation points are chosen to be the Chebyshev-Gauss-Lobatto points associated with the interval [0, L], t k = L 2 (1 -cos(kπ/N)), k = 0, 1, 2, . . . , N . Here, the summation symbol with double primes denotes a sum with both the first and last terms halved.

The operational matrix of derivative
A continuous and bounded function y(t) can be approximated in terms of shifted Chebyshev polynomials in the interval [0, L] by the formula Using the discrete orthogonality relation, the coefficient c k in (6) is given by the explicit formula Applying (6), (7), the function y N (t) can be written collectively in a matrix form where The derivative y N (t) is as follows: We know that in which M is the (N + 1) × (N + 1) operational matrix of derivative given by , so that m 1 , m 2 , and m 3 are respectively N , 0, 2N for odd N and 0, 2N , 0 for even N . Then, we substitute equation (10) into (9) to get Therefore y N (t) can be expressed in a discretized form as follows: So, we can get the operational matrix of derivative Furthermore, the operational matrix of derivative of the n-order derivative can be completely determined from those of the first derivative

Calculation of the operational matrix of fractional order derivatives
According to the definition of Liouville-Caputo fractional derivative, we can write where α > 0. Applying (5), the Liouville-Caputo fractional derivative of the vector T(x) in (13) can be expressed as where Using (4) we have where , .
Employing (13), (14) and (15), we get Then we get the operational matrix of fractional derivative For simplicity, the operational matrix of the fractional derivative in (16) can be written collectively in the following form:

Applications to the Bagley-Torvik equation
To show the fundamental importance of the operational matrix of fractional order derivatives, we apply it for solving the Bagley-Torvik equation. To solve the problem, we first consider incorporating boundary conditions By substituting the approximation (18) in (19) and by using the boundary conditions (2), we get a system of algebraic equations: Solving the system of algebraic equations, we can obtain the vector Y . Then, using (8), we can get the output response

Some useful lemmas
In this section, we give some useful lemmas, which play a significant role in the convergence analysis later. We first introduce some notations that will be used. Let I := (-1, 1) and L 2 ω α,β (I) be the space of measurable functions whose square is Lebesgue integrable in I relative to the weight function ω α,β (x). The inner produce and norm of L 2 ω α,β (I) are defined by For a nonnegative integer m, define with the seminorm and the norm as follows: For a given positive integer N , we denote the points by {x i } N i=0 , which is a set of N + 1 Gauss-Lobatto points, corresponding to the weight ω(x). By P N we denote the space of all polynomials of degree not exceeding N . For all v ∈ C[-1, 1], we define the Lagrange interpolating polynomial I N v ∈ P N , satisfying The Lagrange interpolating polynomial can be written in the form where F i (x) is the Lagrange interpolation basis function associated with {x i } N i=0 .

Lemma 3 ([34]) Assume that v ∈ H m ω , and denote I N v its interpolation polynomial associated with the Gauss-Lobatto points {x
Then the following estimates hold:

Convergence analysis
In this section, an error estimate of the applied method for the solutions of the Bagley-Torvik equation is provided. For the sake of applying the theory of orthogonal polynomials, we use the variable transformations t = T(1 + x)/2, x ∈ [-1, 1] to rewrite (1), (2) as follows: and (22), which is assumed to be sufficiently smooth. Let the approximate solution u N (x) be obtained by using the proposed method. If u(x) ∈ H m ω (I), then for sufficiently large N the following error estimate holds:

Theorem 4 Let u(x) be the exact solution of the Bagley-Torvik equation differential equation
where F j , j = 0, 1, 2, . . . , N , is the Lagrange interpolation basis function. Consider equation (22). By using Chebyshev-Gauss-Lobatto collocation points {x i } N i=0 , we have u(x i ) = Then the numerical scheme (20) can be written as We now subtract (24) from (26) and subtract (25) from (27) to get the error equation Multiplying by F i (x) both sides of (28), (29) and summing from 0 to N yield Consequently, It follows from the Gronwall inequality and [35] that Using Lemma 3, we have We now estimate J 4 . By virtue of Lemma() with m = 2, we have Therefore, a combination of (30), (31), (32), and (33) yields estimate (23).

Illustrative examples
To illustrate the effectiveness of the proposed method in the present paper, some test examples are carried out in this section. The results obtained by the present methods reveal that the present method is very effective and convenient for fractional differential equations.
Example 9.1 As the first example, we consider the following Bagley-Torvik differential equation [36][37][38]: Example 9.2 In this example we consider the following equation: where The analytical solution can be found in [15]. The problem is considered in [18,20,22,28,39]. First, we consider the boundary conditions y(0) = 0, y(20) = -1.48433 and apply the present method to solve the problem with N = 8, 16, 32, 64, 128. In Table 1, we list the L ∞ , L 2 errors and CPU time for the differential values of N . The numerical solutions obtained by the present method and some other numerical methods, such as the wavelet method [22] and the hybrid functions method [28], are given in Tables 2 and 3. Clearly, numerical results show that the present method is working well and the accuracy is comparable     with the existing methods. Also, the numerical results with N = 16, 32 and the exact solution are plotted in Fig. 1. For Example 9.2, Fig. 1 shows that the approximate solutions using the present method are in high agreement with the exact solutions. Second, we solve this problem with the boundary conditions y(0) = 0, y(1) = 2.95179355. We compare the absolute errors of the present method, the Taylor collocation method [18], and the fractional Taylor method [20] in Table 4. This indicates that our results are better than given by [18,20].

Conclusion
In this work, the shifted Chebyshev operational matrix of fractional derivatives has been derived. Also, the operational matrix in combination with a collocation method is used to approximate the unknown function of the Bagley-Torvik equation. Moreover, a convergence analysis was performed under the L ∞ norm. Finally, numerical examples were presented to demonstrate the validity and applicability of the proposed numerical scheme. From examples, we observed that our scheme is simple and accurate. We believe that the ideas introduced in this study can be extended for systems of nonlinear fractional differential equations and fractional integro-differential equations.