Boundary layers in a two-point boundary value ... - Semantic Scholar
Boundary layers in a two-point boundary value problem with a Caputo fractional derivative∗ Martin Stynes† Applied Mathematics Division, Beijing Computational Science Research Center, China Department of Mathematics, National University of Ireland, Cork and Jos´e Luis Gracia‡ Institute of Mathematics and Applications, University of Zaragoza, Spain Department of Applied Mathematics, University of Zaragoza, Spain August 4, 2014
Abstract A two-point boundary value problem is considered on the interval [0, 1], where the leading term in the differential operator is a Caputo fractional derivative of order δ with 1 < δ < 2. Writing u for the solution of the problem, it is known that typically u00 (x) blows up as x → 0. A numerical example demonstrates the possibility of a further phenomenon that imposes difficulties on numerical methods: u may exhibit a boundary layer at x = 1 when δ is near 1. The conditions on the data of the problem under which this layer appears are investigated by first solving the constant-coefficient case using Laplace transforms, determining precisely when a layer is present in this special case, then using this information to enlighten our examination of the general variable-coefficient case (in particular, in the construction of a barrier function for u). This analysis proves that usually no boundary layer can occur in the solution u at x = 0, and that the quantity M = maxx∈[0,1] b(x), where b is the coefficient of the first-order term in the differential operator, is critical: when M < 1, no boundary layer is present when δ is near 1, but when M ≥ 1 then a boundary layer at x = 1 is possible. Numerical results illustrate the sharpness of most of our results.
Keywords: Fractional differential equation, Caputo fractional derivative, boundary value problem, boundary layer. MSC (2010) Classification: Primary 34B08 Abbreviated Title: A two-point boundary value problem with a Caputo derivative
1
Introduction
Boundary value problems whose differential operators involve fractional derivatives are of great interest, as these non-classical derivatives can model some physical processes where integer-order ∗
This research was partly supported by the Institute of Mathematics and Applications (IUMA), project MEC/FEDER MTM 2010-16917 and the Diputaci´ on General de Arag´ on. † Email: [email protected] ‡ Email: [email protected]
1
derivatives are unsuitable; see Jin et al. (2013); Machado et al. (2011) for an extensive list of recent applications and mathematical developments in this area. Thus the precise behaviour of solutions to fractional-derivative boundary value problems is of fundamental importance. Let δ ∈ (1, 2). Let g ∈ C 1 [0, 1] with g 0 absolutely continuous on [0, 1]. Then the Caputo fractional derivative D∗δ g associated with the point x = 0 is defined by Z x 1 δ D∗ g(x) := (x − t)1−δ g 00 (t) dt for 0 < x ≤ 1; Γ(2 − δ) t=0 see Diethelm (2010); Pedas and Tamme (2012); Stynes and Gracia (2014). The Riemannδ Liouville fractional derivative DRL of order δ associated with the point x = 0 is defined by Z x d2 1 δ DRL g(x) = 2 (x − t)1−δ g(t) dt for 0 < x ≤ 1; dx Γ(2 − δ) t=0 see Diethelm (2010). These fractional derivatives are related by the formula δ D∗δ g(x) = DRL g(x) −
g 0 (0) g(0) x−δ − x1−δ ; Γ(1 − δ) Γ(2 − δ)
(1.1)
see Diethelm (2010, Lemma 3.4). In the present paper we shall consider the two-point boundary value problem Lu(x) := −D∗δ u(x) + b(x)u0 (x) + c(x)u(x) = f (x) for x ∈ (0, 1), 0
u(0) − α0 u (0) = γ0 ,
0
u(1) + α1 u (1) = γ1 ,
(1.2a) (1.2b)
where 1 < δ < 2. The constants α0 , α1 , γ0 , γ1 and the functions b, c and f are given. We assume that b, c, f ∈ C 1 [0, 1] with c ≥ 0 in [0, 1]. We assume also that α1 ≥ 0
and α0 ≥
1 . δ−1
(1.3)
The conditions on c, α0 and α1 guarantee that (1.2) satisfies a comparison/maximum principle; see Theorem 3.1 below. The problem (1.2) models superdiffusion of particle motion when convection is present; see the discussion and references in Jin et al. (2013, Section 1). It is a member of the general class of boundary value problems that is analysed in Pedas and Tamme (2012). It is also discussed in Al-Refai (2012). Numerical methods for its solution are presented in Gracia and Stynes (2015); Kopteva and Stynes (2014); Stynes and Gracia (2014) and their references. Existence and uniqueness of a classical solution to (1.2) is established in Stynes and Gracia (2014). It is proved that u ∈ C 1 [0, 1] ∩ C 2 (0, 1], and for some constants C˜i one has the sharp bounds ( C˜i if i = 0, 1, (i) |u (x)| ≤ (1.4) δ−i ˜ C2 x if i = 2, for all x ∈ (0, 1]. Thus u00 (x) may blow up at the interval endpoint x = 0.
These results tell us a lot about the nature of the solution u, but in one respect they are seriously deficient: all constants C˜i that appear above depend on the parameter δ, but they can be extremely large when δ is near 1, because in certain cases (as we shall see) the solution u develops a boundary layer at x = 1 (i.e., |u0 (1)| becomes very large) when δ is near 1. It is well known that, when computing numerical solution of problems with integer-order differential 2
operators, such layers can cause a deterioration in accuracy (Roos et al., 2008). This is also the case for the fractional-derivative problem (1.2): see Gracia and Stynes (2015); Jin et al. (2013); Stynes and Gracia (2014), where computed solutions of (1.2) become less accurate when δ is near 1. This loss of accuracy appears in only some numerical examples in these papers, and no explanation is given there, but it is in fact confined to problems whose solutions exhibit a boundary layer at x = 1. In the present paper we shall cast light on when such layers appear in solutions of (1.2). For a concrete example exhibiting a boundary layer, consider (1.2) with b ≡ 1.9, c ≡ 0 and f ≡ 1. Take α0 = 1/(δ − 1), α1 = 0, γ0 = 0.4 and γ1 = 1.7. For constant-coefficient problems like this, an explicit formula for u(x) will be derived in Section 2 below. Using this formula, we plot the solution for the values δ = 1.6, 1.4, 1.2, 1.1 in Figure 1. It is clear from this figure that a boundary layer at x = 1 develops in the solution when δ is near 1. 2.05 2
1.7
1.95
Exact solution u
Exact solution u
1.6
1.5
1.4
1.9 1.85 1.8 1.75 1.7
1.3
1.65 1.2
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
x
0.6
0.8
1
0.6
0.8
1
x
3.6 3.4
6
3.2
5.5 5 Exact solution u
Exact solution u
3 2.8 2.6 2.4 2.2
4.5 4 3.5 3
2
2.5
1.8
2
1.6 0
0.2
0.4
0.6
0.8
1
0
x
0.2
0.4 x
Figure 1: Exact solution for b ≡ 1.9, c ≡ 0, f ≡ 1, α0 = 1/(δ − 1), α1 = 0, γ0 = 0.4 and γ1 = 1.7, with δ = 1.6 (1st row, left), δ = 1.4 (1st row, right), δ = 1.2 (2nd row, left) and δ = 1.1 (2nd row, right), showing development of a boundary layer when δ is near 1 Figure 1 of Roop (2008) also shows a boundary layer developing as the order of the fractional derivative approaches a certain limiting value, though the boundary value problem under discussion there is not the same as (1.2). Furthermore, there is currently great interest in the construction of numerical methods for problems with variable order of fractional derivative (see, e.g., Chen et al. (2013); Samko (2013)); this implies that one should pay close attention to how the solution of (1.2) depends on the value of the parameter δ. For these two reasons (loss of accuracy in numerical solution only for certain values of δ; 3
design of numerical methods for variable δ) our aim in the present paper is to investigate in detail how u behaves as a function of δ (as well as a function of x). The structure and main results of the paper are as follows. In Section 2 Laplace transforms are used to derive an explicit formula for the solution u of (1.2) in the special case when b is a nonzero constant and c ≡ 0. From this formula we deduce that, when δ is near 1, a boundary layer in u at x = 1 can appear only if b ≥ 1, while u never has a boundary layer at x = 0 even though max[0,1] |u(x)| blows up as δ → 1+ if b ≥ 1. In Section 3 the general case of (1.2) is considered and a comparison/maximum principle is used, with some guidance from Section 2, to explore how u depends on δ. In this general case we find that a boundary layer in u at x = 1 when δ is near 1 is possible only when max[0,1] b(x) ≥ 1. Note here that it is the maximum of b, not of |b|, that is the significant quantity. It is shown also that a boundary layer at x = 0 (when δ is near 1) is possible only if max[0,1] b(x) > 1 and min[0,1] b(x) ≤ 0 and min[0,1] c(x) = 0.
Notation. We use the standard notation C k (I) to denote the space of real-valued functions whose derivatives up to order k are continuous on an interval I, and write C(I) for C 0 (I). For each g ∈ C[0, 1], set kgk∞ = maxx∈[0,1] |g(x)|. In several inequalities C denotes a generic positive constant that depends on the data b, c, f, γ0 , γ1 , α1 of the boundary value problem (1.2) but is independent of δ and x; note that C can take different values in different places. A subscripted C (e.g., C1 ) denotes a fixed positive constant that can depend on all the data of the boundary value problem (1.2) except δ and x.
The case where b is constant and c ≡ 0
2
In Section 2 we consider the special case of problem (1.2) where c ≡ 0, and b is constant with b 6= 0. Our results could be extended to the case where b and c are arbitrary constants, but when c 6= 0 the details become much more complex; see Remark 2.1. We shall use Laplace transforms to derive an explicit formula for the solution u of (1.2) in terms of Mittag-Leffler functions. Our examination of this special case gives useful and penetrating insights into the properties of the solution u of (1.2). Furthermore, the precise form of the solution that we find in Section 2.2 is very helpful when constructing a barrier function in Section 3 to analyse the solution of (1.2) when b and c are no longer constants.
2.1
General right-hand side f
Extend the domain of f to [0, ∞) in such a way that the extension (which we also call f ) is smooth and has support in [0, 2]. To solve (1.2), we treat it as an initial-value problem on [0, ∞) with the initial condition u(0) − α0 u0 (0) = Rγ0 from (1.2b), and apply the standard Laplace ∞ transform operator L, defined by Lv(s) = t=0 e−st v(t) dt. The second boundary condition u(1) + α1 u0 (1) = γ1 of (1.2b) will be invoked later. Define the two-parameter Mittag-Leffler function by Eα,β (z) =
∞ X k=0
zk Γ(αk + β)
for α, β, z ∈ R with α > 0.
(2.1)
We shall need the properties (Podlubny, 1999, (1.80),(1.82)) that for constant α, β, γ and λ one has n o sα−β L xβ−1 Eα,β (±λxα ) = α (2.2) s ∓λ 4
and
γ DRL xβ−1 Eα,β (λxα ) = xβ−γ−1 Eα,β−γ (λxα ).
(2.3)
γ Note that when γ = 1 then DRL ≡ d/dx (Diethelm, 2010, p.27); thus (2.3) implies that
d β−1 x Eα,β (λxα ) = xβ−2 Eα,β−1 (λxα ). (2.4) dx Applying L to (1.2a) and observing that L{D∗δ u} = sδ−2 s2 L{u} − su(0) − u0 (0) from Podlubny (1999, (2.253)), we obtain (−sδ + bs)L{u} + sδ−1 u(0) + sδ−2 u0 (0) − bu(0) = L{f }. Hence L{u} = −
L{f } u(0) sδ−3 u0 (0) L{f } γ0 + α0 u0 (0) sδ−3 u0 (0) + + δ−1 + = − + s s s(sδ−1 − b) sδ−1 − b s(sδ−1 − b) s −b
(2.5)
using the boundary condition (1.2b) at x = 0. To find the inverse Laplace transform of −L{f }/[s(sδ−1 − b)] we imitate Diethelm (2010, p.135). Using the integration theorem for Laplace transforms, R x L r=0 f (r)dr L{f } = . (2.6) s(sδ−1 − b) sδ−1 − b By (2.4) one has d Eδ−1,1 (bxδ−1 ) = x−1 Eδ−1,0 (bxδ−1 ), (2.7) dx where the first term in the Mittag-Leffler series for Eδ−1,0 (bxδ−1 ) vanishes since Γ(0) = ∞. But (2.2) yields sδ−2 L{Eδ−1,1 (bxδ−1 )} = δ−1 s −b so the differentiation theorem for Laplace transforms gives n o n o L x−1 Eδ−1,0 (bxδ−1 ) = sL Eδ−1,1 (bxδ−1 ) − Eδ−1,1 (0) =
b sδ−1 − 1 = δ−1 . δ−1 s −b s −b
Consequently from (2.6) we have o L{f } 1 n −1 δ−1 = L x E (bx ) L δ−1,0 b s(sδ−1 − b)
Z
x
f (r)dr .
r=0
Now the convolution theorem for Laplace transforms yields Z x−t Z L{f } 1 x −1 −1 δ−1 L (x) = t Eδ−1,0 (bt ) f (r) dr dt. b t=0 s(sδ−1 − b) r=0
(2.8)
Taking the inverse transform of (2.5) and invoking (2.8) and (2.2), we get Z x−t Z 1 x −1 t Eδ−1,0 (btδ−1 ) f (r) dr dt. (2.9) u(x) = γ0 + α0 u0 (0) + u0 (0)xEδ−1,2 (bxδ−1 ) − b t=0 r=0 By virtue of (2.4) one can differentiate (2.9) to obtain Z 1 x −1 0 0 δ−1 u (x) = u (0)Eδ−1,1 (bx ) − t Eδ−1,0 (btδ−1 )f (x − t) dt. b t=0 5
(2.10)
Consequently, imposing the boundary condition u(1) + α1 u0 (1) = γ1 of (1.2b), one has Z 1−t Z 1 1 −1 0 0 δ−1 γ1 =γ0 + α0 u (0) + u (0)Eδ−1,2 (b) − t Eδ−1,0 (bt ) f (r) dr dt b t=0 r=0 Z 1 1 −1 + α1 u0 (0)Eδ−1,1 (b) − t Eδ−1,0 (btδ−1 )f (1 − t) dt b t=0 whence u0 (0) =
γ1 − γ0 +
1 b
R1
−1 δ−1 ) t=0 t Eδ−1,0 (bt
hR 1−t r=0
i f (r) dr + α1 f (1 − t) dt
α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b)
.
(2.11)
Substitution of (2.11) into (2.9) and (2.10) yields explicit formulas for u(x) and u0 (x).
2.2
Constant right-hand side f
In this section we simplify the results of Section 2.1 by taking f to be constant so that we can then investigate in detail the solution u. Furthermore, the formulas of Section 2.2 will be of great help in the construction of a barrier function in Section 3 to analyse the structure and behaviour of u in the case of variable b, c and f . First, taking f constant in (2.11) and recalling (2.7), we have R1 d Eδ−1,1 (btδ−1 ) dt γ1 − γ0 + fb t=0 [α1 + (1 − t)] dt 0 u (0) = α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b) i h R1 γ1 − γ0 + fb α1 Eδ−1,1 (b) − α1 − 1 + t=0 Eδ−1,1 (btδ−1 ) dt = α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b) =
γ1 − γ0 + fb [α1 Eδ−1,1 (b) − α1 − 1 + Eδ−1,2 (b)] α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b)
(2.12)
where we integrated by parts then invoked (2.4) with β = 2. Substituting (2.12) into (2.9) gives us u(x); this formula can be written in a variety of ways, the most compact of which seems to be α0 + xEδ−1,2 (bxδ−1 ) [γ1 − γ0 − (1 + α0 + α1 )f /b] (α0 + x)f u(x) = γ0 + + . (2.13) b α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b) The correctness of this formula can be verified by substitution into (1.2), using (2.4) and, by (1.1) and (2.3), D∗δ xEδ−1,2 (bxδ−1 ) = x1−δ Eδ−1,2−δ (bxδ−1 ) − = x1−δ
∞ X k=1
x1−δ Γ(2 − δ)
(bxδ−1 )k Γ(k(δ − 1) + 2 − δ)
= bEδ−1,1 (bxδ−1 ). Remark 2.1. In the case of constant b, c, f with c 6= 0, one can use the above Laplace transform technique to obtain a closed-form representation of the solution by imitating (Mathai et al., 2006, (35)). For example, to invert [s(sδ − bs − c)]−1 , write r X ∞ ∞ X 1 bs br sr−1 1 1 = = = δ δ δ δ s(s − bs − c) s(s − c) s −c (s − c)r+1 s(sδ − c) 1 − sδbs−c r=0 r=0 6
then appeal to (Podlubny, 1999, (1.80)): n o (k) L xαk+β−1 Eα,β (±λxα ) =
k! sα−β (sα ∓ λ)k+1
(2.14)
(k)
where Eα,β (y) = (dk /dy k )Eα,β (y). This will yield an infinite sum of derivatives of Mittag-Leffler functions, which is very complicated. As our main aim in Section 2 is to gain insight into the solution of the general problem (1.2), we do not consider c 6= 0 any further. Remark 2.2. Differentiating (2.13) twice by invoking (2.4), one gets γ1 − γ0 − (1 + α0 + α1 )f /b 1 00 u (x) = Eδ−1,0 (bxδ−1 ). α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b) x
(2.15)
When γ1 − γ0 6= (1 + α0 + α1 )f /b, on recalling the definition (2.1) and Γ(0) = ∞ it follows ˜ δ−2 . This coincides with the bound on |u00 (x)| that is that as x → 0 one has u00 (x) ∼ Cx derived in Stynes and Gracia (2014, Corollary 3.5) for the general case of variable b, c and f . Furthermore, the next term in the series expansion of the Mittag-Leffler function in (2.15) is ˜ −1 (xδ−1 )2 = Cx ˜ 2δ−3 which matches Stynes and Gracia (2014, Example 3.7), where it was Cx shown that the derivative bounds (1.4) are sharp. Remark 2.3. The coefficient of f in the formula (2.13) for u(x) is (α0 + x)[α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b)] − [α0 + xEδ−1,2 (bxδ−1 )](1 + α0 + α1 ) b[α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b)] n 1 = α0 [Eδ−1,2 (b) − 1 − xEδ−1,2 (bxδ−1 ) + x] b[α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b)] o +xEδ−1,2 (b) − xEδ−1,2 (bxδ−1 ) + α1 [(α0 + x)Eδ−1,1 (b) − α0 − xEδ−1,2 (bxδ−1 )] =
n 1 α0 [Eδ−1,δ+1 (b) − xδ Eδ−1,δ+1 (bxδ−1 )] 1 + α0 + bEδ−1,δ+1 (b) + α1 [1 + bEδ−1,δ (b)] o +xEδ−1,δ+1 (b) − xδ Eδ−1,δ+1 (bxδ−1 ) + α1 [(α0 + x)Eδ−1,δ (b) − xδ Eδ−1,δ+1 (bxδ−1 )] ,
(2.16)
where we used the elementary identity Eδ−1,1+i (z) = zEδ−1,δ+i (z) + 1
for i = 0, 1.
(2.17)
The formula (2.16) is the inspiration for our choice of barrier function in Theorem 3.5. We turn now to our main interest: investigating how the solution u changes as δ approaches 1, and in particular determining when boundary layers appear in u. Observe first that, since α0 > 1 by condition (1.3), it follows from (2.12) that |u0 (0)| ≤ C for some constant C, i.e., there is never a boundary layer in u at x = 0. In the subsections that follow, we show that when δ is near 1, the magnitudes of kuk∞ and depend strongly on whether b < 1 or b ≥ 1.
|u0 (1)|
Differentiating (2.13) gives
u0 (1) =
(γ1 − γ0 )Eδ−1,1 (b) f α0 [1 − Eδ−1,1 (b)] + Eδ−1,2 (b) − Eδ−1,1 (b) · + . (2.18) b α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b) α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b) 7
2.2.1
Case b > 1
Assume in Section 2.2.1 that b > 1 and f 6= 0. Suppose that 1/(δ − 1) = K for some positive integer K. Then for i = 1, 2 we have Eδ−1,i (b) =
∞ X k=0
bk bK 1 ≥ ≥ b1/(δ−1) . Γ((δ − 1)k + i) Γ(1 + i) 2
(2.19)
It follows that limδ→1+ Eδ−1,i (b) = ∞ and Eδ−1,i (b) 1/(δ − 1) as δ → 1+ . Consequently, when α0 = 1/(δ − 1), by virtue of (2.12) one has limδ→1+ u0 (0) = f /b and hence, recalling (1.2b), ( ∞ if f > 0, lim u(0) = lim [γ0 + α0 u0 (0)] = + + δ→1 δ→1 −∞ if f < 0. Thus when b > 1 and f 6= 0 we do not have kuk∞ ≤ C independently of δ. Nevertheless, recalling from above that limδ→1+ u0 (0) = f /b, there is no boundary layer at x = 0. To discuss limδ→1+ u0 (1), we use more sophisticated machinery. Set σ = δ − 1 for brevity. Let γ(1, ϕ) be the complex-plane contour of Figure 2 where we choose ϕ = 2σπ/3. This contour divides the complex plane into two regions, which we denote by G− (1, ϕ) and G+ (1, ϕ). The real number b > 1 lies in the region G+ (1, ϕ), so by (Podlubny, 1999, Theorem 1.1), for arbitrary but fixed n > 0 we have Z 1 (1−n)/σ 1 ζ (1−n)/σ exp(ζ 1/σ ) 1/σ Eσ,n (b) = b exp(b ) + dζ, (2.20) σ 2πσı γ(1,ϕ) ζ −b where ı is the imaginary unit.
y
G− (1, ϕ) G+ (1, ϕ)
ϕ 1 b
O γ(1, ϕ)
x
−ϕ
Figure 2: Contour γ(1, ϕ) of equation (2.20) Our interest lies in what happens when δ → 1+ (i.e., σ → 0+ , so ϕ → 0+ ) with b > 1 fixed. For ζ ∈ γ(1, ϕ) and ϕ sufficiently small we have |ζ − b| ≥ min{b − 1, b sin ϕ} = b sin ϕ ≥ bϕ/2 = bσπ/3. 8
Consequently (cf. (Podlubny, 1999, p.33)), for some positive constants C one obtains Z 1 Z (1−n)/σ exp(ζ 1/σ ) ζ C C |ζ (1−n)/σ | exp |ζ|1/σ cos(2π/3) dζ ≤ 2 dζ ≤ 2 2πσı γ(1,ϕ) σ γ(1,ϕ) ζ −b σ as σ → 0+ , where the contour integral is bounded since cos(2π/3) < 0. Thus (2.20) implies that for b > 1 one has 1 (1−n)/(δ−1) 1 1/(δ−1) Eδ−1,n (b) = as δ → 1+ . (2.21) b exp b +O δ−1 (δ − 1)2 It follows from (2.21) that Eδ−1,1 (b)/Eδ−1,2 (b) ≈ b1/(δ−1) when δ is near 1. Consequently, when δ is near 1, if α1 > 0 then (2.18) and (2.19) imply that |u0 (1)| ≈ C min{α0 , Eδ−1,1 (b)},
(2.22)
while if α1 = 0, then (2.18) and (2.19) imply that |u0 (1)| ≈
Cα0 b1/(δ−1) Eδ−1,2 (b) Cα0 Eδ−1,1 (b) ≈ ≈ Cb1/(δ−1) min{α0 , Eδ−1,2 (b)}. α0 + Eδ−1,2 (b) α0 + Eδ−1,2 (b)
(2.23)
Thus when b > 1 there is a boundary layer at x = 1, and it is much stronger when α1 = 0. 2.2.2
Case b = 1
Assume in Section 2.2.2 that b = 1 and f 6= 0. Then (2.13) yields α0 [γ1 − γ0 − (1 + α0 + α1 )f ] u(0) = γ0 + α0 f + α0 + Eδ−1,2 (1) + α1 Eδ−1,1 (1) γ1 − γ0 + f [Eδ−1,2 (1) + α1 Eδ−1,1 (1) − 1 − α1 ] . = γ0 + α0 α0 + Eδ−1,2 (1) + α1 Eδ−1,1 (1)
(2.24)
Suppose that 1/(δ − 1) is an integer. Then Eδ−1,2 (1) =
∞ X k=0
1
2
−1
3
−1
−1
δ−1 δ−1 δ−1 X 1 X 1 X 1 1 ≥ + + + ··· Γ((δ − 1)k + 2) Γ(3) Γ(4) Γ(5) 1 2
k=0
k= δ−1
1 1 1 1 + + + ··· δ − 1 Γ(3) Γ(4) Γ(5) e−2 = . δ−1
k= δ−1
=
(2.25)
Likewise, one has 1 −1 δ−1
Eδ−1,2 (1) ≤
X k=0
1 + Γ(2)
2 −1 δ−1
X 1 k= δ−1
3 −1 δ−1
X 1 1 + + ··· Γ(3) Γ(4) 2 k= δ−1
1 1 1 1 = + + + ··· δ − 1 Γ(2) Γ(3) Γ(4) e−1 = . δ−1 9
(2.26)
One can show similarly that
e−1 e ≤ Eδ−1,1 (1) ≤ . δ−1 δ−1
(2.27)
Invoking (2.25)–(2.27) in (2.24), for f > 0 and δ sufficiently close to 1 one obtains h i 1 (e−1) γ1 − γ0 + f e−2+α − 1 − α 1 δ−1 u(0) ≥ γ0 + α0 1e α0 + e−1+α δ−1 h i 1 (e−1) γ1 − γ0 + f e−2+α − 1 − α1 δ−1 = γ0 + . 1e 1 + e−1+α α0 (δ−1) But α0 ≥ 1/(δ − 1), so it follows that as δ → 1+ with 1/(δ − 1) an integer, then u(0) → ∞ if f > 0. Otherwise, if f < 0 then u(0) → −∞ as δ → 1+ . Thus when b = 1 and f 6= 0 we do not have kuk∞ ≤ C independently of δ.
Furthermore, the condition (1.3) compared with (2.26) and (2.27) shows that α0 ≥ 13 Eδ−1,i (1) for i = 1, 2 when δ is near 1. Consequently (2.18) implies that |u0 (1)| ≈ CEδ−1,1 (1) ≈ C/(δ − 1) when δ is near 1; thus u exhibits a boundary layer at x = 1. 2.2.3
Case 0 ≤ b < 1
Assume in Section 2.2.3 that 0 ≤ b < 1. Then Eδ−1,2 (b) =
∞ X k=0
∞
X bk 1 bk = ≤ . Γ((δ − 1)k + 2) 1−b
(2.28)
k=0
The estimate (2.28) is qualitatively sharp when δ is near 1 because for 0 ≤ k ≤ b1/(δ − 1)c, where bnc denotes the greatest integer less or equal to n, one has 1 1 1 ≥ = Γ((δ − 1)k + 2) Γ(3) 2 so Eδ−1,2 (b) ≥
b1/(δ−1)c k X b k=0
2
=
1 − b1+b1/(δ−1)c 1 ≥ for δ sufficiently close to 1. 2(1 − b) 4(1 − b)
(2.29)
Similarly to (2.28) and (2.29), one has 1 1 ≤ Eδ−1,1 (b) ≤ , 2(1 − b) θ(1 − b)
(2.30)
where θ = min{Γ(x) : 1 ≤ x ≤ 2}.
From (2.28)–(2.30) and (2.12), since α0 ≥ 1/(δ − 1) one obtains |u0 (0)| ≤ C/α0 for some constant C, so u0 (0) → 0 as δ → 1+ . Furthermore, (2.13) and (2.18) now yield kuk∞ + |u0 (1)| ≤ C for 0 ≤ b < 1 for some constant C = C(b). Thus no boundary layers are present in u when 0 ≤ b < 1.
10
(2.31)
Remark 2.4. The analysis above of the cases b > 1, b = 1 and 0 ≤ b < 1 shows that, as a function of δ, the nature of the solution u undergoes a fundamental change when b moves from b ≥ 1 to b < 1. This is shown graphically in Figure 3, where u changes dramatically as b moves from 0.9 to 1.1 (as well as the evident changes in shape in these graphs, note the changes in scale of the y axis). The regime b < 0 is not discussed separately in Section 2 since we are able to handle it easily in Section 3 when we discuss the general variable-coefficient case; see Theorem 3.3.
12 11
6
10
5.5
9 Exact solution u
Exact solution u
7 6.5
5 4.5 4 3.5
8 7 6 5
3
4
2.5
3
2
2 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
x
0.6
0.8
1
x
18 16
Exact solution u
14 12 10 8 6 4 2 0
0.2
0.4
0.6
0.8
1
x
Figure 3: Exact solution for δ = 1.05, c ≡ 0, f ≡ 1, α0 = 1/(δ − 1), α1 = 0, γ0 = 0.4, γ1 = 1.7, and b = 0.9 (top left figure), b = 1.0 (top right figure), b = 1.1 (bottom figure), showing effect on solution u of moving from b < 1 to b > 1
3
Boundary layers in the variable-coefficient problem
In numerical solutions of (1.2) computed by the finite difference method of Stynes and Gracia (2014), when δ is near 1 we have observed boundary layers at x = 1 in certain examples but we have never observed a layer at x = 0. The results of Section 3 will extend those of Section 2, prove in most cases that a boundary layer cannot occur at x = 0, and provide substantial information about when boundary layers can appear at x = 1 when δ is near 1. The key tool in the analysis presented in Section 3 is the following comparison principle, which will be used several times. Theorem 3.1. (Stynes and Gracia, 2014, Theorem 2.1) Let z ∈ C 1 [0, 1] ∩ C 2 (0, 1]. Assume that |z 00 (x)| ≤ Kx−θ for 0 < x ≤ 1, where θ ∈ (0, 1) and K are constants that are independent 11
of x. Let b, c ∈ C[0, 1] with c(x) ≥ 0 for all x ∈ (0, 1). Assume that z satisfies the inequalities − D∗δ z + bz 0 + cz ≥ 0 on (0, 1), 0
(3.1a)
0
z(0) − α0 z (0) ≥ 0 and z(1) + α1 z (1) ≥ 0,
(3.1b)
where α0 and α1 satisfy (1.3). Then z ≥ 0 on [0, 1]. The bounds (1.4) show that the solution u of (1.2) satisfies the regularity hypotheses imposed on z in Theorem 3.1. We shall apply Theorem 3.1 to Bi ± u for various barrier functions Bi , concluding that Bi ± u ≥ 0, i.e., that |u| ≤ Bi on [0, 1]. We begin with a general result that provides a useful relationship between u0 (1) and u0 (0). For any function φ ∈ C[0, 1], set kφk[0,x] = max0≤t≤x |φ(t)| for 0 ≤ x ≤ 1. Lemma 3.1. There exists a constant C such that |u0 (x)| ≤ |u0 (0)| + Cxδ−1 1 + kuk[0,x] + |u0 (0)| Eδ−1,1 (kbk[0,x] xδ−1 )
for 0 ≤ x ≤ 1.
Proof. Set y(x) = u0 (x) − u0 (0) for 0 ≤ x ≤ 1. On multiplying (1.2a) by (x − t)δ−2 /Γ(δ − 1) then integrating from t = 0 to t = x, after some manipulation of the fractional derivative term one obtains (Kopteva and Stynes, 2014) a weakly singular Volterra integral equation of the second kind in the unknown y: for 0 < x ≤ 1, Z x Z x 1 1 y(x) − (x − t)δ−2 b(t)y(t) dt = (x − t)δ−2 {b(t)u0 (0) + c(t)u(t) − f (t)} dt Γ(δ − 1) t=0 Γ(δ − 1) t=0 := g(x), say. (3.2) From Brunner (2004, p.343) the solution of (3.2) can be expressed as Z x y(x) = g(x) + Rα (x, t)g(t) dt for 0 ≤ x ≤ 1,
(3.3)
t=0
where δ−2
Rα (x, t) = (x − t)
∞ X
(x − t)(n−1)(δ−1) Φn (x, t; δ)
(3.4)
n=1
and the Φn are defined iteratively by Φ1 (x, t; δ) := b(t)/Γ(δ − 1) Z 1 1 (1 − z)δ−2 z (n−1)(δ−1)−1 b(t + (x − t)z) Φn−1 (t + (x − t)z, t; δ) dz Φn (x, t; δ) := Γ(δ − 1) z=0 for n = 1, 2, . . . An inductive argument using Euler’s beta function shows easily (cf. Brunner (2004, Lemma 6.1.3)) that |Φn (x, t; δ)| ≤
kbkn[0,x]
Γ(n(δ − 1))
for n = 1, 2, · · · and 0 ≤ t ≤ x ≤ 1.
(3.5)
Consequently the infinite series in (3.4) is uniformly convergent for 0 ≤ t ≤ x ≤ 1. Hence we can move the summation sign of (3.4) outside the integral in (3.3). Again appealing to (3.5),
12
from (3.3) we get |y(x)| ≤ |g(x)| + kgk[0,x] = |g(x)| + kgk[0,x] = |g(x)| + kgk[0,x]
∞ X n=1 ∞ X n=1 ∞ X n=1
kbkn[0,x]
Z
Γ(n(δ − 1))
x
t=0
(x − t)(n−1)(δ−1)+δ−2 dt
kbkn[0,x] xn(δ−1)
n(δ − 1)Γ(n(δ − 1)) (kbk[0,x] xδ−1 )n Γ(n(δ − 1) + 1)
≤ kgk[0,x] Eδ−1,1 (kbk[0,x] xδ−1 )
(3.6)
by the definition (2.1) of the Mittag-Leffler function. Recalling the definition of g in (3.2), it follows that kgk[0,x]
1 Γ(δ − 1) xδ−1 = C 1 + kuk[0,x] + |u0 (0)| Γ(δ) δ−1 0 ≤ C 1 + kuk[0,x] + |u (0)| x 0
≤ C 1 + kuk[0,x] + |u (0)|
Z
x
t=0
(x − t)δ−2 dt
(3.7)
as inf δ∈(1,2) Γ(δ) > 0. Combining (3.6), (3.7) and u0 (x) = u0 (0) + y(x) yields the required result. Inequality (3.6) is sharp if g and b are positive constants; see Brunner (2004, Theorem 6.1.1). First, we show that in the special case where c > 0 on [0, 1], the solution u is uniformly bounded for δ ∈ (1, 2) and no boundary layer appears at x = 0 when δ is near 1; there is also no boundary layer at x = 1 if α1 > 0. If c ≥ c > 0 for some constant c, set kf k∞ C1 = max |γ0 |, |γ1 |, c
.
Theorem 3.2. Assume that c ≥ c > 0 for some constant c. Then (i)
kuk∞ ≤ C1 ,
(3.8)
C1 + |γ0 | ≤ (δ − 1)(C1 + |γ0 |), α0 C1 + |γ1 | (iii) |u0 (1)| ≤ if α1 > 0. α1 (iv) |u0 (1)| ≤ CEδ−1,1 (kbk∞ ) for some constant C if α1 = 0. (ii) |u0 (0)| ≤
(3.9) (3.10) (3.11)
Proof. (i) Define the constant barrier function B1 (x) ≡ C1 . Then B1 ± u ≥ 0 by Theorem 3.1, so (3.8) is valid. (ii) By (1.2b) we have u(0) − α0 u0 (0) = γ0 , so (1.3) and (3.8) yield u(0) − γ0 kuk∞ + |γ0 | C1 + |γ0 | 0 ≤ |u (0)| = ≤ ≤ (δ − 1)(C1 + |γ0 |). α0 α0 α0 (iii) Use the boundary condition u(1) + α1 u0 (1) = γ1 and (3.8) to derive (3.10). (iv) This follows immediately from Lemma 3.1, (3.8) and (3.9). 13
Theorem 3.2 shows that when c is strictly positive on [0, 1] there is no boundary layer at x = 0 and, if one also has α1 > 0, there is no boundary layer at x = 1. But when α1 = 0 one may have a layer at x = 1, as the next example demonstrates. Example 3.1. Consider our boundary value problem (1.2) with δ = 1.01, α0 = 1/(δ − 1), α1 = 0, γ0 = 0.4, γ1 = 1.7 and constant functions b(x) ≡ 1.1, c(x) ≡ 1 and f (x) ≡ 3.25. In Figure 4 we show the solution (values at the mesh points, joined by a piecewise linear curve) computed by the finite difference scheme of Stynes and Gracia (2014) on a uniform mesh of width 1/2048. A 3.2 3
Computed solution
2.8 2.6 2.4 2.2 2 1.8 0
0.1
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
0.9
1
Figure 4: Example showing boundary layer at x = 1 is possible when c > 0 and α1 = 0 boundary layer at x = 1 is clearly visible. Set M = max b(x). x∈[0,1]
We observe that M can have any sign. In Theorems 3.3, 3.4 and 3.5 we shall consider the regimes M < 0, 0 ≤ M ≤ 1 and M > 1 respectively and in each case we shall derive bounds on kuk∞ , |u0 (0)| and |u0 (1)|. The case where b is strictly negative is addressed first. For M < 0, set (1 + α1 )kf k∞ kf k∞ + C2 kck∞ C2 = max |γ0 |, |γ1 | − and C3 = max |γ0 | + C2 , . M −M
Theorem 3.3. Assume that M < 0. Then (i) (ii) (iii)
kuk∞ ≤ C2 , 0
|u (0)| ≤ (δ − 1)(C2 + |γ0 |), 0
|u (1)| ≤ C3 .
(3.12) (3.13) (3.14)
Proof. (i) Set B2 (x) = C2 + xkf k∞ /M for 0 ≤ x ≤ 1. Then B2 ≥ 0 so LB2 (x) = b(x)B20 (x) + c(x)B2 (x) ≥ B2 (0) − α0 B20 (0) = C2 −
b(x)kf k∞ ≥ kf k∞ , M
α0 kf k∞ (1 + α1 )kf k∞ ≥ |γ0 | and B2 (1) + α1 B20 (1) = C2 + ≥ |γ1 | M M
by definition of C2 . It follows from (1.2) and Theorem 3.1 that B2 ± u ≥ 0, i.e., |u(x)| ≤ B2 (x) for 0 ≤ x ≤ 1. As B2 (x) ≤ C2 , we have proved (3.12). 14
(ii) To obtain (3.13), use the boundary condition u(0) − α0 u0 (0) = γ0 , (1.3) and (3.12). (iii) Set w(x) = u(x) − u(1). Then by (1.2) and (3.12) we have
|Lw(x)| = |f (x) − c(x)u(1)| ≤ kf k∞ + C2 kck∞ ,
|w(0) − α0 w0 (0)| = |γ0 − u(1)| ≤ |γ0 | + C2 ,
w(1) = 0.
Set v(x) = C3 (1 − x). Then Lv(x) = −C3 b(x) + C3 (1 − x)c(x) ≥ −C3 M,
v(0) − α0 v 0 (0) = C3 (1 + α0 ) ≥ C3 ,
v(1) = 0.
It follows from the definition of C3 and Theorem 3.1 that v ± w ≥ 0, i.e., |w(x)| ≤ v(x) for all x. Hence |u(1) − u(x)| v(x) |u0 (1)| = lim ≤ lim = C3 , − − 1−x x→1 x→1 1 − x which completes the proof of (3.14). In the proof of Theorem 3.3(iii) we did not use Lemma 3.1 since it will yield only a crude bound on |u0 (1)| when |M | ≥ 1. Thus when M < 0, the solution u of (1.2) is bounded independently of δ and has no boundary layers. For 0 < M ≤ 1 set
(
σ0 (M, δ) = min
) 1.13 δ 0.13 + exp M 1/(δ−1) , 1−M δ−1
(3.15a)
) δ[exp M 1/(δ−1) − 1] 1 , . 1−M (δ − 1)M 1/(δ−1)
(3.15b)
and ( σ1 (M, δ) = min
In these definitions, if M = 1 then each σi equals the second term in {. . . }, while if M → 1− with δ fixed, then the first term blows up so the second term gives the minimum, and if δ → 1+ with M ∈ (0, 1) fixed, then by L’Hˆ opital’s rule the second term is approximately δ/(δ − 1) so the first term gives the minimum. Lemma 3.2. For 0 ≤ M ≤ 1 one has 0 < Eδ−1,δ (M ) ≤ σ0 (M, δ)
and
0 < Eδ−1,δ+1 (M ) ≤ σ1 (M, δ).
(3.16)
Proof. For 0 ≤ M < 1 we have 0 < Eδ−1,δ+1 (M ) ≤ 1/(1 − M ) by the argument used to prove (2.28). For 0 ≤ M ≤ 1, letting dre denote the smallest integer satisfying dre ≥ r, we also have the alternative bound 1 d δ−1 e−1
X
Eδ−1,δ+1 (M ) ≤
k=0
1 + Γ(2)
2 d δ−1 e−1
X 1 e k=d δ−1
M 1/(δ−1) + Γ(3)
X 2 e k=d δ−1
M 2/(δ−1) + ··· Γ(4)
# 2/(δ−1) 3/(δ−1) 1 1 M M ≤ + 1 · 1/(δ−1) M 1/(δ−1) + + + ··· δ−1 2! 3! M h i δ 1/(δ−1) = exp M − 1 . (3.17) (δ − 1) M 1/(δ−1)
"
3 d δ−1 e−1
15
The arguments for Eδ−1,δ (M ) are similar; observe that 1/Γ(δ) < 1.13 since min1 0 on (0, 1) and the definition of M , we have LB3 (x) = −D∗δ B3 (x) + M B30 (x) + [b(x) − M ]B30 (x) + c(x)B3 (x) = kf k∞ + [M − b(x)]φ0M (x)kf k∞ + c(x)B3 (x) ≥ kf k∞
and B3 (0) − α0 B30 (0) = max{|γ0 |, |γ1 |} + kf k∞ [φM (1) + α1 φ0M (1)] ≥ |γ0 |, B3 (1) + α1 B30 (1) = max{|γ0 |, |γ1 |} ≥ |γ1 |.
It follows from Theorem 3.1 that B3 is a barrier function for ±u, i.e., |u(x)| ≤ B3 (x) on [0, 1]. Then (3.21) is immediate from the properties of φM .
16
Theorem 3.4. Assume that 0 ≤ M ≤ 1. Then (i) (ii) (iii)
kuk∞ ≤ max{|γ0 |, |γ1 |} + kf k∞ [σ1 (M, δ) + α1 σ0 (M, δ)],
|u0 (0)| ≤ (δ − 1)[|γ0 | + max{|γ0 |, |γ1 |}] + δ[(e − 1) + α1 (e + 0.13)]kf k∞ , ( (|γ1 | + kuk∞ )/α1 if α1 > 0, |u0 (1)| ≤ σ0 (M, δ) max{kf k∞ + kck∞ |γ1 |, 2[|γ0 | + |γ1 |]} if α1 = 0.
Proof. (i) Combine Lemmas 3.2 and 3.3. (ii) The boundary condition (1.2b) and the condition (1.3) yield |u0 (0)| ≤ (δ−1)(|γ0 |+kuk∞ ). Now invoke part (i) and observe that σ0 (M, δ) ≤ δ(0.13 + e)/(δ − 1), σ1 (M, δ) ≤ δ(e − 1)/(δ − 1).
(iii) If α1 > 0, then the result follows from the boundary condition (1.2b). Thus, assume that α1 = 0. Set w1 (x) = u(x) − u(1) = u(x) − γ1 . Then by (1.2) we have |Lw1 (x)| = |f (x) − c(x)γ1 | ≤ kf k∞ + kck∞ |γ1 |,
|w1 (0) − α0 w10 (0)| = |γ0 − γ1 | ≤ |γ0 | + |γ1 |,
w1 (1) = 0.
Set v1 (x) = K[φM (1) − φM (x)], where φM is defined in (3.18) and K = max {kf k∞ + kck∞ |γ1 |, [|γ0 | + |γ1 |]/φM (1)} . Then v1 ≥ 0 and v10 ≤ 0, so by (3.19) we get Lv1 (x) = −D∗δ v1 (x) + M v10 (x) + [b(x) − M ]v10 (x) + c(x)v1 (x) ≥ K,
v1 (0) − α0 v10 (0) = K[φM (1) − φM (0) + α0 φ0M (0)] = KφM (1),
v1 (1) = 0.
It follows from the definition of K and Theorem 3.1 that v1 ± w1 ≥ 0, i.e., |w1 (x)| ≤ v1 (x) for all x. Hence |u(1) − u(x)| v1 (x) |u0 (1)| = lim ≤ lim = Kφ0M (1). (3.22) − − 1−x x→1 x→1 1 − x
But 0 < φ0M (1) < σ0 (M, δ) by Lemma 3.2, and φM (1) = Eδ−1,δ+1 (M ) ≥ 1/2 on taking the first term in the Mittag-Leffler series. Combining these inequalities with (3.22) and the definition of K completes the proof. If instead of 0 ≤ M ≤ 1 one has the stronger hypothesis that kbk∞ ≤ 1, then a bound similar to that for the case α1 = 0 in Theorem 3.4(iii) can be derived quickly using Lemma 3.1, Theorem 3.4(i)(ii) and (2.27). For 0 ≤ M < 1, Theorem 3.4(i) bounds kuk∞ independently of δ. Theorem 3.4(ii) shows that for 0 ≤ M ≤ 1, when δ is near 1 there is no boundary layer in u at x = 0. The situation in Theorem 3.4(iii) is more complicated; to clarify it, we give now a simpler but slightly less sharp corollary of this part of the theorem. Corollary 3.1. Assume that 0 ≤ M ≤ 1. Then for some constant C we have ( C/(1 − M ) if 0 ≤ M < 1, 0 |u (1)| ≤ C/(δ − 1) if M = 1. Proof. Substitute the bound of Theorem 3.4(i) into Theorem 3.4(iii) and recall the definitions of σ1 (M, δ) and σ0 (M, δ).
17
Theorem 3.4(ii) shows that as δ → 1+ there is no boundary layer at x = 0 when 0 ≤ M ≤ 1. Corollary 3.1 shows that as δ → 1+ , no boundary layer is present at x = 1 when 0 ≤ M < 1, but when M = 1 a boundary layer is possible there and the discussion in Section 2.2.2 shows that then the estimate |u0 (1)| ≤ C/(δ − 1) of Corollary 3.1 is sharp in general.
In Theorem 3.4 and Corollary 3.1 we state no result for the case M > 1, despite the fact that our arguments can be extended to this case, because the bound provided by Lemma 3.3 is then excessively large. The next result gives a satisfactory bound on kuk∞ when M > 1 under the extra hypothesis that b > 0 on [0, 1]. In the definition of B4 below, the expression multiplying kf k∞ imitates closely the coefficient of f that we reported in (2.16). Theorem 3.5. Assume that M > 1 and b ≥ b > 0 for some constant b. For 0 ≤ x ≤ 1, set B4 (x) = max{|γ0 |, |γ1 |} +
kf k∞ {α0 [φM (1) − φM (x)] + xφM (1) − φM (x) + α1 [(α0 + x)φ0M (1) − φM (x)]} . 1 + α0 + bφM (1) + α1 [1 + bφ0M (1)] (3.23)
Then |u(x)| ≤ B4 (x) for x ∈ [0, 1] and for some constant C one has (i) (ii)
kuk∞ ≤ kB4 k∞ ≤ C min{α0 , φM (1)},
(iii) (iv)
|u0 (0)| ≤ C,
|u0 (1)| ≤ C min{α0 , φM (1)} if α1 > 0,
|u0 (1)| ≤ CM 1/(δ−1) min{α0 , φM (1)} if α1 = 0.
Proof. The expression xφM (1) − φM (x) vanishes at x = 0, 1 and its second-order derivative is negative by (2.4), so the expression is positive on (0, 1). It is easy to see that φ0M (1) > φM (1) so we also have xφ0M (1) − φM (x) > 0 on (0, 1). Consequently B4 (x) ≥ 0 for 0 ≤ x ≤ 1. Now
B4 (0) − α0 B40 (0) = max{|γ0 |, |γ1 |} ≥ |γ0 |
and B4 (1) + α1 B40 (1) = max{|γ0 |, |γ1 |} ≥ |γ1 |.
For 0 < x < 1, by (3.19) one has LB4 (x) =c(x)B4 (x) +
kf k∞ {α0 + b(x)φM (1) + 1 + α1 [b(x)φ0M (1) + 1] + (M − b)φ0M (x)(1 + α0 + α1 )} 1 + α0 + bφM (1) + α1 [1 + bφ0M (1)]
≥kf k∞ . Thus B4 is a barrier function for ±u by Theorem 3.1. Hence, recalling that α0 ≥ 1/(δ − 1) > 1, for some C we get kuk∞ ≤ kB4 k∞
kf k∞ {α0 φM (1) + φM (1) + α1 (α0 + 1)φ0M (1)} ≤ max{|γ0 |, |γ1 |} + 1 + α0 + bφM (1) + α1 [1 + bφ0M (1)] α0 φM (1) + α0 α1 φ0M (1) ≤C 1+ α0 + φM (1) + α1 φ0M (1) ≤ C min{α0 , φM (1)}. 18
This proves (i). For (ii) and (iii), use the result of (i) and the boundary conditions (1.2b). To prove (iv), define the function u ˜(x) = u(x) − (1 − x)
γ0 − γ1 − γ1 , 1 + α0
(3.24)
which is the solution of the problem L˜ u(x) = f˜(x) for x ∈ (0, 1),
u ˜(0) − α0 u ˜0 (0) = 0,
u ˜(1) = 0,
with
γ0 − γ1 f˜(x) = f (x) + [b(x) − (1 − x)c(x)] − γ1 c(x). 1 + α0 Applying the barrier function B4 of (i) above to u ˜ yields |˜ u(x)| ≤
kf˜k∞ {α0 [φM (1) − φM (x)] + xφM (1) − φM (x)} , 1 + α0 + bφM (1)
for x ∈ [0, 1].
Hence |˜ u0 (1)| = lim
x→1−
kf˜k∞ [(α0 + 1)φ0M (1) + φM (1)] |˜ u(x)| |˜ u(1) − u ˜(x)| = lim ≤ . 1−x 1 + α0 + bφM (1) x→1− 1 − x
Recalling the identity (2.27), we have
φM (1) = Eδ−1,δ+1 (M ) = [Eδ−1,2 (M ) − 1]/M, φ0M (1) = Eδ−1,δ (M ) = [Eδ−1,1 (M ) − 1]/M
and (2.21) implies that Eδ−1,1 (M )/Eδ−1,2 (M ) ≈ M 1/(δ−1) as δ → 1+ . Consequently |˜ u0 (1)| ≤
Cα0 φ0M (1) Cα0 M 1/(δ−1) φM (1) ≤ ≤ CM 1/(δ−1) min{α0 , φM (1)} . α0 + φM (1) α0 + φM (1)
Part (iv) of the theorem follows from this estimate and (3.24). Thus, in the case M > 1 with b > 0, Theorem 3.5 (ii) shows that there is no boundary layer at x = 0 when δ is near 1, even though kuk∞ may be unbounded as δ → 1+ (see Section 2.2.1).
The accurate approximations (2.22) and (2.23) show that Theorem 3.5 (iii)(iv) is sharp for constant b > 1.
The next example shows numerically that the bounds of Theorem 3.5 (i)(iii) are sharp for a variable-coefficient problem when α0 = 1/(δ − 1) < φM (1). Example 3.2. Consider the test problem −D∗δ u(x) + (x + 0.2)u0 (x) = 2(3 + x) for x ∈ (0, 1), 1 u0 (0) = 0.4, u(1) + u0 (1) = 1.7. u(0) − δ−1
(3.25a) (3.25b)
The exact solution of this problem is unknown. Observe that M = 1.2 so φM (1) 1/(δ−1) = α0 by an inequality similar to (2.19). To check the bounds of Theorem 3.5 (i)(iii), we use the finite difference scheme of Stynes and Gracia (2014) to compute an approximate numerical solution N {uj }N j=0 of (3.25) on a uniform mesh {xj := j/N }j=0 with N = 2048 for various values of δ close to 1, then evaluate (δ − 1) max0≤j≤N |uj | and (δ − 1)|uN − uN −1 |/(xN − xN −1 ) to approximate (δ − 1)kuk∞ and (δ − 1)|u0 (1)|.
As each row of the table is approximately constant as δ → 1+ , the bounds of Theorem 3.5 (i)(iii) are sharp for this example. 19
Table 3.1: Verifying sharpness of Theorem 3.5 (i)(iii) (δ − 1) max0≤j≤N |uj | (δ − 1)|uN − uN −1 |/(xN − xN −1 )
4
δ = 1.1
δ = 1.01
δ = 1.001
δ = 1.0001
δ = 1.00001
δ = 1.000001
δ = 1.0000001
5.8374 5.0045
7.6438 7.5954
7.6043 7.5795
7.6004 7.5777
7.6000 7.5775
7.6000 7.5775
7.6000 7.5775
Conclusions
We considered a two-point boundary value problem whose leading term is a Caputo fractional derivative of order δ with 1 < δ < 2. The dependence of the solution on the parameter δ has not previously been investigated analytically, despite a growing interest in the research literature in the numerical solution of problems with variable fractional derivatives. By considering first the special case of a constant-coefficient operator, for which the solution can be determined explicitly, we showed that when δ is near 1, the solution of the boundary value problem may exhibit a boundary layer at the endpoint x = 1 of the domain. Moving on to the general case of a variable-coefficient differential operator, we then determined conditions on the data of the problem under which boundary layers at each endpoint (x = 0, 1) cannot occur. This analysis showed that a crucial parameter in the presence or absence of a boundary layer at x = 1 is the quantity M := maxx∈[0,1] b(x), where b is the coefficient of the first-order term in the differential operator. In all cases considered, we showed that |u0 (0)| ≤ C, i.e., no boundary layer in u appears at x = 0 when δ is near 1. The only data regime where this bound is not guaranteed by our theory (Theorems 3.2 and 3.5) is when min[0,1] c(x) = 0 and M > 1 without the additional property that b > 0 on [0, 1], but our numerical experience (using the finite difference method of Stynes and Gracia (2014)) is that in this case also no boundary layer appears at x = 0. At x = 1, our theory proves rigorously that when δ is near 1, no boundary layer appears in u if M < 1, but one can have such layers when M ≥ 1. This agrees with our numerical experience. This analysis of the solution u of (1.2) leads naturally to the question: can one construct a numerical method that will yield an accurate approximation of u when it has a boundary layer at x = 1?
References Mohammed Al-Refai. Basic results on nonlinear eigenvalue problems of fractional order. Electron. J. Differ. Equ., 2012(191):1–12, 2012. Hermann Brunner. Collocation methods for Volterra integral and related functional differential equations, volume 15 of Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, Cambridge, 2004. ISBN 0-521-80615-1. doi: 10.1017/CBO9780511543234. Chang-Ming Chen, F. Liu, K. Burrage, and Y. Chen. Numerical methods of the variable-order Rayleigh-Stokes problem for a heated generalized second grade fluid with fractional derivative. IMA J. Appl. Math., 78(5):924–944, 2013. ISSN 0272-4960. doi: 10.1093/imamat/hxr079. Kai Diethelm. The analysis of fractional differential equations, volume 2004 of Lecture Notes in Mathematics. Springer-Verlag, Berlin, 2010. ISBN 978-3-642-14573-5. An application-oriented exposition using differential operators of Caputo type. 20
Jos´e Luis Gracia and Martin Stynes. Upwind and central difference approximation of convection in Caputo fractional derivative two-point boundary value problems. J. Comput. Appl. Math., 273:103–115, 2015. Bangti Jin, Raytcho Lazarov, and Joseph Pasciak. Variational formulation of problems involving fractional order differential operators. ArXiv:1307.4795, 2013. Natalia Kopteva and Martin Stynes. An efficient collocation method for a Caputo two-point boundary value problem. 2014. (Submitted for publication). J. Tenreiro Machado, Virginia Kiryakova, and Francesco Mainardi. Recent history of fractional calculus. Commun. Nonlinear Sci. Numer. Simul., 16(3):1140–1153, 2011. ISSN 1007-5704. A. M. Mathai, R. K. Saxena, and H. J. Haubold. A certain class of Laplace transforms with applications to reaction and reaction-diffusion equations. Astrophys. Space Sci., 305(3):283– 288, 2006. doi: 10.1007/s10509-006-9188-7. Arvet Pedas and Enn Tamme. Piecewise polynomial collocation for linear boundary 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. Igor Podlubny. Fractional differential equations, volume 198 of Mathematics in Science and Engineering. Academic Press Inc., San Diego, CA, 1999. ISBN 0-12-558840-2. An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. John Paul Roop. Numerical approximation of a one-dimensional space fractional advection-dispersion equation with boundary layer. Comput. Math. Appl., 56(7): 1808–1819, 2008. ISSN 0898-1221. doi: 10.1016/j.camwa.2008.04.025. URL http://dx.doi.org/10.1016/j.camwa.2008.04.025. Hans-G¨org Roos, Martin Stynes, and Lutz Tobiska. Robust numerical methods for singularly perturbed differential equations, volume 24 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 2008. ISBN 978-3-540-34466-7. Convection-diffusionreaction and flow problems. Stefan Samko. Fractional integration and differentiation of variable order: an overview. Nonlinear Dynam., 71(4):653–662, 2013. ISSN 0924-090X. doi: 10.1007/s11071-012-0485-0. Martin Stynes and Jos´e Luis Gracia. A finite difference method for a two-point boundary value problem with a Caputo fractional derivative. IMA J. Numer. Anal., 2014. doi: 10.1093/imanum/dru011.
21
Abstract A two-point boundary value problem is considered on the interval [0, 1], where the leading term in the differential operator is a Caputo fractional derivative of order δ with 1 < δ < 2. Writing u for the solution of the problem, it is known that typically u00 (x) blows up as x → 0. A numerical example demonstrates the possibility of a further phenomenon that imposes difficulties on numerical methods: u may exhibit a boundary layer at x = 1 when δ is near 1. The conditions on the data of the problem under which this layer appears are investigated by first solving the constant-coefficient case using Laplace transforms, determining precisely when a layer is present in this special case, then using this information to enlighten our examination of the general variable-coefficient case (in particular, in the construction of a barrier function for u). This analysis proves that usually no boundary layer can occur in the solution u at x = 0, and that the quantity M = maxx∈[0,1] b(x), where b is the coefficient of the first-order term in the differential operator, is critical: when M < 1, no boundary layer is present when δ is near 1, but when M ≥ 1 then a boundary layer at x = 1 is possible. Numerical results illustrate the sharpness of most of our results.
Keywords: Fractional differential equation, Caputo fractional derivative, boundary value problem, boundary layer. MSC (2010) Classification: Primary 34B08 Abbreviated Title: A two-point boundary value problem with a Caputo derivative
1
Introduction
Boundary value problems whose differential operators involve fractional derivatives are of great interest, as these non-classical derivatives can model some physical processes where integer-order ∗
This research was partly supported by the Institute of Mathematics and Applications (IUMA), project MEC/FEDER MTM 2010-16917 and the Diputaci´ on General de Arag´ on. † Email: [email protected] ‡ Email: [email protected]
1
derivatives are unsuitable; see Jin et al. (2013); Machado et al. (2011) for an extensive list of recent applications and mathematical developments in this area. Thus the precise behaviour of solutions to fractional-derivative boundary value problems is of fundamental importance. Let δ ∈ (1, 2). Let g ∈ C 1 [0, 1] with g 0 absolutely continuous on [0, 1]. Then the Caputo fractional derivative D∗δ g associated with the point x = 0 is defined by Z x 1 δ D∗ g(x) := (x − t)1−δ g 00 (t) dt for 0 < x ≤ 1; Γ(2 − δ) t=0 see Diethelm (2010); Pedas and Tamme (2012); Stynes and Gracia (2014). The Riemannδ Liouville fractional derivative DRL of order δ associated with the point x = 0 is defined by Z x d2 1 δ DRL g(x) = 2 (x − t)1−δ g(t) dt for 0 < x ≤ 1; dx Γ(2 − δ) t=0 see Diethelm (2010). These fractional derivatives are related by the formula δ D∗δ g(x) = DRL g(x) −
g 0 (0) g(0) x−δ − x1−δ ; Γ(1 − δ) Γ(2 − δ)
(1.1)
see Diethelm (2010, Lemma 3.4). In the present paper we shall consider the two-point boundary value problem Lu(x) := −D∗δ u(x) + b(x)u0 (x) + c(x)u(x) = f (x) for x ∈ (0, 1), 0
u(0) − α0 u (0) = γ0 ,
0
u(1) + α1 u (1) = γ1 ,
(1.2a) (1.2b)
where 1 < δ < 2. The constants α0 , α1 , γ0 , γ1 and the functions b, c and f are given. We assume that b, c, f ∈ C 1 [0, 1] with c ≥ 0 in [0, 1]. We assume also that α1 ≥ 0
and α0 ≥
1 . δ−1
(1.3)
The conditions on c, α0 and α1 guarantee that (1.2) satisfies a comparison/maximum principle; see Theorem 3.1 below. The problem (1.2) models superdiffusion of particle motion when convection is present; see the discussion and references in Jin et al. (2013, Section 1). It is a member of the general class of boundary value problems that is analysed in Pedas and Tamme (2012). It is also discussed in Al-Refai (2012). Numerical methods for its solution are presented in Gracia and Stynes (2015); Kopteva and Stynes (2014); Stynes and Gracia (2014) and their references. Existence and uniqueness of a classical solution to (1.2) is established in Stynes and Gracia (2014). It is proved that u ∈ C 1 [0, 1] ∩ C 2 (0, 1], and for some constants C˜i one has the sharp bounds ( C˜i if i = 0, 1, (i) |u (x)| ≤ (1.4) δ−i ˜ C2 x if i = 2, for all x ∈ (0, 1]. Thus u00 (x) may blow up at the interval endpoint x = 0.
These results tell us a lot about the nature of the solution u, but in one respect they are seriously deficient: all constants C˜i that appear above depend on the parameter δ, but they can be extremely large when δ is near 1, because in certain cases (as we shall see) the solution u develops a boundary layer at x = 1 (i.e., |u0 (1)| becomes very large) when δ is near 1. It is well known that, when computing numerical solution of problems with integer-order differential 2
operators, such layers can cause a deterioration in accuracy (Roos et al., 2008). This is also the case for the fractional-derivative problem (1.2): see Gracia and Stynes (2015); Jin et al. (2013); Stynes and Gracia (2014), where computed solutions of (1.2) become less accurate when δ is near 1. This loss of accuracy appears in only some numerical examples in these papers, and no explanation is given there, but it is in fact confined to problems whose solutions exhibit a boundary layer at x = 1. In the present paper we shall cast light on when such layers appear in solutions of (1.2). For a concrete example exhibiting a boundary layer, consider (1.2) with b ≡ 1.9, c ≡ 0 and f ≡ 1. Take α0 = 1/(δ − 1), α1 = 0, γ0 = 0.4 and γ1 = 1.7. For constant-coefficient problems like this, an explicit formula for u(x) will be derived in Section 2 below. Using this formula, we plot the solution for the values δ = 1.6, 1.4, 1.2, 1.1 in Figure 1. It is clear from this figure that a boundary layer at x = 1 develops in the solution when δ is near 1. 2.05 2
1.7
1.95
Exact solution u
Exact solution u
1.6
1.5
1.4
1.9 1.85 1.8 1.75 1.7
1.3
1.65 1.2
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
x
0.6
0.8
1
0.6
0.8
1
x
3.6 3.4
6
3.2
5.5 5 Exact solution u
Exact solution u
3 2.8 2.6 2.4 2.2
4.5 4 3.5 3
2
2.5
1.8
2
1.6 0
0.2
0.4
0.6
0.8
1
0
x
0.2
0.4 x
Figure 1: Exact solution for b ≡ 1.9, c ≡ 0, f ≡ 1, α0 = 1/(δ − 1), α1 = 0, γ0 = 0.4 and γ1 = 1.7, with δ = 1.6 (1st row, left), δ = 1.4 (1st row, right), δ = 1.2 (2nd row, left) and δ = 1.1 (2nd row, right), showing development of a boundary layer when δ is near 1 Figure 1 of Roop (2008) also shows a boundary layer developing as the order of the fractional derivative approaches a certain limiting value, though the boundary value problem under discussion there is not the same as (1.2). Furthermore, there is currently great interest in the construction of numerical methods for problems with variable order of fractional derivative (see, e.g., Chen et al. (2013); Samko (2013)); this implies that one should pay close attention to how the solution of (1.2) depends on the value of the parameter δ. For these two reasons (loss of accuracy in numerical solution only for certain values of δ; 3
design of numerical methods for variable δ) our aim in the present paper is to investigate in detail how u behaves as a function of δ (as well as a function of x). The structure and main results of the paper are as follows. In Section 2 Laplace transforms are used to derive an explicit formula for the solution u of (1.2) in the special case when b is a nonzero constant and c ≡ 0. From this formula we deduce that, when δ is near 1, a boundary layer in u at x = 1 can appear only if b ≥ 1, while u never has a boundary layer at x = 0 even though max[0,1] |u(x)| blows up as δ → 1+ if b ≥ 1. In Section 3 the general case of (1.2) is considered and a comparison/maximum principle is used, with some guidance from Section 2, to explore how u depends on δ. In this general case we find that a boundary layer in u at x = 1 when δ is near 1 is possible only when max[0,1] b(x) ≥ 1. Note here that it is the maximum of b, not of |b|, that is the significant quantity. It is shown also that a boundary layer at x = 0 (when δ is near 1) is possible only if max[0,1] b(x) > 1 and min[0,1] b(x) ≤ 0 and min[0,1] c(x) = 0.
Notation. We use the standard notation C k (I) to denote the space of real-valued functions whose derivatives up to order k are continuous on an interval I, and write C(I) for C 0 (I). For each g ∈ C[0, 1], set kgk∞ = maxx∈[0,1] |g(x)|. In several inequalities C denotes a generic positive constant that depends on the data b, c, f, γ0 , γ1 , α1 of the boundary value problem (1.2) but is independent of δ and x; note that C can take different values in different places. A subscripted C (e.g., C1 ) denotes a fixed positive constant that can depend on all the data of the boundary value problem (1.2) except δ and x.
The case where b is constant and c ≡ 0
2
In Section 2 we consider the special case of problem (1.2) where c ≡ 0, and b is constant with b 6= 0. Our results could be extended to the case where b and c are arbitrary constants, but when c 6= 0 the details become much more complex; see Remark 2.1. We shall use Laplace transforms to derive an explicit formula for the solution u of (1.2) in terms of Mittag-Leffler functions. Our examination of this special case gives useful and penetrating insights into the properties of the solution u of (1.2). Furthermore, the precise form of the solution that we find in Section 2.2 is very helpful when constructing a barrier function in Section 3 to analyse the solution of (1.2) when b and c are no longer constants.
2.1
General right-hand side f
Extend the domain of f to [0, ∞) in such a way that the extension (which we also call f ) is smooth and has support in [0, 2]. To solve (1.2), we treat it as an initial-value problem on [0, ∞) with the initial condition u(0) − α0 u0 (0) = Rγ0 from (1.2b), and apply the standard Laplace ∞ transform operator L, defined by Lv(s) = t=0 e−st v(t) dt. The second boundary condition u(1) + α1 u0 (1) = γ1 of (1.2b) will be invoked later. Define the two-parameter Mittag-Leffler function by Eα,β (z) =
∞ X k=0
zk Γ(αk + β)
for α, β, z ∈ R with α > 0.
(2.1)
We shall need the properties (Podlubny, 1999, (1.80),(1.82)) that for constant α, β, γ and λ one has n o sα−β L xβ−1 Eα,β (±λxα ) = α (2.2) s ∓λ 4
and
γ DRL xβ−1 Eα,β (λxα ) = xβ−γ−1 Eα,β−γ (λxα ).
(2.3)
γ Note that when γ = 1 then DRL ≡ d/dx (Diethelm, 2010, p.27); thus (2.3) implies that
d β−1 x Eα,β (λxα ) = xβ−2 Eα,β−1 (λxα ). (2.4) dx Applying L to (1.2a) and observing that L{D∗δ u} = sδ−2 s2 L{u} − su(0) − u0 (0) from Podlubny (1999, (2.253)), we obtain (−sδ + bs)L{u} + sδ−1 u(0) + sδ−2 u0 (0) − bu(0) = L{f }. Hence L{u} = −
L{f } u(0) sδ−3 u0 (0) L{f } γ0 + α0 u0 (0) sδ−3 u0 (0) + + δ−1 + = − + s s s(sδ−1 − b) sδ−1 − b s(sδ−1 − b) s −b
(2.5)
using the boundary condition (1.2b) at x = 0. To find the inverse Laplace transform of −L{f }/[s(sδ−1 − b)] we imitate Diethelm (2010, p.135). Using the integration theorem for Laplace transforms, R x L r=0 f (r)dr L{f } = . (2.6) s(sδ−1 − b) sδ−1 − b By (2.4) one has d Eδ−1,1 (bxδ−1 ) = x−1 Eδ−1,0 (bxδ−1 ), (2.7) dx where the first term in the Mittag-Leffler series for Eδ−1,0 (bxδ−1 ) vanishes since Γ(0) = ∞. But (2.2) yields sδ−2 L{Eδ−1,1 (bxδ−1 )} = δ−1 s −b so the differentiation theorem for Laplace transforms gives n o n o L x−1 Eδ−1,0 (bxδ−1 ) = sL Eδ−1,1 (bxδ−1 ) − Eδ−1,1 (0) =
b sδ−1 − 1 = δ−1 . δ−1 s −b s −b
Consequently from (2.6) we have o L{f } 1 n −1 δ−1 = L x E (bx ) L δ−1,0 b s(sδ−1 − b)
Z
x
f (r)dr .
r=0
Now the convolution theorem for Laplace transforms yields Z x−t Z L{f } 1 x −1 −1 δ−1 L (x) = t Eδ−1,0 (bt ) f (r) dr dt. b t=0 s(sδ−1 − b) r=0
(2.8)
Taking the inverse transform of (2.5) and invoking (2.8) and (2.2), we get Z x−t Z 1 x −1 t Eδ−1,0 (btδ−1 ) f (r) dr dt. (2.9) u(x) = γ0 + α0 u0 (0) + u0 (0)xEδ−1,2 (bxδ−1 ) − b t=0 r=0 By virtue of (2.4) one can differentiate (2.9) to obtain Z 1 x −1 0 0 δ−1 u (x) = u (0)Eδ−1,1 (bx ) − t Eδ−1,0 (btδ−1 )f (x − t) dt. b t=0 5
(2.10)
Consequently, imposing the boundary condition u(1) + α1 u0 (1) = γ1 of (1.2b), one has Z 1−t Z 1 1 −1 0 0 δ−1 γ1 =γ0 + α0 u (0) + u (0)Eδ−1,2 (b) − t Eδ−1,0 (bt ) f (r) dr dt b t=0 r=0 Z 1 1 −1 + α1 u0 (0)Eδ−1,1 (b) − t Eδ−1,0 (btδ−1 )f (1 − t) dt b t=0 whence u0 (0) =
γ1 − γ0 +
1 b
R1
−1 δ−1 ) t=0 t Eδ−1,0 (bt
hR 1−t r=0
i f (r) dr + α1 f (1 − t) dt
α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b)
.
(2.11)
Substitution of (2.11) into (2.9) and (2.10) yields explicit formulas for u(x) and u0 (x).
2.2
Constant right-hand side f
In this section we simplify the results of Section 2.1 by taking f to be constant so that we can then investigate in detail the solution u. Furthermore, the formulas of Section 2.2 will be of great help in the construction of a barrier function in Section 3 to analyse the structure and behaviour of u in the case of variable b, c and f . First, taking f constant in (2.11) and recalling (2.7), we have R1 d Eδ−1,1 (btδ−1 ) dt γ1 − γ0 + fb t=0 [α1 + (1 − t)] dt 0 u (0) = α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b) i h R1 γ1 − γ0 + fb α1 Eδ−1,1 (b) − α1 − 1 + t=0 Eδ−1,1 (btδ−1 ) dt = α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b) =
γ1 − γ0 + fb [α1 Eδ−1,1 (b) − α1 − 1 + Eδ−1,2 (b)] α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b)
(2.12)
where we integrated by parts then invoked (2.4) with β = 2. Substituting (2.12) into (2.9) gives us u(x); this formula can be written in a variety of ways, the most compact of which seems to be α0 + xEδ−1,2 (bxδ−1 ) [γ1 − γ0 − (1 + α0 + α1 )f /b] (α0 + x)f u(x) = γ0 + + . (2.13) b α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b) The correctness of this formula can be verified by substitution into (1.2), using (2.4) and, by (1.1) and (2.3), D∗δ xEδ−1,2 (bxδ−1 ) = x1−δ Eδ−1,2−δ (bxδ−1 ) − = x1−δ
∞ X k=1
x1−δ Γ(2 − δ)
(bxδ−1 )k Γ(k(δ − 1) + 2 − δ)
= bEδ−1,1 (bxδ−1 ). Remark 2.1. In the case of constant b, c, f with c 6= 0, one can use the above Laplace transform technique to obtain a closed-form representation of the solution by imitating (Mathai et al., 2006, (35)). For example, to invert [s(sδ − bs − c)]−1 , write r X ∞ ∞ X 1 bs br sr−1 1 1 = = = δ δ δ δ s(s − bs − c) s(s − c) s −c (s − c)r+1 s(sδ − c) 1 − sδbs−c r=0 r=0 6
then appeal to (Podlubny, 1999, (1.80)): n o (k) L xαk+β−1 Eα,β (±λxα ) =
k! sα−β (sα ∓ λ)k+1
(2.14)
(k)
where Eα,β (y) = (dk /dy k )Eα,β (y). This will yield an infinite sum of derivatives of Mittag-Leffler functions, which is very complicated. As our main aim in Section 2 is to gain insight into the solution of the general problem (1.2), we do not consider c 6= 0 any further. Remark 2.2. Differentiating (2.13) twice by invoking (2.4), one gets γ1 − γ0 − (1 + α0 + α1 )f /b 1 00 u (x) = Eδ−1,0 (bxδ−1 ). α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b) x
(2.15)
When γ1 − γ0 6= (1 + α0 + α1 )f /b, on recalling the definition (2.1) and Γ(0) = ∞ it follows ˜ δ−2 . This coincides with the bound on |u00 (x)| that is that as x → 0 one has u00 (x) ∼ Cx derived in Stynes and Gracia (2014, Corollary 3.5) for the general case of variable b, c and f . Furthermore, the next term in the series expansion of the Mittag-Leffler function in (2.15) is ˜ −1 (xδ−1 )2 = Cx ˜ 2δ−3 which matches Stynes and Gracia (2014, Example 3.7), where it was Cx shown that the derivative bounds (1.4) are sharp. Remark 2.3. The coefficient of f in the formula (2.13) for u(x) is (α0 + x)[α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b)] − [α0 + xEδ−1,2 (bxδ−1 )](1 + α0 + α1 ) b[α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b)] n 1 = α0 [Eδ−1,2 (b) − 1 − xEδ−1,2 (bxδ−1 ) + x] b[α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b)] o +xEδ−1,2 (b) − xEδ−1,2 (bxδ−1 ) + α1 [(α0 + x)Eδ−1,1 (b) − α0 − xEδ−1,2 (bxδ−1 )] =
n 1 α0 [Eδ−1,δ+1 (b) − xδ Eδ−1,δ+1 (bxδ−1 )] 1 + α0 + bEδ−1,δ+1 (b) + α1 [1 + bEδ−1,δ (b)] o +xEδ−1,δ+1 (b) − xδ Eδ−1,δ+1 (bxδ−1 ) + α1 [(α0 + x)Eδ−1,δ (b) − xδ Eδ−1,δ+1 (bxδ−1 )] ,
(2.16)
where we used the elementary identity Eδ−1,1+i (z) = zEδ−1,δ+i (z) + 1
for i = 0, 1.
(2.17)
The formula (2.16) is the inspiration for our choice of barrier function in Theorem 3.5. We turn now to our main interest: investigating how the solution u changes as δ approaches 1, and in particular determining when boundary layers appear in u. Observe first that, since α0 > 1 by condition (1.3), it follows from (2.12) that |u0 (0)| ≤ C for some constant C, i.e., there is never a boundary layer in u at x = 0. In the subsections that follow, we show that when δ is near 1, the magnitudes of kuk∞ and depend strongly on whether b < 1 or b ≥ 1.
|u0 (1)|
Differentiating (2.13) gives
u0 (1) =
(γ1 − γ0 )Eδ−1,1 (b) f α0 [1 − Eδ−1,1 (b)] + Eδ−1,2 (b) − Eδ−1,1 (b) · + . (2.18) b α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b) α0 + Eδ−1,2 (b) + α1 Eδ−1,1 (b) 7
2.2.1
Case b > 1
Assume in Section 2.2.1 that b > 1 and f 6= 0. Suppose that 1/(δ − 1) = K for some positive integer K. Then for i = 1, 2 we have Eδ−1,i (b) =
∞ X k=0
bk bK 1 ≥ ≥ b1/(δ−1) . Γ((δ − 1)k + i) Γ(1 + i) 2
(2.19)
It follows that limδ→1+ Eδ−1,i (b) = ∞ and Eδ−1,i (b) 1/(δ − 1) as δ → 1+ . Consequently, when α0 = 1/(δ − 1), by virtue of (2.12) one has limδ→1+ u0 (0) = f /b and hence, recalling (1.2b), ( ∞ if f > 0, lim u(0) = lim [γ0 + α0 u0 (0)] = + + δ→1 δ→1 −∞ if f < 0. Thus when b > 1 and f 6= 0 we do not have kuk∞ ≤ C independently of δ. Nevertheless, recalling from above that limδ→1+ u0 (0) = f /b, there is no boundary layer at x = 0. To discuss limδ→1+ u0 (1), we use more sophisticated machinery. Set σ = δ − 1 for brevity. Let γ(1, ϕ) be the complex-plane contour of Figure 2 where we choose ϕ = 2σπ/3. This contour divides the complex plane into two regions, which we denote by G− (1, ϕ) and G+ (1, ϕ). The real number b > 1 lies in the region G+ (1, ϕ), so by (Podlubny, 1999, Theorem 1.1), for arbitrary but fixed n > 0 we have Z 1 (1−n)/σ 1 ζ (1−n)/σ exp(ζ 1/σ ) 1/σ Eσ,n (b) = b exp(b ) + dζ, (2.20) σ 2πσı γ(1,ϕ) ζ −b where ı is the imaginary unit.
y
G− (1, ϕ) G+ (1, ϕ)
ϕ 1 b
O γ(1, ϕ)
x
−ϕ
Figure 2: Contour γ(1, ϕ) of equation (2.20) Our interest lies in what happens when δ → 1+ (i.e., σ → 0+ , so ϕ → 0+ ) with b > 1 fixed. For ζ ∈ γ(1, ϕ) and ϕ sufficiently small we have |ζ − b| ≥ min{b − 1, b sin ϕ} = b sin ϕ ≥ bϕ/2 = bσπ/3. 8
Consequently (cf. (Podlubny, 1999, p.33)), for some positive constants C one obtains Z 1 Z (1−n)/σ exp(ζ 1/σ ) ζ C C |ζ (1−n)/σ | exp |ζ|1/σ cos(2π/3) dζ ≤ 2 dζ ≤ 2 2πσı γ(1,ϕ) σ γ(1,ϕ) ζ −b σ as σ → 0+ , where the contour integral is bounded since cos(2π/3) < 0. Thus (2.20) implies that for b > 1 one has 1 (1−n)/(δ−1) 1 1/(δ−1) Eδ−1,n (b) = as δ → 1+ . (2.21) b exp b +O δ−1 (δ − 1)2 It follows from (2.21) that Eδ−1,1 (b)/Eδ−1,2 (b) ≈ b1/(δ−1) when δ is near 1. Consequently, when δ is near 1, if α1 > 0 then (2.18) and (2.19) imply that |u0 (1)| ≈ C min{α0 , Eδ−1,1 (b)},
(2.22)
while if α1 = 0, then (2.18) and (2.19) imply that |u0 (1)| ≈
Cα0 b1/(δ−1) Eδ−1,2 (b) Cα0 Eδ−1,1 (b) ≈ ≈ Cb1/(δ−1) min{α0 , Eδ−1,2 (b)}. α0 + Eδ−1,2 (b) α0 + Eδ−1,2 (b)
(2.23)
Thus when b > 1 there is a boundary layer at x = 1, and it is much stronger when α1 = 0. 2.2.2
Case b = 1
Assume in Section 2.2.2 that b = 1 and f 6= 0. Then (2.13) yields α0 [γ1 − γ0 − (1 + α0 + α1 )f ] u(0) = γ0 + α0 f + α0 + Eδ−1,2 (1) + α1 Eδ−1,1 (1) γ1 − γ0 + f [Eδ−1,2 (1) + α1 Eδ−1,1 (1) − 1 − α1 ] . = γ0 + α0 α0 + Eδ−1,2 (1) + α1 Eδ−1,1 (1)
(2.24)
Suppose that 1/(δ − 1) is an integer. Then Eδ−1,2 (1) =
∞ X k=0
1
2
−1
3
−1
−1
δ−1 δ−1 δ−1 X 1 X 1 X 1 1 ≥ + + + ··· Γ((δ − 1)k + 2) Γ(3) Γ(4) Γ(5) 1 2
k=0
k= δ−1
1 1 1 1 + + + ··· δ − 1 Γ(3) Γ(4) Γ(5) e−2 = . δ−1
k= δ−1
=
(2.25)
Likewise, one has 1 −1 δ−1
Eδ−1,2 (1) ≤
X k=0
1 + Γ(2)
2 −1 δ−1
X 1 k= δ−1
3 −1 δ−1
X 1 1 + + ··· Γ(3) Γ(4) 2 k= δ−1
1 1 1 1 = + + + ··· δ − 1 Γ(2) Γ(3) Γ(4) e−1 = . δ−1 9
(2.26)
One can show similarly that
e−1 e ≤ Eδ−1,1 (1) ≤ . δ−1 δ−1
(2.27)
Invoking (2.25)–(2.27) in (2.24), for f > 0 and δ sufficiently close to 1 one obtains h i 1 (e−1) γ1 − γ0 + f e−2+α − 1 − α 1 δ−1 u(0) ≥ γ0 + α0 1e α0 + e−1+α δ−1 h i 1 (e−1) γ1 − γ0 + f e−2+α − 1 − α1 δ−1 = γ0 + . 1e 1 + e−1+α α0 (δ−1) But α0 ≥ 1/(δ − 1), so it follows that as δ → 1+ with 1/(δ − 1) an integer, then u(0) → ∞ if f > 0. Otherwise, if f < 0 then u(0) → −∞ as δ → 1+ . Thus when b = 1 and f 6= 0 we do not have kuk∞ ≤ C independently of δ.
Furthermore, the condition (1.3) compared with (2.26) and (2.27) shows that α0 ≥ 13 Eδ−1,i (1) for i = 1, 2 when δ is near 1. Consequently (2.18) implies that |u0 (1)| ≈ CEδ−1,1 (1) ≈ C/(δ − 1) when δ is near 1; thus u exhibits a boundary layer at x = 1. 2.2.3
Case 0 ≤ b < 1
Assume in Section 2.2.3 that 0 ≤ b < 1. Then Eδ−1,2 (b) =
∞ X k=0
∞
X bk 1 bk = ≤ . Γ((δ − 1)k + 2) 1−b
(2.28)
k=0
The estimate (2.28) is qualitatively sharp when δ is near 1 because for 0 ≤ k ≤ b1/(δ − 1)c, where bnc denotes the greatest integer less or equal to n, one has 1 1 1 ≥ = Γ((δ − 1)k + 2) Γ(3) 2 so Eδ−1,2 (b) ≥
b1/(δ−1)c k X b k=0
2
=
1 − b1+b1/(δ−1)c 1 ≥ for δ sufficiently close to 1. 2(1 − b) 4(1 − b)
(2.29)
Similarly to (2.28) and (2.29), one has 1 1 ≤ Eδ−1,1 (b) ≤ , 2(1 − b) θ(1 − b)
(2.30)
where θ = min{Γ(x) : 1 ≤ x ≤ 2}.
From (2.28)–(2.30) and (2.12), since α0 ≥ 1/(δ − 1) one obtains |u0 (0)| ≤ C/α0 for some constant C, so u0 (0) → 0 as δ → 1+ . Furthermore, (2.13) and (2.18) now yield kuk∞ + |u0 (1)| ≤ C for 0 ≤ b < 1 for some constant C = C(b). Thus no boundary layers are present in u when 0 ≤ b < 1.
10
(2.31)
Remark 2.4. The analysis above of the cases b > 1, b = 1 and 0 ≤ b < 1 shows that, as a function of δ, the nature of the solution u undergoes a fundamental change when b moves from b ≥ 1 to b < 1. This is shown graphically in Figure 3, where u changes dramatically as b moves from 0.9 to 1.1 (as well as the evident changes in shape in these graphs, note the changes in scale of the y axis). The regime b < 0 is not discussed separately in Section 2 since we are able to handle it easily in Section 3 when we discuss the general variable-coefficient case; see Theorem 3.3.
12 11
6
10
5.5
9 Exact solution u
Exact solution u
7 6.5
5 4.5 4 3.5
8 7 6 5
3
4
2.5
3
2
2 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
x
0.6
0.8
1
x
18 16
Exact solution u
14 12 10 8 6 4 2 0
0.2
0.4
0.6
0.8
1
x
Figure 3: Exact solution for δ = 1.05, c ≡ 0, f ≡ 1, α0 = 1/(δ − 1), α1 = 0, γ0 = 0.4, γ1 = 1.7, and b = 0.9 (top left figure), b = 1.0 (top right figure), b = 1.1 (bottom figure), showing effect on solution u of moving from b < 1 to b > 1
3
Boundary layers in the variable-coefficient problem
In numerical solutions of (1.2) computed by the finite difference method of Stynes and Gracia (2014), when δ is near 1 we have observed boundary layers at x = 1 in certain examples but we have never observed a layer at x = 0. The results of Section 3 will extend those of Section 2, prove in most cases that a boundary layer cannot occur at x = 0, and provide substantial information about when boundary layers can appear at x = 1 when δ is near 1. The key tool in the analysis presented in Section 3 is the following comparison principle, which will be used several times. Theorem 3.1. (Stynes and Gracia, 2014, Theorem 2.1) Let z ∈ C 1 [0, 1] ∩ C 2 (0, 1]. Assume that |z 00 (x)| ≤ Kx−θ for 0 < x ≤ 1, where θ ∈ (0, 1) and K are constants that are independent 11
of x. Let b, c ∈ C[0, 1] with c(x) ≥ 0 for all x ∈ (0, 1). Assume that z satisfies the inequalities − D∗δ z + bz 0 + cz ≥ 0 on (0, 1), 0
(3.1a)
0
z(0) − α0 z (0) ≥ 0 and z(1) + α1 z (1) ≥ 0,
(3.1b)
where α0 and α1 satisfy (1.3). Then z ≥ 0 on [0, 1]. The bounds (1.4) show that the solution u of (1.2) satisfies the regularity hypotheses imposed on z in Theorem 3.1. We shall apply Theorem 3.1 to Bi ± u for various barrier functions Bi , concluding that Bi ± u ≥ 0, i.e., that |u| ≤ Bi on [0, 1]. We begin with a general result that provides a useful relationship between u0 (1) and u0 (0). For any function φ ∈ C[0, 1], set kφk[0,x] = max0≤t≤x |φ(t)| for 0 ≤ x ≤ 1. Lemma 3.1. There exists a constant C such that |u0 (x)| ≤ |u0 (0)| + Cxδ−1 1 + kuk[0,x] + |u0 (0)| Eδ−1,1 (kbk[0,x] xδ−1 )
for 0 ≤ x ≤ 1.
Proof. Set y(x) = u0 (x) − u0 (0) for 0 ≤ x ≤ 1. On multiplying (1.2a) by (x − t)δ−2 /Γ(δ − 1) then integrating from t = 0 to t = x, after some manipulation of the fractional derivative term one obtains (Kopteva and Stynes, 2014) a weakly singular Volterra integral equation of the second kind in the unknown y: for 0 < x ≤ 1, Z x Z x 1 1 y(x) − (x − t)δ−2 b(t)y(t) dt = (x − t)δ−2 {b(t)u0 (0) + c(t)u(t) − f (t)} dt Γ(δ − 1) t=0 Γ(δ − 1) t=0 := g(x), say. (3.2) From Brunner (2004, p.343) the solution of (3.2) can be expressed as Z x y(x) = g(x) + Rα (x, t)g(t) dt for 0 ≤ x ≤ 1,
(3.3)
t=0
where δ−2
Rα (x, t) = (x − t)
∞ X
(x − t)(n−1)(δ−1) Φn (x, t; δ)
(3.4)
n=1
and the Φn are defined iteratively by Φ1 (x, t; δ) := b(t)/Γ(δ − 1) Z 1 1 (1 − z)δ−2 z (n−1)(δ−1)−1 b(t + (x − t)z) Φn−1 (t + (x − t)z, t; δ) dz Φn (x, t; δ) := Γ(δ − 1) z=0 for n = 1, 2, . . . An inductive argument using Euler’s beta function shows easily (cf. Brunner (2004, Lemma 6.1.3)) that |Φn (x, t; δ)| ≤
kbkn[0,x]
Γ(n(δ − 1))
for n = 1, 2, · · · and 0 ≤ t ≤ x ≤ 1.
(3.5)
Consequently the infinite series in (3.4) is uniformly convergent for 0 ≤ t ≤ x ≤ 1. Hence we can move the summation sign of (3.4) outside the integral in (3.3). Again appealing to (3.5),
12
from (3.3) we get |y(x)| ≤ |g(x)| + kgk[0,x] = |g(x)| + kgk[0,x] = |g(x)| + kgk[0,x]
∞ X n=1 ∞ X n=1 ∞ X n=1
kbkn[0,x]
Z
Γ(n(δ − 1))
x
t=0
(x − t)(n−1)(δ−1)+δ−2 dt
kbkn[0,x] xn(δ−1)
n(δ − 1)Γ(n(δ − 1)) (kbk[0,x] xδ−1 )n Γ(n(δ − 1) + 1)
≤ kgk[0,x] Eδ−1,1 (kbk[0,x] xδ−1 )
(3.6)
by the definition (2.1) of the Mittag-Leffler function. Recalling the definition of g in (3.2), it follows that kgk[0,x]
1 Γ(δ − 1) xδ−1 = C 1 + kuk[0,x] + |u0 (0)| Γ(δ) δ−1 0 ≤ C 1 + kuk[0,x] + |u (0)| x 0
≤ C 1 + kuk[0,x] + |u (0)|
Z
x
t=0
(x − t)δ−2 dt
(3.7)
as inf δ∈(1,2) Γ(δ) > 0. Combining (3.6), (3.7) and u0 (x) = u0 (0) + y(x) yields the required result. Inequality (3.6) is sharp if g and b are positive constants; see Brunner (2004, Theorem 6.1.1). First, we show that in the special case where c > 0 on [0, 1], the solution u is uniformly bounded for δ ∈ (1, 2) and no boundary layer appears at x = 0 when δ is near 1; there is also no boundary layer at x = 1 if α1 > 0. If c ≥ c > 0 for some constant c, set kf k∞ C1 = max |γ0 |, |γ1 |, c
.
Theorem 3.2. Assume that c ≥ c > 0 for some constant c. Then (i)
kuk∞ ≤ C1 ,
(3.8)
C1 + |γ0 | ≤ (δ − 1)(C1 + |γ0 |), α0 C1 + |γ1 | (iii) |u0 (1)| ≤ if α1 > 0. α1 (iv) |u0 (1)| ≤ CEδ−1,1 (kbk∞ ) for some constant C if α1 = 0. (ii) |u0 (0)| ≤
(3.9) (3.10) (3.11)
Proof. (i) Define the constant barrier function B1 (x) ≡ C1 . Then B1 ± u ≥ 0 by Theorem 3.1, so (3.8) is valid. (ii) By (1.2b) we have u(0) − α0 u0 (0) = γ0 , so (1.3) and (3.8) yield u(0) − γ0 kuk∞ + |γ0 | C1 + |γ0 | 0 ≤ |u (0)| = ≤ ≤ (δ − 1)(C1 + |γ0 |). α0 α0 α0 (iii) Use the boundary condition u(1) + α1 u0 (1) = γ1 and (3.8) to derive (3.10). (iv) This follows immediately from Lemma 3.1, (3.8) and (3.9). 13
Theorem 3.2 shows that when c is strictly positive on [0, 1] there is no boundary layer at x = 0 and, if one also has α1 > 0, there is no boundary layer at x = 1. But when α1 = 0 one may have a layer at x = 1, as the next example demonstrates. Example 3.1. Consider our boundary value problem (1.2) with δ = 1.01, α0 = 1/(δ − 1), α1 = 0, γ0 = 0.4, γ1 = 1.7 and constant functions b(x) ≡ 1.1, c(x) ≡ 1 and f (x) ≡ 3.25. In Figure 4 we show the solution (values at the mesh points, joined by a piecewise linear curve) computed by the finite difference scheme of Stynes and Gracia (2014) on a uniform mesh of width 1/2048. A 3.2 3
Computed solution
2.8 2.6 2.4 2.2 2 1.8 0
0.1
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
0.9
1
Figure 4: Example showing boundary layer at x = 1 is possible when c > 0 and α1 = 0 boundary layer at x = 1 is clearly visible. Set M = max b(x). x∈[0,1]
We observe that M can have any sign. In Theorems 3.3, 3.4 and 3.5 we shall consider the regimes M < 0, 0 ≤ M ≤ 1 and M > 1 respectively and in each case we shall derive bounds on kuk∞ , |u0 (0)| and |u0 (1)|. The case where b is strictly negative is addressed first. For M < 0, set (1 + α1 )kf k∞ kf k∞ + C2 kck∞ C2 = max |γ0 |, |γ1 | − and C3 = max |γ0 | + C2 , . M −M
Theorem 3.3. Assume that M < 0. Then (i) (ii) (iii)
kuk∞ ≤ C2 , 0
|u (0)| ≤ (δ − 1)(C2 + |γ0 |), 0
|u (1)| ≤ C3 .
(3.12) (3.13) (3.14)
Proof. (i) Set B2 (x) = C2 + xkf k∞ /M for 0 ≤ x ≤ 1. Then B2 ≥ 0 so LB2 (x) = b(x)B20 (x) + c(x)B2 (x) ≥ B2 (0) − α0 B20 (0) = C2 −
b(x)kf k∞ ≥ kf k∞ , M
α0 kf k∞ (1 + α1 )kf k∞ ≥ |γ0 | and B2 (1) + α1 B20 (1) = C2 + ≥ |γ1 | M M
by definition of C2 . It follows from (1.2) and Theorem 3.1 that B2 ± u ≥ 0, i.e., |u(x)| ≤ B2 (x) for 0 ≤ x ≤ 1. As B2 (x) ≤ C2 , we have proved (3.12). 14
(ii) To obtain (3.13), use the boundary condition u(0) − α0 u0 (0) = γ0 , (1.3) and (3.12). (iii) Set w(x) = u(x) − u(1). Then by (1.2) and (3.12) we have
|Lw(x)| = |f (x) − c(x)u(1)| ≤ kf k∞ + C2 kck∞ ,
|w(0) − α0 w0 (0)| = |γ0 − u(1)| ≤ |γ0 | + C2 ,
w(1) = 0.
Set v(x) = C3 (1 − x). Then Lv(x) = −C3 b(x) + C3 (1 − x)c(x) ≥ −C3 M,
v(0) − α0 v 0 (0) = C3 (1 + α0 ) ≥ C3 ,
v(1) = 0.
It follows from the definition of C3 and Theorem 3.1 that v ± w ≥ 0, i.e., |w(x)| ≤ v(x) for all x. Hence |u(1) − u(x)| v(x) |u0 (1)| = lim ≤ lim = C3 , − − 1−x x→1 x→1 1 − x which completes the proof of (3.14). In the proof of Theorem 3.3(iii) we did not use Lemma 3.1 since it will yield only a crude bound on |u0 (1)| when |M | ≥ 1. Thus when M < 0, the solution u of (1.2) is bounded independently of δ and has no boundary layers. For 0 < M ≤ 1 set
(
σ0 (M, δ) = min
) 1.13 δ 0.13 + exp M 1/(δ−1) , 1−M δ−1
(3.15a)
) δ[exp M 1/(δ−1) − 1] 1 , . 1−M (δ − 1)M 1/(δ−1)
(3.15b)
and ( σ1 (M, δ) = min
In these definitions, if M = 1 then each σi equals the second term in {. . . }, while if M → 1− with δ fixed, then the first term blows up so the second term gives the minimum, and if δ → 1+ with M ∈ (0, 1) fixed, then by L’Hˆ opital’s rule the second term is approximately δ/(δ − 1) so the first term gives the minimum. Lemma 3.2. For 0 ≤ M ≤ 1 one has 0 < Eδ−1,δ (M ) ≤ σ0 (M, δ)
and
0 < Eδ−1,δ+1 (M ) ≤ σ1 (M, δ).
(3.16)
Proof. For 0 ≤ M < 1 we have 0 < Eδ−1,δ+1 (M ) ≤ 1/(1 − M ) by the argument used to prove (2.28). For 0 ≤ M ≤ 1, letting dre denote the smallest integer satisfying dre ≥ r, we also have the alternative bound 1 d δ−1 e−1
X
Eδ−1,δ+1 (M ) ≤
k=0
1 + Γ(2)
2 d δ−1 e−1
X 1 e k=d δ−1
M 1/(δ−1) + Γ(3)
X 2 e k=d δ−1
M 2/(δ−1) + ··· Γ(4)
# 2/(δ−1) 3/(δ−1) 1 1 M M ≤ + 1 · 1/(δ−1) M 1/(δ−1) + + + ··· δ−1 2! 3! M h i δ 1/(δ−1) = exp M − 1 . (3.17) (δ − 1) M 1/(δ−1)
"
3 d δ−1 e−1
15
The arguments for Eδ−1,δ (M ) are similar; observe that 1/Γ(δ) < 1.13 since min1 0 on (0, 1) and the definition of M , we have LB3 (x) = −D∗δ B3 (x) + M B30 (x) + [b(x) − M ]B30 (x) + c(x)B3 (x) = kf k∞ + [M − b(x)]φ0M (x)kf k∞ + c(x)B3 (x) ≥ kf k∞
and B3 (0) − α0 B30 (0) = max{|γ0 |, |γ1 |} + kf k∞ [φM (1) + α1 φ0M (1)] ≥ |γ0 |, B3 (1) + α1 B30 (1) = max{|γ0 |, |γ1 |} ≥ |γ1 |.
It follows from Theorem 3.1 that B3 is a barrier function for ±u, i.e., |u(x)| ≤ B3 (x) on [0, 1]. Then (3.21) is immediate from the properties of φM .
16
Theorem 3.4. Assume that 0 ≤ M ≤ 1. Then (i) (ii) (iii)
kuk∞ ≤ max{|γ0 |, |γ1 |} + kf k∞ [σ1 (M, δ) + α1 σ0 (M, δ)],
|u0 (0)| ≤ (δ − 1)[|γ0 | + max{|γ0 |, |γ1 |}] + δ[(e − 1) + α1 (e + 0.13)]kf k∞ , ( (|γ1 | + kuk∞ )/α1 if α1 > 0, |u0 (1)| ≤ σ0 (M, δ) max{kf k∞ + kck∞ |γ1 |, 2[|γ0 | + |γ1 |]} if α1 = 0.
Proof. (i) Combine Lemmas 3.2 and 3.3. (ii) The boundary condition (1.2b) and the condition (1.3) yield |u0 (0)| ≤ (δ−1)(|γ0 |+kuk∞ ). Now invoke part (i) and observe that σ0 (M, δ) ≤ δ(0.13 + e)/(δ − 1), σ1 (M, δ) ≤ δ(e − 1)/(δ − 1).
(iii) If α1 > 0, then the result follows from the boundary condition (1.2b). Thus, assume that α1 = 0. Set w1 (x) = u(x) − u(1) = u(x) − γ1 . Then by (1.2) we have |Lw1 (x)| = |f (x) − c(x)γ1 | ≤ kf k∞ + kck∞ |γ1 |,
|w1 (0) − α0 w10 (0)| = |γ0 − γ1 | ≤ |γ0 | + |γ1 |,
w1 (1) = 0.
Set v1 (x) = K[φM (1) − φM (x)], where φM is defined in (3.18) and K = max {kf k∞ + kck∞ |γ1 |, [|γ0 | + |γ1 |]/φM (1)} . Then v1 ≥ 0 and v10 ≤ 0, so by (3.19) we get Lv1 (x) = −D∗δ v1 (x) + M v10 (x) + [b(x) − M ]v10 (x) + c(x)v1 (x) ≥ K,
v1 (0) − α0 v10 (0) = K[φM (1) − φM (0) + α0 φ0M (0)] = KφM (1),
v1 (1) = 0.
It follows from the definition of K and Theorem 3.1 that v1 ± w1 ≥ 0, i.e., |w1 (x)| ≤ v1 (x) for all x. Hence |u(1) − u(x)| v1 (x) |u0 (1)| = lim ≤ lim = Kφ0M (1). (3.22) − − 1−x x→1 x→1 1 − x
But 0 < φ0M (1) < σ0 (M, δ) by Lemma 3.2, and φM (1) = Eδ−1,δ+1 (M ) ≥ 1/2 on taking the first term in the Mittag-Leffler series. Combining these inequalities with (3.22) and the definition of K completes the proof. If instead of 0 ≤ M ≤ 1 one has the stronger hypothesis that kbk∞ ≤ 1, then a bound similar to that for the case α1 = 0 in Theorem 3.4(iii) can be derived quickly using Lemma 3.1, Theorem 3.4(i)(ii) and (2.27). For 0 ≤ M < 1, Theorem 3.4(i) bounds kuk∞ independently of δ. Theorem 3.4(ii) shows that for 0 ≤ M ≤ 1, when δ is near 1 there is no boundary layer in u at x = 0. The situation in Theorem 3.4(iii) is more complicated; to clarify it, we give now a simpler but slightly less sharp corollary of this part of the theorem. Corollary 3.1. Assume that 0 ≤ M ≤ 1. Then for some constant C we have ( C/(1 − M ) if 0 ≤ M < 1, 0 |u (1)| ≤ C/(δ − 1) if M = 1. Proof. Substitute the bound of Theorem 3.4(i) into Theorem 3.4(iii) and recall the definitions of σ1 (M, δ) and σ0 (M, δ).
17
Theorem 3.4(ii) shows that as δ → 1+ there is no boundary layer at x = 0 when 0 ≤ M ≤ 1. Corollary 3.1 shows that as δ → 1+ , no boundary layer is present at x = 1 when 0 ≤ M < 1, but when M = 1 a boundary layer is possible there and the discussion in Section 2.2.2 shows that then the estimate |u0 (1)| ≤ C/(δ − 1) of Corollary 3.1 is sharp in general.
In Theorem 3.4 and Corollary 3.1 we state no result for the case M > 1, despite the fact that our arguments can be extended to this case, because the bound provided by Lemma 3.3 is then excessively large. The next result gives a satisfactory bound on kuk∞ when M > 1 under the extra hypothesis that b > 0 on [0, 1]. In the definition of B4 below, the expression multiplying kf k∞ imitates closely the coefficient of f that we reported in (2.16). Theorem 3.5. Assume that M > 1 and b ≥ b > 0 for some constant b. For 0 ≤ x ≤ 1, set B4 (x) = max{|γ0 |, |γ1 |} +
kf k∞ {α0 [φM (1) − φM (x)] + xφM (1) − φM (x) + α1 [(α0 + x)φ0M (1) − φM (x)]} . 1 + α0 + bφM (1) + α1 [1 + bφ0M (1)] (3.23)
Then |u(x)| ≤ B4 (x) for x ∈ [0, 1] and for some constant C one has (i) (ii)
kuk∞ ≤ kB4 k∞ ≤ C min{α0 , φM (1)},
(iii) (iv)
|u0 (0)| ≤ C,
|u0 (1)| ≤ C min{α0 , φM (1)} if α1 > 0,
|u0 (1)| ≤ CM 1/(δ−1) min{α0 , φM (1)} if α1 = 0.
Proof. The expression xφM (1) − φM (x) vanishes at x = 0, 1 and its second-order derivative is negative by (2.4), so the expression is positive on (0, 1). It is easy to see that φ0M (1) > φM (1) so we also have xφ0M (1) − φM (x) > 0 on (0, 1). Consequently B4 (x) ≥ 0 for 0 ≤ x ≤ 1. Now
B4 (0) − α0 B40 (0) = max{|γ0 |, |γ1 |} ≥ |γ0 |
and B4 (1) + α1 B40 (1) = max{|γ0 |, |γ1 |} ≥ |γ1 |.
For 0 < x < 1, by (3.19) one has LB4 (x) =c(x)B4 (x) +
kf k∞ {α0 + b(x)φM (1) + 1 + α1 [b(x)φ0M (1) + 1] + (M − b)φ0M (x)(1 + α0 + α1 )} 1 + α0 + bφM (1) + α1 [1 + bφ0M (1)]
≥kf k∞ . Thus B4 is a barrier function for ±u by Theorem 3.1. Hence, recalling that α0 ≥ 1/(δ − 1) > 1, for some C we get kuk∞ ≤ kB4 k∞
kf k∞ {α0 φM (1) + φM (1) + α1 (α0 + 1)φ0M (1)} ≤ max{|γ0 |, |γ1 |} + 1 + α0 + bφM (1) + α1 [1 + bφ0M (1)] α0 φM (1) + α0 α1 φ0M (1) ≤C 1+ α0 + φM (1) + α1 φ0M (1) ≤ C min{α0 , φM (1)}. 18
This proves (i). For (ii) and (iii), use the result of (i) and the boundary conditions (1.2b). To prove (iv), define the function u ˜(x) = u(x) − (1 − x)
γ0 − γ1 − γ1 , 1 + α0
(3.24)
which is the solution of the problem L˜ u(x) = f˜(x) for x ∈ (0, 1),
u ˜(0) − α0 u ˜0 (0) = 0,
u ˜(1) = 0,
with
γ0 − γ1 f˜(x) = f (x) + [b(x) − (1 − x)c(x)] − γ1 c(x). 1 + α0 Applying the barrier function B4 of (i) above to u ˜ yields |˜ u(x)| ≤
kf˜k∞ {α0 [φM (1) − φM (x)] + xφM (1) − φM (x)} , 1 + α0 + bφM (1)
for x ∈ [0, 1].
Hence |˜ u0 (1)| = lim
x→1−
kf˜k∞ [(α0 + 1)φ0M (1) + φM (1)] |˜ u(x)| |˜ u(1) − u ˜(x)| = lim ≤ . 1−x 1 + α0 + bφM (1) x→1− 1 − x
Recalling the identity (2.27), we have
φM (1) = Eδ−1,δ+1 (M ) = [Eδ−1,2 (M ) − 1]/M, φ0M (1) = Eδ−1,δ (M ) = [Eδ−1,1 (M ) − 1]/M
and (2.21) implies that Eδ−1,1 (M )/Eδ−1,2 (M ) ≈ M 1/(δ−1) as δ → 1+ . Consequently |˜ u0 (1)| ≤
Cα0 φ0M (1) Cα0 M 1/(δ−1) φM (1) ≤ ≤ CM 1/(δ−1) min{α0 , φM (1)} . α0 + φM (1) α0 + φM (1)
Part (iv) of the theorem follows from this estimate and (3.24). Thus, in the case M > 1 with b > 0, Theorem 3.5 (ii) shows that there is no boundary layer at x = 0 when δ is near 1, even though kuk∞ may be unbounded as δ → 1+ (see Section 2.2.1).
The accurate approximations (2.22) and (2.23) show that Theorem 3.5 (iii)(iv) is sharp for constant b > 1.
The next example shows numerically that the bounds of Theorem 3.5 (i)(iii) are sharp for a variable-coefficient problem when α0 = 1/(δ − 1) < φM (1). Example 3.2. Consider the test problem −D∗δ u(x) + (x + 0.2)u0 (x) = 2(3 + x) for x ∈ (0, 1), 1 u0 (0) = 0.4, u(1) + u0 (1) = 1.7. u(0) − δ−1
(3.25a) (3.25b)
The exact solution of this problem is unknown. Observe that M = 1.2 so φM (1) 1/(δ−1) = α0 by an inequality similar to (2.19). To check the bounds of Theorem 3.5 (i)(iii), we use the finite difference scheme of Stynes and Gracia (2014) to compute an approximate numerical solution N {uj }N j=0 of (3.25) on a uniform mesh {xj := j/N }j=0 with N = 2048 for various values of δ close to 1, then evaluate (δ − 1) max0≤j≤N |uj | and (δ − 1)|uN − uN −1 |/(xN − xN −1 ) to approximate (δ − 1)kuk∞ and (δ − 1)|u0 (1)|.
As each row of the table is approximately constant as δ → 1+ , the bounds of Theorem 3.5 (i)(iii) are sharp for this example. 19
Table 3.1: Verifying sharpness of Theorem 3.5 (i)(iii) (δ − 1) max0≤j≤N |uj | (δ − 1)|uN − uN −1 |/(xN − xN −1 )
4
δ = 1.1
δ = 1.01
δ = 1.001
δ = 1.0001
δ = 1.00001
δ = 1.000001
δ = 1.0000001
5.8374 5.0045
7.6438 7.5954
7.6043 7.5795
7.6004 7.5777
7.6000 7.5775
7.6000 7.5775
7.6000 7.5775
Conclusions
We considered a two-point boundary value problem whose leading term is a Caputo fractional derivative of order δ with 1 < δ < 2. The dependence of the solution on the parameter δ has not previously been investigated analytically, despite a growing interest in the research literature in the numerical solution of problems with variable fractional derivatives. By considering first the special case of a constant-coefficient operator, for which the solution can be determined explicitly, we showed that when δ is near 1, the solution of the boundary value problem may exhibit a boundary layer at the endpoint x = 1 of the domain. Moving on to the general case of a variable-coefficient differential operator, we then determined conditions on the data of the problem under which boundary layers at each endpoint (x = 0, 1) cannot occur. This analysis showed that a crucial parameter in the presence or absence of a boundary layer at x = 1 is the quantity M := maxx∈[0,1] b(x), where b is the coefficient of the first-order term in the differential operator. In all cases considered, we showed that |u0 (0)| ≤ C, i.e., no boundary layer in u appears at x = 0 when δ is near 1. The only data regime where this bound is not guaranteed by our theory (Theorems 3.2 and 3.5) is when min[0,1] c(x) = 0 and M > 1 without the additional property that b > 0 on [0, 1], but our numerical experience (using the finite difference method of Stynes and Gracia (2014)) is that in this case also no boundary layer appears at x = 0. At x = 1, our theory proves rigorously that when δ is near 1, no boundary layer appears in u if M < 1, but one can have such layers when M ≥ 1. This agrees with our numerical experience. This analysis of the solution u of (1.2) leads naturally to the question: can one construct a numerical method that will yield an accurate approximation of u when it has a boundary layer at x = 1?
References Mohammed Al-Refai. Basic results on nonlinear eigenvalue problems of fractional order. Electron. J. Differ. Equ., 2012(191):1–12, 2012. Hermann Brunner. Collocation methods for Volterra integral and related functional differential equations, volume 15 of Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, Cambridge, 2004. ISBN 0-521-80615-1. doi: 10.1017/CBO9780511543234. Chang-Ming Chen, F. Liu, K. Burrage, and Y. Chen. Numerical methods of the variable-order Rayleigh-Stokes problem for a heated generalized second grade fluid with fractional derivative. IMA J. Appl. Math., 78(5):924–944, 2013. ISSN 0272-4960. doi: 10.1093/imamat/hxr079. Kai Diethelm. The analysis of fractional differential equations, volume 2004 of Lecture Notes in Mathematics. Springer-Verlag, Berlin, 2010. ISBN 978-3-642-14573-5. An application-oriented exposition using differential operators of Caputo type. 20
Jos´e Luis Gracia and Martin Stynes. Upwind and central difference approximation of convection in Caputo fractional derivative two-point boundary value problems. J. Comput. Appl. Math., 273:103–115, 2015. Bangti Jin, Raytcho Lazarov, and Joseph Pasciak. Variational formulation of problems involving fractional order differential operators. ArXiv:1307.4795, 2013. Natalia Kopteva and Martin Stynes. An efficient collocation method for a Caputo two-point boundary value problem. 2014. (Submitted for publication). J. Tenreiro Machado, Virginia Kiryakova, and Francesco Mainardi. Recent history of fractional calculus. Commun. Nonlinear Sci. Numer. Simul., 16(3):1140–1153, 2011. ISSN 1007-5704. A. M. Mathai, R. K. Saxena, and H. J. Haubold. A certain class of Laplace transforms with applications to reaction and reaction-diffusion equations. Astrophys. Space Sci., 305(3):283– 288, 2006. doi: 10.1007/s10509-006-9188-7. Arvet Pedas and Enn Tamme. Piecewise polynomial collocation for linear boundary 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. Igor Podlubny. Fractional differential equations, volume 198 of Mathematics in Science and Engineering. Academic Press Inc., San Diego, CA, 1999. ISBN 0-12-558840-2. An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. John Paul Roop. Numerical approximation of a one-dimensional space fractional advection-dispersion equation with boundary layer. Comput. Math. Appl., 56(7): 1808–1819, 2008. ISSN 0898-1221. doi: 10.1016/j.camwa.2008.04.025. URL http://dx.doi.org/10.1016/j.camwa.2008.04.025. Hans-G¨org Roos, Martin Stynes, and Lutz Tobiska. Robust numerical methods for singularly perturbed differential equations, volume 24 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 2008. ISBN 978-3-540-34466-7. Convection-diffusionreaction and flow problems. Stefan Samko. Fractional integration and differentiation of variable order: an overview. Nonlinear Dynam., 71(4):653–662, 2013. ISSN 0924-090X. doi: 10.1007/s11071-012-0485-0. Martin Stynes and Jos´e Luis Gracia. A finite difference method for a two-point boundary value problem with a Caputo fractional derivative. IMA J. Numer. Anal., 2014. doi: 10.1093/imanum/dru011.
21
