- Email: [email protected]
A new expanded mixed method for parabolic integro-differential equations
Applied Mathematics and Computation 259 (2015) 600–613
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A new expanded mixed method for parabolic integro-differential equations q Yang Liu ⇑, Zhichao Fang, Hong Li *, Siriguleng He, Wei Gao School of Mathematical Sciences, Inner Mongolia University, Hohhot 010021, China
a r t i c l e
i n f o
Keywords: New expanded mixed method Parabolic integro-differential equations New expanded mixed projection Backward Euler scheme Error estimates
a b s t r a c t A new expanded mixed scheme is studied and analyzed for linear parabolic integrodifferential equations. The proposed method’s gradient belongs to the simple square integrable space replacing the classical Hðdiv;XÞ space. The new expanded mixed projection is introduced, the existence and uniqueness of solution for semi-discrete scheme are proved and the fully discrete error estimates based on both backward Euler scheme are obtained. Moreover, the optimal a priori error estimates in L2 and H1 -norm for the scalar unknown u and the error results in L2 ðXÞ-norm for its gradient k, and its flux r (the coefficients times the negative gradient) are derived. Finally, some numerical results are calculated to verify our theoretical analysis. Ó 2015 Elsevier Inc. All rights reserved.
1. Introduction We consider the following linear parabolic integro-differential equations
8 Rt > > < ut r aðx; tÞru þ 0 bðx; s; tÞruds ¼ f ðx; tÞ; ðx; tÞ 2 X J; uðx; tÞ ¼ 0; ðx; tÞ 2 @ X J; > > : uðx; 0Þ ¼ u0 ðxÞ; x 2 X;
ð1:1Þ
where X is a bounded convex polygonal domain in R2 with Lipschitz continuous boundary @ X; J ¼ ð0; T is the time interval with 0 < T < 1. u0 ðxÞ and f ðx; tÞ are given functions. Coefficients a ¼ aðx; tÞ and bðs; tÞ ¼ bðx; s; tÞ are two smooth and bounded functions, which satisfy the property that there exist some positive constants a0 ; a1 ; b0 , and b1 such that
0 < a0 6 aðx; tÞ 6 a1 < þ1; 0 < b0 6 bðs; tÞ 6 b1 < þ1: Parabolic integro-differential equations are a class of very important evolution equations which describe many physical phenomena such as heat conduction in material with memory, compression of viscoelastic media, nuclear reactor dynamics, etc. In recent years, a lot of scholars have studied the numerical methods for parabolic integro-differential equations, such as finite element methods [1–4], mixed finite element methods [5–10], finite volume element method [11] and so forth.
q Foundation item: Supported by the National Natural Science Fund (11301258, 11361035), the Natural Science Fund of Inner Mongolia Autonomous Region (2012MS0108, 2012MS0106, 2014BS0101), the Scientific Research Projection of Higher Schools of Inner Mongolia (NJZZ12011, NJZY14013), the Key Project of Chinese Ministry of Education (12024), the Program of Higher-Level Talents of Inner Mongolia University (30105-135127). ⇑ Corresponding authors. E-mail addresses: [email protected] (Y. Liu), [email protected] (H. Li).
http://dx.doi.org/10.1016/j.amc.2015.02.081 0096-3003/Ó 2015 Elsevier Inc. All rights reserved.
601
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
In recent years, many researchers have studied the Chen’s expanded mixed finite element method [12,13]. Compared to standard mixed element methods, the expanded mixed method is expanded in the sense that three variables are explicitly approximated, namely, the scalar unknown, its gradient, and its flux (the tensor coefficient times the negative gradient). In 1997, Arbogast et al. [14] derived and exploited a connection between the expanded mixed method and a certain cell-centered finite difference method. And Chen proved some mathematical theories for second-order quasilinear elliptic equation [15]. From then on, the expanded mixed method was applied to others evolution equations [16], and some new numerical methods based on the Chen’s expanded mixed method were proposed. For example, two-grid expanded mixed methods [17– 20], expanded upwind-mixed multi-step method [21], expanded characteristic-mixed finite element method [22], H1 Galerkin expanded mixed method [23–25], expanded mixed covolume method [26], expanded mixed hybrid methods [27], and positive definite expanded mixed method [9] and so on. From the above literatures on the study of expanded mixed method, we can find that all papers were studied based on Chen’s expanded mixed method [12,13]. For problem (1.1), we introduce the auxiliary variables k ¼ ru; r ¼ aðx; tÞru Rt Rt bðx; s; tÞruds ¼ aðx; tÞk 0 bðx; s; tÞkds and consider the Chen’s expanded mixed method [12,13], to obtain the following 0 expanded mixed weak formulation
8 > < ðaÞ ðut ; v Þ þ ðr r; v Þ ¼ ðf ðx; tÞ; v Þ; 8v 2 W; ðbÞ ðk; wÞ þ ðu; r wÞ ¼ 0; 8w 2 V; > Rt : ðcÞ ðr; zÞ þ ðak; zÞ þ ð 0 bðs; tÞkds; zÞ ¼ 0; 8z 2 X;
ð1:2Þ
2
2
where W ¼ L2 ðXÞ or W ¼ fw 2 L2 ðXÞjwj@ X ¼ 0g; V ¼ Hðdiv;XÞ ¼ fv 2 ðL2 ðXÞÞ jr v 2 L2 ðXÞg; X ¼ L2 ðXÞ ¼ ðL2 ðXÞÞ . In this article, we develop a new expanded mixed finite element method for linear integro-differential equations based on the new mixed weak formulation [28–32]. The proposed method is different from Chen’s expanded mixed scheme (1.2). Compared to expanded mixed scheme (1.2), our scheme mainly covers three advantages: The gradient for our method avoids the use of the classical Hðdiv;XÞ space and just needs to belong to the simple square integrable space L2 ðXÞ; The number of total degrees of freedom is reduced; The error in H1 -norm for the unknown function u can be derived. We will give the proof for the existence and uniqueness of the solution for semi-discrete scheme and a new expanded mixed projection and the proof of its uniqueness. We derived the fully discrete error estimates based on backward Euler scheme. Moreover, we prove the optimal priori error estimates in L2 and H1 -norm for the scalar unknown u and in L2 -norm for its gradient k, and its flux r (the coefficients times the gradient). Finally, we give some numerical experiments to confirm our theoretical results. Throughout this paper, C will denote a generic positive constant which does not depend on the spatial mesh parameter h 2
and temporal discretization parameter Dt. At the same time, we denote the natural inner product in L2 ðXÞ or L2 ðXÞ ¼ ðL2 ðXÞÞ by ð; Þ with the corresponding norm k k. The others notations and definitions of Sobolev spaces as in Ref. [33,34] are used. 2. New expanded mixed formulation Introduce the auxiliary variables k ¼ ru; r ¼ aðx; tÞru lowing first-order system for (1.1)
Rt 0
bðx; s; tÞruds ¼ aðx; tÞk
Rt 0
bðx; s; tÞkds to obtain the fol-
8 > < ut þ r r ¼ f ðx; tÞ; k ru ¼ 0; > R : r þ ak þ 0t bðs; tÞkds ¼ 0;
ð2:1Þ 2
2
The new expanded mixed weak formulation of (2.1) is to find fu; k; rg : ½0; T # H10 ðL2 ðXÞÞ ðL2 ðXÞÞ such that
8 > ðaÞ ðut ; v Þ ðr; rv Þ ¼ ðf ; v Þ; 8v 2 H10 ; > > < 2 ðbÞ ðk; wÞ ðru; wÞ ¼ 0; 8w 2 ðL2 ðXÞÞ ; > R > > : ðcÞ ðr; zÞ þ ðak; zÞ þ t bðs; tÞkds; z ¼ 0; 8z 2 ðL2 ðXÞÞ2 : 0
ð2:2Þ
Then, the semi-discrete mixed finite element scheme for (2.2) is to determine fuh ; kh ; rh g : ½0; T # V h Wh Wh such that
8 > < ðaÞ ðuht ; v h Þ ðrh ; rv h Þ ¼ ðf ; v h Þ; 8v h 2 V h ; ðbÞ ðkh ; wh Þ ðruh ; wh Þ ¼ 0; 8wh 2 Wh ; > Rt : ðcÞ ðrh ; zh Þ þ ðakh ; zh Þ þ ð 0 bðs; tÞkh ds; zh Þ ¼ 0; 8zh 2 Wh ;
ð2:3Þ
where ðV h ; Wh Þ is chosen as the finite element pair P 1 P 20 as follows
V h ¼ fv h 2 C 0 ðXÞ \ H10 jv h 2 P1 ðKÞ; 8K 2 Kh g; 2
Wh ¼ fwh ¼ ðw1h ; w2h Þ 2 ðL2 ðXÞÞ jw1h ; w2h 2 P0 ðKÞ; 8K 2 Kh g:
ð2:4Þ
602
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
We will prove the existence and uniqueness of solution for semi-discrete scheme (2.3). Theorem 2.1. There exists a unique discrete solution to semi-discrete scheme (2.3). 1 2 Proof. Let fwi ðxÞgni¼1 and fuj ðxÞgnj¼1 be bases of V h , Wh , respectively. Let
uh ¼
n1 X
ui ðtÞwi ðxÞ;
kh ¼
i¼1
n2 X
rh ¼
kj ðtÞuj ðxÞ;
j¼1
n2 X
rk ðtÞuk ðxÞ;
k¼1
and
uh ð0Þ ¼
n1 X i wi ðxÞ; u
kh ð0Þ ¼
n2 X kj u ðxÞ; j
i¼1
rh ð0Þ ¼
j¼1
n2 X
r k uk ðxÞ;
k¼1
and substitute these expressions into (2.3) and choose
8 ~ 0 ðtÞ Cr ~ ðtÞ ¼ FðtÞ; ðaÞ Au > > > > < ðbÞ D~kðtÞ Eu ~ ðtÞ ¼ 0; R > ~ ðtÞ þ H~kðtÞ þ 0t J~kðsÞds ¼ 0; ðcÞ Dr > > > : ~ r ~ ð0Þ ¼ u ; kð0Þ ~ ð0Þ ¼ r ; ðdÞu ¼ k;
v h ¼ wm ; wh ¼ ul ; zh ¼ ul to obtain the following equations ð2:5Þ
where
A ¼ ððwi ; wm ÞÞn1 n1 ; C ¼ ððus ; rwm ÞÞn1 n2 ; D ¼ ððuj ; ul ÞÞn
2 n2
H ¼ ððauj ; us ÞÞn
; E ¼ ððrwi ; ul ÞÞn2 n1 ;
2 n2
; J ¼ ððbuj ; us ÞÞn
2 n2
;
~ ðtÞ ¼ ðu1 ðtÞ; u2 ðtÞ; ; un1 ðtÞÞT ; u ~ ¼ ðk1 ðtÞ; k2 ðtÞ; ; kn ðtÞÞT ; kðtÞ 2
r~ ðtÞ ¼ ðr1 ðtÞ; r2 ðtÞ; ; rn2 ðtÞÞT ; ¼ ðu 1 ; u 2 u n1 ÞT ; FðtÞ ¼ ððf ; wm ÞÞn1 1 ; u k ¼ ðk1 ; k2 ; ; kn ÞT ; r 1; r 2; ; r n2 ÞT : ¼ ðr 2 It is easy to see that A and D are invertible matrixes. From (2.5), the initial value problems can be written as follows
8 ~ 0 ðtÞ A1 Cr ~ ðtÞ ¼ A1 FðtÞ; ðaÞ u > > > > < ðbÞ ~kðtÞ ¼ D1 Eu ~ ðtÞ; R 1 ~ > ~ ðtÞ þ D HkðtÞ þ D1 0t J~kðsÞds ¼ 0; > ðcÞ r > > : ~ r ~ ð0Þ ¼ u ; kð0Þ ~ ð0Þ ¼ r : ¼ k; ðdÞ u
ð2:6Þ
Substitute (2.6b) into (2.6c) to get
r~ ðtÞ ¼ D1 HD1 Eu~ ðtÞ D1
Z
t
~ ðsÞds: JD1 Eu
ð2:7Þ
0
Substituting (2.7) into (2.6a), we obtain
~ 0 ðtÞ þ A1 CðD1 HD1 Eu ~ ðtÞ þ D1 u
Z
t
~ ðsÞdsÞ ¼ A1 FðtÞ: JD1 Eu
ð2:8Þ
0
~ ðtÞ, then (2.7) and (2.6b) has a unique soluThus, by the theory of differential equations [35], (2.8) has a unique solution u ~ ~ ðtÞ and kðtÞ, respectively. Equivalently (2.3) has a unique solution. tion r 3. New expanded mixed projection and some Lemmas In order to analyze the convergence of the method, we first introduce the new expanded mixed elliptic projection associated with our equations. ~h ; ~ ~ h Þ : ½0; T ! V h Wh Wh be given by the following mixed relations kh ; r Let ðu
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
8 ~ h ; rv h Þ ¼ 0; 8v h 2 V h ; ðaÞ ðr r > > < ~ ~h Þ; wh Þ ¼ 0; 8wh 2 Wh ; ðbÞ ðk kh ; wh Þ ðrðu u R > > : ðcÞ ðr r ~ h ; zh Þ þ ðaðk ~kh Þ; zh Þ þ 0t bðs; tÞðk ~kh Þds; zh ¼ 0; 8zh 2 Wh :
603
ð3:1Þ
In what follows, we will give some important Lemmas based on new mixed scheme. 2
Lemma 3.1. There exists a linear operator Ph : ðL2 ðXÞÞ ! Wh such that
ðr Ph r; rv h Þ ¼ 0; 8v h 2 V h ;
ð3:2Þ
kr Ph rkL2 ðXÞ 6 ChkrkðH1 ðXÞÞ2 :
ð3:3Þ
and
Lemma 3.2. For the linear operator Ph of Lemma 3.1, we have
ðk Ph k; wh Þ ¼ 0; 8wh 2 Wh ;
ð3:4Þ
kk Ph kkL2 ðXÞ 6 ChkkkðH1 ðXÞÞ2 :
ð3:5Þ
and
Lemma 3.3. There exists a linear operator Ph : H10 ðXÞ ! V h such that
ðrðu Ph uÞ; wh Þ ¼ 0; 8wh 2 Wh ;
ð3:6Þ
and 2
ku Ph ukL2 ðXÞ þ hku Ph ukH1 6 Ch kukH2 ;
ð3:7Þ
2
ð3:8Þ
kut Ph ut kL2 ðXÞ 6 Ch kut kH2 : From Refs. [29,30], we can obtain the proof for Lemma 3.1–3.3. Using the definition of Ph and Ph , we rewrite g; d; q as
g ¼ u u~ h ¼ u Ph u þ Ph u u~ h ¼ gm þ ge ; d ¼ k ~kh ¼ k Ph k þ Ph k ~kh ¼ dm þ de ;
q ¼ r r~ h ¼ r Ph r þ Ph r r~ h ¼ qm þ qe : Since estimates of gm ; dm and qm are known, it is enough to estimate ge ; de and qe . Using Lemma 3.1–3.3, we rewrite (3.1) as
8 ðaÞ ðqe ; rv h Þ ¼ 0; 8v h 2 V h ; > > < ðbÞ ðde ; wh Þ ðrge ; wh Þ ¼ 0; 8wh 2 Wh ; R R > > : ðcÞ ðqe ; zh Þ þ ðade ; zh Þ þ t bðs; tÞde ds; zh ¼ ðadm ; zh Þ t bðs; tÞdm ds; zh ; 8zh 2 Wh : 0 0
ð3:9Þ
We discuss the following approximation properties for system (3.9). Lemma 3.4. There is a constant C independent of h such that
Z t kdk 6 Ch kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds ;
ð3:10Þ
Z t kqk 6 Ch krkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds ;
ð3:11Þ
Z t krgk 6 Ch kukH2 þ kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds :
ð3:12Þ
0
0
0
604
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
Proof. Choose
v h ¼ ge in (3.9a), wh ¼ qe in (3.9b) and zh ¼ de in (3.9c) to obtain
8 ðaÞ ðqe ; rge Þ ¼ 0; > > > > 0; > < ðbÞ ðde ; qe Þ ðrge ; qe Þ ¼ R t ðcÞ ðqe ; de Þ þ ðade ; de Þ þ 0 bðs; tÞde ds; de > > > > > : ¼ ðad ; d Þ R t bðs; tÞd ds; d : m e m e 0
ð3:13Þ
Add the three equations and use Cauchy–Schwarz inequality to get
Z t Z t a0 kde k2 6 ðade ; de Þ ¼ bðs; tÞde ds; de ðadm ; de Þ bðs; tÞdm ds; de 6C
0
Z
t
ðkde k þ kdm kÞds þ kdm k kde k:
0
ð3:14Þ
0
Using (3.14) and Gronwall lemma, we obtain
kde k 6 C
Z
t
kdm kds þ kdm k :
ð3:15Þ
0
Choose zh ¼ qe in (3.9c) and use Cauchy–Schwarz inequality to obtain
Z t kqe k 6 C kde k þ kdm k þ ðkde k þ kdm kÞds :
ð3:16Þ
0
Choose wh ¼ rge in (3.9b) and use (3.15) and Cauchy–Schwarz inequality to obtain
krge k 6 C
Z
t
kdm kds þ kdm k :
ð3:17Þ
0
Combining (3.15)–(3.17) and Lemma 3.1–3.3 and using the triangle inequality, we get the conclusion of Lemma 3.4. Lemma 3.5. There is a constant C independent of h such that 2
kgk 6 Ch kukH2 ;
ð3:18Þ
2
kgt k 6 Ch kut kH2 ;
ð3:19Þ
Z t kgk1 6 Ch kukH2 þ kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds :
ð3:20Þ
0
Proof. To estimate terms kgk; kgt k and kgk1 , we consider the following auxiliary elliptic problem
r ðarvÞ ¼ g; in X; v ¼ 0; on @ X:
ð3:21Þ
Use (3.21), Lemma 3.1–3.3 to obtain
kgk2
~h ; gÞ ¼ ðu u ~ h ; gÞ ¼ ðu Ph u þ Ph u u ~h ; gÞ ¼ ðu Ph u; gÞ þ ðPh u u ~h ; r ðarvÞÞ ¼ ðu Ph u; gÞ þ ðPh u u ~ h Þ; arvÞ ¼ ðu Ph u; gÞ þ ðrðPh u u ~ h Þ; arv Ph ðarvÞÞ ¼ ðu Ph u; gÞ þ ðrðPh u u ~ þðrðPh u uh Þ; Ph ðarvÞÞ
ð3:22Þ
¼ ðu Ph u; gÞ 6 ku Ph ukkgk 2
6 Ch kukH2 kgk: From (3.22), we obtain 2
kgk 6 Ch kukH2 :
ð3:23Þ
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
605
Using the similar method to kgkL2 , we can obtain 2
kgt k 6 Ch kut kH2 :
ð3:24Þ
Combining (3.12) and (3.23), we obtain
kgk1 6 ChðkukH2 þ kkkðH1 ðXÞÞ2 þ
Z
t
kkkðH1 ðXÞÞ2 dsÞ:
0
ð3:25Þ
Using (3.23)–(3.25), we obtain the conclusion of Lemma 3.5.
4. Semi-discrete error estimates In this section, we discuss a priori error estimates. For the error estimates, we decompose the errors as
~h þ u ~ h uh ¼ g þ 1; u uh ¼ u u k kh ¼ k ~kh þ ~kh kh ¼ d þ h;
r rh ¼ r r~ h þ r~ h rh ¼ q þ n: Using (2.2) and (2.3), we can get the error equations
8 ðaÞ ð1t ; v h Þ ðn; rv h Þ ¼ ðgt ; v h Þ; 8v h 2 V h ; > > < ðbÞ ðh; wh Þ ðr1; wh Þ ¼ 0; 8wh 2 Wh ; R > > : ðcÞ ðn; zh Þ þ ðah; zh Þ þ t bhds; zh ¼ 0; 8zh 2 Wh : 0
ð4:1Þ
We will prove the error estimates for semi-discrete scheme. ~ h ð0Þ, then one has the following estimates Theorem 4.1. Assume that uh ð0Þ ¼ u 2
ku uh k þ kut uht k 6 Ch
Z t kukH2 þ kut kH2 þ kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds ;
ð4:2Þ
0
Z t ku uh k1 6 Ch kukH2 þ kut kH2 þ kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds ;
ð4:3Þ
Z t kk kh k 6 Ch kukH2 þ kut kH2 þ kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds ;
ð4:4Þ
Z t kr rh k 6 Ch kukH2 þ kut kH2 þ krkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds :
ð4:5Þ
0
0
0
Proof. Choose
v h ¼ 1 in (4.1a), wh ¼ n in (4.1b) and zh ¼ h in (4.1c) to obtain
8 > ðaÞ 12 dtd k1k2 ðn; r1Þ ¼ ðgt ; 1Þ; > < ðbÞ ðh; nÞ ðr1; nÞ ¼ 0; > > : ðcÞ ðn; hÞ þ ka12 hk2 þ R t bhds; h ¼ 0:
ð4:6Þ
0
Adding the above three equations, and using Cauchy–Schwarz inequality, Young inequality, we have
1 d 1 k1k2 þ ka2 hk2 ¼ ðgt ; 1Þ 2 dt
Z 0
t
Z t bhds; h 6 C kgt k2 þ k1k2 þ khk2 þ khk2 ds :
ð4:7Þ
0
Integrate with respect to time from 0 to t to obtain
k1k2 þ
Z
t
khk2 ds 6 C
0
Z t Z s kgt k2 þ k1k2 þ khk2 þ khðsÞk2 ds ds: 0
ð4:8Þ
0
Using Gronwall Lemma, we obtain
k1k2 þ
Z
t 0
khk2 ds 6 C
Z 0
t
kgt k2 ds:
ð4:9Þ
606
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
Taking wh ¼ r1 in (4.1b), we have
kr1k 6 khk:
ð4:10Þ
Differentiating (4.1b) and taking wh ¼ n, we obtain
ðht ; nÞ ðr1t ; nÞ ¼ 0:
v h ¼ 1t
Choose
ð4:11Þ
in (4.1a) and zh ¼ ht in (4.1c) to obtain
8 < ðaÞ k1t k2 ðn; r1t Þ ¼ ðgt ; 1t Þ; : ðcÞ ðn; ht Þ þ 1 2
1 d ka2 hk2 dt
¼ 12 ðat h; hÞ þ
R t d dt
0
R t bhds; h dtd 0 bhds; h :
ð4:12Þ
Adding (4.11), (4.12a) and (4.12c), and using Cauchy–Schwarz inequality, Young inequality, we have
k1t k2 þ
Z t Z t 1 d 1 2 1 d ka2 hk ¼ ðgt ; 1t Þ þ ðat h; hÞ þ bt hds; h þ bðx; t; tÞh bhds; h 2 dt 2 dt 0 0 Z t Z t 1 d khk2 ds þ k1t k2 bhds; h : 6 C kgt k2 þ khk2 þ 2 dt 0 0
ð4:13Þ
Integrate with respect to time from 0 to t and use Gronwall lemma to obtain
Z
t
0
Z
k1t k2 ds þ khk2 6 C
t
0
kgt k2 ds:
ð4:14Þ
Choosing zh ¼ n in (4.1c) and using (4.9) and (4.14), we have
1 knk2 ¼ ðah; nÞ 2
Z
t
0
Z t Z t 1 bhds; n knk2 6 C khk2 þ khk2 ds 6 C kgt k2 ds: 2 0 0
ð4:15Þ
Combining (3.11), (3.12), (4.9), (4.10), (4.14), (4.15) and the triangle inequality, we obtain the error estimates for Theorem 4.1.
5. Fully discrete error estimates In this section, we get the error estimates of fully discrete schemes based on backward Euler method. For the fully discrete procedure, let 0 ¼ t 0 < t1 < t2 < < t M ¼ T be a given partition of the time interval ½0; T with step length Dt ¼ T=M and nodes tn ¼ nDt, for some positive integer M. For a smooth function / on ½0; T, define /n ¼ /ðtn Þ and @t /n ¼ ð/n /n1 Þ=Dt. For approximating the integrals, we use the composite left rectangle rule
Dt
Z n1 X /j
tn
/ðsÞds:
0
j¼0
Note that / 2 C 1 ½0; T, the quadrature error satisfies
Z tn Z tn X n1 j /ðsÞds 6 C Dt j/s ðsÞjds: Dt / j¼0 0 0 The Eq. (2.2) has the following equivalent formulation
8 n > ðaÞ ð@ t un ; v Þ ðrn ; rv Þ ¼ ðf þ Rn1 ; v Þ; 8v 2 H10 ; > > > > < ðbÞ ðkn ; wÞ ðrun ; wÞ ¼ 0; 8w 2 ðL2 ðXÞÞ2 ; ! n1 > X > 2 j > n n n > ðcÞ ð ¼ ðRn2 ; zÞ; 8z 2 ðL2 ðXÞÞ ; r ; zÞ þ ða k ; zÞ þ D t bðx; t ; t Þk ; z > n j : j¼0
where
Rn1 ¼ @ t un ut ðt n Þ ¼
Rn2 ¼ Dt
n1 X j¼0
1 Dt
Z
bðx; t n ; tj Þk j
tn
ðtn1 sÞutt ds;
t n1
Z
tn
bðx; t n ; sÞkðsÞds:
0
Now we formulate a fully discrete procedure: Find ðunh ; knh ; rnh Þ 2 V h Wh Wh ; ðn ¼ 0; 1; ; MÞ such that
ð5:1Þ
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
8 n ðaÞ ð@ t unh ; v h Þ ðrnh ; rv h Þ ¼ ðf ; v h Þ; 8v h 2 V h ; > > > > < ðbÞ ðknh ; wh Þ ðrunh ; wh Þ ¼ 0; 8wh 2 Wh ; ! n1 X > > j n n n > > : ðcÞ ðrh ; zh Þ þ ða kh ; zh þ Dt bðx; t n ; tj Þkh ; zh ¼ 0; 8zh 2 Wh :
607
ð5:2Þ
j¼0
For the fully discrete error estimates, we now split the errors
~nh þ u ~ nh unh ¼ gn þ 1n ; uðtn Þ unh ¼ uðtn Þ u kðtn Þ knh ¼ kðtn Þ ~knh þ ~knh knh ¼ dn þ hn ;
rðtn Þ rnh ¼ rðtn Þ r~ nh þ r~ nh rnh ¼ qn þ nn : We will prove the theorem for the fully discrete error estimates. ~ h ð0Þ, then there exists a positive constant C independent of h and Dt such that Theorem 5.1. Assume that u0h ¼ u
kuðt J Þ uJh k 6 Ch
2
Z kukL1 ðH2 Þ þ kut kL1 ðH2 Þ þ kkkL1 ððH1 Þ2 Þ þ C Dt
tJ
t0
Z kuðt J Þ uJh k1 6 Ch kukL1 ðH2 Þ þ kut kL1 ðH2 Þ þ kkkL1 ððH1 Þ2 Þ þ C Dt Z
ð5:3Þ
kkkðH1 Þ2 ds;
ð5:4Þ
tJ
t0
kkðt J Þ kJh k 6 ChðkukL1 ðH2 Þ þ kut kL1 ðH2 Þ þ kkkL1 ððH1 Þ2 Þ Þ þ C Dt
kkkðH1 Þ2 ds;
tJ
t0
kkkðH1 Þ2 ds;
ð5:5Þ
Z krðt J Þ rJh k 6 Ch kukL1 ðH2 Þ þ kut kL1 ðH2 Þ þ kkkL1 ððH1 Þ2 Þ þ krkL1 ððH1 Þ2 Þ þ C Dt
tJ
t0
kkkðH1 Þ2 ds:
ð5:6Þ
Proof. Using (3.1), (5.1) and (5.2) at t ¼ tn , we get the error equations
8 ðaÞ ð@ t 1n ; v h Þ ðnn ; rv h Þ ¼ ð@ t gn ; v h Þ þ ðRn1 ; v h Þ; 8v h 2 V h ; > > > > < ðbÞ ðhn ; wh Þ ðr1n ; wh Þ ¼ 0; 8wh 2 Wh ; ! n1 X > > n n j n n > > : ðcÞ ðn ; zh Þ þ ða h ; zh Þ þ Dt bðx; tn ; t j Þh ; zh ¼ ðR2 ; zh Þ; 8zh 2 Wh :
ð5:7Þ
j¼0
Choose
v h ¼ 1n
in (5.7a), wh ¼ nn in (5.7b) and zh ¼ hn in (5.7c) to obtain
8 > ðaÞ 21Dt ðk1n k2 k1n1 k2 þ k1n 1n1 k2 Þ > > > n n > n n n n > > < ðn ; r1 Þ ¼ ð@ t g ; 1 Þ þ ðR1 ; 1 Þ; n n n ðbÞ ðh ; n Þ ðr1n ; n Þ ¼ 0; ! > > n1 > X > 1 n n j n > n 2 n 2 > ðcÞ ðn ¼ ðRn2 ; hn Þ: ; h Þ þ kða Þ h k þ D t bðx; t ; t Þh ; h > n j :
ð5:8Þ
j¼0
Adding the above three equations, we obtain n1 X 1 1 ðk1n k2 k1n1 k2 Þ þ kðan Þ2 hn k2 ¼ ð@ t gn ; 1n Þ Dt bðx; tn ; t j Þh j ; hn 2Dt j¼0 n 2
6 C k@ t g k þ
kRn1 k2
n 2
þ k1 k þ Dt
!
n1 X j¼0
þ ðRn1 ; 1n Þ þ ðRn2 ; hn Þ ! j 2
kh k
þ
a0 n 2 kh k : 2
ð5:9Þ
Multiplying by 2Dt and summing (5.9) from n ¼ 1 to J, the resulting equation becomes
! J J J1 n1 X X X a0 X n 2 n 2 n 2 j 2 0 2 n 2 n 2 ð1 C DtÞk1 k þ Dt kh k 6 k1 k þ C Dt k@ t g k þ kR1 k þ kR2 k þ C Dt k1 k þ Dt kh k : 2 n¼1 n¼1 n¼1 j¼0 J 2
ð5:10Þ
608
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
Choose Dt0 in such a way that for 0 < Dt 6 Dt 0 ; ð1 C DtÞ > 0. Then an application of Gronwall’s lemma to obtain
k1J k2 þ a0 Dt
J J X X khn k2 6 C Dt k@ t gn k2 þ kRn1 k2 þ kRn2 k2 : n¼1
ð5:11Þ
n¼1
Note that
ðaÞ kRn1 k2 6 C Dt
R tn
kutt k2 ds; R t ðbÞ kRn2 k2 6 CðDtÞ2 0n kkk2 þ kkt k2 ds; R tn kgt ðsÞk2 ds: ðcÞ k@ t gn k2 6 C D1t tn1 tn1
ð5:12Þ
Substitute (5.12) into (5.11) to get
k1J k2 þ a0 Dt
Z J X khn k2 6 C
tJ t0
n¼1
kgt ðsÞk2 ds þ CðDtÞ2
Z
tJ
kutt k2 þ kkk2 þ kkt k2 ds:
ð5:13Þ
t0
Taking wh ¼ r1n in (5.7b), we have
kr1n k 6 khn k:
ð5:14Þ
From (5.7b), we get
ð@ t hn ; wh Þ ðr@ t 1n ; wh Þ ¼ 0: Choose wh ¼ nn in (5.15),
ð5:15Þ
v h ¼ @ t 1n in (5.7a) and zh ¼ @ t hn in (5.7c) to obtain
8 ðaÞ k@ t 1n k2 ðnn ; r@ t 1n Þ ¼ ð@ t gn ; @ t 1n Þ þ ðRn1 ; @ t 1n Þ; > > > > < ðbÞ ð@ t hn ; nn Þ ðr@ t 1n ; nn Þ ¼ 0; ! n1 X > > n n n n n j n n > > : ðcÞ ðn ; @ t h Þ þ ðaðx; t n Þ@ t h ; h Þ þ Dt bðx; t n ; tj Þh ; @ t h ¼ ðR2 ; @ t h Þ:
ð5:16Þ
j¼0
Adding the three equations, we obtain
1 an1 an 1 1 n 12 n 2 kða Þ h k kðan1 Þ2 hn1 k2 þ hn1 ; hn1 2 Dt 2 Dt ! n1 X 6 ð@ t gn ; @ t 1n Þ þ ðRn1 ; @ t 1n Þ Dt bðx; t n ; t j Þh j ; @ t hn þ ðRn2 ; @ t hn Þ:
k@ t 1n k2 þ
ð5:17Þ
j¼0
Note that
ðaÞðRn2 ; @ t hn Þ ¼
ðbÞ Dt
ðhn ;Rn2 Þðhn1 ;Rn1 2 Þ Dt
n1 X bðx; tn ; t j Þh j ; @ t hn
!
ð@ t Rn2 ; hn1 Þ; P Dt
n1 j¼0
Dt
!
n2 X
bðx;t n1 ;t j Þh j ;hn1
j¼0
¼
j¼0
bðx;t n ;t j Þh j ;hn
Dt
P bðx;t n ;t j Þbðx;t n1 ;t j Þ j Dt n1 h ; hn1 j¼0 Dt
ð5:18Þ
bðx; t n1 ; tn1 Þhn1 ; hn1 :
Substitute (5.18) into (5.17) to obtain
k@ t 1n k2 þ
1 an1 an 1 1 n 12 n 2 kða Þ h k kðan1 Þ2 hn1 k2 þ hn1 ; hn1 2 Dt 2 Dt n
n
6 ð@ t g ; @ t 1 Þ þ þ
ðRn1 ; @ t n Þ
1 þ Dt
n1 X bðx; t n ; t j Þ bðx; tn1 ; t j Þ j¼0
Dt
Pn1 j¼0
Dt
P j n1 bðx; t n ; tj Þh j ; hn Dt n2 j¼0 bðx; t n1 ; t j Þh ; h
Dt ðhn ; Rn2 Þ ðhn1 ; Rn1 2 Þ ð@ t Rn2 ; hn1 Þ: Dt
j
h ;h
!
n1
þ ðbðx; t n1 ; t n1 Þhn1 ; hn1 Þ ð5:19Þ
609
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
Multiplying by 2Dt and summing (5.19) from n ¼ 1 to J, the resulting equation becomes
Dt
J J J X X X 1 k@ t 1n k2 þ ða0 C DtÞkhJ k2 6 kða0 Þ2 h0 k2 þ C Dt k@ t gn k2 þ kRn1 k2 þ k@ t Rn2 k2 þ C Dt khn1 k2 n¼0
n¼1
n¼1
J X n1
X
bðx; tn ; t j Þ bðx; t n1 ; t j Þ j 2 n 2
þ ðDtÞ3 h
þ kR2 k Dt n¼1 j¼0
1 2
6 kða0 Þ h0 k2 þ C Dt
J J X X k@ t gn k2 þ kRn1 k2 þ k@ t Rn2 k2 þ C Dt khn1 k2 þ kRn2 k2 : n¼1
ð5:20Þ
n¼1
Using (5.12), we obtain
ðaÞ Dt
J X Rt kR1j k2 6 CðDtÞ2 t0J kutt k2 ds; j¼1
ð5:21Þ
J X Rt ðbÞ Dt k@ t g j k2 6 C t0J kgt k2 ds: j¼1
Now, we estimate the term Dt
@ t Rn2 ¼ Dt
PJ
j 2 j¼1 k@ t R2 k .
Z n2 X ½bðtn ; t j Þ bðtn1 ; tj Þk j þ Dtbðt n ; tn1 Þkn1
bðt n ; sÞkðsÞds
t n1
j¼0
¼ ðDtÞ2
tn
n2 X bt ðt n0 ; t j Þk j þ Dt½bðtn ; t n1 Þkn1 bðtn ; t n1 Þkðtn1 Þ j¼0
n2 X ¼ ðDtÞ2 bt ðt n0 ; t j Þk j þ Dt½bðtn ; t n1 Þkt ðt n2 Þðt n1 t n1 Þ þ bt ðt n ; tn3 Þðtn1 t n1 Þkðt n1 Þ; j¼0
where tn0 ; t n1 2 ðtn1 ; t n Þ; t n2 ; tn3 2 ðt n1 ; t n1 Þ: From (5.22), we obtain
j@ t Rn2 j 6 Dt max jbt ðt n0 ; t j Þk j j þ ðDtÞ2 ½jbðt n ; tn1 Þkt ðt n2 Þj þ jbt ðt n ; t n3 Þkðt n1 Þj:
ð5:23Þ
06j6n2
Using (5.23), we obtain
Dt
J X k@ t R2j k2 6 ðDtÞ2 kkk2L1 ðL2 Þ þ kkt k2L1 ðL2 Þ :
ð5:24Þ
j¼1
Substitute (5.21), (5.24) and (5.12) into (5.20) and use Gronwall lemma to get
2Dt
Z J X k@ t 1n k2 þ khJ k2 6 C
tJ
t0
n¼0
kgt ðsÞk2 ds þ CðDtÞ2 kutt k2L2 ðL2 Þ þ kkk2L1 ðL2 Þ þ kkt k2L1 ðL2 Þ þ kkk2L2 ðL2 Þ þ kkt k2L2 ðL2 Þ :
Taking zh ¼ n in (5.7c), we get
knn k2 ¼ ðan hn ; nn Þ ðDt
n1 n1 X X bðx; t n ; tj Þh j ; nn Þ þ ðRn2 ; nn Þ 6 C khn k2 þ ðDtÞ2 kh j k2 j¼0
j¼0
!
1 þ knn k2 : 2
ð5:25Þ
ð5:26Þ
Substitute (5.25) into (5.26) to get
knJ k2 6 C
Z
tJ
t0
kgt ðsÞk2 ds þ CðDtÞ2 kutt k2L2 ðL2 Þ þ kkk2L1 ðL2 Þ þ kkt k2L1 ðL2 Þ þ kkk2L2 ðL2 Þ þ kkt k2L2 ðL2 Þ :
ð5:27Þ
Combining (3.11), (3.12), (5.13), (5.14), (5.25), (5.27) and the triangle inequality, we complete the proof. 6. Numerical experiments For confirming our theoretical analysis, we need to consider a numerical example based on problem (1.1) with initial boundary conditions. In (1.1), the spatial domain and the temporal interval are X ¼ ð0; 1Þ ð0; 1Þ and J ¼ ð0; 1, respectively. The known source term f ðx; tÞ is decided by the chosen coefficients aðx; tÞ ¼ 1 þ x21 þ x22 þ t2 ; bðx; t; sÞ ¼ 1 þ x21 þ x22 þ t 2 þ s and the exact solution
uðx; tÞ ¼ et sinðpx1 Þ sinð2px2 Þ;
610
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
the exact gradient function
k ¼ ru ¼ ðpet cosðpx1 Þ sinð2px2 Þ; 2pet sinðpx1 Þ cosð2px2 ÞÞ; and the exact flux function
r ¼ aru
Z
t
bðx; t; sÞruds ¼ ðdðx; tÞ cosðpx1 Þ sinð2px2 Þ; 2dðx; tÞ sinðpx1 Þ cosð2px2 ÞÞ;
0
where x ¼ ðx1 ; x2 Þ and dðx; tÞ ¼ pð2 þ x21 þ x22 þ t2 tet et Þ. We first divide the spatial domain X ¼ ½0; 1 ½0; 1 into uniform triangulations with spatial mesh size h, use the backward Euler approximation with uniform time step length Dt, choose the piecewise linear space V h with index k ¼ 1 and the piecewise constant space Wh with index k ¼ 0, then obtain the fully discrete procedure. Based on the above consideration, we will obtain some calculated data on orders of convergence and a priori errors. pffiffi pffiffi pffiffi pffiffi With a fixed time step Dt ¼ 1=100 and some changed space meshes length h ¼ 42 ; 82 ; 162 ; 322, we obtain some calculated results of a priori error estimates in L2 and H1 -norms for the scalar unknown u in Table 1. From Table 1, we find that the convergence rate in L2 -norm for u is close to 2, while the one in H1 -norm is close to 1. So we claim that these computed convergence results, which are in correspondence with our theoretical results, are optimal. With the same choices of space–time 2
meshes to the case for u, we also calculate a priori error results and orders of convergence in L2 ðXÞ ¼ ðL2 ðXÞÞ -norm for the gradient k and its flux r in Table 2. From the calculated results, we see that the convergence rate, which is also in correspondence with the results of theorem, is close to 1.
Table 1 Error estimates and convergence rate in L2 and H1 -norms for u with Dt ¼ 1=100. h pffiffi 2 4 pffiffi 2 8 pffiffi 2 16 pffiffi 2 32
ku uh kL1 ðL2 ðXÞÞ
Rate
ku uh kL1 ðH1 ðXÞÞ
1.7778e001
Rate
1.2816e+000
4.7102e002
1.9162
6.2158e001
1.0439
1.1941e002
1.9799
3.0833e001
1.0115
3.0125e003
1.9869
1.5384e001
1.0030
Rate
kr rh kL1 ðL2 ðXÞÞ
Rate
6.2068e001
1.0366
1.8693e+000
1.0416
3.0822e001
1.0099
9.2991e001
1.0073
1.5383e001
1.0026
4.6435e001
1.0019
Table 2 Error estimates and convergence rate in L2 -norm for gradient k and its flux r with Dt ¼ 1=100. h pffiffi 2 4 pffiffi 2 8 pffiffi 2 16 pffiffi 2 32
kk kh kL1 ðL2 ðXÞÞ 1.2733e+000
3.8480e+000
0.4 0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 1
1 0.8
0.6
0.5 0.4
0.2
0
Fig. 1. The exact solution u.
0
611
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
0.4 0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 1
1 0.8
0.6
0.5 0.4
0.2
0
0
Fig. 2. The numerical solution uh .
2
4
1
2
0
0 1
−1 −2 1
1
−2 −4 1
0.5 0.5
0.5 0.5
0 0
0 0
Fig. 3. The exact gradient function k ¼ ðk1 ; k2 Þ.
2
4
1
2
0
0 1
−1 −2 1
0.5 0.5
1
−2 −4 1
0.5 0.5
0 0
Fig. 4. The numerical gradient function kh ¼ ðk1h ; k2h Þ.
0 0
612
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
20
20
10
0
0 1
−10 −20 1
−20
1
−40 1
0.5
0.5
0.5
0.5 0 0
0 0
Fig. 5. The exact flux function
20
r ¼ ðr1 ; r2 Þ.
20
10
0
0 1
−10 −20 1
0.5
−20
1
−40 1
0.5
0.5 0.5
0 0
Fig. 6. The numerical flux function
0 0
rh ¼ ðr1h ; r2h Þ.
Now we will make some comparisons of figures between the exact solution and the numerical solution. We show the exact unknown function u in Fig. 1 at time t ¼ 1 with the chosen time step length Dt ¼ 1=100 and the spatial mesh size pffiffi h ¼ 162. Under the same conditions, we also show the numerical solution uh in Fig. 2. Form Figs. 1 and 2, it is not hard to find that the numerical function uh can primely approximate the exact function u. The similar comparisons are made between the exact gradient solution k and the numerical gradient solution kh in Figs. 3 and 4. Finally, we also give the comparison of figures between the exact solution r of flux function and the numerical solution rh of flux function in Figs. 5,6. Based on the discussion of the calculated results in Tables 1, 2 and Figs. 1–6, we can confirm our theoretical analysis.
7. Concluding remarks In this article, we proposed and analyzed a new expanded mixed finite element scheme, whose gradient belongs to the square integrable space instead of the classical Hðdiv;XÞ, for linear parabolic integro-differential equations. We introduced the new expanded mixed projection and proved some important Lemmas based on new expanded mixed projection. We derived the optimal priori error estimates in L2 for the scalar unknown u and the suboptimal error estimates in L2 -norm for its gradient k, and its flux r (the coefficients times the negative gradient). Moreover, we can obtained the H1 -norm for the scalar unknown u, which cannot be obtained by the previous expanded mixed finite element methods. In the near future, the proposed expanded mixed scheme will be applied to other equations, such as nonlinear parabolic equation, hyperbolic wave equation and so on.
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
613
Acknowledgment Authors thank the anonymous reviewers and editor for their valuable comments and suggestions, which greatly improve this work. 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. 2015.02.081. References [1] T. Zhang, Finite Element Methods For Partial Differential-Integral Equations, Chinese Science Press, Beijing, 2009. [2] A.K. Pani, T.E. Peterson, Finite element methods with numerical quadrature for parabolic integro-differential equations, SIAM J. Numer. Anal. 33 (3) (1996) 1084–1105. [3] A.K. Pani, R.K. Sinha, Error estimates for semidiscrete Galerkin approximation to a time dependent parabolic integro-differential equation with nonsmooth data, Calcolo 37 (4) (2000) 181–205. [4] Y. Lin, V. Thomée, L.B. Wahlbin, Ritz–Volterra projections to finite-element spaces and applications to integro-differential and related equations, SIAM J. Numer. Anal. 28 (4) (1991) 1047–1070. [5] R.K. Sinha, R.E. Ewing, R.D. Lazarov, Mixed finite element approximations of parabolic integro-differential equations with nonsmooth initial data, SIAM J. Numer. Anal. 47 (5) (2009) 3269–3292. [6] Z. Jiang, L1 ðL2 Þ and L1 ðL1 Þ error estimates for mixed methods for integro-differential equations of parabolic type, Math. Model. Numer. Anal. 33 (3) (1999) 531–546. [7] H. Guo, H.X. Rui, Least-squares Galerkin procedures for parabolic integro-differential equations, Appl. Math. Comput. 150 (3) (2004) 749–762. [8] D.Y. Shi, H.H. Wang, An H1 -Galerkin nonconforming mixed finite element method for integro-differential equation of parabolic type, J. Math. Res. Exposition 29 (5) (2009) 871–881. [9] Y. Liu, H. Li, J.F. Wang, W. Gao, A new positive definite expanded mixed finite element method for parabolic integrodifferential equations, J. Appl. Math. 2012, Article ID 391372, 24pages. [10] H. Guo, J.S. Zhang, H.F. Fu, Two splitting positive definite mixed finite element methods for parabolic integro-differential equations, Appl. Math. Comput. 218 (22) (2012) 11255–11268. [11] H.R. Li, Q. Li, Finite volume element methods for nonlinear parabolic integro-differential problems, J. KSIAM 7 (2) (2003) 35–49. [12] Z.X. Chen, Expanded mixed finite element methods for linear second order elliptic problems I, IMA preprint series ] 1219, Institute for Mathematics and its Applications, University of Minnesota, Minneapolis, MN, 1994. [13] Z.X. Chen, Expanded mixed element methods for linear second-order elliptic problems(I), RAIRO Model. Math. Anal. Numer. 32 (1998) 479–499. [14] T. Arbogast, M.F. Wheeler, I. Yotov, Mixed finite elements for elliptic problems with tensor coefficients as cellcentered finite differences, SIAM J. Numer. Anal. 34 (1997) 828–852. [15] Z.X. Chen, Expanded mixed element methods for quasilinear second-order elliptic problems(II), RAIRO Model. Math. Anal. Numer. 32 (1998). 501–420. [16] C.S. Woodward, C.N. Dawson, Analysis of expanded mixed finite element methods for a non-linear parabolic equation modelling flow into variably saturated porous media, SIAM J. Numer. Anal. 37 (2000) 701–724. [17] Y.P. Chen, Y.Q. Huang, D.H. Yu, A two-grid method for expanded mixed finite-element solution of semilinear reaction-diffusion equations, Int. J. Numer. Methods Eng. 57 (2) (2003) 193–209. [18] Y.P. Chen, H.W. Liu, S. Liu, Analysis of two-grid methods for reaction–diffusion equations by expanded mixed finite element methods, Int. J. Numer. Methods Eng. 69 (2) (2007) 408–422. [19] Y.P. Chen, P. Luan, Z.L. Lu, Analysis of two-grid methods for nonlinear parabolic equations by expanded mixed finite element methods, Adv. Appl. Math. Mech. 1 (6) (2009) 830–844. [20] W. Liu, H.X. Rui, H. Guo, A two-grid method with expanded mixed element for nonlinear reaction–diffusion equations, Acta Math. Appl. Sin. 27 (3) (2011) 495–502. [21] H.L. Song, Y.R. Yuan, The expanded upwind-mixed multi-step method for the miscible displacement problem in three dimensions, Appl. Math. Comput. 195 (2008) 100–109. [22] L. Guo, H.Z. Chen, An expanded characteristic-mixed finite element method for a convection-dominated transport problem, J. Comput. Math. 23 (5) (2005) 479–490. [23] H.Z. Chen, H. Wang, An optimal-order error estimate on an H1 -Galerkin mixed method for a nonlinear parabolic equation in porous medium flow, Numer. Methods Partial Differ. Equ. 26 (2010) 188–205. [24] H.T. Che, Y.J. Wang, Z.J. Zhou, An optimal error estimates of H1 -Galerkin expanded mixed finite element methods for nonlinear viscoelasticity-type equation, Math. Prob. Eng. 2011, Article ID 570980, 18 pages doi:10.1155/2011/570980. [25] Y. Liu, Analysis and numerical simulation of nonstandard mixed element methods (Ph.D. thesis), Inner Mongolia University, Hohhot, China, 2011. [26] H.X. Rui, T.C. Lu, An expanded mixed covolume method for elliptic problems, Numer. Methods Partial Differ. Equ. 21 (1) (2005) 8–23. [27] D. Kim, E.J. Park, A posteriori error estimator for expanded mixed hybrid methods, Numer. Methods Partial Differ. Equ. 23 (2007) 330–349. [28] Y. Liu, H. Li, et al., A new expanded mixed finite element method for elliptic equations, submitted for publication. [29] S.C. Chen, H.R. Chen, New mixed element schemes for second order elliptic problem, Math. Numer. Sin. 32 (2) (2010) 213–218. [30] F. Shi, J.P. Yu, K.T. Li, A new stabilized mixed finite-element method for Poisson equation based on two local Gauss integrations for linear element pair, Int. J. Comput. Math. 88 (11) (2011) 2293–2305. [31] D.Y. Shi, Y.D. Zhang, High accuracy analysis of a new nonconforming mixed finite element scheme for Sobolev equations, Appl. Math. Comput. 218 (7) (2011) 3176–3186. [32] Z.F. Weng, X.L. Feng, P.Z. Huang, A new mixed finite element method based on the Crank–Nicolson scheme for the parabolic problems, Appl. Math. Model. (2012), http://dx.doi.org/10.1016/j.apm.2011.12.044. [33] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [34] Z.D. Luo, Mixed Finite Element Methods and Applications, Chinese Science Press, Beijing, China, 2006. [35] H. Brunner, P.J. Houwen, The Numerical Solution of Volterra Equations, North-Holland Publishing Co., Amsterdam, 1986.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A new expanded mixed method for parabolic integro-differential equations q Yang Liu ⇑, Zhichao Fang, Hong Li *, Siriguleng He, Wei Gao School of Mathematical Sciences, Inner Mongolia University, Hohhot 010021, China
a r t i c l e
i n f o
Keywords: New expanded mixed method Parabolic integro-differential equations New expanded mixed projection Backward Euler scheme Error estimates
a b s t r a c t A new expanded mixed scheme is studied and analyzed for linear parabolic integrodifferential equations. The proposed method’s gradient belongs to the simple square integrable space replacing the classical Hðdiv;XÞ space. The new expanded mixed projection is introduced, the existence and uniqueness of solution for semi-discrete scheme are proved and the fully discrete error estimates based on both backward Euler scheme are obtained. Moreover, the optimal a priori error estimates in L2 and H1 -norm for the scalar unknown u and the error results in L2 ðXÞ-norm for its gradient k, and its flux r (the coefficients times the negative gradient) are derived. Finally, some numerical results are calculated to verify our theoretical analysis. Ó 2015 Elsevier Inc. All rights reserved.
1. Introduction We consider the following linear parabolic integro-differential equations
8 Rt > > < ut r aðx; tÞru þ 0 bðx; s; tÞruds ¼ f ðx; tÞ; ðx; tÞ 2 X J; uðx; tÞ ¼ 0; ðx; tÞ 2 @ X J; > > : uðx; 0Þ ¼ u0 ðxÞ; x 2 X;
ð1:1Þ
where X is a bounded convex polygonal domain in R2 with Lipschitz continuous boundary @ X; J ¼ ð0; T is the time interval with 0 < T < 1. u0 ðxÞ and f ðx; tÞ are given functions. Coefficients a ¼ aðx; tÞ and bðs; tÞ ¼ bðx; s; tÞ are two smooth and bounded functions, which satisfy the property that there exist some positive constants a0 ; a1 ; b0 , and b1 such that
0 < a0 6 aðx; tÞ 6 a1 < þ1; 0 < b0 6 bðs; tÞ 6 b1 < þ1: Parabolic integro-differential equations are a class of very important evolution equations which describe many physical phenomena such as heat conduction in material with memory, compression of viscoelastic media, nuclear reactor dynamics, etc. In recent years, a lot of scholars have studied the numerical methods for parabolic integro-differential equations, such as finite element methods [1–4], mixed finite element methods [5–10], finite volume element method [11] and so forth.
q Foundation item: Supported by the National Natural Science Fund (11301258, 11361035), the Natural Science Fund of Inner Mongolia Autonomous Region (2012MS0108, 2012MS0106, 2014BS0101), the Scientific Research Projection of Higher Schools of Inner Mongolia (NJZZ12011, NJZY14013), the Key Project of Chinese Ministry of Education (12024), the Program of Higher-Level Talents of Inner Mongolia University (30105-135127). ⇑ Corresponding authors. E-mail addresses: [email protected] (Y. Liu), [email protected] (H. Li).
http://dx.doi.org/10.1016/j.amc.2015.02.081 0096-3003/Ó 2015 Elsevier Inc. All rights reserved.
601
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
In recent years, many researchers have studied the Chen’s expanded mixed finite element method [12,13]. Compared to standard mixed element methods, the expanded mixed method is expanded in the sense that three variables are explicitly approximated, namely, the scalar unknown, its gradient, and its flux (the tensor coefficient times the negative gradient). In 1997, Arbogast et al. [14] derived and exploited a connection between the expanded mixed method and a certain cell-centered finite difference method. And Chen proved some mathematical theories for second-order quasilinear elliptic equation [15]. From then on, the expanded mixed method was applied to others evolution equations [16], and some new numerical methods based on the Chen’s expanded mixed method were proposed. For example, two-grid expanded mixed methods [17– 20], expanded upwind-mixed multi-step method [21], expanded characteristic-mixed finite element method [22], H1 Galerkin expanded mixed method [23–25], expanded mixed covolume method [26], expanded mixed hybrid methods [27], and positive definite expanded mixed method [9] and so on. From the above literatures on the study of expanded mixed method, we can find that all papers were studied based on Chen’s expanded mixed method [12,13]. For problem (1.1), we introduce the auxiliary variables k ¼ ru; r ¼ aðx; tÞru Rt Rt bðx; s; tÞruds ¼ aðx; tÞk 0 bðx; s; tÞkds and consider the Chen’s expanded mixed method [12,13], to obtain the following 0 expanded mixed weak formulation
8 > < ðaÞ ðut ; v Þ þ ðr r; v Þ ¼ ðf ðx; tÞ; v Þ; 8v 2 W; ðbÞ ðk; wÞ þ ðu; r wÞ ¼ 0; 8w 2 V; > Rt : ðcÞ ðr; zÞ þ ðak; zÞ þ ð 0 bðs; tÞkds; zÞ ¼ 0; 8z 2 X;
ð1:2Þ
2
2
where W ¼ L2 ðXÞ or W ¼ fw 2 L2 ðXÞjwj@ X ¼ 0g; V ¼ Hðdiv;XÞ ¼ fv 2 ðL2 ðXÞÞ jr v 2 L2 ðXÞg; X ¼ L2 ðXÞ ¼ ðL2 ðXÞÞ . In this article, we develop a new expanded mixed finite element method for linear integro-differential equations based on the new mixed weak formulation [28–32]. The proposed method is different from Chen’s expanded mixed scheme (1.2). Compared to expanded mixed scheme (1.2), our scheme mainly covers three advantages: The gradient for our method avoids the use of the classical Hðdiv;XÞ space and just needs to belong to the simple square integrable space L2 ðXÞ; The number of total degrees of freedom is reduced; The error in H1 -norm for the unknown function u can be derived. We will give the proof for the existence and uniqueness of the solution for semi-discrete scheme and a new expanded mixed projection and the proof of its uniqueness. We derived the fully discrete error estimates based on backward Euler scheme. Moreover, we prove the optimal priori error estimates in L2 and H1 -norm for the scalar unknown u and in L2 -norm for its gradient k, and its flux r (the coefficients times the gradient). Finally, we give some numerical experiments to confirm our theoretical results. Throughout this paper, C will denote a generic positive constant which does not depend on the spatial mesh parameter h 2
and temporal discretization parameter Dt. At the same time, we denote the natural inner product in L2 ðXÞ or L2 ðXÞ ¼ ðL2 ðXÞÞ by ð; Þ with the corresponding norm k k. The others notations and definitions of Sobolev spaces as in Ref. [33,34] are used. 2. New expanded mixed formulation Introduce the auxiliary variables k ¼ ru; r ¼ aðx; tÞru lowing first-order system for (1.1)
Rt 0
bðx; s; tÞruds ¼ aðx; tÞk
Rt 0
bðx; s; tÞkds to obtain the fol-
8 > < ut þ r r ¼ f ðx; tÞ; k ru ¼ 0; > R : r þ ak þ 0t bðs; tÞkds ¼ 0;
ð2:1Þ 2
2
The new expanded mixed weak formulation of (2.1) is to find fu; k; rg : ½0; T # H10 ðL2 ðXÞÞ ðL2 ðXÞÞ such that
8 > ðaÞ ðut ; v Þ ðr; rv Þ ¼ ðf ; v Þ; 8v 2 H10 ; > > < 2 ðbÞ ðk; wÞ ðru; wÞ ¼ 0; 8w 2 ðL2 ðXÞÞ ; > R > > : ðcÞ ðr; zÞ þ ðak; zÞ þ t bðs; tÞkds; z ¼ 0; 8z 2 ðL2 ðXÞÞ2 : 0
ð2:2Þ
Then, the semi-discrete mixed finite element scheme for (2.2) is to determine fuh ; kh ; rh g : ½0; T # V h Wh Wh such that
8 > < ðaÞ ðuht ; v h Þ ðrh ; rv h Þ ¼ ðf ; v h Þ; 8v h 2 V h ; ðbÞ ðkh ; wh Þ ðruh ; wh Þ ¼ 0; 8wh 2 Wh ; > Rt : ðcÞ ðrh ; zh Þ þ ðakh ; zh Þ þ ð 0 bðs; tÞkh ds; zh Þ ¼ 0; 8zh 2 Wh ;
ð2:3Þ
where ðV h ; Wh Þ is chosen as the finite element pair P 1 P 20 as follows
V h ¼ fv h 2 C 0 ðXÞ \ H10 jv h 2 P1 ðKÞ; 8K 2 Kh g; 2
Wh ¼ fwh ¼ ðw1h ; w2h Þ 2 ðL2 ðXÞÞ jw1h ; w2h 2 P0 ðKÞ; 8K 2 Kh g:
ð2:4Þ
602
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
We will prove the existence and uniqueness of solution for semi-discrete scheme (2.3). Theorem 2.1. There exists a unique discrete solution to semi-discrete scheme (2.3). 1 2 Proof. Let fwi ðxÞgni¼1 and fuj ðxÞgnj¼1 be bases of V h , Wh , respectively. Let
uh ¼
n1 X
ui ðtÞwi ðxÞ;
kh ¼
i¼1
n2 X
rh ¼
kj ðtÞuj ðxÞ;
j¼1
n2 X
rk ðtÞuk ðxÞ;
k¼1
and
uh ð0Þ ¼
n1 X i wi ðxÞ; u
kh ð0Þ ¼
n2 X kj u ðxÞ; j
i¼1
rh ð0Þ ¼
j¼1
n2 X
r k uk ðxÞ;
k¼1
and substitute these expressions into (2.3) and choose
8 ~ 0 ðtÞ Cr ~ ðtÞ ¼ FðtÞ; ðaÞ Au > > > > < ðbÞ D~kðtÞ Eu ~ ðtÞ ¼ 0; R > ~ ðtÞ þ H~kðtÞ þ 0t J~kðsÞds ¼ 0; ðcÞ Dr > > > : ~ r ~ ð0Þ ¼ u ; kð0Þ ~ ð0Þ ¼ r ; ðdÞu ¼ k;
v h ¼ wm ; wh ¼ ul ; zh ¼ ul to obtain the following equations ð2:5Þ
where
A ¼ ððwi ; wm ÞÞn1 n1 ; C ¼ ððus ; rwm ÞÞn1 n2 ; D ¼ ððuj ; ul ÞÞn
2 n2
H ¼ ððauj ; us ÞÞn
; E ¼ ððrwi ; ul ÞÞn2 n1 ;
2 n2
; J ¼ ððbuj ; us ÞÞn
2 n2
;
~ ðtÞ ¼ ðu1 ðtÞ; u2 ðtÞ; ; un1 ðtÞÞT ; u ~ ¼ ðk1 ðtÞ; k2 ðtÞ; ; kn ðtÞÞT ; kðtÞ 2
r~ ðtÞ ¼ ðr1 ðtÞ; r2 ðtÞ; ; rn2 ðtÞÞT ; ¼ ðu 1 ; u 2 u n1 ÞT ; FðtÞ ¼ ððf ; wm ÞÞn1 1 ; u k ¼ ðk1 ; k2 ; ; kn ÞT ; r 1; r 2; ; r n2 ÞT : ¼ ðr 2 It is easy to see that A and D are invertible matrixes. From (2.5), the initial value problems can be written as follows
8 ~ 0 ðtÞ A1 Cr ~ ðtÞ ¼ A1 FðtÞ; ðaÞ u > > > > < ðbÞ ~kðtÞ ¼ D1 Eu ~ ðtÞ; R 1 ~ > ~ ðtÞ þ D HkðtÞ þ D1 0t J~kðsÞds ¼ 0; > ðcÞ r > > : ~ r ~ ð0Þ ¼ u ; kð0Þ ~ ð0Þ ¼ r : ¼ k; ðdÞ u
ð2:6Þ
Substitute (2.6b) into (2.6c) to get
r~ ðtÞ ¼ D1 HD1 Eu~ ðtÞ D1
Z
t
~ ðsÞds: JD1 Eu
ð2:7Þ
0
Substituting (2.7) into (2.6a), we obtain
~ 0 ðtÞ þ A1 CðD1 HD1 Eu ~ ðtÞ þ D1 u
Z
t
~ ðsÞdsÞ ¼ A1 FðtÞ: JD1 Eu
ð2:8Þ
0
~ ðtÞ, then (2.7) and (2.6b) has a unique soluThus, by the theory of differential equations [35], (2.8) has a unique solution u ~ ~ ðtÞ and kðtÞ, respectively. Equivalently (2.3) has a unique solution. tion r 3. New expanded mixed projection and some Lemmas In order to analyze the convergence of the method, we first introduce the new expanded mixed elliptic projection associated with our equations. ~h ; ~ ~ h Þ : ½0; T ! V h Wh Wh be given by the following mixed relations kh ; r Let ðu
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
8 ~ h ; rv h Þ ¼ 0; 8v h 2 V h ; ðaÞ ðr r > > < ~ ~h Þ; wh Þ ¼ 0; 8wh 2 Wh ; ðbÞ ðk kh ; wh Þ ðrðu u R > > : ðcÞ ðr r ~ h ; zh Þ þ ðaðk ~kh Þ; zh Þ þ 0t bðs; tÞðk ~kh Þds; zh ¼ 0; 8zh 2 Wh :
603
ð3:1Þ
In what follows, we will give some important Lemmas based on new mixed scheme. 2
Lemma 3.1. There exists a linear operator Ph : ðL2 ðXÞÞ ! Wh such that
ðr Ph r; rv h Þ ¼ 0; 8v h 2 V h ;
ð3:2Þ
kr Ph rkL2 ðXÞ 6 ChkrkðH1 ðXÞÞ2 :
ð3:3Þ
and
Lemma 3.2. For the linear operator Ph of Lemma 3.1, we have
ðk Ph k; wh Þ ¼ 0; 8wh 2 Wh ;
ð3:4Þ
kk Ph kkL2 ðXÞ 6 ChkkkðH1 ðXÞÞ2 :
ð3:5Þ
and
Lemma 3.3. There exists a linear operator Ph : H10 ðXÞ ! V h such that
ðrðu Ph uÞ; wh Þ ¼ 0; 8wh 2 Wh ;
ð3:6Þ
and 2
ku Ph ukL2 ðXÞ þ hku Ph ukH1 6 Ch kukH2 ;
ð3:7Þ
2
ð3:8Þ
kut Ph ut kL2 ðXÞ 6 Ch kut kH2 : From Refs. [29,30], we can obtain the proof for Lemma 3.1–3.3. Using the definition of Ph and Ph , we rewrite g; d; q as
g ¼ u u~ h ¼ u Ph u þ Ph u u~ h ¼ gm þ ge ; d ¼ k ~kh ¼ k Ph k þ Ph k ~kh ¼ dm þ de ;
q ¼ r r~ h ¼ r Ph r þ Ph r r~ h ¼ qm þ qe : Since estimates of gm ; dm and qm are known, it is enough to estimate ge ; de and qe . Using Lemma 3.1–3.3, we rewrite (3.1) as
8 ðaÞ ðqe ; rv h Þ ¼ 0; 8v h 2 V h ; > > < ðbÞ ðde ; wh Þ ðrge ; wh Þ ¼ 0; 8wh 2 Wh ; R R > > : ðcÞ ðqe ; zh Þ þ ðade ; zh Þ þ t bðs; tÞde ds; zh ¼ ðadm ; zh Þ t bðs; tÞdm ds; zh ; 8zh 2 Wh : 0 0
ð3:9Þ
We discuss the following approximation properties for system (3.9). Lemma 3.4. There is a constant C independent of h such that
Z t kdk 6 Ch kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds ;
ð3:10Þ
Z t kqk 6 Ch krkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds ;
ð3:11Þ
Z t krgk 6 Ch kukH2 þ kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds :
ð3:12Þ
0
0
0
604
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
Proof. Choose
v h ¼ ge in (3.9a), wh ¼ qe in (3.9b) and zh ¼ de in (3.9c) to obtain
8 ðaÞ ðqe ; rge Þ ¼ 0; > > > > 0; > < ðbÞ ðde ; qe Þ ðrge ; qe Þ ¼ R t ðcÞ ðqe ; de Þ þ ðade ; de Þ þ 0 bðs; tÞde ds; de > > > > > : ¼ ðad ; d Þ R t bðs; tÞd ds; d : m e m e 0
ð3:13Þ
Add the three equations and use Cauchy–Schwarz inequality to get
Z t Z t a0 kde k2 6 ðade ; de Þ ¼ bðs; tÞde ds; de ðadm ; de Þ bðs; tÞdm ds; de 6C
0
Z
t
ðkde k þ kdm kÞds þ kdm k kde k:
0
ð3:14Þ
0
Using (3.14) and Gronwall lemma, we obtain
kde k 6 C
Z
t
kdm kds þ kdm k :
ð3:15Þ
0
Choose zh ¼ qe in (3.9c) and use Cauchy–Schwarz inequality to obtain
Z t kqe k 6 C kde k þ kdm k þ ðkde k þ kdm kÞds :
ð3:16Þ
0
Choose wh ¼ rge in (3.9b) and use (3.15) and Cauchy–Schwarz inequality to obtain
krge k 6 C
Z
t
kdm kds þ kdm k :
ð3:17Þ
0
Combining (3.15)–(3.17) and Lemma 3.1–3.3 and using the triangle inequality, we get the conclusion of Lemma 3.4. Lemma 3.5. There is a constant C independent of h such that 2
kgk 6 Ch kukH2 ;
ð3:18Þ
2
kgt k 6 Ch kut kH2 ;
ð3:19Þ
Z t kgk1 6 Ch kukH2 þ kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds :
ð3:20Þ
0
Proof. To estimate terms kgk; kgt k and kgk1 , we consider the following auxiliary elliptic problem
r ðarvÞ ¼ g; in X; v ¼ 0; on @ X:
ð3:21Þ
Use (3.21), Lemma 3.1–3.3 to obtain
kgk2
~h ; gÞ ¼ ðu u ~ h ; gÞ ¼ ðu Ph u þ Ph u u ~h ; gÞ ¼ ðu Ph u; gÞ þ ðPh u u ~h ; r ðarvÞÞ ¼ ðu Ph u; gÞ þ ðPh u u ~ h Þ; arvÞ ¼ ðu Ph u; gÞ þ ðrðPh u u ~ h Þ; arv Ph ðarvÞÞ ¼ ðu Ph u; gÞ þ ðrðPh u u ~ þðrðPh u uh Þ; Ph ðarvÞÞ
ð3:22Þ
¼ ðu Ph u; gÞ 6 ku Ph ukkgk 2
6 Ch kukH2 kgk: From (3.22), we obtain 2
kgk 6 Ch kukH2 :
ð3:23Þ
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
605
Using the similar method to kgkL2 , we can obtain 2
kgt k 6 Ch kut kH2 :
ð3:24Þ
Combining (3.12) and (3.23), we obtain
kgk1 6 ChðkukH2 þ kkkðH1 ðXÞÞ2 þ
Z
t
kkkðH1 ðXÞÞ2 dsÞ:
0
ð3:25Þ
Using (3.23)–(3.25), we obtain the conclusion of Lemma 3.5.
4. Semi-discrete error estimates In this section, we discuss a priori error estimates. For the error estimates, we decompose the errors as
~h þ u ~ h uh ¼ g þ 1; u uh ¼ u u k kh ¼ k ~kh þ ~kh kh ¼ d þ h;
r rh ¼ r r~ h þ r~ h rh ¼ q þ n: Using (2.2) and (2.3), we can get the error equations
8 ðaÞ ð1t ; v h Þ ðn; rv h Þ ¼ ðgt ; v h Þ; 8v h 2 V h ; > > < ðbÞ ðh; wh Þ ðr1; wh Þ ¼ 0; 8wh 2 Wh ; R > > : ðcÞ ðn; zh Þ þ ðah; zh Þ þ t bhds; zh ¼ 0; 8zh 2 Wh : 0
ð4:1Þ
We will prove the error estimates for semi-discrete scheme. ~ h ð0Þ, then one has the following estimates Theorem 4.1. Assume that uh ð0Þ ¼ u 2
ku uh k þ kut uht k 6 Ch
Z t kukH2 þ kut kH2 þ kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds ;
ð4:2Þ
0
Z t ku uh k1 6 Ch kukH2 þ kut kH2 þ kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds ;
ð4:3Þ
Z t kk kh k 6 Ch kukH2 þ kut kH2 þ kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds ;
ð4:4Þ
Z t kr rh k 6 Ch kukH2 þ kut kH2 þ krkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 þ kkkðH1 ðXÞÞ2 ds :
ð4:5Þ
0
0
0
Proof. Choose
v h ¼ 1 in (4.1a), wh ¼ n in (4.1b) and zh ¼ h in (4.1c) to obtain
8 > ðaÞ 12 dtd k1k2 ðn; r1Þ ¼ ðgt ; 1Þ; > < ðbÞ ðh; nÞ ðr1; nÞ ¼ 0; > > : ðcÞ ðn; hÞ þ ka12 hk2 þ R t bhds; h ¼ 0:
ð4:6Þ
0
Adding the above three equations, and using Cauchy–Schwarz inequality, Young inequality, we have
1 d 1 k1k2 þ ka2 hk2 ¼ ðgt ; 1Þ 2 dt
Z 0
t
Z t bhds; h 6 C kgt k2 þ k1k2 þ khk2 þ khk2 ds :
ð4:7Þ
0
Integrate with respect to time from 0 to t to obtain
k1k2 þ
Z
t
khk2 ds 6 C
0
Z t Z s kgt k2 þ k1k2 þ khk2 þ khðsÞk2 ds ds: 0
ð4:8Þ
0
Using Gronwall Lemma, we obtain
k1k2 þ
Z
t 0
khk2 ds 6 C
Z 0
t
kgt k2 ds:
ð4:9Þ
606
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
Taking wh ¼ r1 in (4.1b), we have
kr1k 6 khk:
ð4:10Þ
Differentiating (4.1b) and taking wh ¼ n, we obtain
ðht ; nÞ ðr1t ; nÞ ¼ 0:
v h ¼ 1t
Choose
ð4:11Þ
in (4.1a) and zh ¼ ht in (4.1c) to obtain
8 < ðaÞ k1t k2 ðn; r1t Þ ¼ ðgt ; 1t Þ; : ðcÞ ðn; ht Þ þ 1 2
1 d ka2 hk2 dt
¼ 12 ðat h; hÞ þ
R t d dt
0
R t bhds; h dtd 0 bhds; h :
ð4:12Þ
Adding (4.11), (4.12a) and (4.12c), and using Cauchy–Schwarz inequality, Young inequality, we have
k1t k2 þ
Z t Z t 1 d 1 2 1 d ka2 hk ¼ ðgt ; 1t Þ þ ðat h; hÞ þ bt hds; h þ bðx; t; tÞh bhds; h 2 dt 2 dt 0 0 Z t Z t 1 d khk2 ds þ k1t k2 bhds; h : 6 C kgt k2 þ khk2 þ 2 dt 0 0
ð4:13Þ
Integrate with respect to time from 0 to t and use Gronwall lemma to obtain
Z
t
0
Z
k1t k2 ds þ khk2 6 C
t
0
kgt k2 ds:
ð4:14Þ
Choosing zh ¼ n in (4.1c) and using (4.9) and (4.14), we have
1 knk2 ¼ ðah; nÞ 2
Z
t
0
Z t Z t 1 bhds; n knk2 6 C khk2 þ khk2 ds 6 C kgt k2 ds: 2 0 0
ð4:15Þ
Combining (3.11), (3.12), (4.9), (4.10), (4.14), (4.15) and the triangle inequality, we obtain the error estimates for Theorem 4.1.
5. Fully discrete error estimates In this section, we get the error estimates of fully discrete schemes based on backward Euler method. For the fully discrete procedure, let 0 ¼ t 0 < t1 < t2 < < t M ¼ T be a given partition of the time interval ½0; T with step length Dt ¼ T=M and nodes tn ¼ nDt, for some positive integer M. For a smooth function / on ½0; T, define /n ¼ /ðtn Þ and @t /n ¼ ð/n /n1 Þ=Dt. For approximating the integrals, we use the composite left rectangle rule
Dt
Z n1 X /j
tn
/ðsÞds:
0
j¼0
Note that / 2 C 1 ½0; T, the quadrature error satisfies
Z tn Z tn X n1 j /ðsÞds 6 C Dt j/s ðsÞjds: Dt / j¼0 0 0 The Eq. (2.2) has the following equivalent formulation
8 n > ðaÞ ð@ t un ; v Þ ðrn ; rv Þ ¼ ðf þ Rn1 ; v Þ; 8v 2 H10 ; > > > > < ðbÞ ðkn ; wÞ ðrun ; wÞ ¼ 0; 8w 2 ðL2 ðXÞÞ2 ; ! n1 > X > 2 j > n n n > ðcÞ ð ¼ ðRn2 ; zÞ; 8z 2 ðL2 ðXÞÞ ; r ; zÞ þ ða k ; zÞ þ D t bðx; t ; t Þk ; z > n j : j¼0
where
Rn1 ¼ @ t un ut ðt n Þ ¼
Rn2 ¼ Dt
n1 X j¼0
1 Dt
Z
bðx; t n ; tj Þk j
tn
ðtn1 sÞutt ds;
t n1
Z
tn
bðx; t n ; sÞkðsÞds:
0
Now we formulate a fully discrete procedure: Find ðunh ; knh ; rnh Þ 2 V h Wh Wh ; ðn ¼ 0; 1; ; MÞ such that
ð5:1Þ
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
8 n ðaÞ ð@ t unh ; v h Þ ðrnh ; rv h Þ ¼ ðf ; v h Þ; 8v h 2 V h ; > > > > < ðbÞ ðknh ; wh Þ ðrunh ; wh Þ ¼ 0; 8wh 2 Wh ; ! n1 X > > j n n n > > : ðcÞ ðrh ; zh Þ þ ða kh ; zh þ Dt bðx; t n ; tj Þkh ; zh ¼ 0; 8zh 2 Wh :
607
ð5:2Þ
j¼0
For the fully discrete error estimates, we now split the errors
~nh þ u ~ nh unh ¼ gn þ 1n ; uðtn Þ unh ¼ uðtn Þ u kðtn Þ knh ¼ kðtn Þ ~knh þ ~knh knh ¼ dn þ hn ;
rðtn Þ rnh ¼ rðtn Þ r~ nh þ r~ nh rnh ¼ qn þ nn : We will prove the theorem for the fully discrete error estimates. ~ h ð0Þ, then there exists a positive constant C independent of h and Dt such that Theorem 5.1. Assume that u0h ¼ u
kuðt J Þ uJh k 6 Ch
2
Z kukL1 ðH2 Þ þ kut kL1 ðH2 Þ þ kkkL1 ððH1 Þ2 Þ þ C Dt
tJ
t0
Z kuðt J Þ uJh k1 6 Ch kukL1 ðH2 Þ þ kut kL1 ðH2 Þ þ kkkL1 ððH1 Þ2 Þ þ C Dt Z
ð5:3Þ
kkkðH1 Þ2 ds;
ð5:4Þ
tJ
t0
kkðt J Þ kJh k 6 ChðkukL1 ðH2 Þ þ kut kL1 ðH2 Þ þ kkkL1 ððH1 Þ2 Þ Þ þ C Dt
kkkðH1 Þ2 ds;
tJ
t0
kkkðH1 Þ2 ds;
ð5:5Þ
Z krðt J Þ rJh k 6 Ch kukL1 ðH2 Þ þ kut kL1 ðH2 Þ þ kkkL1 ððH1 Þ2 Þ þ krkL1 ððH1 Þ2 Þ þ C Dt
tJ
t0
kkkðH1 Þ2 ds:
ð5:6Þ
Proof. Using (3.1), (5.1) and (5.2) at t ¼ tn , we get the error equations
8 ðaÞ ð@ t 1n ; v h Þ ðnn ; rv h Þ ¼ ð@ t gn ; v h Þ þ ðRn1 ; v h Þ; 8v h 2 V h ; > > > > < ðbÞ ðhn ; wh Þ ðr1n ; wh Þ ¼ 0; 8wh 2 Wh ; ! n1 X > > n n j n n > > : ðcÞ ðn ; zh Þ þ ða h ; zh Þ þ Dt bðx; tn ; t j Þh ; zh ¼ ðR2 ; zh Þ; 8zh 2 Wh :
ð5:7Þ
j¼0
Choose
v h ¼ 1n
in (5.7a), wh ¼ nn in (5.7b) and zh ¼ hn in (5.7c) to obtain
8 > ðaÞ 21Dt ðk1n k2 k1n1 k2 þ k1n 1n1 k2 Þ > > > n n > n n n n > > < ðn ; r1 Þ ¼ ð@ t g ; 1 Þ þ ðR1 ; 1 Þ; n n n ðbÞ ðh ; n Þ ðr1n ; n Þ ¼ 0; ! > > n1 > X > 1 n n j n > n 2 n 2 > ðcÞ ðn ¼ ðRn2 ; hn Þ: ; h Þ þ kða Þ h k þ D t bðx; t ; t Þh ; h > n j :
ð5:8Þ
j¼0
Adding the above three equations, we obtain n1 X 1 1 ðk1n k2 k1n1 k2 Þ þ kðan Þ2 hn k2 ¼ ð@ t gn ; 1n Þ Dt bðx; tn ; t j Þh j ; hn 2Dt j¼0 n 2
6 C k@ t g k þ
kRn1 k2
n 2
þ k1 k þ Dt
!
n1 X j¼0
þ ðRn1 ; 1n Þ þ ðRn2 ; hn Þ ! j 2
kh k
þ
a0 n 2 kh k : 2
ð5:9Þ
Multiplying by 2Dt and summing (5.9) from n ¼ 1 to J, the resulting equation becomes
! J J J1 n1 X X X a0 X n 2 n 2 n 2 j 2 0 2 n 2 n 2 ð1 C DtÞk1 k þ Dt kh k 6 k1 k þ C Dt k@ t g k þ kR1 k þ kR2 k þ C Dt k1 k þ Dt kh k : 2 n¼1 n¼1 n¼1 j¼0 J 2
ð5:10Þ
608
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
Choose Dt0 in such a way that for 0 < Dt 6 Dt 0 ; ð1 C DtÞ > 0. Then an application of Gronwall’s lemma to obtain
k1J k2 þ a0 Dt
J J X X khn k2 6 C Dt k@ t gn k2 þ kRn1 k2 þ kRn2 k2 : n¼1
ð5:11Þ
n¼1
Note that
ðaÞ kRn1 k2 6 C Dt
R tn
kutt k2 ds; R t ðbÞ kRn2 k2 6 CðDtÞ2 0n kkk2 þ kkt k2 ds; R tn kgt ðsÞk2 ds: ðcÞ k@ t gn k2 6 C D1t tn1 tn1
ð5:12Þ
Substitute (5.12) into (5.11) to get
k1J k2 þ a0 Dt
Z J X khn k2 6 C
tJ t0
n¼1
kgt ðsÞk2 ds þ CðDtÞ2
Z
tJ
kutt k2 þ kkk2 þ kkt k2 ds:
ð5:13Þ
t0
Taking wh ¼ r1n in (5.7b), we have
kr1n k 6 khn k:
ð5:14Þ
From (5.7b), we get
ð@ t hn ; wh Þ ðr@ t 1n ; wh Þ ¼ 0: Choose wh ¼ nn in (5.15),
ð5:15Þ
v h ¼ @ t 1n in (5.7a) and zh ¼ @ t hn in (5.7c) to obtain
8 ðaÞ k@ t 1n k2 ðnn ; r@ t 1n Þ ¼ ð@ t gn ; @ t 1n Þ þ ðRn1 ; @ t 1n Þ; > > > > < ðbÞ ð@ t hn ; nn Þ ðr@ t 1n ; nn Þ ¼ 0; ! n1 X > > n n n n n j n n > > : ðcÞ ðn ; @ t h Þ þ ðaðx; t n Þ@ t h ; h Þ þ Dt bðx; t n ; tj Þh ; @ t h ¼ ðR2 ; @ t h Þ:
ð5:16Þ
j¼0
Adding the three equations, we obtain
1 an1 an 1 1 n 12 n 2 kða Þ h k kðan1 Þ2 hn1 k2 þ hn1 ; hn1 2 Dt 2 Dt ! n1 X 6 ð@ t gn ; @ t 1n Þ þ ðRn1 ; @ t 1n Þ Dt bðx; t n ; t j Þh j ; @ t hn þ ðRn2 ; @ t hn Þ:
k@ t 1n k2 þ
ð5:17Þ
j¼0
Note that
ðaÞðRn2 ; @ t hn Þ ¼
ðbÞ Dt
ðhn ;Rn2 Þðhn1 ;Rn1 2 Þ Dt
n1 X bðx; tn ; t j Þh j ; @ t hn
!
ð@ t Rn2 ; hn1 Þ; P Dt
n1 j¼0
Dt
!
n2 X
bðx;t n1 ;t j Þh j ;hn1
j¼0
¼
j¼0
bðx;t n ;t j Þh j ;hn
Dt
P bðx;t n ;t j Þbðx;t n1 ;t j Þ j Dt n1 h ; hn1 j¼0 Dt
ð5:18Þ
bðx; t n1 ; tn1 Þhn1 ; hn1 :
Substitute (5.18) into (5.17) to obtain
k@ t 1n k2 þ
1 an1 an 1 1 n 12 n 2 kða Þ h k kðan1 Þ2 hn1 k2 þ hn1 ; hn1 2 Dt 2 Dt n
n
6 ð@ t g ; @ t 1 Þ þ þ
ðRn1 ; @ t n Þ
1 þ Dt
n1 X bðx; t n ; t j Þ bðx; tn1 ; t j Þ j¼0
Dt
Pn1 j¼0
Dt
P j n1 bðx; t n ; tj Þh j ; hn Dt n2 j¼0 bðx; t n1 ; t j Þh ; h
Dt ðhn ; Rn2 Þ ðhn1 ; Rn1 2 Þ ð@ t Rn2 ; hn1 Þ: Dt
j
h ;h
!
n1
þ ðbðx; t n1 ; t n1 Þhn1 ; hn1 Þ ð5:19Þ
609
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
Multiplying by 2Dt and summing (5.19) from n ¼ 1 to J, the resulting equation becomes
Dt
J J J X X X 1 k@ t 1n k2 þ ða0 C DtÞkhJ k2 6 kða0 Þ2 h0 k2 þ C Dt k@ t gn k2 þ kRn1 k2 þ k@ t Rn2 k2 þ C Dt khn1 k2 n¼0
n¼1
n¼1
J X n1
X
bðx; tn ; t j Þ bðx; t n1 ; t j Þ j 2 n 2
þ ðDtÞ3 h
þ kR2 k Dt n¼1 j¼0
1 2
6 kða0 Þ h0 k2 þ C Dt
J J X X k@ t gn k2 þ kRn1 k2 þ k@ t Rn2 k2 þ C Dt khn1 k2 þ kRn2 k2 : n¼1
ð5:20Þ
n¼1
Using (5.12), we obtain
ðaÞ Dt
J X Rt kR1j k2 6 CðDtÞ2 t0J kutt k2 ds; j¼1
ð5:21Þ
J X Rt ðbÞ Dt k@ t g j k2 6 C t0J kgt k2 ds: j¼1
Now, we estimate the term Dt
@ t Rn2 ¼ Dt
PJ
j 2 j¼1 k@ t R2 k .
Z n2 X ½bðtn ; t j Þ bðtn1 ; tj Þk j þ Dtbðt n ; tn1 Þkn1
bðt n ; sÞkðsÞds
t n1
j¼0
¼ ðDtÞ2
tn
n2 X bt ðt n0 ; t j Þk j þ Dt½bðtn ; t n1 Þkn1 bðtn ; t n1 Þkðtn1 Þ j¼0
n2 X ¼ ðDtÞ2 bt ðt n0 ; t j Þk j þ Dt½bðtn ; t n1 Þkt ðt n2 Þðt n1 t n1 Þ þ bt ðt n ; tn3 Þðtn1 t n1 Þkðt n1 Þ; j¼0
where tn0 ; t n1 2 ðtn1 ; t n Þ; t n2 ; tn3 2 ðt n1 ; t n1 Þ: From (5.22), we obtain
j@ t Rn2 j 6 Dt max jbt ðt n0 ; t j Þk j j þ ðDtÞ2 ½jbðt n ; tn1 Þkt ðt n2 Þj þ jbt ðt n ; t n3 Þkðt n1 Þj:
ð5:23Þ
06j6n2
Using (5.23), we obtain
Dt
J X k@ t R2j k2 6 ðDtÞ2 kkk2L1 ðL2 Þ þ kkt k2L1 ðL2 Þ :
ð5:24Þ
j¼1
Substitute (5.21), (5.24) and (5.12) into (5.20) and use Gronwall lemma to get
2Dt
Z J X k@ t 1n k2 þ khJ k2 6 C
tJ
t0
n¼0
kgt ðsÞk2 ds þ CðDtÞ2 kutt k2L2 ðL2 Þ þ kkk2L1 ðL2 Þ þ kkt k2L1 ðL2 Þ þ kkk2L2 ðL2 Þ þ kkt k2L2 ðL2 Þ :
Taking zh ¼ n in (5.7c), we get
knn k2 ¼ ðan hn ; nn Þ ðDt
n1 n1 X X bðx; t n ; tj Þh j ; nn Þ þ ðRn2 ; nn Þ 6 C khn k2 þ ðDtÞ2 kh j k2 j¼0
j¼0
!
1 þ knn k2 : 2
ð5:25Þ
ð5:26Þ
Substitute (5.25) into (5.26) to get
knJ k2 6 C
Z
tJ
t0
kgt ðsÞk2 ds þ CðDtÞ2 kutt k2L2 ðL2 Þ þ kkk2L1 ðL2 Þ þ kkt k2L1 ðL2 Þ þ kkk2L2 ðL2 Þ þ kkt k2L2 ðL2 Þ :
ð5:27Þ
Combining (3.11), (3.12), (5.13), (5.14), (5.25), (5.27) and the triangle inequality, we complete the proof. 6. Numerical experiments For confirming our theoretical analysis, we need to consider a numerical example based on problem (1.1) with initial boundary conditions. In (1.1), the spatial domain and the temporal interval are X ¼ ð0; 1Þ ð0; 1Þ and J ¼ ð0; 1, respectively. The known source term f ðx; tÞ is decided by the chosen coefficients aðx; tÞ ¼ 1 þ x21 þ x22 þ t2 ; bðx; t; sÞ ¼ 1 þ x21 þ x22 þ t 2 þ s and the exact solution
uðx; tÞ ¼ et sinðpx1 Þ sinð2px2 Þ;
610
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
the exact gradient function
k ¼ ru ¼ ðpet cosðpx1 Þ sinð2px2 Þ; 2pet sinðpx1 Þ cosð2px2 ÞÞ; and the exact flux function
r ¼ aru
Z
t
bðx; t; sÞruds ¼ ðdðx; tÞ cosðpx1 Þ sinð2px2 Þ; 2dðx; tÞ sinðpx1 Þ cosð2px2 ÞÞ;
0
where x ¼ ðx1 ; x2 Þ and dðx; tÞ ¼ pð2 þ x21 þ x22 þ t2 tet et Þ. We first divide the spatial domain X ¼ ½0; 1 ½0; 1 into uniform triangulations with spatial mesh size h, use the backward Euler approximation with uniform time step length Dt, choose the piecewise linear space V h with index k ¼ 1 and the piecewise constant space Wh with index k ¼ 0, then obtain the fully discrete procedure. Based on the above consideration, we will obtain some calculated data on orders of convergence and a priori errors. pffiffi pffiffi pffiffi pffiffi With a fixed time step Dt ¼ 1=100 and some changed space meshes length h ¼ 42 ; 82 ; 162 ; 322, we obtain some calculated results of a priori error estimates in L2 and H1 -norms for the scalar unknown u in Table 1. From Table 1, we find that the convergence rate in L2 -norm for u is close to 2, while the one in H1 -norm is close to 1. So we claim that these computed convergence results, which are in correspondence with our theoretical results, are optimal. With the same choices of space–time 2
meshes to the case for u, we also calculate a priori error results and orders of convergence in L2 ðXÞ ¼ ðL2 ðXÞÞ -norm for the gradient k and its flux r in Table 2. From the calculated results, we see that the convergence rate, which is also in correspondence with the results of theorem, is close to 1.
Table 1 Error estimates and convergence rate in L2 and H1 -norms for u with Dt ¼ 1=100. h pffiffi 2 4 pffiffi 2 8 pffiffi 2 16 pffiffi 2 32
ku uh kL1 ðL2 ðXÞÞ
Rate
ku uh kL1 ðH1 ðXÞÞ
1.7778e001
Rate
1.2816e+000
4.7102e002
1.9162
6.2158e001
1.0439
1.1941e002
1.9799
3.0833e001
1.0115
3.0125e003
1.9869
1.5384e001
1.0030
Rate
kr rh kL1 ðL2 ðXÞÞ
Rate
6.2068e001
1.0366
1.8693e+000
1.0416
3.0822e001
1.0099
9.2991e001
1.0073
1.5383e001
1.0026
4.6435e001
1.0019
Table 2 Error estimates and convergence rate in L2 -norm for gradient k and its flux r with Dt ¼ 1=100. h pffiffi 2 4 pffiffi 2 8 pffiffi 2 16 pffiffi 2 32
kk kh kL1 ðL2 ðXÞÞ 1.2733e+000
3.8480e+000
0.4 0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 1
1 0.8
0.6
0.5 0.4
0.2
0
Fig. 1. The exact solution u.
0
611
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
0.4 0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 1
1 0.8
0.6
0.5 0.4
0.2
0
0
Fig. 2. The numerical solution uh .
2
4
1
2
0
0 1
−1 −2 1
1
−2 −4 1
0.5 0.5
0.5 0.5
0 0
0 0
Fig. 3. The exact gradient function k ¼ ðk1 ; k2 Þ.
2
4
1
2
0
0 1
−1 −2 1
0.5 0.5
1
−2 −4 1
0.5 0.5
0 0
Fig. 4. The numerical gradient function kh ¼ ðk1h ; k2h Þ.
0 0
612
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
20
20
10
0
0 1
−10 −20 1
−20
1
−40 1
0.5
0.5
0.5
0.5 0 0
0 0
Fig. 5. The exact flux function
20
r ¼ ðr1 ; r2 Þ.
20
10
0
0 1
−10 −20 1
0.5
−20
1
−40 1
0.5
0.5 0.5
0 0
Fig. 6. The numerical flux function
0 0
rh ¼ ðr1h ; r2h Þ.
Now we will make some comparisons of figures between the exact solution and the numerical solution. We show the exact unknown function u in Fig. 1 at time t ¼ 1 with the chosen time step length Dt ¼ 1=100 and the spatial mesh size pffiffi h ¼ 162. Under the same conditions, we also show the numerical solution uh in Fig. 2. Form Figs. 1 and 2, it is not hard to find that the numerical function uh can primely approximate the exact function u. The similar comparisons are made between the exact gradient solution k and the numerical gradient solution kh in Figs. 3 and 4. Finally, we also give the comparison of figures between the exact solution r of flux function and the numerical solution rh of flux function in Figs. 5,6. Based on the discussion of the calculated results in Tables 1, 2 and Figs. 1–6, we can confirm our theoretical analysis.
7. Concluding remarks In this article, we proposed and analyzed a new expanded mixed finite element scheme, whose gradient belongs to the square integrable space instead of the classical Hðdiv;XÞ, for linear parabolic integro-differential equations. We introduced the new expanded mixed projection and proved some important Lemmas based on new expanded mixed projection. We derived the optimal priori error estimates in L2 for the scalar unknown u and the suboptimal error estimates in L2 -norm for its gradient k, and its flux r (the coefficients times the negative gradient). Moreover, we can obtained the H1 -norm for the scalar unknown u, which cannot be obtained by the previous expanded mixed finite element methods. In the near future, the proposed expanded mixed scheme will be applied to other equations, such as nonlinear parabolic equation, hyperbolic wave equation and so on.
Y. Liu et al. / Applied Mathematics and Computation 259 (2015) 600–613
613
Acknowledgment Authors thank the anonymous reviewers and editor for their valuable comments and suggestions, which greatly improve this work. 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. 2015.02.081. References [1] T. Zhang, Finite Element Methods For Partial Differential-Integral Equations, Chinese Science Press, Beijing, 2009. [2] A.K. Pani, T.E. Peterson, Finite element methods with numerical quadrature for parabolic integro-differential equations, SIAM J. Numer. Anal. 33 (3) (1996) 1084–1105. [3] A.K. Pani, R.K. Sinha, Error estimates for semidiscrete Galerkin approximation to a time dependent parabolic integro-differential equation with nonsmooth data, Calcolo 37 (4) (2000) 181–205. [4] Y. Lin, V. Thomée, L.B. Wahlbin, Ritz–Volterra projections to finite-element spaces and applications to integro-differential and related equations, SIAM J. Numer. Anal. 28 (4) (1991) 1047–1070. [5] R.K. Sinha, R.E. Ewing, R.D. Lazarov, Mixed finite element approximations of parabolic integro-differential equations with nonsmooth initial data, SIAM J. Numer. Anal. 47 (5) (2009) 3269–3292. [6] Z. Jiang, L1 ðL2 Þ and L1 ðL1 Þ error estimates for mixed methods for integro-differential equations of parabolic type, Math. Model. Numer. Anal. 33 (3) (1999) 531–546. [7] H. Guo, H.X. Rui, Least-squares Galerkin procedures for parabolic integro-differential equations, Appl. Math. Comput. 150 (3) (2004) 749–762. [8] D.Y. Shi, H.H. Wang, An H1 -Galerkin nonconforming mixed finite element method for integro-differential equation of parabolic type, J. Math. Res. Exposition 29 (5) (2009) 871–881. [9] Y. Liu, H. Li, J.F. Wang, W. Gao, A new positive definite expanded mixed finite element method for parabolic integrodifferential equations, J. Appl. Math. 2012, Article ID 391372, 24pages. [10] H. Guo, J.S. Zhang, H.F. Fu, Two splitting positive definite mixed finite element methods for parabolic integro-differential equations, Appl. Math. Comput. 218 (22) (2012) 11255–11268. [11] H.R. Li, Q. Li, Finite volume element methods for nonlinear parabolic integro-differential problems, J. KSIAM 7 (2) (2003) 35–49. [12] Z.X. Chen, Expanded mixed finite element methods for linear second order elliptic problems I, IMA preprint series ] 1219, Institute for Mathematics and its Applications, University of Minnesota, Minneapolis, MN, 1994. [13] Z.X. Chen, Expanded mixed element methods for linear second-order elliptic problems(I), RAIRO Model. Math. Anal. Numer. 32 (1998) 479–499. [14] T. Arbogast, M.F. Wheeler, I. Yotov, Mixed finite elements for elliptic problems with tensor coefficients as cellcentered finite differences, SIAM J. Numer. Anal. 34 (1997) 828–852. [15] Z.X. Chen, Expanded mixed element methods for quasilinear second-order elliptic problems(II), RAIRO Model. Math. Anal. Numer. 32 (1998). 501–420. [16] C.S. Woodward, C.N. Dawson, Analysis of expanded mixed finite element methods for a non-linear parabolic equation modelling flow into variably saturated porous media, SIAM J. Numer. Anal. 37 (2000) 701–724. [17] Y.P. Chen, Y.Q. Huang, D.H. Yu, A two-grid method for expanded mixed finite-element solution of semilinear reaction-diffusion equations, Int. J. Numer. Methods Eng. 57 (2) (2003) 193–209. [18] Y.P. Chen, H.W. Liu, S. Liu, Analysis of two-grid methods for reaction–diffusion equations by expanded mixed finite element methods, Int. J. Numer. Methods Eng. 69 (2) (2007) 408–422. [19] Y.P. Chen, P. Luan, Z.L. Lu, Analysis of two-grid methods for nonlinear parabolic equations by expanded mixed finite element methods, Adv. Appl. Math. Mech. 1 (6) (2009) 830–844. [20] W. Liu, H.X. Rui, H. Guo, A two-grid method with expanded mixed element for nonlinear reaction–diffusion equations, Acta Math. Appl. Sin. 27 (3) (2011) 495–502. [21] H.L. Song, Y.R. Yuan, The expanded upwind-mixed multi-step method for the miscible displacement problem in three dimensions, Appl. Math. Comput. 195 (2008) 100–109. [22] L. Guo, H.Z. Chen, An expanded characteristic-mixed finite element method for a convection-dominated transport problem, J. Comput. Math. 23 (5) (2005) 479–490. [23] H.Z. Chen, H. Wang, An optimal-order error estimate on an H1 -Galerkin mixed method for a nonlinear parabolic equation in porous medium flow, Numer. Methods Partial Differ. Equ. 26 (2010) 188–205. [24] H.T. Che, Y.J. Wang, Z.J. Zhou, An optimal error estimates of H1 -Galerkin expanded mixed finite element methods for nonlinear viscoelasticity-type equation, Math. Prob. Eng. 2011, Article ID 570980, 18 pages doi:10.1155/2011/570980. [25] Y. Liu, Analysis and numerical simulation of nonstandard mixed element methods (Ph.D. thesis), Inner Mongolia University, Hohhot, China, 2011. [26] H.X. Rui, T.C. Lu, An expanded mixed covolume method for elliptic problems, Numer. Methods Partial Differ. Equ. 21 (1) (2005) 8–23. [27] D. Kim, E.J. Park, A posteriori error estimator for expanded mixed hybrid methods, Numer. Methods Partial Differ. Equ. 23 (2007) 330–349. [28] Y. Liu, H. Li, et al., A new expanded mixed finite element method for elliptic equations, submitted for publication. [29] S.C. Chen, H.R. Chen, New mixed element schemes for second order elliptic problem, Math. Numer. Sin. 32 (2) (2010) 213–218. [30] F. Shi, J.P. Yu, K.T. Li, A new stabilized mixed finite-element method for Poisson equation based on two local Gauss integrations for linear element pair, Int. J. Comput. Math. 88 (11) (2011) 2293–2305. [31] D.Y. Shi, Y.D. Zhang, High accuracy analysis of a new nonconforming mixed finite element scheme for Sobolev equations, Appl. Math. Comput. 218 (7) (2011) 3176–3186. [32] Z.F. Weng, X.L. Feng, P.Z. Huang, A new mixed finite element method based on the Crank–Nicolson scheme for the parabolic problems, Appl. Math. Model. (2012), http://dx.doi.org/10.1016/j.apm.2011.12.044. [33] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [34] Z.D. Luo, Mixed Finite Element Methods and Applications, Chinese Science Press, Beijing, China, 2006. [35] H. Brunner, P.J. Houwen, The Numerical Solution of Volterra Equations, North-Holland Publishing Co., Amsterdam, 1986.