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

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] 𝜕𝑢()𝜕𝑡𝑥,𝑡=𝐾₀𝐷_𝑡^1−𝛾𝜕²𝜕𝑥²𝑢(𝑥,𝑡), (1) where₀𝐷_𝑡^1−𝛾𝑓(𝑡) is the fractional derivative defined by Riemann-Liouville operator as ₀𝐷_𝑡^1−𝛾𝑓(𝑡)=1Γ𝛾𝜕𝜕𝑡^𝑡₀𝑓(𝜏)(𝑡−𝜏)^1−𝛾𝑑𝜏, (2) 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 𝑢(0,𝑡)=0,𝑢(1,𝑡)=0 (3) and the initial condition as 𝑢(𝑥,0)=𝑥(1−𝑥). (4) The exact analytical solution of (1) is found by the method of separation of variables [9] as 𝑢()=8𝑥,𝑡𝜋^3∞_𝑛=01(2𝑛+1)³sin[(2𝑛+1)𝜋𝑥]×𝐸_𝛾−𝐾(2𝑛+1)²𝜋²𝑡^𝛾, (5) 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 (0<𝛾≤1), we will also discretize the Riemann-Liouville operator [9, 16] as ₀𝐷_𝑡^1−𝛾𝑢𝑥_𝑗,𝑡_𝑛=₀𝐷_𝑡^1−𝛾𝑢^𝑛_𝑗=1(Δ𝑡)^𝑛1−𝛾_𝑘=0𝜔_𝑘^1−𝛾𝑢_𝑗^𝑛−𝑘+𝑂Δ𝑡^𝑝, (6) where 𝜔₀^1−𝛾=1,𝜔_𝑘^1−𝛾=1−2−𝛾𝑘𝜔^1−𝛾_𝑘−1. (7)

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 𝑥_𝑚, 𝑚=0,1,2,…,𝑁 such that 𝑎=𝑥₀<𝑥₁<⋯<𝑥_𝑁=𝑏 and ℎ=𝑥_𝑚+1−𝑥_𝑚. The cubic B-splines 𝜙_𝑚(𝑥), (𝑚=−1(1)𝑁+1), over the interval [𝑎,𝑏] and at the knots 𝑥_𝑚 are defined by [17]𝜙_𝑚(𝑥)=1ℎ³𝑥−𝑥_𝑚−2³,𝑥∈𝑥_𝑚−2,𝑥_𝑚−1,ℎ³+3ℎ²𝑥−𝑥_𝑚−1+3ℎ𝑥−𝑥_𝑚−1²−3𝑥−𝑥_𝑚−1³,𝑥∈𝑥_𝑚−1,𝑥_𝑚,ℎ³+3ℎ²𝑥_𝑚+1−𝑥+3ℎ𝑥_𝑚+1−𝑥²−3𝑥_𝑚+1−𝑥³,𝑥∈𝑥_𝑚,𝑥_𝑚+1,𝑥_𝑚+2−𝑥³,𝑥∈𝑥_𝑚+1,𝑥_𝑚+2,0otherwise. (8) The set of cubic splines {𝜙₋₁(𝑥),𝜙₀(𝑥),…,𝜙_𝑁+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 𝑈_𝑁(𝑥,𝑡)=^𝑁+1_𝑚=−1𝛿_𝑚(𝑡)𝜙_𝑚(𝑥), (9) 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 𝛿_𝑚(𝑡)𝑈_𝑚=𝑈𝑥_𝑚=𝛿_𝑚−1(𝑡)+4𝛿_𝑚(𝑡)+𝛿_𝑚+1(𝑡),𝑈^′_𝑚=𝑈^′𝑥_𝑚=3ℎ−𝛿_𝑚−1(𝑡)+𝛿_𝑚+1(𝑡),𝑈_𝑚^′′=𝑈^′′𝑥_𝑚=6ℎ²𝛿_𝑚−1(𝑡)−2𝛿_𝑚(𝑡)+𝛿_𝑚+1(𝑡), (10) the variation of 𝑈_𝑁(𝑥,𝑡) over the typical element [𝑥_𝑚,𝑥_𝑚+1] can now be given by 𝑈_𝑁(𝑥,𝑡)=^𝑚+2_{𝑗=𝑚−1}𝛿_𝑗(𝑡)𝜙_𝑗(𝑥). (11) 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: 𝛿̇_𝑚−1(𝑡)+4𝛿̇_𝑚(𝑡)+𝛿̇_𝑚+1(𝑡)−6ℎ²⁰𝐷_𝑡^1−𝛾𝛿_𝑚−1(𝑡)−2𝛿_𝑚(𝑡)+𝛿_𝑚+1(𝑡)=0, (12) 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:1𝛿=2𝛿^𝑛+𝛿^𝑛+1,𝛿̇=𝛿^𝑛+1−𝛿^𝑛.Δ𝑡 (13) Then, if we discretize the fractional operator₀𝐷_𝑡^1−𝛾𝛿^𝑛 by the following formula: ₀𝐷_𝑡^1−𝛾𝛿^𝑛=1Δ𝑡^𝑛1−𝛾_𝑘=0𝜔_𝑘^1−𝛾𝛿^𝑛−𝑘, (14) we easily obtain a recurrence relationship between successive time levels relating unknown parameters 𝛿_𝑚^𝑛+1(𝑡) as (1−𝛼)𝛿^𝑛+1_𝑚−1+(4+2𝛼)𝛿_𝑚^𝑛+1+(1−𝛼)𝛿^𝑛+1_𝑚+1=(1+𝛼)𝛿^𝑛_𝑚−1+(4−2𝛼)𝛿^𝑛_𝑚+(1+𝛼)𝛿^𝑛_𝑚+1+𝛼^𝑛_𝑘=1𝜔_𝑘^1−𝛾𝛿^{𝑛+1−𝑘}_𝑚−1+𝛿^𝑛−𝑘_𝑚−1−2𝛿_𝑚^{𝑛+1−𝑘}+𝛿_𝑚^𝑛−𝑘+𝛿^{𝑛+1−𝑘}_𝑚+1+𝛿^𝑛−𝑘_𝑚+1, (15) where 3()𝛼=Δ𝑡^𝛾ℎ². (16) The system (15) is consisted of 𝑁+1 linear equations including 𝑁+3 unknown parameters (𝛿₋₁,…,𝛿_𝑁+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 𝛿₋₁ and 𝛿_𝑁+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 vector 𝐝⁰=(𝛿₀,𝛿₁,𝛿₂,…,𝛿_𝑁−2,𝛿_𝑁−1,𝛿_𝑁)^𝑇 is determined from the initial and boundary conditions. Thus, the approximation (11) can be rewritten for the initial condition as 𝑈_𝑁(𝑥,0)=^𝑁+1_𝑚=−1𝛿_𝑚(0)𝜙_𝑚(𝑥), (17) where the 𝛿_𝑚(0) are unknown element parameters. We force the initial numerical approximation 𝑈_𝑁(𝑥,0) to meet the following conditions: 𝑈_𝑁(𝑥,0)=𝑈𝑥_𝑚𝑈,0,𝑚=0,1,…,𝑁,_𝑁𝑥𝑥(0,0)=−2,𝑈_𝑁𝑥𝑥(1,0)=−2. (18) Using these conditions results in a three-diagonal system of matrix of the form 𝐖𝐝⁰=𝐛, (19) where ⋱,𝑈𝑥𝐖=6014114114106𝐛=₀+ℎ,0²3,𝑈𝑥₁,0,𝑈𝑥₂,0,…,𝑈𝑥_𝑁−2,𝑈𝑥,0_𝑁−1,0,𝑈(𝑥_𝑁ℎ,0)+²3^𝑇. (20)

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

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) where 𝑖=√−1. First of all, substituting the Fourier mode (21) into the recurrence relationship (15) results in the following equation: 𝜉^𝑛+1(1−𝛼)𝑒^−𝑖𝜑+(4+2𝛼)+(1−𝛼)𝑒^𝑖𝜑=𝜉^𝑛(1+𝛼)𝑒^−𝑖𝜑+(4−2𝛼)+(1+𝛼)𝑒^𝑖𝜑+𝛼^𝑛_𝑘=1𝜔_𝑘^1−𝛾𝜉^{𝑛−𝑘+1}+𝜉^𝑛−𝑘𝑒^−𝑖𝜑−2+𝑒^𝑖𝜑. (22) Secondly, if we write 𝜉^𝑛+1=𝜁𝜉^𝑛 (23) and assume that 𝜁≡𝜁𝜑 (24) is independent of time, then we get the following expression for the amplification factor 𝜁 of the subdiffusion mode: 𝜁(1−𝛼)𝑒^−𝑖𝜑+(4+2𝛼)+(1−𝛼)𝑒^𝑖𝜑=(1+𝛼)𝑒^−𝑖𝜑+(4−2𝛼)+(1+𝛼)𝑒^𝑖𝜑+𝛼^𝑛_𝑘=1𝜔_𝑘^1−𝛾𝜁^1−𝑘+𝜁^−𝑘𝑒^−𝑖𝜑−2+𝑒^𝑖𝜑. (25) 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: 8𝛼sin²𝜑2≥0. (26) 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 norm 𝐿₂ as 𝐿₂=𝑈^exact−𝑈_𝑁2=ℎ^𝑁_𝑗=0𝑈^exact_𝑗−𝑈_𝑁𝑗² (27) and the error norm 𝐿_∞ as 𝐿_∞=𝑈^exact−𝑈_𝑁∞=max_𝑗𝑈^exact_𝑗−𝑈_𝑁𝑗. (28)

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 𝐿₂ 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 𝐿₂ 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 𝐿₂ 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 𝐿₂ and 𝐿_∞.

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 𝐿₂ 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 𝐿₂ 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 𝐿₂ 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 𝐿₂ and 𝐿_∞.

Figure 1 

The comparison of the exact (lines) and numerical solutions for 𝛾=0.50, Δ𝑡=0.0001, and 𝑁=40 at 𝑡=0.001 (triangles), 𝑡=0.01 (squares), and 𝑡=0.1 (stars).

Figure 2 

The comparison of the exact and numerical solutions for 𝛾=0.75, Δ𝑡=0.0001, and 𝑁=40 at at 𝑡=0.001 (triangles), 𝑡=0.01 (squares), and 𝑡=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 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.