A Galerkin-Type Fractional Approach for Solutions of Bagley-Torvik Equations

In this study, we present a numerical scheme similar to the Galerkin method in order to obtain numerical solutions of the Bagley Torvik equation of fractional order 3/2. In this approach, the approximate solution is assumed to have the form of a polynomial in the variable t = x, where α is a positive real parameter of our choice. The problem is firstly expressed in vectoral form via substituting the matrix counterparts of the terms present in the equation. After taking inner product of this vector with nonnegative integer powers of t up to a selected positive parameter N, a set of linear algebraic equations is obtained. After incorporation of the boundary conditions, the approximate solution of the problem is then computed from the solution of this linear system. The present method is illustrated with two examples.


Introduction
Historically, the idea of fractional calculus first appeared in a letter written by Leibniz to L'Hôspital in 1695. Curiously, the first rigorous definition of derivatives having noninteger orders only appeared in a work of Lacroix [Lacroix (1819)] in 1819. Over the years, the subject of fractional calculus attracted many mathematicians; as a result, different approaches were adopted to define fractional differential operators. Interested readers may refer to Ross [Ross (1977)] for an account on the history of fractional calculus.
Although the concept of fractional derivative has such a long history, for almost three centuries it remained as a topic which is only of theoretical interest for mathematicians. The realization that it can be used as a tool to explain physical phenomena took place as late as 1980s. Several of first such studies belong to Torvik (1983, 1984)], where they used fractional calculus to describe the behaviour of real materials. The equations that they proposed in order to simulate the motion of a rigid plate immersed in a Newtonain fluid contained derivative of order 3/2. Fractional derivatives have also been used in various fields such as control theory [Manabe (2002)], electrochemistry [Oldham (2010)], oil industry [Fitt, Goodwin, Ronaldson et al. (2009)] and vibrations [Hedrih (2006)].
In this study, our main interest will be to obtain approximate solutions to the Bagley-Torvik equation given by PðxÞy 00 ðxÞ þ RðxÞy ð3=2Þ ðxÞ þ SðxÞyðxÞ ¼ f ðxÞ; a x b; (1) under the boundary conditions Here y(x) is the unknown function to be determined and y (3/2) (x) is the fractional derivative of order 3/2 which will be defined in the next section, P(x), R(x), S(x) and f(x) are real-valued functions defined on a≤x≤b, and a, b, c 0 , c 1 are real numbers. Here, we note that it is far more common in the literature to consider Eq. (1) subject to initial conditions rather than boundary conditions as given in Eq.
(2). The existence and uniqueness of the solution to problem Eqs. (1)-(2) is proved in Al-Mdallal et al. [Al-Mdallal, Syam and Anwar (2010)] in the special case that P(x), R(x) and Q(x) are constants. We refer the reader to Staněk [Staněk (2013)] for a treatment of existence and uniqueness of solutions of nonlinear fractional differential equations of Bagley-Torvik type.
A few words on the physical significance of Eq. (1) might be helpful for the reader. The relationship between the stress field and the transverse fluid velocity field contains a fractional derivative; namely, the stress field is proportional to the fractional derivative of order 1/2 of the transverse fluid velocity field. In view of this phenomenon, if a rigid plate is immersed in a Newtonian fluid and is applied an external force f(x), the displacement of this plate is known to satisfy the Bagley-Torvik Eq. (1). The reader is referred to Esmaeili [Esmaeili (2017)] for a more thorough explanation.
Since the Bagley-Torvik Eq. (1) is of great importance, a substantial amount of literature has been devoted to the examination of its various aspects. Among these studies, a large proportion is related to obtaining its numerical solutions. To name a few of such studies, Yüzbaşı [Yüzbaşı (2013) [Ji and Hou (2020)] and a method based on fractional Taylor vector approximation [Krishnasamy and Razzaghi (2016)]. Lastly, Esmaeili [Esmaeili (2017)] used exponential integrators to solve the Bagley-Torvik equation by firstly converting it to an equation of order 1/2.
In this paper, we present a Galerkin-like approach to solve the Bagley-Torvik Eq. (1) under the boundary conditions Eq.
(2) instead of the more commonly used initial conditions. We claim that the presented method gives fairly accurate results with relatively small computational cost.
The organization of the paper is as follows: In Section 2, some preliminary information is given. The Galerkin-like method is explained in Section 3. Then, we discuss two numerical examples in Section 4. Finally, the conclusions regarding the present scheme are given in Section 5.

Basic definitions
In this section, we define the Caputo derivative for the fractional derivative present in Eq.
Definition 2.1 A real-valued function f(x) defined for x>0 is said to belong to the space C μ , where l 2 R, if there exist a real number p>μ and a function f 1 (x) ∈ C[0,∞] such that f(x)= x p f 1 (x).
A direct consequence of this definition is that C l &C c for γ<μ.
Definition 2.2 A real-valued function f(x) defined for x>0 is said to belong to the space C m l , where m 2 N [ f0g; l 2 R, if f (m) (x) ∈ C μ holds for the m-th derivative of f(x). Now, we are ready for the definition of fractional derivative. Definition 2.3 Let f ∈ C μ , where μ≥−1. The Riemann-Liouville fractional integral operator of order α≥0 of f, denoted by J α f(x) is defined by the following: The Γ that appears in the above definition is the special function defined by ÀðzÞ ¼ R 1 0 x zÀ1 e Àx dx for every complex number z which is not a nonpositive integer. The fractional derivative in the Caputo sense is advantageous over the Riemann-Liouville sense in that it is easier to deal with the integer order initial conditions using fractional derivative in the Caputo sense [Staněk (2013)]. Therefore, the following definition is in order: The following consequences for the Caputo fractional derivative will be important for us in the remaining of the paper: Here, bac stands for the largest integer that is not greater than α. It is important to notice that the above does not define D a Ã x b for 0<β<1. Therefore, for these values of β, we use the identity which holds for the fractional derivative of order α of x β in the Riemann-Liouville sense. Note that this fractional derivative does not exist if b ¼ 1 2 : The proofs of these properties are straightforward using Definition 2.4. The interested reader can find them among other useful properties in a study by Diethelm et al. [Diethelm, Ford, Freed et al. (2005); Podlubny (1998)].

Method of solution
In this section, we will describe the procedure to solve Eq. (1). The same method was employed to obtain approximate solutions of high-order Fredholm integro-differential equations [Türkyılmazoğlu (2014)] and high-order integro differential equations with weakly singular kernel [Yüzbaşı and Karaçayır (2016)].
As the first step of the Galerkin-like scheme, we assume that the unique solution y(x) of Eq. (1) is uniquely expressible in the form of a power series where t=x α for some α>0 such that M a 6 ¼ 1 2 for any integer M (see the last paragraph of Section 2). Truncating this power series after the (N+1)st term then yields Here, the variable row vector X(x) and the column vector A are given by Under this setting, the coefficients a i , i=0, 1, … , N are the unknown constants which will be determined as the output of the method. The approximate solution y N,α , which is a polynomial of degree N of the variable t=x α , can then be obtained from these coefficients.
Since the solution method should be programmable for computer, we would like all the operations to be expressed in terms of matrices. To this end, the second ordinary derivative of y N,α (x) can be expressed as a product of matrices with the help of a special matrix. Namely, if we define B to be the (N+1) × (N+1) matrix such that B k+1, k+1 = kα(kα−1) for k=1,2, …, N and B i,j = 0 elsewhere, then the following equality holds: Here, multiplication by the term x −2 from the left-hand side is to be interpreted as scalar multiplication.
As for the fractional derivative of order 3/2 of the approximate solution y N,α (x), it is useful to rewrite property Eq.
(3) of the Caputo derivative for a ¼ 3 2 as follows: In order to deal with the task of expressing the fractional derivative of y N,α by means of a product of matrices, we define a new auxiliary matrix as follows: Let Γ (3/2),α be the (N+1) × (N+1) diagonal square matrix defined by À Note again that all the entries of Γ (3/2),α that are not on the main diagonal are equal to 0. As an illustration, for the choice of α=1, the matrix Γ (3/2),α is given by N ;a ðxÞ ¼x À3=2 Xðx a ÞΓ ð3=2Þ;a A: The next step is to substitute the matrix expression Eq. (4) for y N,α (x), the matrix expression Eq. (5) for the ordinary second derivative y N,α ′′ (x) and the expression Eq. (6) for the fractional derivative y N ;a ðxÞ into Eq. (1). This gives us the relation where GðxÞ ¼ PðxÞx À2 Xðx a ÞB þ RðxÞx À3=2 Xðx a ÞÀ ð3=2Þ;a þ SðxÞXðx a Þ: Now, it is time to apply the central idea of the present numerical method. Namely, we now apply inner product to Eq. (7) with the elements of the set Ф = {1, x α , x 2α , …, x Nα }. The inner product to be used here is the standard inner product in the Hilbert space L 2 [a, b], which is defined by where f and g are two functions from L 2 [a,b]. For each k=0, 1, … , N, inner product of Eq. ; Before proceeding with the solution of the linear system WA=F, a restriction on the approximate solution y N,α (x) is in order. We would like the error of the approximate solution to be equal to zero on the boundary points; in other words, we demand y N,α to satisfy the boundary conditions Eq. (2) given by y(a)=c 0 and y(b)=c 1 . This restriction we impose on y N,α implies the linear equations P N k¼0 a k a ka ¼ c 0 and P N k¼0 a k b ka ¼ c 1 . In order to form a new linear system including these two equations in a 0 , a 1 , …, a N , we express the initial conditions in vector form. They can be written as With the aim of including these boundary conditions in the algorithm, we sacrifice the equations corresponding to inner product with x 0 =1 and x 1 in favour of Eq. (8) corresponding to the boundary conditions Eq.
(2). This amounts to updating the first two rows of the system matrix W and the first two entries of the right-hand size F, thus yielding a new systemWA ¼F, explicitly given bỹ ; Finally, provided that the modified system matrixW is of full rank, we compute the matrix of unknown coefficients as A ¼W À1F , and thus the approximate solution is obtained by Before moving on to the next section, it will be of benefit to summarize the scheme explained in this section in a step-by-step fashion. Such a description can be as follows: GðxÞ 1;j dx STEP 8: DEFINE THE COLUMN MATRIX F AS FOLLOWS: f ðxÞdx STEP 9: SOLVE THE LINEAR ALGEBRAIC SYSTEM WA = F STEP 10: SET y N,α (x) = a 0 + a 1 x α + a 2 x 2α + … + x Nα END

Numerical examples
In this section, we solve two example problems using the method explained in Section 3. y 00 ðxÞ þ y ð3=2Þ ðxÞ þ yðxÞ ¼ 1 þ x; 0 x 1; yð0Þ ¼ 1; yð1Þ ¼ 2: The exact solution of this problem is y exact (x)=1+x. Using the method explained in Section 3, we obtained approximate solutions of Eq. (9) corresponding to the values 1 and 1/3 for the parameter α and several values for the parameter N. For instance, applying the method with N=4 results in the linear system given by Solving this system yields the unknown coefficients given by a 0 ¼ 1; a 1 ¼ 1; a 2 ¼ À2:051Â10 À15 ; a 3 ¼ À9:335Â10 À15 ; a 4 ¼ À2:907Â10 À15 : Thus, the approximate solution y 4,1 is obtained by y 4;1 ðxÞ ¼1þxÀ2:051Â10 À15 x 2 À9:335Â10 À15 x 3 À2:907Â10 À15 x 4 : In a similar manner, we have carried out the calculations required to compute the approximate solution y 4,1/3 and obtained y 4;1=3 ðxÞ ¼1À1:1102Â10 À15 x 2=3 þxÀ1:0111Â10 À15 x 4=3 : We obtained the approximate solutions corresponding to other values of N for both values of the parameter α. In order to measure their accuracy, we consider their actual absolute error given by e N,α (x)=|y N,α (x) − y exact (x)|. Fig. 1 depicts the absolute actual errors of the approximate solutions obtained by N=4,5,6,9 corresponding to the parameter value α=1. It is seen from the plot that increasing N significantly improves the accuracy of the approximate solutions. In addition, the approximate solutions corresponding to N=6 and N=9 are illustrated together with the exact solution for both values of the parameter α in Figs. 2 and 3. The approximate solutions seem indistinguishable from the exact solution.
Tabs. 1 and 2 give a more detailed comparison of the actual and approximate solutions for selected values of x.
Example 2. Next, let us consider the following constant coefficient Bagley-Torvik equation without ordinary derivative studied in Diethelm et al. [Diethelm, Ford, Freed et al. (2005); Esmaeili and Shamsi (2011)]: The exact solution of this problem is y exact (x) = x 2 − x. As in the previous example, we obtained approximate solutions which are polynomials of x and ffiffi ffi x 3 p ; in other words, we used the parameter values α=1 and α=1/3. Setting N=5 and carrying out the calculations explained in Section 3 gives rise to the linear algebraic system   The actual absolute errors corresponding to several N values are visualized in Fig. 4 for α=1 and in Fig. 5   Let us now compare the solutions we have obtained with Adams-Bashforth-Moulton method (ABM) [Diethelm, Ford, Freed et al. (2005)] and the pseudospectral method (PM) [Esmaeili and Shamsi (2011)]. As only the maximum absolute error over the entire interval [0,1] is considered in these two studies, we calculated the maximum absolute error e N ;a 1 of our approximate solutions. The results are collected in Tab. 4. The values show that the present scheme outperforms the other two methods.

Conclusion
In this paper, we have presented a Galerkin-like approach for the approximate solution of the fractional Bagley-Torvik equation. It turns out that our approach gives fairly good results besides its simplicity. Another advantage of this approach is that the accuracy of the approximate solutions undergo a significant improvement with increasing N values. Simulation results on two example problems show that high levels of accuracy can be attained even for small values of N. Comparison of one example problem with other methods has revealed that the present scheme gives more accurate results for similar parameter values. On the whole, the results make it clear that the numerical scheme presented in this paper can be relied upon when one would like to solve fractional-order differential equations of Bagley-Torvik type.
Acknowledgement: The authors wish to express their appreciation to the reviewers for their helpful suggestions which greatly improved the presentation of this paper.