Analysis of a finite difference scheme for a nonlinear Caputo fractional differential equation on an adaptive grid

Abstract: A nonlinear initial value problem whose differential operator is a Caputo derivative of order α with 0 < α < 1 is studied. By using the Riemann-Liouville fractional integral transformation, this problem is reformulated as a Volterra integral equation, which is discretized by using the right rectangle formula. Both a priori and an a posteriori error analysis are conducted. Based on the a priori error bound and mesh equidistribution principle, we show that there exists a nonuniform grid that gives first-order convergent result, which is robust with respect to α. Then an a posteriori error estimation is derived and used to design an adaptive grid generation algorithm. Numerical results complement the theoretical findings.


Introduction
In this paper, we consider the following nonlinear Caputo fractional initial value problem: u(t)) = 0, t ∈ Ω = (0, 1), where α ∈ (0, 1), D α 0 is a caputo fractional order derivative operator, which is defined by (1.2) and f (t, u(t)) ∈ C 1 (Ω × R) is a known function. Furthermore, throughout this paper we shall assume that there exist two positive constants β i (i = 1, 2) such that It is well known that fractional-order differential equations (FDEs) have widely been applied in many fields like bio-engineering [1], food science [2], electrical engineering [3], Biomedicine [4][5][6], epidemiology [7], etc. Due to the importance of these problems, there has been tremendous interest in developing numerical methods for FDEs, such as finite difference method [8][9][10][11], finite element method [12], collocation method [13][14][15] and spectral method [16,17]. However, many of these methods are based on local operations and the authors don't consider the weak singularities. For this reason, Stynes and Gracia [18,19] first considered a two point boundary value problem with a Caputo derivative of order δ ∈ (1, 2), and analyzed the bounds for the derivatives of exact solution u(x). Then they constructed finite difference methods on a uniform mesh to solve this problem and gave the corresponding convergence analysis.
As far as we known, there is great current interest in the use of some numerical methods on a graded grid for time fractional diffusion equations with a Caputo derivative (see, e.g., [20][21][22][23]). In these literature, the authors derived the bounds for the time derivatives of u(x, t), and constructed a finite difference scheme on a graded grid. Recently, for a two-point boundary value problem with a Riemann-Liouville fractional derivative, Cen, et.al., [24] developed an adaptive grid method based on an a posteriori error estimation. Furthermore, Liu, et.al., [25] and Huang,et.al.,[26] proposed an adaptive grid method based on a priori error estimation and obtained a first-order and a second-order convergence results, respectively. Inspired by literature [24][25][26], the authors in [27] also developed an adaptive grid method by using a backward Euler formula to approximate the first-order derivative of problem (1.1) and obtained an a posteriori error estimation for the presented discretization scheme. Liu and Chen [28] derived a new a posteriori error estimation and gave the corresponding adaptive strategy for problem (1.1).
In this paper, similar to literature [25,26], we first transform problem (1.1) into a Volterra integro equation by using the Riemann-Liouville fractional integral transformation. Then a right rectangle formula on an arbitrary nonuniform mesh is used to approximate this integro equation. Furthermore, a convergence analysis based on a priori error estimation is carried out by using the mesh equidistribution principle. It is shown that there exists a grid that gives the optimal first-order convergence for the presented method. Besides, we also derive an a posteriori error estimation for the presented discretization scheme and design an adaptive grid generation algorithm. Finally, some numerical results are provided to validate the theoretical results.
Notation. Throughout the paper, C will denote a generic positive constant that is independent of the mesh parameter N and α. To simplify the notation we set g i = g(t i ) for any function g defined on the interval [0, 1]. We define the maximum norm by · ∞ .

Preliminary results
In Section 2.1, we first convert problem (1.1) into an equivalent Volterra integral equation and give the bounds for the exact solution u(t) and its derivatives. Then for the transformed Volterra integral equation, a finite difference scheme and the corresponding stability result are given in Section 2.2.

Reformulation the initial value problem
be the Riemann-Liouville fractional integral of order α ∈ (0, 1) (see [29]), which satisfies the following property Then, the above problem (1.1) can be written into the following Volterra integral equation Similar to the argument of Lemma 2.1 of [24], we derive the following results for the exact solution u(t) and its derivatives.
Proof. By using the Taylor series expansion to f (t, u(t)), problem (2.2) can be changed into the following linear formal where a(s) = f u (s, λu(s)) with 0 < λ < 1. Then, the proof of (2.3) can be followed from (2.5) and Lemma 2.1 of [24]. For the proof of (2.4), we may refer to the Theorem 2.1 in [20].

Discretization scheme
For our numerical scheme we consider an arbitrary nonuniform mesh where N is a positive integer. For i = 1, · · · , N, let τ i = t i − t i−1 be the local mesh step. Then, by using the right rectangle formula to approximate the integral given in (2.2), we obtain where u N i is the approximation solution of u(t) at point t = t i , i = 0, 1, · · · , N. Next, the following lemma gives the stability for the discretization scheme (2.6).

Lemma 2.2.
Under the assumption f (t, u(t)) ∈ C 1 (Ω × R), the solution u N i N i=0 of the scheme (2.6) on an arbitrary meshΩ N satisfies Proof. Similar to (2.5), the numerical scheme (2.6) can be written into the following formula where ξ ∈ (0, 1). Furthermore, we have Then, for each i, it follows from the proof of Lemma 3.3 given in [30] that (2.14) Furthermore, we have Finally, by using the modified Grönwall inequality given in [30,Lemma 3.3], we can obtain the desirable result (2.7).

A priori error analysis
Let e N i = u(t i ) − u N i be the error at t i in the computed solution, i = 0, 1, · · · , N. Then it follows from (2.2) and (2.6) that we can obtain the following error equation is the truncation error at t = t i . be the discrete solution of (2.6) on an adaptive grid. Then Thus, the desired result can be followed from (3.1) and Lemma 2.2.
Remark 3.1. It is shown from Corollary 3.1 that there exists a mesh that gives the optimal first-order convergence of the presented discretization scheme (2.6). However, it is difficult to construct this mesh {t i } N i=0 based on (3.3) since this requires the exact solution u(t).

A posteriori error analysis
In this section, we will derive an a posteriori error estimation for the numerical solution u N i , i = 0, 1, · · · , N. Theorem 4.1. Let u(t) be the exact solution of (1.1), u N i be the solution of (2.6) and u N (t) be the piecewise linear interpolation function through knots t i , u N i , i = 0, 1, · · · , N. Then we have Proof. For ∀t ∈ [t i−1 , t i ], it follows from (2.2) and (2.6) that For p, we have u N (s)) .

(4.7)
It follows from the assumption of f (t, u(t)) that Thus, from (4.2), (4.3) and (4.7)-(4.9), we obtain where we have used the bound of u N i and the condition (1.3). Finally, the desired result can be followed by applying the generalized Grönwall's inequality [31, Corollary 2] to the above inequality. Under the assumption of f (t, u(t)) ∈ C 1 (Ω × R), we have which completes the proof by virtue of the above Theorem 4.1.

Numerical experiments and discussion
In Section 5.1, we describe an adaptive grid generation algorithm by equidistributing a monitor function. Then, numerical experiments are presented in Section 5.2 to demonstrate the validity and efficiency of the presented adaptive grid method.

An adaptive algorithm
As is stated in Remark 3.1 that it is hard to obtain a grid {t i } N i=0 satisfying (3.3). Therefore, in practical computation, the key problem is to find a grid {t i } N i=0 and the corresponding numerical solution u N i such that where M(s, u N i ) is a monitor function which is a function about u N i . For a given monitor function M, the adaptive grid generation algorithm based on (5.1) aims to construct a mesh that equidistributes M. Thus, it is very important to choose a suitable monitor function M. Here, in this paper, by using the a posteriori error estimation given in Theorem 4.1, we construct the monitor function as follows: In order to compute the equidistributed mesh {t i } N i=0 and corresponding numerical solution u N i , an adaptive grid algorithm is given as follows: Step 1.
be an initial uniform mesh with mesh size 1 N .
Step 2. For a given meshΩ (k) j is the value of the monitor function (5.2) computed at the current mesh and corresponding numerical solution.
Step 4.Let Y (k) i = iL (k) N /N and φ (k) (t) be a piecewise linear interpolation function through knots i for i = 0, · · · , N. Let k = k + 1 and return to Step 2.

Two test examples
Here, we give some numerical experiments to illustrate the validity of our presented adaptive grid method. The test problem follows [27] by taking (1.1) with f (t, u(t)) = 2u(t) + sin(u(t)) + 0.1u 2 (t) + s(t), where s(t) is chosen such that the exact solution is The initial value condition is u(0) = 1.
Since the analytic solution of this problem is given, the maximum point-wise error E N and the order of convergence r N are calculated as follows: For different values of N and α, we use our presented adaptive grid algorithm to solve this test problem. The maximum errors E N , the orders of convergence r N and the number of iterations K of the above grid generation algorithm are listed in Table 1. One can see from the numerical results in Table 1 that our presented adaptive grid method is first-order convergent. which is robust with respect to the order of fractional derivative α. Meanwhile, to illustrate the advantage of our presented adaptive grid method, we also use the presented discretization scheme (2.6) on a uniform mesh to solve this test problem, see Table 2.
It is shown from these numerical results that the maximum errors calculated on an adaptive grid are much lower than that computed on a uniform mesh with the increase of N. Besides, the order of convergence obtained on an adaptive grid is more higher than that obtained on a uniform mesh.
Furthermore, in order to verify the relationship between numerical solution u N i and the order of fractional derivative α, Figure 1 represents some graphs of numerical solution for different values of N and α. Obviously, the solution of this test problem has a boundary layer at t = 0 with the decrease of α. When α = 0.1, the Figure 2 shows how a mesh with N = 64 intervals evolves through successive iterations of the above grid generation algorithm.
Finally, for α = 0.2, 0.4, 0.6, 0.8 and the same N given in Table 1, Table 3 gives the numerical results calculated by using our presented adaptive grid method, while we also list the results obtained by the method in [28]. Obviously, it is shown from Table 3 that our presented method produces better results than that produced by method in [28].

Conclusions
This work mainly discusses a nonlinear value problem whose the differential operator is a Caputo derivative of order α with 0 < α < 1. By using the Riemann-Liouville integral operator, the problem (1.1) can be changed into a Volterra integral equation, which is approximated by using the integral discrete formula. A priori error estimation and convergence analysis have been given on an adaptive grid. Meanwhile, an a posterior error estimation has been obtained by using the polynomial interpolation technique and the corresponding adaptive grid generation algorithm is constructed. It should be pointed out that the presented adaptive grid method can be extended to the other time fractional differential equations.