- Email: [email protected]
Galerkin spectral method for solving two-dimensional time-space distributed-order weakly singular integro-partial differential equation
Journal Pre-proof Crank–Nicolson/Galerkin spectral method for solving two-dimensional time-space distributed-order weakly singular integro-partial differential equation Mostafa Abbaszadeh, Mehdi Dehghan, Yong Zhou
PII: DOI: Reference:
S0377-0427(20)30030-3 https://doi.org/10.1016/j.cam.2020.112739 CAM 112739
To appear in:
Journal of Computational and Applied Mathematics
Received date : 21 July 2019 Revised date : 9 January 2020 Please cite this article as: M. Abbaszadeh, M. Dehghan and Y. Zhou, Crank–Nicolson/Galerkin spectral method for solving two-dimensional time-space distributed-order weakly singular integro-partial differential equation, Journal of Computational and Applied Mathematics (2020), doi: https://doi.org/10.1016/j.cam.2020.112739. 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.
© 2020 Elsevier B.V. All rights reserved.
Journal Pre-proof
of
Crank–Nicolson/Galerkin spectral method for solving two-dimensional time-space distributed-order weakly singular integro-partial differential equation Mostafa Abbaszadeha∗, Mehdi Dehghana , Yong Zhoub
p ro
a Department of Applied Mathematics, Faculty of Mathematics and Computer Sciences,
Amirkabir University of Technology, No. 424, Hafez Ave.,15914, Tehran, Iran b School of Mathematics and Computational Science,
Xiangtan University, Xiangtan 411105, Hunan, PR China
Pr e-
January 20, 2020
Abstract
al
The fractional PDEs based upon the distributed-order fractional derivative have several applications in physics. The two-dimensional time-space distributed-order weakly singular integro-partial differential model is investigated by a combination of finite difference and Galekrin spectral methods. A second-order finite difference formula is employed to approximate the temporal variable. In this stage, the stability and convergence of the semi-discrete scheme are proved. Then, the Galerkin spectral method based on the modified Jacobi polynomials is applied to discrete the space variable. Also, in this step, the error estimate of the full-discrete scheme is studied. Finally, two test problems have been presented to confirm the theoretical results.
urn
Keywords: Time-space fractional derivatives, convergence analysis and error estimate, finite element method, stability analysis, finite difference method.
Jo
AMS subject Classification: 65M70, 34A34.
1
Introduction
The mathematical models based upon the fractional calculations are developing in most branches of physics or engineering . ∗
Corresponding author. E-mail addresses: [email protected] , [email protected] (M. Dehghan), [email protected] (M.Abbaszadeh), [email protected] (Y. Zhou).
1
Journal Pre-proof
Pr e-
p ro
of
The fractional PDEs have been investigated by several numerical procedures and their modifications such as finite difference methods [4, 5, 12, 20, 21, 22, 37, 49, 54], spectral methods and their improvements [7, 38, 42, 60, 61], the point interpolation technique [59], FEM on complex regions [25, 53], the ADI scheme with high accuracy [62], finite element method [8], etc. Also, the interested readers can refer to [19, 29, 41, 52, 55, 56, 57, 58]. The main aim of [63] is to study the existence and uniqueness of local and global mild solutions of the Navier–Stokes equations with time-fractional derivative of order 0 < α < 1. Authors of [64] studied the time-fractional reaction-diffusion equation with nonlocal boundary condition. The main aim of [65] is to establish sufficient conditions for the existence of globally attractive solutions for the Cauchy problems in cases that semigroup is compact as well as noncompact. An error estimate has been proposed in [1] to find a second-order finite difference scheme for solving the Riesz space distributed-order diffusion equation. Authors of [2] developed a meshless upwind local radial basis function-finite difference technique for solving the time- fractional distributed-order advection–diffusion equation. The main aim of [3] is to solve the multi-dimensional generalized modified anomalous sub-diffusion equation by using alternating direction implicit-spectral element method (ADI-SEM). An error estimate has been proposed in [15] to solve the 2D weakly singular integro-partial differential equation with space and time fractional derivatives based on the finite element/finite difference scheme. Authors of [18] employed the homotopy analysis method to solve some nonlinear fractional models such as fractional KdV, K(2,2), Burgers, BBM-Burgers, cubic Boussinesq, coupled KdV, and Boussinesq-like B(m,n) equations with initial conditions. Authors of [26] found a special point for the interpolation approximation of the linear combination of multi-term fractional derivatives as based on this issue they derived a numerical differentiation formula with second-order accuracy for solving time multi-term and distributed-order fractional sub-diffusion equations. In this paper, the following two-dimensional time-space distributed-order integro-partial differential equation with the weakly singular kernel is studied [14, 35, 40] Z2
al
Z1
ω(α)0 Dαt u(x, y, t)dα + u(x, y, t) = K
1
+
Zt 0
urn
0
$(β)
− 12
(t − ζ)
∂ β u(x, y, t) ∂|x|β
+
∂ β u(x, y, t) ∂|y|β
∂ 2 u(x, y, ζ) ∂ 2 u(x, y, ζ) + dζ + f (x, y, t), ∂x2 ∂y 2
!
dβ
(x, y, t) ∈ Ω × (0, T ], (1.1)
with the initial and boundary conditions
(x, y) ∈ Ω,
(1.2)
u(x, y, t) = g(x, y, t),
(x, y) ∈ ∂Ω,
(1.3)
Jo
u(x, y, 0) = ϕ(x, y),
in which K is Newtonian contribution to the viscosity. Also, 0 Dαt u(x, y, t) is the Riemann-Liouville fractional derivative. Eq. (1.1) may be appeared in the mathematical modeling of various physical
2
Journal Pre-proof
2
p ro
of
phenomena such as the heat transfer in the materials with memory, convection and radiation systems [35] and viscoelastic mechanics [14, 16, 40]. The fractional equations with the distributedorder have been studied by many researchers for example Caputo [9] proposed the application of differential equations with distributed-order derivatives for generalizing stress-strain relations of unelastic media. Also, Caputo in [10, 11] discussed distributed-order time fractional differential equations and distributed-order space fractional differential equations, respectively and derived the solutions with closed form formulae of the classic problems. Authors of [13, 43] introduced the diffusion-like equations with distributed-order time and space fractional derivatives for the kinetic description of anomalous diffusion and relaxation phenomena [28]. Author of [27] applied distributed-order diffusion equation to discuss ultraslow and lateral diffusion processes. Also, the applications of fractional equations with space distributed-order can be found in [6, 23, 33, 44].
Arrangement of fractional derivatives
Definitions of the Riesz fractional derivatives in the x- and y-directions are [24, 36]
∂|x|β ∂ β u(x, y, t) ∂|y|β
=
−1 RL β RL β D u(x, y, t) + D u(x, y, t) , ℘ (y) x x ℘2 (y) 2 cos( βπ ) 1 2
=
−1 RL β RL β D u(x, y, t) + D u(x, y, t) , = (x) y y =2 (x) ) 1 2 cos( βπ 2
in which [24, 36]
Pr e-
∂ β u(x, y, t)
d2 1 = Γ(2 − β) dx2
Zx
d2 1 RL β D u(x, y, t) = x b Γ(2 − β) dx2
Zb
1 d2 = Γ(2 − β) dy 2
Zy
1 d2 RL β D u(x, y, t) = y d Γ(2 − β) dy 2
Zd
urn
al
RL β a Dx u(x, y, t)
RL β c Dy u(x, y, t)
a
x
c
y
u(s, y, t)
(x − s)β−1 u(s, y, t) (s − x)β−1 u(x, s, t) (y − s)β−1 u(x, s, t) (s − y)β−1
ds,
ds,
ds,
ds.
Jo
Lemma 2.1. [24, 36] Let v ∈ JS% (Ω) and let v¯ be the deployment of v outside of Ω then RL % RL a Dx v,x
RL % RL c Dy v,y
D%b v
D%d v
L2 (Ω)
L2 (Ω)
= =
RL % RL ¯,x a Dx v RL % RL ¯,y c DL v
% 2 D%b v¯ L2 (R2 ) = cos(π%) RL ¯ L2 (R2 ) , −∞ Dx v
% 2 D%d v¯ L2 (R2 ) = cos(π%) RL ¯ L2 (R2 ) . −∞ Dy v 3
Journal Pre-proof
Lemma 2.2. [48, 49] Let β ∈ (0, 1) then the following second-order weighted and shifted Gr¨ unwald difference (WSGD) approximation formula holds β 0 Dt u(x, t)
=τ
n X
−β
λβ (k)u(x, tn−k ) + O(τ 2 ),
(2.1)
λβ (k) =
and ω0β
= 1,
ωkβ
β+2 β 2 ω0 ,
k = 0,
α + 2 ωβ − β ωβ , k 2 2 k−1
Γ(k − β) = , Γ(−β)Γ(k + 1)
ωkβ
k > 0,
p ro
in which
of
k=0
β+1 β = 1− ωk−1 , k
k ≥ 1.
Lemma 2.3. [48, 49, 51] According to relation (2.1), for each N ∈ N and real vector u0 , u1 , . . . , uN ∈ RN +1 , the following inequalities hold
Pr e-
N X n X n=0 k=0
λβ (k) un−k , un ≥ 0,
n X
λβ (k) < 0.
k=0
Lemma 2.4. [47] Assuming f (t) ∈ C 1 ([0, T ]) ∩ C 3 ((0, T )) and 1
(tn − ξ)− 2 f (ξ)dξ.
al
I(f, tn ) =
Ztn 0
urn
Then
I(f, tn ) = an f (t0 ) +
an = 2
Jo
where
b0 =
1 tn2
−
t Z1 2
τ
0
1 τ
n X r=0
tZn+1
tn
1 y 2 dy
3
br f (tn−r ) + O(τ 2 ),
1 y 2 dy ,
√ 4 τ + β, 3 4
n = 1, 2, . . . ,
Journal Pre-proof
b1 =
t Z2 2
τ
1 y 2 dy
−
t0
t1
br
Zt1
√ 4 τ − β, 3
1 y 2 dy
t Zr+1 Ztr 1 1 2 = y 2 dy , y 2 dy − τ tr
r ≥ 2.
tr−1
sin(πθ1 ) + sin(πθ2 ) . sin(πα)
We set the operators
p ro
C(α, θ1 , θ2 ) =
of
Definition 2.1. [34] Let 1 < α < 1, 0 < θ1 , θ2 < α, θ1 + θ2 = α, 0 ≤ p ≤ 1 and
1 C(α, θ1 , θ2 ) (I0α+ + I1α− ) , 2
θ1 ,θ2 ,α D−1,1 : =
dk θ1 ,θ2 ,k−α I , dxk −1,1
Pr e-
θ1 ,θ2 ,α I−1,1 : =
Jm−θ,−s (x) = (1 − x)θ (1 + x)s Pmθ,s ,
θ, s > −1.
Lemma 2.5. [34] Let 1 < p < 1. If 0 < θ1 , θ2 < p such that θ1 + θ2 = α and p sin(πθ1 ) = (1 − p) sin(πθ2 ),
so for x ∈ [−1, 1], n = 0, 1, 2, . . . and k = 1, 2, . . . , n + 2 we have
al
θ1 ,θ2 ,2−α −θ1 ,−θ2 θ2 −2,θ1 −2 b α)Pn+2 I−1,1 Jn (x) = C(n, ,
urn
θ1 ,θ2 ,k−2+α −θ1 ,−θ2 e k, α)P θ2 −2+k,θ1 −2+k , D−1,1 Jn (x) = C(n, n+2−k
Also, for l = −1, 0, 1, . . . , min {m, n} we have =
θ2 ,θ1 Kn,l δm,n
=
b k, α) = Γ(n + k + α − 1) . C(n, 2k−2 n!
Z1
Dpα+κ Jm−θ1 ,−θ2 (x)Dpα+κ Jn−θ1 ,−θ2 (x)ω θ1 +l,θ2 +l (x)dx,
Z1
α+κ −θ2 ,−θ1 α+κ −θ2 ,−θ1 D1−p Jm (x)D1−p Jn (x)ω θ1 +l,θ2 +l (x)dx,
−1
Jo
θ1 ,θ2 Kn,l δm,n
b α) = 4Γ(n + α − 1) , C(n, n!
−1
θ2 ,θ1 Kn,l = 2α+1
Γ(n + θ1 + 1)Γ(n + θ2 + 1)Γ(n + l + α + 1) . (2n + α + 1)(n − l)!(n!)2 5
Journal Pre-proof
3
Analytical investigation
In the next, the order of convergence and the stability issue of the new developed methodology have been investigated. For 0 < γ < 1 we define the following functional spaces n o 1+2γ H 2γ (Ω) = u ∈ L2 (Ω) s.t. (1 + ω 2 ) 2 u b ∈ L2 (Ω) , H02γ (Ω) = u ∈ H 2γ (Ω) s.t. u = 0
on ∂Ω .
Study the semi-discrete formulation
p ro
3.1
of
in which u b is the Fourier transformation of u and
At first, we have $(β)
∂ β u(x, y, t) β
∂|x|
1
1
dβ =
MG X
(1)
wj $(βj )
j=0
∂ βj u(x, y, t) βj
∂|x|
Pr e-
Z2
+ O((MG1 )−r ),
Z2
$(β)
Z1
MG X ∂ αj u(x, y, t) ∂ α u(x, y, t) (3) ω(α ) ω(α) dα = w + O((MG2 )−r ), j j αj ∂tα ∂t j=0
∂ β u(x, y, t)
1
0
∂|y|β
1
dβ =
MG X
(1)
wj $(βj )
∂ βj u(x, y, t)
j=0
∂|y|βj
+ O((MG1 )−r ),
2
2
k=0 j=0
(3)
λαj (k)wj ω(αj )u(x, y, tn−k ) + u(x, y, tn )
urn
MG n X X
Jo
τ
−α
al
in which αj and βj are the Legendre-Gauss points on the intervals (0, 1) and (1, 2), respectively. By employing Lemmas 2.2 and 2.4, the time-discrete of Eq. (1.1) is as follows [17]
1
= K
MG X j=0
(1)
wj $(βj )
(3.1) (
∂ βj u(x, y, tn ) ∂|x|βj
+
∂ βj u(x, y, tn ) ∂|y|βj
∂ 2 u(x, y, 0) ∂ 2 u(x, y, 0) + an + ∂x2 ∂y 2 2 n X ∂ u(x, y, tn−r ) ∂ 2 u(x, y, tn−r ) + br + ∂x2 ∂y 2 r=0 + f (x, y, tn ) + O((MG1 )−r ) + O((MG2 )−r ).
6
)
Journal Pre-proof
By omitting the small term, we can get 2
1
(3) λαj (k)wj ω(αj )U n−k
+ Un = K
MG X
(1) wj $(βj )
(
∂ βj U n
∂ βj U n
)
(3.2) ∂|y|βj 2 X 2 n−r n ∂ ϕ ∂ 2ϕ ∂ U ∂ 2 U n−r + an br + f n. + 2 + + 2 2 2 ∂x ∂y ∂x ∂y r=0
k=0 j=0
j=0
+
∂|x|βj
β
of
τ −α
MG n X X
Also, the weak form of Eq. (3.2) is: Find a U n ∈ H02 (Ω) such that τ
MG n X X k=0 j=0
(3) λαj (k)wj ω(αj ) U n−k , v
1
n
+ hU , vi = K
MG X j=0
n
Bj (U , v) = Kj + Kj
βj
RL 2 ℘1 (y) Dx
βj
RL 2 =1 (x) Dy
and
n X r=0
D℘2 (y) v +
U n ,RL y
βj 2
D=2 (x) v + K
β
2 cos( 2j π)
al
Kj =
βj 2
U n ,RL x
(3.3)
br ∇U n−r , ∇v + hf n , vi , ∀v ∈ H0 (Ω),
Pr e-
+ an h∆ϕ, vi − where
(1)
wj $(βj )Bj (U n , v)
p ro
2
−α
βj
RL n RL 2 x D℘2 (y) U ,℘1 (y)
βj
β 2
βj 2
Dx v
RL n RL 2 y D=2 (x) U ,=1 (x)
βj 2
Dy v
.
Lemma 3.1. [30] For an in Eq. (2.2), we have m X n=0
m X
√ |an | ≤ 4τ T ,
urn
τ
n=0
1
an ≤ Cτ 2 .
Jo
Proof. According to definition an , we have tZn+1 tZn+1 m m m X X X 1 1 1 1 1 tn2 − tn2 + 1 τ |an | = 2τ y 2 dy ≤ 2τ y 2 dy τ τ n=0 n=0 n=0 ≈ 2τ
m h X n=0
1 2
1 2
i
tn + y1 ≤ 2τ
tn
m X n=0
−1 tn 2
tn
≤ 2τ
Zmτ 0
7
√ dx √ ≤ 4τ T , x
,
Journal Pre-proof
in which y1 ∈ (tn , tn+1 ). In the current paper, we employ the following inequalities to prove the theorems.
ξ0 ≤ ε0 ,
n−1
of
Proposition 3.1. (Discrete Gronwall inequality) [39] Assume that $n is a non-negative sequence and that the sequence ξn satisfies
n−1
X X νk + $k ξk , ξn ≤ ε0 + k=0
Then, if ε0 ≥ 0 and ν0 ≥ 0, it follows ε0 +
k=0
νk
!
exp
n−1 X
$k
k=0
!
n ≥ 1.
,
Pr e-
ξn ≤
n−1 X
p ro
k=0
n ≥ 1.
Lemma 3.2. (Gronwall’s inequality) [50] Suppose that the discrete function {wn | n = 0, 1, 2, . . . , N ; N τ = T } , satisfies the inequality
wn ≤ γ + τ
θl wl ,
l=1
al
where γ and θl are nonnegative constants. Then
n X
kwn kL2 (Ω) ≤ γ exp 2τ
urn
β
n X l=1
θl
!
.
Theorem 3.2. Let U n ∈ H02 (Ω). The formulation (3.1) is unconditionally stable. e n be the approximate solution of (3.3), then the roundoff error equation is as follows Proof. Let U 2
MG n X X k=0 j=0
(3) λαj (k)wj ω(αj ) E n−k , v
Jo
τ
−α
1
n
+ hE , vi =
MG X
(1)
wj $(βj )Bj (E n , v)
j=0
n β
0 X
+ an ∆E , v − br ∇E n−r , ∇v , ∀v ∈ H02 (Ω), r=0
8
Journal Pre-proof
e n and U e n is approximation solution of Eq. (3.4). If we set v = E n , then in which E n = U n − U τ −α
k=0 j=0
(3) λαj (k)wj ω(αj ) E n−k , E n + hE n , E n i 1
=
MG X j=0
(1) wj $(βj )Bj (E n , E n )
n
0 n X br ∇E n−r , ∇E n . + an ∆E , E − r=0
of
2
MG n X X
p ro
Also from [8] we obtain βj βj βj βj RL RL n n n n n RL n RL 2 2 2 2 Bj (E , E ) = |2Kj | ℘1 (y) Dx E ,x D℘2 (y) E + |2Kj | =1 (x) Dy E ,y D=2 (x) E
βj βj βj βj RL RL n RL n n RL n 2 2 2 2 ≥ min{|2Kj | , |2Kj |} ℘1 (y) Dx E ,x D℘2 (y) E + =1 (x) Dy E ,y D=2 (x) E
Pr e-
≥ Kj kE n k2H β (Ω) .
(3.4)
Summing Eq. (3.4) for n from 0 to N , gives 2
τ −α
MG N X n X X n=0 k=0 j=0
n=0
= +
N X
Using Lemma 2.1, we have
n=0 k=0 j=0
n=0 j=0
an ∆E 0 , E n −
N X n X n=0 r=0
(3.5)
br ∇E n−r , ∇E n .
M2
N X n G X
X
(3) (3) λαj (k)wj ω(αj ) E n−k , E n = wj ω(αj ) λαj (k) E n−k , E n ≥ 0.
urn
2
MG N X n X X
al
n=0
τ −α
1
MG N X N X X
n−k n (1) (3) n n hE , E i + wj $(βj )Bj (E n , E n ) λαj (k)wj ω(αj ) E , E +
n=0 k=0
j=0
Also, form Eq. (3.4), we can write
n=0 j=0
and
1
(1) wj $(βj )Bj (E n , E n )
Jo
1
MG N X X
=
MG X
(1) wj $(βj )
j=0
−
N X n X n=0 r=0
N X n=0
1
n
n
Bj (E , E ) ≥
MG X
br ∇E n−r , ∇E n ≤ 0.
9
j=0
(1) wj $(βj )Kj
N X n=0
kE n k2H β (Ω) , (3.6) (3.7)
Journal Pre-proof
Thus, using the above inequalities, Eq. (3.5) can be rewritten as 1
MG X
(1) wj $(βj )Kj
N X n=0
j=0
kE n k2H β (Ω)
+
N X n=0
n
N X
n
hE , E i =
n=0
an ∆E 0 , E n .
(3.8)
By changing the index, we will have
j=0
1
≤
MG X
i=0
i=0
It can be concluded that kE n k2L2 (Ω)
n n n X X
i 2
i i X
0 i
E β + E , E = a ∆E , E . i H (Ω)
(1)
wj $(βj )Kj
(3.9)
i=0
n n n X X
j 2
j j X
E β = aj ∆E 0 , E j + E , E H (Ω) j=0
j=0
j=0
of
(1)
wj $(βj )Kj
p ro
1
MG X
j=0
0
+
n X j=1
Pr e-
≤ a0 ∆E , E
0
aj ∆E 0 , E j
n X
0 2
0
≤ |a0 | ∇E L2 (Ω) + ∆E L2 (Ω) aj E j L2 (Ω)
Thus, we get the following inequality
j=1
n X
0 2
0
≤ 1 + |a0 | ∇E L2 (Ω) + ∆E L2 (Ω) aj E j L2 (Ω) .
Jo
urn
kE kL2 (Ω)
n X
2
≤ 1 + |a0 | ∇E 0 L2 (Ω) + ∆E 0 L2 (Ω) aj E j L2 (Ω) .
al
n
j=1
10
j=1
(3.10)
Journal Pre-proof
j=1
aj ≤
n X
(aj )2 , we can write Eq. (3.10) as follows
j=1
n
kE kL2 (Ω)
n X
0 2
0
1
≤ 1 + |a0 | ∇E L2 (Ω) + ∆E L2 (Ω) τ (aj )2 E j L2 (Ω) τ j=1 ! n X
0 2
0 1 (aj )2 ≤ 1 + |a0 | ∇E L2 (Ω) exp τ ∆E L2 (Ω) τ j=1 !2 n X
0 2 ≤ 1 + |a0 | ∇E L2 (Ω) exp C aj
≤
and this completes the proof.
p ro
j=1
of
Since
n X
2 1 + |a0 | ∇E 0 L2 (Ω) exp (Cτ ) , β
Pr e-
Theorem 3.3. Let un , U n ∈ H02 (Ω) then relation (3.1) is convergent.
Proof. The exact and approximate solutions can be found from the following relations 2
τ −α
MG n X X
1
(3) λαj (k)wj ω(αj )un−k
k=0 j=0
+ un = K
MG X
(1) wj $(βj )
(
∂ βj un
+
∂ βj un
)
(3.11) ∂|y|βj X 2 n−r 2 n ∂ 2 un−r ∂ u ∂ ϕ ∂ 2ϕ + 2 + br + + an + f n + Rτ , 2 2 2 ∂x ∂y ∂x ∂y r=0 j=0
∂|x|βj
3
al
in which there exists a positive constant C such that |Rτ | ≤ Cτ 2 and 1
2
MG n X X
(3) λαj (k)wj ω(αj )U n−k
+U
n
= K
k=0 j=0
MG X
(1) wj $(βj )
(
∂ βj U n
j=0
∂|x|βj
+
∂ βj U n
)
(3.12) ∂|y|βj 2 X 2 n−r n ∂ ϕ ∂ 2ϕ ∂ U ∂ 2 U n−r + an + f n. + 2 + br + 2 2 ∂x2 ∂y ∂x ∂y r=0
urn
τ
−α
Jo
We set Qn = un − U n then subtracting (3.12) from (3.11) gives 2
τ −α
MG n X X k=0 j=0
(3) λαj (k)wj ω(αj )Qn−k
1
+ Qn = K
+
MG X j=0
n X r=0
11
(1) wj $(βj )
br
2
(
n−r
∂ Q ∂x2
∂ βj Qn ∂|x|βj 2
+
+
n−r
∂ Q ∂y 2
∂ βj Qn
∂|y|βj
) (3.13)
+ Rτ .
Journal Pre-proof
Thus, the weak form of Eq. (3.13) is τ −α
k=0 j=0
(3) λαj (k)wj ω(αj ) Qn−k , ξ =−
n X r=0
1
+ hQn , ξi + K
MG X j=0
br ∇Qn−r , ∇ξ + hRτ , ξi ,
β
ξ ∈ H02 (Ω).
If we set ξ = Qn , then τ −α
k=0 j=0
(3) λαj (k)wj ω(αj ) Qn−k , Qn =−
n X r=0
1
n
n
+ hQ , Q i + K
n=0 k=0 j=0
=−
N X n X n=0 r=0
1
+
1
(1)
N X n=0
n
1
n
hQ , Q i + K
MG N X X n=0 j=0
(1)
wj $(βj )Bj (Qn , Qn )
n=0
wj $(βj )Bj (E n , Qn ) =
MG X
(1)
wj $(βj )
j=0
urn
n=0 j=0
(3.14)
N
X n−r n br ∇Q , ∇Q + hRτ , Qn i.
Now we have MG N X X
j=0
Pr e-
τ −α
(3) λαj (k)wj ω(αj ) Qn−k , Qn
al
2
(1)
wj $(βj )Bj (Qn , Qn )
br ∇Qn−r , ∇Qn + hRτ , Qn i .
Summing Eq. (3.14) for n = 0 : N , yields MG N X n X X
MG X
p ro
2
MG n X X
(1)
wj $(βj )Bj (Qn , ξ)
of
2
MG n X X
N X n=0
(3.15) 1
Bj (Qn , Qn ) ≥
MG X
(1)
wj $(βj )Kj
j=0
N X n=0
kQn k2H β (Ω) .
Also, using Lemma 2.3, concludes 2
MG N X n X X n=0 k=0 j=0
(3) λαj (k)wj ω(αj ) Qn−k , E n
Jo
τ
−α
−
N X n X n=0 r=0
2
=
MG X j=0
(3) wj ω(αj )
N X n X n=0 k=0
br ∇Qn−r , ∇Qn ≤ 0,
12
λα (k) Qn−k , Qn ≥ 0,
Journal Pre-proof
then, Eq. (3.15) can be rewritten as follows 1
≤ τ
i=0
i=0
i=0
i=0
n X
i
i 3
Q 2 2 kRτ kL2 (Ω) Q L2 (Ω) = Cτ L (Ω)
"
i=0
n X
1 2
of
# n n X X
i
i 2 1
Q 2 = Cτ τ Q L2 (Ω) ≤ Cτ τ+ L (Ω) 2 i=0 i=0 i=0 " # n n X X
i 2
i 2 1
Q 2 , Cτ T + ≤ Q L2 (Ω) ≤ Cτ T + Cτ L (Ω) 2 i=0 i=0
kQn k2L2 (Ω)
n X
i 2
Q 2 ≤ Cτ T exp (Cτ (n + 1)) ≤ C∗ τ, ≤ Cτ T + Cτ L (Ω) i=0
in which this totalizes the proof.
3.2
n n n X X
i 2
i i X
Q β + Rτ , Qi Q ,Q = H (Ω)
p ro
thus, we obtain
j=0
n X
≤
(1) wj $(βj )Kj
Pr e-
kQn k2L2 (Ω)
MG X
Analysis of full-discrete scheme
Let weight function ω(x) be positive, thus the weighted Sobolev space can be defined as Z 2 2 Lω (Ω) = u : u (x)ωdx < +∞ ,
with corresponding norm
al
Ω
urn
kukL2ω (Ω) =
Z Ω
21
u2 (x)ωdx .
For known values p, 1 < α < 2, 0 < µ, ν < α, l = −1, 0, 1, . . . , m and m ∈ N, we define m Bα,p,µ,ν = u ∈ L2ω−µ,−ν (Ω) :
Jo
Let Phm be the interpolation operator and
Dpα+1 u ∈ L2ων+l,µ+l (Ω) f or − 1 ≤ l ≤ m .
F−µ,−ν (Ω) = {φ = (1 − x)µ (1 + x)ν ϕ : ϕ ∈ PN } , N
be the finite dimensional fractional polynomial space. m Lemma 3.3. [34] Let 1 < α < 2 and u ∈ Bα,p,µ,ν . For for values p, 1 < α < 2, 0 < µ, ν < α and
13
Journal Pre-proof
l = −1, 0, 1, . . . , m, we have
α+l −µ,−ν
Dp PN u − u L2
ω ν+l,µ+l
−µ,−ν
P u − u L2 N
ω −µ,−ν
−µ,−ν
P u − u L2 N
ω −µ,−ν
(Ω)
(Ω)
≤ cN
l−m 2
s
(N − m + 1)!
Dpα+l u 2 , L ν+m,µ+m (Ω) (N − l + 1)! ω
α+l
Dp u 2 L
ω ν+m,µ+m
(Ω)
s
,
(for fixed value m),
(N − m + 1)!
Dpα+m u 2 , 0 ≤ m ≤ N, L ν+m,µ+m (Ω) (N + m + 1)! ω s
(N − m + 1)!
Dpα+m u 2 ≤ cN −α−m , (fixed value m). L ν+m,µ+m (Ω) (N + m + 1)! ω ≤ cN −α
(Ω)
(Ω)
of
ω ν+l,µ+l
≤ cN
l−m 2
p ro
α+l −µ,−ν
Dp PN u − u L2
We define the projection operator
Πh : H0β (Ω) ∩ C 0 (Ω) → Xh ,
Pr e-
such that
B(Πh z, w) = B(z, w),
∀ w ∈ Xh .
β
Theorem 3.4. Assume un = u(x, y, tn ) ∈ H02 (Ω) and Uhn (x, y, tn ) ∈ Xh are the exact and approximate solutions of (1.1), respectively. Thus 3 kun − Uhn k ≤ C ∗ τ 2 + N −β−m . Proof. We can write τ
MG n X X
1
(3) λαj (k)wj ω(αj )un−k
2
τ −α
MG n X X k=0 j=0
(3) λαj (k)wj ω(αj )UNn−k
Jo
and
+U
n
= K
(1) wj $(βj )
(
∂ βj un
j=0
+ UNn = K
1
MG X
+ an
(1) wj $(βj )
j=0
∂ 2ϕ ∂ 2ϕ + 2 ∂x2 ∂y
14
∂|x|βj
(
+
∂ βj un
)
(3.16) ∂|y|βj 2 X 2 n−r n ∂ ϕ ∂ 2ϕ ∂ 2 un−r ∂ u + an + 2 + br + + f n + Rτ , 2 2 ∂x2 ∂y ∂x ∂y r=0
urn
k=0 j=0
MG X
al
2
−α
∂ βj UNn ∂|x|βj
+
n X r=0
+
br
∂ βj UNn
∂|y|βj
)
∂ 2 UNn−r ∂ 2 UNn−r + ∂x2 ∂y 2
(3.17)
+ f n.
Journal Pre-proof
Subtracting (3.17) from (3.16) yields τ −α
k=0 j=0
(3)
λαj (k)wj ω(αj ) un−k − UN n−k
+ (un − UN n ) 1
= K
MG X
(1) wj $(βj )
j=0
n X
+
br
Thus, the weak form of Eq. (3.18) is 2
τ
MG n X X k=0 j=0
j=0
∂|x|βj n−r
+
∂ βj (un − UN n )
∂ 2 un−r − UN ∂x2
(1)
wj $(βj )Bj (un − UN n , ϑ) = −
n X r=0
br ∇ un−r − UN n−r , ∇ϑ + hRτ , ϑi .
The use of this equation and the notations
ZNn = PN un − UNn , yield τ −α
k=0 j=0 2
(3) λαj (k)wj ω(αj ) XNn−k , ϑ + hXNn , ϑi
urn
MG n X X
(3) λαj (k)wj ω(αj ) ZNn−k , ϑ + hZNn , ϑi
k=0 j=0
1
+K
MG X
(1) wj $(βj )Bj
(ZNn , ϑ)
j=0
Jo
+τ
−α
XNn = un − PN un ,
al
2
MG n X X
15
= −
n X r=0
∂|y|βj
∂ 2 un−r − UN n−r + ∂y 2
Pr e-
+K
∂ βj (un − UN n )
(3) λαj (k)wj ω(αj ) un−k − UN n−k , ϑ + hun − UN n , ϑi
1
MG X
(
p ro
r=0
−α
(3.18)
of
2
MG n X X
br ∇ZNn−r , ∇ϑ + hRτ , ϑi .
!
)
+ Rτ .
Journal Pre-proof
We set ϑ = ZNn , then τ −α
k=0 j=0 2
+τ −α
MG n X X k=0 j=0
(3) λαj (k)wj ω(αj ) ZNn−k , ZNn + hZNn , ZNn i
(3) λαj (k)wj ω(αj ) XNn−k , ZNn + hXNn , ZNn i 1
(1) wj $(βj )Bj
j=0
(ZNn , ZNn )
= −
Summing Eq. (3.19) for n = 0 : N , confirms that τ
MG K X n X X n=0 k=0 j=0 2
+τ
−α
MG N X n X X n=0 k=0 j=0
(3) λαj (k)wj ω(αj ) ZNn−k , ZNn
MG N X X
n=0 k=0 j=0
kZNn k2L2 (Ω)
(3.20)
(1)
N X n X n=0 r=0
N
X hRτ , ZNn i. br ∇ZNn−r , ∇ZNn + n=0
(3) λαj (k)wj ω(αj ) ZNn−k , ZNn ≥ 0,
(3.21)
al
τ −α
n=0
wj $(βj )Bj (ZNn , ZNn ) = −
Now, since 2
+
n=0
n=0 j=0
MG K X n X X
N X
N X
(3) λαj (k)wj ω(αj ) XNn−k , ZNn + hXNn , ZNn i 1
+K
r=0
br ∇ZNn−r , ∇ZNn + hRτ , ZNn i . (3.19)
Pr e-
2
−α
n X
p ro
+K
MG X
of
2
MG n X X
−
urn
ZN0 = 0,
N X n X n=0 r=0 1
MG N X X
(1) wj $(βj )Bj (ZNn , ZNn )
Jo
n=0 j=0
br ∇ZNn−r , ∇ZNn
(3.22)
≤ 0,
(3.23) 1
=
MG X j=0
(1) wj $(βj )
N X n=0
1
≥
MG X j=0
(1) wj $(βj )Kj
Bj (ZNn , ZNn )
N X n=0
kZNn k2H β (Ω) ,
! 2 2 MG N n MG N n X X 1 −α X X X 1 −α X (3) (3) n 2 n 2 τ λαj (k)wj ω(αj ) kZN kL2 (Ω) = τ wj ω(αj ) kZN kL2 (Ω) λαj (k) < 0, 2 2 n=0 k=0 j=0 n=0 j=0 k=0 16
Journal Pre-proof
Eq. (3.20) can be rewritten as follows N X n=0
1
kZNn k2L2 (Ω)
+
MG X
N X
(1) wj $(βj )Kj
n=0
j=0
kZNn k2H β (Ω)
2
of
N n MG
2 1 −α X X X (3) ≤ τ λαj (k)wj ω(αj ) XNn−k L2 (Ω) 2 n=0 k=0 j=0
N N N N X X 1X n 2 1X n 2 2 n 2 + kRτ kL2 (Ω) + kX k 2 + kZ k 2 + kZ k 2 , 4 n=0 N L (Ω) n=0 N L (Ω) n=0 4 n=0 N L (Ω)
p ro
or 1
MG n n X X X
i 2
i 2 (1)
ZN 2
ZN β + 2 w $(β )K j j j L (Ω) H (Ω) i=0
j=0
i=0
2
MG n X i X X
+ 2
2 (3) λαj (k)wj ω(αj ) XNi−k L2 (Ω)
Pr e-
≤ τ −α
i=0 k=0 j=0
n n X X
k 2
XN 2 + 2 kRτ k2L2 (Ω) . L (Ω) k=0
k=0
Changing the indices of summations arrives at
1
MG n n X X X
i 2
i 2 (1)
ZN β
ZN L2 (Ω) + 2 ≤ wj $(βj )Kj H (Ω) i=0
i=0
j=0
al
kZNn k2L2 (Ω)
2
MG n X i X X i=0 k=0 j=0
urn
≤ τ
−α
+ 2
2 (3) λαj (k)wj ω(αj ) XNi−k L2 (Ω)
n n X X
k 2
XN 2 + 2 kRτ k2L2 (Ω) L (Ω) k=0
k=0
M2
Jo
≤ τ −α C2 N + 2C2
n X k=0
n X i X G X −β−m 2
N
(3)
λαj (k)wj ω(αj )
i=0 k=0 j=0
−β−m 2
17
+ 2C1
n X k=0
9
τ 4.
Journal Pre-proof
Finally, we obtain N
k=0
−β−m 2
≤ 2τ (C1 + C2 ) ≤ 2 (C1 + C2 ) τ
k=0
3
τ 2 + N −β−m
k=0
n X
9
τ4
3
2
τ 2 + N −β−m
2
k=0
3 ≤ C exp (T ) τ + N −β−m ≤ C τ 2 + N −β−m ,
which completes the proof.
4
n X
+ 2τ C1
n X
of
≤ 2τ C2
n X
3 2
p ro
τ
kZNn k2L2 (Ω)
Numerical investigation
Pr e-
To verify the theoretical results, we apply the new numerical procedure for two examples. The used computational order is E1 E
log10 2
Rate =
h1 h2
log10
,
such that E1 and E2 are the absolute errors related to h1 and h2 mesh sizes, respectively.
4.1
Example 1.
ω(α)0 Dαt u(x, y, t)dα
0
+ u(x, y, t) =
Z2
$(β)
∂ β u(x, y, t)
urn
Z1
al
For the first model, we consider the following example
1
+
Zt
(t − ζ)
+
∂|y|β
∂ 2 u(x, y, ζ) ∂ 2 u(x, y, ζ) dζ + f (x, y, t), + ∂x2 ∂y 2
Jo
0
− 21
∂|x|β
∂ β u(x, y, t)
18
!
dβ
(x, y, t) ∈ Ω × (0, T ],
Journal Pre-proof
in which ω = Γ(2 − α) and $(β) = −2Γ(5 − β) cos f (x, y, t) = 10x2 (1 − x2 )y 2 (1 − y 2 )
(t − 1) ln(t)
πβ 2
and also
− 10ty 2 (1 − y 2 ) [φ1 (x) + φ1 (1 − x) + φ2 (x) + φ2 (1 − x)] 3 40 (1 − x2 )y 2 (1 − y 2 )t 2 , 3
of
− 10tx2 (1 − x2 ) [φ1 (y) + φ1 (1 − y) + φ2 (y) + φ2 (1 − y)] −
Pr e-
p ro
2 ϑ3 − ϑ2 3ϑ − 2ϑ ϑ2 − ϑ φ1 (ϑ) = Γ (5) − 2Γ (4) − ln ϑ ln ϑ (ln ϑ)2 Γ (3) 2 5ϑ 3 2ϑ + − , 6ϑ − 2 − + + ln ϑ ln ϑ ln ϑ (ln ϑ)2 (ln ϑ)2 3 ϑ4 − ϑ3 4ϑ − 3ϑ2 ϑ3 − ϑ2 φ2 (ϑ) = Γ (5) − 2Γ (4) − ln ϑ ln ϑ (ln ϑ)2 Γ (3) 2ϑ 1 2ϑ2 2 2 + + 12ϑ − 6ϑ − 7ϑ − 5ϑ − , ln ϑ ln ϑ ln ϑ ln ϑ with the exact solution
Jo
urn
al
u(x, y, t) = 10tx2 (1 − x2 )y 2 (1 − y 2 ).
Figure 1: Error obtained for Example 1.
19
Journal Pre-proof
Table 1 Errors and computational orders at T = 0.5 for Example 1 N = 10 Rate
3.7826 × 10
1/10
7.6644 × 10
1/20
1.8202 × 10
1/40
−4
−
−5
2.0022
6.9621 × 10−8
2.0009
2.7865 × 10
1/320
Rate
CPU time
−
5
−6
2.0664
25
−7
2.0234
46
−8
2.0092
77
−8
2.0039
93
8.2710 × 10
2.0621 × 10
2.0064
−7
−6
3.3295 × 10
2.0208
1.1163 × 10
1/640
1.3536 × 10
2.0741
−6
−6
1/160
2.6694 × 10
2.3031
−5
4.4853 × 10
1/80
L∞
of
L∞
−9
2.0018
126
1.2865 × 10−9
2.0008
178
5.1489 × 10
p ro
τ
N = 20
4.2
Pr e-
Table 1 is based on the N = 10, N = 20 and 16 Gauss–Legendre-Lobatto points to integration. In Table 1, we can observe that the theoretical results are complied to the computational conclusions. Figure 1 illustrates the error obtained as a function based upon the different values of N and τ = 10−4 for left panel and different values of time step and N = 20 in the right panel for Example 1.
Example 2.
For the second example, we consider a test problem without exact solution Z1
ω(α)0 Dαt u(x, y, t)dα + u(x, y, t) =
0
Z2
$(β)
1
− 21
(t − ζ)
with initial data
∂ 2 u(x, y, ζ) ∂ 2 u(x, y, ζ) + dζ, ∂x2 ∂y 2
urn
0
∂|x|β
al
+
Zt
∂ β u(x, y, t)
u(x, y, 0) = exp (−10 [(x − 0.5) + (y − 0.5)]) ,
+
∂ β u(x, y, t) ∂|y|β
!
dβ
(x, y, t) ∈ Ω × (0, T ],
(x, y) ∈ Ω,
and homogeneous boundary condition
and also
Jo
u(x, y, t) = 0,
(x, y) ∈ ∂Ω,
ω(α) = Γ (2 − α) ,
0 ≤ t ≤ 1,
$(β) = exp (−β) .
To check the convergence of the present numerical technique, we adopt the following strategy: • Select a N ∗ and τ small enough and obtain the numerical solution UNτ ∗ using two parameters. 20
Journal Pre-proof
• Set N1 , N2 , . . . in which N ∗ > N i for i = 1, 2, . . . and obtain the related numerical solution UNτ 1 , UNτ 2 , . . . based on the N1 , N2 , . . .. • Put
∗
LN ∞ =
max |Uhτ∗ (xj ) − Uhτi (xj )| ,
i = 1, 2, . . . .
1≤j≤M −1
of
Note that the dimensions of UNτ ∗ and UNτ i are different but we consider the common nodes in these solutions. Table 2
Errors and computational orders at T = 0.5 for Example 2 τ = 10−5
∗
∗
p ro
τ = 10−3 N
LN ∞
Rate
LN ∞
5
1.0378 × 10−3
−
1.1754 × 10−6
1.9749 × 10
1.2057 × 10
−8
7.6627
1.0694 × 10−9
8.4209
2.6951 × 10
20 30 40
5.7156
−7
6.1953
CPU time
−
11
2.2442 × 10
−6
5.7108
71
−8
7.7133
131
3.2966 × 10
−10
8.5812
241
11.6660
388
1.0694 × 10
1.1496 × 10−11
Pr e-
10
−5
Rate
Table 3
Errors and computational orders at T = 0.5 for Example 2 N ∗ = 40, τ = 10−4 τ
L∞
1/10
1.8781 × 10−4
1/40 1/80 1/160
L∞
−
1.6384 × 10−6
1.9921
−5
1.9985
2.9543 × 10−6
1.9997
1.1814 × 10 7.3862 × 10 1.8466 × 10
−7 −7
1.9999 2.0000
urn
1/320
4.7210 × 10
Rate
al
1/20
−5
N ∗ = 80, τ = 10−4 Rate −
−7
4.4409 × 10
1.8833
−7
1.9711
2.8459 × 10−8
1.9928
−9
1.9981
−9
1.9991
1.1327 × 10 7.1242 × 10 1.7821 × 10
Jo
Table 2 presents the errors and computational orders obtained at T = 0.5 for Example 2. In Table 2, we assume that the reference solution is based on N ∗ = 80 and τ = 10−5 . In this table, two values are fixed for time step τ and then the number of basis N is increasing. To construct the Table 3, we consider two reference solutions such that one of them is based on N ∗ = 40 and τ = 10−4 and another is N ∗ = 80 and τ = 10−4 . Table 3 shows that the proposed numerical procedure has second-order accuracy in the temporal direction. It must be noted that an extrapolation approach is used to improve the numerical results. For this reason, computational order of temporal variable has been increased.
21
Journal Pre-proof
5
Conclusion
of
A high-order numerical scheme is proposed to solve the two-dimensional time-space distributedorder weakly singular integro-partial differential equation. The temporal direction is discretized by using a second-order difference formula. Furthermore, to get a numerical procedure with high accuracy, derivatives of spatial directions are approximated by a Galerkin spectral idea based upon the modified Jacobi polynomials. The stability, convergence and error estimation of the semi- and full-discrete plan have been studied in details. In the end, the theoretical results have been checked by using two test problems.
p ro
Acknowledgments
The authors are very grateful to both reviewers for carefully reading this paper and for their comments and suggestions which have improved the paper.
Pr e-
References
[1] M. Abbaszadeh, Error estimate of second-order finite difference scheme for solving the Riesz space distributed-order diffusion equation, Applied Mathematics Letters 88, (9)2019) 179-185. [2] M. Abbaszadeh, M. Dehghan, Meshless upwind local radial basis function-finite difference technique to simulate the time- fractional distributed-order advection–diffusion equation, Engineering with Computers, In press (2019).
al
[3] M. Abbaszadeh, M. Dehghan, Y. Zhou, Alternating direction implicit-spectral element method (ADI-SEM) for solving multi-dimensional generalized modified anomalous sub-diffusion equation, Computers & Mathematics with Applications, 78 (2019) 1772-1792.
urn
[4] A. A. Alikhanov, A new difference scheme for the time fractional diffusion equation, J. Comput. Phys., 280 (2015) 424-438. [5] A. A. Alikhanov, Numerical methods of solutions of boundary value problems for the multiterm variable-distributed order diffusion equation, Appl. Math. Comput., 268 (2015) 12-22.
Jo
[6] T. M. Atanackovic, S. Pilipovi, D. Zorica, Distributed-order fractional wave equation on a finite domain, stress relaxation in a rod, International Journal of Engineering Science, 49 (2011) 175-190. [7] A. H. Bhrawy, M. A. Zaky, An improved collocation method for multi-dimensional space-time variable-order fractional Schrodinger equations, Appl. Numer. Math., 111 (2017) 197-218. [8] W. Bu, Y. Tang, Y. Wu, J. Yang, Finite difference/finite element method for two-dimensional space and time fractional Bloch-Torrey equations, J. Comput. Phys. 293 (2015) 264-279. 22
Journal Pre-proof
[9] M. Caputo, Mean fractional-order-derivatives differential equations and filters, Annali dellUniversita di Ferrara, 41 (1995) 73–84. [10] M. Caputo, Distributed order differential equations modelling dielectric induction and diffusion, Fractional Calculus and Applied Analysis, 4 (2001) 421-442.
of
[11] M. Caputo, Diffusion with space memory modelled with distributed order space fractional differential equations, Annals of Geophysics, 46 (2003). [12] C. Celik, M. Duman, Crank-Nicolson method for the fractional diffusion equation with the Riesz fractional derivative, Journal of Computational Physics, 231 (2012) 1743-1750.
p ro
[13] A. Chechkin, R. Gorenflo, I. Sokolov, Retarding subdiffusion and accelerating superdiffusion governed by distributed-order fractional diffusion equations, Physical Review E, 66 (2002) 046129. [14] R. M. Christensen, Theory of Viscoelasticity, Academic Press, New York, 1971.
Pr e-
[15] M. Dehghan, M. Abbaszadeh, Error estimate of finite element/finite difference technique for solution of two-dimensional weakly singular integro-partial differential equation with space and time fractional derivatives, Journal of Computational and Applied Mathematics 356 (2019) 314-328. [16] M. Dehghan, Solution of a partial integro–differential equation arising from viscoelasticity, International Journal of Computer Mathematics 83 (1), (2006) 123–129.
al
[17] M. Dehghan, Finite difference procedures for solving a problem arising in modelling and design of certain optoelectronic devices, Mathematics and Computers in Simulation, 71 (1), (2006) 16-30.
urn
[18] M. Dehghan, J. Manafian, A. Saadatmandi, Solving nonlinear fractional partial differential equations using the homotopy analysis method, Numerical Methods for Partial Differential Equations, 26(2), (2010), 448-479. [19] W. Deng, J. S. Hesthaven, Local discontinuous Galerkin methods for fractional ordinary differential equations, BIT Numerical Mathematics 55 (4), (2015) 967-985.
Jo
[20] H. Ding, General Pade approximation method for time–space fractional diffusion equation, Journal of Computational and Applied Mathematics, 299 (2016) 221-228. [21] H. Ding, C. P. Li, Y. Q. Chen, High-order algorithms for Riesz derivative and their applications (II), Journal of Computational Physics, 293 (2015) 218-237. [22] H. Ding, C. P. Li, High-order numerical algorithms for Riesz derivatives via constructing new generating functions, Journal of Scientific Computing, 71 (2017) 759–784.
23
Journal Pre-proof
[23] C. H. Eab, S. C. Lim, Fractional Langevin equations of distributed order, Physical Review E., 83 (2011), 031136. [24] V. J. Ervin, J. P. Roop, Variational formulation for the stationary fractional advection dispersion equation, Numer. Methods Partial Differ. Equ. 22 (2006) 558-576.
of
[25] W. Fan, F. Liu, X. Jiang, I. Turner, A novel unstructured mesh finite element method for solving the time-space fractional wave equation on a two-dimensional irregular convex domain, Fractional Calculus and Applied Analysis 20.2 (2017) 352-383.
p ro
[26] G. H. Gao, A. A. Alikhanov, Z. Z. Sun, The temporal second order difference schemes based on the interpolation approximation for solving the time multi-term and distributed order fractional sub-diffusion equations, J. Sci. Comput., 73 (2017) 93-121. [27] A. N. Kochubei, Distributed order calculus and equations of ultraslow diffusion, Journal of Mathematical Analysis and Applications, 340 (2008) 252-281.
Pr e-
[28] J. Li, F. Liu, L. Feng, I. Turner, A novel finite volume method for the Riesz space distributedorder diffusion equation, Computers and Mathematics with Applications, 46 (2017) 536-553. [29] H. L. Liao, P. Lyu, S. Vong, Second-order BDF time approximation for Riesz space-fractional diffusion equations, Int. J. Comput. Math., 95 (2018) 144–158. [30] M. Luo, D. Xu, L. Li, A compact difference scheme for a partial integro-differential equation with a weakly singular kernel, Appl. Math. Model. 39 (2015) 947-954.
al
[31] J. E. Macias-Diaz, Numerical study of the process of nonlinear supratransmission in Riesz space-fractional sine-Gordon equations, Commun. NonLinear Sci. Numer. Simul. 46 (2017) 89-102.
urn
[32] J. E. Macias-Diaz, A structure-preserving method for a class of nonlinear dissipative wave equations with Riesz space-fractional derivatives, Journal of Computational Physics, 351 (2017) 40-58. [33] F. Mainardi, A. Mura, R. Gorenflo, M. Stojanovic, The two forms of fractional relaxation of distributed order, Journal of Vibration and Control, 13 (2007) 1249–1268.
Jo
[34] Z. Mao, G. Em Karniadakis, A spectral method (of exponential convergence) for singular solutions of the diffusion equation with general two–sided fractional derivative, SIAM J. Numer. Anal., 56(1) (2018) 24–49. [35] R. K. Miller, An integro-differential equation for grid heat conductors with memory, J. Math. Anal. Appl. 66 (1978) 313–332. [36] J. P. Roop, Variational solution of the fractional advection dispersion equation, PhD thesis, Clemson University, 2004. 24
Journal Pre-proof
[37] H. K. Pang, H. W. Sun, Fourth order finite difference schemes for time-space fractional subdiffusion equations, Comput. Math. Appl., 71 (2016) 1287-1302. [38] E. Pindza, K. M. Owolabi, Fourier spectral method for higher order space fractional reactiondiffusion equations, Communications in Nonlinear Science and Numerical Simulation 40 (2016) 112-128.
of
[39] A. Quarteroni, A. Valli, Numerical Approximation of Partial Differential Equations, SpringerVerlag, New York, 1997.
p ro
[40] M. Rcnardy, Mathematical analysis of viscoelastic flows, Ann. Rev. Fluid Mech., 21 (1989) 21–36. [41] A. Saadatmandi, M. Dehghan, A Legendre collocation method for fractional integro–differential equations, Journal of Vibration and Control, 17(13), (2011) 2050–2058
Pr e-
[42] S. Shen, F. Liu, V. Anh, I. Turner, J. Chen, A novel numerical approximation for the space fractional advection–dispersion equation, IMA Journal of Applied Mathematics, 79 (3), (2014) 431-444. [43] I. Sokolov, A. Chechkin, J. Klafter, Distributed-order fractional kinetics, Acta Physica Polonica B, 35 (2004) 1323-1341. [44] N. Su, P. N. Nelson, S. Connor, The distributed-order fractional diffusion-wave equation of groundwater flow: theory and application to pumping and slug tests, Journal of Hydrology, 529 (2015) 1262-1273.
al
[45] H. Sun, Z. Z. Sun, G. H. Gao, Some high order difference schemes for the space and time fractional Bloch-Torrey equations, Appl. Math. Comput. 281 (2016) 356-380.
urn
[46] Z. Z. Sun, X.N. Wu, A fully discrete difference scheme for a diffusion-wave system, Appl. Numer. Math. 56 (2006) 193-209. [47] T. Tang, A finite difference scheme for a partial integro-differential equations with a weakly singular kernel, Appl. Numer. Math. 11 (1993) 309–319. [48] W. Tian, H. Zhou, W. Deng, A class of second order difference approximations for solving space fractional diffusion equations, Mathematics of Computation 84 (294) (2015) 1703–1727.
Jo
[49] Z. Wang, S. Vong, Compact difference schemes for the modified anomalous fractional sub– diffusion equation and the fractional diffusion–wave equation, J. Comput. Phys., 277 (2014) 1–15. [50] T. Wang, B. Guo, L. Zhang, New conservative difference schemes for a coupled nonlinear Schrodinger system, Appl. Math. Comput., 217 (2010) 1604-1619.
25
Journal Pre-proof
[51] J. Wang, T. Liu, H. Li, Y. Liu, S. He, Second-order approximation scheme combined with H 1 Galerkin MFE method for nonlinear time fractional convection–diffusion equation, Computers & Mathematics with Applications 73 (2017) 1182–1196. [52] Y. Yan, K. Pal, N.. J. Ford, Higher order numerical methods for solving fractional differential equations, BIT Numerical Mathematics, 54 (2014) 555–584.
of
[53] Z. Yang, Z. Yuan, Y. Nie, J. Wang, X. Zhu, F. Liu, Finite element method for nonlinear Riesz space fractional diffusion equations on irregular domains, J. Comput. Phys., 330 (2017) 863–883.
p ro
[54] Y. Yu, W. Deng, Y. Wu, J. Wu, Third order difference schemes (without using points outside of the domain) for one sided space tempered fractional partial differential equations, Appl. Numer. Math., 112 (2017) 126–145.
Pr e-
[55] Q. Yu, F. Liu, I. Turner, K. Burrage, A computationally effective alternating direction method for the space and time fractional Bloch-Torrey equation in 3–D, Applied Mathematics and Computation, 219 (8), (2012) 4082–4095. [56] Q. Yu, F. Liu, I. Turner, K. Burrage, Numerical investigation of three types of space and time fractional Bloch–Torrey equations in 2D, Central European Journal of Physics, 11 (6), (2013) 646–665. [57] Q. Yu, F. Liu, I. Turner, K. Burrage, Numerical simulation of the fractional Bloch equations, Journal of Computational and Applied Mathematics, 255(1), (2014) 635–651.
al
[58] Q. Yu, F. Liu, I. Turner, K. Burrage, Stability and convergence of an implicit numerical method for the space and time fractional Bloch-Torrey equation, Phil. Trans. R. Soc. A 371 (1990), 20120150.
urn
[59] Z. B. Yuan, Y. F. Nie, F. Liu, I. Turner, G. Y. Zhang, Y. T. Gu, An advanced numerical modeling for Riesz space fractional advection-dispersion equations by a meshfree approach, Appl. Math. Model., 40 (2016) 7816-7829. [60] M. Zayernouri, G. E. Karniadakis, Exponentially accurate spectral and spectral element methods for fractional ODEs, J. Comput. Phys., 257 (2014) 460–480.
Jo
[61] M. Zayernouri, G. E. Karniadakis, Discontinuous spectral element methods for time-and spacefractional advection equations, SIAM J. Sci. Comput., 36 (2014) B684–B707. [62] X. Zhao, Z. Z. Sun, Z. P. Hao, A fourth-order compact ADI scheme for two-dimensional nonlinear space fractional Schrodinger equation, SIAM J. Sci. Comput., 36 (2014) A2865A2886. [63] Y. Zhou, Li Peng, On the time–fractional Navier-Stokes equations, Computers & Mathematics with Applications, 73 (2017) 874–891. 26
Journal Pre-proof
[64] Y. Zhou, L. Shangerganesh, J. Manimaran, A. Debbouche, A class of time–fractional reaction– diffusion equation with nonlocal boundary condition, Math. Meth. Appl. Sci. 41 (2018) 2987– 2999.
Jo
urn
al
Pr e-
p ro
of
[65] Y. Zhou, Attractively for fractional evolution equations with almost sectorial operators, Fract. Calc. Appl. Anal. 21.3 (2018) 786–800.
27
PII: DOI: Reference:
S0377-0427(20)30030-3 https://doi.org/10.1016/j.cam.2020.112739 CAM 112739
To appear in:
Journal of Computational and Applied Mathematics
Received date : 21 July 2019 Revised date : 9 January 2020 Please cite this article as: M. Abbaszadeh, M. Dehghan and Y. Zhou, Crank–Nicolson/Galerkin spectral method for solving two-dimensional time-space distributed-order weakly singular integro-partial differential equation, Journal of Computational and Applied Mathematics (2020), doi: https://doi.org/10.1016/j.cam.2020.112739. 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.
© 2020 Elsevier B.V. All rights reserved.
Journal Pre-proof
of
Crank–Nicolson/Galerkin spectral method for solving two-dimensional time-space distributed-order weakly singular integro-partial differential equation Mostafa Abbaszadeha∗, Mehdi Dehghana , Yong Zhoub
p ro
a Department of Applied Mathematics, Faculty of Mathematics and Computer Sciences,
Amirkabir University of Technology, No. 424, Hafez Ave.,15914, Tehran, Iran b School of Mathematics and Computational Science,
Xiangtan University, Xiangtan 411105, Hunan, PR China
Pr e-
January 20, 2020
Abstract
al
The fractional PDEs based upon the distributed-order fractional derivative have several applications in physics. The two-dimensional time-space distributed-order weakly singular integro-partial differential model is investigated by a combination of finite difference and Galekrin spectral methods. A second-order finite difference formula is employed to approximate the temporal variable. In this stage, the stability and convergence of the semi-discrete scheme are proved. Then, the Galerkin spectral method based on the modified Jacobi polynomials is applied to discrete the space variable. Also, in this step, the error estimate of the full-discrete scheme is studied. Finally, two test problems have been presented to confirm the theoretical results.
urn
Keywords: Time-space fractional derivatives, convergence analysis and error estimate, finite element method, stability analysis, finite difference method.
Jo
AMS subject Classification: 65M70, 34A34.
1
Introduction
The mathematical models based upon the fractional calculations are developing in most branches of physics or engineering . ∗
Corresponding author. E-mail addresses: [email protected] , [email protected] (M. Dehghan), [email protected] (M.Abbaszadeh), [email protected] (Y. Zhou).
1
Journal Pre-proof
Pr e-
p ro
of
The fractional PDEs have been investigated by several numerical procedures and their modifications such as finite difference methods [4, 5, 12, 20, 21, 22, 37, 49, 54], spectral methods and their improvements [7, 38, 42, 60, 61], the point interpolation technique [59], FEM on complex regions [25, 53], the ADI scheme with high accuracy [62], finite element method [8], etc. Also, the interested readers can refer to [19, 29, 41, 52, 55, 56, 57, 58]. The main aim of [63] is to study the existence and uniqueness of local and global mild solutions of the Navier–Stokes equations with time-fractional derivative of order 0 < α < 1. Authors of [64] studied the time-fractional reaction-diffusion equation with nonlocal boundary condition. The main aim of [65] is to establish sufficient conditions for the existence of globally attractive solutions for the Cauchy problems in cases that semigroup is compact as well as noncompact. An error estimate has been proposed in [1] to find a second-order finite difference scheme for solving the Riesz space distributed-order diffusion equation. Authors of [2] developed a meshless upwind local radial basis function-finite difference technique for solving the time- fractional distributed-order advection–diffusion equation. The main aim of [3] is to solve the multi-dimensional generalized modified anomalous sub-diffusion equation by using alternating direction implicit-spectral element method (ADI-SEM). An error estimate has been proposed in [15] to solve the 2D weakly singular integro-partial differential equation with space and time fractional derivatives based on the finite element/finite difference scheme. Authors of [18] employed the homotopy analysis method to solve some nonlinear fractional models such as fractional KdV, K(2,2), Burgers, BBM-Burgers, cubic Boussinesq, coupled KdV, and Boussinesq-like B(m,n) equations with initial conditions. Authors of [26] found a special point for the interpolation approximation of the linear combination of multi-term fractional derivatives as based on this issue they derived a numerical differentiation formula with second-order accuracy for solving time multi-term and distributed-order fractional sub-diffusion equations. In this paper, the following two-dimensional time-space distributed-order integro-partial differential equation with the weakly singular kernel is studied [14, 35, 40] Z2
al
Z1
ω(α)0 Dαt u(x, y, t)dα + u(x, y, t) = K
1
+
Zt 0
urn
0
$(β)
− 12
(t − ζ)
∂ β u(x, y, t) ∂|x|β
+
∂ β u(x, y, t) ∂|y|β
∂ 2 u(x, y, ζ) ∂ 2 u(x, y, ζ) + dζ + f (x, y, t), ∂x2 ∂y 2
!
dβ
(x, y, t) ∈ Ω × (0, T ], (1.1)
with the initial and boundary conditions
(x, y) ∈ Ω,
(1.2)
u(x, y, t) = g(x, y, t),
(x, y) ∈ ∂Ω,
(1.3)
Jo
u(x, y, 0) = ϕ(x, y),
in which K is Newtonian contribution to the viscosity. Also, 0 Dαt u(x, y, t) is the Riemann-Liouville fractional derivative. Eq. (1.1) may be appeared in the mathematical modeling of various physical
2
Journal Pre-proof
2
p ro
of
phenomena such as the heat transfer in the materials with memory, convection and radiation systems [35] and viscoelastic mechanics [14, 16, 40]. The fractional equations with the distributedorder have been studied by many researchers for example Caputo [9] proposed the application of differential equations with distributed-order derivatives for generalizing stress-strain relations of unelastic media. Also, Caputo in [10, 11] discussed distributed-order time fractional differential equations and distributed-order space fractional differential equations, respectively and derived the solutions with closed form formulae of the classic problems. Authors of [13, 43] introduced the diffusion-like equations with distributed-order time and space fractional derivatives for the kinetic description of anomalous diffusion and relaxation phenomena [28]. Author of [27] applied distributed-order diffusion equation to discuss ultraslow and lateral diffusion processes. Also, the applications of fractional equations with space distributed-order can be found in [6, 23, 33, 44].
Arrangement of fractional derivatives
Definitions of the Riesz fractional derivatives in the x- and y-directions are [24, 36]
∂|x|β ∂ β u(x, y, t) ∂|y|β
=
−1 RL β RL β D u(x, y, t) + D u(x, y, t) , ℘ (y) x x ℘2 (y) 2 cos( βπ ) 1 2
=
−1 RL β RL β D u(x, y, t) + D u(x, y, t) , = (x) y y =2 (x) ) 1 2 cos( βπ 2
in which [24, 36]
Pr e-
∂ β u(x, y, t)
d2 1 = Γ(2 − β) dx2
Zx
d2 1 RL β D u(x, y, t) = x b Γ(2 − β) dx2
Zb
1 d2 = Γ(2 − β) dy 2
Zy
1 d2 RL β D u(x, y, t) = y d Γ(2 − β) dy 2
Zd
urn
al
RL β a Dx u(x, y, t)
RL β c Dy u(x, y, t)
a
x
c
y
u(s, y, t)
(x − s)β−1 u(s, y, t) (s − x)β−1 u(x, s, t) (y − s)β−1 u(x, s, t) (s − y)β−1
ds,
ds,
ds,
ds.
Jo
Lemma 2.1. [24, 36] Let v ∈ JS% (Ω) and let v¯ be the deployment of v outside of Ω then RL % RL a Dx v,x
RL % RL c Dy v,y
D%b v
D%d v
L2 (Ω)
L2 (Ω)
= =
RL % RL ¯,x a Dx v RL % RL ¯,y c DL v
% 2 D%b v¯ L2 (R2 ) = cos(π%) RL ¯ L2 (R2 ) , −∞ Dx v
% 2 D%d v¯ L2 (R2 ) = cos(π%) RL ¯ L2 (R2 ) . −∞ Dy v 3
Journal Pre-proof
Lemma 2.2. [48, 49] Let β ∈ (0, 1) then the following second-order weighted and shifted Gr¨ unwald difference (WSGD) approximation formula holds β 0 Dt u(x, t)
=τ
n X
−β
λβ (k)u(x, tn−k ) + O(τ 2 ),
(2.1)
λβ (k) =
and ω0β
= 1,
ωkβ
β+2 β 2 ω0 ,
k = 0,
α + 2 ωβ − β ωβ , k 2 2 k−1
Γ(k − β) = , Γ(−β)Γ(k + 1)
ωkβ
k > 0,
p ro
in which
of
k=0
β+1 β = 1− ωk−1 , k
k ≥ 1.
Lemma 2.3. [48, 49, 51] According to relation (2.1), for each N ∈ N and real vector u0 , u1 , . . . , uN ∈ RN +1 , the following inequalities hold
Pr e-
N X n X n=0 k=0
λβ (k) un−k , un ≥ 0,
n X
λβ (k) < 0.
k=0
Lemma 2.4. [47] Assuming f (t) ∈ C 1 ([0, T ]) ∩ C 3 ((0, T )) and 1
(tn − ξ)− 2 f (ξ)dξ.
al
I(f, tn ) =
Ztn 0
urn
Then
I(f, tn ) = an f (t0 ) +
an = 2
Jo
where
b0 =
1 tn2
−
t Z1 2
τ
0
1 τ
n X r=0
tZn+1
tn
1 y 2 dy
3
br f (tn−r ) + O(τ 2 ),
1 y 2 dy ,
√ 4 τ + β, 3 4
n = 1, 2, . . . ,
Journal Pre-proof
b1 =
t Z2 2
τ
1 y 2 dy
−
t0
t1
br
Zt1
√ 4 τ − β, 3
1 y 2 dy
t Zr+1 Ztr 1 1 2 = y 2 dy , y 2 dy − τ tr
r ≥ 2.
tr−1
sin(πθ1 ) + sin(πθ2 ) . sin(πα)
We set the operators
p ro
C(α, θ1 , θ2 ) =
of
Definition 2.1. [34] Let 1 < α < 1, 0 < θ1 , θ2 < α, θ1 + θ2 = α, 0 ≤ p ≤ 1 and
1 C(α, θ1 , θ2 ) (I0α+ + I1α− ) , 2
θ1 ,θ2 ,α D−1,1 : =
dk θ1 ,θ2 ,k−α I , dxk −1,1
Pr e-
θ1 ,θ2 ,α I−1,1 : =
Jm−θ,−s (x) = (1 − x)θ (1 + x)s Pmθ,s ,
θ, s > −1.
Lemma 2.5. [34] Let 1 < p < 1. If 0 < θ1 , θ2 < p such that θ1 + θ2 = α and p sin(πθ1 ) = (1 − p) sin(πθ2 ),
so for x ∈ [−1, 1], n = 0, 1, 2, . . . and k = 1, 2, . . . , n + 2 we have
al
θ1 ,θ2 ,2−α −θ1 ,−θ2 θ2 −2,θ1 −2 b α)Pn+2 I−1,1 Jn (x) = C(n, ,
urn
θ1 ,θ2 ,k−2+α −θ1 ,−θ2 e k, α)P θ2 −2+k,θ1 −2+k , D−1,1 Jn (x) = C(n, n+2−k
Also, for l = −1, 0, 1, . . . , min {m, n} we have =
θ2 ,θ1 Kn,l δm,n
=
b k, α) = Γ(n + k + α − 1) . C(n, 2k−2 n!
Z1
Dpα+κ Jm−θ1 ,−θ2 (x)Dpα+κ Jn−θ1 ,−θ2 (x)ω θ1 +l,θ2 +l (x)dx,
Z1
α+κ −θ2 ,−θ1 α+κ −θ2 ,−θ1 D1−p Jm (x)D1−p Jn (x)ω θ1 +l,θ2 +l (x)dx,
−1
Jo
θ1 ,θ2 Kn,l δm,n
b α) = 4Γ(n + α − 1) , C(n, n!
−1
θ2 ,θ1 Kn,l = 2α+1
Γ(n + θ1 + 1)Γ(n + θ2 + 1)Γ(n + l + α + 1) . (2n + α + 1)(n − l)!(n!)2 5
Journal Pre-proof
3
Analytical investigation
In the next, the order of convergence and the stability issue of the new developed methodology have been investigated. For 0 < γ < 1 we define the following functional spaces n o 1+2γ H 2γ (Ω) = u ∈ L2 (Ω) s.t. (1 + ω 2 ) 2 u b ∈ L2 (Ω) , H02γ (Ω) = u ∈ H 2γ (Ω) s.t. u = 0
on ∂Ω .
Study the semi-discrete formulation
p ro
3.1
of
in which u b is the Fourier transformation of u and
At first, we have $(β)
∂ β u(x, y, t) β
∂|x|
1
1
dβ =
MG X
(1)
wj $(βj )
j=0
∂ βj u(x, y, t) βj
∂|x|
Pr e-
Z2
+ O((MG1 )−r ),
Z2
$(β)
Z1
MG X ∂ αj u(x, y, t) ∂ α u(x, y, t) (3) ω(α ) ω(α) dα = w + O((MG2 )−r ), j j αj ∂tα ∂t j=0
∂ β u(x, y, t)
1
0
∂|y|β
1
dβ =
MG X
(1)
wj $(βj )
∂ βj u(x, y, t)
j=0
∂|y|βj
+ O((MG1 )−r ),
2
2
k=0 j=0
(3)
λαj (k)wj ω(αj )u(x, y, tn−k ) + u(x, y, tn )
urn
MG n X X
Jo
τ
−α
al
in which αj and βj are the Legendre-Gauss points on the intervals (0, 1) and (1, 2), respectively. By employing Lemmas 2.2 and 2.4, the time-discrete of Eq. (1.1) is as follows [17]
1
= K
MG X j=0
(1)
wj $(βj )
(3.1) (
∂ βj u(x, y, tn ) ∂|x|βj
+
∂ βj u(x, y, tn ) ∂|y|βj
∂ 2 u(x, y, 0) ∂ 2 u(x, y, 0) + an + ∂x2 ∂y 2 2 n X ∂ u(x, y, tn−r ) ∂ 2 u(x, y, tn−r ) + br + ∂x2 ∂y 2 r=0 + f (x, y, tn ) + O((MG1 )−r ) + O((MG2 )−r ).
6
)
Journal Pre-proof
By omitting the small term, we can get 2
1
(3) λαj (k)wj ω(αj )U n−k
+ Un = K
MG X
(1) wj $(βj )
(
∂ βj U n
∂ βj U n
)
(3.2) ∂|y|βj 2 X 2 n−r n ∂ ϕ ∂ 2ϕ ∂ U ∂ 2 U n−r + an br + f n. + 2 + + 2 2 2 ∂x ∂y ∂x ∂y r=0
k=0 j=0
j=0
+
∂|x|βj
β
of
τ −α
MG n X X
Also, the weak form of Eq. (3.2) is: Find a U n ∈ H02 (Ω) such that τ
MG n X X k=0 j=0
(3) λαj (k)wj ω(αj ) U n−k , v
1
n
+ hU , vi = K
MG X j=0
n
Bj (U , v) = Kj + Kj
βj
RL 2 ℘1 (y) Dx
βj
RL 2 =1 (x) Dy
and
n X r=0
D℘2 (y) v +
U n ,RL y
βj 2
D=2 (x) v + K
β
2 cos( 2j π)
al
Kj =
βj 2
U n ,RL x
(3.3)
br ∇U n−r , ∇v + hf n , vi , ∀v ∈ H0 (Ω),
Pr e-
+ an h∆ϕ, vi − where
(1)
wj $(βj )Bj (U n , v)
p ro
2
−α
βj
RL n RL 2 x D℘2 (y) U ,℘1 (y)
βj
β 2
βj 2
Dx v
RL n RL 2 y D=2 (x) U ,=1 (x)
βj 2
Dy v
.
Lemma 3.1. [30] For an in Eq. (2.2), we have m X n=0
m X
√ |an | ≤ 4τ T ,
urn
τ
n=0
1
an ≤ Cτ 2 .
Jo
Proof. According to definition an , we have tZn+1 tZn+1 m m m X X X 1 1 1 1 1 tn2 − tn2 + 1 τ |an | = 2τ y 2 dy ≤ 2τ y 2 dy τ τ n=0 n=0 n=0 ≈ 2τ
m h X n=0
1 2
1 2
i
tn + y1 ≤ 2τ
tn
m X n=0
−1 tn 2
tn
≤ 2τ
Zmτ 0
7
√ dx √ ≤ 4τ T , x
,
Journal Pre-proof
in which y1 ∈ (tn , tn+1 ). In the current paper, we employ the following inequalities to prove the theorems.
ξ0 ≤ ε0 ,
n−1
of
Proposition 3.1. (Discrete Gronwall inequality) [39] Assume that $n is a non-negative sequence and that the sequence ξn satisfies
n−1
X X νk + $k ξk , ξn ≤ ε0 + k=0
Then, if ε0 ≥ 0 and ν0 ≥ 0, it follows ε0 +
k=0
νk
!
exp
n−1 X
$k
k=0
!
n ≥ 1.
,
Pr e-
ξn ≤
n−1 X
p ro
k=0
n ≥ 1.
Lemma 3.2. (Gronwall’s inequality) [50] Suppose that the discrete function {wn | n = 0, 1, 2, . . . , N ; N τ = T } , satisfies the inequality
wn ≤ γ + τ
θl wl ,
l=1
al
where γ and θl are nonnegative constants. Then
n X
kwn kL2 (Ω) ≤ γ exp 2τ
urn
β
n X l=1
θl
!
.
Theorem 3.2. Let U n ∈ H02 (Ω). The formulation (3.1) is unconditionally stable. e n be the approximate solution of (3.3), then the roundoff error equation is as follows Proof. Let U 2
MG n X X k=0 j=0
(3) λαj (k)wj ω(αj ) E n−k , v
Jo
τ
−α
1
n
+ hE , vi =
MG X
(1)
wj $(βj )Bj (E n , v)
j=0
n β
0 X
+ an ∆E , v − br ∇E n−r , ∇v , ∀v ∈ H02 (Ω), r=0
8
Journal Pre-proof
e n and U e n is approximation solution of Eq. (3.4). If we set v = E n , then in which E n = U n − U τ −α
k=0 j=0
(3) λαj (k)wj ω(αj ) E n−k , E n + hE n , E n i 1
=
MG X j=0
(1) wj $(βj )Bj (E n , E n )
n
0 n X br ∇E n−r , ∇E n . + an ∆E , E − r=0
of
2
MG n X X
p ro
Also from [8] we obtain βj βj βj βj RL RL n n n n n RL n RL 2 2 2 2 Bj (E , E ) = |2Kj | ℘1 (y) Dx E ,x D℘2 (y) E + |2Kj | =1 (x) Dy E ,y D=2 (x) E
βj βj βj βj RL RL n RL n n RL n 2 2 2 2 ≥ min{|2Kj | , |2Kj |} ℘1 (y) Dx E ,x D℘2 (y) E + =1 (x) Dy E ,y D=2 (x) E
Pr e-
≥ Kj kE n k2H β (Ω) .
(3.4)
Summing Eq. (3.4) for n from 0 to N , gives 2
τ −α
MG N X n X X n=0 k=0 j=0
n=0
= +
N X
Using Lemma 2.1, we have
n=0 k=0 j=0
n=0 j=0
an ∆E 0 , E n −
N X n X n=0 r=0
(3.5)
br ∇E n−r , ∇E n .
M2
N X n G X
X
(3) (3) λαj (k)wj ω(αj ) E n−k , E n = wj ω(αj ) λαj (k) E n−k , E n ≥ 0.
urn
2
MG N X n X X
al
n=0
τ −α
1
MG N X N X X
n−k n (1) (3) n n hE , E i + wj $(βj )Bj (E n , E n ) λαj (k)wj ω(αj ) E , E +
n=0 k=0
j=0
Also, form Eq. (3.4), we can write
n=0 j=0
and
1
(1) wj $(βj )Bj (E n , E n )
Jo
1
MG N X X
=
MG X
(1) wj $(βj )
j=0
−
N X n X n=0 r=0
N X n=0
1
n
n
Bj (E , E ) ≥
MG X
br ∇E n−r , ∇E n ≤ 0.
9
j=0
(1) wj $(βj )Kj
N X n=0
kE n k2H β (Ω) , (3.6) (3.7)
Journal Pre-proof
Thus, using the above inequalities, Eq. (3.5) can be rewritten as 1
MG X
(1) wj $(βj )Kj
N X n=0
j=0
kE n k2H β (Ω)
+
N X n=0
n
N X
n
hE , E i =
n=0
an ∆E 0 , E n .
(3.8)
By changing the index, we will have
j=0
1
≤
MG X
i=0
i=0
It can be concluded that kE n k2L2 (Ω)
n n n X X
i 2
i i X
0 i
E β + E , E = a ∆E , E . i H (Ω)
(1)
wj $(βj )Kj
(3.9)
i=0
n n n X X
j 2
j j X
E β = aj ∆E 0 , E j + E , E H (Ω) j=0
j=0
j=0
of
(1)
wj $(βj )Kj
p ro
1
MG X
j=0
0
+
n X j=1
Pr e-
≤ a0 ∆E , E
0
aj ∆E 0 , E j
n X
0 2
0
≤ |a0 | ∇E L2 (Ω) + ∆E L2 (Ω) aj E j L2 (Ω)
Thus, we get the following inequality
j=1
n X
0 2
0
≤ 1 + |a0 | ∇E L2 (Ω) + ∆E L2 (Ω) aj E j L2 (Ω) .
Jo
urn
kE kL2 (Ω)
n X
2
≤ 1 + |a0 | ∇E 0 L2 (Ω) + ∆E 0 L2 (Ω) aj E j L2 (Ω) .
al
n
j=1
10
j=1
(3.10)
Journal Pre-proof
j=1
aj ≤
n X
(aj )2 , we can write Eq. (3.10) as follows
j=1
n
kE kL2 (Ω)
n X
0 2
0
1
≤ 1 + |a0 | ∇E L2 (Ω) + ∆E L2 (Ω) τ (aj )2 E j L2 (Ω) τ j=1 ! n X
0 2
0 1 (aj )2 ≤ 1 + |a0 | ∇E L2 (Ω) exp τ ∆E L2 (Ω) τ j=1 !2 n X
0 2 ≤ 1 + |a0 | ∇E L2 (Ω) exp C aj
≤
and this completes the proof.
p ro
j=1
of
Since
n X
2 1 + |a0 | ∇E 0 L2 (Ω) exp (Cτ ) , β
Pr e-
Theorem 3.3. Let un , U n ∈ H02 (Ω) then relation (3.1) is convergent.
Proof. The exact and approximate solutions can be found from the following relations 2
τ −α
MG n X X
1
(3) λαj (k)wj ω(αj )un−k
k=0 j=0
+ un = K
MG X
(1) wj $(βj )
(
∂ βj un
+
∂ βj un
)
(3.11) ∂|y|βj X 2 n−r 2 n ∂ 2 un−r ∂ u ∂ ϕ ∂ 2ϕ + 2 + br + + an + f n + Rτ , 2 2 2 ∂x ∂y ∂x ∂y r=0 j=0
∂|x|βj
3
al
in which there exists a positive constant C such that |Rτ | ≤ Cτ 2 and 1
2
MG n X X
(3) λαj (k)wj ω(αj )U n−k
+U
n
= K
k=0 j=0
MG X
(1) wj $(βj )
(
∂ βj U n
j=0
∂|x|βj
+
∂ βj U n
)
(3.12) ∂|y|βj 2 X 2 n−r n ∂ ϕ ∂ 2ϕ ∂ U ∂ 2 U n−r + an + f n. + 2 + br + 2 2 ∂x2 ∂y ∂x ∂y r=0
urn
τ
−α
Jo
We set Qn = un − U n then subtracting (3.12) from (3.11) gives 2
τ −α
MG n X X k=0 j=0
(3) λαj (k)wj ω(αj )Qn−k
1
+ Qn = K
+
MG X j=0
n X r=0
11
(1) wj $(βj )
br
2
(
n−r
∂ Q ∂x2
∂ βj Qn ∂|x|βj 2
+
+
n−r
∂ Q ∂y 2
∂ βj Qn
∂|y|βj
) (3.13)
+ Rτ .
Journal Pre-proof
Thus, the weak form of Eq. (3.13) is τ −α
k=0 j=0
(3) λαj (k)wj ω(αj ) Qn−k , ξ =−
n X r=0
1
+ hQn , ξi + K
MG X j=0
br ∇Qn−r , ∇ξ + hRτ , ξi ,
β
ξ ∈ H02 (Ω).
If we set ξ = Qn , then τ −α
k=0 j=0
(3) λαj (k)wj ω(αj ) Qn−k , Qn =−
n X r=0
1
n
n
+ hQ , Q i + K
n=0 k=0 j=0
=−
N X n X n=0 r=0
1
+
1
(1)
N X n=0
n
1
n
hQ , Q i + K
MG N X X n=0 j=0
(1)
wj $(βj )Bj (Qn , Qn )
n=0
wj $(βj )Bj (E n , Qn ) =
MG X
(1)
wj $(βj )
j=0
urn
n=0 j=0
(3.14)
N
X n−r n br ∇Q , ∇Q + hRτ , Qn i.
Now we have MG N X X
j=0
Pr e-
τ −α
(3) λαj (k)wj ω(αj ) Qn−k , Qn
al
2
(1)
wj $(βj )Bj (Qn , Qn )
br ∇Qn−r , ∇Qn + hRτ , Qn i .
Summing Eq. (3.14) for n = 0 : N , yields MG N X n X X
MG X
p ro
2
MG n X X
(1)
wj $(βj )Bj (Qn , ξ)
of
2
MG n X X
N X n=0
(3.15) 1
Bj (Qn , Qn ) ≥
MG X
(1)
wj $(βj )Kj
j=0
N X n=0
kQn k2H β (Ω) .
Also, using Lemma 2.3, concludes 2
MG N X n X X n=0 k=0 j=0
(3) λαj (k)wj ω(αj ) Qn−k , E n
Jo
τ
−α
−
N X n X n=0 r=0
2
=
MG X j=0
(3) wj ω(αj )
N X n X n=0 k=0
br ∇Qn−r , ∇Qn ≤ 0,
12
λα (k) Qn−k , Qn ≥ 0,
Journal Pre-proof
then, Eq. (3.15) can be rewritten as follows 1
≤ τ
i=0
i=0
i=0
i=0
n X
i
i 3
Q 2 2 kRτ kL2 (Ω) Q L2 (Ω) = Cτ L (Ω)
"
i=0
n X
1 2
of
# n n X X
i
i 2 1
Q 2 = Cτ τ Q L2 (Ω) ≤ Cτ τ+ L (Ω) 2 i=0 i=0 i=0 " # n n X X
i 2
i 2 1
Q 2 , Cτ T + ≤ Q L2 (Ω) ≤ Cτ T + Cτ L (Ω) 2 i=0 i=0
kQn k2L2 (Ω)
n X
i 2
Q 2 ≤ Cτ T exp (Cτ (n + 1)) ≤ C∗ τ, ≤ Cτ T + Cτ L (Ω) i=0
in which this totalizes the proof.
3.2
n n n X X
i 2
i i X
Q β + Rτ , Qi Q ,Q = H (Ω)
p ro
thus, we obtain
j=0
n X
≤
(1) wj $(βj )Kj
Pr e-
kQn k2L2 (Ω)
MG X
Analysis of full-discrete scheme
Let weight function ω(x) be positive, thus the weighted Sobolev space can be defined as Z 2 2 Lω (Ω) = u : u (x)ωdx < +∞ ,
with corresponding norm
al
Ω
urn
kukL2ω (Ω) =
Z Ω
21
u2 (x)ωdx .
For known values p, 1 < α < 2, 0 < µ, ν < α, l = −1, 0, 1, . . . , m and m ∈ N, we define m Bα,p,µ,ν = u ∈ L2ω−µ,−ν (Ω) :
Jo
Let Phm be the interpolation operator and
Dpα+1 u ∈ L2ων+l,µ+l (Ω) f or − 1 ≤ l ≤ m .
F−µ,−ν (Ω) = {φ = (1 − x)µ (1 + x)ν ϕ : ϕ ∈ PN } , N
be the finite dimensional fractional polynomial space. m Lemma 3.3. [34] Let 1 < α < 2 and u ∈ Bα,p,µ,ν . For for values p, 1 < α < 2, 0 < µ, ν < α and
13
Journal Pre-proof
l = −1, 0, 1, . . . , m, we have
α+l −µ,−ν
Dp PN u − u L2
ω ν+l,µ+l
−µ,−ν
P u − u L2 N
ω −µ,−ν
−µ,−ν
P u − u L2 N
ω −µ,−ν
(Ω)
(Ω)
≤ cN
l−m 2
s
(N − m + 1)!
Dpα+l u 2 , L ν+m,µ+m (Ω) (N − l + 1)! ω
α+l
Dp u 2 L
ω ν+m,µ+m
(Ω)
s
,
(for fixed value m),
(N − m + 1)!
Dpα+m u 2 , 0 ≤ m ≤ N, L ν+m,µ+m (Ω) (N + m + 1)! ω s
(N − m + 1)!
Dpα+m u 2 ≤ cN −α−m , (fixed value m). L ν+m,µ+m (Ω) (N + m + 1)! ω ≤ cN −α
(Ω)
(Ω)
of
ω ν+l,µ+l
≤ cN
l−m 2
p ro
α+l −µ,−ν
Dp PN u − u L2
We define the projection operator
Πh : H0β (Ω) ∩ C 0 (Ω) → Xh ,
Pr e-
such that
B(Πh z, w) = B(z, w),
∀ w ∈ Xh .
β
Theorem 3.4. Assume un = u(x, y, tn ) ∈ H02 (Ω) and Uhn (x, y, tn ) ∈ Xh are the exact and approximate solutions of (1.1), respectively. Thus 3 kun − Uhn k ≤ C ∗ τ 2 + N −β−m . Proof. We can write τ
MG n X X
1
(3) λαj (k)wj ω(αj )un−k
2
τ −α
MG n X X k=0 j=0
(3) λαj (k)wj ω(αj )UNn−k
Jo
and
+U
n
= K
(1) wj $(βj )
(
∂ βj un
j=0
+ UNn = K
1
MG X
+ an
(1) wj $(βj )
j=0
∂ 2ϕ ∂ 2ϕ + 2 ∂x2 ∂y
14
∂|x|βj
(
+
∂ βj un
)
(3.16) ∂|y|βj 2 X 2 n−r n ∂ ϕ ∂ 2ϕ ∂ 2 un−r ∂ u + an + 2 + br + + f n + Rτ , 2 2 ∂x2 ∂y ∂x ∂y r=0
urn
k=0 j=0
MG X
al
2
−α
∂ βj UNn ∂|x|βj
+
n X r=0
+
br
∂ βj UNn
∂|y|βj
)
∂ 2 UNn−r ∂ 2 UNn−r + ∂x2 ∂y 2
(3.17)
+ f n.
Journal Pre-proof
Subtracting (3.17) from (3.16) yields τ −α
k=0 j=0
(3)
λαj (k)wj ω(αj ) un−k − UN n−k
+ (un − UN n ) 1
= K
MG X
(1) wj $(βj )
j=0
n X
+
br
Thus, the weak form of Eq. (3.18) is 2
τ
MG n X X k=0 j=0
j=0
∂|x|βj n−r
+
∂ βj (un − UN n )
∂ 2 un−r − UN ∂x2
(1)
wj $(βj )Bj (un − UN n , ϑ) = −
n X r=0
br ∇ un−r − UN n−r , ∇ϑ + hRτ , ϑi .
The use of this equation and the notations
ZNn = PN un − UNn , yield τ −α
k=0 j=0 2
(3) λαj (k)wj ω(αj ) XNn−k , ϑ + hXNn , ϑi
urn
MG n X X
(3) λαj (k)wj ω(αj ) ZNn−k , ϑ + hZNn , ϑi
k=0 j=0
1
+K
MG X
(1) wj $(βj )Bj
(ZNn , ϑ)
j=0
Jo
+τ
−α
XNn = un − PN un ,
al
2
MG n X X
15
= −
n X r=0
∂|y|βj
∂ 2 un−r − UN n−r + ∂y 2
Pr e-
+K
∂ βj (un − UN n )
(3) λαj (k)wj ω(αj ) un−k − UN n−k , ϑ + hun − UN n , ϑi
1
MG X
(
p ro
r=0
−α
(3.18)
of
2
MG n X X
br ∇ZNn−r , ∇ϑ + hRτ , ϑi .
!
)
+ Rτ .
Journal Pre-proof
We set ϑ = ZNn , then τ −α
k=0 j=0 2
+τ −α
MG n X X k=0 j=0
(3) λαj (k)wj ω(αj ) ZNn−k , ZNn + hZNn , ZNn i
(3) λαj (k)wj ω(αj ) XNn−k , ZNn + hXNn , ZNn i 1
(1) wj $(βj )Bj
j=0
(ZNn , ZNn )
= −
Summing Eq. (3.19) for n = 0 : N , confirms that τ
MG K X n X X n=0 k=0 j=0 2
+τ
−α
MG N X n X X n=0 k=0 j=0
(3) λαj (k)wj ω(αj ) ZNn−k , ZNn
MG N X X
n=0 k=0 j=0
kZNn k2L2 (Ω)
(3.20)
(1)
N X n X n=0 r=0
N
X hRτ , ZNn i. br ∇ZNn−r , ∇ZNn + n=0
(3) λαj (k)wj ω(αj ) ZNn−k , ZNn ≥ 0,
(3.21)
al
τ −α
n=0
wj $(βj )Bj (ZNn , ZNn ) = −
Now, since 2
+
n=0
n=0 j=0
MG K X n X X
N X
N X
(3) λαj (k)wj ω(αj ) XNn−k , ZNn + hXNn , ZNn i 1
+K
r=0
br ∇ZNn−r , ∇ZNn + hRτ , ZNn i . (3.19)
Pr e-
2
−α
n X
p ro
+K
MG X
of
2
MG n X X
−
urn
ZN0 = 0,
N X n X n=0 r=0 1
MG N X X
(1) wj $(βj )Bj (ZNn , ZNn )
Jo
n=0 j=0
br ∇ZNn−r , ∇ZNn
(3.22)
≤ 0,
(3.23) 1
=
MG X j=0
(1) wj $(βj )
N X n=0
1
≥
MG X j=0
(1) wj $(βj )Kj
Bj (ZNn , ZNn )
N X n=0
kZNn k2H β (Ω) ,
! 2 2 MG N n MG N n X X 1 −α X X X 1 −α X (3) (3) n 2 n 2 τ λαj (k)wj ω(αj ) kZN kL2 (Ω) = τ wj ω(αj ) kZN kL2 (Ω) λαj (k) < 0, 2 2 n=0 k=0 j=0 n=0 j=0 k=0 16
Journal Pre-proof
Eq. (3.20) can be rewritten as follows N X n=0
1
kZNn k2L2 (Ω)
+
MG X
N X
(1) wj $(βj )Kj
n=0
j=0
kZNn k2H β (Ω)
2
of
N n MG
2 1 −α X X X (3) ≤ τ λαj (k)wj ω(αj ) XNn−k L2 (Ω) 2 n=0 k=0 j=0
N N N N X X 1X n 2 1X n 2 2 n 2 + kRτ kL2 (Ω) + kX k 2 + kZ k 2 + kZ k 2 , 4 n=0 N L (Ω) n=0 N L (Ω) n=0 4 n=0 N L (Ω)
p ro
or 1
MG n n X X X
i 2
i 2 (1)
ZN 2
ZN β + 2 w $(β )K j j j L (Ω) H (Ω) i=0
j=0
i=0
2
MG n X i X X
+ 2
2 (3) λαj (k)wj ω(αj ) XNi−k L2 (Ω)
Pr e-
≤ τ −α
i=0 k=0 j=0
n n X X
k 2
XN 2 + 2 kRτ k2L2 (Ω) . L (Ω) k=0
k=0
Changing the indices of summations arrives at
1
MG n n X X X
i 2
i 2 (1)
ZN β
ZN L2 (Ω) + 2 ≤ wj $(βj )Kj H (Ω) i=0
i=0
j=0
al
kZNn k2L2 (Ω)
2
MG n X i X X i=0 k=0 j=0
urn
≤ τ
−α
+ 2
2 (3) λαj (k)wj ω(αj ) XNi−k L2 (Ω)
n n X X
k 2
XN 2 + 2 kRτ k2L2 (Ω) L (Ω) k=0
k=0
M2
Jo
≤ τ −α C2 N + 2C2
n X k=0
n X i X G X −β−m 2
N
(3)
λαj (k)wj ω(αj )
i=0 k=0 j=0
−β−m 2
17
+ 2C1
n X k=0
9
τ 4.
Journal Pre-proof
Finally, we obtain N
k=0
−β−m 2
≤ 2τ (C1 + C2 ) ≤ 2 (C1 + C2 ) τ
k=0
3
τ 2 + N −β−m
k=0
n X
9
τ4
3
2
τ 2 + N −β−m
2
k=0
3 ≤ C exp (T ) τ + N −β−m ≤ C τ 2 + N −β−m ,
which completes the proof.
4
n X
+ 2τ C1
n X
of
≤ 2τ C2
n X
3 2
p ro
τ
kZNn k2L2 (Ω)
Numerical investigation
Pr e-
To verify the theoretical results, we apply the new numerical procedure for two examples. The used computational order is E1 E
log10 2
Rate =
h1 h2
log10
,
such that E1 and E2 are the absolute errors related to h1 and h2 mesh sizes, respectively.
4.1
Example 1.
ω(α)0 Dαt u(x, y, t)dα
0
+ u(x, y, t) =
Z2
$(β)
∂ β u(x, y, t)
urn
Z1
al
For the first model, we consider the following example
1
+
Zt
(t − ζ)
+
∂|y|β
∂ 2 u(x, y, ζ) ∂ 2 u(x, y, ζ) dζ + f (x, y, t), + ∂x2 ∂y 2
Jo
0
− 21
∂|x|β
∂ β u(x, y, t)
18
!
dβ
(x, y, t) ∈ Ω × (0, T ],
Journal Pre-proof
in which ω = Γ(2 − α) and $(β) = −2Γ(5 − β) cos f (x, y, t) = 10x2 (1 − x2 )y 2 (1 − y 2 )
(t − 1) ln(t)
πβ 2
and also
− 10ty 2 (1 − y 2 ) [φ1 (x) + φ1 (1 − x) + φ2 (x) + φ2 (1 − x)] 3 40 (1 − x2 )y 2 (1 − y 2 )t 2 , 3
of
− 10tx2 (1 − x2 ) [φ1 (y) + φ1 (1 − y) + φ2 (y) + φ2 (1 − y)] −
Pr e-
p ro
2 ϑ3 − ϑ2 3ϑ − 2ϑ ϑ2 − ϑ φ1 (ϑ) = Γ (5) − 2Γ (4) − ln ϑ ln ϑ (ln ϑ)2 Γ (3) 2 5ϑ 3 2ϑ + − , 6ϑ − 2 − + + ln ϑ ln ϑ ln ϑ (ln ϑ)2 (ln ϑ)2 3 ϑ4 − ϑ3 4ϑ − 3ϑ2 ϑ3 − ϑ2 φ2 (ϑ) = Γ (5) − 2Γ (4) − ln ϑ ln ϑ (ln ϑ)2 Γ (3) 2ϑ 1 2ϑ2 2 2 + + 12ϑ − 6ϑ − 7ϑ − 5ϑ − , ln ϑ ln ϑ ln ϑ ln ϑ with the exact solution
Jo
urn
al
u(x, y, t) = 10tx2 (1 − x2 )y 2 (1 − y 2 ).
Figure 1: Error obtained for Example 1.
19
Journal Pre-proof
Table 1 Errors and computational orders at T = 0.5 for Example 1 N = 10 Rate
3.7826 × 10
1/10
7.6644 × 10
1/20
1.8202 × 10
1/40
−4
−
−5
2.0022
6.9621 × 10−8
2.0009
2.7865 × 10
1/320
Rate
CPU time
−
5
−6
2.0664
25
−7
2.0234
46
−8
2.0092
77
−8
2.0039
93
8.2710 × 10
2.0621 × 10
2.0064
−7
−6
3.3295 × 10
2.0208
1.1163 × 10
1/640
1.3536 × 10
2.0741
−6
−6
1/160
2.6694 × 10
2.3031
−5
4.4853 × 10
1/80
L∞
of
L∞
−9
2.0018
126
1.2865 × 10−9
2.0008
178
5.1489 × 10
p ro
τ
N = 20
4.2
Pr e-
Table 1 is based on the N = 10, N = 20 and 16 Gauss–Legendre-Lobatto points to integration. In Table 1, we can observe that the theoretical results are complied to the computational conclusions. Figure 1 illustrates the error obtained as a function based upon the different values of N and τ = 10−4 for left panel and different values of time step and N = 20 in the right panel for Example 1.
Example 2.
For the second example, we consider a test problem without exact solution Z1
ω(α)0 Dαt u(x, y, t)dα + u(x, y, t) =
0
Z2
$(β)
1
− 21
(t − ζ)
with initial data
∂ 2 u(x, y, ζ) ∂ 2 u(x, y, ζ) + dζ, ∂x2 ∂y 2
urn
0
∂|x|β
al
+
Zt
∂ β u(x, y, t)
u(x, y, 0) = exp (−10 [(x − 0.5) + (y − 0.5)]) ,
+
∂ β u(x, y, t) ∂|y|β
!
dβ
(x, y, t) ∈ Ω × (0, T ],
(x, y) ∈ Ω,
and homogeneous boundary condition
and also
Jo
u(x, y, t) = 0,
(x, y) ∈ ∂Ω,
ω(α) = Γ (2 − α) ,
0 ≤ t ≤ 1,
$(β) = exp (−β) .
To check the convergence of the present numerical technique, we adopt the following strategy: • Select a N ∗ and τ small enough and obtain the numerical solution UNτ ∗ using two parameters. 20
Journal Pre-proof
• Set N1 , N2 , . . . in which N ∗ > N i for i = 1, 2, . . . and obtain the related numerical solution UNτ 1 , UNτ 2 , . . . based on the N1 , N2 , . . .. • Put
∗
LN ∞ =
max |Uhτ∗ (xj ) − Uhτi (xj )| ,
i = 1, 2, . . . .
1≤j≤M −1
of
Note that the dimensions of UNτ ∗ and UNτ i are different but we consider the common nodes in these solutions. Table 2
Errors and computational orders at T = 0.5 for Example 2 τ = 10−5
∗
∗
p ro
τ = 10−3 N
LN ∞
Rate
LN ∞
5
1.0378 × 10−3
−
1.1754 × 10−6
1.9749 × 10
1.2057 × 10
−8
7.6627
1.0694 × 10−9
8.4209
2.6951 × 10
20 30 40
5.7156
−7
6.1953
CPU time
−
11
2.2442 × 10
−6
5.7108
71
−8
7.7133
131
3.2966 × 10
−10
8.5812
241
11.6660
388
1.0694 × 10
1.1496 × 10−11
Pr e-
10
−5
Rate
Table 3
Errors and computational orders at T = 0.5 for Example 2 N ∗ = 40, τ = 10−4 τ
L∞
1/10
1.8781 × 10−4
1/40 1/80 1/160
L∞
−
1.6384 × 10−6
1.9921
−5
1.9985
2.9543 × 10−6
1.9997
1.1814 × 10 7.3862 × 10 1.8466 × 10
−7 −7
1.9999 2.0000
urn
1/320
4.7210 × 10
Rate
al
1/20
−5
N ∗ = 80, τ = 10−4 Rate −
−7
4.4409 × 10
1.8833
−7
1.9711
2.8459 × 10−8
1.9928
−9
1.9981
−9
1.9991
1.1327 × 10 7.1242 × 10 1.7821 × 10
Jo
Table 2 presents the errors and computational orders obtained at T = 0.5 for Example 2. In Table 2, we assume that the reference solution is based on N ∗ = 80 and τ = 10−5 . In this table, two values are fixed for time step τ and then the number of basis N is increasing. To construct the Table 3, we consider two reference solutions such that one of them is based on N ∗ = 40 and τ = 10−4 and another is N ∗ = 80 and τ = 10−4 . Table 3 shows that the proposed numerical procedure has second-order accuracy in the temporal direction. It must be noted that an extrapolation approach is used to improve the numerical results. For this reason, computational order of temporal variable has been increased.
21
Journal Pre-proof
5
Conclusion
of
A high-order numerical scheme is proposed to solve the two-dimensional time-space distributedorder weakly singular integro-partial differential equation. The temporal direction is discretized by using a second-order difference formula. Furthermore, to get a numerical procedure with high accuracy, derivatives of spatial directions are approximated by a Galerkin spectral idea based upon the modified Jacobi polynomials. The stability, convergence and error estimation of the semi- and full-discrete plan have been studied in details. In the end, the theoretical results have been checked by using two test problems.
p ro
Acknowledgments
The authors are very grateful to both reviewers for carefully reading this paper and for their comments and suggestions which have improved the paper.
Pr e-
References
[1] M. Abbaszadeh, Error estimate of second-order finite difference scheme for solving the Riesz space distributed-order diffusion equation, Applied Mathematics Letters 88, (9)2019) 179-185. [2] M. Abbaszadeh, M. Dehghan, Meshless upwind local radial basis function-finite difference technique to simulate the time- fractional distributed-order advection–diffusion equation, Engineering with Computers, In press (2019).
al
[3] M. Abbaszadeh, M. Dehghan, Y. Zhou, Alternating direction implicit-spectral element method (ADI-SEM) for solving multi-dimensional generalized modified anomalous sub-diffusion equation, Computers & Mathematics with Applications, 78 (2019) 1772-1792.
urn
[4] A. A. Alikhanov, A new difference scheme for the time fractional diffusion equation, J. Comput. Phys., 280 (2015) 424-438. [5] A. A. Alikhanov, Numerical methods of solutions of boundary value problems for the multiterm variable-distributed order diffusion equation, Appl. Math. Comput., 268 (2015) 12-22.
Jo
[6] T. M. Atanackovic, S. Pilipovi, D. Zorica, Distributed-order fractional wave equation on a finite domain, stress relaxation in a rod, International Journal of Engineering Science, 49 (2011) 175-190. [7] A. H. Bhrawy, M. A. Zaky, An improved collocation method for multi-dimensional space-time variable-order fractional Schrodinger equations, Appl. Numer. Math., 111 (2017) 197-218. [8] W. Bu, Y. Tang, Y. Wu, J. Yang, Finite difference/finite element method for two-dimensional space and time fractional Bloch-Torrey equations, J. Comput. Phys. 293 (2015) 264-279. 22
Journal Pre-proof
[9] M. Caputo, Mean fractional-order-derivatives differential equations and filters, Annali dellUniversita di Ferrara, 41 (1995) 73–84. [10] M. Caputo, Distributed order differential equations modelling dielectric induction and diffusion, Fractional Calculus and Applied Analysis, 4 (2001) 421-442.
of
[11] M. Caputo, Diffusion with space memory modelled with distributed order space fractional differential equations, Annals of Geophysics, 46 (2003). [12] C. Celik, M. Duman, Crank-Nicolson method for the fractional diffusion equation with the Riesz fractional derivative, Journal of Computational Physics, 231 (2012) 1743-1750.
p ro
[13] A. Chechkin, R. Gorenflo, I. Sokolov, Retarding subdiffusion and accelerating superdiffusion governed by distributed-order fractional diffusion equations, Physical Review E, 66 (2002) 046129. [14] R. M. Christensen, Theory of Viscoelasticity, Academic Press, New York, 1971.
Pr e-
[15] M. Dehghan, M. Abbaszadeh, Error estimate of finite element/finite difference technique for solution of two-dimensional weakly singular integro-partial differential equation with space and time fractional derivatives, Journal of Computational and Applied Mathematics 356 (2019) 314-328. [16] M. Dehghan, Solution of a partial integro–differential equation arising from viscoelasticity, International Journal of Computer Mathematics 83 (1), (2006) 123–129.
al
[17] M. Dehghan, Finite difference procedures for solving a problem arising in modelling and design of certain optoelectronic devices, Mathematics and Computers in Simulation, 71 (1), (2006) 16-30.
urn
[18] M. Dehghan, J. Manafian, A. Saadatmandi, Solving nonlinear fractional partial differential equations using the homotopy analysis method, Numerical Methods for Partial Differential Equations, 26(2), (2010), 448-479. [19] W. Deng, J. S. Hesthaven, Local discontinuous Galerkin methods for fractional ordinary differential equations, BIT Numerical Mathematics 55 (4), (2015) 967-985.
Jo
[20] H. Ding, General Pade approximation method for time–space fractional diffusion equation, Journal of Computational and Applied Mathematics, 299 (2016) 221-228. [21] H. Ding, C. P. Li, Y. Q. Chen, High-order algorithms for Riesz derivative and their applications (II), Journal of Computational Physics, 293 (2015) 218-237. [22] H. Ding, C. P. Li, High-order numerical algorithms for Riesz derivatives via constructing new generating functions, Journal of Scientific Computing, 71 (2017) 759–784.
23
Journal Pre-proof
[23] C. H. Eab, S. C. Lim, Fractional Langevin equations of distributed order, Physical Review E., 83 (2011), 031136. [24] V. J. Ervin, J. P. Roop, Variational formulation for the stationary fractional advection dispersion equation, Numer. Methods Partial Differ. Equ. 22 (2006) 558-576.
of
[25] W. Fan, F. Liu, X. Jiang, I. Turner, A novel unstructured mesh finite element method for solving the time-space fractional wave equation on a two-dimensional irregular convex domain, Fractional Calculus and Applied Analysis 20.2 (2017) 352-383.
p ro
[26] G. H. Gao, A. A. Alikhanov, Z. Z. Sun, The temporal second order difference schemes based on the interpolation approximation for solving the time multi-term and distributed order fractional sub-diffusion equations, J. Sci. Comput., 73 (2017) 93-121. [27] A. N. Kochubei, Distributed order calculus and equations of ultraslow diffusion, Journal of Mathematical Analysis and Applications, 340 (2008) 252-281.
Pr e-
[28] J. Li, F. Liu, L. Feng, I. Turner, A novel finite volume method for the Riesz space distributedorder diffusion equation, Computers and Mathematics with Applications, 46 (2017) 536-553. [29] H. L. Liao, P. Lyu, S. Vong, Second-order BDF time approximation for Riesz space-fractional diffusion equations, Int. J. Comput. Math., 95 (2018) 144–158. [30] M. Luo, D. Xu, L. Li, A compact difference scheme for a partial integro-differential equation with a weakly singular kernel, Appl. Math. Model. 39 (2015) 947-954.
al
[31] J. E. Macias-Diaz, Numerical study of the process of nonlinear supratransmission in Riesz space-fractional sine-Gordon equations, Commun. NonLinear Sci. Numer. Simul. 46 (2017) 89-102.
urn
[32] J. E. Macias-Diaz, A structure-preserving method for a class of nonlinear dissipative wave equations with Riesz space-fractional derivatives, Journal of Computational Physics, 351 (2017) 40-58. [33] F. Mainardi, A. Mura, R. Gorenflo, M. Stojanovic, The two forms of fractional relaxation of distributed order, Journal of Vibration and Control, 13 (2007) 1249–1268.
Jo
[34] Z. Mao, G. Em Karniadakis, A spectral method (of exponential convergence) for singular solutions of the diffusion equation with general two–sided fractional derivative, SIAM J. Numer. Anal., 56(1) (2018) 24–49. [35] R. K. Miller, An integro-differential equation for grid heat conductors with memory, J. Math. Anal. Appl. 66 (1978) 313–332. [36] J. P. Roop, Variational solution of the fractional advection dispersion equation, PhD thesis, Clemson University, 2004. 24
Journal Pre-proof
[37] H. K. Pang, H. W. Sun, Fourth order finite difference schemes for time-space fractional subdiffusion equations, Comput. Math. Appl., 71 (2016) 1287-1302. [38] E. Pindza, K. M. Owolabi, Fourier spectral method for higher order space fractional reactiondiffusion equations, Communications in Nonlinear Science and Numerical Simulation 40 (2016) 112-128.
of
[39] A. Quarteroni, A. Valli, Numerical Approximation of Partial Differential Equations, SpringerVerlag, New York, 1997.
p ro
[40] M. Rcnardy, Mathematical analysis of viscoelastic flows, Ann. Rev. Fluid Mech., 21 (1989) 21–36. [41] A. Saadatmandi, M. Dehghan, A Legendre collocation method for fractional integro–differential equations, Journal of Vibration and Control, 17(13), (2011) 2050–2058
Pr e-
[42] S. Shen, F. Liu, V. Anh, I. Turner, J. Chen, A novel numerical approximation for the space fractional advection–dispersion equation, IMA Journal of Applied Mathematics, 79 (3), (2014) 431-444. [43] I. Sokolov, A. Chechkin, J. Klafter, Distributed-order fractional kinetics, Acta Physica Polonica B, 35 (2004) 1323-1341. [44] N. Su, P. N. Nelson, S. Connor, The distributed-order fractional diffusion-wave equation of groundwater flow: theory and application to pumping and slug tests, Journal of Hydrology, 529 (2015) 1262-1273.
al
[45] H. Sun, Z. Z. Sun, G. H. Gao, Some high order difference schemes for the space and time fractional Bloch-Torrey equations, Appl. Math. Comput. 281 (2016) 356-380.
urn
[46] Z. Z. Sun, X.N. Wu, A fully discrete difference scheme for a diffusion-wave system, Appl. Numer. Math. 56 (2006) 193-209. [47] T. Tang, A finite difference scheme for a partial integro-differential equations with a weakly singular kernel, Appl. Numer. Math. 11 (1993) 309–319. [48] W. Tian, H. Zhou, W. Deng, A class of second order difference approximations for solving space fractional diffusion equations, Mathematics of Computation 84 (294) (2015) 1703–1727.
Jo
[49] Z. Wang, S. Vong, Compact difference schemes for the modified anomalous fractional sub– diffusion equation and the fractional diffusion–wave equation, J. Comput. Phys., 277 (2014) 1–15. [50] T. Wang, B. Guo, L. Zhang, New conservative difference schemes for a coupled nonlinear Schrodinger system, Appl. Math. Comput., 217 (2010) 1604-1619.
25
Journal Pre-proof
[51] J. Wang, T. Liu, H. Li, Y. Liu, S. He, Second-order approximation scheme combined with H 1 Galerkin MFE method for nonlinear time fractional convection–diffusion equation, Computers & Mathematics with Applications 73 (2017) 1182–1196. [52] Y. Yan, K. Pal, N.. J. Ford, Higher order numerical methods for solving fractional differential equations, BIT Numerical Mathematics, 54 (2014) 555–584.
of
[53] Z. Yang, Z. Yuan, Y. Nie, J. Wang, X. Zhu, F. Liu, Finite element method for nonlinear Riesz space fractional diffusion equations on irregular domains, J. Comput. Phys., 330 (2017) 863–883.
p ro
[54] Y. Yu, W. Deng, Y. Wu, J. Wu, Third order difference schemes (without using points outside of the domain) for one sided space tempered fractional partial differential equations, Appl. Numer. Math., 112 (2017) 126–145.
Pr e-
[55] Q. Yu, F. Liu, I. Turner, K. Burrage, A computationally effective alternating direction method for the space and time fractional Bloch-Torrey equation in 3–D, Applied Mathematics and Computation, 219 (8), (2012) 4082–4095. [56] Q. Yu, F. Liu, I. Turner, K. Burrage, Numerical investigation of three types of space and time fractional Bloch–Torrey equations in 2D, Central European Journal of Physics, 11 (6), (2013) 646–665. [57] Q. Yu, F. Liu, I. Turner, K. Burrage, Numerical simulation of the fractional Bloch equations, Journal of Computational and Applied Mathematics, 255(1), (2014) 635–651.
al
[58] Q. Yu, F. Liu, I. Turner, K. Burrage, Stability and convergence of an implicit numerical method for the space and time fractional Bloch-Torrey equation, Phil. Trans. R. Soc. A 371 (1990), 20120150.
urn
[59] Z. B. Yuan, Y. F. Nie, F. Liu, I. Turner, G. Y. Zhang, Y. T. Gu, An advanced numerical modeling for Riesz space fractional advection-dispersion equations by a meshfree approach, Appl. Math. Model., 40 (2016) 7816-7829. [60] M. Zayernouri, G. E. Karniadakis, Exponentially accurate spectral and spectral element methods for fractional ODEs, J. Comput. Phys., 257 (2014) 460–480.
Jo
[61] M. Zayernouri, G. E. Karniadakis, Discontinuous spectral element methods for time-and spacefractional advection equations, SIAM J. Sci. Comput., 36 (2014) B684–B707. [62] X. Zhao, Z. Z. Sun, Z. P. Hao, A fourth-order compact ADI scheme for two-dimensional nonlinear space fractional Schrodinger equation, SIAM J. Sci. Comput., 36 (2014) A2865A2886. [63] Y. Zhou, Li Peng, On the time–fractional Navier-Stokes equations, Computers & Mathematics with Applications, 73 (2017) 874–891. 26
Journal Pre-proof
[64] Y. Zhou, L. Shangerganesh, J. Manimaran, A. Debbouche, A class of time–fractional reaction– diffusion equation with nonlocal boundary condition, Math. Meth. Appl. Sci. 41 (2018) 2987– 2999.
Jo
urn
al
Pr e-
p ro
of
[65] Y. Zhou, Attractively for fractional evolution equations with almost sectorial operators, Fract. Calc. Appl. Anal. 21.3 (2018) 786–800.
27