1.
Introduction
Fractional differential equations have been applied in modelling of many fields of physics and processes such as earthquakes, optics, finance, hydrology, traffic flow, fluid mechanics, fractional kinetics, mathematical biology, measurement of visco-elastic material properties, electrical network, electro-chemistry electro-magnetic, signal processing, control theory, acoustics, material sciences [2,5,8,9,17,22,39,41,42,43,44,45,46]. The nonlocal property is the most important advantage of these equations showing state of a complex system does not depend only on its current state but also depends on its all-previous states. Due to this reason, the fractional calculus has becoming more and more popular.
The fractional diffusion equations (FDEs) are generalized forms of classical diffusion equations. These equations have been solved by many researchers. For example, authors of [6,11,24,26,32] have been applied finite element method, compact FDM, Crank-Nicolson FDM, B-spline based method, implicit difference approximation, respectively, to solve FDEs. Khader [14] considered an efficient method based upon Chebyshev approximations. Murio [20] developed an implicit unconditionally stable method. Sun et al. [25] applied a semi-analytical FEM to solve a class of these equations. Tadjeran et al. [27] applied aforesaid method with spatial extrapolation to solve a class of variable coefficients FDEs. Lin and Xu [15] constructed a method based on FDM and Legendre spectral method while Ghanbari and Atangana [7] presented a method based on the product-integration rule for solving aforesaid equations. Çelik and Duman [3] examined a Crank-Nicolson method with the Riesz fractional derivative. Zhai and Feng [35] introduced a block-centered FDM for solving these equations.
Additionally, Zhai et al. [34,36] constructed unconditionally stable FDMs to solve the time-space fractional diffusion and three-dimensional time-fractional subdiffusion equations, respectively. Wu and Zhai [33] proposed a high order FDM to solve the 2D time-fractional convection-dominated diffusion equation while Zhai et al. [37,38] proposed a high-order compact FDM and ADI method to solve the 3D fractional convection-diffusion equation. Furthermore, Hanert [10] proposed a flexible numerical method to discretize the space-time FDEs, where pseudo-spectral expansion has been used for discretization of the time derivative and either high order pseudo-spectral or low order FE and FD have been used for space derivative. Murillo and Yuste [19] considered explicit difference scheme, where three-point centered formula has been used to approximate the spatial derivative. Verma et al. [31] studied nonlinear diffusion equations analytically and numerically using the classical Lie symmetry method.
In this paper, we consider a Caputo time-fractional diffusion equation
subject to the initial condition
the Dirichlet boundary conditions
where cDαtu(x,tn) denotes the Caputo time-fractional derivative [40] as
2.
Discretization of the problem
First, the Caputo time-fractional derivative is discretized. For this, we partition the time domain [0,T] uniformly as 0=t0<t1,...,<tN=T, where interval Δt=tn−tn−1=TN for n=1,2,…,N. Here N is the number of time mesh. By the definition of Caputo derivative at time t=tn, the time-fractional derivative is approximated as
where ˉbl=(l+1)1−α−l1−α ∀ 0≤l≤N and TrnΔt gives the truncation error ≤ˉcuΔt2−α, where the constant ˉcu is only related to u(x,t).
Lemma 3.1 The below relations should be hold for the coefficients ˉbl
(1)ˉb0=1,
(2)ˉbl>0,∀ 0≤l≤N,
(3)ˉbl−1>ˉbl, ∀ 1≤l≤N.
Proof. For the proof of Lemma 3.1, see [16].
Next, we use CExpB-spline collocation technique to discretize the space derivatives. First, we partition the space domain [a,b] uniformly as a=x0<x1,...,xM=b with space size Δx=h=xm+1−xm=b−aM where m=0,1,…,M. For discrete form, we denote unm=u(xm,tn) for m=0,1,...,M and n=0,1,...,N. The CExpB-spline functions Epm(x) for m=−1, 0, 1,…,M+1 are defined as [28,29,30]:
where
In Eq (2.2), the free parameter ˆκ have been used for obtaining the different forms of CExpB-spline functions. The set of Epm(x) ∀ m=−1, 0, 1,…,M+1 forms a basis over the problem domain. We assume that the approximation uM to the exact u(x,t) at the point (xm,tn) is expressed in terms of linear combinations of the CExpB-spline functions and unknown time-dependent quantities as follows:
where Ci(t) are the unknown quantities which we have to evaluate for the approximated solution uM(x,t) at (xi,tj). Since each CExpB-spline covers four elements, so each element is covered by four CExpB-splines. So the variation of the uM(x,t), over the element, can be written as:
Using Eq (2.4), the u(x,t) and its first two derivatives at the knots in terms of Cnm are given as:
Now, we modify the CExpB-spline basis functions which generate a new set of CExpB-spline basis functions. The procedure, for modifying the basis functions, is given as [1,12,13]:
where {^Ep0(x),^Ep1(x)……^EpM(x)} forms a basis over the problem domain. Next, the modified form of the approximated solutions, as a linear combination of modified CExpB-spline functions, is given by
Next, the initial vector (C0m, m=0, 1,…, M) can be obtained by initial condition and boundary values of its derivatives as:
Eq (2.10) yields a (M+1)×(M+1) system as:
Now, discretizing the Eq (1.1) by using (2.1) and Crank-Nicolson method, we have
Using Eqs (2.4) and (2.6) in Eq (2.12) and simplifying the terms, we have
where
For m=0, from Eq (2.13), we get
For m = M, from Eq (2.13), we get
At time step tn, n=1,2,...,N, the Eqs (2.13)–(2.15) can be formulated into a (M+1)×(M+1) linear system as follows:
where
3.
Stability analysis
The von-Neumann stability analysis [4,18,21,23] have been done in this section. Taking f(x,t)=0 in Eq (1.1) and discretizing, we get
Simplifying and rearranging the terms, we get
which can be written as
where
and
Now we consider one Fourier mode out of the full solution Cnm=δnekiϕ as trial solutions at a given point xi, where ϕ=θh. The θ and h are the mode number and element size respectively, and k=√−1. Substituting the trial solution in above equation and simplifying the terms, we get
where ˉF=n−1∑l=1(ˉbl−1−ˉbl)δn−l+ˉbn−1δ0.
Now we define
Using (3.4), we get
Now, using above equation and the values of A, B, D and E in Eq (3.3), we have
From (3.6), we get |δ|≤1 , and hence the method is unconditionally stable for the discretized system of the Caputo time-fractional diffusion equation.
4.
Computational results and discussions
In this section, we consider numerical examples of the Caputo time-fractional diffusion equation in order to check the accuracy and efficiency of the method. We use the following error norms :
The ROC is calculated by
where, Eh1 and Eh2 represent the errors at space mesh sizes h1 and h2, respectively and Ek1 and Ek2 represent the errors at time mesh sizes k1 and k2, respectively.
4.1. Example 1
First, we consider the Caputo time-fractional diffusion equation
with the exact solution
where
We solve the Caputo time-fractional diffusion Eq (4.3) for different values of M and N. Table 1 shows the comparison of the present method with the finite element method [6] and cubic B-spline collocation method [24] for different values of M, fixed Δt=0.001 at the time-fractional order α=0.5 in terms of L2 and Lmax error norms. This table shows that the present method gives better results than the results obtained by those available in [6,24]. Moreover, one can notice from Table 1 that the error norms are decreasing as we increase the space as well as time mesh sizes. The analytical and approximate solutions together with the absolute errors are presented in Figure 1 for ˆκ=11, Δx=0.001, and α=0.5 at Δt=0.2, Δt=0.1, and Δt=0.05, respectively. Figure 2 shows the comparison between analytical and approximate solutions for M=64, Δt=0.001, and α=0.5 at different t. From these figures, we notice that there is an excellent agreement between analytical and approximate solutions. Also, the absolute errors are decreasing on increasing the time interval.
4.2. Example 2
Next, we consider the Caputo time-fractional diffusion equation
with the exact solution
where
We solve the Caputo time-fractional diffusion Eq (4.4) for various values of M. Table 2 shows the maximum absolute error norms Lmax for different values of time-fractional orders i. e. α=0.2, 0.5 and 0.8 at fixed time mesh size of N=50 for different space mesh sizes. One can see from Tables 2 and 3 that the error norms are decreasing as we increase the mesh sizes. As we can notice from this table that the present method is of order O((Δx)2,(Δt)1+α). The approximate solutions for time levels t= 0.1, 0.4, 0.6 and 0.8 are plotted in Figure 3 for M=32, N=50 and time-fractional order α=0.5. One can notice that the amplitude of the approximate solutions is increasing on incresing the time t.
4.3. Example 3
Now, we consider the Caputo time-fractional diffusion equation
where
subject to the initial and boundary conditions
and with the exact solution
Now, we solve the Caputo time-fractional diffusion Eq (4.5) for different values of M and N. Table 4 shows the comparison in term of Lmax error norms obtained by the present method and the method based on the product-integration (PI) rule presented in [7] for Δx=1101, β=6, different values of N and α at t=1. As we can see that the present method gives more accurate results than the method presented in [7]. Also, we notice that the error norms are decreasing as we increase the time mesh sizes. Figure 4 shows the analytical and approximate solutions together with absolute error norms for Δx=1101, α=0.85, β=5, Δt=0.01 at T=0.35 showing that errors are decreasing on increasing the grid sizes.
5.
Conclusions
A modified CExpB-spline collocation technique has been presented for solving the Caputo time-fractional diffusion equation. The modified CExpB-spline collocation technique is used to discretize the space derivatives. The three examples of the Caputo time-fractional diffusion equation have been considered. The obtained results show that the present method gives more accurate results than the results obtained in [6,7,24]. It is observed numerically that the method is second-order accurate in space and (1+α) order in time. The stability analysis shows that the method is unconditionally stable. Moreover, the implementation of the present method is easy and needs low memory storage which is the advantage. The present method can easily be extended for solving higher dimensional fractional PDEs.
Conflict of interest
The authors have no conflict of interest.