- Email: [email protected]
A mixed finite element method for a time-fractional fourth-order partial differential equation
Applied Mathematics and Computation 243 (2014) 703–717
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A mixed finite element method for a time-fractional fourth-order partial differential equation q Yang Liu ⇑, Zhichao Fang, Hong Li *, Siriguleng He School of Mathematical Sciences, Inner Mongolia University, Hohhot 010021, China
a r t i c l e
i n f o
Keywords: Time-fractional fourth-order PDE Mixed finite element method Caputo-fractional derivative Finite difference scheme Stability A priori error estimates
a b s t r a c t In this paper, a numerical theory based on the mixed finite element method for a timefractional fourth-order partial differential equation (PDE) is presented and analyzed. An auxiliary variable r ¼ Du is introduced, then the fourth-order equation can be split into the coupled system of two second-order equations. The time Caputo-fractional derivative is discretized by a finite difference method and the spatial direction is approximated by the mixed finite element method. The stabilities based on a priori analysis for two variables are discussed and some a priori error estimates in L2 -norm for the scalar unknown u and the variable r ¼ Du, are derived, respectively. Moreover, an a priori error result in H1 -norm for the scalar unknown u also is proved. For verifying the theoretical analysis, a numerical test is made by using Matlab procedure. Ó 2014 Elsevier Inc. All rights reserved.
1. Introduction Fractional partial differential equations (FPDEs), whose theoretical analysis and numerical methods have been paid close attention by more and more math researchers, include many types based on the fractional derivative in different directions, such as space FPDEs, time FPDEs and space–time FPDEs. So far, we have found a large number of numerical methods for hunting for the numerical solutions of FPDEs. These methods include finite difference methods [1,4–6,2,9,10,12– 14,16,22,23,30], spectral methods [8], finite element methods [3,11,17–21,24], mixed finite element method [7], finite volume element method [9], DG method [15] and so forth. In recent years, finite element methods for helping people to obtain the numerical solutions for FPDEs have been increasingly concerned by most people. In the recent literatures, we find that the study of mixed finite element methods for FPDEs is very limited. So far, only a paper [7] has been studied and analyzed for a mixed finite element method of time-FPDE with second-order space derivative. However the theoretical analysis of the (mixed) finite element methods for solving the fractional fourth-order PDEs have not been mentioned and reported. In this article, our goal is to give some detailed numerical analysis of a mixed finite element method for studying the following time-FPDE with fourth-order derivative term
@ a uðx; tÞ Du þ D2 u ¼ f ðx; tÞ; ðx; tÞ 2 X J; @t a
ð1:1Þ
q Foundation item: Supported by the National Natural Science Fund (11301258, 11361035), the Natural Science Fund of Inner Mongolia Autonomous Region (2012MS0108, 2012MS0106), the Scientific Research Projection of Higher Schools of Inner Mongolia (NJZZ12011, NJZY13199, and NJZY14013). ⇑ Corresponding authors. E-mail addresses: [email protected] (Y. Liu), [email protected] (H. Li).
http://dx.doi.org/10.1016/j.amc.2014.06.023 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.
704
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
with boundary condition
uðx; tÞ ¼ Duðx; tÞ ¼ 0; ðx; tÞ 2 @ X J;
ð1:2Þ
and initial condition
uðx; 0Þ ¼ u0 ðxÞ; x 2 X;
ð1:3Þ
d
where X R ðd 6 2Þ and J ¼ ð0; T are a bounded convex polygonal domain with Lipschitz continuous boundary @ X and the a uðx;tÞ time interval with 0 < T < 1, respectively. f ðx; tÞ and u0 ðxÞ are given known functions and @ @t is defined by the following a Caputo fractional derivative
@ a uðx; tÞ 1 ¼ @t a Cð1 aÞ
Z 0
t
@uðx; sÞ ds ; 0 < a < 1: @ s ðt sÞa
ð1:4Þ
Here, we approximate the Caputo fractional derivative by a finite difference method, formulate the mixed weak formulation and fully discrete scheme, prove the stability of the fully discrete scheme and derive the theoretical analysis of some a priori error results in detail. The remaining parts of the article is as follow. In Section 2, we introduce a finite difference method for approximating the Caputo time-fractional derivative, and then discuss the detailed proof of the truncation error. In Section 3, we formulate a fully discrete mixed scheme for the fractional fourth-order PDE (1.1) and derive the stable results for two important variables in detail. Moreover, we prove some a priori error estimates in L2 and H1 -norms. In Section 4, we show some numerical results to illustrate the rationality and effectiveness of our method. In Section 5, we make a brief summary about the presented method and the future development. Throughout this paper, we will denote C > 0 as a generic constant free of the space–time step parameters h and d. At the 2 same time, we define the natural inner product in L2 ðXÞ or ðL2 ðXÞÞ by ð; Þ with the corresponding norm k k. The other notations and definitions of Sobolev spaces can be easily followed in Ref. [29]. 2. Approximation of time-fractional derivative For the discretization for time-fractional derivative, let 0 ¼ t 0 < t 1 < t 2 < < t M ¼ T be a given partition of the time interval ½0; T with step length d ¼ T=M and nodes t n ¼ nd, for some positive integer M. For a smooth function / on ½0; T, define /n ¼ /ðtn Þ. Lemma 2.1. Assuming that u 2 C 2 ð½0; TÞ, then the time fractional derivative 0
@ a uðx;tÞ @t a
at t ¼ t nþ1 can be approximated by, for
i u1 u0 @ a uðx; tnþ1 Þ d1a h ¼ ðn þ 1Þ1a ðnÞ1a a @t Cð2 aÞ d n h 1a i kþ1 4uk þ uk1 X d 1a 1a 3u þ ðn k þ 1Þ ðn kÞ þ Enþ1 0 ; Cð2 aÞ k¼1 2d where Enþ1 ¼ E01 þ Enþ1 0 2 ,
Z
t1
"
# 2 t1 þ t 0 @ uðx; t12 Þ ds 2 2 þ Oðð s t ; 1 Þ Þ þ Oðd Þ 2 2 ðt nþ1 sÞa @t 2
E01 ¼
1 Cð1 aÞ
Enþ1 2
" # n Z tkþ1 X 1 @ 2 uðx; tkþ1 Þ ds 2 2 ¼ ðs t kþ1 Þ þ Oððs tkþ1 Þ Þ þ Oðd Þ : Cð1 aÞ k¼1 tk ðt nþ1 sÞa @t2
s
t0
ð2:1Þ
ð2:2Þ
and
ð2:3Þ
Proof. Writing the integral (1.4) into two parts, we get
@ a uðx; tnþ1 Þ 1 ¼ @t a Cð1 aÞ
Z
t1
t0
n Z t kþ1 X @uðx; sÞ ds 1 @uðx; sÞ ds : ¼I þ II: aþ @ s ðtnþ1 sÞ @ s ðtnþ1 sÞa Cð1 aÞ k¼1 tk
ð2:4Þ
Then discretizing parts I and II by the similar approximated schemes to [8,7], respectively, we get the discrete formula (2.1) with truncation errors (2.2) and (2.3). h
705
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
Lemma 2.2. With u 2 C 2 ð½0; TÞ, the truncation error Enþ1 is bounded by 0
! M D d2a C a d2 MD d þ Cd2 1a : þ ln T ¼CðaÞðd2a þ d2 þ ln dÞ; þ Cð3 aÞ Cð1 aÞ Cð2 aÞ
jEnþ1 0 j 6
where M 4 ¼ sup t2½0;T
o n 2 @ uðx;tÞ @t2 ;
ð2:5Þ
l0 ¼ 0; ln ¼ 1; for any n P 1.
Proof. Now we estimate the truncation error E01 . Firstly, we only need to consider the following term. Using the Taylor formula, we have
Z
t1
s
t0
Z t1 Z t1 t1 þ t0 ds t1 þ t0 ds 1a ¼ ðt s Þ d s þ t nþ1 nþ1 a 2 2 ðt nþ1 sÞa t0 t 0 ðt nþ1 sÞ 1 2n þ 1 ¼ d2a ðn2a ðn þ 1Þ2a Þ þ ððn þ 1Þ1a n1a Þ 2a 2ð1 aÞ a 1 2a 2a ¼d ððn þ 1Þ n2a Þ þ ðnðn þ 1Þ1a ðn þ 1Þn1a Þ 2ð1 aÞ 2ð1 aÞð2 aÞ " ! !# 1a 1a a 1 1 1 ¼ d2a n þ n 1þ ðn þ 1Þ n1a ðn þ 1Þ 1 þ n 2ð1 aÞ n 2ð1 aÞð2 aÞ a 1 1 n ¼ d2a ðn þ 1Þ 1 þ ð1 aÞ þ O 2 n n 2ð1 aÞð2 aÞ 1 1 1 nð1 þ ð1 aÞ þ O 2 Þ ðn þ 1Þ n1a þ 2ð1 aÞ n n a 1 1 1 Oð1Þd2a na ¼ d2a 2aþO : þ n1a ¼ a þ O n 2ð1 aÞ n 2ð1 aÞð2 aÞ ð1 aÞð2 aÞ ð2:6Þ
By the simple calculation for (2.2), we have
Z t1 h i @ 2 uðx; t1 Þ Oð1Þd2a na 1 ds 2 O ðs t1 Þ2 þ Oðd2 Þ þ a 2 2 C ð3 a Þ C ð1 a Þ ðt @t t0 nþ1 sÞ h i a Oððth t 1 Þ2 Þ þ Oðd2 Þ na @ 2 uðx; t1 Þ Oð1Þd2a na n 2 2 ¼ n þ ðn þ 1Þ 2 nþ1 Cð3 aÞ Cð2 aÞ @t h i " # 2 2 a 2 Oððth t 1 Þ Þ þ Oðd Þ n @ uðx; t1 Þ Oð1Þd2a na ð1 aÞ 2 2 ; þ ¼ Cð3 aÞ Cð2 aÞ 1 þ a 1n þ Oðn12 Þ @t2
E01 ¼
ð2:7Þ
where t 0 < t h < t1 . Noting that ðnÞa 6 1, we get
jE01 j 6
M D d2a C a d2 ; þ Cð3 aÞ Cð1 aÞ
ð2:8Þ
Secondly, we now start to estimate the error bounds for jEnþ1 2 j. By using the integral calculation, we easily obtain
jEnþ1 2 j 6
n M D d þ Cd2 X Cð1 aÞ k¼1
Z
t kþ1
tk
ds MD d2a þ Cd3a 1a M D d þ Cd2 1a n 6 T : a ¼ Cð2 aÞ Cð2 aÞ ðt nþ1 sÞ
ð2:9Þ
Combining (2.8) and (2.9), we get 0 nþ1 jEnþ1 0 j 6 jE1 j þ ln jE2 j 6
where
MD d2a C a d2 M D d þ Cd2 1a T ¼ CðaÞðd2a þ d2 þ ln dÞ; þ þ ln Cð3 aÞ Cð1 aÞ Cð2 aÞ
l0 ¼ 0; ln ¼ 1; for any n P 1. h
3. Mixed method and the analysis of stability 3.1. Weak formulation and semi-discrete scheme Introducing an auxiliary variable
r ¼ Du, we first rewrite Eq. (1.1) into the following coupled system
ð2:10Þ
706
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
@ a uðx; tÞ Du þ Dr ¼ f ðx; tÞ; @t a
ð3:1Þ
and
r Du ¼ 0:
ð3:2Þ
By using Green’s formula, we get the following mixed weak formulation of (3.1) and (3.2): find fu; rg : ½0; T # such that
H10
H10
a @ uðx; tÞ ; v þ ðru; rv Þ ðrr; rv Þ ¼ ðf ; v Þ; 8v 2 H10 ; a @t
ð3:3Þ
ðr; wÞ þ ðru; rwÞ ¼ 0; 8w 2 H10 :
ð3:4Þ
and
We now define the finite element space V h as
V h ¼ fv h 2 H10 ðXÞ \ C 0 ðXÞjv h je 2 Pm ðeÞ; 8e 2 Kh g;
ð3:5Þ
where Kh is a quasiuniform partition of the domain X. Then the corresponding semi-discrete mixed finite element scheme of (3.3) and (3.4) is to find fuh ; rh g : ½0; T # V h V h such that
a @ uh ðx; tÞ ; v h þ ðruh ; rv h Þ ðrrh ; rv h Þ ¼ ðf ; v h Þ; 8v h 2 V h ; a @t
ð3:6Þ
ðrh ; wh Þ þ ðruh ; rwh Þ ¼ 0; 8wh 2 V h :
ð3:7Þ
and
Remark 3.1 (i) So far, we have not found any related studies on the mixed finite element method for solving fractional fourth-order PDEs. In this article, we derive some detailed theoretical process on a mixed finite element method for solving the fractional fourth-order PDEs (1.1). (ii) We can also get another mixed weak formulation
a @ uðx; tÞ ; v ðr; v Þ ðrr; rv Þ ¼ ðf ; v Þ; 8v 2 H10 ; a @t
ð3:8Þ
ðr; wÞ þ ðru; rwÞ ¼ 0; 8w 2 H10 :
ð3:9Þ
and
(iii) For the system (3.8) and (3.9), we can get some similar results by the similar analysis to that based on the system (3.3) and (3.4), so we omit the corresponding discussions in this article. 3.2. The analysis of stability for discrete scheme In this subsection, we will give the analysis of stability for the case 0 < a < 1. For the convenience of theoretical analysis, we now denote
(
Gak ¼ ðk þ 1Þ
1a
k
1a
and dt ukþ1 ¼
u1 u0 d
;
3ukþ1 4uk þuk1 2d
k¼0 ; k P 1:
Based on the discrete formula (2.1), the system (3.3) and (3.4) can be rewritten as the following equivalent weak formulation
d1a
n X 1 Gank ðdt ukþ1 ; v Þ þ ðrunþ1 ; rv Þ ðrrnþ1 ; rv Þ ¼ ðf nþ1 þ Enþ1 0 ; v Þ; 8v 2 H 0 ;
Cð2 aÞ k¼0
ð3:10Þ
and
ðrnþ1 ; wÞ þ ðrunþ1 ; rwÞ ¼ 0; 8w 2 H10 : Now, a fully discrete finite element procedure is to determine 1a
d
n X
Cð2 aÞ k¼0
ð3:11Þ ðunþ1 h ;
r
nþ1 h Þ
2 V h V h ; ðn ¼ 0; 1; . . . ; M 1Þ such that
nþ1 nþ1 nþ1 Gank ðdt ukþ1 ; v h Þ; 8v h 2 V h ; h ; v h Þ þ ðruh ; rv h Þ ðrrh ; rv h Þ ¼ ðf
ð3:12Þ
707
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
and nþ1 ðrnþ1 h ; wh Þ þ ðruh ; rwh Þ ¼ 0; 8wh 2 V h :
ð3:13Þ
In the next process, we firstly consider the stability for scheme (3.12) and (3.13). Theorem 3.2. For the scheme (3.12) and (3.13), the following stable inequality for unh 2 V h holds, for sufficiently small d
kunh k 6 Cðku0h k þ ka ðdÞmax kf k kÞ;
ð3:14Þ
06k6n
: where ka ðdÞ¼2da Cð2 aÞ; f 2 L2 ðXÞ and C > 0 is a constant free of two mesh parameters d and h. Proof. We take
v h ¼ unþ1 ; h
wh ¼ rnþ1 in (3.12) and (3.13), respectively, and add the two equations to get h
n d1a X nþ1 2 nþ1 2 nþ1 Ga ðd ukþ1 ; unþ1 ; unþ1 h Þ þ kruh k þ krh k ¼ ðf h Þ: Cð2 aÞ k¼0 nk t h
ð3:15Þ
Now we handle the first term on the left hand side of equality (3.15). We easily get the following equality by the simple verification n n1 d1a X d1a d1a X Gank ðdt uhkþ1 ; unþ1 Gan ðdt u1h ; unþ1 Ga ðd unkþ1 ; unþ1 h Þ ¼ h Þþ h Þ Cð2 aÞ k¼0 Cð2 aÞ Cð2 aÞ k¼0 k t h
¼
da
Cð2 aÞ
Gan ðu1h u0h ; unþ1 h Þþ
n1 X da Ga ð3uhnkþ1 4unk þ unk1 ; unþ1 h h h Þ: 2Cð2 aÞ k¼0 k
Noting that Ga0 ¼ 1, it follows that n1 X
Gak 3unkþ1 4unk þ uhnk1 ; unþ1 ¼ h h h
n2 X
3
k¼0
a
Gkþ1 unk h
k¼1
¼
2 3kunþ1 h k
n1 n X X nþ1 4 Gak unk þ Gak1 unk h h ; uh k¼0
a n
a
þ 3G1 4G0
uh ; unþ1 h
ð3:16Þ
!
k¼1
n1 X nþ1 þ ð3Gakþ1 4Gak þ Gak1 Þunk h ; uh
!
3Gan u1h ; unþ1 h
k¼1
þ Gn1 u0h ; unþ1 : h a
ð3:17Þ
: Apply a combination for (3.17), (3.16) and (3.15), multiply (3.15) by ka ðdÞ¼2da Cð2 aÞ and use Cauchy–Schwarz inequality to arrive at 2 nþ1 2 nþ1 2 nþ1 3kunþ1 ;unþ1 h k þ ka ðdÞðkruh k þ krh k Þ ¼ ka ðdÞðf h Þ
n1 X nþ1 ð3Gakþ1 4Gak þ Gak1 Þunk h ;uh
!
k¼1
a 1 nþ1 a a 0 nþ1 þ ð4Ga0 3Ga1 Þðunh ; unþ1 h Þ þ Gn ðuh ;uh Þ þ ð2Gn Gn1 Þðuh ;uh Þ
6
n1 X a a a a a 1 n 0 j3Gakþ1 4Gak þ Gak1 jkunk h k þ j4G0 3G1 jkuh k þ j2Gn Gn1 jkuh k þ Gn kuh k k¼1
þka ðdÞkf nþ1 k kunþ1 h k: 2 unþ1 h k
Remove the second positive term ka ðdÞðkr
nþ1 u 6 h
nþ1 2 k Þ, h
þ kr
ð3:18Þ
then divide kunþ1 h k in inequality (3.18) to obtain
! n1 X 3Ga 4Ga þ Ga kunk k þ j4Ga 3Ga jkun k þ j2Ga Ga jku0 k þ Ga ku1 k þ ka ðdÞkf nþ1 k : h h h h kþ1 k k1 0 1 n n1 n
ð3:19Þ
k¼1
We apply triangle inequality j3Gakþ1 4Gak þ Gak1 j 6 3jGakþ1 Gak j þ jGak Gak1 j; j2Gan Gan1 j 6 jGan1 Gan j þ Gan and note that inequality 1 > Gak1 > Gak > Gakþ1 > 0 in (3.19) to get
kunþ1 h k 6
! n1 n1 X X 1 a a a a a a 1 nk n 0 nþ1 3 ðGak Gakþ1 Þkunk k þ ðG G Þku k þ ð4G 3G Þku k þ G ku k þ G ku k þ k ðdÞkf k : a h h h h h k1 k 0 1 n1 n 3 k¼1 k¼1
ð3:20Þ ku1h k
ku0h k.
Now we derive the relationship between and When n ¼ 0 in (3.12) and (3.13), we take (3.12) and (3.13), respectively, and add the two equations to get
da ðu1 u0h ; u1h Þ þ kr1h k2 þ kru1h k2 ¼ ðf 1 ; u1h Þ: Cð2 aÞ h
vh ¼
u1h
and wh ¼ r1h in
ð3:21Þ
708
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
Multiply (3.21) by da Cð2 aÞ and use Cauchy–Schwarz inequality to get
ku1h k2 þ da Cð2 aÞ kr1h k2 þ kru1h k2 6 ðku0h k þ da Cð2 aÞkf 1 kÞku1h k:
ð3:22Þ
By the simple calculation, we easily get
ku1h k 6 ku0h k þ da Cð2 aÞkf 1 k:
ð3:23Þ
Make a combination for (3.23) and (3.20) to get
! n1 n1 X X 1 a a a a a a a a nk nk n 0 nþ1 1 3 ðGk Gkþ1 Þkuh k þ ðGk1 Gk Þkuh k þ ð4G0 3G1 Þkuh k þ ðGn1 þ Gn Þkuh k þ ka ðdÞkf k þ ka ðdÞkf k : 3 2 k¼1 k¼1
1 kunþ1 h k6
ð3:24Þ From (3.24), we can prove that 0 k kunþ1 h k 6 Cðkuh k þ ka ðdÞ max kf kÞ:
ð3:25Þ
06k6nþ1
In the following discussion, we use mathematical induction to prove the inequality (3.25). When n ¼ 1 in (3.12) and (3.13), we take v h ¼ u2h in (3.12) and (3.13) and sum the two equations to derive
da da Ga1 ðu1h u0h ; u2h Þ þ ð3u2h 4u1h þ u0h ; u2h Þ þ kr2h k2 þ kru2h k2 ¼ ðf 2 ; u2h Þ: Cð2 aÞ 2Cð2 aÞ
ð3:26Þ
Multiply (3.26) by ka ðdÞ and use Cauchy–Schwarz inequality to get
3ku2h k2 6 ka ðdÞðf 2 ; u2h Þ þ 2ðu1h ; u2h Þ þ ðu0h ; u2h Þ 6 ðka ðdÞkf 2 k þ 2ku1h k þ ku0h kÞku2h k: Divide
ku2h k
ð3:27Þ
in (3.27) and substitute (3.23) into (3.27) to get
ku2h k 6
1 1 ðka ðdÞkf 2 k þ 2ku1h k þ ku0h kÞ 6 3ku0h k þ ka ðdÞðkf 2 k þ kf 1 kÞ ¼ C ku0h k þ ka ðdÞmax kf k k : 06k61 3 3
ð3:28Þ
It is not hard to see from (3.28) that (3.25) holds for the case n ¼ 1. We now assume that
kunh k 6 C ku0h k þ ka ðdÞmax kf k k ;
ð3:29Þ
06k6n
holds for the cases n ¼ 1; 2; . . . ; j, then we prove that (3.25) holds based on the induction hypothesis (3.29). We now consider (3.24) with the simple calculation to get for the case n ¼ j þ 1
kujþ1 h k 6
1 ð½3ðGa1 Gaj Þ þ ðGa0 Gaj1 ÞC ku0h k þ ka ðdÞ max kf k k þ ð4Ga0 3Ga1 ÞC ku0h k þ ka ðdÞmaxkf k k 06k6j1 06k6j 3 1 a a þ Gj1 þ Gj ku0h k þ ka ðdÞkf jþ1 k þ ka ðdÞkf 1 kÞ: 2
ð3:30Þ
By the simplification for (3.30), we have for n ¼ j þ 1
0 k kujþ1 h k 6 C kuh k þ ka ðdÞ max kf k :
ð3:31Þ
06k6jþ1
Based on (3.28) and (3.31), we easily see that the conclusion (3.25) of Theorem 3.2 holds. In the next theorem, we will mainly discuss the stability of variable rnh . h Theorem 3.3. For the scheme (3.12) and (3.13), the following stable inequality for
pffiffiffiffiffiffiffiffiffiffiffi krnh k þ krunh k 6 C kr0h k þ kru0h k þ ka ðdÞ max kf k k ;
rnh 2 V h holds, for sufficiently small d
06k6n
ð3:32Þ
where ka ðdÞ is defined as Theorem 3.2 and C is a positive constant free of space–time step h and d. Proof. In order to analyze the stable result for following
rnh , we have to handle the scheme (3.13). Now we rewrite (3.13) as the
n n d1a X d1a X Gank ðdt rkþ1 Ga ðd rukþ1 ; rwh Þ ¼ 0: h ; wh Þ þ Cð2 aÞ k¼0 Cð2 aÞ k¼0 nk t h
Set wh ¼ rnþ1 in (3.33) and h
v h ¼ Cdð2aaÞ 1
Pn
a kþ1 k¼0 Gnkdt uh
in (3.12) and take the sum of the two equations to arrive at
ð3:33Þ
709
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
2 d1a X n n
d1a X a kþ1 kþ1 nþ1 G du þ Ga ðd rkþ1 ; rnþ1 h Þ þ ðdt ruh ; ruh Þ Cð2 aÞ k¼0 nk t h Cð2 aÞ k¼0 nk t h ¼
n X
d1a
Cð2 aÞ k¼0
nþ1 Gank ðdt ukþ1 Þ; h ;f
ð3:34Þ
Use Cauchy–Schwarz inequality and Young inequality to get n
1 nþ1 2 d1a X kþ1 nþ1 Ga ðd rkþ1 ; rnþ1 kf k : h Þ þ ðdt ruh ; ruh Þ 6 2 Cð2 aÞ k¼0 nk t h
ð3:35Þ
Using the similar process of calculation to (3.20), we have 2 jk!nþ1 h kj 6
! n1 n1 X 2 X nk n 0 1 a a a a a a jk!nþ1 3 ðGak Gakþ1 Þjk!nk kj þ ðG G Þjk ! kj þ ð4G 3G Þjk ! kj þ G jk ! kj þ G jk ! kj k1 k 0 1 n1 n h h h h h h kj 3 k¼1 k¼1 1 þ ka ðdÞkf nþ1 k2 ; 3 ð3:36Þ
: j h kj¼k
j hk
ujh k.
where jk! r þ kr Applying the inequalities we have
2 jk!nþ1 h kj 6 C 3
1 2
P
k j¼1 aj
2
6
Pk
2 j¼1 aj
6
P k
j¼1 aj
2
;2
Pk
j;m¼1;j–m aj am
6
P
k j¼1 aj
2
; aj P 0; ðj; m ¼ 0; 1; 2; . . . ; k),
n1 n1 X X pffiffiffiffiffiffiffiffiffiffiffi nþ1 n 0 1 a a a a ka ðdÞkf k ðGak Gakþ1 Þjk!nk ðGak1 Gak Þjk!nk h kj þ h kj þ ð4G0 3G1 Þjk!h kj þ Gn1 jk!h kj þ Gn jk!h kj þ k¼1
!2
k¼1
1 nþ1 2 þ jk!h kj 2 ð3:37Þ
By the simple calculation for (3.37), we easily get jk!
nþ1 h kj 6 C
! n1 n1 X X pffiffiffiffiffiffiffiffiffiffiffi nþ1 nk nk n 0 1 a a a a a a a a 3 ðGk Gkþ1 Þjk!h kj þ ðGk1 Gk Þjk!h kj þ ð4G0 3G1 Þjk!h kj þ Gn1 jk!h kj þ Gn jk!h kj þ ka ðdÞkf k : k¼1
k¼1
ð3:38Þ
In the next consideration, we use a similar process of proof to (3.25) based on mathematical induction to get 0 jk!nþ1 h kj 6 Cðjk!h kj þ
pffiffiffiffiffiffiffiffiffiffiffi ka ðdÞ max kf k kÞ:
ð3:39Þ
06k6n
Noting that the notation for jk!jh kj, we get the stability on
rnh . h
3.3. Fully discrete error analysis In this subsection, we will give some fully discrete analysis of a priori error estimates for the case 0 < a < 1. Now a Ritzprojection associated with our equations is introduced to help us to analyze the convergence of the method. Lemma 3.4. Define a Ritz-projection operator {h : H10 ðXÞ ! V h by
ðrðu {h uÞ; rwh Þ ¼ 0; 8wh 2 V h ;
ð3:40Þ
and holds
ku {h uk þ hku {h uk1 6 Ch
mþ1
kukmþ1 ; 8u 2 H10 ðXÞ \ Hmþ1 ðXÞ;
P
R
1 2
ð3:41Þ
where the norms is defined by kzkk ¼ ½ 06jbj6k X jDb zj2 dx with the degree k of the polynomial. The contents of Lemma 3.4 can be found in many literatures, such as Refs. [29,28]. In the next analysis, we will discuss the proof of fully discrete a priori error estimates based on the Ritz-projection Lemma 3.4 in detail. To simplify the process of writing in the proof, we now rewrite the errors as
uðtn Þ unh ¼ ðuðt n Þ {h un Þ þ ð{h un unh Þ ¼ fn þ /n ;
rðtn Þ rnh ¼ ðrðtn Þ {h rn Þ þ ð{h rn rnh Þ ¼ pn þ vn :
710
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
Theorem 3.5. Supposing that un 2 H10 ðXÞ \ Hmþ1 ðXÞ; unh 2 V h and u0h ¼ {h uð0Þ, then the error estimate for u in L2 -norm holds with a parameter 0 < a < 1
kuðt n Þ unh k 6
Cðu; aÞT a mþ1 mþ1 ðln d þ d2 þ d2a þ ð1 þ da Þh Þ þ CðuÞh ; 1a
ð3:42Þ
and when a ! 1 mþ1
kuðt n Þ unh k 6 CðuÞðln d þ d2 þ d þ ð1 þ d1 Þh
Þ;
ð3:43Þ
where both Cðu; aÞ (with a) and CðuÞ depend on the norm kukmþ1 . Proof. Subtract (3.12) from (3.10) and apply the Ritz-projection (3.40) to obtain the error equations n n d1a X d1a X Gank ðdt /kþ1 ; v h Þ þ ðr/nþ1 ; rv h Þ ðrvnþ1 ; rv h Þ ¼ Ga ðd fkþ1 ; v h Þ þ ðEnþ1 0 ; v h Þ; 8v h 2 V h ; Cð2 aÞ k¼0 Cð2 aÞ k¼0 nk t
ð3:44Þ
ðvnþ1 ; wh Þ þ ðr/nþ1 ; rwh Þ ¼ ðpnþ1 ; wh Þ; 8wh 2 V h :
ð3:45Þ
and
For the first term on the left hand side of (3.44), we take a similar process to (3.16) to get n n1 X d1a X da da Gank ðdt /kþ1 ; v h Þ ¼ Gan ð/1 /0 ; v h Þ þ Ga ð3/nkþ1 4/nk þ /nk1 ; v h Þ: Cð2 aÞ k¼0 Cð2 aÞ 2Cð2 aÞ k¼0 k
ð3:46Þ
We add (3.44) and (3.45), substitute (3.46), take ðv h ; wh Þ ¼ ð/nþ1 ; vnþ1 Þ and multiply by ka ðdÞ ¼ 2da Cð2 aÞ to arrive at n1 X Gak ð3/nkþ1 4/nk þ /nk1 ; /nþ1 Þ þ ka ðdÞðkvnþ1 k2 þ kr/nþ1 k2 Þ k¼0 n X nþ1 ¼ 2d Gank ðdt fkþ1 ; /nþ1 Þ 2Gan ð/1 /0 ; /nþ1 Þ þ ka ðdÞðEnþ1 Þ ka ðpnþ1 ; vnþ1 Þ: 0 ;/
ð3:47Þ
k¼0
Now we deal with the first term on the left hand side of (3.47). Noting that Ga0 ¼ 1, we use the similar process to (3.17) to arrive at n1 X Gak ð3/nkþ1 4/nk þ /nk1 ; /nþ1 ¼ 3k/nþ1 k2 þ ð3Ga1 4Ga0 Þð/n ; /nþ1 Þ k¼0
þ
! n1 X ð3Gakþ1 4Gak þ Gak1 Þ/nk ; /nþ1 3Gan ð/1 ; /nþ1 Þ k¼1
þ Gan1 ð/0 ; /nþ1 Þ:
ð3:48Þ
By substituting (3.48) into (3.47), we obtain n X nþ1 3k/nþ1 k2 þ ka ðdÞðkvnþ1 k2 þ kr/nþ1 k2 Þ ¼ 2d Gank ðdt fkþ1 ; /nþ1 Þ þ ka ðdÞðEnþ1 Þ 0 ;/ k¼0
ð3Ga1 4Ga0 Þð/n ; /nþ1 Þ
n1 X
! ð3Gakþ1 4Gak þ Gak1 Þ/nk ; /nþ1
k¼1
þ Gan ð/1 ; /nþ1 Þ þ ð2Gan Gan1 Þð/0 ; /nþ1 Þ ka ðpnþ1 ; vnþ1 Þ:
ð3:49Þ
For the next proof, we have to deal with the first term on the right hand side of (3.49). Now we use the similar proof to (3.46) and (3.48) to obtain n n1 X X 2d Gank ðdt fkþ1 ; /nþ1 Þ ¼ 2Gan ðf1 f0 ; /nþ1 Þ Gak ð3fnkþ1 4fnk þ fnk1 ; /nþ1 Þ k¼0
k¼0 n1 X ð3Gakþ1 4Gak þ Gak1 Þfnk ; /nþ1 ¼ 3ðfnþ1 ; /nþ1 Þ
! ð3Ga1 4Ga0 Þðfn ; /nþ1 Þ
k¼1
þ Gan ðf1 ; /nþ1 Þ þ ð2Gan Gan1 Þðf0 ; /nþ1 Þ:
ð3:50Þ
Applying (3.41), Cauchy–Schwarz inequality and the inequality
j3Gakþ1 4Gak þ Gak1 j 6 3jGakþ1 Gak j þ jGak Gak1 j ¼ 3ðGak Gakþ1 Þ þ ðGak1 Gak Þ;
ð3:51Þ
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
711
we have n n1 X X 2d Gank ðdt fkþ1 ; /nþ1 Þ ¼ ð3kfnþ1 k þ j3Gakþ1 4Gak þ Gak1 jkfnk k þ ð4Ga0 3Ga1 Þkfn k þ Gan kf1 k k¼0
k¼1
þ Gan1 kf0 kÞk/nþ1 k 6 CðuÞð3 þ 3ðGa1 Gan Þ þ ðGa0 Gan1 Þ þ ð4Ga0 3Ga1 Þ þ Gan þ Gan1 Þh
mþ1
k/nþ1 k ¼ Cðu; aÞh
mþ1
k/nþ1 k:
ð3:52Þ
Substitute (3.52) into (3.49) and use Cauchy–Schwarz inequality to get
3k/nþ1 k2 þ ka ðdÞðkvnþ1 k2 þ kr/nþ1 k2 Þ 6 ðj4Ga0 3Ga1 jk/n k n1 X
þ
j3Gakþ1 4Gak þ Gak1 jk/nk k þ Gan k/1 k þ ðjGan1 Gan j þ Gan Þk/0 k
k¼1
!
þ ka ðdÞkEnþ1 0 k þ Cðu; aÞh
mþ1
k/nþ1 k þ ka kpnþ1 kkvnþ1 k:
ð3:53Þ
Noting that the following inequalities
1 2 2 2 ða þ bÞ 6 ka þ b ; 2
ac þ bd 6 aðc þ dÞ þ bðc þ dÞ ¼ ða þ bÞðc þ dÞ;
ð3:54Þ
for any a; b; c; d > 0; k 2 N þ , we have
pffiffiffiffiffiffiffiffiffiffiffi 2 1 ðk/nþ1 k þ ka ðdÞkvnþ1 kÞ 6 ðj4Ga0 3Ga1 jk/n k 2 n1 X j3Gakþ1 4Gak þ Gak1 jk/nk k þ Gan k/1 k þ ðjGan1 Gan j þ Gan Þk/0 k þ k¼1
þ
ka ðdÞkEnþ1 0 k
þ Cðu; aÞh
mþ1
! pffiffiffiffiffi pffiffiffiffiffi nþ1 þ ka kp k ðk/nþ1 k þ ka kvnþ1 kÞ:
ð3:55Þ
Use inequality (3.51) and the inequality 1 > Gak1 > Gak > Gakþ1 > 0 to get
k/nþ1 k 6 2ðð4Ga0 3Ga1 Þk/n k þ ½3ðGan1 Gan Þ þ Ga0 Ga1 þ 4
n1 X ðGak1 Gak Þk/nk k þ Gan k/1 k þ Gan1 k/0 k k¼2
þ ka ðdÞkEnþ1 0 k þ Cðu; aÞh
mþ1
Þ:
ð3:56Þ
Based on (3.56), we can derive
k/nþ1 k 6
Cðu; aÞ mþ1 ðln d1þa þ d2þa þ d2 þ h Þ: Gan
ð3:57Þ
In the following analysis and discussion, we use mathematical induction to discuss the conclusion (3.57). When n ¼ 0 in (3.44) and (3.45), taking v h ¼ /1 and wh ¼ v1 in (3.44) and (3.45), respectively, and summing the two equations, we use Cauchy–Schwarz inequality to obtain by the simplification
k/1 k2 k/0 k2 1 1 þ k/ /0 k2 þ da Cð2 aÞðkv1 k2 þ kr/1 k2 Þ 2 2 ¼ ðf1 f0 ; /1 Þ þ da Cð2 aÞðE10 ; /1 Þ da Cð2 aÞðp1 ; v1 Þ 6 ðkf1 k þ kf0 k þ CðaÞðd2þa þ d2 ÞCð2 aÞÞk/1 k þ da Cð2 aÞkp1 kkv1 k:
ð3:58Þ
Noting that /0 ¼ 0 and using two important inequalities (3.54), we use the similar discussion to (3.55) to obtain
k/1 k 6 Cðu; aÞðh
mþ1
þ d2þa þ d2 Þ:
ð3:59Þ
Based on the case n ¼ 0, we now discuss the case n ¼ 1, which is the first step of mathematical induction. Under the condition of n ¼ 1, setting v h ¼ /2 in (3.44) and wh ¼ v2 in (3.45), and taking the sum of the two equations, we arrive at
1 ðk/2 k2 þ k2/2 /1 k2 þ k/2 2/1 þ /0 k2 Þ þ ka ðdÞðkv2 k2 þ kr/2 k2 Þ 2 1 ¼ ðk/1 k2 þ k2/1 /0 k2 Þ Ga1 ð/1 /0 ; /2 Þ Ga1 ðf1 f0 ; /2 Þ þ ka ðdÞðE20 ; /2 Þ ka ðdÞðp2 ; v2 Þ: 2 By an application of Cauchy–Schwarz inequality and Young inequality, we can get
ð3:60Þ
712
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
k/2 k2 þ k2/2 /1 k2 þ k/2 2/1 þ /0 k2 6 CðaÞðk/1 k2 þ k/0 k2 þ kf1 k2 þ kf0 k2 þ ka ðdÞkp2 k2 þ ln d2þ2a þ d4þ2a þ d4 Þ:
ð3:61Þ
Combining (3.59) and (3.41), we arrive at
k/2 k 6 Cðu; aÞðh ¼
mþ1
þ ln d1þa þ d2þa þ d2 Þ ¼ ð21a 1Þ
Cðu; aÞ mþ1 ðh þ ln d1þa þ d2þa þ d2 Þ Ga1
Cðu; aÞ mþ1 ðh þ ln d1þa þ d2þa þ d2 Þ; Ga1
ð3:62Þ
which is to tell that (3.57) holds for the first case n ¼ 1. We firstly assume that (3.57) holds, for 1 6 n 6 j
k/jþ1 k 6
Cðu; aÞ mþ1 ðln d1þa þ d2þa þ d2 þ h Þ: Gaj
ð3:63Þ
Now, we consider the case n ¼ j þ 1. By (3.56) and the induction hypothesis (3.63), we arrive at n1 X 1 mþ1 ðjþ1Þþ1 k/ðjþ1Þþ1 k 6 ðð4Ga0 3Ga1 Þk/jþ1 k þ ½3ðGan1 Gan Þ þ Ga0 Ga1 þ 4 ðGak1 Gak Þk/ðjþ1Þk k þ Gan k/1 k þ Gan1 k/0 k þ ka ðdÞkE0 k þ Cðu; aÞh Þ 3 k¼2
6 Cðu; aÞðð4Ga0 3Ga1 Þ
Noting that
1 Gajk
ðjþ1Þþ1
k/
1 1 mþ1 þ ½3ðGan1 Gan Þ þ Ga0 Ga1 þ 4ðGa1 Gan1 Þ a þ Gan þ 2Þðln d1þa þ d2þa þ d2 þ h Þ: Gaj Gjk
ð3:64Þ
< G1a < Ga1 and 0 < Gajþ1 < Gaj < 1, we have j
jþ1
mþ1
k 6 Cðu; aÞð4Ga0 3Ga1 þ 3ðGan1 Gan Þ þ Ga0 Ga1 þ 4ðGa1 Gan1 Þ þ Gan þ 2Þðln d1þa þ d2þa þ d2 þ h 6
Cðu; aÞ mþ1 ðln d1þa þ d2þa þ d2 þ h Þ: Gajþ1
Þ
1 Gajþ1 ð3:65Þ
that is to say that (3.57) holds for n ¼ j þ 1. By using induction based on (3.62) and (3.65), we announce that (3.57) holds for any n P 1. a In order to get the final result based on (3.65), we have to give the estimate for ðnþ1Þ . We now consider Taylor formula to Gan get
ðn þ 1Þa ðn þ 1Þa 1 1 ¼ ¼ a a ¼ a2 2 Gn ðn þ 1Þ1a n1a ðn þ 1Þ nð1 þ 1nÞ ðn þ 1Þ nð1 þ a 1n þ aða2!1Þ ð1 þ h 1nÞ ð1nÞ Þ ¼
1 1aþ
að1aÞ 2!
ð1 þ
a2 h 1nÞ 1nÞ
<
1 ; 1a
ð3:66Þ
where 0 < h < 1 and 0 < a < 1. Noting that the result of (3.66) and ðn þ 1Þd 6 Md ¼ T, we have
Cðu; aÞ Cðu; aÞðn þ 1Þa Cðu; aÞT a da ¼ ðn þ 1Þa da da 6 : a a 1a Gn Gn
ð3:67Þ
Combining (3.57) and (3.67), we get
k/nþ1 k 6
Cðu; aÞT a mþ1 ðln d þ d2 þ d2a þ da h Þ: 1a
ð3:68Þ
By a combination for (3.41) and (3.68) with triangle inequality, we get the error result (3.42) of theorem. Now we briefly discuss the case a ! 1. We note that the error result (3.42) in this case has no meaning since the aÞT a coefficient Cðu; 1a ! 1 as a ! 1. So, we have to seek another error estimate’s form. Considering the fact nd 6 T ðn ¼ 1; 2; . . . ; MÞ, we can get the error result
k/n k 6 CðuÞndðln d þ d2 þ d þ d1 h
mþ1
Þ:
ð3:69Þ
Now we can use induction to derive the proof. The detail process is similar to the above proof, so we do not derive that again. Combining (3.69) and (3.41) with triangle inequality, we can arrive at the estimate (3.43). In the next analysis, our aim is to focus on deriving a priori error estimates in L2 -norm for variable r. At the same time, we also analyze the error result in H1 -norm for u. h Theorem 3.6. With un ; rn 2 H10 ðXÞ \ Hmþ1 ðXÞ; unh ; estimates with 0 < a < 1
rnh 2 V h ; u0h ¼ {h uð0Þ and r0h ¼ {h rð0Þ, one has the following a priori error
713
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
krðt n Þ rnh k 6 kuðt n Þ unh k1 6
3a mþ1 Cðu; r; aÞT a mþ1 mþ1 ðln d þ d2 þ d2a þ d 2 h þ h Þ þ CðrÞh ; 1a
ð3:70Þ
3a mþ1 Cðu; r; aÞT a mþ1 m ðln d þ d2 þ d2a þ d 2 h þ h Þ þ CðuÞh ; 1a
ð3:71Þ
and when a ! 1 3
krðt n Þ rnh k 6 Cðu; rÞTðln d þ d2 þ d þ d2 h 3
mþ1
kuðt n Þ unh k1 6 Cðu; rÞTðln d þ d2 þ d þ d2 h
mþ1
þh
mþ1
þh
mþ1
Þ þ CðrÞh
mþ1
ð3:72Þ
;
m
Þ þ CðuÞh :
ð3:73Þ
where Cðu; r; aÞ (with a) and Cðu; rÞ are dependent on the norms kukmþ1 and krkmþ1 ; CðrÞ only depends on krkmþ1 and CðuÞ is described as the statement in Theorem 3.5. Proof. For the need of proof, we have to rewrite (3.13) as n n n d1a X d1a X d1a X Gank ðdt vkþ1 ; wh Þ þ Gank ðrdt /kþ1 ; rwh Þ ¼ Ga ðd pkþ1 ; wh Þ: Cð2 aÞ k¼0 Cð2 aÞ k¼0 Cð2 aÞ k¼0 nk t 1a
We add (3.74) and (3.12), take wh ¼ vnþ1 and v h ¼ Cdð2aÞ ity to get
Pn
a kþ1 , k¼0 Gnkdt /
ð3:74Þ
multiply by ka ðdÞ and use Cauchy–Schwarz inequal-
2 d1a X n n X ka ðdÞ Gankdt /kþ1 þ 2d Gank ½ðdt vkþ1 ; vnþ1 Þ þ ðrdt /kþ1 ; r/kþ1 Þ Cð2 aÞ k¼0 k¼0 ! n n n X a X d1a X a kþ1 nþ1 kþ1 2d Gank ðdt pkþ1 ; vnþ1 Þ Gnkdt / ¼ 2d Gnkdt f þ ka ðdÞE0 ; Cð2 aÞ k¼0 k¼0 k¼0 ! X d1a X n n n X Ga d /kþ1 þ k2d Gankdt pkþ1 kkvnþ1 k: 6 2d Gankdt fkþ1 þ ka ðdÞkEnþ1 0 k k¼0 Cð2 aÞ k¼0 nk t k¼0
ð3:75Þ
Using Young inequality and the first inequality in (3.54), we have
2 d1a X n n X 1 ka ðdÞ Gankdt /kþ1 þ 2d Gank ½ðdt vkþ1 ; vnþ1 Þ þ ðrdt /kþ1 ; r/kþ1 Þ Cð2 aÞ k¼0 2 k¼0 2 2 X n n 2 1 X a kþ1 2 d k þ 2d G p 6 ðka ðdÞÞ1 2d Gankdt fkþ1 þ ka ðdÞkEnþ1 þ e vnþ1 : nk t 0 k¼0 4e k¼0
ð3:76Þ
2
u
1 0 1
−1 −2 1
0.8 0.6 0.8
0.6
0.4 0.4
0.2
0.2 0
0
Fig. 1. Surface for exact solution u.
714
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
Table 1 The errors and convergence orders for u and
a
M
pffiffiffi
r with P1 element and space–time mesh d ¼ h= 2 ¼ 1=M. Order
L2 -norm u
Order
H1 -norm u
L2 -norm
Order
r
0.3
16 32 64
6.7904e002 1.7450e002 4.3932e003
1.9603 1.9899
1.6779e+000 8.3477e001 4.1677e001
1.0072 1.0021
2.9915e+000 7.5426e001 1.8901e001
1.9877 1.9966
0.5
16 32 64
6.7903e002 1.7450e002 4.3932e003
1.9603 1.9899
1.6779e+000 8.3477e001 4.1677e001
1.0072 1.0021
2.9915e+000 7.5424e001 1.8900e001
1.9878 1.9966
0.7
16 32 64
6.7902e002 1.7450e002 4.3930e003
1.9603 1.9898
1.6779e+000 8.3477e001 4.1677e001
1.0072 1.0021
2.9914e+000 7.5421e001 1.8899e001
1.9878 1.9967
0.9
16 32 64
6.7900e002 1.7448e002 4.3924e003
1.9603 1.9900
1.6779e+000 8.3476e001 4.1677e001
1.0072 1.0021
2.9912e+000 7.5413e001 1.8895e001
1.9879 1.9968
2
u
h
1 0 1
−1 −2 1
0.8 0.6 0.8
0.6
0.4 0.4
0.2
0.2 0
0
Fig. 2. Surface for numerical solution uh .
Now we discuss the results for any n P 0. Firstly, taking n ¼ 0 and
e ¼ 12 in (3.76), we have
2 da 1 1 a 1 a a 1 1 2 0 1 0 1 0 2 1 0 1 ka ðdÞ G ð/ / Þ Cð2 aÞ 0 þ 2G0 ½ðv v ; v Þ þ ðr/ r/ ; r/ Þ 6 ðka ðdÞÞ k2G0 ðf f Þk þ ka ðdÞkE0 k 2 1 1 þ k2Ga0 ðp1 p0 Þk2 þ kv1 k2 : 2 2 Noting that Ga0 ¼ 1;
ð3:77Þ
v0 ¼ 0; /0 ¼ 0 and (3.41), we have by the simple simplification
2 da 1 3 1 0 1 2 þ 2 r/1 2 6 CðuÞðka ðdÞÞ1 h2þ2m þ CðaÞka ðdÞðd42a þ d4 Þ þ CðrÞh2þ2m : ka ðdÞ / Þ ð/ Cð2 aÞ þ2 v 2
ð3:78Þ
From (3.78), we easily find that
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi da Cð2 aÞ
1 da 1þm 1 a2 1þm þ d2 þ d2þa þ h Þ: ð/1 /0 Þ þ v þ r/ 6 Cðu; r; aÞðd h Cð2 aÞ
For the n P 1 in (3.76), we use mathematical induction and a similar discussion to Theorem 3.5 to get
ð3:79Þ
715
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
200
σ
100 0 1
−100 −200 1
0.8 0.6 0.8
0.4
0.6
0.4
0.2
0.2 0
0
Fig. 3. Surface for exact solution
r.
200
σh
100 0 1
−100 −200 1
0.8 0.6 0.8
0.6
0.4 0.4
0.2
0.2 0
Fig. 4. Surface for numerical solution
0 rh .
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 1a X 3a 1þm Cðu; r; aÞT a d 1þm a kþ1 a d Cð2 aÞ Gnk dt / þ kvn k þ kr/n k 6 ðln d þ d2 þ d2a þ d 2 h þ h Þ: Cð2 aÞ k¼0 1a
ð3:80Þ
Combining (3.80) and (3.41) with triangle inequality, we get the estimates (3.70) and (3.71) for the case 0 < a < 1. For the case a ! 1, use a similar analysis as that in Theorem 3.5 to get the error estimates (3.72) and (3.73). h Remark 3.7. From the error results in Theorems 3.5 and 3.6, we can see that the results of convergence in the spatial direction is affected by the inverse of the time mesh d. However the time step’s inverse would not affect the global accuracy, mþ1 because h can be much smaller than d for sufficiently regular solutions on the spatial variable. The similar discussions can be found in the Remark 4.2 of Lin and Xu Ref. [8].
716
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
4. Some numerical results In this section, we provide a numerical test to confirm the theoretical analysis of our method. Here, we choose h i 2 2t 2a 2 4 Cð3aÞ þ ð8p þ 64p Þðt þ 1Þ sinð2pxÞ sinð2pyÞ; 8ðx; yÞ 2 ½0; 1 ½0; 1; t 2 ½0; 1 and initial value function
function f ðx; tÞ ¼
u0 ðx; y; 0Þ ¼ sinð2pxÞ sinð2pyÞ. We now get the exact solution
uðx; y; tÞ ¼ ðt2 þ 1Þ sinð2pxÞ sinð2pyÞ;
ð4:1Þ
r ¼ Du ¼ 8p2 ðt2 þ 1Þ sinð2pxÞ sinð2pyÞ:
ð4:2Þ
and
We now take the finite element space V h as the linear space with continuous piecewise linear functions, and divide the qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 1 2ffi pffiffiffi 1 2 2 spatial domain ½0; 1 ½0; 1 into uniform triangulations with mesh size h ¼ hx þ hy ¼ þ N ¼ 2 N. For formulating N fully discrete procedure, we now consider the discrete time step d ¼ M1 for the time interval ½0; 1. For the convenience of calculation, we now take M ¼ N ¼ 16; 32; 64, respectively, then compute some a priori errors and orders of convergence of u in both L2 and H1 -norms for different a ¼ 0:3; 0:5; 0:7; 0:9. From the results of convergence in Table 1, we can see that the rates of convergence of u in L2 and H1 -norms are close to 2 and 1, respectively, which are in mþ1 accord with the spatial convergence order Oðh Þ with index m ¼ 1. Moreover, the error in L2 -norm for the variable r is also computed and the rate of convergence is close to 2. In Fig. 1, the surface of the exact solution u is shown. With the chosen pffiffiffi pffiffiffi space–time steps h ¼ 2d ¼ 2=16, the surface of the numerical solution uh is displayed at time t ¼ 1 in Fig. 2. It not hard to see that the numerical uh can approximate well the exact solution u. The similar comparison between rh and r is also made in Figs. 3 and 4. Based on the results of errors and rates of convergence, we know our method is feasible and effective. 5. Some remarks and extensions So far, we have not found any reports on the (mixed) finite element method for solving the fractional fourth-order PDE (1.1). In this article, our aim is propose a mixed finite element method for the fractional fourth-order PDE with the case 0 < a < 1. Based on a finite difference method for time fractional derivative, we give a fully discrete mixed finite element scheme, and then analyze the detailed proofs’ process for both the stabilities in L2 -norm and some a priori error results in L2 and H1 -norms. Finally, we present a numerical example to show the feasibility of the studied method. For the case 1 < a < 2, we will give the numerical analysis and simulation based on the mixed finite element method in another work. In the near future, we will discuss a moving (mixed) finite element schemes based on the discussion in [19,21] and obtain the optimal convergence analysis of mixed finite element method for fractional fourth-order PDE by using a special interpolation operator in Ref. [25]. At the same time, we will apply some other (mixed) finite element methods [26,27] to analyze the fractional fourth-order PDE. Acknowledgments Authors thank the reviewers and editor for their valuable comments and suggestions. This work is supported by the National Natural Science Fund (11301258, 11361035), the Natural Science Fund of Inner Mongolia Autonomous Region (2012MS0108, 2012MS0106), the Scientific Research Projection of Higher Schools of Inner Mongolia (NJZZ12011, NJZY13199, and NJZY14013). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.amc.2014.06.023. References [1] D. Baleanu, A.H. Bhrawy, T.M. Taha, Two efficient generalized Laguerre spectral algorithms for fractional initial value problems, Abstract Appl. Anal. 2013 (2013). Article ID 546502, 10 pages. [2] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for two-sided space-fractional partial differential equations, Appl. Math. Comput. 56 (1) (2006) 80–90. [3] N. Zhang, W. Deng, Y. Wu, Finite difference/element method for a two-dimensional modified fractional diffusion equation, Adv. Appl. Math. Mech. 4 (2012) 496–518. [4] C. Tadjeran, M. Meerschaert, A second-order accurate numerical method for the two-dimensional fractional diffusion equation, J. Comput. Phys. 220 (2007) 813–823. [5] C. Tadjeran, M. Meerschaert, H. Scheffler, A second-order accurate numerical approximation for the fractional diffusion equation, J. Comput. Phys. 213 (2006) 205–213. [6] C.M. Chen, F. Liu, K. Burrage, Finite difference methods and a Fourier analysis for the fractional reaction-subdiffusion equation, Appl. Math. Comput. 198 (2) (2008) 754–769.
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
717
[7] Y. Liu, H. Li, W. Gao, S. He, Z.C. Fang, A new mixed element method for a class of time-fractional partial differential equations, Sci. World J. 2014 (2014). Article ID 141467, 8 pages. [8] Y.M. Lin, C.J. Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, J. Comput. Phys. 225 (2) (2007) 1533–1552. [9] P. Zhuang, F. Liu, I. Turner, Y.T. Gu, Finite volume and finite element methods for solving a one-dimensional space-fractional Boussinesq equation, Appl. Math. Model. (2013). in press, http://dx.doi.org/10.1016/j.apm.2013.10.008. [10] C.M. Chen, F. Liu, V. Anh, I. Turner, Numerical methods for solving a two-dimensional variable-order anomalous subdiffusion equation, Math. Comput. 81 (2012) 345–366. [11] H. Zhang, F. Liu, V. Anh, Galerkin finite element approximation of symmetric space-fractional partial differential equations, Appl. Math. Comput. 217 (6) (2010) 2534–2545. [12] R. Lin, F. Liu, V. Anh, I. Turner, Stability and convergence of a new explicit finite-difference approximation for the variable-order nonlinear fractional diffusion equation, Appl. Math. Comput. 212 (2) (2009) 435–445. [13] S.B. Yuste, Weighted average finite difference methods for fractional diffusion equations, J. Comput. Phys. 216 (1) (2006) 264–274. [14] E. Sousa, Finite difference approximations for a fractional advection diffusion problem, J. Comput. Phys. 228 (2009) 4038–4054. [15] L.L. Wei, Y.N. He, Analysis of a fully discrete local discontinuous Galerkin method for time-fractional fourth-order problems, Appl. Math. Model. 38 (4) (2014) 1511–1522. [16] C.P. Li, F.H. Zeng, Finite difference methods for fractional differential equations, Int. J. Bifurcation Chaos 22 (2012). 28 pages. [17] C.P. Li, Z.G. Zhao, Y.Q. Chen, Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion, Comput. Math. Appl. 62 (3) (2011) 855–875. [18] Z.G. Zhao, C.P. Li, Fractional difference/finite element approximations for the time-space fractional telegraph equation, Appl. Math. Comput. 219 (6) (2012) 2975–2988. [19] J.T. Ma, J.Q. Liu, Z.Q. Zhou, Convergence analysis of moving finite element methods for space fractional differential equations, J. Comput. Appl. Math. 255 (2014) 661–670. [20] Y.J. Jiang, J.T. Ma, High-order finite element methods for time-fractional partial differential equations, J. Comput. Appl. Math. 235 (11) (2011) 3285– 3290. [21] Y.J. Jiang, J.T. Ma, Moving finite element methods for time fractional partial differential equations, Sci. China Math. 56 (6) (2013) 1287–1300. [22] G. Gao, Z. Sun, A compact finite difference scheme for the fractional sub-diffusion equations, J. Comput. Phys. 230 (2011) 586–595. [23] Y. Zhang, Z. Sun, H. Wu, Error estimate of Crank–Nicolson-type difference for subdiffusion equation, SIAM J. Numer. Anal. 49 (2011) 2302–2322. [24] N.J. Ford, J.Y. Xiao, Y.B. Yan, A finite element method for time fractional partial differential equations, Fract. Calc. Appl. Anal. 14 (3) (2011) 454–474. [25] J.C. Li, Optimal convergence analysis of mixed finite element methods for fourth-order elliptic and parabolic problems, Numer. Methods Partial Differ. Equ. 22 (2006) 884–896. [26] D.Y. Shi, S.P. Mao, S.C. Chen, On the anisotropic accuracy analysis of ACM’s nonconforming finite element, J. Comput. Math. 23 (6) (2005) 635–646. [27] Y. Liu, Z.C. Fang, H. Li, S. He, W. Gao, A coupling method based on new MFE and FE for fourth-order parabolic equation, J. Appl. Math. Comput. 43 (2013) 249–269. [28] Z.D. Luo, Mixed Finite Element Methods and Applications, Chinese Science Press, Beijing, China, 2006. [29] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, The Netherlands, 1978. [30] K.X. Wang, H. Wang, A fast characteristic finite difference method for fractional advection-diffusion equations, Adv. Water Resour. 34 (7) (2011) 810–816.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A mixed finite element method for a time-fractional fourth-order partial differential equation q Yang Liu ⇑, Zhichao Fang, Hong Li *, Siriguleng He School of Mathematical Sciences, Inner Mongolia University, Hohhot 010021, China
a r t i c l e
i n f o
Keywords: Time-fractional fourth-order PDE Mixed finite element method Caputo-fractional derivative Finite difference scheme Stability A priori error estimates
a b s t r a c t In this paper, a numerical theory based on the mixed finite element method for a timefractional fourth-order partial differential equation (PDE) is presented and analyzed. An auxiliary variable r ¼ Du is introduced, then the fourth-order equation can be split into the coupled system of two second-order equations. The time Caputo-fractional derivative is discretized by a finite difference method and the spatial direction is approximated by the mixed finite element method. The stabilities based on a priori analysis for two variables are discussed and some a priori error estimates in L2 -norm for the scalar unknown u and the variable r ¼ Du, are derived, respectively. Moreover, an a priori error result in H1 -norm for the scalar unknown u also is proved. For verifying the theoretical analysis, a numerical test is made by using Matlab procedure. Ó 2014 Elsevier Inc. All rights reserved.
1. Introduction Fractional partial differential equations (FPDEs), whose theoretical analysis and numerical methods have been paid close attention by more and more math researchers, include many types based on the fractional derivative in different directions, such as space FPDEs, time FPDEs and space–time FPDEs. So far, we have found a large number of numerical methods for hunting for the numerical solutions of FPDEs. These methods include finite difference methods [1,4–6,2,9,10,12– 14,16,22,23,30], spectral methods [8], finite element methods [3,11,17–21,24], mixed finite element method [7], finite volume element method [9], DG method [15] and so forth. In recent years, finite element methods for helping people to obtain the numerical solutions for FPDEs have been increasingly concerned by most people. In the recent literatures, we find that the study of mixed finite element methods for FPDEs is very limited. So far, only a paper [7] has been studied and analyzed for a mixed finite element method of time-FPDE with second-order space derivative. However the theoretical analysis of the (mixed) finite element methods for solving the fractional fourth-order PDEs have not been mentioned and reported. In this article, our goal is to give some detailed numerical analysis of a mixed finite element method for studying the following time-FPDE with fourth-order derivative term
@ a uðx; tÞ Du þ D2 u ¼ f ðx; tÞ; ðx; tÞ 2 X J; @t a
ð1:1Þ
q Foundation item: Supported by the National Natural Science Fund (11301258, 11361035), the Natural Science Fund of Inner Mongolia Autonomous Region (2012MS0108, 2012MS0106), the Scientific Research Projection of Higher Schools of Inner Mongolia (NJZZ12011, NJZY13199, and NJZY14013). ⇑ Corresponding authors. E-mail addresses: [email protected] (Y. Liu), [email protected] (H. Li).
http://dx.doi.org/10.1016/j.amc.2014.06.023 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.
704
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
with boundary condition
uðx; tÞ ¼ Duðx; tÞ ¼ 0; ðx; tÞ 2 @ X J;
ð1:2Þ
and initial condition
uðx; 0Þ ¼ u0 ðxÞ; x 2 X;
ð1:3Þ
d
where X R ðd 6 2Þ and J ¼ ð0; T are a bounded convex polygonal domain with Lipschitz continuous boundary @ X and the a uðx;tÞ time interval with 0 < T < 1, respectively. f ðx; tÞ and u0 ðxÞ are given known functions and @ @t is defined by the following a Caputo fractional derivative
@ a uðx; tÞ 1 ¼ @t a Cð1 aÞ
Z 0
t
@uðx; sÞ ds ; 0 < a < 1: @ s ðt sÞa
ð1:4Þ
Here, we approximate the Caputo fractional derivative by a finite difference method, formulate the mixed weak formulation and fully discrete scheme, prove the stability of the fully discrete scheme and derive the theoretical analysis of some a priori error results in detail. The remaining parts of the article is as follow. In Section 2, we introduce a finite difference method for approximating the Caputo time-fractional derivative, and then discuss the detailed proof of the truncation error. In Section 3, we formulate a fully discrete mixed scheme for the fractional fourth-order PDE (1.1) and derive the stable results for two important variables in detail. Moreover, we prove some a priori error estimates in L2 and H1 -norms. In Section 4, we show some numerical results to illustrate the rationality and effectiveness of our method. In Section 5, we make a brief summary about the presented method and the future development. Throughout this paper, we will denote C > 0 as a generic constant free of the space–time step parameters h and d. At the 2 same time, we define the natural inner product in L2 ðXÞ or ðL2 ðXÞÞ by ð; Þ with the corresponding norm k k. The other notations and definitions of Sobolev spaces can be easily followed in Ref. [29]. 2. Approximation of time-fractional derivative For the discretization for time-fractional derivative, let 0 ¼ t 0 < t 1 < t 2 < < t M ¼ T be a given partition of the time interval ½0; T with step length d ¼ T=M and nodes t n ¼ nd, for some positive integer M. For a smooth function / on ½0; T, define /n ¼ /ðtn Þ. Lemma 2.1. Assuming that u 2 C 2 ð½0; TÞ, then the time fractional derivative 0
@ a uðx;tÞ @t a
at t ¼ t nþ1 can be approximated by, for
i u1 u0 @ a uðx; tnþ1 Þ d1a h ¼ ðn þ 1Þ1a ðnÞ1a a @t Cð2 aÞ d n h 1a i kþ1 4uk þ uk1 X d 1a 1a 3u þ ðn k þ 1Þ ðn kÞ þ Enþ1 0 ; Cð2 aÞ k¼1 2d where Enþ1 ¼ E01 þ Enþ1 0 2 ,
Z
t1
"
# 2 t1 þ t 0 @ uðx; t12 Þ ds 2 2 þ Oðð s t ; 1 Þ Þ þ Oðd Þ 2 2 ðt nþ1 sÞa @t 2
E01 ¼
1 Cð1 aÞ
Enþ1 2
" # n Z tkþ1 X 1 @ 2 uðx; tkþ1 Þ ds 2 2 ¼ ðs t kþ1 Þ þ Oððs tkþ1 Þ Þ þ Oðd Þ : Cð1 aÞ k¼1 tk ðt nþ1 sÞa @t2
s
t0
ð2:1Þ
ð2:2Þ
and
ð2:3Þ
Proof. Writing the integral (1.4) into two parts, we get
@ a uðx; tnþ1 Þ 1 ¼ @t a Cð1 aÞ
Z
t1
t0
n Z t kþ1 X @uðx; sÞ ds 1 @uðx; sÞ ds : ¼I þ II: aþ @ s ðtnþ1 sÞ @ s ðtnþ1 sÞa Cð1 aÞ k¼1 tk
ð2:4Þ
Then discretizing parts I and II by the similar approximated schemes to [8,7], respectively, we get the discrete formula (2.1) with truncation errors (2.2) and (2.3). h
705
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
Lemma 2.2. With u 2 C 2 ð½0; TÞ, the truncation error Enþ1 is bounded by 0
! M D d2a C a d2 MD d þ Cd2 1a : þ ln T ¼CðaÞðd2a þ d2 þ ln dÞ; þ Cð3 aÞ Cð1 aÞ Cð2 aÞ
jEnþ1 0 j 6
where M 4 ¼ sup t2½0;T
o n 2 @ uðx;tÞ @t2 ;
ð2:5Þ
l0 ¼ 0; ln ¼ 1; for any n P 1.
Proof. Now we estimate the truncation error E01 . Firstly, we only need to consider the following term. Using the Taylor formula, we have
Z
t1
s
t0
Z t1 Z t1 t1 þ t0 ds t1 þ t0 ds 1a ¼ ðt s Þ d s þ t nþ1 nþ1 a 2 2 ðt nþ1 sÞa t0 t 0 ðt nþ1 sÞ 1 2n þ 1 ¼ d2a ðn2a ðn þ 1Þ2a Þ þ ððn þ 1Þ1a n1a Þ 2a 2ð1 aÞ a 1 2a 2a ¼d ððn þ 1Þ n2a Þ þ ðnðn þ 1Þ1a ðn þ 1Þn1a Þ 2ð1 aÞ 2ð1 aÞð2 aÞ " ! !# 1a 1a a 1 1 1 ¼ d2a n þ n 1þ ðn þ 1Þ n1a ðn þ 1Þ 1 þ n 2ð1 aÞ n 2ð1 aÞð2 aÞ a 1 1 n ¼ d2a ðn þ 1Þ 1 þ ð1 aÞ þ O 2 n n 2ð1 aÞð2 aÞ 1 1 1 nð1 þ ð1 aÞ þ O 2 Þ ðn þ 1Þ n1a þ 2ð1 aÞ n n a 1 1 1 Oð1Þd2a na ¼ d2a 2aþO : þ n1a ¼ a þ O n 2ð1 aÞ n 2ð1 aÞð2 aÞ ð1 aÞð2 aÞ ð2:6Þ
By the simple calculation for (2.2), we have
Z t1 h i @ 2 uðx; t1 Þ Oð1Þd2a na 1 ds 2 O ðs t1 Þ2 þ Oðd2 Þ þ a 2 2 C ð3 a Þ C ð1 a Þ ðt @t t0 nþ1 sÞ h i a Oððth t 1 Þ2 Þ þ Oðd2 Þ na @ 2 uðx; t1 Þ Oð1Þd2a na n 2 2 ¼ n þ ðn þ 1Þ 2 nþ1 Cð3 aÞ Cð2 aÞ @t h i " # 2 2 a 2 Oððth t 1 Þ Þ þ Oðd Þ n @ uðx; t1 Þ Oð1Þd2a na ð1 aÞ 2 2 ; þ ¼ Cð3 aÞ Cð2 aÞ 1 þ a 1n þ Oðn12 Þ @t2
E01 ¼
ð2:7Þ
where t 0 < t h < t1 . Noting that ðnÞa 6 1, we get
jE01 j 6
M D d2a C a d2 ; þ Cð3 aÞ Cð1 aÞ
ð2:8Þ
Secondly, we now start to estimate the error bounds for jEnþ1 2 j. By using the integral calculation, we easily obtain
jEnþ1 2 j 6
n M D d þ Cd2 X Cð1 aÞ k¼1
Z
t kþ1
tk
ds MD d2a þ Cd3a 1a M D d þ Cd2 1a n 6 T : a ¼ Cð2 aÞ Cð2 aÞ ðt nþ1 sÞ
ð2:9Þ
Combining (2.8) and (2.9), we get 0 nþ1 jEnþ1 0 j 6 jE1 j þ ln jE2 j 6
where
MD d2a C a d2 M D d þ Cd2 1a T ¼ CðaÞðd2a þ d2 þ ln dÞ; þ þ ln Cð3 aÞ Cð1 aÞ Cð2 aÞ
l0 ¼ 0; ln ¼ 1; for any n P 1. h
3. Mixed method and the analysis of stability 3.1. Weak formulation and semi-discrete scheme Introducing an auxiliary variable
r ¼ Du, we first rewrite Eq. (1.1) into the following coupled system
ð2:10Þ
706
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
@ a uðx; tÞ Du þ Dr ¼ f ðx; tÞ; @t a
ð3:1Þ
and
r Du ¼ 0:
ð3:2Þ
By using Green’s formula, we get the following mixed weak formulation of (3.1) and (3.2): find fu; rg : ½0; T # such that
H10
H10
a @ uðx; tÞ ; v þ ðru; rv Þ ðrr; rv Þ ¼ ðf ; v Þ; 8v 2 H10 ; a @t
ð3:3Þ
ðr; wÞ þ ðru; rwÞ ¼ 0; 8w 2 H10 :
ð3:4Þ
and
We now define the finite element space V h as
V h ¼ fv h 2 H10 ðXÞ \ C 0 ðXÞjv h je 2 Pm ðeÞ; 8e 2 Kh g;
ð3:5Þ
where Kh is a quasiuniform partition of the domain X. Then the corresponding semi-discrete mixed finite element scheme of (3.3) and (3.4) is to find fuh ; rh g : ½0; T # V h V h such that
a @ uh ðx; tÞ ; v h þ ðruh ; rv h Þ ðrrh ; rv h Þ ¼ ðf ; v h Þ; 8v h 2 V h ; a @t
ð3:6Þ
ðrh ; wh Þ þ ðruh ; rwh Þ ¼ 0; 8wh 2 V h :
ð3:7Þ
and
Remark 3.1 (i) So far, we have not found any related studies on the mixed finite element method for solving fractional fourth-order PDEs. In this article, we derive some detailed theoretical process on a mixed finite element method for solving the fractional fourth-order PDEs (1.1). (ii) We can also get another mixed weak formulation
a @ uðx; tÞ ; v ðr; v Þ ðrr; rv Þ ¼ ðf ; v Þ; 8v 2 H10 ; a @t
ð3:8Þ
ðr; wÞ þ ðru; rwÞ ¼ 0; 8w 2 H10 :
ð3:9Þ
and
(iii) For the system (3.8) and (3.9), we can get some similar results by the similar analysis to that based on the system (3.3) and (3.4), so we omit the corresponding discussions in this article. 3.2. The analysis of stability for discrete scheme In this subsection, we will give the analysis of stability for the case 0 < a < 1. For the convenience of theoretical analysis, we now denote
(
Gak ¼ ðk þ 1Þ
1a
k
1a
and dt ukþ1 ¼
u1 u0 d
;
3ukþ1 4uk þuk1 2d
k¼0 ; k P 1:
Based on the discrete formula (2.1), the system (3.3) and (3.4) can be rewritten as the following equivalent weak formulation
d1a
n X 1 Gank ðdt ukþ1 ; v Þ þ ðrunþ1 ; rv Þ ðrrnþ1 ; rv Þ ¼ ðf nþ1 þ Enþ1 0 ; v Þ; 8v 2 H 0 ;
Cð2 aÞ k¼0
ð3:10Þ
and
ðrnþ1 ; wÞ þ ðrunþ1 ; rwÞ ¼ 0; 8w 2 H10 : Now, a fully discrete finite element procedure is to determine 1a
d
n X
Cð2 aÞ k¼0
ð3:11Þ ðunþ1 h ;
r
nþ1 h Þ
2 V h V h ; ðn ¼ 0; 1; . . . ; M 1Þ such that
nþ1 nþ1 nþ1 Gank ðdt ukþ1 ; v h Þ; 8v h 2 V h ; h ; v h Þ þ ðruh ; rv h Þ ðrrh ; rv h Þ ¼ ðf
ð3:12Þ
707
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
and nþ1 ðrnþ1 h ; wh Þ þ ðruh ; rwh Þ ¼ 0; 8wh 2 V h :
ð3:13Þ
In the next process, we firstly consider the stability for scheme (3.12) and (3.13). Theorem 3.2. For the scheme (3.12) and (3.13), the following stable inequality for unh 2 V h holds, for sufficiently small d
kunh k 6 Cðku0h k þ ka ðdÞmax kf k kÞ;
ð3:14Þ
06k6n
: where ka ðdÞ¼2da Cð2 aÞ; f 2 L2 ðXÞ and C > 0 is a constant free of two mesh parameters d and h. Proof. We take
v h ¼ unþ1 ; h
wh ¼ rnþ1 in (3.12) and (3.13), respectively, and add the two equations to get h
n d1a X nþ1 2 nþ1 2 nþ1 Ga ðd ukþ1 ; unþ1 ; unþ1 h Þ þ kruh k þ krh k ¼ ðf h Þ: Cð2 aÞ k¼0 nk t h
ð3:15Þ
Now we handle the first term on the left hand side of equality (3.15). We easily get the following equality by the simple verification n n1 d1a X d1a d1a X Gank ðdt uhkþ1 ; unþ1 Gan ðdt u1h ; unþ1 Ga ðd unkþ1 ; unþ1 h Þ ¼ h Þþ h Þ Cð2 aÞ k¼0 Cð2 aÞ Cð2 aÞ k¼0 k t h
¼
da
Cð2 aÞ
Gan ðu1h u0h ; unþ1 h Þþ
n1 X da Ga ð3uhnkþ1 4unk þ unk1 ; unþ1 h h h Þ: 2Cð2 aÞ k¼0 k
Noting that Ga0 ¼ 1, it follows that n1 X
Gak 3unkþ1 4unk þ uhnk1 ; unþ1 ¼ h h h
n2 X
3
k¼0
a
Gkþ1 unk h
k¼1
¼
2 3kunþ1 h k
n1 n X X nþ1 4 Gak unk þ Gak1 unk h h ; uh k¼0
a n
a
þ 3G1 4G0
uh ; unþ1 h
ð3:16Þ
!
k¼1
n1 X nþ1 þ ð3Gakþ1 4Gak þ Gak1 Þunk h ; uh
!
3Gan u1h ; unþ1 h
k¼1
þ Gn1 u0h ; unþ1 : h a
ð3:17Þ
: Apply a combination for (3.17), (3.16) and (3.15), multiply (3.15) by ka ðdÞ¼2da Cð2 aÞ and use Cauchy–Schwarz inequality to arrive at 2 nþ1 2 nþ1 2 nþ1 3kunþ1 ;unþ1 h k þ ka ðdÞðkruh k þ krh k Þ ¼ ka ðdÞðf h Þ
n1 X nþ1 ð3Gakþ1 4Gak þ Gak1 Þunk h ;uh
!
k¼1
a 1 nþ1 a a 0 nþ1 þ ð4Ga0 3Ga1 Þðunh ; unþ1 h Þ þ Gn ðuh ;uh Þ þ ð2Gn Gn1 Þðuh ;uh Þ
6
n1 X a a a a a 1 n 0 j3Gakþ1 4Gak þ Gak1 jkunk h k þ j4G0 3G1 jkuh k þ j2Gn Gn1 jkuh k þ Gn kuh k k¼1
þka ðdÞkf nþ1 k kunþ1 h k: 2 unþ1 h k
Remove the second positive term ka ðdÞðkr
nþ1 u 6 h
nþ1 2 k Þ, h
þ kr
ð3:18Þ
then divide kunþ1 h k in inequality (3.18) to obtain
! n1 X 3Ga 4Ga þ Ga kunk k þ j4Ga 3Ga jkun k þ j2Ga Ga jku0 k þ Ga ku1 k þ ka ðdÞkf nþ1 k : h h h h kþ1 k k1 0 1 n n1 n
ð3:19Þ
k¼1
We apply triangle inequality j3Gakþ1 4Gak þ Gak1 j 6 3jGakþ1 Gak j þ jGak Gak1 j; j2Gan Gan1 j 6 jGan1 Gan j þ Gan and note that inequality 1 > Gak1 > Gak > Gakþ1 > 0 in (3.19) to get
kunþ1 h k 6
! n1 n1 X X 1 a a a a a a 1 nk n 0 nþ1 3 ðGak Gakþ1 Þkunk k þ ðG G Þku k þ ð4G 3G Þku k þ G ku k þ G ku k þ k ðdÞkf k : a h h h h h k1 k 0 1 n1 n 3 k¼1 k¼1
ð3:20Þ ku1h k
ku0h k.
Now we derive the relationship between and When n ¼ 0 in (3.12) and (3.13), we take (3.12) and (3.13), respectively, and add the two equations to get
da ðu1 u0h ; u1h Þ þ kr1h k2 þ kru1h k2 ¼ ðf 1 ; u1h Þ: Cð2 aÞ h
vh ¼
u1h
and wh ¼ r1h in
ð3:21Þ
708
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
Multiply (3.21) by da Cð2 aÞ and use Cauchy–Schwarz inequality to get
ku1h k2 þ da Cð2 aÞ kr1h k2 þ kru1h k2 6 ðku0h k þ da Cð2 aÞkf 1 kÞku1h k:
ð3:22Þ
By the simple calculation, we easily get
ku1h k 6 ku0h k þ da Cð2 aÞkf 1 k:
ð3:23Þ
Make a combination for (3.23) and (3.20) to get
! n1 n1 X X 1 a a a a a a a a nk nk n 0 nþ1 1 3 ðGk Gkþ1 Þkuh k þ ðGk1 Gk Þkuh k þ ð4G0 3G1 Þkuh k þ ðGn1 þ Gn Þkuh k þ ka ðdÞkf k þ ka ðdÞkf k : 3 2 k¼1 k¼1
1 kunþ1 h k6
ð3:24Þ From (3.24), we can prove that 0 k kunþ1 h k 6 Cðkuh k þ ka ðdÞ max kf kÞ:
ð3:25Þ
06k6nþ1
In the following discussion, we use mathematical induction to prove the inequality (3.25). When n ¼ 1 in (3.12) and (3.13), we take v h ¼ u2h in (3.12) and (3.13) and sum the two equations to derive
da da Ga1 ðu1h u0h ; u2h Þ þ ð3u2h 4u1h þ u0h ; u2h Þ þ kr2h k2 þ kru2h k2 ¼ ðf 2 ; u2h Þ: Cð2 aÞ 2Cð2 aÞ
ð3:26Þ
Multiply (3.26) by ka ðdÞ and use Cauchy–Schwarz inequality to get
3ku2h k2 6 ka ðdÞðf 2 ; u2h Þ þ 2ðu1h ; u2h Þ þ ðu0h ; u2h Þ 6 ðka ðdÞkf 2 k þ 2ku1h k þ ku0h kÞku2h k: Divide
ku2h k
ð3:27Þ
in (3.27) and substitute (3.23) into (3.27) to get
ku2h k 6
1 1 ðka ðdÞkf 2 k þ 2ku1h k þ ku0h kÞ 6 3ku0h k þ ka ðdÞðkf 2 k þ kf 1 kÞ ¼ C ku0h k þ ka ðdÞmax kf k k : 06k61 3 3
ð3:28Þ
It is not hard to see from (3.28) that (3.25) holds for the case n ¼ 1. We now assume that
kunh k 6 C ku0h k þ ka ðdÞmax kf k k ;
ð3:29Þ
06k6n
holds for the cases n ¼ 1; 2; . . . ; j, then we prove that (3.25) holds based on the induction hypothesis (3.29). We now consider (3.24) with the simple calculation to get for the case n ¼ j þ 1
kujþ1 h k 6
1 ð½3ðGa1 Gaj Þ þ ðGa0 Gaj1 ÞC ku0h k þ ka ðdÞ max kf k k þ ð4Ga0 3Ga1 ÞC ku0h k þ ka ðdÞmaxkf k k 06k6j1 06k6j 3 1 a a þ Gj1 þ Gj ku0h k þ ka ðdÞkf jþ1 k þ ka ðdÞkf 1 kÞ: 2
ð3:30Þ
By the simplification for (3.30), we have for n ¼ j þ 1
0 k kujþ1 h k 6 C kuh k þ ka ðdÞ max kf k :
ð3:31Þ
06k6jþ1
Based on (3.28) and (3.31), we easily see that the conclusion (3.25) of Theorem 3.2 holds. In the next theorem, we will mainly discuss the stability of variable rnh . h Theorem 3.3. For the scheme (3.12) and (3.13), the following stable inequality for
pffiffiffiffiffiffiffiffiffiffiffi krnh k þ krunh k 6 C kr0h k þ kru0h k þ ka ðdÞ max kf k k ;
rnh 2 V h holds, for sufficiently small d
06k6n
ð3:32Þ
where ka ðdÞ is defined as Theorem 3.2 and C is a positive constant free of space–time step h and d. Proof. In order to analyze the stable result for following
rnh , we have to handle the scheme (3.13). Now we rewrite (3.13) as the
n n d1a X d1a X Gank ðdt rkþ1 Ga ðd rukþ1 ; rwh Þ ¼ 0: h ; wh Þ þ Cð2 aÞ k¼0 Cð2 aÞ k¼0 nk t h
Set wh ¼ rnþ1 in (3.33) and h
v h ¼ Cdð2aaÞ 1
Pn
a kþ1 k¼0 Gnkdt uh
in (3.12) and take the sum of the two equations to arrive at
ð3:33Þ
709
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
2 d1a X n n
d1a X a kþ1 kþ1 nþ1 G du þ Ga ðd rkþ1 ; rnþ1 h Þ þ ðdt ruh ; ruh Þ Cð2 aÞ k¼0 nk t h Cð2 aÞ k¼0 nk t h ¼
n X
d1a
Cð2 aÞ k¼0
nþ1 Gank ðdt ukþ1 Þ; h ;f
ð3:34Þ
Use Cauchy–Schwarz inequality and Young inequality to get n
1 nþ1 2 d1a X kþ1 nþ1 Ga ðd rkþ1 ; rnþ1 kf k : h Þ þ ðdt ruh ; ruh Þ 6 2 Cð2 aÞ k¼0 nk t h
ð3:35Þ
Using the similar process of calculation to (3.20), we have 2 jk!nþ1 h kj 6
! n1 n1 X 2 X nk n 0 1 a a a a a a jk!nþ1 3 ðGak Gakþ1 Þjk!nk kj þ ðG G Þjk ! kj þ ð4G 3G Þjk ! kj þ G jk ! kj þ G jk ! kj k1 k 0 1 n1 n h h h h h h kj 3 k¼1 k¼1 1 þ ka ðdÞkf nþ1 k2 ; 3 ð3:36Þ
: j h kj¼k
j hk
ujh k.
where jk! r þ kr Applying the inequalities we have
2 jk!nþ1 h kj 6 C 3
1 2
P
k j¼1 aj
2
6
Pk
2 j¼1 aj
6
P k
j¼1 aj
2
;2
Pk
j;m¼1;j–m aj am
6
P
k j¼1 aj
2
; aj P 0; ðj; m ¼ 0; 1; 2; . . . ; k),
n1 n1 X X pffiffiffiffiffiffiffiffiffiffiffi nþ1 n 0 1 a a a a ka ðdÞkf k ðGak Gakþ1 Þjk!nk ðGak1 Gak Þjk!nk h kj þ h kj þ ð4G0 3G1 Þjk!h kj þ Gn1 jk!h kj þ Gn jk!h kj þ k¼1
!2
k¼1
1 nþ1 2 þ jk!h kj 2 ð3:37Þ
By the simple calculation for (3.37), we easily get jk!
nþ1 h kj 6 C
! n1 n1 X X pffiffiffiffiffiffiffiffiffiffiffi nþ1 nk nk n 0 1 a a a a a a a a 3 ðGk Gkþ1 Þjk!h kj þ ðGk1 Gk Þjk!h kj þ ð4G0 3G1 Þjk!h kj þ Gn1 jk!h kj þ Gn jk!h kj þ ka ðdÞkf k : k¼1
k¼1
ð3:38Þ
In the next consideration, we use a similar process of proof to (3.25) based on mathematical induction to get 0 jk!nþ1 h kj 6 Cðjk!h kj þ
pffiffiffiffiffiffiffiffiffiffiffi ka ðdÞ max kf k kÞ:
ð3:39Þ
06k6n
Noting that the notation for jk!jh kj, we get the stability on
rnh . h
3.3. Fully discrete error analysis In this subsection, we will give some fully discrete analysis of a priori error estimates for the case 0 < a < 1. Now a Ritzprojection associated with our equations is introduced to help us to analyze the convergence of the method. Lemma 3.4. Define a Ritz-projection operator {h : H10 ðXÞ ! V h by
ðrðu {h uÞ; rwh Þ ¼ 0; 8wh 2 V h ;
ð3:40Þ
and holds
ku {h uk þ hku {h uk1 6 Ch
mþ1
kukmþ1 ; 8u 2 H10 ðXÞ \ Hmþ1 ðXÞ;
P
R
1 2
ð3:41Þ
where the norms is defined by kzkk ¼ ½ 06jbj6k X jDb zj2 dx with the degree k of the polynomial. The contents of Lemma 3.4 can be found in many literatures, such as Refs. [29,28]. In the next analysis, we will discuss the proof of fully discrete a priori error estimates based on the Ritz-projection Lemma 3.4 in detail. To simplify the process of writing in the proof, we now rewrite the errors as
uðtn Þ unh ¼ ðuðt n Þ {h un Þ þ ð{h un unh Þ ¼ fn þ /n ;
rðtn Þ rnh ¼ ðrðtn Þ {h rn Þ þ ð{h rn rnh Þ ¼ pn þ vn :
710
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
Theorem 3.5. Supposing that un 2 H10 ðXÞ \ Hmþ1 ðXÞ; unh 2 V h and u0h ¼ {h uð0Þ, then the error estimate for u in L2 -norm holds with a parameter 0 < a < 1
kuðt n Þ unh k 6
Cðu; aÞT a mþ1 mþ1 ðln d þ d2 þ d2a þ ð1 þ da Þh Þ þ CðuÞh ; 1a
ð3:42Þ
and when a ! 1 mþ1
kuðt n Þ unh k 6 CðuÞðln d þ d2 þ d þ ð1 þ d1 Þh
Þ;
ð3:43Þ
where both Cðu; aÞ (with a) and CðuÞ depend on the norm kukmþ1 . Proof. Subtract (3.12) from (3.10) and apply the Ritz-projection (3.40) to obtain the error equations n n d1a X d1a X Gank ðdt /kþ1 ; v h Þ þ ðr/nþ1 ; rv h Þ ðrvnþ1 ; rv h Þ ¼ Ga ðd fkþ1 ; v h Þ þ ðEnþ1 0 ; v h Þ; 8v h 2 V h ; Cð2 aÞ k¼0 Cð2 aÞ k¼0 nk t
ð3:44Þ
ðvnþ1 ; wh Þ þ ðr/nþ1 ; rwh Þ ¼ ðpnþ1 ; wh Þ; 8wh 2 V h :
ð3:45Þ
and
For the first term on the left hand side of (3.44), we take a similar process to (3.16) to get n n1 X d1a X da da Gank ðdt /kþ1 ; v h Þ ¼ Gan ð/1 /0 ; v h Þ þ Ga ð3/nkþ1 4/nk þ /nk1 ; v h Þ: Cð2 aÞ k¼0 Cð2 aÞ 2Cð2 aÞ k¼0 k
ð3:46Þ
We add (3.44) and (3.45), substitute (3.46), take ðv h ; wh Þ ¼ ð/nþ1 ; vnþ1 Þ and multiply by ka ðdÞ ¼ 2da Cð2 aÞ to arrive at n1 X Gak ð3/nkþ1 4/nk þ /nk1 ; /nþ1 Þ þ ka ðdÞðkvnþ1 k2 þ kr/nþ1 k2 Þ k¼0 n X nþ1 ¼ 2d Gank ðdt fkþ1 ; /nþ1 Þ 2Gan ð/1 /0 ; /nþ1 Þ þ ka ðdÞðEnþ1 Þ ka ðpnþ1 ; vnþ1 Þ: 0 ;/
ð3:47Þ
k¼0
Now we deal with the first term on the left hand side of (3.47). Noting that Ga0 ¼ 1, we use the similar process to (3.17) to arrive at n1 X Gak ð3/nkþ1 4/nk þ /nk1 ; /nþ1 ¼ 3k/nþ1 k2 þ ð3Ga1 4Ga0 Þð/n ; /nþ1 Þ k¼0
þ
! n1 X ð3Gakþ1 4Gak þ Gak1 Þ/nk ; /nþ1 3Gan ð/1 ; /nþ1 Þ k¼1
þ Gan1 ð/0 ; /nþ1 Þ:
ð3:48Þ
By substituting (3.48) into (3.47), we obtain n X nþ1 3k/nþ1 k2 þ ka ðdÞðkvnþ1 k2 þ kr/nþ1 k2 Þ ¼ 2d Gank ðdt fkþ1 ; /nþ1 Þ þ ka ðdÞðEnþ1 Þ 0 ;/ k¼0
ð3Ga1 4Ga0 Þð/n ; /nþ1 Þ
n1 X
! ð3Gakþ1 4Gak þ Gak1 Þ/nk ; /nþ1
k¼1
þ Gan ð/1 ; /nþ1 Þ þ ð2Gan Gan1 Þð/0 ; /nþ1 Þ ka ðpnþ1 ; vnþ1 Þ:
ð3:49Þ
For the next proof, we have to deal with the first term on the right hand side of (3.49). Now we use the similar proof to (3.46) and (3.48) to obtain n n1 X X 2d Gank ðdt fkþ1 ; /nþ1 Þ ¼ 2Gan ðf1 f0 ; /nþ1 Þ Gak ð3fnkþ1 4fnk þ fnk1 ; /nþ1 Þ k¼0
k¼0 n1 X ð3Gakþ1 4Gak þ Gak1 Þfnk ; /nþ1 ¼ 3ðfnþ1 ; /nþ1 Þ
! ð3Ga1 4Ga0 Þðfn ; /nþ1 Þ
k¼1
þ Gan ðf1 ; /nþ1 Þ þ ð2Gan Gan1 Þðf0 ; /nþ1 Þ:
ð3:50Þ
Applying (3.41), Cauchy–Schwarz inequality and the inequality
j3Gakþ1 4Gak þ Gak1 j 6 3jGakþ1 Gak j þ jGak Gak1 j ¼ 3ðGak Gakþ1 Þ þ ðGak1 Gak Þ;
ð3:51Þ
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
711
we have n n1 X X 2d Gank ðdt fkþ1 ; /nþ1 Þ ¼ ð3kfnþ1 k þ j3Gakþ1 4Gak þ Gak1 jkfnk k þ ð4Ga0 3Ga1 Þkfn k þ Gan kf1 k k¼0
k¼1
þ Gan1 kf0 kÞk/nþ1 k 6 CðuÞð3 þ 3ðGa1 Gan Þ þ ðGa0 Gan1 Þ þ ð4Ga0 3Ga1 Þ þ Gan þ Gan1 Þh
mþ1
k/nþ1 k ¼ Cðu; aÞh
mþ1
k/nþ1 k:
ð3:52Þ
Substitute (3.52) into (3.49) and use Cauchy–Schwarz inequality to get
3k/nþ1 k2 þ ka ðdÞðkvnþ1 k2 þ kr/nþ1 k2 Þ 6 ðj4Ga0 3Ga1 jk/n k n1 X
þ
j3Gakþ1 4Gak þ Gak1 jk/nk k þ Gan k/1 k þ ðjGan1 Gan j þ Gan Þk/0 k
k¼1
!
þ ka ðdÞkEnþ1 0 k þ Cðu; aÞh
mþ1
k/nþ1 k þ ka kpnþ1 kkvnþ1 k:
ð3:53Þ
Noting that the following inequalities
1 2 2 2 ða þ bÞ 6 ka þ b ; 2
ac þ bd 6 aðc þ dÞ þ bðc þ dÞ ¼ ða þ bÞðc þ dÞ;
ð3:54Þ
for any a; b; c; d > 0; k 2 N þ , we have
pffiffiffiffiffiffiffiffiffiffiffi 2 1 ðk/nþ1 k þ ka ðdÞkvnþ1 kÞ 6 ðj4Ga0 3Ga1 jk/n k 2 n1 X j3Gakþ1 4Gak þ Gak1 jk/nk k þ Gan k/1 k þ ðjGan1 Gan j þ Gan Þk/0 k þ k¼1
þ
ka ðdÞkEnþ1 0 k
þ Cðu; aÞh
mþ1
! pffiffiffiffiffi pffiffiffiffiffi nþ1 þ ka kp k ðk/nþ1 k þ ka kvnþ1 kÞ:
ð3:55Þ
Use inequality (3.51) and the inequality 1 > Gak1 > Gak > Gakþ1 > 0 to get
k/nþ1 k 6 2ðð4Ga0 3Ga1 Þk/n k þ ½3ðGan1 Gan Þ þ Ga0 Ga1 þ 4
n1 X ðGak1 Gak Þk/nk k þ Gan k/1 k þ Gan1 k/0 k k¼2
þ ka ðdÞkEnþ1 0 k þ Cðu; aÞh
mþ1
Þ:
ð3:56Þ
Based on (3.56), we can derive
k/nþ1 k 6
Cðu; aÞ mþ1 ðln d1þa þ d2þa þ d2 þ h Þ: Gan
ð3:57Þ
In the following analysis and discussion, we use mathematical induction to discuss the conclusion (3.57). When n ¼ 0 in (3.44) and (3.45), taking v h ¼ /1 and wh ¼ v1 in (3.44) and (3.45), respectively, and summing the two equations, we use Cauchy–Schwarz inequality to obtain by the simplification
k/1 k2 k/0 k2 1 1 þ k/ /0 k2 þ da Cð2 aÞðkv1 k2 þ kr/1 k2 Þ 2 2 ¼ ðf1 f0 ; /1 Þ þ da Cð2 aÞðE10 ; /1 Þ da Cð2 aÞðp1 ; v1 Þ 6 ðkf1 k þ kf0 k þ CðaÞðd2þa þ d2 ÞCð2 aÞÞk/1 k þ da Cð2 aÞkp1 kkv1 k:
ð3:58Þ
Noting that /0 ¼ 0 and using two important inequalities (3.54), we use the similar discussion to (3.55) to obtain
k/1 k 6 Cðu; aÞðh
mþ1
þ d2þa þ d2 Þ:
ð3:59Þ
Based on the case n ¼ 0, we now discuss the case n ¼ 1, which is the first step of mathematical induction. Under the condition of n ¼ 1, setting v h ¼ /2 in (3.44) and wh ¼ v2 in (3.45), and taking the sum of the two equations, we arrive at
1 ðk/2 k2 þ k2/2 /1 k2 þ k/2 2/1 þ /0 k2 Þ þ ka ðdÞðkv2 k2 þ kr/2 k2 Þ 2 1 ¼ ðk/1 k2 þ k2/1 /0 k2 Þ Ga1 ð/1 /0 ; /2 Þ Ga1 ðf1 f0 ; /2 Þ þ ka ðdÞðE20 ; /2 Þ ka ðdÞðp2 ; v2 Þ: 2 By an application of Cauchy–Schwarz inequality and Young inequality, we can get
ð3:60Þ
712
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
k/2 k2 þ k2/2 /1 k2 þ k/2 2/1 þ /0 k2 6 CðaÞðk/1 k2 þ k/0 k2 þ kf1 k2 þ kf0 k2 þ ka ðdÞkp2 k2 þ ln d2þ2a þ d4þ2a þ d4 Þ:
ð3:61Þ
Combining (3.59) and (3.41), we arrive at
k/2 k 6 Cðu; aÞðh ¼
mþ1
þ ln d1þa þ d2þa þ d2 Þ ¼ ð21a 1Þ
Cðu; aÞ mþ1 ðh þ ln d1þa þ d2þa þ d2 Þ Ga1
Cðu; aÞ mþ1 ðh þ ln d1þa þ d2þa þ d2 Þ; Ga1
ð3:62Þ
which is to tell that (3.57) holds for the first case n ¼ 1. We firstly assume that (3.57) holds, for 1 6 n 6 j
k/jþ1 k 6
Cðu; aÞ mþ1 ðln d1þa þ d2þa þ d2 þ h Þ: Gaj
ð3:63Þ
Now, we consider the case n ¼ j þ 1. By (3.56) and the induction hypothesis (3.63), we arrive at n1 X 1 mþ1 ðjþ1Þþ1 k/ðjþ1Þþ1 k 6 ðð4Ga0 3Ga1 Þk/jþ1 k þ ½3ðGan1 Gan Þ þ Ga0 Ga1 þ 4 ðGak1 Gak Þk/ðjþ1Þk k þ Gan k/1 k þ Gan1 k/0 k þ ka ðdÞkE0 k þ Cðu; aÞh Þ 3 k¼2
6 Cðu; aÞðð4Ga0 3Ga1 Þ
Noting that
1 Gajk
ðjþ1Þþ1
k/
1 1 mþ1 þ ½3ðGan1 Gan Þ þ Ga0 Ga1 þ 4ðGa1 Gan1 Þ a þ Gan þ 2Þðln d1þa þ d2þa þ d2 þ h Þ: Gaj Gjk
ð3:64Þ
< G1a < Ga1 and 0 < Gajþ1 < Gaj < 1, we have j
jþ1
mþ1
k 6 Cðu; aÞð4Ga0 3Ga1 þ 3ðGan1 Gan Þ þ Ga0 Ga1 þ 4ðGa1 Gan1 Þ þ Gan þ 2Þðln d1þa þ d2þa þ d2 þ h 6
Cðu; aÞ mþ1 ðln d1þa þ d2þa þ d2 þ h Þ: Gajþ1
Þ
1 Gajþ1 ð3:65Þ
that is to say that (3.57) holds for n ¼ j þ 1. By using induction based on (3.62) and (3.65), we announce that (3.57) holds for any n P 1. a In order to get the final result based on (3.65), we have to give the estimate for ðnþ1Þ . We now consider Taylor formula to Gan get
ðn þ 1Þa ðn þ 1Þa 1 1 ¼ ¼ a a ¼ a2 2 Gn ðn þ 1Þ1a n1a ðn þ 1Þ nð1 þ 1nÞ ðn þ 1Þ nð1 þ a 1n þ aða2!1Þ ð1 þ h 1nÞ ð1nÞ Þ ¼
1 1aþ
að1aÞ 2!
ð1 þ
a2 h 1nÞ 1nÞ
<
1 ; 1a
ð3:66Þ
where 0 < h < 1 and 0 < a < 1. Noting that the result of (3.66) and ðn þ 1Þd 6 Md ¼ T, we have
Cðu; aÞ Cðu; aÞðn þ 1Þa Cðu; aÞT a da ¼ ðn þ 1Þa da da 6 : a a 1a Gn Gn
ð3:67Þ
Combining (3.57) and (3.67), we get
k/nþ1 k 6
Cðu; aÞT a mþ1 ðln d þ d2 þ d2a þ da h Þ: 1a
ð3:68Þ
By a combination for (3.41) and (3.68) with triangle inequality, we get the error result (3.42) of theorem. Now we briefly discuss the case a ! 1. We note that the error result (3.42) in this case has no meaning since the aÞT a coefficient Cðu; 1a ! 1 as a ! 1. So, we have to seek another error estimate’s form. Considering the fact nd 6 T ðn ¼ 1; 2; . . . ; MÞ, we can get the error result
k/n k 6 CðuÞndðln d þ d2 þ d þ d1 h
mþ1
Þ:
ð3:69Þ
Now we can use induction to derive the proof. The detail process is similar to the above proof, so we do not derive that again. Combining (3.69) and (3.41) with triangle inequality, we can arrive at the estimate (3.43). In the next analysis, our aim is to focus on deriving a priori error estimates in L2 -norm for variable r. At the same time, we also analyze the error result in H1 -norm for u. h Theorem 3.6. With un ; rn 2 H10 ðXÞ \ Hmþ1 ðXÞ; unh ; estimates with 0 < a < 1
rnh 2 V h ; u0h ¼ {h uð0Þ and r0h ¼ {h rð0Þ, one has the following a priori error
713
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
krðt n Þ rnh k 6 kuðt n Þ unh k1 6
3a mþ1 Cðu; r; aÞT a mþ1 mþ1 ðln d þ d2 þ d2a þ d 2 h þ h Þ þ CðrÞh ; 1a
ð3:70Þ
3a mþ1 Cðu; r; aÞT a mþ1 m ðln d þ d2 þ d2a þ d 2 h þ h Þ þ CðuÞh ; 1a
ð3:71Þ
and when a ! 1 3
krðt n Þ rnh k 6 Cðu; rÞTðln d þ d2 þ d þ d2 h 3
mþ1
kuðt n Þ unh k1 6 Cðu; rÞTðln d þ d2 þ d þ d2 h
mþ1
þh
mþ1
þh
mþ1
Þ þ CðrÞh
mþ1
ð3:72Þ
;
m
Þ þ CðuÞh :
ð3:73Þ
where Cðu; r; aÞ (with a) and Cðu; rÞ are dependent on the norms kukmþ1 and krkmþ1 ; CðrÞ only depends on krkmþ1 and CðuÞ is described as the statement in Theorem 3.5. Proof. For the need of proof, we have to rewrite (3.13) as n n n d1a X d1a X d1a X Gank ðdt vkþ1 ; wh Þ þ Gank ðrdt /kþ1 ; rwh Þ ¼ Ga ðd pkþ1 ; wh Þ: Cð2 aÞ k¼0 Cð2 aÞ k¼0 Cð2 aÞ k¼0 nk t 1a
We add (3.74) and (3.12), take wh ¼ vnþ1 and v h ¼ Cdð2aÞ ity to get
Pn
a kþ1 , k¼0 Gnkdt /
ð3:74Þ
multiply by ka ðdÞ and use Cauchy–Schwarz inequal-
2 d1a X n n X ka ðdÞ Gankdt /kþ1 þ 2d Gank ½ðdt vkþ1 ; vnþ1 Þ þ ðrdt /kþ1 ; r/kþ1 Þ Cð2 aÞ k¼0 k¼0 ! n n n X a X d1a X a kþ1 nþ1 kþ1 2d Gank ðdt pkþ1 ; vnþ1 Þ Gnkdt / ¼ 2d Gnkdt f þ ka ðdÞE0 ; Cð2 aÞ k¼0 k¼0 k¼0 ! X d1a X n n n X Ga d /kþ1 þ k2d Gankdt pkþ1 kkvnþ1 k: 6 2d Gankdt fkþ1 þ ka ðdÞkEnþ1 0 k k¼0 Cð2 aÞ k¼0 nk t k¼0
ð3:75Þ
Using Young inequality and the first inequality in (3.54), we have
2 d1a X n n X 1 ka ðdÞ Gankdt /kþ1 þ 2d Gank ½ðdt vkþ1 ; vnþ1 Þ þ ðrdt /kþ1 ; r/kþ1 Þ Cð2 aÞ k¼0 2 k¼0 2 2 X n n 2 1 X a kþ1 2 d k þ 2d G p 6 ðka ðdÞÞ1 2d Gankdt fkþ1 þ ka ðdÞkEnþ1 þ e vnþ1 : nk t 0 k¼0 4e k¼0
ð3:76Þ
2
u
1 0 1
−1 −2 1
0.8 0.6 0.8
0.6
0.4 0.4
0.2
0.2 0
0
Fig. 1. Surface for exact solution u.
714
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
Table 1 The errors and convergence orders for u and
a
M
pffiffiffi
r with P1 element and space–time mesh d ¼ h= 2 ¼ 1=M. Order
L2 -norm u
Order
H1 -norm u
L2 -norm
Order
r
0.3
16 32 64
6.7904e002 1.7450e002 4.3932e003
1.9603 1.9899
1.6779e+000 8.3477e001 4.1677e001
1.0072 1.0021
2.9915e+000 7.5426e001 1.8901e001
1.9877 1.9966
0.5
16 32 64
6.7903e002 1.7450e002 4.3932e003
1.9603 1.9899
1.6779e+000 8.3477e001 4.1677e001
1.0072 1.0021
2.9915e+000 7.5424e001 1.8900e001
1.9878 1.9966
0.7
16 32 64
6.7902e002 1.7450e002 4.3930e003
1.9603 1.9898
1.6779e+000 8.3477e001 4.1677e001
1.0072 1.0021
2.9914e+000 7.5421e001 1.8899e001
1.9878 1.9967
0.9
16 32 64
6.7900e002 1.7448e002 4.3924e003
1.9603 1.9900
1.6779e+000 8.3476e001 4.1677e001
1.0072 1.0021
2.9912e+000 7.5413e001 1.8895e001
1.9879 1.9968
2
u
h
1 0 1
−1 −2 1
0.8 0.6 0.8
0.6
0.4 0.4
0.2
0.2 0
0
Fig. 2. Surface for numerical solution uh .
Now we discuss the results for any n P 0. Firstly, taking n ¼ 0 and
e ¼ 12 in (3.76), we have
2 da 1 1 a 1 a a 1 1 2 0 1 0 1 0 2 1 0 1 ka ðdÞ G ð/ / Þ Cð2 aÞ 0 þ 2G0 ½ðv v ; v Þ þ ðr/ r/ ; r/ Þ 6 ðka ðdÞÞ k2G0 ðf f Þk þ ka ðdÞkE0 k 2 1 1 þ k2Ga0 ðp1 p0 Þk2 þ kv1 k2 : 2 2 Noting that Ga0 ¼ 1;
ð3:77Þ
v0 ¼ 0; /0 ¼ 0 and (3.41), we have by the simple simplification
2 da 1 3 1 0 1 2 þ 2 r/1 2 6 CðuÞðka ðdÞÞ1 h2þ2m þ CðaÞka ðdÞðd42a þ d4 Þ þ CðrÞh2þ2m : ka ðdÞ / Þ ð/ Cð2 aÞ þ2 v 2
ð3:78Þ
From (3.78), we easily find that
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi da Cð2 aÞ
1 da 1þm 1 a2 1þm þ d2 þ d2þa þ h Þ: ð/1 /0 Þ þ v þ r/ 6 Cðu; r; aÞðd h Cð2 aÞ
For the n P 1 in (3.76), we use mathematical induction and a similar discussion to Theorem 3.5 to get
ð3:79Þ
715
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
200
σ
100 0 1
−100 −200 1
0.8 0.6 0.8
0.4
0.6
0.4
0.2
0.2 0
0
Fig. 3. Surface for exact solution
r.
200
σh
100 0 1
−100 −200 1
0.8 0.6 0.8
0.6
0.4 0.4
0.2
0.2 0
Fig. 4. Surface for numerical solution
0 rh .
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 1a X 3a 1þm Cðu; r; aÞT a d 1þm a kþ1 a d Cð2 aÞ Gnk dt / þ kvn k þ kr/n k 6 ðln d þ d2 þ d2a þ d 2 h þ h Þ: Cð2 aÞ k¼0 1a
ð3:80Þ
Combining (3.80) and (3.41) with triangle inequality, we get the estimates (3.70) and (3.71) for the case 0 < a < 1. For the case a ! 1, use a similar analysis as that in Theorem 3.5 to get the error estimates (3.72) and (3.73). h Remark 3.7. From the error results in Theorems 3.5 and 3.6, we can see that the results of convergence in the spatial direction is affected by the inverse of the time mesh d. However the time step’s inverse would not affect the global accuracy, mþ1 because h can be much smaller than d for sufficiently regular solutions on the spatial variable. The similar discussions can be found in the Remark 4.2 of Lin and Xu Ref. [8].
716
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
4. Some numerical results In this section, we provide a numerical test to confirm the theoretical analysis of our method. Here, we choose h i 2 2t 2a 2 4 Cð3aÞ þ ð8p þ 64p Þðt þ 1Þ sinð2pxÞ sinð2pyÞ; 8ðx; yÞ 2 ½0; 1 ½0; 1; t 2 ½0; 1 and initial value function
function f ðx; tÞ ¼
u0 ðx; y; 0Þ ¼ sinð2pxÞ sinð2pyÞ. We now get the exact solution
uðx; y; tÞ ¼ ðt2 þ 1Þ sinð2pxÞ sinð2pyÞ;
ð4:1Þ
r ¼ Du ¼ 8p2 ðt2 þ 1Þ sinð2pxÞ sinð2pyÞ:
ð4:2Þ
and
We now take the finite element space V h as the linear space with continuous piecewise linear functions, and divide the qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 1 2ffi pffiffiffi 1 2 2 spatial domain ½0; 1 ½0; 1 into uniform triangulations with mesh size h ¼ hx þ hy ¼ þ N ¼ 2 N. For formulating N fully discrete procedure, we now consider the discrete time step d ¼ M1 for the time interval ½0; 1. For the convenience of calculation, we now take M ¼ N ¼ 16; 32; 64, respectively, then compute some a priori errors and orders of convergence of u in both L2 and H1 -norms for different a ¼ 0:3; 0:5; 0:7; 0:9. From the results of convergence in Table 1, we can see that the rates of convergence of u in L2 and H1 -norms are close to 2 and 1, respectively, which are in mþ1 accord with the spatial convergence order Oðh Þ with index m ¼ 1. Moreover, the error in L2 -norm for the variable r is also computed and the rate of convergence is close to 2. In Fig. 1, the surface of the exact solution u is shown. With the chosen pffiffiffi pffiffiffi space–time steps h ¼ 2d ¼ 2=16, the surface of the numerical solution uh is displayed at time t ¼ 1 in Fig. 2. It not hard to see that the numerical uh can approximate well the exact solution u. The similar comparison between rh and r is also made in Figs. 3 and 4. Based on the results of errors and rates of convergence, we know our method is feasible and effective. 5. Some remarks and extensions So far, we have not found any reports on the (mixed) finite element method for solving the fractional fourth-order PDE (1.1). In this article, our aim is propose a mixed finite element method for the fractional fourth-order PDE with the case 0 < a < 1. Based on a finite difference method for time fractional derivative, we give a fully discrete mixed finite element scheme, and then analyze the detailed proofs’ process for both the stabilities in L2 -norm and some a priori error results in L2 and H1 -norms. Finally, we present a numerical example to show the feasibility of the studied method. For the case 1 < a < 2, we will give the numerical analysis and simulation based on the mixed finite element method in another work. In the near future, we will discuss a moving (mixed) finite element schemes based on the discussion in [19,21] and obtain the optimal convergence analysis of mixed finite element method for fractional fourth-order PDE by using a special interpolation operator in Ref. [25]. At the same time, we will apply some other (mixed) finite element methods [26,27] to analyze the fractional fourth-order PDE. Acknowledgments Authors thank the reviewers and editor for their valuable comments and suggestions. This work is supported by the National Natural Science Fund (11301258, 11361035), the Natural Science Fund of Inner Mongolia Autonomous Region (2012MS0108, 2012MS0106), the Scientific Research Projection of Higher Schools of Inner Mongolia (NJZZ12011, NJZY13199, and NJZY14013). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.amc.2014.06.023. References [1] D. Baleanu, A.H. Bhrawy, T.M. Taha, Two efficient generalized Laguerre spectral algorithms for fractional initial value problems, Abstract Appl. Anal. 2013 (2013). Article ID 546502, 10 pages. [2] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for two-sided space-fractional partial differential equations, Appl. Math. Comput. 56 (1) (2006) 80–90. [3] N. Zhang, W. Deng, Y. Wu, Finite difference/element method for a two-dimensional modified fractional diffusion equation, Adv. Appl. Math. Mech. 4 (2012) 496–518. [4] C. Tadjeran, M. Meerschaert, A second-order accurate numerical method for the two-dimensional fractional diffusion equation, J. Comput. Phys. 220 (2007) 813–823. [5] C. Tadjeran, M. Meerschaert, H. Scheffler, A second-order accurate numerical approximation for the fractional diffusion equation, J. Comput. Phys. 213 (2006) 205–213. [6] C.M. Chen, F. Liu, K. Burrage, Finite difference methods and a Fourier analysis for the fractional reaction-subdiffusion equation, Appl. Math. Comput. 198 (2) (2008) 754–769.
Y. Liu et al. / Applied Mathematics and Computation 243 (2014) 703–717
717
[7] Y. Liu, H. Li, W. Gao, S. He, Z.C. Fang, A new mixed element method for a class of time-fractional partial differential equations, Sci. World J. 2014 (2014). Article ID 141467, 8 pages. [8] Y.M. Lin, C.J. Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, J. Comput. Phys. 225 (2) (2007) 1533–1552. [9] P. Zhuang, F. Liu, I. Turner, Y.T. Gu, Finite volume and finite element methods for solving a one-dimensional space-fractional Boussinesq equation, Appl. Math. Model. (2013). in press, http://dx.doi.org/10.1016/j.apm.2013.10.008. [10] C.M. Chen, F. Liu, V. Anh, I. Turner, Numerical methods for solving a two-dimensional variable-order anomalous subdiffusion equation, Math. Comput. 81 (2012) 345–366. [11] H. Zhang, F. Liu, V. Anh, Galerkin finite element approximation of symmetric space-fractional partial differential equations, Appl. Math. Comput. 217 (6) (2010) 2534–2545. [12] R. Lin, F. Liu, V. Anh, I. Turner, Stability and convergence of a new explicit finite-difference approximation for the variable-order nonlinear fractional diffusion equation, Appl. Math. Comput. 212 (2) (2009) 435–445. [13] S.B. Yuste, Weighted average finite difference methods for fractional diffusion equations, J. Comput. Phys. 216 (1) (2006) 264–274. [14] E. Sousa, Finite difference approximations for a fractional advection diffusion problem, J. Comput. Phys. 228 (2009) 4038–4054. [15] L.L. Wei, Y.N. He, Analysis of a fully discrete local discontinuous Galerkin method for time-fractional fourth-order problems, Appl. Math. Model. 38 (4) (2014) 1511–1522. [16] C.P. Li, F.H. Zeng, Finite difference methods for fractional differential equations, Int. J. Bifurcation Chaos 22 (2012). 28 pages. [17] C.P. Li, Z.G. Zhao, Y.Q. Chen, Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion, Comput. Math. Appl. 62 (3) (2011) 855–875. [18] Z.G. Zhao, C.P. Li, Fractional difference/finite element approximations for the time-space fractional telegraph equation, Appl. Math. Comput. 219 (6) (2012) 2975–2988. [19] J.T. Ma, J.Q. Liu, Z.Q. Zhou, Convergence analysis of moving finite element methods for space fractional differential equations, J. Comput. Appl. Math. 255 (2014) 661–670. [20] Y.J. Jiang, J.T. Ma, High-order finite element methods for time-fractional partial differential equations, J. Comput. Appl. Math. 235 (11) (2011) 3285– 3290. [21] Y.J. Jiang, J.T. Ma, Moving finite element methods for time fractional partial differential equations, Sci. China Math. 56 (6) (2013) 1287–1300. [22] G. Gao, Z. Sun, A compact finite difference scheme for the fractional sub-diffusion equations, J. Comput. Phys. 230 (2011) 586–595. [23] Y. Zhang, Z. Sun, H. Wu, Error estimate of Crank–Nicolson-type difference for subdiffusion equation, SIAM J. Numer. Anal. 49 (2011) 2302–2322. [24] N.J. Ford, J.Y. Xiao, Y.B. Yan, A finite element method for time fractional partial differential equations, Fract. Calc. Appl. Anal. 14 (3) (2011) 454–474. [25] J.C. Li, Optimal convergence analysis of mixed finite element methods for fourth-order elliptic and parabolic problems, Numer. Methods Partial Differ. Equ. 22 (2006) 884–896. [26] D.Y. Shi, S.P. Mao, S.C. Chen, On the anisotropic accuracy analysis of ACM’s nonconforming finite element, J. Comput. Math. 23 (6) (2005) 635–646. [27] Y. Liu, Z.C. Fang, H. Li, S. He, W. Gao, A coupling method based on new MFE and FE for fourth-order parabolic equation, J. Appl. Math. Comput. 43 (2013) 249–269. [28] Z.D. Luo, Mixed Finite Element Methods and Applications, Chinese Science Press, Beijing, China, 2006. [29] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, The Netherlands, 1978. [30] K.X. Wang, H. Wang, A fast characteristic finite difference method for fractional advection-diffusion equations, Adv. Water Resour. 34 (7) (2011) 810–816.