A finite difference method for a two-point boundary value problem with a Caputo fractional derivative
A two-point boundary value problem whose highest-order term is a Caputo fractional derivative of order is considered. Al-Refai’s comparison principle is improved and modified to fit our problem. Sharp a priori bounds on derivatives of the solution of the boundary value problem are established, showing that may be unbounded at the interval endpoint . These bounds and a discrete comparison principle are used to prove pointwise convergence of a finite difference method for the problem, where the convective term is discretized using simple upwinding to yield stability on coarse meshes for all values of . Numerical results are presented to illustrate the performance of the method. Fractional differential equation; Caputo fractional derivative; boundary value problem; derivative bounds; finite difference method; convergence proof.
Fractional derivatives are used in an ever-widening range of models of physical processes, and as a consequence the last decade has seen an explosive growth in the number of numerical analysis papers examining differential equations with fractional-order derivatives (see the references in Machado et al., 2011). While the analysis of some of these papers (e.g., Mustapha and McLean, 2012; Pedas and Tamme, 2012) takes account of the possibly singular behaviour of solutions near some domain boundaries, most fractional-derivative numerical analysis papers work only with very special cases by assuming (explicitly or implicitly) that the solutions they approximate are smooth on the closure of the domain where the problem is posed. In particular, we know of no paper where a finite difference method for a fractional-derivative boundary value problem posed on a bounded domain is analysed rigorously under reasonably general and realistic hypotheses on the behaviour of the solution near the boundaries of that domain. In the present paper we provide the first such rigorous analysis.
Even though we deal with the one-dimensional case—a two-point boundary value problem—the analysis is nevertheless lengthy and requires the development of various techniques that do not appear in the context of “classical” problems (i.e., problems with integer-order derivatives).
Let satisfy for some positive integer . The Riemann-Liouville fractional derivative is defined by
for all functions such that exists. Our interest centres on the Caputo fractional derivative , which is defined (Diethelm, 2010, Definition 3.2) in terms of by
where denotes the Taylor polynomial of degree of the function expanded around . If and is absolutely continuous on , then (Diethelm, 2010, Theorem 3.1) one also has the equivalent formulation
Since the integrals in and are associated in a special way with the point , many authors write instead and , but for simplicity of notation we omit the extra subscript .
Let the parameter satisfy . Throughout the paper we consider the two-point boundary value problem
where the constants and the functions and are given. We assume that and
The condition (1.4) comes from Al-Refai (2012a); it will be used in Sections 2 and 4 below to ensure that (1.3) and its discretization each satisfy a suitable comparison principle. For the moment we assume that ; further hypotheses will be placed later on the regularity of these functions.
The problem (1.3) is discussed by Al-Refai (2012a) and is a particular case of the wide class of boundary value problems considered in Pedas and Tamme (2012). It is a steady-state version of the time-dependent problems discussed in Saadatmandi and Dehghan (2011); Shen and Liu (2004/05); Zheng et al. (2010) and Ji and Tang (2012)—who describe some advantages of the Caputo fractional derivative over the Riemann-Liouville fractional derivative.
Our paper is structured as follows. Section 2 obtains a comparison principle for the differential operator and boundary operators in (1.3). In Section 3 existence and uniqueness of a solution to (1.3) is shown, and sharp pointwise bounds on the integer-order derivatives of this solution are derived. The finite difference discretization of (1.3) on a uniform mesh of width is described and analysed in Section 4, and it is proved to be convergent at the mesh points. Two numerical examples are presented in Section 5.
Notation. We use the standard notation to denote the space of real-valued functions whose derivatives up to order are continuous on an interval , and write for . For each , set . As in Diethelm (2010), for each positive integer define
2 Comparison principle
We begin with a basic result.
(Al-Refai, 2012b, Theorem 2.1) Let achieve a global minimum at . Then
A careful inspection of the argument used to prove this lemma in Al-Refai (2012b) shows that it remains valid under the weaker regularity hypothesis that
where and are some fixed constants with . Observe that any function that satisfies Regularity Hypothesis 1 can be extended to a function (which we also call ) lying in . We shall see in Section 3 that the solution of the boundary value problem (1.3) satisfies Regularity Hypothesis 1 but does not in general lie in .
Lemma 2.1 is the key tool needed to prove the following comparison principle.
Recalling our observation above that Lemma 2.1 is still true when the hypothesis is replaced by Regularity Hypothesis 1, one sees quickly from the proof of Al-Refai (2012a, Lemma 3.3) that Lemma 2.2 remains valid when the assumption is replaced by Regularity Hypothesis 1. In fact, one can go further: Lemma 2.1 shows immediately that when one has at the global minimum, and invoking this observation in the proof of Al-Refai (2012a, Lemma 3.3) and changing a few inequalities there from strict to weak or vice versa, the hypothesis can be weakened to . That is, one has the following more general version of Lemma 2.2.
The next example shows that, for our Caputo differential operator, one does not have a comparison principle for the simplest case of Dirichlet boundary conditions, unlike the situation for classical second-order boundary value problems. Thus one cannot permit in Theorem 2.1.
3 A priori bounds on derivatives of the solution
The only source we know for bounds on (certain) derivatives of the solution of (1.3) is Pedas and Tamme (2012), who prove a very general existence result for two-point boundary value problems with differential operators involving fractional-order derivatives. Their analysis is based on Brunner et al. (2001).
Notation. For integer and , define to be the set of continuous functions that are times continuously differentiable in and satisfy the bounds
for and , where is some constant. Observe that as increases, the smoothness of functions in decreases. Clearly for and .
For our problem (1.3), the Pedas and Tamme result is as follows.
While Theorem 3.1 bounds and the integer-order derivatives of , it gives no bound on the derivatives for , but these derivatives will be needed in the consistency analysis of our finite difference method. Thus we now deduce bounds on the integer-order derivatives of from Theorem 3.1. Our elementary argument can be regarded as interpolating between the integer-order derivatives of ; it relies only on the derivative bounds stated in Theorem 3.1 and makes no use of the differential equation (1.3a).
We shall prove this bound on the integer-order derivatives in a general setting that is suited to fractional-derivative boundary value problems of arbitrary order—such as those considered in Pedas and Tamme (2012)—since the general proof is essentially the same as the proof for . At various places in our calculations we shall need the formula
then applying the fundamental theorem of calculus.
The first step is the following technical result.
Let be a positive integer and let satisfy . Suppose that
where with for and some constant . Then
where is Euler’s Beta function.
The essential property of the bound (3.4) is that it takes the form with a constant that is independent of .
Now we can proceed with the main result of this section.
Let be a positive integer with . Assume that and for some integer . Then and for all there exists a constant , which is independent of , such that
First, implies (3.5) for .
Next, we show that . Set . Then and . By definition
after integrations by parts followed by differentiations using (3.2). Consequently
This is an Abel integral equation for the function . Thus from Samko et al. (1993, Section 2) it follows that
But by hypothesis, so we can integrate by parts then use (3.2) to get
As the integrand here lies in the space of Lebesgue integrable functions, it follows that is absolutely continuous on . Hence is absolutely continuous on , i.e., .
We come now to the main part of the proof. Set . As , by Diethelm (2010, Corollary 3.9) one has
Integration by parts yields
Differentiating this formula times using and (3.2), we obtain
Hence, since , for some constants one obtains
which proves (3.5) for , and
from which one can see that .
Comparing (3.6) and (3.8), the relationship between their leading terms is simple, while the integrals in both are but the integral in (3.8) is multiplied by . One now proceeds to differentiate (3.8), invoking Lemma 3.1 with ; this will yield a rather complicated formula for that involves two integrals, but one sees readily that these integrals are
which proves (3.5) for . Continuing in this way, each higher derivative of introduces a further factor in the estimates, and we can derive successively the bounds of (3.5) for . The calculation must stop when one reaches an integral involving , i.e., when . ∎
We now apply this result to our boundary value problem (1.3).
Observe that any function in satisfies Regularity Hypothesis 1 of Section 2 since . Consequently Theorem 2.1 implies that if and , then the problem (1.3) has in only the trivial solution . Hence Theorem 3.1 yields existence and uniqueness of a solution of (1.3) with and . An appeal to Theorem 3.2 completes the proof. ∎
If we impose the additional hypothesis that on , it is then possible to give a simple proof of Corollary 3.1 directly from Theorem 3.1 without using Theorem 3.2: differentiating (1.3a) then solving for yields
Finally, we give an example to show that the bounds of Theorem 3.2 are sharp.
4 Discretization and convergence
4.1 The discretization of the boundary value problem
where denotes the computed approximation to , and we set
Note that for .
Set for each mesh point , where can be or . To discretize the convective term we shall use simple upwinding (Roos et al., 2008, p.47), because the standard approximation may yield a non-monotone difference scheme when is near 1; see Gracia and Stynes (2013) for details. Thus we use the approximation
This difference approximation can also be written as
The full discretization of (1.3a) is
for . The boundary conditions (1.3b) are discretized by approximating by and by .
Let denote the matrix corresponding to this discretization of (1.3), i.e., where
and the superscript denotes transpose.
Thus the row of is and its
. For , the entries of the row of satisfy
where we set
Hence the matrix is lower Hessenberg.
Observe that (4.5) implies that
so the sign of depends on and .
One has for .