- Email: [email protected]
Preprocessing schemes for fractional-derivative problems to improve their convergence rates
Applied Mathematics Letters 74 (2017) 187–192
Contents lists available at ScienceDirect
Applied Mathematics Letters www.elsevier.com/locate/aml
Preprocessing schemes for fractional-derivative problems to improve their convergence rates Martin Stynesa , José Luis Graciab, * a
Applied and Computational Mathematics Division, Beijing Computational Science Research Center, Haidian District, Beijing 100193, China b Department of Applied Mathematics, Torres Quevedo Building, Campus Rio Ebro, University of Zaragoza, 50018 Zaragoza, Spain
article
info
abstract
Article history: Received 1 April 2017 Received in revised form 25 May 2017 Accepted 25 May 2017 Available online 15 June 2017
A simple and inexpensive preprocessing of an initial–boundary value problem with a Caputo time derivative is shown theoretically and numerically to yield an enhanced convergence rate for the L1 scheme. The same preprocessing can also be used with other methods for time-dependent problems. © 2017 Elsevier Ltd. All rights reserved.
Keywords: Caputo derivative Weak singularity Preprocessing of solution L1 scheme
1. Introduction Consider the initial–boundary value problem Lu(x, t) := Dtα u − puxx + c(x)u = f (x, t)
(1a)
for (x, t) ∈ Q := (0, l) × (0, T ], with u(0, t) = u(l, t) = 0 u(x, 0) = ϕ(x)
for t ∈ (0, T ],
for x ∈ [0, l],
(1b) (1c)
¯ where Q ¯ := [0, l] × [0, T ], and where 0 < α < 1, p is a positive constant, c ∈ C[0, l] with c ≥ 0, f ∈ C(Q) α ϕ ∈ C[0, l]. In (1a), Dt u denotes a Caputo fractional derivative of order α with 0 < α < 1, which is defined by ∫ t 1 ∂u(x, s) Dtα u(x, t) = (t − s)−α ds for (x, t) ∈ Q. Γ (1 − α) s=0 ∂s
*
Corresponding author. E-mail addresses: [email protected] (M. Stynes), [email protected] (J.L. Gracia).
http://dx.doi.org/10.1016/j.aml.2017.05.016 0893-9659/© 2017 Elsevier Ltd. All rights reserved.
M. Stynes, J.L. Gracia / Applied Mathematics Letters 74 (2017) 187–192
188
It is known [1,2] that under reasonable hypotheses on its data, the problem (1) has a unique solution which typically exhibits a weak singularity at t = 0; roughly speaking, u(x, t) behaves like the function tα near t = 0. Consequently the temporal derivatives ut and utt are unbounded at t = 0, which presents challenges for numerical methods for (1), and for their analysis. Many different numerical methods have been suggested for (1); see [3] for a list of these. In the present paper we show that a simple preprocessing of the solution u of (1) will improve the rate of convergence of the computed solution when the well-known L1 scheme is used to approximate the fractional derivative Dtα u. While this scheme is the only one that we discuss here, nevertheless the preprocessing idea discussed in this paper can clearly be used to improve the rate of convergence of many methods for (1) and for other problems such as the fractional wave-diffusion equation (for which 1 < α < 2 in (1a)). Notation. By C we denote a generic constant that depends on the data of the problem (1), i.e., C = C(α, p, c, f, ϕ, l, T ). It can take different values in different places but is independent of any mesh used. 2. The basic idea As in [1], we set D(Lγ ) =
⎧ ⎨
ω ∈ L2 (0, l) :
⎩
∞ ∑ i=1
⎫ ⏐2 ⏐ ⎬ ⏐ ω(s)ψi (s) ds⏐ < ∞ ⏐ s=0 ⏐ ⎭ ⏐∫ ⏐
l
⏐ λ2γ i ⏐
for γ ≥ 0,
where {(λi , ψi ) : i = 1, 2, . . .} are the eigenvalues and normalised eigenfunctions of the Sturm–Liouville two-point boundary value problem Lψi := −pψi′′ + cψi = λi ψi on (0, l),
ψi (0) = ψi (l) = 0.
(2)
The norm associated with the space D(Lγ ) is denoted by ∥ · ∥Lγ . As in [2, Theorem 2.1], assume that ϕ ∈ D(L5/2 ), f (·, t) ∈ D(L5/2 ), ft (·, t) and ftt (·, t) are in D(L1/2 ) for each t ∈ (0, T ] with ∥f (·, t)∥L5/2 + ∥ft (·, t)∥L1/2 + tρ ∥ftt (·, t)∥L1/2 ≤ C1 for all t ∈ (0, T ] and some constant ρ < 1, where C1 is a constant independent of t. These hypotheses imply in particular that 0 = ϕ(0) = ϕ′′ (0) = ϕ(l) = ϕ′′ (l) = f (0, t) = f (l, t) for 0 ≤ t ≤ T.
(3)
To solve (1) numerically, use the rectangular grid (xm , tn ) where the uniform spatial mesh is xm = mh for m = 0, 1, . . . , M , where h = l/M for some positive integer M , while the temporal mesh is graded: tn = T (n/N )r for n = 0, 1, . . . , N and the mesh grading r ≥ 1 is chosen by the user. As in [2], we approximate Dtα u by the L1 scheme and −puxx + c(x)u by a standard O(h2 ) finite difference approximation; we shall refer to this time and space discretisation simply as “the L1 scheme”. It is shown in [2, Theorem 5.2] that the solution {unm } of (1) of this scheme satisfies ( ) max |u(xm , tn ) − unm | ≤ C h2 + N − min{2−α, rα} (4) ¯ (xm ,tn )∈Q
for some constant C. ¯ and it In [4, Lemma 1] it is assumed that ϕ ∈ C 4 [0, l], (3) holds true, c ∈ C 2 [0, l] and f, fx , fxx ∈ C(Q), is shown that the solution u of (1) can be decomposed as ¯ u(x, t) = z(x)tα + ϕ(x) + v(x, t) for (x, t) ∈ Q,
(5)
M. Stynes, J.L. Gracia / Applied Mathematics Letters 74 (2017) 187–192
189
where z(x) =
1 [f (x, 0) + pϕ′′ (x) − c(x)ϕ(x)] Γ (α + 1)
(6)
can be computed explicitly from the given data, and the remainder term v is the solution of the initial– boundary value problem ( ) (Dtα v − pvxx + cv) (x, t) = f (x, t) − L z(x)tα + ϕ(x) = f (x, t) − f (x, 0) + pz ′′ (x)tα − c(x)z(x)tα
(7a)
in Q, with (observe that (3) implies that z(0) = z(l) = 0) v(0, t) = v(l, t) = 0 for t ∈ (0, T ], and v(x, 0) = 0
for x ∈ [0, l].
(7b)
The correctness of (5)–(7) can be verified easily by a direct calculation. Now in [4, Theorem 1] it is proved, under the extra regularity assumptions fttt (·, t) ∈ D(L1/2 ) and ∥ftt (·, t)∥L1/2 + tρ˜∥fttt (·, t)∥L1/2 ≤ C2 , where C2 is a constant independent of t and 0 < ρ˜ < 1, that v satisfies ⏐ ⏐ k ⏐∂ v ⏐ ⏐ ≤ C for k = 0, 1, 2, 3, 4, ⏐ (8a) (x, t) ⏐ ⏐ ∂xk ⏐ ℓ ⏐ ⏐∂ v ⏐ ⏐ ⏐ ≤ Ct2α−ℓ for ℓ = 1, 2, (8b) (x, t) ⏐ ∂tℓ ⏐ for all (x, t) ∈ [0, l] × (0, T ] and some constant C. That is, v is smoother than u (for which one has similar bounds but with 2α replaced by α in (8b); see [2, Theorem 2.1]). Consequently when the L1 scheme is used to solve (7) one obtains, instead of (4), the improved bound ( ) n max |v(xm , tn ) − vm | ≤ C h2 + N − min{2−α, 2rα} (9) ¯ (xm ,tn )∈Q
n is the computed solution. This observation is the inspiration for our preprocessing idea. where vm
2.1. Preprocessed L1 method n Compute z from (6), then solve (7) numerically using the L1 scheme to get {vm }. Our approximation u ˜nm of u(xm , tn ) for each m and n is then—based on (5)—defined by n u ˜nm := z(xm )tα n + ϕ(xm ) + vm .
As z and ϕ can be computed exactly, one has n u(xm , tn ) − u ˜nm = v(xm , tn ) − vm for all m and n,
so from (9) we deduce that max
( ) |u(xm , tn ) − u ˜nm | ≤ C h2 + N − min{2−α, 2rα} .
¯ (xm ,tn )∈Q
(10)
Comparing this bound with (4) shows that: (i) when the mesh is uniform (r = 1), the preprocessed scheme is much more accurate (to confirm this, see the r = 1 rows in Tables 1 and 2); (ii) the optimal mesh for the preprocessed scheme is much less severely graded than the optimal mesh for the L1 scheme (which reduces the risk of rounding error—see Remark 1). As the preprocessed scheme needs negligible extra computation, it is a clear improvement on the L1 scheme.
3.704E−2 0.099
3.601E−4 1.655
r=1
r = (2 − α)/α
N = M = 64
1.144E−4 1.672
3.458E−2 0.109
N = M = 128
3.589E−5 1.689
3.207E−2 0.118
N = M = 256
Table 1 L1 scheme: Maximum errors and orders of convergence for α = 0.2. N = M = 512
1.113E−5 1.704
2.955E−2 0.127
N = M = 1024
3.418E−6 1.717
2.706E−2 0.135
1.040E−6 1.729
2.465E−2 0.142
N = M = 2048
3.138E−7
2.233E−2
N = M = 4096
190 M. Stynes, J.L. Gracia / Applied Mathematics Letters 74 (2017) 187–192
1.184E−2 0.301
1.397E−4 1.632
r=1
r = (2 − α)/(2α)
N = M = 64
4.509E−5 1.658
9.610E−3 0.312
N = M = 128
1.429E−5 1.677
7.741E−3 0.322
N = M = 256
Table 2 Preprocessed L1 scheme: Maximum errors and orders of convergence for α = 0.2. N = M = 512
4.469E−6 1.693
6.194E−3 0.330
N = M = 1024
1.382E−6 1.707
4.928E−3 0.338
4.235E−7 1.719
3.899E−3 0.345
N = M = 2048
1.287E−7
3.070E−3
N = M = 4096
M. Stynes, J.L. Gracia / Applied Mathematics Letters 74 (2017) 187–192 191
192
M. Stynes, J.L. Gracia / Applied Mathematics Letters 74 (2017) 187–192
−r Remark 1. The finest mesh interval in our graded mesh tn = T (n/N )r has diameter O(N ). The change ( ) from rα in (4) to 2rα in (10) implies that to attain the optimal rate of convergence O h2 + N −(2−α) , the √ diameters mu , mv of the finest mesh intervals needed for u, v respectively will satisfy mv ≈ mu . Thus the optimal v-mesh is much coarser than the optimal u-mesh near t = 0, which reduces the risk of rounding errors caused by very fine meshes.
Remark 2. If a convective term ∫ x q(x)(∂u/∂x) also appears in (1a), it can be removed by the change of variable y(x, t) := u(x, t) exp[−( 0 q(s) ds)/(2p)], which yields a differential equation for y of the form (1a). 3. Numerical example Consider the test problem [2, Example 6.1] Dtα u − uxx = f (x, t) for (x, t) ∈ (0, π) × (0, 1],
(11)
with initial condition u(x, 0) = sin x for 0 < x < π, and boundary conditions u(0, t) = u(π, t) = 0 for 0 ≤ t[ ≤ 1. The function f in (11) is chosen such that the exact solution of the problem is ] u(x, t) = Eα (−tα ) + t3 sin x, where Eα (·) is the classical Mittag-Leffler function. This solution u has exactly the regularity at t = 0 that typical solutions of (1) enjoy; see [2]. The maximum error in each computed solution and the corresponding rate of convergence are denoted by ( ) eM,N := max |u(xm , tn ) − unm |, rateM,N := log2 eM,N /e2M,2N . ¯ (xm ,tn )∈Q
We take α = 0.2. This value is chosen because it gives slow convergence of the L1 solution. Table 1 gives the maximum errors and the rates of convergence for the L1 scheme of [2], where we used the mesh gradings r = 1 (uniform mesh) and r = (2 − α)/α = 9 (optimal choice). The numerical results for the preprocessed L1 scheme are given in Table 2 with r = 1 and r = (2 − α)/(2α) = 4.5 (optimal choice). The numerical results in both tables are in agreement with our theoretical convergence results. If one approximates the problem using a uniform mesh, the preprocessed L1 method provides more accurate approximations to the solution and better orders of convergence than the L1 scheme. ( If the meshes ) are graded in an optimal way, then although both methods have theoretical accuracy O h2 + N −(2−α) , the preprocessed method gives smaller errors because its grid points are not excessively condensed near t = 0. Acknowledgments The research of Martin Stynes was supported in part by the National Natural Science Foundation of China under grants 91430216 and NSAF-U1530401. The research of Jos´e Luis Gracia was partly supported by the Institute of Mathematics and Applications (IUMA), the project MTM2016-75139-R and the Diputaci´on General de Arag´ on. References
[1] K. Sakamoto, M. Yamamoto, Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems, J. Math. Anal. Appl. 382 (1) (2011) 426–447. http://dx.doi.org/10.1016/j.jmaa.2011.04.058. [2] M. Stynes, E. O’Riordan, J.L. Gracia, Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation, SIAM J. Numer. Anal. 55 (2) (2017) 1057–1079. http://dx.doi.org/10.1137/16M1082329. [3] B. Jin, R. Lazarov, Z. Zhou, Two fully discrete schemes for fractional diffusion and diffusion-wave equations with nonsmooth data, SIAM J. Sci. Comput. 38 (1) (2016) A146–A170. http://dx.doi.org/10.1137/140979563. [4] J.L. Gracia, E. O’Riordan, M. Stynes, A fitted scheme for a Caputo initial–boundary value problem, in preparation.
Contents lists available at ScienceDirect
Applied Mathematics Letters www.elsevier.com/locate/aml
Preprocessing schemes for fractional-derivative problems to improve their convergence rates Martin Stynesa , José Luis Graciab, * a
Applied and Computational Mathematics Division, Beijing Computational Science Research Center, Haidian District, Beijing 100193, China b Department of Applied Mathematics, Torres Quevedo Building, Campus Rio Ebro, University of Zaragoza, 50018 Zaragoza, Spain
article
info
abstract
Article history: Received 1 April 2017 Received in revised form 25 May 2017 Accepted 25 May 2017 Available online 15 June 2017
A simple and inexpensive preprocessing of an initial–boundary value problem with a Caputo time derivative is shown theoretically and numerically to yield an enhanced convergence rate for the L1 scheme. The same preprocessing can also be used with other methods for time-dependent problems. © 2017 Elsevier Ltd. All rights reserved.
Keywords: Caputo derivative Weak singularity Preprocessing of solution L1 scheme
1. Introduction Consider the initial–boundary value problem Lu(x, t) := Dtα u − puxx + c(x)u = f (x, t)
(1a)
for (x, t) ∈ Q := (0, l) × (0, T ], with u(0, t) = u(l, t) = 0 u(x, 0) = ϕ(x)
for t ∈ (0, T ],
for x ∈ [0, l],
(1b) (1c)
¯ where Q ¯ := [0, l] × [0, T ], and where 0 < α < 1, p is a positive constant, c ∈ C[0, l] with c ≥ 0, f ∈ C(Q) α ϕ ∈ C[0, l]. In (1a), Dt u denotes a Caputo fractional derivative of order α with 0 < α < 1, which is defined by ∫ t 1 ∂u(x, s) Dtα u(x, t) = (t − s)−α ds for (x, t) ∈ Q. Γ (1 − α) s=0 ∂s
*
Corresponding author. E-mail addresses: [email protected] (M. Stynes), [email protected] (J.L. Gracia).
http://dx.doi.org/10.1016/j.aml.2017.05.016 0893-9659/© 2017 Elsevier Ltd. All rights reserved.
M. Stynes, J.L. Gracia / Applied Mathematics Letters 74 (2017) 187–192
188
It is known [1,2] that under reasonable hypotheses on its data, the problem (1) has a unique solution which typically exhibits a weak singularity at t = 0; roughly speaking, u(x, t) behaves like the function tα near t = 0. Consequently the temporal derivatives ut and utt are unbounded at t = 0, which presents challenges for numerical methods for (1), and for their analysis. Many different numerical methods have been suggested for (1); see [3] for a list of these. In the present paper we show that a simple preprocessing of the solution u of (1) will improve the rate of convergence of the computed solution when the well-known L1 scheme is used to approximate the fractional derivative Dtα u. While this scheme is the only one that we discuss here, nevertheless the preprocessing idea discussed in this paper can clearly be used to improve the rate of convergence of many methods for (1) and for other problems such as the fractional wave-diffusion equation (for which 1 < α < 2 in (1a)). Notation. By C we denote a generic constant that depends on the data of the problem (1), i.e., C = C(α, p, c, f, ϕ, l, T ). It can take different values in different places but is independent of any mesh used. 2. The basic idea As in [1], we set D(Lγ ) =
⎧ ⎨
ω ∈ L2 (0, l) :
⎩
∞ ∑ i=1
⎫ ⏐2 ⏐ ⎬ ⏐ ω(s)ψi (s) ds⏐ < ∞ ⏐ s=0 ⏐ ⎭ ⏐∫ ⏐
l
⏐ λ2γ i ⏐
for γ ≥ 0,
where {(λi , ψi ) : i = 1, 2, . . .} are the eigenvalues and normalised eigenfunctions of the Sturm–Liouville two-point boundary value problem Lψi := −pψi′′ + cψi = λi ψi on (0, l),
ψi (0) = ψi (l) = 0.
(2)
The norm associated with the space D(Lγ ) is denoted by ∥ · ∥Lγ . As in [2, Theorem 2.1], assume that ϕ ∈ D(L5/2 ), f (·, t) ∈ D(L5/2 ), ft (·, t) and ftt (·, t) are in D(L1/2 ) for each t ∈ (0, T ] with ∥f (·, t)∥L5/2 + ∥ft (·, t)∥L1/2 + tρ ∥ftt (·, t)∥L1/2 ≤ C1 for all t ∈ (0, T ] and some constant ρ < 1, where C1 is a constant independent of t. These hypotheses imply in particular that 0 = ϕ(0) = ϕ′′ (0) = ϕ(l) = ϕ′′ (l) = f (0, t) = f (l, t) for 0 ≤ t ≤ T.
(3)
To solve (1) numerically, use the rectangular grid (xm , tn ) where the uniform spatial mesh is xm = mh for m = 0, 1, . . . , M , where h = l/M for some positive integer M , while the temporal mesh is graded: tn = T (n/N )r for n = 0, 1, . . . , N and the mesh grading r ≥ 1 is chosen by the user. As in [2], we approximate Dtα u by the L1 scheme and −puxx + c(x)u by a standard O(h2 ) finite difference approximation; we shall refer to this time and space discretisation simply as “the L1 scheme”. It is shown in [2, Theorem 5.2] that the solution {unm } of (1) of this scheme satisfies ( ) max |u(xm , tn ) − unm | ≤ C h2 + N − min{2−α, rα} (4) ¯ (xm ,tn )∈Q
for some constant C. ¯ and it In [4, Lemma 1] it is assumed that ϕ ∈ C 4 [0, l], (3) holds true, c ∈ C 2 [0, l] and f, fx , fxx ∈ C(Q), is shown that the solution u of (1) can be decomposed as ¯ u(x, t) = z(x)tα + ϕ(x) + v(x, t) for (x, t) ∈ Q,
(5)
M. Stynes, J.L. Gracia / Applied Mathematics Letters 74 (2017) 187–192
189
where z(x) =
1 [f (x, 0) + pϕ′′ (x) − c(x)ϕ(x)] Γ (α + 1)
(6)
can be computed explicitly from the given data, and the remainder term v is the solution of the initial– boundary value problem ( ) (Dtα v − pvxx + cv) (x, t) = f (x, t) − L z(x)tα + ϕ(x) = f (x, t) − f (x, 0) + pz ′′ (x)tα − c(x)z(x)tα
(7a)
in Q, with (observe that (3) implies that z(0) = z(l) = 0) v(0, t) = v(l, t) = 0 for t ∈ (0, T ], and v(x, 0) = 0
for x ∈ [0, l].
(7b)
The correctness of (5)–(7) can be verified easily by a direct calculation. Now in [4, Theorem 1] it is proved, under the extra regularity assumptions fttt (·, t) ∈ D(L1/2 ) and ∥ftt (·, t)∥L1/2 + tρ˜∥fttt (·, t)∥L1/2 ≤ C2 , where C2 is a constant independent of t and 0 < ρ˜ < 1, that v satisfies ⏐ ⏐ k ⏐∂ v ⏐ ⏐ ≤ C for k = 0, 1, 2, 3, 4, ⏐ (8a) (x, t) ⏐ ⏐ ∂xk ⏐ ℓ ⏐ ⏐∂ v ⏐ ⏐ ⏐ ≤ Ct2α−ℓ for ℓ = 1, 2, (8b) (x, t) ⏐ ∂tℓ ⏐ for all (x, t) ∈ [0, l] × (0, T ] and some constant C. That is, v is smoother than u (for which one has similar bounds but with 2α replaced by α in (8b); see [2, Theorem 2.1]). Consequently when the L1 scheme is used to solve (7) one obtains, instead of (4), the improved bound ( ) n max |v(xm , tn ) − vm | ≤ C h2 + N − min{2−α, 2rα} (9) ¯ (xm ,tn )∈Q
n is the computed solution. This observation is the inspiration for our preprocessing idea. where vm
2.1. Preprocessed L1 method n Compute z from (6), then solve (7) numerically using the L1 scheme to get {vm }. Our approximation u ˜nm of u(xm , tn ) for each m and n is then—based on (5)—defined by n u ˜nm := z(xm )tα n + ϕ(xm ) + vm .
As z and ϕ can be computed exactly, one has n u(xm , tn ) − u ˜nm = v(xm , tn ) − vm for all m and n,
so from (9) we deduce that max
( ) |u(xm , tn ) − u ˜nm | ≤ C h2 + N − min{2−α, 2rα} .
¯ (xm ,tn )∈Q
(10)
Comparing this bound with (4) shows that: (i) when the mesh is uniform (r = 1), the preprocessed scheme is much more accurate (to confirm this, see the r = 1 rows in Tables 1 and 2); (ii) the optimal mesh for the preprocessed scheme is much less severely graded than the optimal mesh for the L1 scheme (which reduces the risk of rounding error—see Remark 1). As the preprocessed scheme needs negligible extra computation, it is a clear improvement on the L1 scheme.
3.704E−2 0.099
3.601E−4 1.655
r=1
r = (2 − α)/α
N = M = 64
1.144E−4 1.672
3.458E−2 0.109
N = M = 128
3.589E−5 1.689
3.207E−2 0.118
N = M = 256
Table 1 L1 scheme: Maximum errors and orders of convergence for α = 0.2. N = M = 512
1.113E−5 1.704
2.955E−2 0.127
N = M = 1024
3.418E−6 1.717
2.706E−2 0.135
1.040E−6 1.729
2.465E−2 0.142
N = M = 2048
3.138E−7
2.233E−2
N = M = 4096
190 M. Stynes, J.L. Gracia / Applied Mathematics Letters 74 (2017) 187–192
1.184E−2 0.301
1.397E−4 1.632
r=1
r = (2 − α)/(2α)
N = M = 64
4.509E−5 1.658
9.610E−3 0.312
N = M = 128
1.429E−5 1.677
7.741E−3 0.322
N = M = 256
Table 2 Preprocessed L1 scheme: Maximum errors and orders of convergence for α = 0.2. N = M = 512
4.469E−6 1.693
6.194E−3 0.330
N = M = 1024
1.382E−6 1.707
4.928E−3 0.338
4.235E−7 1.719
3.899E−3 0.345
N = M = 2048
1.287E−7
3.070E−3
N = M = 4096
M. Stynes, J.L. Gracia / Applied Mathematics Letters 74 (2017) 187–192 191
192
M. Stynes, J.L. Gracia / Applied Mathematics Letters 74 (2017) 187–192
−r Remark 1. The finest mesh interval in our graded mesh tn = T (n/N )r has diameter O(N ). The change ( ) from rα in (4) to 2rα in (10) implies that to attain the optimal rate of convergence O h2 + N −(2−α) , the √ diameters mu , mv of the finest mesh intervals needed for u, v respectively will satisfy mv ≈ mu . Thus the optimal v-mesh is much coarser than the optimal u-mesh near t = 0, which reduces the risk of rounding errors caused by very fine meshes.
Remark 2. If a convective term ∫ x q(x)(∂u/∂x) also appears in (1a), it can be removed by the change of variable y(x, t) := u(x, t) exp[−( 0 q(s) ds)/(2p)], which yields a differential equation for y of the form (1a). 3. Numerical example Consider the test problem [2, Example 6.1] Dtα u − uxx = f (x, t) for (x, t) ∈ (0, π) × (0, 1],
(11)
with initial condition u(x, 0) = sin x for 0 < x < π, and boundary conditions u(0, t) = u(π, t) = 0 for 0 ≤ t[ ≤ 1. The function f in (11) is chosen such that the exact solution of the problem is ] u(x, t) = Eα (−tα ) + t3 sin x, where Eα (·) is the classical Mittag-Leffler function. This solution u has exactly the regularity at t = 0 that typical solutions of (1) enjoy; see [2]. The maximum error in each computed solution and the corresponding rate of convergence are denoted by ( ) eM,N := max |u(xm , tn ) − unm |, rateM,N := log2 eM,N /e2M,2N . ¯ (xm ,tn )∈Q
We take α = 0.2. This value is chosen because it gives slow convergence of the L1 solution. Table 1 gives the maximum errors and the rates of convergence for the L1 scheme of [2], where we used the mesh gradings r = 1 (uniform mesh) and r = (2 − α)/α = 9 (optimal choice). The numerical results for the preprocessed L1 scheme are given in Table 2 with r = 1 and r = (2 − α)/(2α) = 4.5 (optimal choice). The numerical results in both tables are in agreement with our theoretical convergence results. If one approximates the problem using a uniform mesh, the preprocessed L1 method provides more accurate approximations to the solution and better orders of convergence than the L1 scheme. ( If the meshes ) are graded in an optimal way, then although both methods have theoretical accuracy O h2 + N −(2−α) , the preprocessed method gives smaller errors because its grid points are not excessively condensed near t = 0. Acknowledgments The research of Martin Stynes was supported in part by the National Natural Science Foundation of China under grants 91430216 and NSAF-U1530401. The research of Jos´e Luis Gracia was partly supported by the Institute of Mathematics and Applications (IUMA), the project MTM2016-75139-R and the Diputaci´on General de Arag´ on. References
[1] K. Sakamoto, M. Yamamoto, Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems, J. Math. Anal. Appl. 382 (1) (2011) 426–447. http://dx.doi.org/10.1016/j.jmaa.2011.04.058. [2] M. Stynes, E. O’Riordan, J.L. Gracia, Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation, SIAM J. Numer. Anal. 55 (2) (2017) 1057–1079. http://dx.doi.org/10.1137/16M1082329. [3] B. Jin, R. Lazarov, Z. Zhou, Two fully discrete schemes for fractional diffusion and diffusion-wave equations with nonsmooth data, SIAM J. Sci. Comput. 38 (1) (2016) A146–A170. http://dx.doi.org/10.1137/140979563. [4] J.L. Gracia, E. O’Riordan, M. Stynes, A fitted scheme for a Caputo initial–boundary value problem, in preparation.