A finite difference method for a two-point boundary value problem with a Caputo fractional derivative
aa r X i v : . [ m a t h . NA ] D ec A finite difference method for a two-point boundary valueproblem with a Caputo fractional derivative ∗ Martin Stynes † Department of Mathematics, National University of Ireland, Cork, IrelandandJos´e Luis Gracia ‡ IUMA and Department of Applied Mathematics, University of Zaragoza, SpainDecember 19, 2013
Abstract
A two-point boundary value problem whose highest-order term is a Caputo fractionalderivative of order δ ∈ (1 ,
2) is considered. Al-Refai’s comparison principle is improved andmodified to fit our problem. Sharp a priori bounds on derivatives of the solution u of theboundary value problem are established, showing that u ′′ ( x ) may be unbounded at the in-terval endpoint x = 0. These bounds and a discrete comparison principle are used to provepointwise convergence of a finite difference method for the problem, where the convectiveterm is discretized using simple upwinding to yield stability on coarse meshes for all valuesof δ . Numerical results are presented to illustrate the performance of the method. Frac-tional differential equation; Caputo fractional derivative; boundary value problem; derivativebounds; finite difference method; convergence proof. Fractional derivatives are used in an ever-widening range of models of physical processes, and asa consequence the last decade has seen an explosive growth in the number of numerical analysispapers examining differential equations with fractional-order derivatives (see the references inMachado 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 nearsome domain boundaries, most fractional-derivative numerical analysis papers work only withvery special cases by assuming (explicitly or implicitly) that the solutions they approximate aresmooth on the closure of the domain where the problem is posed. In particular, we know of nopaper where a finite difference method for a fractional-derivative boundary value problem posedon a bounded domain is analysed rigorously under reasonably general and realistic hypotheseson the behaviour of the solution near the boundaries of that domain. In the present paper weprovide the first such rigorous analysis. ∗ This research was partly supported by the Instituto Universitario de Matem´aticas y Aplicaciones (IUMA),the Ministerio de Educaci´on, Cultura y Deporte (Programa Nacional de Movilidad de Recursos Humanos del PlanNacional de I+D+i 2008-2011), project MEC/FEDER MTM 2010-16917 and the Diputaci´on General de Arag´on. † Corresponding author. Email: [email protected] ‡ Email: [email protected] n ∈ R satisfy m − < n < m for some positive integer m . The Riemann-Liouvillefractional derivative D n is defined by D n g ( x ) = (cid:18) ddx (cid:19) m (cid:20) m − n ) Z xt =0 ( x − t ) m − n − g ( t ) dt (cid:21) for 0 < x ≤ g such that D n g ( x ) exists. Our interest centres on the Caputo fractionalderivative D n ∗ , which is defined (Diethelm, 2010, Definition 3.2) in terms of D n by D n ∗ g = D n [ g − T m − [ g ; 0]] , (1.1)where T m − [ g ; 0] denotes the Taylor polynomial of degree m − g expandedaround x = 0. If g ∈ C m − [0 ,
1] and g ( m − is absolutely continuous on [0 , D n ∗ g ( x ) := 1Γ( m − n ) Z xt =0 ( x − t ) m − n − g ( m ) ( t ) dt for 0 < x ≤ . (1.2)Our work relies heavily on Pedas and Tamme (2012), who use the definition (1.1) of D n ∗ .Since the integrals in D n g ( x ) and D n ∗ g ( x ) are associated in a special way with the point x = 0, many authors write instead D n g ( x ) and D n ∗ g ( x ), but for simplicity of notation we omitthe extra subscript 0.Let the parameter δ satisfy 1 < δ <
2. Throughout the paper we consider the two-pointboundary value problem − D δ ∗ u ( x ) + b ( x ) u ′ ( x ) + c ( x ) u ( x ) = f ( x ) for x ∈ (0 , , (1.3a) u (0) − α u ′ (0) = γ , u (1) + α u ′ (1) = γ , (1.3b)where the constants α , α , γ , γ and the functions b, c and f are given. We assume that α ≥ α ≥ δ − . (1.4)The condition (1.4) comes from Al-Refai (2012a); it will be used in Sections 2 and 4 belowto ensure that (1.3) and its discretization each satisfy a suitable comparison principle. For themoment we assume that b, c, f ∈ C [0 , h is described and analysed in Section 4, and it is proved to be O ( h δ − ) convergent at the meshpoints. Two numerical examples are presented in Section 5. Notation.
We use the standard notation C k ( I ) to denote the space of real-valued functionswhose derivatives up to order k are continuous on an interval I , and write C ( I ) for C ( I ). For2ach g ∈ C [0 , k g k ∞ = max x ∈ [0 , | g ( x ) | . As in Diethelm (2010), for each positive integer m define A m [0 ,
1] = { g ∈ C m − [0 ,
1] : g ( m − is absolutely continuous on [0 , } . In several inequalities C denotes a generic constant that depends on the data of the boundaryvalue problem (1.3) but is independent of any mesh used to solve (1.3) numerically; note that C can take different values in different places. We begin with a basic result.
Lemma 2.1. (Al-Refai, 2012b, Theorem 2.1) Let g ∈ C [0 , achieve a global minimum at x ∈ (0 , . Then D δ ∗ g ( x ) ≥ x − δ Γ(2 − δ ) (cid:8) ( δ −
1) [ g (0) − g ( x )] − x g ′ (0) (cid:9) . (2.1)A careful inspection of the argument used to prove this lemma in Al-Refai (2012b) showsthat it remains valid under the weaker regularity hypothesis that(Regularity Hypothesis 1) g ∈ C (0 ,
1] and | g ′′ ( x ) | ≤ Cx − θ for 0 < x ≤ , where C and θ are some fixed constants with 0 < θ <
1. Observe that any function g thatsatisfies Regularity Hypothesis 1 can be extended to a function (which we also call g ) lying in C [0 , ∩ A [0 , C [0 , Lemma 2.2. (Al-Refai, 2012a, Lemma 3.3) Let g ∈ C [0 , . Let b, c ∈ C [0 , with c ( x ) > for all x ∈ (0 , . Assume that g satisfies the inequalities − D δ ∗ g + bg ′ + cg ≥ on (0 , , (2.2a) g (0) − α g ′ (0) ≥ and g (1) + α g ′ (1) ≥ , (2.2b) where α satisfies (1.4) and α ≥ . Then g ≥ on [0 , . Recalling our observation above that Lemma 2.1 is still true when the hypothesis g ∈ C [0 , g ∈ C [0 ,
1] is replaced byRegularity Hypothesis 1. In fact, one can go further: Lemma 2.1 shows immediately thatwhen g ′ (0) < D δ ∗ g ( x ) > c > c ≥
0. That is, one has thefollowing more general version of Lemma 2.2.
Theorem 2.1.
Let g satisfy Regularity Hypothesis 1. Let b, c ∈ C [0 , with c ( x ) ≥ for all x ∈ (0 , . Assume that g satisfies the inequalities (2.2) , where α satisfies (1.4) and α ≥ .Then g ≥ on [0 , . α = 0 inTheorem 2.1. Example 2.1.
Take δ = 1 . . From Diethelm (2010, Appendix B) we have D . ∗ x = Γ(3)Γ(1 . x . and D . ∗ x = Γ(4)Γ(2 . x . = 3Γ(3)(1 . . x . ≤ (1 . . x . for < x < . Also D . ∗ x = 0 . Set g ( x ) = x − . x + 0 . x . Then − D . ∗ g ( x ) ≥ for < x < , g (0) = g (1) = 0 , (2.3) but g (0 .
8) = 0 . − . . < , so the Dirichlet boundary conditions in (2.3) do notjustify a comparison principle for − D . ∗ on [0 , .In this example one has g (0) = 0 and g ′ (0) > , so the condition g (0) − α g ′ (0) ≥ of (2.2b) is not satisfied. 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 dif-ferential operators involving fractional-order derivatives. Their analysis is based on Brunner et al.(2001).
Notation.
For integer q ≥ ν ∈ ( −∞ , C q,ν (0 ,
1] to be the set of continuousfunctions y : [0 , → R that are q times continuously differentiable in (0 ,
1] and satisfy thebounds | y ( i ) ( x ) | ≤ C if i < − ν,C (1 + | ln x | ) if i = 1 − ν,Cx − ν − i if i > − ν, (3.1)for 0 < x ≤ i = 1 , , . . . , q , where C is some constant. Observe that as ν increases, thesmoothness of functions in C q,ν (0 ,
1] decreases. Clearly C q [0 , ⊂ C q,ν (0 , ⊂ C m,µ (0 , ⊂ C [0 ,
1] for q ≥ m ≥ ν ≤ µ < Theorem 3.1. (Pedas and Tamme, 2012, Theorem 2.1) Let b, c, f ∈ C q,µ (0 , for some integer q ≥ and µ ∈ ( −∞ , . Set ν = max { µ, − δ } . Let S denote the set of functions w defined on [0 , for which D δ ∗ w ∈ C q,ν (0 , . Assume that the problem (1.3) with f ≡ , γ = 0 and γ = 0 has in S only the trivial solution w ≡ . Then (1.3) has a unique solution u ∈ S ; furthermore, u ∈ C [0 , . Remark 3.1.
In Pedas and Tamme (2012, Theorem 2.1) there is the additional assumptionthat the only linear polynomial y ( x ) that satisfies the boundary conditions (1.3b) is y ≡ , butit is straightforward to check that this condition is implied by α ≥ and (1.4) . While Theorem 3.1 bounds u ′ and the integer-order derivatives of D δ ∗ u , it gives no boundon the derivatives u ( i ) for i = 2 , , . . . , but these derivatives will be needed in the consistencyanalysis of our finite difference method. Thus we now deduce bounds on the integer-order4erivatives of u from Theorem 3.1. Our elementary argument can be regarded as interpolatingbetween the integer-order derivatives of D δ ∗ u ; it relies only on the derivative bounds stated inTheorem 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 suitedto fractional-derivative boundary value problems of arbitrary order—such as those consideredin Pedas and Tamme (2012)—since the general proof is essentially the same as the proof for D δ ∗ .At various places in our calculations we shall need the formula ddx (cid:20)Z xs =0 ( x − s ) θ r ( s ) ds (cid:21) = Z xs =0 θ ( x − s ) θ − r ( s ) ds for 0 < x < , (3.2)when | r ( s ) | ≤ Cs − θ for 0 < s ≤ θ , θ lie in (0 , Z xs =0 ( x − s ) θ r ( s ) ds = Z xs =0 r ( s ) Z xt = s θ ( t − s ) θ − dt ds = Z xt =0 Z ts =0 θ ( t − s ) θ − r ( s ) ds dt then applying the fundamental theorem of calculus.The first step is the following technical result. Lemma 3.1.
Let m be a positive integer and let σ ∈ R satisfy m − < σ < m . Suppose that w ( x ) = Z xs =0 ( x − s ) σ − m ψ ( s ) ds for < x < , where ψ ∈ C (0 , with | ψ ( s ) | + s | ψ ′ ( s ) | ≤ C s σ − m for < s < and some constant C . Then w ′ ( x ) = 1 x Z xs =0 ( x − s ) σ − m [ sψ ′ ( s ) + ( σ − m + 1) ψ ( s )] ds (3.3) and | w ′ ( x ) | ≤ C β ( σ − m + 1 , σ − m + 1) x σ − m ) for < x < , (3.4) where β ( · , · ) is Euler’s Beta function.Proof. For 0 < x < xw ( x ) = Z xs =0 [( x − s ) σ − m +1 + ( x − s ) σ − m s ] ψ ( s ) ds = Z xs =0 (cid:26) ( x − s ) σ − m +1 ψ ( s ) + ( x − s ) σ − m +1 σ − m + 1 [ sψ ′ ( s ) + ψ ( s )] (cid:27) ds, after an integration by parts. Applying (3.2) one gets (cid:0) xw ( x ) (cid:1) ′ = Z xs =0 ( x − s ) σ − m [ sψ ′ ( s ) + ( σ − m + 2) ψ ( s )] ds, and hence w ′ ( x ) = 1 x h(cid:0) xw ( x ) (cid:1) ′ − w ( x ) i = 1 x Z xs =0 ( x − s ) σ − m [ sψ ′ ( s ) + ( σ − m + 1) ψ ( s )] ds, as desired. Furthermore, the hypotheses of the lemma imply that | w ′ ( x ) | ≤ C x Z xs =0 ( x − s ) σ − m s σ − m ds = C β ( σ − m + 1 , σ − m + 1) x σ − m ) , where the value of the integral is given by Euler’s Beta function (Diethelm, 2010, Theorem D.6).5he essential property of the bound (3.4) is that it takes the form Cx σ − m ) with a constant C that is independent of x .Now we can proceed with the main result of this section. Theorem 3.2.
Let m be a positive integer with m − < σ < m . Assume that r ∈ C m − [0 , and D σ ∗ r ∈ C q,m − σ (0 , for some integer q ≥ . Then r ∈ C q + m − (0 , and for all x ∈ (0 , there exists a constant C , which is independent of x , such that | r ( i ) ( x ) | ≤ ( C if i = 0 , , . . . , m − ,Cx σ − i if i = m, m + 1 , . . . , q + m − . (3.5) Proof.
First, r ∈ C m − [0 ,
1] implies (3.5) for i = 0 , , . . . , m − r ∈ A m [0 , w = r − T m − [ r ; 0]. Then w ∈ C m − [0 ,
1] and0 = w (0) = w ′ (0) = · · · = w ( m − (0). By definition D σ w ( x ) = (cid:18) ddx (cid:19) m (cid:20) m − σ ) Z xt =0 ( x − t ) m − σ − w ( t ) dt (cid:21) = (cid:18) ddx (cid:19) m (cid:20) m − σ − Z xt =0 ( x − t ) m − σ − w ( m − ( t ) dt (cid:21) = ddx (cid:20) m − σ ) Z xt =0 ( x − t ) m − σ − w ( m − ( t ) dt (cid:21) , after m − m − Z xs =0 D σ w ( s ) ds = 1Γ( m − σ ) Z xt =0 ( x − t ) m − σ − w ( m − ( t ) dt. This is an Abel integral equation for the function w ( m − . Thus from Samko et al. (1993, Section2) it follows that w ( m − ( x ) = ddx (cid:26) σ + 1 − m ) Z xt =0 ( x − t ) σ − m (cid:20)Z ts =0 D σ w ( s ) ds (cid:21) dt (cid:27) . But D σ w = D σ ∗ r ∈ C [0 ,
1] by hypothesis, so we can integrate by parts then use (3.2) to get w ( m − ( x ) = ddx (cid:20) σ + 2 − m ) Z xt =0 ( x − t ) σ +1 − m D σ w ( t ) dt (cid:21) = 1Γ( σ + 1 − m ) Z xt =0 ( x − t ) σ − m D σ w ( t ) dt. As the integrand here lies in the space L [0 ,
1] of Lebesgue integrable functions, it follows that w ( m − is absolutely continuous on [0 , r ( m − = ( w + T m − [ r ; 0]) ( m − is absolutelycontinuous on [0 , r ∈ A m [0 , φ = D σ ∗ r . As r ∈ A m [0 , r ( x ) = m − X j =0 r ( j ) (0) j ! x j + 1Γ( σ ) Z xs =0 ( x − s ) σ − φ ( s ) ds . Integration by parts yields r ( x ) = m − X j =0 r ( j ) (0) j ! x j + φ (0) σ Γ( σ ) x σ + 1 σ Γ( σ ) Z xs =0 ( x − s ) σ φ ′ ( s ) ds for 0 ≤ x ≤ . m times using Γ( n + 1) = n Γ( n ) and (3.2), we obtain r ( m ) ( x ) = φ (0)Γ( σ − m + 1) x σ − m + 1Γ( σ − m + 1) Z xs =0 ( x − s ) σ − m φ ′ ( s ) ds. (3.6)Hence, since φ = D σ ∗ r ∈ C q,m − σ (0 , C one obtains | r ( m ) ( x ) | ≤ | φ (0) | Γ( σ − m + 1) x σ − m + C Γ( σ − m + 1) Z xs =0 ( x − s ) σ − m s σ − m ds ≤ Cx σ − m + Cx σ − m )+1 ≤ Cx σ − m (3.7)as x σ − m > x σ − m x σ +1 − m = x σ − m )+1 , and Diethelm (2010, Theorem D.6) was invoked to boundthe integral (Euler’s Beta function). This is the desired bound (3.5) for i = m . Furthermore, itis easy to see from (3.6) that r ( m ) ∈ C (0 , i = m + 1 , m + 2 , . . . from (3.6). Applying Lemma 3.1 with ψ ( s ) = φ ′ ( s ) to differentiate (3.6), one gets | r ( m +1) ( x ) | ≤ C [ x σ − m − + x σ − m ) ] ≤ Cx σ − m − , which proves (3.5) for i = m + 1, and r ( m +1) ( x ) = φ (0)Γ( σ − m ) x σ − m − + 1Γ( σ − m + 1) · x Z xs =0 ( x − s ) σ − m [ sφ ′′ ( s ) + ( σ − m + 1) φ ′ ( s )] ds, (3.8)from which one can see that r ( m +1) ∈ C (0 , O ( x σ − m )+1 ) but the integral in (3.8) is multiplied by 1 /x . One nowproceeds to differentiate (3.8), invoking Lemma 3.1 with ψ ( s ) = sφ ′′ ( s ) + ( σ − m + 1) φ ′ ( s ); thiswill yield a rather complicated formula for r ( m +2) ( x ) that involves two integrals, but one seesreadily that these integrals are1 x · O ( x σ − m ) ) + 1 x · O ( x σ − m )+1 ) , whence | r ( m +2) ( x ) | ≤ C [ x σ − m − + x σ − m ) − ] ≤ Cx σ − m − , which proves (3.5) for i = m + 2. Continuing in this way, each higher derivative of r introducesa further factor 1 /x in the estimates, and we can derive successively the bounds of (3.5) for i = m + 2 , m + 3 , . . . . The calculation must stop when one reaches an integral involving φ ( q ) ( s ),i.e., when i = q + m − Corollary 3.1.
Let b, c, f ∈ C q,µ (0 , for some integer q ≥ and µ ≤ − δ . Assume that c ≥ , α ≥ and the condition (1.4) is satisfied. Then (1.3) has a unique solution u with u ∈ C [0 , ∩ C q +1 (0 , , and for all x ∈ (0 , there exists a constant C such that | u ( i ) ( x ) | ≤ ( C if i = 0 , ,Cx δ − i if i = 2 , , . . . , q + 1 . (3.9)7 roof. Observe that any function in C q,µ (0 ,
1] satisfies Regularity Hypothesis 1 of Section 2since µ ≤ − δ <
1. Consequently Theorem 2.1 implies that if f ≡ , γ = 0 and γ = 0, thenthe problem (1.3) has in C q,µ (0 ,
1] only the trivial solution u ≡
0. Hence Theorem 3.1 yieldsexistence and uniqueness of a solution u of (1.3) with u ∈ C [0 ,
1] and D δ ∗ u ∈ C q, − δ (0 , Remark 3.2.
If we impose the additional hypothesis that | b ( x ) | ≥ C > on [0 , , it is thenpossible to give a simple proof of Corollary 3.1 directly from Theorem 3.1 without using Theo-rem 3.2: differentiating (1.3a) then solving for u ′′ yields u ′′ ( x ) = 1 b ( x ) (cid:20) f ′ + (cid:16) D δ ∗ u (cid:17) ′ − c ′ u − ( c + b ′ ) u ′ (cid:21) ( x ) , (3.10) and an appeal to the bounds of Theorem 3.1 yields (3.9) immediately for the case i = 2 . Onecan then differentiate (3.10) iteratively and use Theorem 3.1 to prove (3.9) for i = 3 , , . . . .If instead b ( x ) ≡ and | c ( x ) | ≥ C > on [0 , , a similar technique will work—note that b ≡ enables the conclusion of Theorem 3.1 to be strengthened to D δ ∗ u ∈ C q,ν (0 , where ν = max { µ, − δ } 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 m be a positive integer with m − < σ < m . Set r ( x ) = x σ + x σ − m +1 for x ∈ [0 , . Clearly r ∈ C m − [0 , ∩ C ∞ (0 , . Then from (Diethelm, 2010, Appendix B) one gets D σ ∗ r ( x ) = Γ( σ + 1) + Γ(2 σ − m + 2)Γ( σ − m + 2) x σ − m +1 . Hence D σ ∗ r ∈ C [0 , ∩ C ∞ (0 , and | D σ ∗ r ( x ) | ≤ Γ( σ + 1) + [Γ(2 σ − m + 2) / Γ( σ − m + 2)] , 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 ) x σ − m +1 − i for i = 1 , , . . . Thus the derivatives of D σ ∗ r satisfy the hypotheses of Theorem 3.2 and the derivatives of r agreewith (3.5) for i ≥ , i.e., the outcome of Theorem 3.2 cannot be sharpened for i ≥ . Assume the hypotheses of Corollary 3.1. Let N be a positive integer. Subdivide [0 ,
1] bythe uniform mesh x j = j/N =: jh , for j = 0 , , . . . , N . Then the standard discretization of − D δ ∗ u ( x j ) for j = 1 , , . . . , N − − D δ ∗ u ( x j ) ≈ − h δ Γ(3 − δ ) j − X k =0 d j − k ( u k +2 − u k +1 + u k ) , (4.1)where u k denotes the computed approximation to u ( x k ), and we set d r = r − δ + − ( r − − δ + for all integers r, (4.2)8ith s + = ( s if s ≥ , s < . Note that d r = 0 for r ≤ g j = g ( x j ) for each mesh point x j , where g can be b, c or f . To discretize the convectiveterm bu ′ we shall use simple upwinding (Roos et al., 2008, p.47), because the standard approxi-mation u ′ ( x j ) ≈ ( u j +1 − u j − ) / (2 h ) may yield a non-monotone difference scheme when δ is near1; see Gracia and Stynes (2013) for details. Thus we use the approximation( bu ′ )( x j ) ≈ b j D ℵ u j := ( b j ( u j +1 − u j ) /h if b j < ,b j ( u j − u j − ) /h if b j ≥ . (4.3)This difference approximation can also be written as b j D ℵ u j = − ( b j + | b j | ) u j − h + | b j | u j h + ( b j − | b j | ) u j +1 h . (4.4)The full discretization of (1.3a) is − h δ Γ(3 − δ ) j − X k =0 d j − k ( u k +2 − u k +1 + u k ) + b j D ℵ u j + c j u j = f j , (4.5)for j = 1 , , . . . , N −
1. The boundary conditions (1.3b) are discretized by approximating u ′ (0)by ( u − u ) /h and u ′ (1) by ( u N − u N − ) /h .Let A = ( a jk ) Nj,k =0 denote the ( N + 1) × ( N + 1) matrix corresponding to this discretization of(1.3), i.e., A~u = ~f where ~u := ( u u . . . u N ) T , ~f := ( γ f f . . . f N − γ ) T and the superscript T denotes transpose. Thus the 0 th row of A is ((1 + α h − ) − α h − . . .
0) and its N th row is(0 0 . . . − α h − (1 + α h − )). For j = 1 , , . . . , N −
1, the entries of the j th row of A satisfy a j = − d j h δ Γ(3 − δ ) − ǫ j b + | b | h , (4.6a) a j = − d j − + 2 d j h δ Γ(3 − δ ) − ǫ j b + | b | h + ǫ j (cid:18) | b | h + c (cid:19) , (4.6b) a jk = − d j − k + 2 d j − k +1 − d j − k +2 h δ Γ(3 − δ ) − ǫ j,k +1 b j + | b j | h + ǫ jk (cid:18) | b j | h + c j (cid:19) + ǫ j,k − b j − | b j | h for k = 2 , , . . . , N, (4.6c)where we set ǫ jk = ( j = k, . Hence the matrix A is lower Hessenberg.Observe that (4.5) implies that N X k =0 a jk = c j for j = 1 , , . . . , N − . (4.7)We shall prove various inequalities for the non-zero entries of A . First, a jj > j by(4.6), (4.2) and c ≥
0. From (4.6a), one has a j < j = 1 , , . . . , N − . (4.8)9y (4.6b), a = 2 − δ − h δ Γ(3 − δ ) − b + | b | h , (4.9)so the sign of a depends on δ, h and b . Lemma 4.1.
One has a j > for j = 3 , , . . . , N − .Proof. For j = 3 , , . . . , N −
1, equations (4.6b) and (4.2) yield h δ Γ(3 − δ ) a j = − ( j − − δ + ( j − − δ + 2 j − δ − j − − δ = 2 h ( j − − δ − ( j − − δ i (cid:20) j − δ − ( j − − δ ( j − − δ − ( j − − δ − (cid:21) . By Cauchy’s mean value theorem, for some θ ∈ ( j − , j ) we have j − δ − ( j − − δ ( j − − δ − ( j − − δ = (2 − δ ) θ − δ (2 − δ )( θ − − δ = (cid:18) − θ (cid:19) δ − > (cid:18) − (cid:19) δ − > , where we used θ > j − ≥
2. It follows that a j > Lemma 4.2.
One has a jk < for j = 1 , , . . . , N − and k ∈ { , , . . . , j − , j − , j + 1 } .Proof. If k = j + 1, the inequality a j,j +1 < k ∈ { , , . . . , j − } . Taylor expansions imply that for some η ∈ ( j − k, j − k + 2) wehave − d j − k + 2 d j − k +1 − d j − k +2 = − d dr h r − δ − ( r − − δ i (cid:12)(cid:12)(cid:12)(cid:12) r = η = − (2 − δ )(1 − δ ) h η − δ − ( η − − δ i < . Hence (4.6c) implies that a jk < k ∈ { , , . . . , j − } . This completes the description of the entries a jk defined in (4.6). The sign pattern of A is + − . . . . . . . − + − . . . . . . . − (4 .
9) + − . . . . . . . − + − + − . . . . . . . − + − − + − . . . . . . . − + − − . . . . . . . . . − + − − + − − . . . . . . . . . . . . − + − . . . . . . . . . . . . . . . . . . . . . . − + . A We shall show that A − exists and A − ≥
0. Here and subsequently, an inequality betweentwo matrices or vectors means that this inequality holds true for all the corresponding pairs ofentries in those matrices or vectors. 10he positive off-diagonal entries in column 1 of A are inconvenient for our analysis. We shallchange their signs, while simultaneously simplifying column 0 of A , by the following device. Set A ′ = E ( N − E ( N − . . . E (1) A where the elementary matrix E ( k ) := ( e ( k ) ij ) Ni,j =0 with e ( k ) ij := δ ij − a k a δ ik δ j . This multiplication of matrix A on the left by elementary matrices adds a positive multiple ofrow 0 of A to each lower row to reduce to zero the off-diagonal entries of column 0. Write A ′ = ( a ′ jk ) Nj,k =0 .Row 0 of A ′ is ( a a . . . a = 1 + α h − and a = − α h − .By construction a ′ j = 0 for j = 1 , . . . , N . For k > j we clearly have a ′ jk = a jk . Theremaining entries of column 1 of A ′ will be examined below. Lemma 4.3.
The entries of column 1 of the matrix A ′ satisfy a ′ > and a ′ j < for j = 2 , , . . . , N − . Proof.
For j = 1 , , . . . , N −
1, from (4.6) one has a ′ j = a j + α h − α h − a j = 1 h δ Γ(3 − δ ) (cid:20) − d j − + 2 d j − α h − α h − d j (cid:21) − ǫ j b + | b | h + ǫ j (cid:18) | b | h − α h − α h − · b + | b | h + c (cid:19) = 1 h δ Γ(3 − δ ) (cid:20) d j (cid:18) hh + α (cid:19) − d j − (cid:21) − ǫ j b + | b | h + ǫ j (cid:18) | b | h − α h − α h − · b + | b | h + c (cid:19) . (4.10)Hence a ′ > d = 1 , d = 0 and c ≥ d j (cid:18) hh + α (cid:19) − d j − < j = 2 , , . . . , N − . (4.11)This inequality is equivalent to d j − d j > hh + α = 1 + 1 N α + 1 . (4.12)By Cauchy’s mean value theorem, for some η ∈ ( j − , j −
1) we have d j − d j = ( j − − δ − ( j − − δ j − δ − ( j − − δ = (2 − δ ) η − δ (2 − δ )( η + 1) − δ = (cid:18) η (cid:19) δ − > (cid:18) N − (cid:19) δ − (4.13)because j ≤ N −
1. Hence, using the well-known inequality t t < ln(1 + t ) < t for t > , (cid:18) d j − d j (cid:19) > ( δ −
1) ln (cid:18) N − (cid:19) > ( δ −
1) 1 / ( N − / ( N −
2) = δ − N − (cid:18) N α + 1 (cid:19) < N α + 1 . Consequently (4.12) is proved if δ − N − ≥ N α + 1 , i.e., if N − ≤ ( δ − N α + 1) = α ( δ − N + δ −
1, which is true by (1.4). Thus (4.11) isvalid.Combining (4.10) and (4.11), one obtains immediately a ′ j < j = 2 , , . . . , N − ~ . . . T . If all the off-diagonal entries of a square matrix G are non-positive andthere exists a vector ~v ≥ ~ G~v > ~
0, then (see, e.g., Fiedler, 1986) G is an M-matrix, i.e., G − exists with G − ≥ Theorem 4.1.
The matrix A is invertible and A − ≥ .Proof. The properties of A , the construction of A ′ , and Lemma 4.3 imply that the entries of the( N + 1) × ( N + 1) matrix A ′ are positive on its main diagonal and non-positive off this diagonal.We see immediately that P Nk =0 a ′ k = a + a = 1; for j = 1 , , . . . , N − N X k =0 a ′ jk = 0 + (cid:18) a j + α h − α h − a j (cid:19) + N X k =2 a jk = (cid:18) α h − α h − − (cid:19) a j + N X k =0 a jk , but P Nk =0 a jk = c j by (4.7) and a j < N X k =0 a ′ jk = (cid:18) α h − α h − − (cid:19) a j + c j > A ′ is (0 0 . . . − α h − (1 + α h − )). Hence A ′ (1 , , . . . , T > ~ A ′ is an M-matrix. Thus ( A ′ ) − exists with ( A ′ ) − ≥ A ′ = E ( N − E ( N − . . . E (1) A, which implies that A − = ( A ′ ) − E ( N − E ( N − . . . E (1) exists, and, since it is a product of matrices with non-negative entries, A − ≥ A is said to be monotone because A − ≥
0; see, e.g., Fiedler (1986).12 .3 Error analysis
Define the truncation error ~τ := ( τ τ . . . τ N ) T by A ( u − ~u ) = ~τ , (4.14)where by “ Au ” we mean that A multiplies the restriction of u to the mesh. Then τ = ( Au ) − γ = α u ′ ( x ) − α h [ u ( x ) − u ( x )] , (4.15) τ N = ( Au ) N − γ = − α u ′ ( x N ) + α h [ u ( x N ) − u ( x N − )] , (4.16)and for j = 1 , , . . . , N − τ j = ( Au ) j − f ( x j ) = ( Au ) j + D δ ∗ u ( x j ) − b ( x j ) u ′ ( x j ) − c ( x j ) u ( x j ) . (4.17) Lemma 4.4.
There exists a constant C , which is independent of j , such that j − X k =1 (cid:2) ( j − k ) − δ − ( j − k − − δ (cid:3) k δ − ≤ Cj − δ for all integers j ≥ . (4.18) Proof.
For x ∈ R , let ⌈ x ⌉ denote the smallest integer that is greater than or equal to x . By themean value theorem applied to the function x x − δ one gets ⌈ j/ ⌉− X k =1 (cid:2) ( j − k ) − δ − ( j − k − − δ (cid:3) k δ − ≤ ⌈ j/ ⌉− X k =1 (2 − δ )( j − k − − δ k δ − ≤ (2 − δ )( j − ⌈ j/ ⌉ ) − δ ⌈ j/ ⌉− X k =1 k δ − ≤ C (2 − δ )( j/ − δ ∞ X k =1 k δ − ≤ Cj − δ , as δ < < δ < j − X k = ⌈ j/ ⌉ (cid:2) ( j − k ) − δ − ( j − k − − δ (cid:3) k δ − ≤ ⌈ j/ ⌉ δ − j − X k = ⌈ j/ ⌉ (cid:2) ( j − k ) − δ − ( j − k − − δ (cid:3) = ⌈ j/ ⌉ δ − (cid:0) j − ⌈ j/ ⌉ (cid:1) − δ ≤ ( j/ δ − ( j/ − δ = 2 j − ≤ j − δ . Adding these two bounds yields (4.18).We can now bound the truncation errors of the finite difference scheme.13 emma 4.5.
There exists a constant C such that the truncation errors in the discretization A~u = ~f of (1.3) satisfy | τ j | ≤ Ch δ − if j = 0 ,C if j = 1 ,C ( j − − δ if j = 2 , . . . , N − ,Cα h if j = N. (4.19) Proof.
We shall consider separately the four cases in (4.19).
Case j = 0 : By the mean value theorem, for some η ∈ (0 , h ) we have | τ | = α (cid:12)(cid:12)(cid:12)(cid:12) u ′ (0) − u ( h ) − u (0) h (cid:12)(cid:12)(cid:12)(cid:12) = α | u ′ (0) − u ′ ( η ) | = α (cid:12)(cid:12)(cid:12)(cid:12)Z η t =0 u ′′ ( t ) dt (cid:12)(cid:12)(cid:12)(cid:12) ≤ α Z ht =0 Ct δ − dt ≤ Cα h δ − , where we used Corollary 3.1. Case j = N : Similarly to the case j = 0, one gets | τ N | ≤ α Z t =1 − h Ct δ − dt ≤ Cα h. Case < j < N : Fix j ∈ { , , . . . , N − } . By virtue of (4.17) and (4.1) we can write τ j = R j + P j − k =0 τ j,k , where R j = b j D ℵ u j − ( bu ′ )( x j )and τ j,k = 1Γ(2 − δ ) Z x k +1 x k ( x j − s ) − δ u ′′ ( s ) ds − d j − k h δ Γ(3 − δ ) [ u ( x k +2 ) − u ( x k +1 ) + u ( x k )] . (4.20)By Taylor expansions, for some constant C one has | R j | ≤ Ch k b k ∞ max x ∈ [ x j − ,x j +1 ] | u ′′ ( x ) | . Hence Corollary 3.1, hj ≤ δ > | R j | ≤ Ch [ h ( j − δ − = C [ h ( j − δ − ( j − − ≤ C ( j − − δ . (4.21)Returning to the other component of τ j , we have τ j,k = u ′′ ( η )Γ(2 − δ ) Z x k +1 x k ( x j − s ) − δ ds − d j − k h δ Γ(3 − δ ) [ u ( x k +2 ) − u ( x k +1 ) + u ( x k )]= h − δ d j − k Γ(3 − δ ) (cid:20) u ′′ ( η ) − u ( x k +2 ) − u ( x k +1 ) + u ( x k ) h (cid:21) for some η ∈ ( x k , x k +1 ), by the mean value theorem for integrals. Taylor expansions show thatfor some η ∈ ( x k , x k +2 ) and k ≥ (cid:12)(cid:12)(cid:12)(cid:12) u ′′ ( η ) − u ( x k +2 ) − u ( x k +1 ) + u ( x k ) h (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12) u ′′ ( η ) − u ′′ ( η ) (cid:12)(cid:12) = (cid:12)(cid:12) ( η − η ) u ′′′ ( η ) (cid:12)(cid:12) ≤ Chx δ − k , η ∈ ( x k , x k +2 ), then the bound on u ′′′ givenby Corollary 3.1. Consequently | τ j,k | ≤ Ch − δ d j − k x δ − k Γ(3 − δ ) = Cd j − k k δ − Γ(3 − δ ) for k ≥ . (4.22)Now an appeal to Lemma 4.4 gives j − X k =1 | τ j,k | ≤ Cj − δ . (4.23)To complete the case 1 < j < N , it remains to bound τ j, . By (4.20) and a triangle inequalityone has | τ j, | ≤ (cid:12)(cid:12)(cid:12)(cid:12) − δ ) Z x x ( x j − s ) − δ u ′′ ( s ) ds (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) d j h δ Γ(3 − δ ) [ u ( x ) − u ( x ) + u ( x )] (cid:12)(cid:12)(cid:12)(cid:12) . (4.24)We bound these two terms separately. First, by Corollary 3.1 and j > (cid:12)(cid:12)(cid:12)(cid:12) − δ ) Z x x ( x j − s ) − δ u ′′ ( s ) ds (cid:12)(cid:12)(cid:12)(cid:12) ≤ C [( j − h ] − δ Γ(2 − δ ) Z x x s δ − ds = C ( j − − δ ( δ − − δ ) . For the second term in (4.24), the mean value theorem and Corollary 3.1 give (cid:12)(cid:12)(cid:12)(cid:12) d j h δ Γ(3 − δ ) [ u ( x ) − u ( x ) + u ( x )] (cid:12)(cid:12)(cid:12)(cid:12) = d j h | u ′ ( η ) − u ′ ( η ) | h δ Γ(3 − δ ) ≤ d j hh δ Γ(3 − δ ) Z η t = η | u ′′ ( t ) | dt ≤ Cd j hh δ Γ(3 − δ ) Z ht =0 t δ − dt = Cd j ≤ C ( j − − δ , where η ∈ ( x , x ) and η ∈ ( x , x ). Combining these inequalities with (4.24) yields | τ j, | ≤ C ( j − − δ . Add this bound to (4.21) and (4.23) to obtain finally | τ j | ≤ j − X k =0 | τ j,k | + | R j | ≤ C ( j − − δ for 1 < j < N. (4.25) Case j = 1 : This resembles the analysis above of τ j, ; one starts from (4.24) with j = 1there. The only change is that now one invokes a standard bound on Euler’s Beta function(Diethelm, 2010, Theorem D.6) to see that (cid:12)(cid:12)(cid:12)(cid:12) − δ ) Z x x ( x − s ) − δ u ′′ ( s ) ds (cid:12)(cid:12)(cid:12)(cid:12) ≤ C Γ(2 − δ ) Z x x ( x − s ) − δ s δ − ds = C Γ(2 − δ )Γ( δ − − δ ) ≤ C, (cid:12)(cid:12)(cid:12)(cid:12) d h δ Γ(3 − δ ) [ u ( x ) − u ( x ) + u ( x )] (cid:12)(cid:12)(cid:12)(cid:12) ≤ Cd = C. By the mean value theorem, for some η ∈ ( x , x ) one has | R | = (cid:12)(cid:12)(cid:12) b D ℵ u − ( bu ′ )( x ) (cid:12)(cid:12)(cid:12) ≤ k b k ∞ | u ′ ( η ) − u ′ ( x ) | ≤ C (4.26)for some constant C , because u ∈ C [0 , | τ | ≤ | τ , | + | R | ≤ C. (4.27) Remark 4.1.
In Shen and Liu (2004/05, Lemma 3) Taylor expansions are used to prove that | τ j | ≤ Ch for j = 1 , , . . . , N − under the implicit assumptions that u ′′′ and u (4) are boundedon [0 , . An inspection of the argument shows that it can be modified slightly to yield the sameresult under the assumption that u ′′′ is bounded on [0 , ; no assumption on u (4) is needed. Nev-ertheless the assumption that u ′′′ is bounded on [0 , is very strong and restricts the applicabilityof this result to special cases of (1.3) whose solutions are exceptionally smooth. We can now prove that our finite difference method is O (cid:0) h δ − (cid:1) accurate in the discretemaximum norm. Theorem 4.2.
Let b, c, f ∈ C q,µ (0 , for some integer q ≥ and µ ≤ − δ . Assume that c ≥ , α ≥ and the condition (1.4) is satisfied. Then the error in the discretization A~u = ~f of (1.3) satisfies k u − ~u k ∞ ,d ≤ Ch δ − , where k u − ~u k ∞ ,d := max j =0 , ,...,N | u ( x j ) − u j | .Proof. Recalling (4.2) and (4.6a), the construction of Section 4.2 added the multiple − a j α h − = 11 + α h − (cid:20) d j h δ Γ(3 − δ ) + ǫ j b + | b | h (cid:21) of row 0 of A to row j for j = 1 , , . . . , N −
1, yielding an M-matrix A ′ . When this constructionis applied to the system of equations (4.14), one modifies τ j to τ ′ j := τ j + τ α h − (cid:20) d j h δ Γ(3 − δ ) + ǫ j b + | b | h (cid:21) for j = 1 , , . . . , N − , (4.28)and then A ′ ( u − ~u ) = ~τ ′ , with ~τ ′ := ( τ τ ′ τ ′ . . . τ ′ N − τ N ) T . (4.29)For j = 1 , , . . . , N −
1, the proof of Theorem 4.1 shows that the j th row sum of A ′ is (cid:18) α h − α h − − (cid:19) a j + c j = − a j α h − + c j = 11 + α h − (cid:20) d j h δ Γ(3 − δ ) + ǫ j b + | b | h (cid:21) + c j . (4.30)Thus the value of the j th row sum depends strongly on j . We rescale rows 1 , , . . . , N − j th equation by(1 + α h − ) h δ Γ(3 − δ ) d j ,
16o that by (4.30) each of these row sums is now O (1) and equals1 + h δ Γ(3 − δ ) d j (cid:20) ǫ j b + | b | h + (1 + α h − ) c j (cid:21) ≥ . (4.31)Write ˜ A for the rescaled matrix of the system of equations and ~ ˜ τ for the rescaled right-handside, so now we have ˜ A ( u − ~u ) = ~ ˜ τ, with ~ ˜ τ := ( τ ˜ τ ˜ τ . . . ˜ τ N − τ N ) T (4.32)and for j = 1 , , . . . , N − τ j = (1 + α h − ) h δ Γ(3 − δ ) d j τ ′ j = (1 + α h − ) h δ Γ(3 − δ ) d j τ j + (cid:20) ǫ j b + | b | h · h δ Γ(3 − δ ) d j (cid:21) τ by (4.28). Hence, Lemma 4.5 and d j ≥ (2 − δ ) j − δ ≥ (2 − δ )2 − δ ( j − − δ for j ≥ | ˜ τ j | ≤ Ch δ − for j = 1 , , . . . , N − . (4.33)But the off-diagonal entries of ˜ A are non-positive and ˜ A (1 1 . . . T ≥ (1 1 . . . T > ~ A is an M-matrix; furthermore, it follows that in the standard matrix normnotation k · k ∞ one has k ( ˜ A ) − k ∞ ≤ | τ | ≤ Ch δ − and | τ N | ≤ Cα h fromLemma 4.5, imply that k u − ~u k ∞ ,d ≤ Ch δ − , as desired. Remark 4.2.
The order of convergence proved in Theorem 4.2 is the same order (in k · k ∞ ,d ) asthat proved in the norm k · k ∞ for the simplest collocation method ( m = 1 ) in Pedas and Tamme(2012, Theorem 4.1) on a uniform mesh, but Theorem 4.2 places no condition on the meshdiameter, while the results of Pedas and Tamme (2012) implicitly require h δ − to be smaller thansome fixed constant—this may be restrictive when δ is near 1. This mesh condition arises becausethe proofs of the main convergence results of Pedas and Tamme (2012) rely on the property that,for sufficiently large N , the operator I − P N T is invertible and k ( I − P N T ) − k is bounded in acertain operator norm; to verify this property, the authors appeal to a standard argument from(Brunner et al., 2001, Lemma 3.2), but on a uniform mesh this relies on inequality (3.12) ofBrunner et al. (2001) with r = 1 , ν = 2 − δ and m ≥ , and consequently “for sufficiently large N ” is equivalent to “for h δ − sufficiently small”. Note however that the mesh restriction is lessdemanding for the graded meshes that are also considered in Pedas and Tamme (2012). We first consider a problem whose solution u lies in C [0 , ∩ C ∞ (0 , u / ∈ C [0 ,
1] and thebehaviour of u mimics exactly the behaviour of the estimates of the solution in Corollary 3.1;cf. Example 3.1. Test Problem 1. − D δ ∗ u ( x ) + x u ′ ( x ) + (1 + x ) u ( x ) = f ( x ) on (0 , , (5.1a) u (0) − δ − u ′ (0) = γ and u (1) + u ′ (1) = γ , (5.1b)where the function f and the constants γ and γ are chosen such that the exact solution of (5.1)is u ( x ) = x δ + x δ − + 1 + 3 x − x + 4 x + x .17 E x a c t and c o m pu t ed s o l u t i on E x a c t and c o m pu t ed s o l u t i on Figure 1: Exact (solid) and computed (dashed) solutions of (5.1) for N = 64 with δ = 1 . δ = 1 . D e r i v a t i v e o f e x a c t s o l u t i on −0.02 0 0.02 0.04 0.06 0.08 0.133.23.43.63.84 x D e r i v a t i v e o f e x a c t s o l u t i on Figure 2: Derivative u ′ ( x ) of the solution of (5.1) with δ = 1 .
1, for x ∈ [0 ,
1] (left figure) and x near 0 (right figure)The numerical solution { u j } Nj =0 of problem (5.1) is computed on a uniform mesh of width h = 1 /N as described in Section 4.1. Figure 1 shows the exact solution u for δ = 1 . δ = 1 . N = 64. Figures 2 (left) and 3 (left) show the derivative u ′ of the exactsolution for δ = 1 . δ = 1 .
4, respectively. We observe that the solution u and its derivative u ′ are bounded functions. In Figures 2 (right) and 3 (right) a zoom of u ′ ( x ) in the vicinity of x = 0 is displayed, and we can observe a vertical tangent at this point in both cases.Table 5.1 presents numerical results for several values of the parameters δ and N . Each tableentry shows first the maximum pointwise error e δN := max ≤ j ≤ N | u ( x j ) − u j | , and then the order of convergence, which is computed in the standard way: p δN = log (cid:18) e δN e δ N (cid:19) . To show that the numerical results do not depend strongly on the value of δ , we have also18 D e r i v a t i v e o f e x a c t s o l u t i on −0.02 0 0.02 0.04 0.06 0.08 0.12.52.62.72.82.933.13.2 x D e r i v a t i v e o f e x a c t s o l u t i on Figure 3: Derivative u ′ ( x ) of the solution of (5.1) with δ = 1 .
4, for x ∈ [0 ,
1] (left figure) and x near 0 (right figure)computed the uniform errors for the set of values of the parameter δ considered in Table 5.1 andthe corresponding orders of convergence. These values are defined for each N by e N = max δ e δN and p N = log (cid:18) e N e N (cid:19) and they appear in the last row of Table 5.1.These numerical results show that one obtains first-order convergence for each value of δ considered in Table 5.1, and this convergence is uniform in δ . Figure 4 exhibits the pointwiseerrors in the computed solution for δ = 1 . , . N = 64 , , O ( h δ − ) convergence guar-anteed by Theorem 4.2, but our second numerical example will show that the rate of convergencecan indeed deteriorate when δ is close to 1. E rr o r s E rr o r s Figure 4: Pointwise errors in computed solutions for test problem (5.1) with δ = 1 . δ = 1 . N = 64 (solid line)and N = 128 (dashed line). 19able 5.1: Test Problem (5.1): Maximum and uniform errors e δN , e N and their orders of conver-gence p δN , p N N=64 N=128 N=256 N=512 N=1024 N=2048 δ = 1 . δ = 1 . δ = 1 . δ = 1 . δ = 1 . δ = 1 . δ = 1 . δ = 1 . δ = 1 . e N p N Test Problem 2.
Consider the constant-coefficient problem − D δ ∗ u ( x ) + 2 u ′ ( x ) + 3 u ( x ) = 1 .
25 on (0 , , (5.2a) u (0) − δ − u ′ (0) = 0 . u (1) = 1 . . (5.2b)As the exact solution of (5.2) is unknown, in order to estimate the errors of the solution { u j } Nj =0 computed on a uniform mesh of width h = 1 /N by our finite difference method, we use thetwo-mesh principle (Farrell et al., 2000, Section 5.6): on a uniform mesh of width h/
2, computethe numerical solution { z j } Nj =0 with the same method and hence the two-mesh differences d δN := max ≤ j ≤ N | u j − z j | ;from these values we estimate the orders of convergence by q δN = log (cid:18) d δN d δ N (cid:19) . The uniform two-mesh differences and their corresponding uniform orders of convergence arecomputed analogously to Table 5.1 and denoted by d N and q N , respectively. The numericalresults obtained are displayed in Table 5.2 and we observe that the finite difference method isagain convergent but a significant decrease in the order of convergence occurs for values of δ close to 1 and practical values of N . In this paper we discussed a two-point boundary value problem whose highest-order derivativewas a Caputo fractional derivative of order δ ∈ (1 , d δN , d N and theirorders of convergence q δN , q N N=64 N=128 N=256 N=512 N=1024 N=2048 δ = 1 . δ = 1 . δ = 1 . δ = 1 . δ = 1 . δ = 1 . δ = 1 . δ = 1 . δ = 1 . d N q N this differential operator on its domain [0 ,
1] provided that the boundary conditions satisfiedcertain restrictions. Then we derived sharp a priori bounds on the integer-order derivatives ofthe solution u of the boundary value problem, using elementary analytical techniques to extractthis information from previously-known bounds on the integer-order derivatives of D δ ∗ u . Thesenew bounds were used to analyse a finite difference scheme for this problem via a truncationerror analysis, but this analysis was complicated by the awkward fact that u ′′ ( x ) and u ′′′ ( x ) blowup at the boundary x = 0 of the domain [0 , O ( h δ − ) accurate at the nodes of our mesh ( h is the mesh width and the mesh isuniform), but our numerical experience has been that the method is often more accurate; wehave observed first-order convergence for all values of δ in several numerical examples, thoughone can have some deterioration in the rate of convergence when δ is near 1, as we saw in TestProblem 2.In future work in Gracia and Stynes (2013) and other papers we shall discuss the use ofalternative difference approximations of the convective term of (1.3a), investigate why the rateof convergence of (4.5) is sometimes first order for all values of δ , and extend our approach tohigher-order difference schemes and to graded meshes (cf. Pedas and Tamme (2012)). References
Mohammed Al-Refai. Basic results on nonlinear eigenvalue problems of fractional order.
Elec-tron. J. Differ. Equ. , 2012(191):1–12, 2012a.Mohammed Al-Refai. On the fractional derivatives at extreme points.
Electron. J. Qual. TheoryDiffer. Equ. , 2012(55):1–5, 2012b. ISSN 1417-3875.O. Axelsson and L. Kolotilina. Monotonicity and discretization error estimates.
SIAM J. Numer.Anal. , 27(6):1591–1611, 1990. ISSN 0036-1429. doi: 10.1137/0727093.Hermann Brunner, Arvet Pedas, and Gennadi Vainikko. Piecewise polynomial colloca-21ion methods for linear Volterra integro-differential equations with weakly singular ker-nels.
SIAM J. Numer. Anal. , 39(3):957–982 (electronic), 2001. ISSN 0036-1429. doi:10.1137/S0036142900376560. URL http://dx.doi.org/10.1137/S0036142900376560 .Kai Diethelm.
The analysis of fractional differential equations , volume 2004 of
Lecture Notes inMathematics . Springer-Verlag, Berlin, 2010. ISBN 978-3-642-14573-5. An application-orientedexposition using differential operators of Caputo type.P. A. Farrell, A. F. Hegarty, J. J. H. Miller, E. O’Riordan, and G. I. Shishkin.
Robust com-putational techniques for boundary layers , volume 16 of
Applied Mathematics (Boca Raton) .Chapman & Hall/CRC, Boca Raton, FL, 2000. ISBN 1-58488-192-5.Miroslav Fiedler.
Special matrices and their applications in numerical mathematics . MartinusNijhoff Publishers, Dordrecht, 1986. ISBN 90-247-2957-2. doi: 10.1007/978-94-009-4335-3.Translated from the Czech by Petr Pˇrikryl and Karel Segeth.Jos Luis Gracia and Martin Stynes. Upwind and central difference approximation of convectionin Caputo fractional derivative two-point boundary value problems. (In preparation) , 2013.Xia Ji and Huazhong Tang. High-order accurate Runge-Kutta (local) discontinuous Galerkinmethods for one- and two-dimensional fractional diffusion equations.
Numer. Math. TheoryMethods Appl. , 5(3):333–358, 2012. ISSN 1004-8979. doi: 10.4208/nmtma.2012.m1107.J. Tenreiro Machado, Virginia Kiryakova, and Francesco Mainardi. Recent history of fractionalcalculus.
Commun. Nonlinear Sci. Numer. Simul. , 16(3):1140–1153, 2011. ISSN 1007-5704.doi: 10.1016/j.cnsns.2010.05.027.Kassem Mustapha and William McLean. Uniform convergence for a discontinuous Galerkin,time-stepping method applied to a fractional diffusion equation.
IMA J. Numer. Anal. , 32(3):906–925, 2012. ISSN 0272-4979. doi: 10.1093/imanum/drr027.Arvet Pedas and Enn Tamme. Piecewise polynomial collocation for linear bound-ary value problems of fractional differential equations.
J. Comput. Appl. Math. ,236(13):3349–3359, 2012. ISSN 0377-0427. doi: 10.1016/j.cam.2012.03.002. URL http://dx.doi.org/10.1016/j.cam.2012.03.002 .Hans-G¨org Roos, Martin Stynes, and Lutz Tobiska.
Robust numerical methods for singularlyperturbed differential equations , volume 24 of
Springer Series in Computational Mathematics .Springer-Verlag, Berlin, second edition, 2008. ISBN 978-3-540-34466-7. Convection-diffusion-reaction and flow problems.Abbas Saadatmandi and Mehdi Dehghan. A tau approach for solution of the space fractionaldiffusion equation.
Comput. Math. Appl. , 62(3):1135–1142, 2011. ISSN 0898-1221. doi:10.1016/j.camwa.2011.04.014.Stefan G. Samko, Anatoly A. Kilbas, and Oleg I. Marichev.
Fractional integrals and derivatives .Gordon and Breach Science Publishers, Yverdon, 1993. ISBN 2-88124-864-0. Theory andapplications, Edited and with a foreword by S. M. Nikol ′ ski˘ı, Translated from the 1987 Russianoriginal, Revised by the authors.S. Shen and F. Liu. Error analysis of an explicit finite difference approximation for the spacefractional diffusion equation with insulated ends. ANZIAM J. , 46((C)):C871–C887, 2004/05.ISSN 1446-1811. 22rcilia Sousa. How to approximate the fractional derivative of order 1 < α ≤ Internat.J. Bifur. Chaos Appl. Sci. Engrg. , 22(4):1250075, 13, 2012. ISSN 0218-1274. doi: 10.1142/S0218127412500757.Erik Talvila. Necessary and sufficient conditions for differentiating under the integral sign.
Amer.Math. Monthly , 108(6):544–548, 2001. ISSN 0002-9890. doi: 10.2307/2695709.Yunying Zheng, Changpin Li, and Zhengang Zhao. A note on the finite element method for thespace-fractional advection diffusion equation.