- Email: [email protected]
An efficient numerical method for q -fractional differential equations
Journal Pre-proof An efficient numerical method for q-fractional differential equations
Pin Lyu, Seakweng Vong
PII: DOI: Reference:
S0893-9659(19)30482-3 https://doi.org/10.1016/j.aml.2019.106156 AML 106156
To appear in:
Applied Mathematics Letters
Received date : 23 August 2019 Revised date : 18 November 2019 Accepted date : 18 November 2019 Please cite this article as: P. Lyu and S. Vong, An efficient numerical method for q-fractional differential equations, Applied Mathematics Letters (2019), doi: https://doi.org/10.1016/j.aml.2019.106156. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Elsevier Ltd. All rights reserved.
Journal Pre-proof
An efficient numerical method for q-fractional differential equations Pin Lyua , Seakweng Vongb,∗ a
School of Economic Mathematics, Southwestern University of Finance and Economics, Chengdu, China b Department of Mathematics, University of Macau, Macao, China
Abstract In this work, we propose and analyze an efficient numerical method for solving nonlinear Caputo q-fractional differential equations with non-smooth solutions. The numerical method is a linearized one and can be implemented on the graded mesh. For equations with fractional derivatives of order α ∈ (1, 2), the method converges with second-order on uniform mesh. For α ∈ (0, 1), second-order convergence can be obtained when the grading parameter of the nonuniform mesh is set appropriately.
1. Introduction
na lP repr oo f
Keywords: Caputo q-fractional derivative, nonuniform mesh, fractional nonlinear equation
Recently, q-fractional differential equations have attracted lots of interests in their analytical and numerical studies, see also [1–3, 22, 25, 26] and references therein. It is well known that the solutions of fractional (ordinary and partial) differential equations typically exhibit weak singularities near the initial/boundary point. This has become a major concern in the numerical analysis of fractional differential equations [8, 9, 13–19, 23, 24]. Similar weak singularities may also appear in q-fractional differential equations, referring to our brief discussion in the next section. A method to solve the q-fractional differential equations is proposed in a very recent work [26], in which the convergence analysis is absent. Furthermore, the truncation error analysis in [26] requires to assume that the solution is smooth enough on the whole closed domain. In view of this, based on some techniques of our recent work [18], we propose and analyze another efficient numerical method for the q-fractional differential equation: c
Dqα y(t) = f (t, y),
t ∈ (0, T ],
(1.1) dαe−1
subject to the initial condition y(0) = y0 if α ∈ (0, 1); and y(0) = y0 , Dq y(0) = yˆ0 if α ∈ (1, 2). Definitions of q-analogy (0 < |q| < 1) derivatives, integrals and related functions which will be utilized in this paper are listed below:
ur
• Caputo q-fractional derivative [1]:
c
• Riemann-Liouville q-fractional integral [5]:
Jo
• q-Gamma function [5]:
dαe−α
Dqα y(t) := Iq
dαe
Dq y(t). Rt Iqδ g(t) := Γq1(δ) 0 (t − qs)(δ−1) g(s)dq s, δ > 0.
Γq (δ) := (1 − q)1−δ (1 − q)(δ−1) .
• q-analogy of (t − s)(ν) [5]: (t − s)(ν) := tν
Q∞
t−q i s i=0 t−q ν+i s
for a non-positive integer ν.
I The first author is supported by the Fundamental Research Funds for the Central Universities (JBK1901027), the second author is supported by The Science and Technology Development Fund, Macau SAR (File no. 0005/2019/A) and the grants MYRG2018-00047-FST and MYRG2017-00098-FST from University of Macau. ∗ Corresponding author. Email addresses: [email protected] (Pin Lyu), [email protected] (Seakweng Vong)
Preprint submitted to Applied Mathematics Letters
December 2, 2019
Journal Pre-proof
• q-derivative [7] :
Dq g(t) :=
• q-integral [6]: Iq (g) :=
Rb a
g(qt)−g(t) (q−1)t ,
g(s)dq s =
Rb 0
• q-Beta function [4]: Bq (x, y) = (1 − q) 2. Discussion on regularity of solutions
t 6= 0. Ra
g(s)dq s − P∞
n=0 q
0
g(s)dq s,
Q∞
Rt
0 g(s)dq s := (1 − q)
1−q n+1+i i=0 1−q n+y+i .
nx
∞ P
tq n g(tq n ).
n=0
We notice that the solution of the initial value problem (1.1) is not smooth enough for some general cases of f (t, y(t)). For instance, consider the solutions of the two examples given in [26], [26, Example Γq (2) 3 Γ (13/6) 7 1]: x(t) = Γq (5/2) t 2 ; [26, Example 2]: x(t) = Γqq (11/6) t 6 + 0.0001. For both of them, the parameter M1 = sup0≤t≤1 |Dq2 x(t)| fails to be bounded which was assumed in the truncation error analysis of [26, 1
5
Theorem 3.2]. Precisely, one can see that M1 = C1 t− 2 for [26, Example 1] and M1 = C2 t− 6 for [26, Example 2]. They all fail to be bounded at t = 0. Referring to some usual assumptions on regularities of solutions to nonlinear fractional ordinary differential equations [16, 17, 19], the convergence analysis in the next section will be based on the following: k = 0, 1, 2,
t ∈ (0, 1],
na lP repr oo f
|z (k) (t)| ≤ Ctα−k ,
(2.2)
where z(t) = c Dqα y(t) and C is a generic positive constant which is independent of t. To demonstrate the assumption (2.2) may be generated from solutions with typical weak singularities, we recall that every solution of the q-fractional initial value problem (1.1) is also a solution of the following q-Volterra equation and vice versa [25]: y(t) =
dαe−1
X j=0
Z
1 t(j) Dj y(0) + Γq (j + 1) q Γq (α)
t
0
(t − qs)(α−1) f (s, y(s)) dq s.
(2.3)
Consider an usual and simplified case of (1.1): c
Dqα y(t) = f (t, y(t)) =
X
fij ti+jα + fm (t),
0≤i+jα
i, j, m ∈ N, t ∈ (0, 1],
(2.4)
where fij are constants, fm ∈ C m [0, 1] and (i, j) 6= (1, 0) when α ∈ (1, 2). We can see that condition (2.2) holds for (2.4). From (2.3) and (2.4), we have y(t) =
X
fij Γ (α) −α
Z
t
0
(t − qs)(α−1) si+jα dq s +
1 Γq (α)
Z
0
t
(t − qs)(α−1) fm (s) dq s +
dαe−1
X
k=0
t(k) Dk y(0). Γq (k + 1) q (2.5)
The definitions of q-integral and q-beta function yield Z
0
t
(t − qs)(α−1) si+jα dq s =
Jo
ur
1 Γq (α)
= =
∞ X 1 (1 − q) tq n (t − tq n+1 )(α−1) (tq n )i+jα Γq (α) n=0
∞ ∞ X Y ti+(j+1)α 1 − q n+1+r (1 − q) q n(1+i+jα) Γq (α) 1 − q n+α+r n=0 r=0
ti+(j+1)α Γq (1 + i + jα) Bq (α, 1 + i + jα) = ti+(j+1)α . Γq (α) Γq (1 + i + (j + 1)α)
(2.6)
where the property [4] Bq (x, y) = Γq (x)Γq (y)/Γq (x + y) has been utilized. As fm (t) ∈ C m [0, 1], it has the Taylor expansion 0 fm (t) = fm (0) + fm (0)t + · · · +
(m−1)
(m)
fm (0) m−1 fm (ξ) m t + t , (m − 1)! m!
ξ ∈ [0, t].
(2.7)
P Combining (2.5)–(2.7), we obtain that y(t) = 0≤i+jα
Journal Pre-proof
3. The second-order method Implementation of numerical schemes on nonuniform mesh was verified to be an extremely effective way to deal with fractional differential equations with non-smooth solutions, e.g. [10, 14–18, 21, 23, 24]. In this section, we will propose and analyze a scheme on nonuniform mesh for solving (2.3) which provides the same solution to (1.1). Numerical approximation will be conducted on the graded mesh: γ tk = Nk , γ ≥ 1 and 0 ≤ k ≤ N . It is clear that the step size satisfies γ(k − 1)γ−1 N −γ ≤ τk ≤ γk γ−1 N −γ , where τk = tk − tk−1 . Consider the equation (2.3) at the grid point tn (1 ≤ n ≤ N ): y(tn ) =
dαe−1
X j=0
t(j) 1 Dj y(0) + Γq (j + 1) q Γq (α)
Z
tn
0
(tn − qs)(α−1) f (s, y(s)) dq s.
(3.8)
Define the linear interpolation operators: t − tk t − tk−1 z(tk−1 ) + z(tk ), t ∈ [tk−1 , tk ], 1 ≤ k ≤ n, tk−1 − tk tk − tk−1 ˜ 1,n z(t) = t − tn−1 z(tn−2 ) + t − tn−2 z(tn−1 ), t ∈ (tn−1 , tn ]. Π tn−2 − tn−1 tn−1 − tn−2
Π1,k z(t) =
By using the Taylor formula with integral remainder, it is easy to find that
na lP repr oo f
Z 1 00 rk (t) , z(t) − Π1,k z(t) ≤ z (tξ)ξ dξ , 0 Z 1 2 00 ˜ r˜n (t) , z(t) − Π1,n z(t) ≤ Cτn z (tξ)ξ dξ , Cτk2
0
t ∈ (tk−1 , tk ),
(3.9)
t ∈ (tn−1 , tn ].
(3.10)
To construct an efficient numerical method, without loss of generality, we suppose f (t, y(t)) = ω(t)y(t) + g(t) + h(t)F (y(t)), where F (y) 6= c1 + c2 y and |ω(t)| ≤ κ for a positive constant κ. Denote ωk = ω(tk ), k ≥ 0, and similar definitions are applied to notations gk and hk . Throughout the paper, the sums and integrals are always set to zero if the upper index is less than the lower one. We approximate the integration term of (3.8) as follows: 1 Γq (α)
Z
tn
0
(tn − qs)(α−1) f (s, y(s)) dq s =
Z
0
tn
(tn − qs)(α−1) [ω(s)y(s) + g(s) + h(s)F (y(s))] dq s
Z t1 1 (tn − qs)(α−1) [h(t0 )F (y(t0 ))] dq s Γ (α) q t0 k=1 tk−1 Z tk Z tn n−1 X 1 1 ˜ 1,n [h(s)F (y(s))] dq s + (tn − qs)(α−1) Π (tn − qs)(α−1) Π1,k [h(s)F (y(s))] dq s + Γq (α) Γ q (α) tmax{n−1,1} tk−1 ≈
n
1 X Γq (α)
Z
1 Γq (α)
tk
(tn − qs)(α−1) Π1,k [ω(s)y(s) + g(s)] dq s +
k=2
,
n X
(q)
Pk,n [ωk y(tk ) + gk ] +
k=0
n−1 X
(q)
Qk,n [hk F (y(tk ))],
k=0
(q)
1 ≤ n ≤ N,
(3.11)
(q)
1 Γq (α)
Z
tj
tj−1
(tn − qs)(α−1)
Jo
aj,n =
ur
n , 0 ≤ k ≤ n − 1, with where Pk,n = ak+1,n + bk,n , 0 ≤ k ≤ n; and Qk,n = lkn + ck+1,n + dk,n + lk−n n n aj,n , bj,n , cj,n , dj,n , lj , lj−n are given by the following integrals:
s − tj dq s, −τj
bj,n =
1 Γq (α) ljn
Z
tj
tj−1
(tn − qs)(α−1)
1 = Γq (α)
Z
s − tj−1 dq s τj
for 1 ≤ j ≤ n;
t1
cj,n = aj,n , dj,n = bj,n for 2 ≤ j ≤ n − 1, n ≥ 3; (tn − qs)(α−1) dq s for j = 0; t0 R tn (α−1) s−tn−2 1 j = n − 1, Γq (α) tn−1 (tn − qs) τn−1 dq s, n lj−n = for n ≥ 2; R tn s−t 1 n−1 (α−1) j = n − 2, Γq (α) tn−1 (tn − qs) −τn−1 dq s,
otherwise they are all assigned to be 0. We remark that the Newton linearized method (see [11, 12]) can also be employed for the approximation of the nonlinear function F (y) in (3.11) on the interval 3
Journal Pre-proof
[tmax{n−1,1} , tn ] to obtain a linearized numerical scheme which can achieve similar convergence based on our framework. Let yk be the numerical approximation of y(tk ), k ≥ 1. Therefore, we obtain a linearized finite difference scheme to solve equation (1.1): yn =
dαe−1
X j=0
n
n−1
k=0
k=0
X (q) X (q) t(j) Dqj y(0) + Pk,n (ωk yk + gk ) + Qk,n [hk F (yk )], Γq (j + 1)
1 ≤ n ≤ N.
(3.12)
Denote en = u(tn ) − un , n ≥ 0. Subtracting (3.12) from (3.8), we obtain the following error equation: en =
n X
(q)
Pk,n ωk ek +
k=0
n−1 X k=0
(q)
Qk,n hk [F (y(tk )) − F (yk )] + Rn ,
1 ≤ n ≤ N,
(3.13)
where the truncation error follows from (3.11): n Z 1 X tk (tn − qs)(α−1) (κy(s) + g(s) − Π1,k [κy(s) + g(s)]) dq s Γq (α) t k=1 k−1 Z t1 1 (tn − qs)α−1 [h(s)F (y(s)) − h(t0 )F (y(t0 ))] dq s + Γq (α) 0 n−1 Z 1 X tk + (tn − qs)α−1 [h(s)F (y(s)) − Π2,k h(s)F (y(s)))] dq s Γq (α) t k−1 k=2 Z tn h i 1 ˜ 2,n h(s)F (y(s)) dq s , I1 + I2 + I3 + I4 . + (tn − s)α−1 h(s)F (y(s)) − Π Γq (α) tn−1
na lP repr oo f
Rn =
Lemma 3.1. For α ∈ (0, 1) ∪ (1, 2), and 0 ≤ s ≤ t ≤ 1, it holds that (t − qs)(α−1) < Proof
For 0 < α < 1, we have
(q i+α − q i+1 )s q i+α − q i+1 1 − q i+1 t − q i+1 s = 1 + ≤ max{1, 1 + } = , t − q i+α s t − q i+α s 1 − q i+α 1 − q i+α
Then (t − qs)(α−1) = tα−1
Q∞
For 1 < α < 2, since t −
t−q i+1 s tα−1 i=0 t−q i+α s < 1−q α . q i+1 s < t − q i+α s, we
obtain (t − qs)(α−1) = tα−1
Lemma 3.2. For α ∈ (0, 1) ∪ (1, 2), and 0 ≤ s ≤ t ≤ 1, it holds that Proof
tα−1 1−q α .
We have
Rt 0
sα dq s = tα+1 (1 − q)
P∞
i=0 q
i(α+1)
Rt 0
i ≥ 0.
Q∞
t−q i+1 s i=0 t−q i+α s
< tα−1 .
sα dq s ≤ tα+1 .
1−q α+1 . = tα+1 1−q α+1 ≤ t
Lemma 3.3. [20] Let f (x) be a continuous function on a segment [0, a] (a > 0). Then for any Ra q ∈ (0, 1), it holds that ∃ ξ ∈ [0, a] : 0 f (s)dq s = af (ξ).
ur
Lemma 3.4. [20] Let f (x) and g(x) be some continuous functions on a segment [a, b]. Then there exists qˆ ∈ (0, 1), for any q ∈ (ˆ q , 1), it holds that ∃ ξ ∈ (a, b) : Iq (f g) = g(ξ)Iq (f ).
Jo
Denote the maximum norm error e∞ (N ) = max0≤n≤N {|en |}. Next theorem concludes the convergence. Theorem 3.1. Based on the assumption (2.2), for N ≥ 3Γ4κ , if the grading parameter fulfills γ = q (α) −2 max{1/α, 1} with α ∈ (0, 1) ∪ (1, 2), we have e∞ (N ) ≤ CN . Proof By Lemmas 3.1, 3.2 and Theorems 3.3 and 3.4, with (3.9)–(3.10) and the assumption (2.2), the terms in the truncation error Rn can be estimated as |I1 | ≤C
n Z X
k=1
tk
tk−1
(tn −
qs)(α−1) τk2 sα−2
dq s ≤ C
4
n X
k=1
τk2 tkα−2
Z
tk
tk−1
(tn − qs)(α−1) dq s
Journal Pre-proof
n X
≤Ctα−1 n
τk3 tα−2 k
≤
(
CN − min{γ(2α),2} ,
0 < α < 1,
−2
CN , 1 < α < 2; k=1 Z s Z t1 1 (α−1) [ht (ξ)F (y(ξ)) + h(ξ)Ft (y(ξ))] dξ dq s (tn − qs) |I2 | = Γq (α) 0 0 ( Z t1 CN −γ(2α) , 0 < α < 1, α+1 (tn − qs)(α−1) sα dq s ≤ Ctnα−1 t1 ≤ ≤C −γ(α+1) CN , 1 < α < 2; 0 ( n−1 X Z tk CN − min{γ(2α),2} , 0 < α < 1, |I3 | ≤ (tn − qs)(α−1) τk2 sα−2 dq s ≤ −2 CN , 1 < α < 2.; k=2 tk−1 ( Z tn CN − min{γ(2α),3} , 0 < α < 1, (tn − qs)(α−1) sα−2 dq s ≤ |I4 | ≤Cτn2 − min{γ(α+1),3} CN , 1 < α < 2. tn−1
Consequently, |Rn | ≤ |I1 | + |I2 | + |I3 | + |I4 | ≤ CN − min{γ(2α),γ(α+1),2} ,
α ∈ (0, 1) ∪ (1, 2).
(3.14)
(q)
na lP repr oo f
n , we may go Meanwhile, by Lemma 3.1, Theorems 3.3 and 3.4, with the definition of ak,n , bk,n , lkn , lk−n though derivations similar to that of [18, Lemma 2.1] to get (q)
{|Pk,n |, |Qk,n |} ≤ Cτk+1 (tn − tk )α−dαe , for α ∈ (0, 1) ∪ (1, 2), (q)
0 ≤ k ≤ n − 1.
(3.15)
α−1
1 and Pn,n = bn,n ≤ τΓnqtn(α) ≤ Γq (α)N . This yields that 1 − ωn Pn,n ≥ 1/4 provided N ≥ 3Γ4κ . q (α) Therefore, from (3.13)–(3.15), and utlizing the modified Gr¨onwall inequality [18, Lemma 3.2], we finally get
|en | ≤ CN −2
for γ = max{1/α, 1},
n = 1, 2, · · · , N.
4. Numerical experiments
Numerical examples are now carried out to illustrate the efficiency of the proposed scheme (3.12). All our tests were done in MATLAB R2017a with a laptop. Example 4.1. We consider problem (1.1) for
f (t, y(t)) = Γq (α + 1) + tα − t2α − y + y 2 ,
y(0) = Dqdαe−1 y(0) = 0.
ur
The exact solution of Example 4.1 is y = tα . Table 1 lists the numerical accuracy of scheme (3.12) for Example 4.1. We can see that the proposed scheme works very efficiently with second-order convergence rate for different α ∈ (0, 1) ∪ (1, 2). Particularly the scheme is implemented on uniform grid for α ∈ (1, 2). 1
Jo
We also apply our proposed scheme to solve the Example 1 in [26]: c Dq2 y(t) = t for 0 < t ≤ 1 and y(0) = 0. The exact solution is y(t) = Γq (2)t3/2 Γq (5/2). Numerical results listed in Table 2 imply that our method is third-order convergent for this example. This may be owing to its source term (f (t, y(t)) = t) being smoother than those restricted by (2.2). Meanwhile, we note that the absolute errors at the final time point are less than those listed in Table 1 of [26]. 5. Acknowledgements The authors would like to thank anonymous reviewers for their valuable comments and suggestions, which helped them to improve greatly the presentation and quality of the paper. 5
Journal Pre-proof
Table 1: Numerical accuracy of scheme (3.12) with q = 2/5, and γ = max{1/α, 1}. α = 0.2 e∞ (N ) Order 1.1048e-01 ∗ 2.9135e-02 1.92 8.3895e-03 1.80 2.3782e-03 1.82 6.4236e-04 1.89 α = 1.01 e∞ (N ) Order 6.7331e-03 ∗ 1.7238e-03 1.97 4.3329e-04 1.99 1.0813e-04 2.01 2.6961e-05 2.01
N 16 32 64 128 256 N 16 32 64 128 256
α = 0.6 e∞ (N ) 3.9331e-03 1.0746e-03 2.8238e-04 7.1951e-05 1.8072e-05 α = 1.5 e∞ (N ) 1.0157e-02 2.7234e-03 6.9959e-04 1.7650e-04 4.4255e-05
Order ∗ 1.87 1.93 1.97 1.99 Order ∗ 1.90 1.96 1.99 2.00
α = 0.9 e∞ (N ) 6.6738e-03 1.7303e-03 4.3817e-04 1.0970e-04 2.7387e-05 α = 1.9 e∞ (N ) 1.2051e-02 3.3065e-03 8.5872e-04 2.1784e-04 5.4773e-05
Order ∗ 1.95 1.98 2.00 2.00 Order ∗ 1.87 1.95 1.98 1.99
Table 2: Numerical accuracy of scheme (3.12) for solving Example 1 in [26].
References
Maximum-norm e∞ (N ) 2.1515e-02 2.6894e-03 3.3618e-04 4.2022e-05 5.2527e-06
errors Order ∗ 3.00 3.00 3.00 3.00
Errors at the terminal point |eN | 1.1102e-16 9.9920e-16 7.4385e-15 1.0769e-14 3.9968e-15
na lP repr oo f
N 8 16 32 64 128
Jo
ur
[1] T. Abdeljawad, D. Baleanu, Caputo q-fractional initial value problems and a q-analogue Mittag-Leffler function, Commun. Nonlinear Sci. Numer. Simulat., 16 (2011) 4682–4688. [2] T. Abdeljawad, B. Benli, D. Baleanu, A generalized q-Mittag-Leffler function by q-Caputo fractional linear equations, Abstr. Appl. Anal., 2012 (2012) 1–11. [3] R. Almeida, N. Martins, Existence results for fractional q-difference equations of order α ∈]2, 3[ with three-point boundary conditions, Commun. Nonlinear Sci. Numer. Simul., 19 (2014) 1675–1685. [4] R. Askey, The q-gamma and q-beta functions, Appl. Anan., 8 (1978) 125–141. [5] F.M. Atici, P.W. Eloe, Fractional q-calculus on a time scale, J. Nonlinear Math. Phys. 14 (2007) 341–352. [6] F. H. Jackson, On q-definite integrals, Quart. J. Pure Appl. Math., 41 (1910)193–203. [7] F. H. Jackson, On q-functions and a certain difference operator, Trans. Roy. Soc. Edin., 46 (1908) 253–281. [8] B. Jin, R. Lazarov, Z. Zhou, An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data. IMA J. Numer. Anal., 36 (2016) 197–221. [9] B. Jin, B. Li, Z. Zhou, Subdiffusion with a time-dependent coefficient: Analysis and numerical solution, Math. Comput., 88 (2019) 2157–2186. [10] C.P. Li, Q. Yi, A. Chen, Finite difference methods with non-uniform meshes for nonlinear fractional differential equations. J. Comput. Phys., 316 (2016) 614–631. [11] D. Li, H.L. Liao, W. Sun, J. Wang, J. Zhang, Analysis of L1-Galerkin FEMs for time-fractional nonlinear parabolic problems, Commun. Comput. Phys., 2018, in press. [12] D. Li, C. Wu, Z. Zhang, Linearized Galerkin FEMs for nonlinear time fractional parabolic problems with non-smooth solutions in time direction, J. Sci. Comput., 80 (2019) 403–419. [13] L. Li, D. Li, Exact solutions and numerical study of time fractional Burgers’ equations, Appl. Math. Lett., 100 (2020) 106011. [14] H.L. Liao, D. Li, J. Zhang, Sharp error estimate of a nonuniform L1 formula for time-fractional reactionsubdiffusion equations, SIAM J. Numer. Anal., 56 (2018) 1112–1133.
6
Journal Pre-proof
na lP repr oo f
[15] H.L. Liao, W. McLean, J. Zhang, A discrete Gr¨ onwall inequality with applications to numerical schemes for subdiffusion problems, SIAM J. Numer. Anal., 57 (2019) 218–237. [16] Y. Liu, J. Roberts, Y., Yan, A note on finite difference methods for nonlinear fractional differential equations with non-uniform meshes, Int. J. Comput. Math., 95 (2018)1151–1169. [17] Y. Liu, J. Roberts, Y., Yan, Detailed error analysis for a fractional Adams method with graded meshes, Numer. Algorithms, 78 (2018) 1195–1216. [18] P. Lyu, S. Vong, A high-order method with a temporal nonuniform mesh for a time-fractional BenjaminBona-Mahony equation, J. Sci. Comput, 80 (2019) 1607–1628. [19] A. Pedas, E. Tamme, Numerical solution of nonlinear fractional differential equations by spline collocation methods, J. Comput. Appl. Math., 255 (2014) 216–230. [20] P.M. Rajkovi´c, M.S. Stankovi´c, S.D. Marinkovi´c, Mean value theorems in q-Calculus, Math. Vesnic, 54 (2002) 171–178. [21] J. Ren, H. Chen, A numerical method for distributed order time fractional diffusion equation with weakly singular solutions, Appl. Math. Lett., 96 (2019) 159–165. [22] S. Salahshour, A. Ahmadian, C.S. Chan, Successive approximation method for Caputo q-fractional IVPs, Commun. Nonlinear Sci. Numer. Simul., 24 (2015) 153–158. [23] M. Stynes, J.L. Gracia, Preprocessing schemes for fractional-derivative problems to improve their convergence rate, Appl. Math. Lett., 74 (2017) 187–192. [24] 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 (2017) 1057–1079. [25] Y. Tang, T. Zhang, A remark on the q-fractional order differential equations, App. Math. Comput., 350 (2019) 198–208.
Jo
ur
[26] T. Zhang, Y. Tang, A difference method for solving the q-fractional differential equations, App. Math. Lett., 98 (2019) 292–299.
7
Journal Pre-proof
Jo
ur
na lP repr oo f
Both authors of the manuscript have engaged in the preparation of this manuscript. All of us have worked our best on analyzing, composing and reviewing the draft and the revision version.
Pin Lyu, Seakweng Vong
PII: DOI: Reference:
S0893-9659(19)30482-3 https://doi.org/10.1016/j.aml.2019.106156 AML 106156
To appear in:
Applied Mathematics Letters
Received date : 23 August 2019 Revised date : 18 November 2019 Accepted date : 18 November 2019 Please cite this article as: P. Lyu and S. Vong, An efficient numerical method for q-fractional differential equations, Applied Mathematics Letters (2019), doi: https://doi.org/10.1016/j.aml.2019.106156. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Elsevier Ltd. All rights reserved.
Journal Pre-proof
An efficient numerical method for q-fractional differential equations Pin Lyua , Seakweng Vongb,∗ a
School of Economic Mathematics, Southwestern University of Finance and Economics, Chengdu, China b Department of Mathematics, University of Macau, Macao, China
Abstract In this work, we propose and analyze an efficient numerical method for solving nonlinear Caputo q-fractional differential equations with non-smooth solutions. The numerical method is a linearized one and can be implemented on the graded mesh. For equations with fractional derivatives of order α ∈ (1, 2), the method converges with second-order on uniform mesh. For α ∈ (0, 1), second-order convergence can be obtained when the grading parameter of the nonuniform mesh is set appropriately.
1. Introduction
na lP repr oo f
Keywords: Caputo q-fractional derivative, nonuniform mesh, fractional nonlinear equation
Recently, q-fractional differential equations have attracted lots of interests in their analytical and numerical studies, see also [1–3, 22, 25, 26] and references therein. It is well known that the solutions of fractional (ordinary and partial) differential equations typically exhibit weak singularities near the initial/boundary point. This has become a major concern in the numerical analysis of fractional differential equations [8, 9, 13–19, 23, 24]. Similar weak singularities may also appear in q-fractional differential equations, referring to our brief discussion in the next section. A method to solve the q-fractional differential equations is proposed in a very recent work [26], in which the convergence analysis is absent. Furthermore, the truncation error analysis in [26] requires to assume that the solution is smooth enough on the whole closed domain. In view of this, based on some techniques of our recent work [18], we propose and analyze another efficient numerical method for the q-fractional differential equation: c
Dqα y(t) = f (t, y),
t ∈ (0, T ],
(1.1) dαe−1
subject to the initial condition y(0) = y0 if α ∈ (0, 1); and y(0) = y0 , Dq y(0) = yˆ0 if α ∈ (1, 2). Definitions of q-analogy (0 < |q| < 1) derivatives, integrals and related functions which will be utilized in this paper are listed below:
ur
• Caputo q-fractional derivative [1]:
c
• Riemann-Liouville q-fractional integral [5]:
Jo
• q-Gamma function [5]:
dαe−α
Dqα y(t) := Iq
dαe
Dq y(t). Rt Iqδ g(t) := Γq1(δ) 0 (t − qs)(δ−1) g(s)dq s, δ > 0.
Γq (δ) := (1 − q)1−δ (1 − q)(δ−1) .
• q-analogy of (t − s)(ν) [5]: (t − s)(ν) := tν
Q∞
t−q i s i=0 t−q ν+i s
for a non-positive integer ν.
I The first author is supported by the Fundamental Research Funds for the Central Universities (JBK1901027), the second author is supported by The Science and Technology Development Fund, Macau SAR (File no. 0005/2019/A) and the grants MYRG2018-00047-FST and MYRG2017-00098-FST from University of Macau. ∗ Corresponding author. Email addresses: [email protected] (Pin Lyu), [email protected] (Seakweng Vong)
Preprint submitted to Applied Mathematics Letters
December 2, 2019
Journal Pre-proof
• q-derivative [7] :
Dq g(t) :=
• q-integral [6]: Iq (g) :=
Rb a
g(qt)−g(t) (q−1)t ,
g(s)dq s =
Rb 0
• q-Beta function [4]: Bq (x, y) = (1 − q) 2. Discussion on regularity of solutions
t 6= 0. Ra
g(s)dq s − P∞
n=0 q
0
g(s)dq s,
Q∞
Rt
0 g(s)dq s := (1 − q)
1−q n+1+i i=0 1−q n+y+i .
nx
∞ P
tq n g(tq n ).
n=0
We notice that the solution of the initial value problem (1.1) is not smooth enough for some general cases of f (t, y(t)). For instance, consider the solutions of the two examples given in [26], [26, Example Γq (2) 3 Γ (13/6) 7 1]: x(t) = Γq (5/2) t 2 ; [26, Example 2]: x(t) = Γqq (11/6) t 6 + 0.0001. For both of them, the parameter M1 = sup0≤t≤1 |Dq2 x(t)| fails to be bounded which was assumed in the truncation error analysis of [26, 1
5
Theorem 3.2]. Precisely, one can see that M1 = C1 t− 2 for [26, Example 1] and M1 = C2 t− 6 for [26, Example 2]. They all fail to be bounded at t = 0. Referring to some usual assumptions on regularities of solutions to nonlinear fractional ordinary differential equations [16, 17, 19], the convergence analysis in the next section will be based on the following: k = 0, 1, 2,
t ∈ (0, 1],
na lP repr oo f
|z (k) (t)| ≤ Ctα−k ,
(2.2)
where z(t) = c Dqα y(t) and C is a generic positive constant which is independent of t. To demonstrate the assumption (2.2) may be generated from solutions with typical weak singularities, we recall that every solution of the q-fractional initial value problem (1.1) is also a solution of the following q-Volterra equation and vice versa [25]: y(t) =
dαe−1
X j=0
Z
1 t(j) Dj y(0) + Γq (j + 1) q Γq (α)
t
0
(t − qs)(α−1) f (s, y(s)) dq s.
(2.3)
Consider an usual and simplified case of (1.1): c
Dqα y(t) = f (t, y(t)) =
X
fij ti+jα + fm (t),
0≤i+jα
i, j, m ∈ N, t ∈ (0, 1],
(2.4)
where fij are constants, fm ∈ C m [0, 1] and (i, j) 6= (1, 0) when α ∈ (1, 2). We can see that condition (2.2) holds for (2.4). From (2.3) and (2.4), we have y(t) =
X
fij Γ (α) −α
Z
t
0
(t − qs)(α−1) si+jα dq s +
1 Γq (α)
Z
0
t
(t − qs)(α−1) fm (s) dq s +
dαe−1
X
k=0
t(k) Dk y(0). Γq (k + 1) q (2.5)
The definitions of q-integral and q-beta function yield Z
0
t
(t − qs)(α−1) si+jα dq s =
Jo
ur
1 Γq (α)
= =
∞ X 1 (1 − q) tq n (t − tq n+1 )(α−1) (tq n )i+jα Γq (α) n=0
∞ ∞ X Y ti+(j+1)α 1 − q n+1+r (1 − q) q n(1+i+jα) Γq (α) 1 − q n+α+r n=0 r=0
ti+(j+1)α Γq (1 + i + jα) Bq (α, 1 + i + jα) = ti+(j+1)α . Γq (α) Γq (1 + i + (j + 1)α)
(2.6)
where the property [4] Bq (x, y) = Γq (x)Γq (y)/Γq (x + y) has been utilized. As fm (t) ∈ C m [0, 1], it has the Taylor expansion 0 fm (t) = fm (0) + fm (0)t + · · · +
(m−1)
(m)
fm (0) m−1 fm (ξ) m t + t , (m − 1)! m!
ξ ∈ [0, t].
(2.7)
P Combining (2.5)–(2.7), we obtain that y(t) = 0≤i+jα
Journal Pre-proof
3. The second-order method Implementation of numerical schemes on nonuniform mesh was verified to be an extremely effective way to deal with fractional differential equations with non-smooth solutions, e.g. [10, 14–18, 21, 23, 24]. In this section, we will propose and analyze a scheme on nonuniform mesh for solving (2.3) which provides the same solution to (1.1). Numerical approximation will be conducted on the graded mesh: γ tk = Nk , γ ≥ 1 and 0 ≤ k ≤ N . It is clear that the step size satisfies γ(k − 1)γ−1 N −γ ≤ τk ≤ γk γ−1 N −γ , where τk = tk − tk−1 . Consider the equation (2.3) at the grid point tn (1 ≤ n ≤ N ): y(tn ) =
dαe−1
X j=0
t(j) 1 Dj y(0) + Γq (j + 1) q Γq (α)
Z
tn
0
(tn − qs)(α−1) f (s, y(s)) dq s.
(3.8)
Define the linear interpolation operators: t − tk t − tk−1 z(tk−1 ) + z(tk ), t ∈ [tk−1 , tk ], 1 ≤ k ≤ n, tk−1 − tk tk − tk−1 ˜ 1,n z(t) = t − tn−1 z(tn−2 ) + t − tn−2 z(tn−1 ), t ∈ (tn−1 , tn ]. Π tn−2 − tn−1 tn−1 − tn−2
Π1,k z(t) =
By using the Taylor formula with integral remainder, it is easy to find that
na lP repr oo f
Z 1 00 rk (t) , z(t) − Π1,k z(t) ≤ z (tξ)ξ dξ , 0 Z 1 2 00 ˜ r˜n (t) , z(t) − Π1,n z(t) ≤ Cτn z (tξ)ξ dξ , Cτk2
0
t ∈ (tk−1 , tk ),
(3.9)
t ∈ (tn−1 , tn ].
(3.10)
To construct an efficient numerical method, without loss of generality, we suppose f (t, y(t)) = ω(t)y(t) + g(t) + h(t)F (y(t)), where F (y) 6= c1 + c2 y and |ω(t)| ≤ κ for a positive constant κ. Denote ωk = ω(tk ), k ≥ 0, and similar definitions are applied to notations gk and hk . Throughout the paper, the sums and integrals are always set to zero if the upper index is less than the lower one. We approximate the integration term of (3.8) as follows: 1 Γq (α)
Z
tn
0
(tn − qs)(α−1) f (s, y(s)) dq s =
Z
0
tn
(tn − qs)(α−1) [ω(s)y(s) + g(s) + h(s)F (y(s))] dq s
Z t1 1 (tn − qs)(α−1) [h(t0 )F (y(t0 ))] dq s Γ (α) q t0 k=1 tk−1 Z tk Z tn n−1 X 1 1 ˜ 1,n [h(s)F (y(s))] dq s + (tn − qs)(α−1) Π (tn − qs)(α−1) Π1,k [h(s)F (y(s))] dq s + Γq (α) Γ q (α) tmax{n−1,1} tk−1 ≈
n
1 X Γq (α)
Z
1 Γq (α)
tk
(tn − qs)(α−1) Π1,k [ω(s)y(s) + g(s)] dq s +
k=2
,
n X
(q)
Pk,n [ωk y(tk ) + gk ] +
k=0
n−1 X
(q)
Qk,n [hk F (y(tk ))],
k=0
(q)
1 ≤ n ≤ N,
(3.11)
(q)
1 Γq (α)
Z
tj
tj−1
(tn − qs)(α−1)
Jo
aj,n =
ur
n , 0 ≤ k ≤ n − 1, with where Pk,n = ak+1,n + bk,n , 0 ≤ k ≤ n; and Qk,n = lkn + ck+1,n + dk,n + lk−n n n aj,n , bj,n , cj,n , dj,n , lj , lj−n are given by the following integrals:
s − tj dq s, −τj
bj,n =
1 Γq (α) ljn
Z
tj
tj−1
(tn − qs)(α−1)
1 = Γq (α)
Z
s − tj−1 dq s τj
for 1 ≤ j ≤ n;
t1
cj,n = aj,n , dj,n = bj,n for 2 ≤ j ≤ n − 1, n ≥ 3; (tn − qs)(α−1) dq s for j = 0; t0 R tn (α−1) s−tn−2 1 j = n − 1, Γq (α) tn−1 (tn − qs) τn−1 dq s, n lj−n = for n ≥ 2; R tn s−t 1 n−1 (α−1) j = n − 2, Γq (α) tn−1 (tn − qs) −τn−1 dq s,
otherwise they are all assigned to be 0. We remark that the Newton linearized method (see [11, 12]) can also be employed for the approximation of the nonlinear function F (y) in (3.11) on the interval 3
Journal Pre-proof
[tmax{n−1,1} , tn ] to obtain a linearized numerical scheme which can achieve similar convergence based on our framework. Let yk be the numerical approximation of y(tk ), k ≥ 1. Therefore, we obtain a linearized finite difference scheme to solve equation (1.1): yn =
dαe−1
X j=0
n
n−1
k=0
k=0
X (q) X (q) t(j) Dqj y(0) + Pk,n (ωk yk + gk ) + Qk,n [hk F (yk )], Γq (j + 1)
1 ≤ n ≤ N.
(3.12)
Denote en = u(tn ) − un , n ≥ 0. Subtracting (3.12) from (3.8), we obtain the following error equation: en =
n X
(q)
Pk,n ωk ek +
k=0
n−1 X k=0
(q)
Qk,n hk [F (y(tk )) − F (yk )] + Rn ,
1 ≤ n ≤ N,
(3.13)
where the truncation error follows from (3.11): n Z 1 X tk (tn − qs)(α−1) (κy(s) + g(s) − Π1,k [κy(s) + g(s)]) dq s Γq (α) t k=1 k−1 Z t1 1 (tn − qs)α−1 [h(s)F (y(s)) − h(t0 )F (y(t0 ))] dq s + Γq (α) 0 n−1 Z 1 X tk + (tn − qs)α−1 [h(s)F (y(s)) − Π2,k h(s)F (y(s)))] dq s Γq (α) t k−1 k=2 Z tn h i 1 ˜ 2,n h(s)F (y(s)) dq s , I1 + I2 + I3 + I4 . + (tn − s)α−1 h(s)F (y(s)) − Π Γq (α) tn−1
na lP repr oo f
Rn =
Lemma 3.1. For α ∈ (0, 1) ∪ (1, 2), and 0 ≤ s ≤ t ≤ 1, it holds that (t − qs)(α−1) < Proof
For 0 < α < 1, we have
(q i+α − q i+1 )s q i+α − q i+1 1 − q i+1 t − q i+1 s = 1 + ≤ max{1, 1 + } = , t − q i+α s t − q i+α s 1 − q i+α 1 − q i+α
Then (t − qs)(α−1) = tα−1
Q∞
For 1 < α < 2, since t −
t−q i+1 s tα−1 i=0 t−q i+α s < 1−q α . q i+1 s < t − q i+α s, we
obtain (t − qs)(α−1) = tα−1
Lemma 3.2. For α ∈ (0, 1) ∪ (1, 2), and 0 ≤ s ≤ t ≤ 1, it holds that Proof
tα−1 1−q α .
We have
Rt 0
sα dq s = tα+1 (1 − q)
P∞
i=0 q
i(α+1)
Rt 0
i ≥ 0.
Q∞
t−q i+1 s i=0 t−q i+α s
< tα−1 .
sα dq s ≤ tα+1 .
1−q α+1 . = tα+1 1−q α+1 ≤ t
Lemma 3.3. [20] Let f (x) be a continuous function on a segment [0, a] (a > 0). Then for any Ra q ∈ (0, 1), it holds that ∃ ξ ∈ [0, a] : 0 f (s)dq s = af (ξ).
ur
Lemma 3.4. [20] Let f (x) and g(x) be some continuous functions on a segment [a, b]. Then there exists qˆ ∈ (0, 1), for any q ∈ (ˆ q , 1), it holds that ∃ ξ ∈ (a, b) : Iq (f g) = g(ξ)Iq (f ).
Jo
Denote the maximum norm error e∞ (N ) = max0≤n≤N {|en |}. Next theorem concludes the convergence. Theorem 3.1. Based on the assumption (2.2), for N ≥ 3Γ4κ , if the grading parameter fulfills γ = q (α) −2 max{1/α, 1} with α ∈ (0, 1) ∪ (1, 2), we have e∞ (N ) ≤ CN . Proof By Lemmas 3.1, 3.2 and Theorems 3.3 and 3.4, with (3.9)–(3.10) and the assumption (2.2), the terms in the truncation error Rn can be estimated as |I1 | ≤C
n Z X
k=1
tk
tk−1
(tn −
qs)(α−1) τk2 sα−2
dq s ≤ C
4
n X
k=1
τk2 tkα−2
Z
tk
tk−1
(tn − qs)(α−1) dq s
Journal Pre-proof
n X
≤Ctα−1 n
τk3 tα−2 k
≤
(
CN − min{γ(2α),2} ,
0 < α < 1,
−2
CN , 1 < α < 2; k=1 Z s Z t1 1 (α−1) [ht (ξ)F (y(ξ)) + h(ξ)Ft (y(ξ))] dξ dq s (tn − qs) |I2 | = Γq (α) 0 0 ( Z t1 CN −γ(2α) , 0 < α < 1, α+1 (tn − qs)(α−1) sα dq s ≤ Ctnα−1 t1 ≤ ≤C −γ(α+1) CN , 1 < α < 2; 0 ( n−1 X Z tk CN − min{γ(2α),2} , 0 < α < 1, |I3 | ≤ (tn − qs)(α−1) τk2 sα−2 dq s ≤ −2 CN , 1 < α < 2.; k=2 tk−1 ( Z tn CN − min{γ(2α),3} , 0 < α < 1, (tn − qs)(α−1) sα−2 dq s ≤ |I4 | ≤Cτn2 − min{γ(α+1),3} CN , 1 < α < 2. tn−1
Consequently, |Rn | ≤ |I1 | + |I2 | + |I3 | + |I4 | ≤ CN − min{γ(2α),γ(α+1),2} ,
α ∈ (0, 1) ∪ (1, 2).
(3.14)
(q)
na lP repr oo f
n , we may go Meanwhile, by Lemma 3.1, Theorems 3.3 and 3.4, with the definition of ak,n , bk,n , lkn , lk−n though derivations similar to that of [18, Lemma 2.1] to get (q)
{|Pk,n |, |Qk,n |} ≤ Cτk+1 (tn − tk )α−dαe , for α ∈ (0, 1) ∪ (1, 2), (q)
0 ≤ k ≤ n − 1.
(3.15)
α−1
1 and Pn,n = bn,n ≤ τΓnqtn(α) ≤ Γq (α)N . This yields that 1 − ωn Pn,n ≥ 1/4 provided N ≥ 3Γ4κ . q (α) Therefore, from (3.13)–(3.15), and utlizing the modified Gr¨onwall inequality [18, Lemma 3.2], we finally get
|en | ≤ CN −2
for γ = max{1/α, 1},
n = 1, 2, · · · , N.
4. Numerical experiments
Numerical examples are now carried out to illustrate the efficiency of the proposed scheme (3.12). All our tests were done in MATLAB R2017a with a laptop. Example 4.1. We consider problem (1.1) for
f (t, y(t)) = Γq (α + 1) + tα − t2α − y + y 2 ,
y(0) = Dqdαe−1 y(0) = 0.
ur
The exact solution of Example 4.1 is y = tα . Table 1 lists the numerical accuracy of scheme (3.12) for Example 4.1. We can see that the proposed scheme works very efficiently with second-order convergence rate for different α ∈ (0, 1) ∪ (1, 2). Particularly the scheme is implemented on uniform grid for α ∈ (1, 2). 1
Jo
We also apply our proposed scheme to solve the Example 1 in [26]: c Dq2 y(t) = t for 0 < t ≤ 1 and y(0) = 0. The exact solution is y(t) = Γq (2)t3/2 Γq (5/2). Numerical results listed in Table 2 imply that our method is third-order convergent for this example. This may be owing to its source term (f (t, y(t)) = t) being smoother than those restricted by (2.2). Meanwhile, we note that the absolute errors at the final time point are less than those listed in Table 1 of [26]. 5. Acknowledgements The authors would like to thank anonymous reviewers for their valuable comments and suggestions, which helped them to improve greatly the presentation and quality of the paper. 5
Journal Pre-proof
Table 1: Numerical accuracy of scheme (3.12) with q = 2/5, and γ = max{1/α, 1}. α = 0.2 e∞ (N ) Order 1.1048e-01 ∗ 2.9135e-02 1.92 8.3895e-03 1.80 2.3782e-03 1.82 6.4236e-04 1.89 α = 1.01 e∞ (N ) Order 6.7331e-03 ∗ 1.7238e-03 1.97 4.3329e-04 1.99 1.0813e-04 2.01 2.6961e-05 2.01
N 16 32 64 128 256 N 16 32 64 128 256
α = 0.6 e∞ (N ) 3.9331e-03 1.0746e-03 2.8238e-04 7.1951e-05 1.8072e-05 α = 1.5 e∞ (N ) 1.0157e-02 2.7234e-03 6.9959e-04 1.7650e-04 4.4255e-05
Order ∗ 1.87 1.93 1.97 1.99 Order ∗ 1.90 1.96 1.99 2.00
α = 0.9 e∞ (N ) 6.6738e-03 1.7303e-03 4.3817e-04 1.0970e-04 2.7387e-05 α = 1.9 e∞ (N ) 1.2051e-02 3.3065e-03 8.5872e-04 2.1784e-04 5.4773e-05
Order ∗ 1.95 1.98 2.00 2.00 Order ∗ 1.87 1.95 1.98 1.99
Table 2: Numerical accuracy of scheme (3.12) for solving Example 1 in [26].
References
Maximum-norm e∞ (N ) 2.1515e-02 2.6894e-03 3.3618e-04 4.2022e-05 5.2527e-06
errors Order ∗ 3.00 3.00 3.00 3.00
Errors at the terminal point |eN | 1.1102e-16 9.9920e-16 7.4385e-15 1.0769e-14 3.9968e-15
na lP repr oo f
N 8 16 32 64 128
Jo
ur
[1] T. Abdeljawad, D. Baleanu, Caputo q-fractional initial value problems and a q-analogue Mittag-Leffler function, Commun. Nonlinear Sci. Numer. Simulat., 16 (2011) 4682–4688. [2] T. Abdeljawad, B. Benli, D. Baleanu, A generalized q-Mittag-Leffler function by q-Caputo fractional linear equations, Abstr. Appl. Anal., 2012 (2012) 1–11. [3] R. Almeida, N. Martins, Existence results for fractional q-difference equations of order α ∈]2, 3[ with three-point boundary conditions, Commun. Nonlinear Sci. Numer. Simul., 19 (2014) 1675–1685. [4] R. Askey, The q-gamma and q-beta functions, Appl. Anan., 8 (1978) 125–141. [5] F.M. Atici, P.W. Eloe, Fractional q-calculus on a time scale, J. Nonlinear Math. Phys. 14 (2007) 341–352. [6] F. H. Jackson, On q-definite integrals, Quart. J. Pure Appl. Math., 41 (1910)193–203. [7] F. H. Jackson, On q-functions and a certain difference operator, Trans. Roy. Soc. Edin., 46 (1908) 253–281. [8] B. Jin, R. Lazarov, Z. Zhou, An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data. IMA J. Numer. Anal., 36 (2016) 197–221. [9] B. Jin, B. Li, Z. Zhou, Subdiffusion with a time-dependent coefficient: Analysis and numerical solution, Math. Comput., 88 (2019) 2157–2186. [10] C.P. Li, Q. Yi, A. Chen, Finite difference methods with non-uniform meshes for nonlinear fractional differential equations. J. Comput. Phys., 316 (2016) 614–631. [11] D. Li, H.L. Liao, W. Sun, J. Wang, J. Zhang, Analysis of L1-Galerkin FEMs for time-fractional nonlinear parabolic problems, Commun. Comput. Phys., 2018, in press. [12] D. Li, C. Wu, Z. Zhang, Linearized Galerkin FEMs for nonlinear time fractional parabolic problems with non-smooth solutions in time direction, J. Sci. Comput., 80 (2019) 403–419. [13] L. Li, D. Li, Exact solutions and numerical study of time fractional Burgers’ equations, Appl. Math. Lett., 100 (2020) 106011. [14] H.L. Liao, D. Li, J. Zhang, Sharp error estimate of a nonuniform L1 formula for time-fractional reactionsubdiffusion equations, SIAM J. Numer. Anal., 56 (2018) 1112–1133.
6
Journal Pre-proof
na lP repr oo f
[15] H.L. Liao, W. McLean, J. Zhang, A discrete Gr¨ onwall inequality with applications to numerical schemes for subdiffusion problems, SIAM J. Numer. Anal., 57 (2019) 218–237. [16] Y. Liu, J. Roberts, Y., Yan, A note on finite difference methods for nonlinear fractional differential equations with non-uniform meshes, Int. J. Comput. Math., 95 (2018)1151–1169. [17] Y. Liu, J. Roberts, Y., Yan, Detailed error analysis for a fractional Adams method with graded meshes, Numer. Algorithms, 78 (2018) 1195–1216. [18] P. Lyu, S. Vong, A high-order method with a temporal nonuniform mesh for a time-fractional BenjaminBona-Mahony equation, J. Sci. Comput, 80 (2019) 1607–1628. [19] A. Pedas, E. Tamme, Numerical solution of nonlinear fractional differential equations by spline collocation methods, J. Comput. Appl. Math., 255 (2014) 216–230. [20] P.M. Rajkovi´c, M.S. Stankovi´c, S.D. Marinkovi´c, Mean value theorems in q-Calculus, Math. Vesnic, 54 (2002) 171–178. [21] J. Ren, H. Chen, A numerical method for distributed order time fractional diffusion equation with weakly singular solutions, Appl. Math. Lett., 96 (2019) 159–165. [22] S. Salahshour, A. Ahmadian, C.S. Chan, Successive approximation method for Caputo q-fractional IVPs, Commun. Nonlinear Sci. Numer. Simul., 24 (2015) 153–158. [23] M. Stynes, J.L. Gracia, Preprocessing schemes for fractional-derivative problems to improve their convergence rate, Appl. Math. Lett., 74 (2017) 187–192. [24] 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 (2017) 1057–1079. [25] Y. Tang, T. Zhang, A remark on the q-fractional order differential equations, App. Math. Comput., 350 (2019) 198–208.
Jo
ur
[26] T. Zhang, Y. Tang, A difference method for solving the q-fractional differential equations, App. Math. Lett., 98 (2019) 292–299.
7
Journal Pre-proof
Jo
ur
na lP repr oo f
Both authors of the manuscript have engaged in the preparation of this manuscript. All of us have worked our best on analyzing, composing and reviewing the draft and the revision version.