A Galerkin Finite Element Method to Solve Fractional Diffusion and Fractional Diffusion-Wave Equations

Abstract In the present study, numerical solutions of the fractional diffusion and fractional diffusion-wave equations where fractional derivatives are considered in the Caputo sense have been obtained by a Galerkin finite element method using quadratic B-spline base functions. For the fractional diffusion equation, the L1 discretizaton formula is applied, whereas the L2 discretizaton formula is applied for the fractional diffusion-wave equation. The error norms L 2 and L ∞ are computed to test the accuracy of the proposed method. It is shown that the present scheme is unconditionally stable by applying a stability analysis to the approximation obtained by the proposed scheme.


Introduction
Recently, it has turned out that many phenomena in engineering, physics, chemistry and other sciences can be described very successfully by models using mathematical tools from fractional calculus, i.e. the theory of derivatives and integrals of fractional (non-integer) order [10]. In fact, the concept of differentiation and integration to non-integer order is by no means new. Interest in this subject was evident almost as soon as the ideas of the classical calculus were known [12]. However, in the last few decades many authors have pointed out that derivatives and integrals of non-integer order are very suitable for the description of the behavior of various materials, e.g. polymers. It has been shown that new fractional-order models are more adequate than previously used integer-order models. The growing number of fractional derivative applications in various fields of science and engineering indicates that there is a significant demand for better mathematical models of real objects, and that the fractional calculus provides one possible approach on the way to more adequate mathematical modelling of real objects and processes. For example, the modelling of diffusion in a specific type of porous medium (in fractal media) is one of the most significant applications of fractional-order derivatives [13], and the fractional diffusion-wave equations have been proposed to deal with viscoelastic problems such as propagation of stress waves in viscoelastic solids [9]. Although there are several analytical techniques [6,16] for dealing with these fractional equations, as also happens with ordinary (non-fractional) partial differential equations, in many cases the initial condition, and/or the external force are such that the only reasonable option is to resort to numerical methods. However, although there have been an increasing number of works on this topic during the last few years [2,4,5,7,8,11], this field of applied mathematics is by far much less developed and understood than its non-fractional counterpart [20]. Although there have been many methods applied to solve fractional partial differential equations, there is still a long way to go in this field.
Many of the numerical methods for solving fractional partial differential equations that have been proposed differ essentially in the way in which the normal and fractional derivatives are discretized [20]. In this paper, the finite element method is applied to solve fractional differential equations, namely fractional diffusion and fractional diffusion-wave equations.
There are several studies about the fractional equations in the literature. For example, Sun et al. [18] have used a semi-analytical finite element method for a class of time-fractional diffusion equations. Sweilam et al. [19] solved timefractional diffusion equation by using Crank-Nicolson finite difference method. Monami and Adibat [10] have implemented relatively new analytical techniques, the variational iteration method and the Adomian decomposition method, for solving linear fractional partial differential equations arising in fluid mechanics. Celik and Duman [1] have used Crank-Nicolson method for the fractional diffusion equation with the Riesz fractional derivative and obtained numerical results using fractional centered difference approach. Murillo and Yuste [15] constructed the difference schemes using the L1 discretization formula for the fractional diffusion equation and the L2 discretization for the fractional diffusion wave equation.
The general form of the fractional diffusion and fractional diffusion-wave equations which are going to be used as a model is given by ∂τ n dτ, n − 1 < γ < n is the fractional derivative in the Caputo's sense [6,13], K is the diffusion coefficient and n is an integer. In the present paper, diffusion coefficient K is going to be taken as 1. For 0 < γ ≤ 1, Eq. (1.1) is called the fractional diffusion equation or sub-diffusion equation, and for 1 < γ ≤ 2 it is called the fractional diffusion-wave equation. Although various homogeneous and nonhomogeneous boundary conditions can be taken, in the present study, for the diffusion equation, we will take the boundary conditions of the model problem (1.1) given in the interval 0 ≤ x ≤ π as u(0, t) = 0, u(π, t) = 0 (1.2) and the initial condition as and for the diffusion-wave equation, along with the above boundary and the initial conditions, the following additional initial condition will be taken. The exact solution of the problems is found as follows [10,16] u( where E γ is the Mittag-Leffler function [13]. Fujita [3] has shown the existence and uniqueness of the solution of the Cauchy problem of the following form If we take β = 2, it results in fractional diffusion and fractional diffusion-wave equations given by Eq. (1.1). In our numerical solutions, to obtain a finite element scheme for solving the fractional diffusion equation (0 < γ ≤ 1), we will discretize the Caputo derivative by means of the so-called L1 formula [12] where b γ k = (k + 1) 1−γ − k 1−γ and to solve the fractional diffusion-wave equation (1 < γ ≤ 2), we will discretize the Caputo derivative by means of the so-called L2 formula [12] ∂ γ f ∂t γ

Quadratic B-Spline Galerkin Finite Element Solutions
Before starting to solve Eq. (1.1) with the boundary conditions (1.2) and the initial conditions (1.3)-(1.4) using Galerkin finite element method, we firstly define quadratic B-spline functions. Let us partition the interval [a, b] into M finite elements of uniformly equal length by the knots x m , m = 0, 1, 2, . . . , M such that a = x 0 < x 1 · · · < x M = b and h = x m+1 − x m . Since a quadratic Bspline covers three elements, we also introduce four additional grid points x −2 , } forms a base for the functions defined over [a, b]. Therefore, an approximation solution U M (x, t) can be written in terms of the quadratic B-splines trial functions as follows where δ m (t)'s are unknown, time dependent parameters to be determined from the boundary and weighted residual conditions. Each quadratic B-spline function covers three elements so that each element [x m , x m+1 ] is covered by three quadratic B-spline functions. For this problem, the finite elements are identified with the interval [x m , x m+1 ] and the element knots x m , x m+1 . Using the nodal values U m and U m given in terms of the parameter δ m (t) Before starting to apply the Galerkin method to the Eq. (1.1) with the appropriate boundary conditions, first of all we need to construct the weak form of the Eq. (1.1). For this purpose, all terms in Eq. (1.1) are taken to the right hand side of the equation and then multiplied by the weight function Ψ (x). Finally, by integrating the resulting equation over the region [0, π] and setting it to zero, we get where Ψ (x) is the weighted function taken as quadratic B-spline functions. When the differentiation is distributed between the approximate solution U and the weight function Ψ , it results in an integral form requiring weaker continuity conditions on trial base functions, and we obtain Since the weak form (2.2) is valid on the whole region, particularly it is valid over the typical element [x m , x m+1 ], thus Eq. (2.2) can be rewritten as follows To change from the global coordinate system to the local one, we use the The newly obtained Eq. (2.4) is the element equation for a typical element "e". Now, Eqs. (2.1) can be rewritten as follows (2.5) where dot denotes γ th fractional derivative with respect to time and A e ij , B e ij and C e ij are element matrices given by the following statements: By combining the coefficient matrix for each element in terms of global element parameters, Eq. (2.6) yields the system For the fractional diffusion equation (0 < γ ≤ 1), if time parameters δ(t)'s and its fractional time derivativesδ(t)'s in Eq. (2.7) are discretized by the Crank-Nicolson formula and L1 formula, respectively: we obtain a recurrence relationship between successive time levels relating unknown element parameters δ n+1 m (t) Next, for fractional diffusion-wave equation (1 < γ ≤ 2), if element parameters δ m (t)'s and its fractional time derivativesδ m (t)'s in Eq. (2.6) are discretized by the Crank-Nicolson formula and L2 formula, respectively we obtain a recurrence relationship between successive time levels relating element parameters δ n+1 m (t)

Stability analysis
The investigation of the stability of the approximation obtained by the scheme will be based on the von Neumann stability analysis. In the stability analysis, the growth factor of a typical Fourier mode is defined as: Next, by writing ξ n+1 = ζξ n and assuming that ζ ≡ ζ(ϕ) is independent of time, we easily get the following expression for the amplification factor ζof the sub-diffusion mode: According to the Fourier stability analysis, for the given scheme to be stable, the condition |ζ| ≤ 1 must be satisfied. Considering the extreme value ζ = 1, from the Eqs. (2.11) and (2.12), we obtain the following inequality α(2 sin 2 ϕ/2 + sin 2 (ϕ)) ≥ 0. (2.13) Substituting the Fourier mode (2.10) into the recurrence relationship (2.12), and following a similar way, we again obtain (2.13). Since α ≥ 0, both of the schemes are unconditionally stable.

Numerical Examples and Results
Numerical results of the diffusion and diffusion-wave problems are obtained by Galerkin finite element method using quadratic B-spline base functions. The accuracy of the method is measured by the discrete norm L 2 and the maximum nodal norm L ∞ In Fig. 1 a, we present numerical solutions of the fractional diffusion problem at the midpoint for various values of γ and M = 40.
The comparison of the analytical solution and numerical solutions obtained for diffusion equation for values of γ = 0.25, γ = 0.50 and γ = 0.75 is given in Table 1. As it is clearly seen from the table, the analytical and numerical solutions obtained by the present scheme are in good agreement with each other. In Table 2, we demonstrate the numerical results for γ = 0.5, ∆t = 0.00007 and t f = 0.35 and for different number of divisions of the region. Table 2 clearly shows that as the number of division increases, the obtained numerical results become more accurate. We see this from the decreasing values of the discrete norm L 2 and maximum nodal norm L ∞ . In Fig. 1b, we demonstrate the graphs of numerical solutions obtained for γ = 0.50 and M = 40 at different time levels. In Table 3, it is clearly seen that as the values of time steps decrease, the agreement between the approximate solutions and analytical solutions becomes better, and the values of the discrete norm L 2 and maximum nodal norm L ∞ become smaller. In Table 4, we demonstrate the discrete norm L 2 and maximum nodal norm L ∞ for γ = 0.25, γ = 0.50 and γ = 0.75 at various time levels.
In Fig. 2a, the numerical solutions of the fractional diffusion-wave problem at the midpoint for various values of γ and M = 40 are illustrated. The comparison of the analytical solution and numerical solutions obtained by the scheme for diffusion-wave equation for values of γ = 1.25, γ = 1.50 and γ = 1.75 is given in Table 5. The table clearly demonstrates that the obtained numerical results are satisfactorily in good agreement with the analytical ones. As the value of γ increases, so the values of the discrete error norm L 2 and maximum nodal norm L ∞ . In Table 6, the numerical results for γ = 1.5, ∆t = 0.0015 and t f = 3.75 for various values of M are given. As it is seen from the table, as the number of division increases, the values of error decrease. Table 7 shows that  as the values of time steps decrease, the agreement between the approximate solutions and analytical solutions becomes better, and the values of the discrete   error norm L 2 and maximum nodal norm L ∞ become smaller. In Fig. 2b, the numerical solutions for the values of γ = 1.50 and M = 40 at various time levels are presented. In Table 8, we show the discrete error norm L 2 and maximum nodal norm L ∞ for M = 40, ∆t = 0.0015 and t f = 3.75 at different values of γ and t.

Conclusions
In the present study, a Galerkin finite element method has been successfully used to obtain the numerical solutions of diffusion and diffusion-wave equations. In these equations, the fractional derivative is considered of the Caputo form.
The fractional derivative appearing in the fractional diffusion and diffusionwave equations is approximated, respectively, by means of the so-called L1 and L2 formulae the same as used by Ref. [15] in the explicit finite difference method solution. One can easily conclude from the presented results that the applied method is a highly good one to obtain numerical solutions of this kind fractional partial differential equations.