Abstract and Applied Analysis

Volume 2013 (2013), Article ID 187383, 6 pages

http://dx.doi.org/10.1155/2013/187383

Research Article

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

O. Tasbozan, A. Esen, N. M. Yagmurlu, and Y. Ucar

Department of Mathematics, Faculty of Art and Sciences, Inonu University, 44280 Malatya, Turkey

Received 8 February 2013; Accepted 27 April 2013

Academic Editor: Mustafa Bayram

Copyright © 2013 O. Tasbozan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

A collocation finite element method for solving fractional diffusion equation for force-free case is considered. In this paper, we develop an approximation method based on collocation finite elements by cubic B-spline functions to solve fractional diffusion equation for force-free case formulated with Riemann-Liouville operator. Some numerical examples of interest are provided to show the accuracy of the method. A comparison between exact analytical solution and a numerical one has been made.

1. 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–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 is the fractional derivative defined by Riemann-Liouville operator as

where is the diffusion coefficient and 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 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 force-free case , we will also discretize the Riemann-Liouville operator [9, 16] as

where

2. 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 , such that and . The cubic B-splines , (), over the interval and at the knots are defined by [17]

The set of cubic splines 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 is covered by four cubic B-splines. For this problem, the finite elements are identified with the interval and the elements knots . Using the nodal values , , and given in terms of the parameter

the variation of over the typical element 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 by the following formula:

we easily obtain a recurrence relationship between successive time levels relating unknown parameters as

where

The system (15) is consisted of linear equations including unknown parameters . 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 and from the system.

2.1. Initial State

To start iteration, we do need to evaluate the initial vector at starting time level. The initial vector is determined from the initial and boundary conditions. Thus, the approximation (11) can be rewritten for the initial condition as

where the are unknown element parameters. We force the initial numerical approximation to meet the following conditions:

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 . Now, it is time to find out the values of element parameters at different time levels using the iterative system (15).

2.2. 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 . 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

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 must be satisfied. Considering the extreme value , from (22) and (25), we obtain the following inequality:

Since , we can say that the scheme is unconditionally stable.

3. Numerical Examples and Results

Numerical results for the diffusion and diffusion-wave problems are obtained by collocation method using cubic B-spline base functions. The accuracy of the present method is measured by the error norm as

and the error norm as

Figures 1 and 2 show the graphs of the exact (denoted by lines) solutions and the numerical ones for and at (denoted by triangles), (denoted by squares), and (denoted by stars) for two different values of and , respectively. Table 1 shows the comparison of the exact solutions with the numerical ones with , , and for different values of . The calculated error norms 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 , and for different values of is illustrated and then the error norms and are computed and presented in the table. In Table 3, we have listed the numerical and exact solutions of the problem for , , and for different values of and the error norms and . Table 4 illustrates the comparison of the exact solutions with the numerical solutions for , , and for different values of and the error norms and .

Table 1: The comparison of the exact solutions with the numerical solutions with , , and for different values of and the error norms and .

Exact

0.0
0.000000
0.000000
0.000000
0.000000
0.000000
0.000000

0.1
0.013806
0.013886
0.013920
0.013925
0.013926
0.013765

0.2
0.026128
0.026284
0.026325
0.026335
0.026336
0.026183

0.3
0.035777
0.035982
0.036034
0.036046
0.036048
0.036037

0.4
0.041903
0.042137
0.042195
0.042210
0.042212
0.042364

0.5
0.044001
0.044245
0.044306
0.044321
0.044323
0.044544

0.6
0.041903
0.042137
0.042195
0.042210
0.042212
0.042364

0.7
0.035777
0.035982
0.036034
0.036046
0.036048
0.036037

0.8
0.026128
0.026284
0.026325
0.026335
0.026336
0.026183

0.9
0.013806
0.013886
0.013920
0.013925
0.013926
0.013765

1.0
0.000000
0.000000
0.000000
0.000000
0.000000
0.000000

0.294172
0.159576
0.142246
0.139756
0.139501

0.543026
0.299690
0.238942
0.223758
0.221936

Table 2: The comparison of the exact solutions with the numerical solutions with , , and for different values of and the error norms and .

Exact

0.0
0.000000
0.000000
0.000000
0.000000
0.000000

0.1
0.015039
0.013920
0.013929
0.013930
0.013765

0.2
0.026111
0.026325
0.026347
0.026349
0.026183

0.3
0.035392
0.036034
0.036063
0.036065
0.036037

0.4
0.041638
0.042195
0.042229
0.042233
0.042364

0.5
0.043830
0.044306
0.044341
0.044345
0.044544

0.6
0.041638
0.042195
0.042229
0.042233
0.042364

0.7
0.035392
0.036034
0.036063
0.036065
0.036037

0.8
0.026111
0.026325
0.026347
0.026349
0.026183

0.9
0.015039
0.013920
0.013929
0.013930
0.013765

1.0
0.000000
0.000000
0.000000
0.000000
0.000000

0.836889
0.142246
0.136548
0.136264

1.589369
0.238942
0.203235
0.199723

Table 3: The comparison of the exact solutions with the numerical solutions with , , and for different values of and the error norms and .

Exact

0.0
0.000000
0.000000
0.000000
0.000000
0.000000
0.000000

0.1
0.018544
0.018677
0.018710
0.018718
0.018719
0.018574

0.2
0.035151
0.035393
0.035453
0.035468
0.035470
0.035331

0.3
0.048210
0.048531
0.048611
0.048631
0.048633
0.048629

0.4
0.056530
0.056900
0.056992
0.057015
0.057018
0.057166

0.5
0.048210
0.059771
0.059868
0.059892
0.059895
0.060108

0.6
0.056530
0.056900
0.056992
0.057015
0.057018
0.057166

0.7
0.048210
0.048531
0.048611
0.048631
0.048633
0.048629

0.8
0.035151
0.035393
0.035453
0.035468
0.035470
0.035331

0.9
0.018544
0.018677
0.018710
0.018718
0.018719
0.018574

1.0
0.000000
0.000000
0.000000
0.000000
0.000000
0.000000

0.418167
0.174101
0.136253
0.130959
0.130466

0.722901
0.336912
0.240523
0.216433
0.213542

Table 4: The comparison of the exact solutions with the numerical solutions with , , and for different values of and the error norms and .

Exact

0.0
0.000000
0.000000
0.000000
0.000000
0.000000

0.1
0.018556
0.018710
0.018722
0.018723
0.018574

0.2
0.035089
0.035453
0.035475
0.035477
0.035331

0.3
0.048219
0.048611
0.048642
0.048644
0.048629

0.4
0.056522
0.056992
0.057028
0.057031
0.057166

0.5
0.059360
0.059868
0.059905
0.059909
0.060108

0.6
0.056522
0.056992
0.057028
0.057031
0.057166

0.7
0.048219
0.048611
0.048642
0.048644
0.048629

0.8
0.035089
0.035453
0.035475
0.035477
0.035331

0.9
0.018556
0.018710
0.018722
0.018723
0.018574

1.0
0.000000
0.000000
0.000000
0.000000
0.000000

0.449135
0.136253
0.128451
0.127966

0.747758
0.240523
0.203042
0.199418

Figure 1: The comparison of the exact (lines) and numerical solutions for , , and at (triangles), (squares), and (stars).

Figure 2: The comparison of the exact and numerical solutions for , , and at at (triangles), (squares), and (stars).

4. 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 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.

References

R. Hilfer, Ed., Applications of Fractional Calculus in Physics, World Scientific Publishing, River Edge, NJ, USA, 2000.

I. M. Sokolov, J. Klafter, and A. Blumen, “Fractional kinetics,” Physics Today, vol. 55, no. 11, pp. 48–54, 2002.

A. Iomin, S. Dorfman, and L. Dorfman, “On tumor development: fractional transport approach,” http://arxiv.org/abs/q-bio/0406001

R. Metzler and J. Klafter, “The random walk's guide to anomalous diffusion: a fractional dynamics approach,” Physics Reports, vol. 339, no. 1, p. 77, 2000.

R. Metzler and J. Klafter, “The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics,” Journal of Physics A, vol. 37, no. 31, pp. R161–R208, 2004.

R. Metzler, E. Barkai, and J. Klafter, “Anomalous diffusion and relaxation close to thermal equilibrium: a fractional Fokker-Planck equation approach,” Physical Review Letters, vol. 82, no. 18, pp. 3563–3567, 1999.

E. Barkai, R. Metzler, and J. Klafter, “From continuous time random walks to the fractional Fokker-Planck equation,” Physical Review E, vol. 61, no. 1, pp. 132–138, 2000.

K. Diethelm and A. D. Freed, “On the solution of nonlinear fractional order differential equations used in the modelling of viscoplasticity,” in Scientific Computing in Chemical Engineering II:, Computational Fluid Dynamics, Reaction Engineering and Molecular Properties, pp. 217–224, Springer, Heidelberg, Germany, 1999.

S. B. Yuste, “Weighted average finite difference methods for fractional diffusion equations,” Journal of Computational Physics, vol. 216, no. 1, pp. 264–274, 2006.

T. A. M. Langlands and B. I. Henry, “The accuracy and stability of an implicit solution method for the fractional diffusion equation,” Journal of Computational Physics, vol. 205, no. 2, pp. 719–736, 2005.

D. A. Murio, “Implicit finite difference approximation for time fractional diffusion equations,” Computers & Mathematics with Applications, vol. 56, no. 4, pp. 1138–1145, 2008.

S. B. Yuste and L. Acedo, “An explicit finite difference method and a new von Neumann-type stability analysis for fractional diffusion equations,” SIAM Journal on Numerical Analysis, vol. 42, no. 5, pp. 1862–1874, 2005.

V. Balakrishman, “Anomalous diffusion in one dimension,” Phsica A, vol. 132, pp. 569–580, 1985.

W. Wyss, “The fractional diffusion equation,” Journal of Mathematical Physics, vol. 27, no. 11, pp. 2782–2785, 1986.

L. Podlubny, Fractional Differantial Equations, Academic Press, London, UK, 1999.

S. B. Yuste and L. Acedo, “An explicit finite difference method and a new von Neumann-type stability analysis for fractional diffusion equations,” SIAM Journal on Numerical Analysis, vol. 42, no. 5, pp. 1862–1874, 2005.

P. M. Prenter, Splines and Variational Methods, Pure and Applied Mathematics, John Wiley & Sons, New York, NY, USA, 1975.