A finite difference method for a two-point boundary value problem with a Caputo fractional derivativeThis research was partly supported by the Instituto Universitario de Matemáticas y Aplicaciones (IUMA), the Ministerio de Educación, Cultura y Deporte (Programa Nacional de Movilidad de Recursos Humanos del Plan Nacional de I+D+i 2008-2011), project MEC/FEDER MTM 2010-16917 and the Diputación General de Aragón.

# A finite difference method for a two-point boundary value problem with a Caputo fractional derivative1

## Abstract

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.

## 1 Introduction

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

 Dng(x)=(ddx)m[1Γ(m−n)∫xt=0(x−t)m−n−1g(t)dt]for  0

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

 Dn∗g=Dn[g−Tm−1[g;0]], (1.1)

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

 Dn∗g(x):=1Γ(m−n)∫xt=0(x−t)m−n−1g(m)(t)dtfor  0

Our work relies heavily on Pedas and Tamme (2012), who use the definition (1.1) of .

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

 −Dδ∗ u(x)+b(x)u′(x)+c(x)u(x)=f(x)  for x∈(0,1), (1.3a) u(0)−α0u′(0)=γ0,u(1)+α1u′(1)=γ1, (1.3b)

where the constants and the functions and are given. We assume that and

 α0≥1δ−1. (1.4)

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

 Am[0,1]={g∈Cm−1[0,1]:g(m−1) is absolutely continuous on [0,1]}.

In several inequalities denotes a generic constant that depends on the data of the boundary value problem (1.3) but is independent of any mesh used to solve (1.3) numerically; note that can take different values in different places.

## 2 Comparison principle

We begin with a basic result.

###### Lemma 2.1.

(Al-Refai, 2012b, Theorem 2.1) Let achieve a global minimum at . Then

 Dδ∗g(x0)≥x−δ0Γ(2−δ){(δ−1)[g(0)−g(x0)]−x0g′(0)}. (2.1)

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

 (Regularity Hypothesis 1)g∈C2(0,1] and |g′′(x)|≤Cx−θ for 0

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.

###### Lemma 2.2.

(Al-Refai, 2012a, Lemma 3.3) Let . Let with for all . Assume that satisfies the inequalities

 −Dδ∗g+bg′+cg≥0  on (0,1), (2.2a) g(0) −α0g′(0)≥0  and  g(1)+α1g′(1)≥0, (2.2b)

where satisfies (1.4) and . Then on .

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.

###### Theorem 2.1.

Let satisfy Regularity Hypothesis 1. Let with for all . Assume that satisfies the inequalities (2.2), where satisfies (1.4) and . Then on .

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.

###### Example 2.1.

Take . From Diethelm (2010, Appendix B) we have

 D1.2∗x2=Γ(3)Γ(1.8)x0.8  and  D1.2∗x3=Γ(4)Γ(2.8)x1.8=3Γ(3)(1.8)Γ(1.8)x1.8≤(1.67)Γ(3)Γ(1.8)x0.8

for . Also . Set . Then

 −D1.2∗g(x)≥0  for  0

but , so the Dirichlet boundary conditions in (2.3) do not justify a comparison principle for on .

In this example one has and , so the condition of (2.2b) is not satisfied.

## 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

 |y(i)(x)|≤⎧⎨⎩Cif i<1−ν,C(1+|lnx|)if i=1−ν,Cx1−ν−iif i>1−ν, (3.1)

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.

###### Theorem 3.1.

(Pedas and Tamme, 2012, Theorem 2.1) Let for some integer and . Set . Let denote the set of functions defined on for which . Assume that the problem (1.3) with and has in only the trivial solution . Then (1.3) has a unique solution ; furthermore, .

###### Remark 3.1.

In Pedas and Tamme (2012, Theorem 2.1) there is the additional assumption that the only linear polynomial that satisfies the boundary conditions (1.3b) is , but it is straightforward to check that this condition is implied by  and (1.4).

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

 ddx[∫xs=0(x−s)θ1r(s)ds]=∫xs=0θ1(x−s)θ1−1r(s)dsfor  0

when for and the constants lie in ; one can justify (3.2) from the results of Talvila (2001) or by writing

 ∫xs=0(x−s)θ1r(s)ds=∫xs=0r(s)∫xt=sθ1(t−s)θ1−1dtds=∫xt=0∫ts=0θ1(t−s)θ1−1r(s)dsdt

then applying the fundamental theorem of calculus.

The first step is the following technical result.

###### Lemma 3.1.

Let be a positive integer and let satisfy . Suppose that

 w(x)=∫xs=0(x−s)σ−mψ(s)dsfor  0

where with for and some constant . Then

 w′(x) =1x∫xs=0(x−s)σ−m[sψ′(s)+(σ−m+1)ψ(s)]ds (3.3) and |w′(x)| ≤C1β(σ−m+1,σ−m+1)x2(σ−m)% for  0

where is Euler’s Beta function.

###### Proof.

For ,

 xw(x) =∫xs=0[(x−s)σ−m+1+(x−s)σ−ms]ψ(s)ds =∫xs=0{(x−s)σ−m+1ψ(s)+(x−s)σ−m+1σ−m+1[sψ′(s)+ψ(s)]}ds,

after an integration by parts. Applying (3.2) one gets

 (xw(x))′ =∫xs=0(x−s)σ−m[sψ′(s)+(σ−m+2)ψ(s)]ds, and hence w′(x) =1x[(xw(x))′−w(x)]=1x∫xs=0(x−s)σ−m[sψ′(s)+(σ−m+1)ψ(s)]ds,

as desired. Furthermore, the hypotheses of the lemma imply that

 |w′(x)|≤C1x∫xs=0(x−s)σ−msσ−mds=C1β(σ−m+1,σ−m+1)x2(σ−m),

where the value of the integral is given by Euler’s Beta function (Diethelm, 2010, Theorem D.6). ∎

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.

###### Theorem 3.2.

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

 |r(i)(x)|≤{Cif  i=0,1,…,m−1,Cxσ−iif  i=m,m+1,…,q+m−1. (3.5)
###### Proof.

First, implies (3.5) for .

Next, we show that . Set . Then and . By definition

 Dσw(x) =(ddx)m[1Γ(m−σ)∫xt=0(x−t)m−σ−1w(t)dt] =(ddx)m[1Γ(2m−σ−1)∫xt=0(x−t)2m−σ−2w(m−1)(t)dt] =ddx[1Γ(m−σ)∫xt=0(x−t)m−σ−1w(m−1)(t)dt],

after integrations by parts followed by differentiations using (3.2). Consequently

 ∫xs=0Dσw(s)ds=1Γ(m−σ)∫xt=0(x−t)m−σ−1w(m−1)(t)dt.

This is an Abel integral equation for the function . Thus from Samko et al. (1993, Section 2) it follows that

 w(m−1)(x)=ddx{1Γ(σ+1−m)∫xt=0(x−t)σ−m[∫ts=0Dσw(s)ds]dt}.

But by hypothesis, so we can integrate by parts then use (3.2) to get

 w(m−1)(x) =ddx[1Γ(σ+2−m)∫xt=0(x−t)σ+1−mDσw(t)dt] =1Γ(σ+1−m)∫xt=0(x−t)σ−mDσw(t)dt.

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

 r(x)=m−1∑j=0r(j)(0)j!xj+1Γ(σ)∫xs=0(x−s)σ−1ϕ(s)ds.

Integration by parts yields

 r(x) =m−1∑j=0r(j)(0)j!xj+ϕ(0)σΓ(σ)xσ+1σΓ(σ)∫xs=0(x−s)σϕ′(s)ds  for  0≤x≤1.

Differentiating this formula times using and (3.2), we obtain

 r(m)(x)=ϕ(0)Γ(σ−m+1)xσ−m+1Γ(σ−m+1)∫xs=0(x−s)σ−mϕ′(s)ds. (3.6)

Hence, since , for some constants one obtains

 |r(m)(x)| ≤|ϕ(0)|Γ(σ−m+1)xσ−m+CΓ(σ−m+1)∫xs=0(x−s)σ−msσ−mds ≤Cxσ−m+Cx2(σ−m)+1 ≤Cxσ−m (3.7)

as , and Diethelm (2010, Theorem D.6) was invoked to bound the integral (Euler’s Beta function). This is the desired bound (3.5) for . Furthermore, it is easy to see from (3.6) that .

We now deduce (3.5) for from (3.6). Applying Lemma 3.1 with to differentiate (3.6), one gets

 |r(m+1)(x)|≤C[xσ−m−1+x2(σ−m)]≤Cxσ−m−1,

which proves (3.5) for , and

 r(m+1)(x)=ϕ(0)Γ(σ−m)xσ−m−1+1Γ(σ−m+1)⋅1x∫xs=0(x−s)σ−m[sϕ′′(s)+(σ−m+1)ϕ′(s)]ds, (3.8)

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

 1x⋅O(x2(σ−m))+1x2⋅O(x2(σ−m)+1),

whence

 |r(m+2)(x)|≤C[xσ−m−2+x2(σ−m)−1]≤Cxσ−m−2,

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).

###### Corollary 3.1.

Let for some integer and . Assume that and the condition (1.4) is satisfied. Then (1.3) has a unique solution with , and for all there exists a constant  such that

 |u(i)(x)|≤{Cif  i=0,1,Cxδ−iif  i=2,3,…,q+1. (3.9)
###### Proof.

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. ∎

###### Remark 3.2.

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

 u′′(x)=1b(x)[f′+(Dδ∗u)′−c′u−(c+b′)u′](x), (3.10)

and an appeal to the bounds of Theorem 3.1 yields (3.9) immediately for the case . One can then differentiate (3.10) iteratively and use Theorem 3.1 to prove (3.9) for .

If instead and on , a similar technique will work—note that enables the conclusion of Theorem 3.1 to be strengthened to where by Pedas and Tamme (2012, Remark 2.2).

Finally, we give an example to show that the bounds of Theorem 3.2 are sharp.

###### Example 3.1.

Let be a positive integer with . Set for . Clearly . Then from (Diethelm, 2010, Appendix B) one gets

 Dσ∗r(x)=Γ(σ+1)+Γ(2σ−m+2)Γ(σ−m+2)xσ−m+1.

Hence and , while

 (Dσ∗r)(i)(x) =Γ(2σ−m+2)Γ(σ−m+2−i)xσ−m+1−i and r(i)(x) =Γ(σ+1)Γ(σ+1−i)xσ−i+Γ(2σ−m+2)Γ(2σ−m+2−i)x2σ−m+1−i  for  i=1,2,…

Thus the derivatives of satisfy the hypotheses of Theorem 3.2 and the derivatives of agree with (3.5) for , i.e., the outcome of Theorem 3.2 cannot be sharpened for .

## 4 Discretization and convergence

### 4.1 The discretization of the boundary value problem

Assume the hypotheses of Corollary 3.1. Let be a positive integer. Subdivide by the uniform mesh , for . Then the standard discretization of for is (see, e.g., Sousa, 2012) given by

 −Dδ∗u(xj)≈−1hδΓ(3−δ)j−1∑k=0dj−k(uk+2−2uk+1+uk), (4.1)

where denotes the computed approximation to , and we set

 dr=r2−δ+−(r−1)2−δ+  for all integers r, (4.2)

with

 s+={sif s≥0,0if s<0.

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

 (bu′)(xj)≈bjDℵuj:={bj(uj+1−uj)/hif bj<0,bj(uj−uj−1)/hif bj≥0. (4.3)

This difference approximation can also be written as

 bjDℵuj=−(bj+|bj|)uj−12h+|bj|ujh+(bj−|bj|)uj+12h. (4.4)

The full discretization of (1.3a) is

 −1hδΓ(3−δ)j−1∑k=0dj−k(uk+2−2uk+1+uk)+bjDℵuj+cjuj=fj, (4.5)

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 row is
. For , the entries of the row of satisfy

 aj0= −djhδΓ(3−δ)−ϵj1b1+|b1|2h, (4.6a) aj1= −dj−1+2djhδΓ(3−δ)−ϵj2b2+|b2|2h+ϵj1(|b1|h+c1), (4.6b) ajk= −dj−k+2dj−k+1−dj−k+2hδΓ(3−δ)−ϵj,k+1bj+|bj|2h+ϵjk(|bj|h+cj)+ϵj,k−1bj−|bj|2h %fork=2,3,…,N, (4.6c)

where we set

 ϵjk={1if j=k,0otherwise.

Hence the matrix is lower Hessenberg.

Observe that (4.5) implies that

 N∑k=0ajk=cj  for j=1,2,…,N−1. (4.7)

We shall prove various inequalities for the non-zero entries of . First, for all  by (4.6), (4.2) and . From (4.6a), one has

 aj0<0for j=1,2,…,N−1. (4.8)

By (4.6b),

 a21=23−δ−3hδΓ(3−δ)−b2+|b2|2h, (4.9)

so the sign of depends on and .

One has for .

###### Proof.

For , equations (4.6b) and (4.2) yield