An exponential spline approximation for fractional Bagley–Torvik equation

In this paper, we approximate the solution of fractional Bagley–Torvik equation by using the exponential spline function and the shifted Grünwald difference operator. The proposed methods reduce to the system of algebraic equations. The convergence analysis of the methods has been discussed. The numerical examples are presented to illustrate the applications of the methods and to compare the computed results with the other methods.


Introduction
Fractional calculus is an old topic in mathematical analysis, which goes back to Leibniz (1695) and Euler (1730) (see [15,16]). In recent years, the numerical solution of fractional equations has become a popular topic in applied sciences control and engineering. Bagley-Torvik equation appears in the modeling of the motion of a rigid plate submerged in a Newtonian fluid [6]. Existence and uniqueness theorem for Bagley-Torvik equation with Dirichlet boundary condition is given in [5]. In this article, the exponential spline will be employed to obtain the approximate solution of Bagley-Torvik equation with Caputo derivative subject to boundary conditions Here, D α is the Caputo derivative, f (x) is a continuous function, ω i (i = 1, 2), η, μ are real constants, and m = 1 or 2. In general, it is difficult to solve most of the fractional differential equations analytically. Therefore, numerical methods to find an approximate solution and qualitative behaviors of the solution for fractional differential equation have been investigated by authors in [1-9, 11, 14, 19-25, 27-29], and some references therein. The reproducing kernel method is employed for the fractional order differential equations in [1][2][3]. In [21], numerical solution of boundary value problem of fractional Bagley-Torvik equation is given in the reproducing kernel space. In [14], the authors study the numerical approach based on operational matrices of fractional differential equations with a hybrid of block-pulse functions and Chebyshev polynomials. The existence of positive and negative solutions and properties of their derivatives for the generalized Bagley-Torvik fractional differential equation is given in [24]. Numerical solution of the fractional Bagley-Torvik equation arising in fluid mechanics based on Taylor matrix method is given in [11].
In [29] the numerical solutions for fractional boundary value problem have been found by cubic spline polynomials. The numerical scheme for solving two-point fractional Bagley-Torvik equation using the Chebyshev collocation method has been solved in [22]. In [23] the Bagley-Torvik equation as a prototype fractional differential equation with two derivatives is investigated by means of homotopy perturbation method. The numerical solution to the Bagley-Torvik equation by exponential integrators is discussed in [9]. Also Adomian decomposition method for solving the initial value problem of Bagley-Torvik equation is discussed in [19], fractional linear multistep method and a predictor-corrector method of Adams type based on finite difference methods for initial value problem of Bagley-Torvik equation are discussed in [8]; Legendre operational matrix method for fractional differential equation is applied in [20]; and a combination of collocation points and first-order Bessel functions, which is called Bessel-collocation method for boundary value problem of Bagley-Torvik equation, is discussed in [27]. Quadratic spline solution for boundary value problem of fractional order is applied in [28], and an exponential spline technique for solving fractional boundary value problem is employed in [4]. The first aim of the present work is to explore exponential spline interpolation with multiple parameters and to produce the error of approximate exponential spline. The second aim is to introduce a new approximate technique to find solutions of fractional boundary value problem, and we demonstrate the convergence analysis for this technique.
In [10], the authors tried to approximate the solution of nonlinear fractional differential pantograph equations by sinc interpolation. At first, they have transformed the problem into a nonlinear integral equation with some delay terms and the kernel of this integral equation is weakly singular for the case 0 < α < 1, thus the solution is weakly singular and the numerical methods cannot achieve high accuracy in approximating solutions. The main advantage of our algorithm is that it can be used directly without using assumption or transformation formulae. This paper is organized into four sections. In Sect. 2, we describe basic definitions and the nonpolynomial spline method to approximate the solutions of fractional Bagley-Torvik equation. Convergence analysis is proved in Sect. 3. In Sect. 4, the numerical examples are given to illustrate the applications of the method, and also the computed results are compared with another known method in [4,9,17,21,28,29].

Basic definitions and description of the methods
In this section, we recall some definitions and properties of the fractional calculus theory, which are used in this paper. There are several definitions of a fractional derivative of order α > 0, such as Riemann-Liouville, Grunwald-Letnikov, and Caputo. In the present work, Caputo and Grunwald-Letnikov fractional derivatives are used for the formulation of the problem.

Definition 1
Let u(x) be a function defined on (a, b), then the Riemann-Liouville fractional derivative is of the following form [5]: where Γ is the gamma function.

Definition 3
Let u(x) be a function defined on (a, b), then the Caputo fractional derivative is of the following form [5]: ( 3 )
Let us consider a mesh with nodal points x i on [a, b] such that where h = b-a n , x i = a + ih for i = 0(1)n. Let u(x) be the exact solution of (1) and S i be an approximation to u i = u(x i ) obtained by the exponential spline function Q i (x) ∈ C ∞ [a, b] passing through the points (x i , S i ) and (x i+1 , S i+1 ). Then in each subinterval the parametric spline segment Q i (x) has the following form (see [12,18,26]): where β is a free parameter of the spline functions which can be real or pure imaginary and which will be used to raise the accuracy of the method, see [26].
For the development of consistency relations between the exponential spline approximation and its derivatives at the nodal points, we consider the following four rela- After a simple calculation, we obtain the values of coefficients, and using the second-order derivative continuity at the knots x i , for i = 1, . . . , n -1, we get where In the limiting case, when θ → 0, relations (9) and (11) reduce into the ordinary cubic spline relation: The proposed differential Eq. (1) in the mesh point (x i ) may be discretized by

Lemma 1
The local truncation error x i associated with equations (9) and (11) for i = 1, . . . , n -1 in the limiting case when θ → 0 is given by Proof The above expressions can be obtained by expanding the terms M i+1 , M i-1 , m i+1 , m i-1 , u i+1 , and u i-1 about the points x i in relations (12) using Taylor series respectively. Moreover, In a similar manner, we get

Cubic and exponential splines method for approximate fractional Bagley-Torvik equation
In this section, we give some methods to approximate D α (u(x))| x=x i by using spline function.
Method I. The discrete approximation of the Caputo fractional derivative D α (u(x)) can be obtained by a cubic spline (in the limiting case when θ → 0, the relations of exponential spline reduce into ordinary cubic spline relation) formula as follows (see [12] and [13]): Using the Caputo fractional derivative for 1 < α < 2, we get Using a piecewise technique, the following equation is obtained using equation (21) and Lemma 1: Since (x iη) 1-α does not change sign on [(j -1)h, jh], by the weighted mean value theorem for integrals and by applying to each integration of the last summation, we get After simple calculations, equations (6), (10), and (22) become Using [12] and [13], we have Qu ∞ = O(h 2 ). Also we obtain the following relation for i = 1, 2, . . . , n -1: where a in = c in = b i1 = d i1 = 0 for i = 1, 2, . . . , n. Also we approximate u i byû i and m i bym i such thatû i andm i for i = 1, 2, . . . , n satisfy system (11). We get Finally, we approximate the exact solution u i by the natural cubic spline function Q i (x) for i = 1, 2, . . . , n. In the matrix notation, we get where A = (a ij ), B = (b ij ), C = (c ij ), and D = (d ij ). Now, the valuesM i andm i are determined as the solutions of linear systems (12)(I) and (12)(II). We approximatem 0 = -3û 0 +4û 1 -û 2 2h , m n = 3û n-2 -4û n-1 -û n 2h , and also, by using boundary conditions, we approximateM i for i = 0, n. We need the following lemma.

Lemma 2
The matrices W and Z are obtained with the help of systems (9) and (11) invertible.
Proof The valuesM i are determined as the solutions of the following linear system: Also the valuesm i are determined as follows: Also, for determining the valuesM i andm i , in the limiting case when θ → 0, by using relations (12), the matrices W and Z are strictly diagonally-dominant matrices, then the matrices W and Z are invertible. Hence, Therefore, from (26) and (29) we obtain Method II. Suppose that M j = Q (x j ) is approximated u (x j ) in the subintervals [(j -1)h, jh] for i = 1, 2, 3, . . . , n -1 and j = 1, 2, . . . , i. Also, by using equation (23), the following recurrence relation is obtained: Bagley-Torvik equation (1)-(2) for 1 < α < 2 can be discretized as follows: Finally, we approximate the exact solution u i by the natural cubic spline function Q i (x) for i = 1, 2, . . . , n. In the matrix notation, we get where ρ = 1 (29) and (32) we also obtain Method III. In this section, we approximate the exact solution by use of the Caputo fractional derivative for 0 < α < 1 as follows: Using a piecewise technique, equation (34) becomes Since (x iη) -α does not change sign on [(j -1)h, jh], by the weighted mean value theorem for integrals and by applying to each integration of the last summation, we get After simple calculations, from equations (34) and (36) we get where a i k for i = 1, 2, 3, 4 are given in relation (8).

The weighted and shifted Grünwald difference operator and exponential spline function
Method IV. In this section, we would like to develop a numerical method based on the methods in references [4,25,29], and [28]. Also we investigate the convergence analysis of this method. Let U = (u i ), S = (s i ), C = (c i ), T = (t i ), and E = (e i ) = U -S = U -Q i (x) be (n -1)-dimensional column vectors. We used the weighted and shifted Grünwald difference operator and exponential spline function. By using consistency relation (12)(I) and using the boundary condition, we get the system of algebraic equations |i -j| = 1, (45) Page 12 of 20 We assume that F = (f 1 , f 2 , f 3 , . . . , f n-1 ) t , S = (s 1 , s 2 , s 3 , . . . , s n- where 3 . . . Substituting equation (46) into equation (43), we obtain and Now the error equation can be written as follows: In order to get a bound on E (the infinite norm), we need the following lemma.

Convergence analysis
In this section, we discuss the convergence analysis of exponential spline Method III. Convergence analyses of Method I and Method II are similar. So, first we write equation (41) in the points of x i , i = 1, 2, . . . , n -1: Now, using the results obtained in (18), we have Equations (60) construct a nonlinear system, which can be solved by Newton's iterations method. Let u(x) be the exact solution of the problem and Q(x) ∈ C ∞ [0, T] be the exponential spline approximation to u(x) satisfied in Q(x i ) = u(x i ), i = 1, 2, . . . , n -1, and Q (x i ) = u (x i ), i = 0, n. We should approximate the error u(x) -Q(x) . Let us assume thatQ(x) is the computed spline approximation to Q(x). To estimate u(x) -Q(x) , we will estimate u(x) -Q(x) and Q (x) -Q(x) separately.

Lemma 4 LetQ(x) be the unique spline interpolation to Q(x)
, and also suppose that partial derivatives of F exist and | ∂F ∂u | ≤ k 1 , | ∂F ∂u | ≤ k 2 for some constants k 1 and k 2 . Then, for 0 ≤ i ≤ n, we have Proof For 1 ≤ i ≤ n -1, we get Now, using the mean value theorem for two parts of the above relation, there exist ξ i and ν i such that Using relation (57), we have |Q( , and taking the absolute value, we obtain (62) T] be the exact solution (1) and Q(x) be the exponential spline approximation to u(x), then we have Proof SinceQ(x) is an interpolation to u(x), thus there is a finite constant 1 independent of h that we get where 1 is a finite constant. Now, using the triangular inequality and Lemma 1, we can obtain the results as follows: We can prove the convergence analysis for Method I and Method II in the same manner.

Numerical results
In this section, we have implemented our methods for solving some of the Bagley-Torvik differential equations with different values of h = 1 8 ,  [4,9,17,21,28,29]. The convergence order (C.O.) is obtained by where E(h) is the maximum absolute error. Numerical results can be derived by using MATHEMATICA 9.
Example 2 Consider the following boundary value problem: where the exact solution is given by the relation u(x) = x 6 (1x 2 ). The maximum absolute errors of Method III and Method IV are presented in Tables 7 and 4. Also compare the computed results with the methods [28] and [4] in Tables 5 and 6.

Conclusion
Computational methods for solving the fractional Bagley-Torvik equation were proposed. The fractional differential equation term in the fractional Bagley-Torvik equation was discretized using the exponential spline function and the shifted Grünwald difference operator. Also we obtain the four numerical schemes based on the exponential spline. The convergence analyses of the shifted Grünwald difference and the exponential spline are discussed. The feasibility of the numerical algorithms was illustrated with four examples, and the approximated results were compared with the methods in [4,9,17,21,28,29].