Quadratic B-Spline Galerkin Method for Numerical Solutions of Fractional Telegraph Equations

Letters which were written by two famous mathematicians G.W. Leibnitz and L’Hospital to each other in 1695 can be considered as the beginning of the fractional calculus’ taking part in literature. After that time, a huge contribution has been made for the development of arbitrary order differentiation and integration by a lot of famous mathematicians such as Euler, Laplace, Fourier, Lacroix, Abel, Riemann, Liouville, Caputo [7]. Fractional order differentiation concept helps to state the many physical problems and lead to huge amount of applications [8]. Recently scientists have shown the effectiveness of the fractional order differential equations on expressing the complex events and modelling many physical, engineering phenomenons in their studies [4]. The huge amount of applications on these type equations can bee seen frequently in various fields such as viscoelastic, biology, signal process, electromagnetic, chaos and fractals, traffic system, chemistry, control system, economics, finance and etc. [9]. Riemann-Liouville fractional concept which proposed by two famous mathematicians Riemann and Liouville is composed of fractional integral and fractional derivative. In this concept, in the fractional initial value problems, the initial conditions are comprised of limit values of Rieman-Liouville derivative in the initial point. This situation causes a big problem in solving fractional initial value problems. There isn’t any physical meaning of these initial points. Caputo’s fractional derivative which is presented byM. Caputo in 1967 involves the limit values in the initial points of integer order derivatives with initial values that are given with fractional order equation. Due to this user friendliness of Caputo’s definition, many scientists prefer Caputo’s derivative as fractional derivative operator in many in their studies [4]. In last decades, the importance of fractional order differential equations increased rapidly. Due to this increasement investigating the analytical and numerical solutions of fractional order differential equations(FDEs) becomes an attractive topic for scientists. So a lot of methods are used to obtain analytical and numerical solutions of FDEs such as Laplace transform method [4], power series method [4], finite element method [10, 11], Adomian decomposition method [12], variational iteration method [13], differential transform method [14], homotopy perturbation method [15], homotopy analysis method [16, 17], finite difference methods [18], and etc. Finite elements method which is used commonly in various fields of physics and engineering, first arise in 1960. Argyris, Clough and Zienkiewicz made contribution to this method [19]. After the development of computer technology in recent 50 years, this method has been occupied in an important place in solving many problems which arise in physics and engineering [20]. Bulletin of Mathematical Sciences and Applications Submitted: 2016-09-21 ISSN: 2278-9634, Vol. 18, pp 23-39 Revised: 2016-12-09 doi:10.18052/www.scipress.com/BMSA.18.23 Accepted: 2016-12-09 2017 SciPress Ltd, Switzerland Online: 2017-05-31


Introduction
Letters which were written by two famous mathematicians G.W. Leibnitz and L'Hospital to each other in 1695 can be considered as the beginning of the fractional calculus' taking part in literature. After that time, a huge contribution has been made for the development of arbitrary order differentiation and integration by a lot of famous mathematicians such as Euler, Laplace, Fourier, Lacroix, Abel, Riemann, Liouville, Caputo [7]. Fractional order differentiation concept helps to state the many physical problems and lead to huge amount of applications [8]. Recently scientists have shown the effectiveness of the fractional order differential equations on expressing the complex events and modelling many physical, engineering phenomenons in their studies [4]. The huge amount of applications on these type equations can bee seen frequently in various fields such as viscoelastic, biology, signal process, electromagnetic, chaos and fractals, traffic system, chemistry, control system, economics, finance and etc. [9].
Riemann-Liouville fractional concept which proposed by two famous mathematicians Riemann and Liouville is composed of fractional integral and fractional derivative. In this concept, in the fractional initial value problems, the initial conditions are comprised of limit values of Rieman-Liouville derivative in the initial point. This situation causes a big problem in solving fractional initial value problems. There isn't any physical meaning of these initial points. Caputo's fractional derivative which is presented by M. Caputo in 1967 involves the limit values in the initial points of integer order derivatives with initial values that are given with fractional order equation. Due to this user friendliness of Caputo's definition, many scientists prefer Caputo's derivative as fractional derivative operator in many in their studies [4].
In last decades, the importance of fractional order differential equations increased rapidly. Due to this increasement investigating the analytical and numerical solutions of fractional order differential equations(FDEs) becomes an attractive topic for scientists. So a lot of methods are used to obtain analytical and numerical solutions of FDEs such as Laplace transform method [4], power series method [4], finite element method [10,11], Adomian decomposition method [12], variational iteration method [13], differential transform method [14], homotopy perturbation method [15], homotopy analysis method [16,17], finite difference methods [18], and etc.
Finite elements method which is used commonly in various fields of physics and engineering, first arise in 1960. Argyris, Clough and Zienkiewicz made contribution to this method [19]. After the development of computer technology in recent 50 years, this method has been occupied in an important place in solving many problems which arise in physics and engineering [20].
In this study, let us consider the time fractional telegraph equations [1,2] as follows where λ, s 1 , s 2 and s 3 are constants and denote Caputo fractional derivatives [3,4]. Considered fractional order telegraph Eqs.(1)-(2) are solved by quadratic B-spline Galerkin method with the initial and boundary conditions. The fractional order derivatives are taken into account in the Caputo's sense.
where b γ k = (k + 1) 1−γ − k 1−γ and for 1 < γ ≤ 2, L2 approximation [6] ∂ γ f (t) where b γ k = (k + 1) 2−γ − k 2−γ are subrogated into the equations instead of the fractional derivatives. Now, let us to define quadratic B-spline base functions. Partitioning interval [a, b] into N finite elements of uniformly equal length by the knots x m , m = 0, 1, 2, ..., N such that a = x 0 < x 1 · · · < x N = b and h = x m+1 − x m . The quadratic B-splines Q m (x) at the knots x m are defined over the interval [a, b] by [21] Q m (x) = 1 where m = −1(1)N . The set of splines {Q −1 (x), Q 0 (x), . . . , Q N (x)} forms a basis for the functions defined over [a, b]. Therefore, we can write an approximation solution U N (x, t) with regard to the quadratic B-splines trial functions as: where δ m (t)'s are unknown, time dependent parameters to be determined from the boundary and weighted residual conditions. For these problems, we define the finite elements with the interval

Quadratic B-Spline Finite Element Galerkin Solutions
We take all the terms in Eq. (1) to the right hand side of the equation and then multiply by the weight function Ψ(x) to apply the Galerkin method to equation with the appropriate boundary conditions. Then, integrating the resulting equation over the region [0, 1] and setting it to zero, we get where Ψ(x) is the weighted function that taken as quadratic B-spline functions. As Eq. (9) is valid on the whole region, specially it is valid over the typical element [x m , x m+1 ] as follows When the differentiation is scattered between the U and Ψ, it results in an integral form which needs weaker continuity conditions on trial base functions and we have The transformation ξ = x − x m is used for changing the global coordinate system into the local one. Thus, Eq. (11) turns into the form The Eq. (12) is the element equation for a typical element "e". Eq. (5) turns into the form

Bulletin of Mathematical Sciences and Applications Vol. 18 25
Subrogating the Eqs. (13) into Eq. (12), we handle whose matrix form is where¨states γ th fractional derivative and˙states (γ − 1) th fractional derivative with respect to time and A e ij , B e ikj , C e ij , D e ij are element matrices given as follows: The element matrices are evaluated as Unifying the accretion values coming from all the elements, Eq. (15) produces the system where δ(t)'s are unknown parameters and A, B, C and D are (N + 2) × (N + 2) global matrices with generalized m th row. If time parameters δ(t)'s and its fractional time derivativesδ(t)'s andδ(t)'s in Eq. (16) are discretized with the aid of the Crank-Nicolson, L1 and L2 formulaes, respectively:

BMSA Volume 18
we have recurrence relationship between successive time levels relating unknown parameters δ n+1 Applying the same procedure to Eq.
(2), we get the matrix equation where The element matrices are evaluated as Unifying the accretion values coming from all the elements, Eq. (21) produces the system where δ(t)'s are unknown parameters and A, B, C and D are (N + 2) × (N + 2) global matrices with generalized m th row. Subrogating the Eqs. (17)- (19) into Eq. (22), we handle a recurrence relationship between successive time levels relating unknown parameters δ n+1 The unknown values of δ −1 arise for n = 0 and k = n in the systems (20) and (23). By using central difference approximation, we can write instead of the derivative term in the second initial value U t (x, 0) = g 2 (x). Hence we can take δ −1 as δ 1 − 2∆tg 2 (x). By using the initial values and vanishing the terms δ −1 and δ N from the systems (20) and (23) which are (N + 2) × (N + 2) type we acquire N × N type algebraic equation systems.

Initial state
The initial vector d 0 = (δ −1 , δ 0 , δ 1 , . . . , δ N −2 , δ N −1 , δ N ) T is determined from the initial and boundary conditions. So we can rewrite the approximation (6) for the initial condition as where the δ m (0)'s are unknown parameters. It is necessary to have the initial numerical approximation U N (x, 0) performs the following conditions: . Thus, using the these conditions cause to a two-diagonal system of matrix of the form

Numerical Examples and Results
In this section the numerical solutions of discussed problems are evaluated by using B-spline Galerkin Method. Problem 1: Consider the time fractional telegraph equation (1) with s 1 = 1, s 2 = 1 and s 3 = π with boundary conditions U (0, t) = 0 , U (1, t) = t 3 sin 2 (1), t ≥ 0 and the initial conditions as and the exact solution of this problem is given by [2] U (x, t) = t 3 sin 2 (x).
In Table 1 numerical solutions, L 2 and L ∞ are given for Problem 1. It is seen from Table 1 that when partition number N increases for results that evaluated to different values of N , ∆t = 0.001, t = 1 and γ = 1.50, L 2 and L ∞ error norms decrease as it expected. In addition to this from Table 2, it is seen that L 2 and L ∞ error norms decrease when ∆t time step decreases for results that evaluated to different values of ∆t, N = 30, t = 1 and γ = 1.50 as it is expected. The comparisons of exact and numerical solutions to Problem 1 for ∆t = 0.001, t = 1, N = 30 and different values of γ are given in Table 3. It is clearly seen that the results are better for the smaller values of γ. Absolute error graphics are given in Figure 1 for numerical results which are obtained by Galerkin method in the values ∆t = 0.001, N = 30, t = 1 and different values of γ. When the absolute error graphics are compared with each other, it is understood that the numerical results which are obtained by using quadratic B-spline Galerkin method are better then the results which are obtained by the aid of radial basis functions taken from [2]. Table 1: Error norms and numerical solutions of Problem 1 for γ = 1.50, ∆t = 0.001, t = 1.
The exact solution of the problem is given by [1] U (x, t) = t 3 sin(2πx).
The results for Problem 2 are presented in Table 4-6. In Table 4 the exact solutions and the numerical solutions for ∆t = 0.0005, t = 1, γ = 1.50 and different values of N are given. It is clearly seen from the decreasement in the L 2 and L ∞ error norms, the numerical solutions are closer to the analytical solution while the partition number N increases. Also, from Table 5 it is deduced that L 2 and L ∞ error norms decreases when the ∆t time step decreased for the N = 100, t = 1, γ = 1.50. The error norms and the numerical results that obtained for different values of γ are given in Table  6. The error norms of the exact solutions and the numerical solutions that obtained for ∆t = 0.0005, N = 100, t = 1, ν = 1 and different values of γ by using Galerkin method are given in Figure 2. The comparison of the error norms that obtained by the implemented method for γ = 1.10, 1.50, 1.90, t = 1 with the error norms which are obtained by local discontinuous Galerkin method in ref. [1] are given in Table 7-9 respectively. Regarding that the results get better by increasing the partition number N in Galerkin method, these results are compared with the results that obtained by local discontinuous galerkin method and it is concluded that the results are better for the implemented method. Table 4: Error norms and numerical solutions of Problem 2 for γ = 1.50, ∆t = 0.0005, t = 1.   Table 6: Error norms and numerical solutions of Problem 2 for ∆t = 0.0005, N = 100, t = 1.

BMSA Volume 18
and the initial conditions as The exact solution of the problem is given by [1] U (x, t) = t 3 cos(x).
Lastly, the numerical results that obtained by applying quadratic B-spline Galerkin method to Problem 3 and the L 2 , L ∞ error norms are given in Table 10-12. As the N partition number increases, a significant amount of decreasement is seen in the error norms which are obtained for ∆t = 0.0002, t = 1, γ = 1.50 and different values of partition number N shown in Table 10. From Table 11 it is seen that the error norms which are obtained for N = 150, t = 1, γ = 1.50 and different values of ∆t decrease when ∆t time step decreases for ∆t, N = 30, t = 1 and γ = 1.50 as it is expected. In Table  12 the evaluated results which are found for ∆t = 0.0002, t = 1, N = 150, different values of γ are given and absolute errors obtained for the same values are given in Figure 3 graphically.   Table 12: Error norms and numerical solutions of Problem 3 for ∆t = 0.0002, N = 150, t = 1.

Conclusion
In this paper, the numerical solutions of three problems for the time fractional telegraph equation are achieved with aid of quadratic B-spline Galerkin method. The time fractional derivative operator is considered by means Caputo fractional derivative in these problems. It can be easily viewed from the numerical solutions, error norms and the tables that implemented method is an efficient tool to achieve numerical solutions of time fractional partial differential equations arising in physics and engineering.