A Characteristic Difference Scheme for Time-Fractional Heat Equations Based on the Crank-Nicholson Difference Schemes

Abstract

We consider the numerical solution of a time-fractional heat equation, which is obtained from the standard diffusion equation by replacing the first-order time derivative with Riemann-Liouville fractional derivative of order α, where 0<𝛼<1. The main purpose of this work is to extend the idea on Crank-Nicholson method to the time-fractional heat equations. We prove that the proposed method is unconditionally stable, and the numerical solution converges to the exact one with the order 𝑂(𝜏²+ℎ²). Numerical experiments are carried out to support the theoretical claims.

1. Introduction

Fractional calculus is one of the most popular subjects in many scientific areas for decades. Many problems in applied science, physics and engineering are modeled mathematically by the fractional partial differential equations (FPDEs). We can see these models adoption in viscoelasticity [1, 2], finance [3, 4], hydrology [5, 6], engineering [7, 8], and control systems [9–11]. FPDEs may be investigated into two fundamental types: time-fractional differential equations and space-fractional differential equations.

Several different methods have been used for solving FPDEs. For the analytical solutions to problems, some methods have been proposed: the variational iteration method [12, 13], the Adomian decomposition method [13–16], as well as the Laplace transform and Fourier transform methods [17, 18].

On the other hand, numerical methods which based on a finite-difference approximation to the fractional derivative, for solving FDPEs [19–24], have been proposed. A practical numerical method for solving multidimensional fractional partial differential equations, using a variation on the classical alternating-directions implicit (ADI) Euler method, is presented in [25]. Many finite-difference approximations for the FPDEs are only first-order accurate. Some second-order accurate numerical approximations for the space-fractional differential equations were presented in [26–28]. Here, we propose a Crank-Nicholson-type method for time-fractional differential heat equations with the accuracy of order 𝑂(𝜏²+ℎ²).

In this work, we consider the following time-fractional heat equation: 𝜕^𝛼_𝑀𝑢(𝑡,𝑥)𝜕𝑡^𝛼=𝜕²𝑢(𝑡,𝑥)𝜕𝑥²(),(),𝑢()(𝑥)𝑢()()+𝑓𝑡,𝑥0<𝑥<1,0<𝑡<10,𝑥=𝑟,0≤𝑥≤1,𝑡,0=0,𝑢𝑡,1=0,0≤𝑡≤1. (1.1) Here, the term 𝜕^𝛼_𝑀𝑢(𝑡,𝑥)/𝜕𝑡^𝛼 denotes 𝛼-order-modified Riemann-Liouville fractional derivative [29] given with the formula: 𝜕^𝛼_𝑀𝑢(𝑡,𝑥)𝜕𝑡^𝛼=1Γ()𝜕1−𝛼∫𝜕𝑡^𝑡0𝑢(𝑠,𝑥)−𝑢(0,𝑥)(𝑡−𝑠)^𝛼𝜕𝑑𝑠,if0<𝛼<1,𝑢()𝜕𝑡𝑡,𝑥,if𝛼=1, (1.2) where Γ(⋅) is the Gamma function.

Remark 1.1. If 𝑟(𝑥)=0, then the Riemann-Liouville and the modified Riemann-Liouville fractional derivatives are identical, since the Riemann-Liouville derivative is given by the following formula: 𝜕^𝛼𝑢(𝑡,𝑥)𝜕𝑡^𝛼=1Γ()𝜕1−𝛼∫𝜕𝑡^𝑡0𝑢(𝑠,𝑥)(𝑡−𝑠)^𝛼𝜕𝑑𝑠,if0<𝛼<1,𝑢()𝜕𝑡𝑡,𝑥,if𝛼=1. (1.3) If 𝑟(𝑥) is nonzero, then there are some problems about the existence of the solutions for the heat equation (1.1). To rectify the situation, two main approaches can be used: the modified Riemann-Liouville fractional derivative can be used [29] or the initial condition should be modified [30]. We chose the first approach in our work.

2. Discretization of the Problem

In this section, we introduce the basic ideas for the numerical solution of the time-fractional heat equation (1.1) by Crank-Nicholson difference scheme.

For some positive integers 𝑀 and 𝑁, the grid sizes in space and time for the finite-difference algorithm are defined by ℎ=1/𝑀 and 𝜏=1/𝑁, respectively. The grid points in the space interval [0,1] are the numbers 𝑥_𝑖=𝑖ℎ, 𝑖=0,1,2,…,𝑀, and the grid points in the time interval [0,1] are labeled 𝑡_𝑛=𝑛𝜏, 𝑛=0,1,2,…,𝑁. The values of the functions 𝑈 and 𝑓 at the grid points are denoted 𝑈^𝑛_𝑖=𝑈(𝑡_𝑛,𝑥_𝑖) and 𝑓^𝑛_𝑖=𝑓(𝑡_𝑛,𝑥_𝑖), respectively.

As in the classical Crank-Nicholson difference scheme, we will obtain a discrete approximation to the fractional derivative 𝜕^𝛼𝑈(𝑡,𝑥)/𝜕𝑡^𝛼 at (𝑡_𝑛+(1/2),𝑥_𝑖). Let 𝐻()=1𝑡,𝑥Γ(1−𝛼)^𝑡₀𝑢(𝑠,𝑥)−𝑢(0,𝑥)(𝑡−𝑠)^𝛼𝑑𝑠. (2.1) Then, we have 𝜕^𝛼𝑈𝑡_𝑛+1/2,𝑥_𝑖𝜕𝑡^𝛼=𝜕𝐻𝑡𝜕𝑡^𝑛+1/2,𝑥_𝑖=𝐻𝑡_𝑛+1,𝑥_𝑖−𝐻𝑡_𝑛,𝑥_𝑖𝜏+𝑂𝜏². (2.2) Now, we will find the approximations for 𝐻(𝑡_𝑛+1,𝑥_𝑖) and 𝐻(𝑡_𝑛,𝑥_𝑖): 𝐻𝑡_𝑛+1,𝑥_𝑖=1Γ(1−𝛼)^𝑡𝑛+1₀𝑢𝑠,𝑥_𝑖−𝑢0,𝑥_𝑖𝑡_𝑛+1−𝑠^𝛼=1𝑑𝑠Γ(1−𝛼)^𝑛+1_𝑗=1^𝑗𝜏_{(𝑗−1)𝜏}𝑢𝑠,𝑥_𝑖𝑡_𝑛+1−𝑠^𝛼𝑑𝑠−𝑢0,𝑥_𝑖((𝑛+1)𝜏)^1−𝛼Γ()=12−𝛼Γ(1−𝛼)^𝑛+1_𝑗=1^𝑗𝜏_{(𝑗−1)𝜏}𝑠−𝑡_𝑗𝑈−𝜏^𝑖𝑗−1+𝑠−𝑡_𝑗−1𝜏𝑈^𝑗_𝑖+𝑂𝜏²1𝑡_𝑛+1−𝑠^𝛼𝑑𝑠−𝑈⁰_𝑖((𝑛+1)𝜏)^1−𝛼Γ(2−𝛼)=𝜏^𝑛_𝑗=0𝑎_𝑗−𝑗𝑏_𝑗𝑈_𝑖^𝑛−𝑗−𝜏^𝑛_𝑗=0𝑎_𝑗−𝑗+1𝑏_𝑗𝑈_𝑖^{𝑛−𝑗+1}−𝑈⁰_𝑖((𝑛+1)𝜏)^1−𝛼Γ(2−𝛼)+𝑅_𝑛+1, (2.3) where 𝑅_𝑛+1=1Γ(1−𝛼)^𝑛+1_𝑗=1^𝑗𝜏_{(𝑗−1)𝜏}𝑂𝜏²𝑑𝑠𝑡_𝑛+1−𝑠^𝛼=1()Γ()𝑂𝜏1−𝛼1−𝛼^2𝑛+1_𝑗=1𝑛−𝑗+2^1−𝛼−𝑛−𝑗+1^1−𝛼𝜏^1−𝛼=1Γ()()2−𝛼𝑛+1^1−𝛼𝑂𝜏^3−𝛼. (2.4)

Similarly, we can obtain 𝐻𝑡_𝑛,𝑥_𝑖=1Γ(1−𝛼)^𝑡𝑛₀𝑢𝑠,𝑥_𝑖−𝑢0,𝑥_𝑖𝑡_𝑛−𝑠^𝛼𝑑𝑠=𝜏^𝑛_𝑗=1𝑎_𝑗−1−𝑗−1𝑏_𝑗−1𝑈_𝑖^𝑛−𝑗−𝜏^𝑛_𝑗=1𝑎_𝑗−1−𝑗𝑏_𝑗−1𝑈_𝑖^{𝑛−𝑗+1}−𝑈⁰_𝑖(𝑛𝜏)^1−𝛼Γ(2−𝛼)+𝑅_𝑛, (2.5) where 𝑅_𝑛=(1/Γ(2−𝛼))𝑛^1−𝛼𝑂(𝜏^3−𝛼) and 𝑎_𝑗=𝜏^−𝛼(2−𝛼)Γ(1−𝛼)𝑗+1^2−𝛼−𝑗^2−𝛼,𝑏_𝑗=𝜏^−𝛼(1−𝛼)Γ(1−𝛼)𝑗+1^1−𝛼−𝑗^1−𝛼. (2.6)

Then, we can write the following approximation: 𝜕^𝛼𝑈𝑡_𝑛+1/2,𝑥_𝑖𝜕𝑡^𝛼=𝐻𝑡_𝑛+1,𝑥_𝑖−𝐻𝑡_𝑛,𝑥_𝑖𝜏+𝑂𝜏²=𝑞_𝑛𝑈⁰_𝑖+^𝑛_𝑗=0𝑝_𝑗𝑈_𝑖^{𝑛+1−𝑗}+𝑅^𝑛+1−𝑅^𝑛𝜏+𝑂𝜏²=𝑞_𝑛𝑈⁰_𝑖+^𝑛_𝑗=0𝑝_𝑗𝑈_𝑖^{𝑛+1−𝑗}+1Γ())2−𝛼(𝑛+1^1−𝛼−𝑛^1−𝛼𝑂𝜏^2−𝛼+𝑂𝜏²=𝑞_𝑛𝑈⁰_𝑖+^𝑛_𝑗=0𝑝_𝑗𝑈_𝑖^{𝑛+1−𝑗}+1Γ()()2−𝛼𝑛+1^1−𝛼−𝑛^1−𝛼𝜏𝑂𝜏^3−𝛼+𝑂𝜏²=𝑞_𝑛𝑈⁰_𝑖+^𝑛_𝑗=0𝑝_𝑗𝑈_𝑖^{𝑛+1−𝑗}+1Γ()(𝜏(2−𝛼𝑛+1))^1−𝛼−(𝜏𝑛)^1−𝛼𝜏𝑂𝜏²+𝑂𝜏², (2.7) where 𝑞₀=3𝑎₀−𝑎₁+2𝑏₁−2𝑏₀,𝑞_𝑛=𝑎_𝑛−𝑎_𝑛−1+(𝑛−1)𝑏_𝑛−1−(𝑛+1)𝑏_𝑛𝑝,for1≤𝑛≤𝑁−1,₀=𝑏₀−𝑎₀,𝑝₁=2𝑎₀−𝑎₁+2𝑏₁−𝑏₀,𝑝_𝑗=−𝑎_𝑗−2+2𝑎_𝑗−1−𝑎_𝑗+𝑗−2𝑏_𝑗−2−2𝑗−1𝑏_𝑗−1+𝑗+1𝑏_𝑗,for𝑗≥2. (2.8)

On the other hand, using the mean-value theorem, we get (𝜏(𝑛+1))^1−𝛼−(𝜏𝑛)^1−𝛼𝜏=𝑓^′(𝑐)=constant, (2.9) where 𝑓(𝑥)=𝑥^1−𝛼 and 𝑡_𝑛<𝑐<𝑡_𝑛+1. So, we obtain the following second-order approximation for the modified Riemann-Liouville derivative: 𝜕^𝛼𝑈𝑡_𝑛+(1/2),𝑥_𝑖𝜕𝑡^𝛼=𝐻𝑡_𝑛+1,𝑥_𝑖−𝐻𝑡_𝑛,𝑥_𝑖𝜏+𝑂𝜏²=𝑞_𝑛𝑈⁰_𝑖+^𝑛_𝑗=0𝑝_𝑗𝑈_𝑖^{𝑛+1−𝑗}+𝑂𝜏². (2.10)

3. Crank-Nicholson Difference Scheme

Using the approximation above, we obtain the following difference scheme which is accurate of order 𝑂(𝜏²+ℎ²): 𝑞_𝑛𝑈⁰_𝑖+^𝑛_𝑗=0𝑝_𝑗𝑈_𝑖^{𝑛+1−𝑗}−𝑈^𝑛+1_𝑖+1−2𝑈_𝑖^𝑛+1+𝑈^𝑛+1_𝑖−12ℎ²+𝑈^𝑛_𝑖+1−2𝑈^𝑛_𝑖+𝑈^𝑛_𝑖−12ℎ²=𝑓𝑡_𝑛+𝜏2,𝑥_𝑖𝑈,0≤𝑛≤𝑁−1,1≤𝑖≤𝑀−1,_0𝑖=𝑟𝑥_𝑖𝑈,1≤𝑖≤𝑀−1,_𝑛0=0,𝑈^𝑛_𝑀=0,0≤𝑛≤𝑁. (3.1)

We can arrange the system above to obtain −12ℎ²𝑈^𝑛+1_𝑖+1+𝑈^𝑛_𝑖+1+𝑞_𝑛𝑈⁰_𝑖+^𝑛_𝑗=0𝑝_𝑗𝑈_𝑖^{𝑛+1−𝑗}+−12ℎ²𝑈^𝑛+1_𝑖−1+𝑈^𝑛_𝑖−1=𝑓𝑡_𝑛+𝜏2,𝑥_𝑖𝑈,0≤𝑛≤𝑁−1,1≤𝑖≤𝑀−1,_0𝑖=𝑟𝑥_𝑖𝑈,1≤𝑖≤𝑀−1,_𝑛0=0,𝑈^𝑛_𝑀=0,0≤𝑛≤𝑁. (3.2)

The difference scheme above can be written in matrix form: 𝐴𝑈_𝑖+1+𝐵𝑈_𝑖+𝐴𝑈_𝑖−1=𝜑_𝑖, (3.3) where 𝜑_𝑖=[𝜑⁰_𝑖,𝜑¹_𝑖,𝜑²_𝑖,…,𝜑^𝑁_𝑖]^𝑇, 𝜑⁰_𝑖=𝑟(𝑥_𝑖), 𝜑^𝑛_𝑖=𝑓(𝑡_𝑛+1/2,𝑥_𝑖), 1≤𝑛≤𝑁, 1≤𝑖≤𝑀, and 𝑈_𝑖=[𝑈⁰_𝑖,𝑈¹_𝑖,𝑈²_𝑖,…,𝑈^𝑁_𝑖]^𝑇.

Here, 𝐴_{(𝑁+1)×(𝑁+1)} and 𝐵_{(𝑁+1)×(𝑁+1)} are the matrices of the form −1𝐴=2ℎ²0,1𝑞1111⋱⋱11𝐵=⁰+1ℎ²𝑝₀+1ℎ²𝑞₁𝑝₁+1ℎ²𝑝₀+1ℎ²𝑞₂𝑝₂𝑝₁+1ℎ²𝑝₀+1ℎ²𝑞⋮⋱⋱⋱_𝑁−1𝑝_𝑁−1…𝑝₂𝑝₁+1ℎ²𝑝₀+1ℎ². (3.4)

We note that the unspecified entries are zero at the matrices above.

Using the idea on the modified Gauss-Elimination method, we can convert (3.3) into the following form: 𝑈_𝑖=𝛼_𝑖+1𝑈_𝑖+1+𝛽_𝑖+1,𝑖=𝑀−1,…,2,1,0. (3.5)

This way, the two-step form of difference schemes in (3.3) is transformed to one-step method as in (3.5).

Now, we need to determine the matrices 𝛼_𝑖+1 and 𝛽_𝑖+1 satisfying the last equality. Since 𝑈₀=𝛼₁𝑈₁+𝛽₁=0, we can select 𝛼₁=𝑂_{(𝑁+1)×(𝑁+1)} and 𝛽₁=𝑂_(𝑁+1)×1. Combining the equalities 𝑈_𝑖=𝛼_𝑖+1𝑈_𝑖+1+𝛽_𝑖+1 and 𝑈_𝑖−1=𝛼_𝑖𝑈_𝑖+𝛽_𝑖 and the matrix equation (3.3), we have 𝐴+𝐵𝛼_𝑖+1+𝐴𝛼_𝑖𝛼_𝑖+1𝑈_𝑖+1+𝐵𝛽_𝑖+1+𝐴𝛼_𝑖𝛽_𝑖+1+𝐴𝛽_𝑖=𝜑_𝑖. (3.6) Then, we write 𝐴+𝐵𝛼_𝑖+1+𝐴𝛼_𝑖𝛼_𝑖+1=0,𝐵𝛽_𝑖+1+𝐴𝛼_𝑖𝛽_𝑖+1+𝐴𝛽_𝑖=𝜑_𝑖, (3.7) where  1≤𝑖≤𝑀−1.

So, we obtain the following pair of formulas: 𝛼_𝑖+1=−𝐵+𝐴𝛼_𝑖⁻¹𝛽𝐴,_𝑖+1=𝐵+𝐴𝛼_𝑖⁻¹𝜑_𝑖−𝐴𝛽_𝑖, (3.8) where  1≤𝑖≤𝑀−1.

4. Stability of the Method

The stability analysis is done by using the analysis of the eigenvalues of the iteration matrix 𝛼_𝑖 (1≤𝑖≤𝑀) of the scheme (3.5).

Let  𝜌(𝐴) denote the spectral radius of a matrix 𝐴, that is, the maximum of the absolute value of the eigenvalues of the matrix 𝐴.

We will prove that 𝜌(𝛼_𝑖)<1, (1≤𝑖≤𝑀), by induction.

Since 𝛼₁ is a zero matrix 𝜌(𝛼₁)=0<1.

Moreover, 𝛼₂=−𝐵⁻¹𝐴,𝜌(𝛼₂)=𝜌(−𝐵⁻¹𝐴)=−11/ℎ²+𝑝₀⋅−12ℎ²=1/ℎ²2(1/ℎ²+𝑝₀), since 𝛼₂ is of the form 𝛼₂=0∗1/ℎ²21/ℎ²+𝑝₀∗∗1/ℎ²21/ℎ²+𝑝₀⋱∗∗∗1/ℎ²21/ℎ²+𝑝_{0(𝑁+1)×(𝑁+1)},𝑝₀=𝑏₀−𝑎₀=𝜏^−𝛼()Γ()−𝜏1−𝛼1−𝛼()Γ^−𝛼()=𝜏2−𝛼1−𝛼Γ(3−^−𝛼𝛼)>0, (4.1) therefore, 𝜌(𝛼₂)<1.

Now, assume 𝜌(𝛼_𝑖)<1. After some calculations, we find that 𝛼_𝑖+1=−𝐵+𝐴𝛼_𝑖⁻¹𝐴=12ℎ²0∗1𝐵_2,2−1/2ℎ²𝛼_𝑖2,21∗∗𝐵_3,3−1/2ℎ²𝛼_𝑖3,31∗∗∗⋱∗∗∗𝐵_{𝑁+1,𝑁+1}−1/2ℎ²𝛼_{𝑖𝑁+1,𝑁+1} (4.2) and we already know that 𝐵_𝑗,𝑗=1/ℎ²+𝑤₀ and 𝛼_{𝑖𝑗,𝑗}=𝜌(𝛼_𝑖) for 2≤𝑗≤𝑁+1: 𝜌𝛼_𝑖+1=1/2ℎ²1/ℎ²+𝑝₀−1/2ℎ²𝜌𝛼_𝑖=𝑀²2𝑀²1−𝜌𝛼_𝑖/2+𝑝₀. (4.3)

Since 0≤𝜌(𝛼_𝑖)<1, it follows that 𝜌(𝛼_𝑖+1)<1. So, 𝜌(𝛼_𝑖)<1 for any 𝑖, where 1≤𝑖≤𝑀.

Remark 4.1. The convergence of the method follows from the Lax equivalence theorem [31] because of the stability and consistency of the proposed scheme.

5. Numerical Analysis

Example 5.1. Consider 𝜕^𝛼𝑢(𝑡,𝑥)𝜕𝑡^𝛼=𝜕²𝑢(𝑡,𝑥)𝜕𝑥²+2𝑡^(2−𝛼)Γ(3−𝛼)sin𝑥−𝑥²+𝑡²sin𝑥−𝑥²(1−2𝑥)²+2cos𝑥−𝑥²,(),𝑢()𝑢()()0<𝑥<1,0<𝑡<10,𝑥=0,0≤𝑥≤1,𝑡,0=0,𝑢𝑡,1=0,0≤𝑡≤1. (5.1)
Exact solution of this problem is 𝑈(𝑡,𝑥)=𝑡²sin(1−𝑥)𝑥. The solution by the Crank-Nicholson scheme is given in Figure 1. The errors when solving this problem are listed in the Table 1 for various values of time and space nodes.
The errors in the table are calculated by the formula max_{0≤𝑛≤𝑀,0≤𝑘≤𝑁}|𝑢(𝑡_𝑘,𝑥_𝑛)−𝑈^𝑘_𝑛| and the error rate formula is |𝐸_𝑘|/|𝐸_𝑘+1|.

Table 1 

Error table for Example 5.1.

(a)

(b)

Figure 1 

(a) The approximate solutions of Example 5.1 by the proposed method when 𝑁=32,  𝑀=32, and 𝛼=0.5. (b) The errors for some values of 𝑀 and 𝑁 when 𝑡=1 and 𝛼=0.5.

Example 5.2. Consider 𝜕^𝛼𝑢(𝑡,𝑥)𝜕𝑡^𝛼=𝜕²𝑢(𝑡,𝑥)𝜕𝑥²+24𝑡^(2−𝛼)Γ()𝑥5−𝛼²−𝑥−2𝑡⁴,(),𝑢()𝑢()()0<𝑥<1,0<𝑡<10,𝑥=0,0≤𝑥≤1,𝑡,0=0,𝑢𝑡,1=0,0≤𝑡≤1. (5.2)
Exact solution of this problem is 𝑈(𝑡,𝑥)=𝑡⁴𝑥(𝑥−1). The solution by the Crank-Nicholson scheme is given in Figure 2. The errors when solving this problem are listed in Table 2 for various values of time and space nodes and several values of 𝛼.
It can be concluded from the tables and the figures that when the step size is reduced by a factor of 1/2, the error decreases by about 1/4. The numerical results support the claim about the order of the convergence.

Table 2 

The errors for some values of 𝑀, 𝑁, and 𝛼.

(a)

(b)

Figure 2 

(a) The approximate solutions of Example 5.2 by the proposed method when 𝑁=32,  𝑀=32, and 𝛼=0.5. (b) The errors for some values of 𝑀 and 𝑁 when 𝑡=1 and 𝛼=0.5.

6. Conclusion

In this work, the Crank-Nicholson difference scheme was successfully extended to solve the time-fractional heat equations. A second-order approximation for the Riemann-Liouville fractional derivative is obtained. It is proven that the time-fractional Crank-Nicholson difference scheme is unconditionally stable and convergent. Numerical results are in good agreement with the theoretical results.