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

Authors: O. Tasbozan; A. Esen; N. M. Yagmurlu; Y. Ucar
Journal: Abstract and Applied Analysis
 (2013)

---

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


---

## Body

## 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]
(1)∂∂tu(x,t)=K0Dt1-γ∂2∂x2u(x,t),
whereD0t1-γf(t) is the fractional derivative defined by Riemann-Liouville operator as
(2)0Dt1-γf(t)=1Γ(γ)∂∂t∫0t‍f(τ)(t-τ)1-γdτ,
where K is the diffusion coefficient and γ∈(0,1) is anomalous diffusion exponent. In all numerical computations, diffusion coefficient K is going to be taken as 1. In this paper, we will take the boundary conditions of (1) given in the interval 0≤x≤1 as
(3)u(0,t)=0,u(1,t)=0
and the initial condition as
(4)u(x,0)=x(1-x).
The exact analytical solution of (1) is found by the method of separation of variables [9] as
(5)u(x,t)=8π3∑‍∞n=01(2n+1)3sin[(2n+1)πx]×Eγ(-K(2n+1)2π2tγ),
where Eγ 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(0<γ≤1), we will also discretize the Riemann-Liouville operator [9, 16] as
(6)0Dt1-γu(xj,tn)=0Dt1-γujn=1(Δt)1-γ∑‍nk=0ωk1-γujn-k+O(Δtp),
where
(7)ω01-γ=1,ωk1-γ=(1-2-γk)ωk-11-γ.


## 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 [a,b] is partitioned into N finite elements of uniformly equal length by the knots xm, m=0,1,2,…,N such that a=x0<x1<⋯<xN=b and h=xm+1-xm. The cubic B-splines ϕm(x), (m=-1(1)N+1), over the interval [a,b] and at the knots xm are defined by [17](8)ϕm(x)=1h3{(x-xm-2)3,x∈[xm-2,xm-1],h3+3h2(x-xm-1)+3h(x-xm-1)2-3(x-xm-1)3,x∈[xm-1,xm],h3+3h2(xm+1-x)+3h(xm+1-x)2-3(xm+1-x)3,x∈[xm,xm+1],(xm+2-x)3,x∈[xm+1,xm+2],0otherwise.
The set of cubic splines {ϕ-1(x),ϕ0(x),…,ϕN+1(x)} constitutes a basis for the functions defined over the interval [a,b]. Thus, an approximate solution UN(x,t) over the interval can be written in terms of the cubic B-splines trial functions as
(9)UN(x,t)=∑m=-1N+1‍δm(t)ϕm(x),
where δm(t)'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 [xm,xm+1] is covered by four cubic B-splines. For this problem, the finite elements are identified with the interval [xm,xm+1] and the elements knots xm,xm+1. Using the nodal values Um, Um′, and Um′′ given in terms of the parameter δm(t)(10)Um=U(xm)=δm-1(t)+4δm(t)+δm+1(t),Um′=U′(xm)=3h(-δm-1(t)+δm+1(t)),Um′′=U′′(xm)=6h2(δm-1(t)-2δm(t)+δm+1(t)),
the variation of UN(x,t) over the typical element [xm,xm+1] can now be given by
(11)UN(x,t)=∑j=m-1m+2‍δj(t)ϕj(x).
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:
(12)δ˙m-1(t)+4δ˙m(t)+δ˙m+1(t)-6h2D0t1-γ[δm-1(t)-2δm(t)+δm+1(t)]=0,
where dot stands for derivative with respect to time. In the first place, time parameters δm(t) and their time derivatives δ˙m(t) in (12) are discretized by the following Crank-Nicolson formula and first order difference formula, respectively:(13)δ=12(δn+δn+1),δ˙=δn+1-δnΔt.
Then, if we discretize the fractional operatorD0t1-γδn by the following formula:
(14)0Dt1-γδn=1Δt1-γ∑‍nk=0ωk1-γδn-k,
we easily obtain a recurrence relationship between successive time levels relating unknown parameters δmn+1(t) as
(15)(1-α)δm-1n+1+(4+2α)δmn+1+(1-α)δm+1n+1=(1+α)δm-1n+(4-2α)δmn+(1+α)δm+1n+α∑‍nk=1ωk1-γ[(δm-1n+1-k+δm-1n-k)-2(δmn+1-k+δmn-k)+(δm+1n+1-k+δm+1n-k)],
where
(16)α=3(Δt)γh2.
The system (15) is consisted of N+1 linear equations including N+3 unknown parameters (δ-1,…,δN+1)T. 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 δN+1 from the system.


### 2.1. Initial State

To start iteration, we do need to evaluate the initial vector at starting time level. The initial vectord0=(δ0,δ1,δ2,…,δN-2,δN-1,δN)T is determined from the initial and boundary conditions. Thus, the approximation (11) can be rewritten for the initial condition as
(17)UN(x,0)=∑m=-1N+1‍δm(0)ϕm(x),
where the δm(0) are unknown element parameters. We force the initial numerical approximation UN(x,0) to meet the following conditions:
(18)UN(x,0)=U(xm,0),m=0,1,…,N,(UN)xx(0,0)=-2,(UN)xx(1,0)=-2.
Using these conditions results in a three-diagonal system of matrix of the form
(19)Wd0=b,
where
(20)W=[60141141⋱14106],b=(U(x0,0)+h23,U(x1,0),U(x2,0),…,U(xN-2,0),U(xN-1,0),U(xN,0)+h23)T.Solving this system yields the values of element parameters att=0. 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(21)δmn=ξneimφ,
where i=-1. First of all, substituting the Fourier mode (21) into the recurrence relationship (15) results in the following equation:
(22)ξn+1((1-α)e-iφ+(4+2α)+(1-α)eiφ)=ξn((1+α)e-iφ+(4-2α)+(1+α)eiφ)+α∑‍nk=1ωk1-γ[(ξn-k+1+ξn-k)(e-iφ-2+eiφ)].
Secondly, if we write
(23)ξn+1=ζξn
and assume that
(24)ζ≡ζ(φ)
is independent of time, then we get the following expression for the amplification factor ζ of the subdiffusion mode:
(25)ζ((1-α)e-iφ+(4+2α)+(1-α)eiφ)=((1+α)e-iφ+(4-2α)+(1+α)eiφ)+α∑‍nk=1ωk1-γ[(ζ1-k+ζ-k)(e-iφ-2+eiφ)].
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:
(26)8αsin2(φ2)≥0.
Since α>0, we can say that the scheme is unconditionally stable.


## 2.1. Initial State

To start iteration, we do need to evaluate the initial vector at starting time level. The initial vectord0=(δ0,δ1,δ2,…,δN-2,δN-1,δN)T is determined from the initial and boundary conditions. Thus, the approximation (11) can be rewritten for the initial condition as
(17)UN(x,0)=∑m=-1N+1‍δm(0)ϕm(x),
where the δm(0) are unknown element parameters. We force the initial numerical approximation UN(x,0) to meet the following conditions:
(18)UN(x,0)=U(xm,0),m=0,1,…,N,(UN)xx(0,0)=-2,(UN)xx(1,0)=-2.
Using these conditions results in a three-diagonal system of matrix of the form
(19)Wd0=b,
where
(20)W=[60141141⋱14106],b=(U(x0,0)+h23,U(x1,0),U(x2,0),…,U(xN-2,0),U(xN-1,0),U(xN,0)+h23)T.Solving this system yields the values of element parameters att=0. 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(21)δmn=ξneimφ,
where i=-1. First of all, substituting the Fourier mode (21) into the recurrence relationship (15) results in the following equation:
(22)ξn+1((1-α)e-iφ+(4+2α)+(1-α)eiφ)=ξn((1+α)e-iφ+(4-2α)+(1+α)eiφ)+α∑‍nk=1ωk1-γ[(ξn-k+1+ξn-k)(e-iφ-2+eiφ)].
Secondly, if we write
(23)ξn+1=ζξn
and assume that
(24)ζ≡ζ(φ)
is independent of time, then we get the following expression for the amplification factor ζ of the subdiffusion mode:
(25)ζ((1-α)e-iφ+(4+2α)+(1-α)eiφ)=((1+α)e-iφ+(4-2α)+(1+α)eiφ)+α∑‍nk=1ωk1-γ[(ζ1-k+ζ-k)(e-iφ-2+eiφ)].
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:
(26)8αsin2(φ2)≥0.
Since α>0, 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 normL2 as
(27)L2=∥Uexact-UN∥2=h∑j=0N‍|Ujexact-(UN)j|2
and the error norm L∞ as
(28)L∞=∥Uexact-UN∥∞=maxj|Ujexact-(UN)j|.Figures1 and 2 show the graphs of the exact (denoted by lines) solutions and the numerical ones for Δt=0.0001 and N=40 at t=0.001 (denoted by triangles), t=0.01 (denoted by squares), and t=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, Δt=0.001, and tf=0.1 for different values of N. The calculated error norms L2 and L∞ 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, N=40 and tf=0.1 for different values of Δt is illustrated and then the error norms L2 and L∞ are computed and presented in the table. In Table 3, we have listed the numerical and exact solutions of the problem for γ=0.75, Δt=0.001, and tf=0.1 for different values of N and the error norms L2 and L∞. Table 4 illustrates the comparison of the exact solutions with the numerical solutions for γ=0.75, N=40, and tf=0.1 for different values of Δt and the error norms L2 and L∞.Table 1
The comparison of the exact solutions with the numerical solutions withγ=0.5, Δt=0.001, and tf=0.1 for different values of N and the error norms L2 and L∞.
x
N
=
10
N
=
20
N
=
40
N
=
80
N
=
100
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
L
2
×
1
0
3
0.294172
0.159576
0.142246
0.139756
0.139501
L
∞
×
1
0
3
0.543026
0.299690
0.238942
0.223758
0.221936Table 2
The comparison of the exact solutions with the numerical solutions withγ=0.5, N=40, and tf=0.1 for different values of Δt and the error norms L2 and L∞.
x
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.0001
Δ
t
=
0.00001
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
L
2
×
1
0
3
0.836889
0.142246
0.136548
0.136264
L
∞
×
1
0
3
1.589369
0.238942
0.203235
0.199723Table 3
The comparison of the exact solutions with the numerical solutions withγ=0.75, Δt=0.001, and tf=0.1 for different values of N and the error norms L2 and L∞.
x
N
=
10
N
=
20
N
=
40
N
=
80
N
=
100
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
L
2
×
1
0
3
0.418167
0.174101
0.136253
0.130959
0.130466
L
∞
×
1
0
3
0.722901
0.336912
0.240523
0.216433
0.213542Table 4
The comparison of the exact solutions with the numerical solutions withγ=0.75, N=40, and tf=0.1 for different values of Δt and the error norms L2 and L∞.
x
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.0001
Δ
t
=
0.00001
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
L
2
×
1
0
3
0.449135
0.136253
0.128451
0.127966
L
∞
×
1
0
3
0.747758
0.240523
0.203042
0.199418Figure 1
The comparison of the exact (lines) and numerical solutions forγ=0.50, Δt=0.0001, and N=40 at t=0.001 (triangles), t=0.01 (squares), and t=0.1 (stars).Figure 2
The comparison of the exact and numerical solutions forγ=0.75, Δt=0.0001, and N=40 at at t=0.001 (triangles), t=0.01 (squares), and t=0.1 (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 normsL2 and L∞. 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.



---

*Source: 187383-2013-05-16.xml*