Abstract
In this paper we introduce a new numerical method for solving time fractional partial differential equation. The time discretization is based on Diethelm’s method where the Hadamard finite-part integral is approximated by using the piecewise quadratic interpolation polynomials. The space discretization is based on the standard finite element method. The error estimates with the convergence order \(O(\tau ^{3-\alpha } + h^2), 0<\alpha <1\) are proved in detail by using the argument developed recently by Lv and Xu (SIAM J Sci Comput 38:A2699–A2724, 2016), where \(\tau \) and h denote the time and space step sizes, respectively. Numerical examples in both one- and two-dimensional cases are given.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
1 Introduction
In this paper, we will consider the numerical methods for solving the following time fractional partial differential equation
where \(\varOmega \subset \mathbb {R}^{d}, d=1, 2, 3\) is a convex polygonal/polyhedral domain and \( \varDelta \) denotes the Laplacian. Here f and q are two given functions and \(_0^C D_{t}^{\alpha } x(t), 0< \alpha <1\) denotes the Caputo fractional order derivative. Many application problems can be modelled by (1)–(3), for example, the thermaldiffusion in media with fractional geometry [28], highly heterogeneous aquifer [1], underground environmental problem [10], random walk [9, 25], etc.
To solve (1)–(3) numerically one needs to approximate the time fractional order derivative. There are three predominant approximations in literature: finite difference method (L1 scheme) [13, 20], the Grünwald–Letnikov method, [35, 36], Diethelm’s method [4, 7]. The L1 scheme is obtained by approximating the first order derivative with the finite difference quotients in the definition of the fractional derivative. The Grünwald–Letnikov method is based on the convolution quadrature and finally the Diethelm’s method is based on the approximation of the Hadamard finite-part integral.
Langlands and Henry [13] considered the L1 scheme for the Riemann–Liouville derivative and proved that the convergence order is \(O(\tau ^{2- \alpha })\) if \(u \in C^{2}[0, T]\). Lin and Xu [20] studied the L1 scheme for the Caputo fractional derivative and proved that the convergence order is also \(O(\tau ^{2- \alpha })\) if \(u \in C^{2}[0, T]\), see also [17, 21, 31]. Li and Ding [14] obtained a finite difference method with order \(O(\tau ^{2})\) if \( \, _{-\infty }^{R} D_{t}^{3- \alpha } u \in L^{1}(0, T)\), see also [12, 24, 26, 27, 29, 30], etc.
Yuste and Acedo [35] considered a Grünwald–Letnikov discretization of the Riemann–Liouville derivative and provided a von Neumann type stability analysis. Zeng et al. [36] introduced two fully discrete schemes with convergence order \(O(\tau ^{2- \alpha })\) if \(u \in C^{2}[0, T]\) by using fractional linear multistep method in time to approximate the convolution integral, see also [2, 11, 18, 22, 34], etc.
Diethelm [4] introduced a finite difference scheme to approximate the Riemann–Liouville fractional derivative by using the Hadamard finite-part integral and showed that the truncation error is \(O(\tau ^{2- \alpha })\) if \(u \in C^{2}[0, T]\). The scheme in [4] is obtained by approximating the Hadamard finite-part integral with the linear interpolation polynomials. This scheme is actually equivalent to the L1 scheme [13, 20] since the weights are the same. Ford et al. [6] applied the Diethelm’s method for solving time fractional partial differential equation and proved that the convergence order is \(O(\tau ^{2- \alpha })\) if \(u \in C^{2}[0, T]\). Higher order Diethelm’s schemes are also available in the literature, see [5, 7, 19, 33], etc.
Recently, Gao et al. [8] obtained a high order numerical differentiation formula with \(O(\tau ^{3- \alpha }), 0<\alpha <1\) for the Caputo fractional derivative by discretizing fractional derivative directly and applied this formula for solving a time fractional diffusion equation. But there are no error estimates in [8]. Li et al. [15] also introduced a high order \(O(\tau ^{3- \alpha }), 0<\alpha <1\) numerical method to approximate the Caputo fractional derivative and applied this method for solving time fractional advection–diffusion equations, see also [3, 16]. However the error estimates and stability analysis in [15] are provided only for \(\alpha \in (0, \alpha _{0})\) with some positive \(\alpha _{0} \in (0, 1)\), see also [3, 16]. More recently, Lv and Xu [23] proposed a higher order numerical method which is slightly different from the method in Gao et al. [8] for solving time fractional diffusion equation and proved that the scheme has the convergence order \(O(\tau ^{3- \alpha })\) for all \(\alpha \in (0, 1)\).
Yan et al. [33] introduced a numerical method for solving linear fractional differential equation with convergence order \(O(\tau ^{3- \alpha }), 0< \alpha <1\) by approximating the Hadamard finite-part integral with the quadratic interpolation polynomials following Diethelm’s idea in [4]. They obtained an asymptotic expansion of the error, but there are no error estimates proved in [33]. Recently, Li et al. [19] gave the detailed and thorough error estimates for the numerical method in [33] for solving the linear fractional differential equation.
In this paper, we will consider the numerical method for solving time fractional partial differential equation. The time discretization is based on the numerical method in [33] and the space discretization is based on the standard finite element method. The error estimates with convergence order \(O( \tau ^{3- \alpha } + h^2)\) are proved in detail by using the argument developed in [23] (see also [19]). The assumption of the regularity for the exact solution in our paper is \(u \in C^{3}([0, T], H^{2}(\varOmega )) = C^{3}(H^2)\) which is much weaker than the assumption for the solution u in Li et al. [15] where \(u \in C^{5}(H^{2})\). In [15], the exact solution u needs to satisfy \(\frac{\partial u(x, 0)}{\partial t} = \frac{\partial ^2 u(x, 0)}{\partial t^2} =0\) in order to obtain the required convergence order \(O(\tau ^{3- \alpha }), 0< \alpha <1\). Our numerical method has no such requirements for the solution u.
The paper is organized as follows. In Sect. 2, we consider the time discretization of the time fractional partial differential equations and prove that the numerical method has the convergence order \(O(\tau ^{3- \alpha })\) for all \(0< \alpha <1\). In Sect. 3, we consider the error estimates for solving time fractional partial differential equation in the fully discrete case where the spatial variables are discretized by using standard Galerkin finite element method. Finally in Sect. 4, we give some numerical examples in both one-dimensional and two-dimensional cases.
By C we denote a positive constant independent of the functions and parameters concerned, but not necessarily the same at different occurrences. By \( c_{0}, c_{1}, c_{2}\) we denote some particular positive constants independent of the functions and parameters concerned.
2 Time Discretization
In this section, we will consider the time discretization of (1)–(3). Recall that the Riemann–Liouville fractional derivative is defined by, with \(0< \alpha <1\),
Let \(\mathbb {N}\) denote the set of all natural numbers then, for \(p\notin \mathbb {N}\), on a general interval [a, b] Hadamard finite-part integral is defined in [4] as follows: with \(p >1\),
where
and \(\oint \) denotes the Hadamard finite-part integral. \(\lfloor {p}\rfloor \) denotes the largest integer not exceeding p, where \(p\not \in \mathbb {N}\).
It is easy to show that the Riemann–Liouville fractional derivative in (4) can be written as, [4]
We will approximate the Hadamard finite-part integral by using piecewise quadratic interpolation polynomial. Let \(n=2M\), where M denotes a fixed positive integer. Let \(0 = t_{0}< t_{1}< t_{2}< \dots< t_{2j}< t_{2j+1}<\dots < t_{2M} =T\) be a partition of [0, T] and \(\tau \) the step size. For simplicity of notation, we assume that \(T=1\) below. At the point \(t_{2j} = \frac{2j}{2M} T = \frac{2j}{2M}\), the Eq. (1) can be written as
and at the point \(t_{2j+1} = \frac{2j+1}{2M}\), the equation can be written as
Let us first consider the discretization of (8). Note that
where the integral denotes the Hadamard finite-part integral.
We will approximate the integral by a piecewise quadratic interpolation polynomial with the equispaced nodes \( 0, \frac{1}{2j}, \frac{2}{2j}, \dots , \frac{2j}{2j}, \, j=1, 2, \dots , M\). More precisely, for the sufficiently smooth function g(w), we have
where \(g_{2}(w)\) is the piecewise quadratic interpolation polynomial of g(w) defined on the nodes \( 0< \frac{1}{2j}< \frac{2}{2j}< \dots < \frac{2j}{2j} =1, \, j=1, 2, \dots , M\), and \(E_{2j} (g)\) is the remainder term.
Lemma 1
[19, Lemma 3.1] Let \(0< \alpha <1\). Assume that \( g \in C^{3}[0, 1]\). Then, with \( j=1, 2, \dots , M\),
where, with \(j=1\),
and, with \(j=2, 3, \dots , M\),
Here
and
Next we consider the discretization of (9). At the point \(t_{2j+1}= \frac{2j+1}{2M}, \, j=1,2, \dots , M-1\) we have
We approximate this Hadamard finite-part integral by a piecewise quadratic interpolation polynomial with the equispaced nodes \( 0, \frac{1}{2j+1}, \frac{2}{2j+1}, \dots , \frac{2j}{2j+1}, \, j=1, 2, \dots , M-1\). More precisely, we have, for the sufficiently smooth function g(w),
where \(g_{2}(w)\) is the piecewise quadratic interpolation polynomial of g(w) defined on the nodes \(0, \frac{1}{2j+1}, \frac{2}{2j+1}, \dots , \frac{2j}{2j+1}, \, j=1, 2, \dots , M-1\) and \(E_{2j+1} (g)\) is the remainder term.
Lemma 2
[19, Lemma 3.2] Let \(0< \alpha <1\). Assume that \(g \in C^{3}[0, 1]\). Then
where \( \alpha _{k, 2j+1} = \alpha _{k, 2j}, \, k=0, 1, 2, \dots , 2j, \, j=1, 2, \dots , M-1\) and \(\alpha _{k, 2j}\) are given in Lemma 1.
By using (10)–(12), we obtain the following approximation of the Riemann–Liouville fractional derivative \(_0^R D_{t}^{\alpha } u(x,t)\) at \(t = t_{2j}, \, j=1, 2, \dots , M\)
where \(R_{2}^{2j} \le C \tau ^{3- \alpha } \big ( \max _{0 \le s \le 1} | u^{\prime \prime \prime } (x,s) | \big ) \) and we denote \( R_{2}^{2j} = O(\tau ^{3- \alpha })\) [33] and the weights \(w_{k, 2j}, k=0, 1, 2, \dots , 2j, j=1, 2, \dots , M\) satisfy
Similarly, we have at \(t = t_{2j+1}, \, j=1, 2, \dots , M - 1\),
where
For the Caputo fractional derivative \(_0^C D_{t}^{\alpha } u(x,t)\) at \(t = t_{2j}, \, j=1, 2, \dots , M\), we have, noting that \(_0^R D_{t}^{\alpha } u(x,0) = u(x,0) \, _0^R D_{t}^{\alpha } (1) = \frac{u(x,0)}{\varGamma (1- \alpha )} t^{-\alpha }\),
where the weights, with \(k=0, 1, 2, \dots , 2j-1, j=1, 2, \dots , M\),
Similarly, we have at \(t = t_{2j+1}, \, j=1, 2, \dots , M-1\),
where, with \(k=0, 1, 2, \dots , 2j, \, j=1, 2, \dots , M-1\),
The exact solution u in (8)–(9) then satisfies, with \( l=2, 3, \dots , 2M\),
or
where \(d_{k, l} = -\bar{w}_{k, l}/\bar{w}_{0, l}, k=1, 2, \dots , l, \; l=2, 3, \dots , 2M\), where \(I_{l}\) is defined by
For the simplification of notation, we omit the dependence of the exact solution u(x, t) on x below. Let \( u^{l} \approx u(x,t_{l}), \, l=0, 1, 2, \dots , 2 M\) denote the approximate solution of \( u(x,t_{l})\). We define the following numerical method to approximate the exact solutions in (23), with \(l=2,3, \dots , 2M\),
where \(\tilde{I}_{l}\) is some approximation of \(I_{l}\) discussed below in (26). Here we assume that \(u^{0} = u_{0}\) and \(u^{1}\) will be approximated below in (25).
To approximate \(u (x,t_{1}) \) with the required accuracy \(O(\tau ^{3- \alpha })\) which will be the convergence order of our numerical method (24), we divide the interval \([0, t_{1}]\) by the equispaced nodes \(0= t_{1}^{(0)}< t_{1}^{(1)}< \dots < t_{1}^{(n_{1})}= t_{1}\) with step size \( \widetilde{\tau } \) such that \(\widetilde{\tau }^{2- \alpha } \approx \tau ^{3 - \alpha }\), where \(n_{1}\) is some positive integer. We then apply the numerical method with the convergence order \( O ( \tilde{\tau }^{2- \alpha })\) in [4] to get the approximate value \(u^{1} \approx u(x,t_{1})\) such that, with \(e^{1} = u^{1} - u(x, t_{1})\),
where \( \Vert \cdot \Vert \) denotes the \(L^{2}\) norm.
We also need to approximate the integral \(I_{l}\) in (23) with the required accuracy \(O(\tau ^{3})\) which we shall use in (27). Let \(n_{2}\) be some positive integer, we divide the interval \([0, t_{1}]\) by the equispaced nodes \(0= t_{1}^{(0)}< t_{1}^{(1)}< \dots < t_{1}^{(n_{2})}= t_{1}\) with step size \(\overline{\tau }\) such that \( \overline{\tau }^{2} \approx \tau ^{3-\alpha }\). We then apply the composite trapezoidal quaduature rule on \([0, t_{1}]\) which has the convergence order \( O ( \overline{\tau }^2)\). More precisely, we have
where \(\tilde{u} (x,s)\) is the piecewise linear interpolation polynomial of u(x, s) on \([0, t_{1}]\), which implies that \( I_{l} - \tilde{I}_{l} = O(\tau ^{3})\). We need this approximation below in (27).
Let \( e^{l}= u^{l}- u(x,t_{l}), l=0, 1, \dots , 2M\). Subtracting (23) from (24), we have, by (26),
where \(e^{0} =0\) and \(e^{1}\) is approximated in (25) and \(d_{k, l} = -\bar{w}_{k, l}/\bar{w}_{0, l}, k=1, 2, \dots , l, \, l =2, 3, \dots , 2M\) are defined as in (21) and (22).
Note that, by (21) and (22), \(d_{1, l}= -\bar{w}_{1, l}/\bar{w}_{0, l}, l=2, 3, \dots , 2M\) is a constant which is independent on \(l=2, 3, \dots , 2M\). Further we define
We now denote
where \(\eta \) is a constant independent on \(l=1, 2, \dots , 2M\). We have, with \(l=2, 3, \dots , 2M\),
Denote
we have, with \(\bar{d}_{1, l } = \eta \),
Lemma 3
[19, Lemma 3.3] For \(0< \alpha <1\), the coefficients in (30) satisfy, with \(l =2, 3, \dots , 2M\),
Next we will consider the error estimates in Theorem 1. To prove the error estimates in Theorem 1 below, we need the further assumption for the initial approximation \(u^{1}\). Let \(\bar{e}^{1} = e^{1} -\eta e^{0}=e^{1}\) be defined as in (28), where \(e^{1} = u^{1} - u(x, t_{1})\). By (32)–(34), we see that \(1 \le (\bar{d}_{l, l})^{-1} \le c_{0} \tau ^{-\alpha }\). Thus, with some fixed \(l=2, 3, \dots , 2M\), choosing \(\tilde{\tau }\) sufficiently small as in (25) , we may find
and a constant \(c_{1}>0\) such that,
which we need in (40) below. By (34), we have, noting that \(\bar{d}_{l, l} <1\) by (32) and (33),
Now it is ready to introduce the following error estimates.
Theorem 1
Let \(u(x,t_{l})\) and \(u^{l}, l=0, 1, \dots , 2M\) be the exact and approximate solutions of (23) and (24), respectively. Assume that \(u(x,t) \in C^{3}[0, T]\). Further assume that \(u^{0}= u_{0}\) and \(u^{1}\) satisfies (36). Then there exists a constant \(C= C(\alpha , f, T)\) such that
Proof
Multiplying \(2 \bar{e}^{l}\) in both sides of (30), we have, with \(l =2, 3, \dots , 2M\),
Note that
We have, with \( l =2, 3, \dots , 2M\),
We write, with \( l =2, 3, \dots , 2M\),
By Cauchy-Schwarz inequality, we have, with \(l=2, 3, \dots , 2M\), noting that \(\bar{d}_{l, l} >0\) from (32),
By (33) and noting that \(\bar{d}_{1, l} = \eta \) , we have
By (31), we have, noting that \(\bar{d}_{1, l} = \eta \),
Note that, by (21) and (22), \(\bar{w}_{0, l} = \frac{2^{-\alpha } (\alpha +2)}{\varGamma (3- \alpha )}, \, l=2, 3, \dots , 2M\) which is independent on \(l=2, 3, \dots , 2M\). Further we define
We now denote the norm, with \(l=1, 2, \dots , 2M\),
We then have, with \( l =2, 3, \dots , 2M\),
By the initial approximation estimate (36), we have, with \(l=2, 3, \dots , 2M\),
We next prove the following by the mathematical induction, with \(k=2,3, \dots , l, \; l =2, 3, \dots , 2M \),
Assume that (40) and (41) hold true for \(k=1, 2, \dots , l-1, \; l =2, 3, \dots , 2M\), we have, for \(k=l\), by (32),
which is (41).
By (33), we have, with \( l =2, 3, \dots , 2M\),
where the last inequality is due to the fact \(R_{2}^{s} = O(\tau ^{3- \alpha }), 1 \le s \le l\) by (35), (18), (20). Hence we have \( \Vert \bar{e}^{l}\Vert _{1} \le C \tau ^{3-\alpha }, \; l=2, 3, \dots , 2M. \)
Further we have, by (31), with \(l=2,3, \dots , 2M\),
Together these estimates complete the proof of Theorem 1.
3 The Fully Discrete Scheme
In this section, we shall consider the fully discretization scheme for solving (1)–(3). Here we only consider the error estimates for the homogeneous Dirichlet boundary condition, i.e., \(q=0\) in (3). But the error estimates in Theorem 2 is also true for the non homogeneous Dirichlet boundary condition. Let \(S_{h}\subseteq H_{0}^{1} (\varOmega )\) denote the standard linear finite element space and h the space step size. Let \( R_{h}:H_{0}^{1}(\varOmega )\rightarrow S_{h}\) denote the Ritz projection defined by, for \(\forall \varphi \in H_{0}^{1}(\varOmega )\),
It is well known that, see [32],
We now consider the fully discrete scheme for solving (23). Let \(U_{h}^{l} \in S_{h}, l=0, 1, 2, \dots , {2M}\) denote the approximation of \(u(x, t_{l})\). We choose \(U_{h}^{0} = R_{h} u_{0}\) and we assume that \(U_{h}^{1} \in S_{h}\) is a suitable approximation of \(u(t_{1})\) which can be obtained by using some special numerical methods and satisfies the condition (45) below. For \(l=2, 3, \dots , 2M\), we define the following finite element method: find \(U_{h}^{l} \in S_{h}\), such that
Theorem 2
Let \(u(x,t_{l})\) and \(U_{h}^{l}, l=0, 1, \dots , 2M\) be the exact and approximate solutions of (23) and (44), respectively. Assume that \(u \in C^{3} ([0, T];H^{2}(\varOmega ))\) and the initial condition \(U_{h}^{0}=R_{h}u_{0}\). Assume that there exists a constant \(c_{2}\) such that
where \(\delta _{2}^{1}=( R_{h} - I) \big ( \, _0^C D_{t}^{\alpha } u (t_{1})-R_{2}^{1} \big )+R_{2}^{1}\) and \(R_{2}^{1}\) is defined in (36). Then there exists a constant \(C= C(\alpha , f,T)\) such that, with \(l=2, 3, \dots , 2M\),
Further we have
Proof
Denote, with \(l =1, 2, \dots , N\),
By (43), we have
We now estimate \(\theta ^{l}, l=2, 3, \dots , N\). We have, for \(\forall \, \chi \in S_{h}\),
Hence we get
where \(\delta _{2}^{l}=( R_{h} - I) \big ( \, _0^C D_{t}^{\alpha } u (t_{l})-R_{2}^{l} \big )+R_{2}^{l}.\)
From the triangle inequality, we have
Note that
We have
Let
Applying (45), we have, noting that \(R_{2}^{1} = O(\tau ^{3- \alpha })\) by (35),
We next consider the case for \(l=2, 3, \dots , N\). We may write the equation (47) as
Following the same argument as in the proof of Theorem 1, we may get, with \(l=2, 3, \dots , 2M\),
Together these estimates with (46) completes the proof of Theorem 2. \(\square \)
4 Numerical Simulations
In this section, we will consider four examples. Examples 1 and 2 are one-dimensional problems and Examples 3 and 4 are two-dimensional problems.
Example 1
Consider
where \(f(x, t) =\frac{\varGamma (m+1)}{\varGamma (m+1 -\alpha )} t^{m- \alpha } \sin (2 \pi x) + 4 \pi ^2 t^{m} \sin (2 \pi x) \). The exact solution is \(u(x, t)= t^{m} \sin (2 \pi x)\). In Table 1, we choose \(m=3.5\). In this case the exact solution \(u(x, \cdot ) \in C^{3}[0, T]\).
The main purpose is to check the order of convergence of the numerical method with respect to the time step size \( \tau \) for the different fractional orders \(\alpha \). For various choices of \(\alpha \in (0,1)\), we compute the errors in the \(L^{2}\) norm at \(T=1\). We use the linear finite element space \(S_{h}\) with the space step size \(h = 1/2^6\) which is sufficiently small such that the error will be dominated by the time discretization of the method. We choose the time step size \( \tau = 1/2^{l}, l=3, 4, 5, 6,7\), i.e, and we divide the interval [0, T] into \(N=1/\tau \) subintervals with nodes \(0=t_{0}< t_{1}< \dots < t_{N}=1\). Then we compute the error \(e(t_{N})=\Vert u(x, t_{N}) - U_{h}^{N} \Vert \). By Theorem 2, we have
To observe the order of convergence we shall compute the error \(\Vert e(t_{N})\Vert \) at \(t_{N}=1\) with respect to the different values of \(\tau \). Denote \(\Vert e_{\tau _{l}}(t_{N})\Vert \) the error at \(t_{N}=1\) with respect to the time step size \(\tau _{l}\). Let \(\tau _{l}= \tau = 1/ 2^{l} \) for a fixed \(l =3, 4, 5, 6, 7\). We then have
which implies that the order of convergence satisfies \( 3- \alpha \approx \text{ log }_{2} \Big ( \frac{\Vert e_{\tau _{l}}(t_{N})\Vert }{\Vert e_{\tau _{l+1}}(t_{N})\Vert } \Big )\). In Table 1, we compute the orders of convergence for the different values of \(\alpha \). The numerical results are consistent with the theoretical results in Theorem 2. We also compare the numerical results of our scheme with the results obtained by the scheme in [8]. We observe that the convergence rate of our method is slightly higher than the convergence rate of the method in [8] for small \(\alpha \in (0, 1)\). We will investigate this interesting issue in our future work. We have the similar observations in other examples below.
Example 2
Consider
where \( u(x, t) = e^{x} t^{4+\alpha }\) and \(f(x, t) = e^{x} \frac{\varGamma (5+ \alpha )}{24} t^{4} - t^{4+ \alpha } e^{x}\).
We use the same notations as in the experiments in Example 1. In Table 2, we observe that the numerical results are consistent with the theoretical results.
Next we will consider two examples in two-dimensional cases. Let us first introduce the algorithm for solving the following time fractional partial differential equations in two-dimensional case by using finite element method.
where \(\varOmega = (0, 1) \times (0, 1)\) and \( \frac{\partial u(x, t)}{\partial n}\) denotes the normal derivative on the boundary \(\partial \varOmega \) and q(x, t) is some function defined on the boundary \(x \in \partial \varOmega \) and \(t \in (0, T)\). Here \(\kappa \) is a constant. When \(\kappa \) is sufficiently large, (58) is reduced to the Dirichlet boundary condition. In our numerical examples below, we will only consider the Dirichlet boundary conditions.
The variational form of (56)–(58) is to find \(u(t) \in L^{2}(\varOmega ) \), such that
Let \(x= (x_{1}, x_{2})\). Let \(0= x_{1}^{0}< x_{1}^{1}< \dots < x_{1}^{M1}=1\) be a partition of [0, 1] on the \(x_{1}\) axis and h the step size. Similarly we let \(0= x_{2}^{0}< x_{2}^{1}< \dots < x_{2}^{M1}=1\) be a partition of [0, 1] on the \(x_{2}\) axis and h the step size. For simplicity of notation, we use the same step size on both \(x_{1}\) and \(x_{2}\) axes. We divide the domain \(\varOmega \) into the small triangles which have the same sizes.
Let \(S_{h}\) denote the linear finite element space defined on \(\varOmega \). Let \(P_{0}, P_{1}, \dots P_{N_{h}}\) denote all the nodes on the triangulation of \(\varOmega \). Let \(\varphi _{j}, j=0, 1, 2, \dots , N_{h}\) be the linear basis function corresponding to the node \(P_{j}, j=0, 1, 2, \dots , N_{h}\). The finite element method is to find \(u_{h} \in S_{h}\) such that
Let \(0< t_{0}< t_{1}< \dots < t_{2M}=T\) be the time partition of [0, T] and \(\tau \) be the step size. Below we only consider the time discretization at the even node \(t_{n}, n=2, 4, \dots , 2M\) to get the idea of the algorithm of our method. (Similarly one can consider the discretization at the odd node \(t_{n}, n=1, 3, \dots , 2M-1\)). By (18), we have
with some suitable weights \(w_{k, n}, k=0, 1, \dots , n\).
Define the following time discretization scheme: find \(U^{n} \approx u_{h}(t_{n}), n=2, 4, \dots , t_{2M}=T\) such that, noting that \(\, _0^C D_{t}^{\alpha } u(t_{n}) =\, _0^R D_{t}^{\alpha } ( u(t_{n}) - u_{0})\),
or
Let \(U^{n} = \sum _{j=0}^{N_{h}} \alpha _{j}^{n} \varphi _{j}\) be the approximate solution of \(u_{h}(t_{n}), n=1, 2, \dots , N\). Choose \(\chi = \varphi _{l}, l=0, 1, \dots , N_{h}\), we have
Denote
We have the matrix form
Solving this system, we get the finite element solution \(\alpha ^{n}, n=2, 4, \dots , 2M\).
Example 3
Consider
where \(\varOmega = (0, 1) \times (0, 1)\). The exact solution is \( u(x, t) = t^{m} \sin (2 \pi x_{1}) \sin (2 \pi x_{2})\) for some \(m >0\) and
Here \( q(x, t) = t^{m} \sin (2 \pi x_{1}) \sin (2 \pi x_{2})\) and \(u_{0} (x) =0\). In Table 3, we choose \(m=3.5\).
We will consider the Dirichlet boundary condition and therefore we choose \(\kappa =10000\) in our numerical simulation. For various choices of \(\alpha \in (0,1)\), we compute the errors at \(T=1\). We use the linear finite element space \(S_{h}\) with the space step size \(h = 1/2^6\) which is sufficiently small such that the error will be dominated by the time discretization of the method. We choose the time step size \( \tau = 1/2^{l}, l=3, 4, 5, 6,7\), i.e, we divide the interval [0, T] into \(2M=1/\tau \) subintervals with nodes \(0=t_{0}< t_{1}< \dots < t_{2M}=1\). In Table 3, we compute the orders of convergence for the different values of \(\alpha \). The numerical results are consistent with the theoretical results in Theorem 2.
Example 4
Consider
where \(\varOmega = (0, 1) \times (0, 1)\). The exact solution is \( u(x, t) = e^{x_{1} + x_{2}} t^{4+\alpha }\) and \(f(x, t) = e^{x_{1} +x_{2}} \big ( \frac{\varGamma (5+ \alpha )}{24} t^{4} \big ) - t^{4+ \alpha } \big (2 e^{x_{1} + x_{2}}\big )\) Here \( q(x, t) =e^{x_{1} +x_{2}} t^{4+\alpha }\) and \(u_{0} (x) =0\).
We use the same notations as in the experiments in Example 3. In Table 4, we observe that the numerical results are consistent with the theoretical results.
References
Adams, E.E., Gelhar, L.W.: Field study of dispersion in a heterogeneous aquifer: 2. Spatial moments analysis. Water Resour. Res. 28, 3293–3307 (1992)
Cuesta, E., Lubich, Ch., Palencia, C.: Convolution quadrature time discretization of fractional diffusion-wave equations. Math. Comput. 75, 673–696 (2006)
Cao, J., Li, C., Chen, Y.: High-order approximation to Caputo derivatives and Caputo-type advection–diffusion equation II. Fract. Calc. Appl. Anal. 18, 735–761 (2015)
Diethelm, K.: An algorithm for the numerical solution of differential equations of fractional order. Electron. Trans. Numer. Anal. 5, 1–6 (1997)
Ford, N.J., Pal, K., Yan, Y.: An algorithm for the numerical solution of space-fractional partial differential equations. Comput. Methods Appl. Math. 15, 497–514 (2015)
Ford, N.J., Xiao, J., Yan, Y.: Stability of a numerical method for a space–time-fractional telegraph equation. Comput. Methods Appl. Math. 12, 1–16 (2012)
Ford, N.J., Xiao, J., Yan, Y.: A finite element method for time fractional partial differential equations. Fract. Calc. Appl. Anal. 14, 454–474 (2011)
Gao, G.-H., Sun, Z.-Z., Zhang, H.-W.: A new fractional numerical differentiation formula to approximate the Caputo fractional derivative and its applications. J. Comput. Phys. 259, 33–50 (2014)
Gorenflo, R., Mainardi, F.: Random walk models for space fractional diffusion processes. Fract. Calc. Appl. Anal. 1, 167–191 (1998)
Hatano, Y., Hatano, N.: Dispersive transport of ions in column experiments: an explanation of long-tailed profiles. Water Resour. Res. 34, 1027–1033 (1998)
Jin, B., Lazarov, R., Zhou, Z.: On two schemes for fractional diffusion and diffusion-wave equations. Preprint arXiv:1404.3800 (2014)
Jin, B., Lazarov, R., Zhou, Z.: An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data. IMA J. Numer. Anal. 36, 197–221 (2016)
Langlands, T.A.M., Henry, B.I.: The accuracy and stability of an implicit solution method for the fractional diffusion equation. J. Comput. Phys. 205, 719–736 (2005)
Li, C., Ding, H.: Higher order finite difference method for the reaction and anomalous-diffusion equation. Appl. Math. Model. 38, 3802–3821 (2014)
Li, C., Wu, R., Ding, H.: High-order approximation to Caputo derivatives and Caputo-type advection–diffusion equations. Commun. Appl. Ind. Math. 6, e536 (2014)
Li, H., Cao, J., Li, C.: High-order approximation to Caputo derivatives and Caputo-type advection–diffusion equation III. J. Comput. Appl. Math. 299, 159–175 (2016)
Li, X., Xu, C.: Existence and uniqueness of the weak solution of the space–time fractional diffusion equation and a spectral method approximation. Commun. Comput. Phys. 8, 1016–1051 (2010)
Li, W., Xu, D.: Finite central difference/finite element approximations for parabolic integro-differential equations. Computing 90, 89–111 (2010)
Li, Z., Yan, Y., Ford, N.J.: Error estimates of a high order numerical method for solving linear fractional differential equation. Appl. Numer. Math. (2016). doi:10.1016/j.apnum.2016.04.010
Lin, Y., Xu, C.: Finite difference/spectral approximations for the time-fractional diffusion equation. J. Comput. Phys. 225, 1533–1552 (2007)
Lin, Y., Li, X., Xu, C.: Finite difference/spectral approximations for the fractional cable equation. Math. Comput. 80, 1369–1396 (2011)
Lubich, Ch., Sloan, I.H., Thomée, V.: Nonsmooth data error estimates for approximations of an evolution equation with a positive-type memory term. Math. Comput. 65, 1–17 (1996)
Lv, C., Xu, C.: Error analysis of a high order method for time-fractional diffusion equations. SIAM J. Sci. Comput. 38, A2699–A2724 (2016)
McLean, W., Mustapha, K.: Time-stepping error bounds for fractional diffusion problems with non-smooth initial data. J. Comput. Phys. 293, 201–217 (2015)
Metzler, R., Klafter, J.: The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics. J. Phys. A Math. Gen. 37, 161–208 (2004)
Mustapha, K., McLean, W.: Uniform convergence for a discontinuous Galerkin, time-stepping method applied to a fractional diffusion equation. IMA J. Numer. Anal. 32, 906–925 (2012)
Mustapha, K., Schötzau, D.: Well-posedness of hp-version discontinuous Galerkin methods for fractional diffusion wave equations. IMA J. Numer. Anal. 34, 1426–1446 (2014)
Nigmatulin, R.: The realization of the generalized transfer equation in a medium with fractal geometry. Phys. Stat. Sol. B 133, 425–430 (1986)
Sousa, E.: Numerical approximations for fractional diffusion equation via splines. Comput. Math. Appl. 62, 938–944 (2011)
Sousa, E.: A second order explicit finite difference method for the fractional advection diffusion equation. Comput. Math. Appl. 64, 3143–3152 (2012)
Sun, Z.-Z., Wu, X.: A fully discrete scheme for a diffusion wave system. Appl. Numer. Math. 56, 193–209 (2011)
Thomée, V.: Galerkin Finite Element Methods for Parabolic Problems. Springer, Berlin (2007)
Yan, Y., Pal, K., Ford, N.J.: Higher order numerical methods for solving fractional differential equations. BIT Numer. Math. 54, 555–584 (2014)
Yuste, S.B.: Weighted average finite difference methods for fractional diffusion equations. J. Comput. Phys. 216, 264–274 (2006)
Yuste, S.B., Acedo, L.: An explicit finite difference method and a new von Neumann-type stability analysis for fractional diffusion equations. SIAM J. Numer. Anal. 42, 1862–1874 (2005)
Zeng, F., Li, C., Liu, F., Turner, I.: The use of finite difference/element approaches for solving the time-fractional subdiffusion equation. SIAM J. Sci. Comput. 35, A2976–A3000 (2013)
Author information
Authors and Affiliations
Corresponding author
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Li, Z., Liang, Z. & Yan, Y. High-Order Numerical Methods for Solving Time Fractional Partial Differential Equations. J Sci Comput 71, 785–803 (2017). https://doi.org/10.1007/s10915-016-0319-1
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10915-016-0319-1