1.
Introduction
In this paper, the following one-sided tempered fractional diffusion equations are considered:
and
where 1<α<2, λ≥0, the diffusion coefficient K is positive, and f(x,t) is the source term. Dα,λa,xu(x,t) and Dα,λx,bu(x,t) represent the left and right Riemann-Liouville tempered fractional derivatives, respectively, and are defined as
The concept of fractional derivatives appeared almost at the same time as integer derivatives, but the lack of application background of fractional derivatives at the beginning of their appearance made fractional models not widely developed until recent decades. Fractional models are widely used in physics [1,2,3,4], finance [5], biology [6], and hydrology [7,8,9]. Recently, the study of fractional diffusion equations has attracted a lot of attention, in which the integral derivatives of diffusion equations are replaced by the fractional derivatives to obtain fractional diffusion equations. Usually the time fractional derivatives describe anomalous sub-diffusion, and the space fractional derivatives describe anomalous super-diffusion [10].
The analytical solutions of fractional differential equations usually cannot be obtained because of the nonlocality of fractional derivatives. Therefore, a lot of attention has been paid to the development of high-precision numerical methods, and a lot of results have been obtained [11,12,13,14,15]. For the space fractional advection-dispersion equations, Meerschaert and Tadjeran [11] point out that the standard Grünwald difference operator to approximate the Riemann-Liouville fractional derivative is unconditionally unstable regardless of the implicit and explicit Euler methods, so a modified Grünwald difference operator, called the shifted Grünwald difference operator, is proposed to solve this problem. Based on the shift Grünwald difference operator and its idea, more research has been done [16,17,18,19,20]. While developing numerical methods for fractional diffusion equations, some researchers have found that if the power law (waiting time or jump length) is tempered by an exponential factor, it has practical advantages [21,22,23]. Therefore, the time and space tempered fractional derivatives are obtained, and the time tempered fractional derivative yields the time tempered fractional diffusion equation [22]. Luo et al. [24] proposed an effective Lagrange-quadratic spline optimal collocation method to solve it. The space tempered fractional derivative yields the space tempered fractional diffusion equation [23]. It is very important to develop numerical methods for tempered fractional diffusion equations.
A common way to solve the space tempered fractional diffusion equation is to develop a modified Grünwald difference operator based on the idea of the shifted Grünwald difference operator, called the tempered shifted Grünwald difference operator by Baeumer and Meerschaert [25], which is used to solve the tempered fractional advection-dispersion equation, but the space direction has only first-order precision. Based on the idea of the tempered shifted Grünwald difference operator, Li and Deng [26] obtained a class of second-order approximations of left and right Riemann-Liouville tempered fractional derivatives by using the weighted idea, called the tempered weighted and shifted Grünwald difference operators, which combined with the Crank-Nicolson method is used to solve the two-sided tempered fractional diffusion equations. See more studies on space-time tempered fractional differential equations [27,28], time tempered fractional differential equations [29,30], and space tempered fractional differential equations [31,32,33,34,35]. Similar to the quasi-compact scheme of fractional diffusion equations [36,37,38], it is natural to think of applying the quasi-compact technique to the numerical solutions of two-sided tempered fractional diffusion equations. However, due to the incompatibility of the left and right Riemann-Liouville tempered fractional quasi-compact operators, Yu et al. [39] considered the quasi-compact technique on the one-sided tempered fractional diffusion equations and obtained that the numerical schemes are stable and convergent with order O(τ+h3) and the schemes don't depend on time and space steps. The high-order quasi-compact scheme of the two-sided tempered fractional diffusion equation needs to be studied.
The novelty of this paper is that the fourth-order numerical schemes of the one-sided tempered fractional diffusion equations are given, and the convergence accuracy is one order higher than that of the existing literature. In the process of proposing the numerical schemes, the left and right third-order Riemann-Liouville tempered derivatives will be generated, so it is necessary to study their effective second-order numerical approximations, and we obtain them through the definition of Riemann-Liouville tempered integer derivatives. Finally, the effective fourth-order numerical schemes for solving these two classes of problems are obtained. In addition, there is a simpler numerical scheme that will be available to solve the two-sided fractional advection-diffusion equations [38] by using the idea of this paper.
The remaining sections of this paper are arranged as follows: Section 2 gives the fourth-order quasi-compact approximations of tempered fractional derivatives. In Section 3, the numerical scheme for solving problems (1.1) and (1.2) are derived. The stability and convergence of the numerical schemes are proved in Section 4. Section 5 shows by examples that the numerical schemes are effective. Finally, a brief summary of this work is given in Section 6.
2.
Fourth-order quasi-compact approximations of tempered fractional derivatives
Define a fractional Sobolev space Sn+αλ(R),
where ˆν(w)=∫Rν(x)e−iwxdx is the Fourier transform of ν(x).
Lemma 2.1. [25,26] Let 1<α<2, λ≥0. The shift number p is an integer, h is the step size, ν(x) is defined on the bounded interval [a,b], and belongs to Sn+αλ(R) after zero extension on the interval x∈(−∞,a)∪(b,+∞). The tempered and shifted Grünwald type difference operators are defined as
then
where g(α)k=(−1)k(αk)(k≥0) denotes the normalized Grünwald weights, that is,
cα,pk are the power series expansion coefficients of the function Wp(s)=eps(1−e−s)αs=+∞∑k=0cα,pksk, and the first four coefficients are given as
Remark 2.1. Let the equations be
then
and there are the following approximations:
where
in particular, when α=1, the following approximations are given as
For the first left and right Riemann-Liouville tempered derivatives D1,λa,xν(x) and D1,λx,bν(x), we know
so the normalized left and right Riemann-Liouville tempered fractional derivatives can be rewritten as
and
Let Λαx:=(I+ch2D2,λa,x) and Δαx:=(I+ch2D2,λx,b), where c=−α2+α+424∈(112,16]. Now, we give the following theorem as the contribution of this paper.
Theorem 2.1. Let ν(x)∈S4+αλ(R),1<α<2, λ≥0, and the continuous operators Λαx and Δαx are given to operate on the normalized left and right Riemann-Liouville tempered fractional derivatives, respectively. Then,
and
have a fourth-order approximation, respectively. The details are as follows:
and
where
Proof. From Remark 2.1, the conclusion of (i) and (ii) in Eqs (2.12) and (2.13) can be obtained.
For the conclusion (iii) in Eqs (2.12) and (2.13), it is easy to obtain the following:
Finally, D2,λa,x(D1,λa,xν(x)) and D2,λx,b(D1,λx,bν(x)) can be represented as
For derivatives of order 1 to 3, we adopt the following discretizations, respectively,
Thus, the whole conclusion of the theorem is proved. □
3.
Derivation of the numerical schemes
In this section, we consider the numerical scheme for solving problems (1.1) and (1.2). Here, we always assume that the function u(x,⋅) in problem (1.1) and (1.2) belongs to S4+αλ(R) after zero extension.
We make space grid {xi=a+ih}M0 and time grid {tn=nτ}N0, where h=b−aM and τ=TN represent the space step size and the time step size, respectively. Let uni=u(xi,tn), fni=f(xi,tn), and Uni represent the numerical solution at the point (xi,tn).
For problem (1.1), the operator Λαx is applied to both the left and right ends of equation. We obtain
then the backward Euler method is carried out at the point (xi,tn) to discrete the time partial derivative, and by Theorem 2.1, we get
where Rni=O(τ+h4) is the local truncation error.
Removing the local truncation error in Eq (3.2) yields the numerical scheme as
the matrix form of the numerical scheme (3.3) can be written as
where P=B(α)+(α−1)λαC(α)−αλα−1B(1)+(16−c)αλα+2h2E, Q=(16−c)αλα−1h2(3λ2D(1)+3λD(2)+D(3)), Un=(Un1,Un2,...,UnM−2,UnM−1)T, fn=(fn1,fn2,...,fnM−2,fnM−1)T, E is the identity matrix of order M−1,
and
The same idea is used to solve problem (1.2). The operator Δαx is applied to both the left and right ends of the equation, and we obtain
then the backward Euler method is carried out at the point (xi,tn) to discrete the time partial derivative. By Theorem 2.1, we get
where ˆRni=O(τ+h4) is the local truncation error.
Removing the local truncation error in Eq (3.14) yields the numerical scheme as
the matrix form of the numerical scheme (3.15) can be written as
where ˆP=ˆB(α)+(α−1)λαˆC(α)−αλα−1ˆB(1)+(16−c)αλα+2h2E, ˆQ=(16−c)αλα−1h2(−3λ2D(1)+3λD(2)−ˆD(3)), ˆB(α)=(B(α))T, ˆB(1)=(B(1))T, ˆC(α)=(C(α))T, ˆP=PT, ˆD(3)=−(D(3))T, ˆQ=QT,
4.
Stability and convergence analysis of the numerical schemes
In this section, we use the energy method to prove in detail that numerical schemes (3.3) and (3.15) are valid. Before doing so, we first introduce some lemmas that will be used.
Lemma 4.1. [40] A real matrix A of order M is positive definite if and only if D=A+AT2 is positive definite.
Lemma 4.2. [41] Let T be a Toeplitz matrix with the generating function f(x) being a 2π-periodic continuous real-valued function. Denote λmin(T) and λmax(T) as the smallest and largest eigenvalues of T, respectively. Then, we have
where fmin(x) and fmax(x) denote the minimum and maximum values of f(x), respectively. In particular, if f(x) is a nonpositive function and is not always zero, and fmin(x)≠fmax(x),
Lemma 4.3. [41] (Weyl's theorem). Let A,E∈Cn×n be a Hermitian matrix and the eigenvalues λi(A),λi(E),λi(A+E) be arranged in an increasing order. Then, for each k = 1, 2, ..., n, we have
Lemma 4.4. For 1<α<2,0<λh≤1, then the matrices P, Q, and C(α) in Eq (3.4) have the following properties:
and the matrix P is negative definite.
(ii). The matrix Q is negative definite.
(iii). For all given nonzero real column vectors ϵ, C(α) satisfies that
Proof. (i). By simple calculation, it is easy to obtain the properties of the elements of the matrix P; From the Gerschgorin disk theorem [42], we get that the matrix P is negative definite.
(ii). Note that Q=(16−c)αλα−1h2(3λ2D(1)+3λD(2)+D(3)), so we just consider 3λ2D(1)+3λD(2)+D(3). First, it is easy to show that the matrix 3λ2D(1)+3λD(2) is negative definite; and then for the matrix D(3), let D=D(3)+(D(3))T2, which is a Toeplitz matrix, and the generating function [41] is
From Lemma 4.2, we can see that the matrix D is negative definite. Further from Lemma 4.1, the matrix Q is negative definite.
(iii). Let H=C(α)+(C(α))T2, and the generating function of the matrix H is
From Lemma 4.2, we obtain 1−c(eλh+e−λh+2)≤λ(H)≤1+c(eλh+e−λh−2), and it is easy to check 112<λ(H)<43, which means 112ϵTϵ<ϵTHϵ<43ϵTϵ, that is, 112ϵTϵ<ϵTC(α)ϵ<43ϵTϵ.
Thus, all the conclusions of the lemma are proved. □
Lemma 4.5. [43] Assume that {kn} and {pn} are nonnegative sequences, and the sequence {ϕn} satisfies
where g0≥0. Then, the sequence {ϕn} satisfies
To prove the stability and convergence of numerical schemes by the energy method, we define Uh={u|u={ui} as a grid function defined on {xi=a+ih}M−1i=1}, for u∈Uh, and the corresponding discrete L2-norm is defined as ‖u‖L2=(hM−1∑i=1u2i)1/2. Next, we present the theoretical analysis.
Theorem 4.1. For α∈(1,2), let 0<λh≤1, then the numerical scheme (3.3) is stable for solving problem (1.1).
Proof. Let Uni and Vni represent the solution obtained by solving problem (1.1) using scheme (3.3) from different initial values. Denoting εni=Uni−Vni,εn=(εn1,εn2,...,εnM−1)T, it can be seen from Eq (3.3) that
Multiplying the left side of both sides of Eq (4.1) by h(εn)T, we obtain
From Lemma 4.3, we know that matrix P+Q is negative definite, so
which implies that
Finally, it follows from Lemma 4.4 that
that is,
At this point, we have proved that the numerical scheme (3.3) is stable. □
Theorem 4.2. For α∈(1,2), let 0<λh≤1, then the numerical scheme (3.3) is convergent for solving problem (1.1), that is, the following relationship is satisfied:
where en=(en1,en2,...,enM−1)T, eni=uni−Uni, and C1 is a constant independent of n, τ, and h.
Proof. Subtracting (3.3) from (3.2) and from scheme (3.4), we obtain
where Rn=(Rn1,Rn2,...,RnM−1)T.
Let's multiply both sides of (4.6) by h(en)T, which can be written as
Denoting En=h(en)TC(α)en, similar to the proof of Theorem 4.1,
Noticing that ‖e0‖2L2=0, the Eq (4.8) can be rearranged as
By the Lemma 4.5 (discrete Gronwall inequalities), we obtain
that is,
Thus, the theorem is proved. □
By the same idea, we can obtain the following theorem, whose proof we skip here.
Theorem 4.3. For α∈(1,2), let 0<λh≤1, then the numerical scheme (3.15) is stable for solving problem (1.2).
Theorem 4.4. For α∈(1,2), let 0<λh≤1, then the numerical scheme (3.15) is convergent for solving problem (1.2), that is, the following relationship is satisfied:
where en=(en1,en2,...,enM−1)T, eni=uni−Uni, and C2 is a constant independent of n, τ, and h.
5.
Numerical experiments
In this section, we present some numerical experiments to verify the convergence accuracy of the numerical schemes. Let
be the order of observation.
Example 5.1. Consider the following left Riemann-Liouville tempered fractional diffusion equation with initial-boundary value problem
where 1<α<2, and the linear source term is
and the exact solution is u(x,t)=et−λxx6.
Set the relationship between time step size and space step size as τ=h4 for all numerical experiments. Choose different α, λ, and space step size h, and use numerical scheme (3.3) to solve Example 5.1. The error and observation order results obtained are shown in Table 1. From Table 1, it can be seen that the space convergence order reaches the fourth order, and the numerical results are in perfect agreement with the theoretical analysis.
Example 5.2. Consider the following right Riemann-Liouville tempered fractional diffusion equation with initial-boundary value problem
where 1<α<2, and the linear source term is
and the exact solution is u(x,t)=et+λx(1−x)8.
The numerical scheme (3.15) is used to solve Example 5.2, and the obtained numerical results are presented in Table 2 for different α,λ, and space step size h. From Table 2, we find that the space convergence order is fourth order, and the numerical results are consistent with the theoretical results.
6.
Conclusions
In this paper, we focus on fourth-order numerical algorithms for one-sided tempered fractional diffusion equations. By using the quasi-compact technique, the fourth-order quasi-compact approximations of the tempered fractional derivatives and the effective second-order approximations of the third-order tempered derivatives are given. It is proved that the numerical schemes are stable and convergent. Finally, some experiments are given to demonstrate the effectiveness of the numerical schemes.
Acknowledgments
The author would like to express thanks to the editors and anonymous reviewers for their valuable comments and suggestions. This research is supported by the Xinjiang University of Political Science and Law.
Conflict of interest
The author has no conflict of interest.