A Numerical Solution to Fractional Diffusion Equation for Force-Free Case

and Applied Analysis 3


Introduction
Scientific and engineering problems including fractional derivatives have become more important in recent years.Since the description of physical and chemical processes by means of equations including fractional derivatives is more accurate and precise, their numerical solutions have been the primary interest of many recently published articles.The applications are so wide that they include such diverse areas as control theory [1], transport problems [2], tumor development [3], subdiffusive anomalous transport in the presence of an external field [4][5][6][7], and viscoelastic and viscoplastic flow [8].These diverse areas of applications have led to an increase in the number of studies on fractional differential equations and have caused it to be an important topic in mathematics and science.Yuste [9] has used weighted average finite difference methods for fractional diffusion equations and provided some examples in which the new methods' numerical solutions are obtained and compared against exact solutions.Langlands and Henry [10] have investigated the accuracy and stability of an implicit numerical scheme for solving the fractional diffusion equation.Murio [11] has developed an implicit unconditionally stable numerical method to solve the one-dimensional linear time fractional diffusion equation, formulated with Caputo's fractional derivative, on a finite slab.Yuste and Acedo [12] have combined the forward time centered space (FTCS) method, well known for the numerical integration of ordinary diffusion equations, with the Grünwald-Letnikov discretization of the Riemann-Liouville derivative to obtain an FTCS scheme for solving the fractional diffusion equation.
The general form of the fractional diffusion equation for force-free case is given by [4,13,14] where 0  1−  () is the fractional derivative defined by Riemann-Liouville operator as where  is the diffusion coefficient and  ∈ (0, 1) is anomalous diffusion exponent.In all numerical computations, diffusion coefficient  is going to be taken as 1.In this paper, we will take the boundary conditions of (1) given in the interval 0 ≤  ≤ 1 as and the initial condition as The exact analytical solution of (1) is found by the method of separation of variables [9] as where   is the Mittag-Leffler function [15].
In our numerical solutions, to obtain a finite element scheme for solving the fractional diffusion equation for forcefree case (0 <  ≤ 1), we will also discretize the Riemann-Liouville operator [9,16] as where

Cubic B-Spline Finite Element Collocation Solutions
To solve (1) with the boundary conditions (3) and the initial condition (4) using collocation finite element method, first of all, we define cubic B-spline base functions.Let us consider that the interval [, ] is partitioned into  finite elements of uniformly equal length by the knots   ,  = 0, 1, 2, . . .,  such that  =  0 <  1 < ⋅ ⋅ ⋅ <   =  and ℎ =  +1 −   .The cubic B-splines   (), ( = −1(1) + 1), over the interval [, ] and at the knots   are defined by [17]   () The set of cubic splines { −1 (),  0 (), . . .,  +1 ()} constitutes a basis for the functions defined over the interval [, ].Thus, an approximate solution   (, ) over the interval can be written in terms of the cubic B-splines trial functions as where   ()'s are unknown, time dependent quantities to be determined from the boundary and cubic B-spline collocation conditions.Due to the fact that each cubic B-spline covers four consecutive elements, each element [  ,  +1 ] is covered by four cubic B-splines.For this problem, the finite elements are identified with the interval [  ,  +1 ] and the elements knots   ,  +1 .Using the nodal values   ,    , and    given in terms of the parameter   () the variation of   (, ) over the typical element [  ,  +1 ] can now be given by If we substitute the global approximation (11) and its necessary derivatives ( 10) into (1), we directly obtain the following set of the first-order ordinary differential equations: where dot stands for derivative with respect to time.In the first place, time parameters   () and their time derivatives δ  () in ( 12) are discretized by the following Crank-Nicolson formula and first order difference formula, respectively: Then, if we discretize the fractional operator 0  1−    by the following formula: we easily obtain a recurrence relationship between successive time levels relating unknown parameters  +1  () as where The system ( 15) is consisted of +1 linear equations including  + 3 unknown parameters ( −1 , . . .,  +1 )  .To obtain a unique solution to this system, we need two additional constraints.These are obtained from the boundary conditions and can be used to eliminate  −1 and  +1 from the system.

Initial State.
To start iteration, we do need to evaluate the initial vector at starting time level.The initial vector d 0 = ( 0 ,  1 ,  2 , . . .,  −2 ,  −1 ,   )  is determined from the initial and boundary conditions.Thus, the approximation (11) can be rewritten for the initial condition as where the   (0) are unknown element parameters.We force the initial numerical approximation   (, 0) to meet the following conditions: (, 0) =  (  , 0) ,  = 0, 1, . . ., , Using these conditions results in a three-diagonal system of matrix of the form where Solving this system yields the values of element parameters at  = 0. Now, it is time to find out the values of element parameters at different time levels using the iterative system (15).

Stability Analysis.
The study of the stability of the approximation obtained by the present scheme will be based on the von Neumann stability analysis.In this analysis, the growth factor of a typical Fourier mode is defined as where  = √ −1.First of all, substituting the Fourier mode (21) into the recurrence relationship (15) results in the following equation: Secondly, if we write and assume that  ≡  () (24) is independent of time, then we get the following expression for the amplification factor  of the subdiffusion mode: If we want the given scheme to be stable in terms of Fourier stability analysis, then the condition || ≤ 1 must be satisfied.
Considering the extreme value  = 1, from ( 22) and (25), we obtain the following inequality: Since  > 0, we can say that the scheme is unconditionally stable.

Numerical Examples and Results
Numerical results for the diffusion and diffusion-wave problems are obtained by collocation method using cubic Bspline base functions.
Figures 1 and 2 show the graphs of the exact (denoted by lines) solutions and the numerical ones for Δ = 0.0001 and  = 40 at  = 0.001 (denoted by triangles),  = 0.01 (denoted by squares), and  = 0.1 (denoted by stars) for two different values of  = 0.50 and  = 0.75, respectively.Table 1 shows the comparison of the exact solutions with the numerical ones with  = 0.5, Δ = 0.001, and   = 0.1 for different values of .The calculated error norms  2 and  ∞ at those time levels are also presented in the table.In Table 2, the comparison of the exact solutions with the numerical ones with  = 0.5,  = 40 and   = 0.1 for different values of Δ is illustrated and then the error norms  2 and  ∞ are computed and presented in the table.In Table 3, we have listed the numerical and exact solutions of the problem for  = 0.75, Δ = 0.001, and   = 0.1 for different values of  and the error norms  2 and  ∞ .Table 4 illustrates the comparison of the exact solutions with the numerical solutions for  = 0.75,  = 40, and   = 0.1 for different values of Δ and the error norms  2 and  ∞ .

Conclusion
In the present study, first of all, a collocation finite element method has been constructed.Then, the method has been applied using cubic B-spline base functions.During the implementation of the method, Crank-Nicolson formula and first-order difference formula have been applied for discretization process.The stability of the method presented in the paper has been tested using the von Neumann stability analysis in which the growth factor of a typical Fourier mode is used.The accuracy of the method is also measured by the error norms  2 and  ∞ .The successful application of the present method prompts the probability of extending it to other finite element methods and other kinds of fractional differential equations.The available results suggest that this is highly probable.

Table 1 :
The comparison of the exact solutions with the numerical solutions with  = 0.5, Δ = 0.001, and   = 0.1 for different values of  and the error norms  2 and  ∞ .

Table 2 :
The comparison of the exact solutions with the numerical solutions with  = 0.5,  = 40, and   = 0.1 for different values of Δ and the error norms  2 and  ∞ .

Table 3 :
The comparison of the exact solutions with the numerical solutions with  = 0.75, Δ = 0.001, and   = 0.1 for different values of  and the error norms  2 and  ∞ .

Table 4 :
The comparison of the exact solutions with the numerical solutions with  = 0.75,  = 40, and   = 0.1 for different values of Δ and the error norms  2 and  ∞ .
The accuracy of the present method is measured by the error norm  2 as ∞ =       exact −       ∞ = max         exact  − (  )        .