Error estimates of high-order numerical methods for solving time fractional partial differential equations

Abstract Error estimates of some high-order numerical methods for solving time fractional partial differential equations are studied in this paper. We first provide the detailed error estimate of a high-order numerical method proposed recently by Li et al. [21] for solving time fractional partial differential equation. We prove that this method has the convergence order O(τ3−α) for all α ∈ (0, 1) when the first and second derivatives of the solution are vanish at t = 0, where τ is the time step size and α is the fractional order in the Caputo sense. We then introduce a new time discretization method for solving time fractional partial differential equations, which has no requirements for the initial values as imposed in Li et al. [21]. We show that this new method also has the convergence order O(τ3−α) for all α ∈ (0, 1). The proofs of the error estimates are based on the energy method developed recently by Lv and Xu [26]. We also consider the space discretization by using the finite element method. Error estimates with convergence order O(τ3−α + h2) are proved in the fully discrete case, where h is the space step size. Numerical examples in both one- and two-dimensional cases are given to show that the numerical results are consistent with the theoretical results.


Introduction
Consider the following time fractional partial differential equation u(x, 0) = u 0 , x ∈ Ω, (1.2) where 0 < α < 1 and ∆ denotes the Laplacian and Ω ⊂ R d , d = 1, 2, 3 is a bounded and regular domain.Here f is a given function and C 0 D α t v(t) denotes the Caputo fractional order derivative defined by where v (s) = dv(s) ds denotes the derivative of the first order.Many engineering and physical problems can be modeled by using (1.1)-(1.3),for example, thermal diffusion in media with fractional geometry [29], highly heterogeneous aquifer [1], underground environmental problems [13], random walks [12,28], etc.
Since we use the different approaches to obtain (1.6)-(1.7)and (1.10), respectively, the weights in (1.6)-(1.7)and (1.10) are completely different.Therefore it is natural to consider the error estimates of the numerical methods developed in Li et al. [21] for solving (1.1)-(1.3)by using (1.10) which is one of the objects in the current paper.
Note that the discretization scheme (1.10) must satisfy the initial conditions v (t 0 ) = v (t 0 ) = 0 in order to get the convergence order O(τ 3−α ).In this paper we will introduce a new scheme which does not require v satisfy v (t 0 ) = v (t 0 ) = 0. We can also show that this scheme has the convergence order O(τ 3−α ) for all α ∈ (0, 1).Let us see how to construct this new scheme.For k = 2, 3, . . ., N , we will use the following way to approximate the Caputo derivative at t k , where P j 2 (s) denotes the quadratic interpolation polynomials defined on the nodes s = t j−1 , t j , t j+1 as above.To apply this scheme to solve the time fractional partial differential equation (1.1)-(1.3),we have to obtain the approximate value of the solution at t 1 with the required accuracy by using other numerical methods, then we use this scheme to approximate the solutions at t k with k = 2, 3, . . ., N .In (1.11), we use the same notations wi,k as in (1.10), but they denote the different values.It is easy to show that the truncation error of (1.11) satisfies We remark that the new scheme (1.12) is different from (1.6)-(1.7)and (1.10).In section 3, we apply the scheme (1.12) to solve the time fractional partial differential equation (1.1)-(1.3)and prove that this scheme has the convergence order O(τ 3−α ) for all α ∈ (0, 1).
The main contributions of this paper are as follows: (i).We obtain the properties of the weights in (1.10), i.e., Lemmas 2.1 and 2.2 for the approximation of the Caputo derivatives, which are important to prove the stability and error estimates of the numerical methods for solving (1.1)- (1.3).
(iii).We introduce a new time discretization scheme for solving (1.1)-(1.3)by using (1.12) and prove that this scheme has the convergence rate O(τ 3−α ) for all α ∈ (0, 1).This scheme does not require that the first and second time derivatives of the solutions are vanish at t = 0.
The paper is organized as follows.In Section 1, we consider the time discretization for solving time fractional partial differential equations and prove that the numerical method proposed in [21] has the convergence order O(τ 3−α ) for all 0 < α < 1.In Section 2, we consider error estimates of the scheme proposed in [21] for solving time fractional partial differential equations in the fully discrete case where the spatial variables are discretized by using the standard Galerkin finite element method.In Section 3, we consider the error estimate of the new time discretization method constructed by using (1.12) for solving time fractional partial differential equation in the fully discrete case.Finally in Section 4, we give numerical examples to show that the numerical results are consistent with our theoretical results.
By C we denote a positive constant independent of the functions and parameters concerned, but not necessarily the same at different occurrences.
Finally we estimate (2.34).It is easy to consider the case for k = 2.Here we only consider the case for k = 3, 4, . . ., N , we then have We shall show that Assume (2.35) holds at the moment, we then have Similarly we may show (2.36) holds also for k = 2, which we need in the proof of (2.39).Thus, by (2.21), (2.23), (2.24), with k = 3, 4, . . ., N , ) which implies that ( dk,k ) −1 < d −1 k,k < c 0 τ −α for some constant c 0 .It remains to show (2.35).In fact, we have In the estimates below, we will frequently use the binomial expansion, with β ∈ R, (2.38) By using (2.38), we have After some tedious but simple calculation, we obtain 1) .
After some simple calculation, we may get Z. Li, Y. Yan Hence (2.35) follows from the fact, with m ≥ 1, Together these estimates complete the proof of Lemma 2.2.
Proof of Theorem 2.1.We first consider the case for k = 1.Multiplying 2u 1 in both sides of (2.20) and integrating on Ω, we have, By using the Cauchy-Schwarz inequality and Young inequality, we have which implies that where in the last inequality, we use the fact, by (2.18) and (2.36), ( d2,2 ) (2.39) Thus we get

.40)
Here we choose the bound 4 √ 6 in order to make sure that this bound is the same bound as obtained in (2.42) for k = 2, 3, . . ., N below.
We now turn to the case for k = 2, 3, . . .N .Multiplying 2ū k in both sides of (2.28) and integrating on Ω, we have Using the inequality (2.33), we have Denote the norm, with k = 1, 2, . . ., N , We then have, with k = 2, 3, . . ., N , We will use mathematical induction to prove, with k = 2, 3, . . ., N , Let us first consider the case for k = 2.We have

The fully Discretization scheme
In this section, we shall consider the fully discretization scheme for solving (1.1)-(1.3).Here we only consider the error estimates for the homogeneous Dirichlet boundary condition.But the error estimates in Theorem 3.1 is also true for the nonhomogeneous Dirichlet boundary condition.

A new time discretization scheme
In this section, we shall consider the new approximation scheme (1.12).The weights wi,k , i = 0, 1, . . ., k, k = 2, 3, . . ., N in (1.12) satisfy the following: For k ≥ 5, we have Let us now turn to the finite difference method for solving (1.1)-(1.3)by using the scheme (1.12).For simplicity, we only consider the case for k ≥ 4 since the weights wi,k has the same expressions for all k ≥ 4. Similarly we can consider the case for k = 2, 3.
At t k , k = 4, 5, . . ., N , the equation (1.1)-( 1.3) has the form The equation (4.49) can be written as (4.51) Let u k ≈ u(x, t k ), k = 1, 2, . . ., N denote the approximate solution of u(x, t k ).Assume that we obtain the starting approximation u 1 , u 2 , u 3 with the required accuracy by using other methods.To calculate u k , k = 4, 5, . . ., N , we define the following time discretization scheme for solving (1.1)-(1.3):Let us now turn to the stability proof of u k , k = 1, 2, . . ., N .As in Section 2, the stability of the starting approximations u 1 , u 2 , u 3 can be considered separately which can be done easily.We may use the same ways to prove the corresponding lemmas and theorems as in Section 3. Therefore we only state the corresponding lemmas and theorems below and leave the proofs for the readers.Now we come to the main theorem in this section.
Theorem 4.1.Let 0 < α < 1.Let u k , k = 0, 1, . . ., N be the approximate solutions defined by (4.52).Assume that u(•, x) ∈ C 3 [0, T ].Assume that the starting approximations u 1 , u 2 , u 3 satisfy the required stability estimates.Then there exists a constant C = C(α, T ) such that For the error estimates of (4.52), we have the following Theorem 4.2.Let 0 < α < 1.Let u(x, t k ) and u k , k = 0, 1, . . ., N be the exact and the approximate solutions of (4.51) and (4.52), respectively.Assume that u(x, •) ∈ C 3 [0, T ].Further we assume that the starting errors e l = u l − u(x, t l ), l = 1, 2, 3 satisfy the required accuracy.Then there exists a constant C = C(α, f, T ) such that Remark 4.2.We can also consider the error estimates of the fully discrete scheme for solving (4.52) and obtain the similar results as in Theorem 3.1.We omit it here.5. Numerical simulations 5.1.Time discretization scheme (2.19) proposed in Li et al. [21] In this subsection, we will consider two numerical examples to investigate the time convergence rates of the numerical method (2.19) proposed in Li et al. [21].The first example is for the one-dimensional case and the second example is for the two-dimensional case.
For the various choices of α ∈ (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 τ = 1/2 l , l = 3, 4, 5, 6, 7, i.e, we divide the interval [0, T ] into N = 1/τ subintervals with nodes 0 = t 0 < t 1 < • • • < t N = 1.In Table 2, we compute the orders of convergence for the different values of α.The numerical results are consistent with the theoretical results in Theorem 2.2. .

= 2 3
.67) To observe the order of convergence we shall compute the error e(t N ) at t N = 1 with respect to the different values of τ .Denote e τ l (t N ) the error at t N = 1 for the time step size τ l .Let τ l = 1/2 l for a fixed l = 3, 4, 5, 6, 7. We then have e τ l (t N ) e τ l+1 (t N ) −α , which implies that the order of convergence satisfies 3−α ≈ log 2 eτ l (t N )

Table 1 .
The time convergence rate with h = 1/2 6 at T = 1 in Example 5.1

Table 3 .
The time convergence rate with h = 1/2 6 at T = 1 in Example 5.3