Numerical approach for solving time fractional diffusion equation

Article history: Received:16 June 2017 Accepted: 13 November 2017 Available Online: 21 November 2017 In this article one of the fractional partial differential equations was solved by finite difference scheme based on five point and three point central space method with discretization in time. We use between the Caputo and the Riemann-Liouville derivative definition and the Grünwald-Letnikov operator for the fractional calculus. The stability analysis of this scheme is examined by using von-Neumann method. A comparison between exact solutions and numerical solutions is made. Some figures and tables are included.


Introduction
Anomalous equation is a diffusion process with a nonlinear relationship to time contrary to typical diffusion process, connected with the interactions within the complex and non-homogeneous background.In contrast to typical diffusion, anomalous diffusion is described by a power law [1].Anomalous diffusion has an important role in the literature to describe many physical phenomena, crowded systems or diffusion through porous media.On the other hand this phenomena is observed in heat baths [2], diffusion through also porous material [3,4], nuclear magnetic resonance diffusometry in disordered materials [4], behaviour of polymers in a glass transition [5,6]and also particle dynamics inside polymer network [7].Fractional order linear and nonlinear differential equations were examined by some researchers by using different methods [8][9][10][11][12].In this study we consider the fractional subdiffusion equation in the following form [13]: Recently, because of their practical applications, anomalous subdiffusion equation has received much attention.The common methods for solving fractional-order equations are purely mathematical [14], in terms of Mittag Lefler function [15], Green's function solution [16], numerical methods [17], variational iteration method [18] and differential transformations [19].In the present work we concentrate on one of the numerical methods which is fractional finite difference method for solutions of equation (1).In section 2, we give some information about fractional calculus, some definition of fractional derivative which are used in the next sections.We talk about finite difference method, stability of this explicit scheme and calculate the maximum absolute error ( max E ) at the time t and the root mean square error RMS E at the same time in section 3.In the last section 4, we apply the method to fractional diffusion equation with the given initial and boundary conditions and compared the numerical solutions with the exact results.

Basic concepts of fractional calculus
In this paper we consider the fractional partial differential Eq.( 1), with the boundary and initial-value conditions: Fractional calculus includes different definition of fractional derivative operator such as Riemann-Liouville derivative, the Caputo derivative, the Grünvald Letnikov derivative, the Riesz derivative and also the Weyl-Marchaud derivative.
We introduced the definition of Caputo operator [1], , . denotes an integer part of a real number.The operator 0 t D  is defined in the Riemann-Liouville sense as [1], According to fractional calculus between the Caputo and the Riemann-Liouville derivative is following form [1], We divide the spatial domain   0, L into the uniform mesh with N subintervals.

j M 
We mention another definition of the fractional derivative which is Grünwald-Letnikov sense as the following [1]: This definition is equivalent to the Riemann-Liouville definition (3).However the Grünwald-Letnikov operator is more suitable and more obvious in numerical calculations.The Grünwald-Letnikov operator Eq.( 6) is approximated with the interval [0,T] with the subinterval length t  as [1], where By using the recurrence relationship we obtain, These coefficients come from the generating function [22], If we put Eq.( 7) in place of Eq.( 5), we obtain the following equation, We will use this fractional derivative definition for the left-side of Eq (1).

Finite difference scheme of the fractional diffusion equation
In this work we take initial value conditions as in the Eq.( 2).If we consider them as discrete forms . The second derivative at the space in the right side of Eq.( 1) in terms of 5-point central differences can be written as, In case of 1 i  and 9 i  for the second derivative at the space we can use the following 3-point central difference formula, ,  , .
and substituting Eq.( 10) and Eq.( 11) into Eq.( 1) we obtain, Disintegrating partially into terms the first sum in Eq.( 13) we have, Numerical approach for solving time fractional diffusion equation

283
If we write more clearly, ) ) To obtain the formula of , ij u for 1 and have to use 3-point central differences.Substituting Eq.( 10) and Eq.( 12) into Eq.( 1) we obtain, Disintegrating partially into terms the first sum in Eq.( 17) and forming we obtain the following equation; and we have, As usual for explicit methods, the present explicit difference schemes Eq.( 16) and Eq.( 20) is not unconditionally stable.Therefore it is important to determine the conditions under which these schemes are stable.

Stability of method
The present explicit difference schemes, Eq.( 16) and Eq.( 20) is not unconditionally stable.Therefore it is important to determine the conditions under which the method is stable.We will employ that fractional von Neumann stability (or Fourier stability) analysis put forward in [22].We start by assuming a solution of our problem with the form () j q i ph i u e    , where q is a real spatial number.Inserting this expression into Eq.(16) we get, In this study, because of the order of fractional differential equation 0 m  then, 0 j p and the source term is neglected (   , f x t =0) and the stability of the process is determined by the behaviour of and assuming  is independent of time , we get, where If 1   for some q , the temporal factor of the solution grows to infinity according to Eq.( 22) and the mode is unstable.Considering the extreme value 1   , we obtain from Eq.( 23) the following stability bound on S [22]: x j S approaches lim . The value of x S can be concluded from Eq.( 24) by taking into consideration that [23], and since we use the sufficient condition for the present method to be stable is x SS  [22].Therefore one is inside the stable region and gets reasonable numerical solution.
To verify the numerical accuracy in the event of initial value problems on the time interval   0,T we calculate the max and the corresponding RMS error, and by the end of calculations we have, max 0.0020 and 0.0024

Numerical example
Let us consider the fractional subdiffusion equation (1) with the boundary conditions The exact solution of this problem is given by [8], For computational work, we have taken    In Fig. 2 we observe numerical solutions of fractional diffusion equation (Eq.1) at different times 0. These solutions given by Table 4.In Figure 3 we observe numerical solutions of fractional diffusion equation (Eq.( 1)) for different orders 0.85  In Figure 4 we observe absolute errors of the numerical method for the Eq.( 1

Conclusion
We introduced numerical schemes for solving time fractional diffusion equation a new finite difference method which is an explicit method.This method involves fractional Riemann-Liouville derivatives, which were used with Grünwald Letnikov formula.In this study, we have used 5-point and 3-point difference scheme while discretizing of equation.For the all numerical algorithm we used Maple 15 algebraic system with acceptable CPU time and for the graphics of the solution Matlab 2009 has been used.On the other hand it was also used that a fractional von-Neumann stability analysis which leads to correct estimate of the stability conditions for diffusion equations.According to this analysis we found a stability bound on   2 12h t S    .Furthermore we investigated absolute and RMS errors, they are compatible.The numerical and exact solutions of Eq.(1) at several times and for different orders have been shown in tables and solutions were supported by the graphics.These results show that the present method gives us new solutions in agreement with the analytical solutions.Hence, we conclude that the proposed approach is very powerful and efficient to investigate fractional differential equation numerically.

1 
is a field variable (probability density of diffusive displacements x in a time t), K  is the diffusion coefficient, 0 c t D  is the fractional derivative in the Caputo sense, α is a real order of this operator and   , f x t is a suitable right hand side function.The time fractional subdiffusion equation has been obtained by replacing the time derivative in ordinary diffusion by a fractional order.In a typical diffusion process super diffusion equation.

Figure 2 .
Figure 2. Numerical solution of Eq.(1) for various values of t .

Figure 4 .
Figure 4. Absolute errors for various values of  .

Table 2 .
Numerical results for 0.90

Table 3 .
Numerical results for 0.95