# A New Positive Definite Expanded Mixed Finite Element Method for Parabolic Integrodifferential Equations

Authors: Yang Liu; Hong Li; Jinfeng Wang; Wei Gao
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

A new positive definite expanded mixed finite element method is proposed for parabolic partial integrodifferential equations. Compared to expanded mixed scheme, the new expanded mixed element system is symmetric positive definite and both the gradient equation and the flux equation are separated from its scalar unknown equation. The existence and uniqueness for semidiscrete scheme are proved and error estimates are derived for both semidiscrete and fully discrete schemes. Finally, some numerical results are provided to confirm our theoretical analysis.


---

## Body

## 1. Introduction

In this paper, we consider the following initial-boundary value problem of parabolic partial integrodifferential equations:(1.1)ut-∇⋅(a(x,t)∇u+b(x,t)∫0t‍∇uds)=f(x,t),(x,t)∈Ω×J,u(x,t)=0,(x,t)∈∂Ω×J-,u(x,0)=u0(x),x∈Ω,
where Ω is a bounded convex polygonal domain in ℝd,(d=1,2,3) with a smooth boundary ∂Ω, J=(0,T] is the time interval with 0<T<∞. The coefficients a=a(x,t), b=b(x,t) are two functions, which satisfy the property that there exist some positive constants amin, amax, bmin, and bmax such that 0<amin≤a(x,t)≤amax and 0<bmin≤b(x,t)≤bmax.Parabolic integrodifferential 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, and nuclear reactor dynamics. In recent years, a lot of researchers have studied the numerical methods for parabolic integrodifferential equations, such as finite element methods [1–5], mixed finite element methods [6–9], and finite volume element method [10] and so forth.In 1994, Chen [11, 12] proposed a new mixed method, which is called a expanded mixed finite element method and proved some mathematical theories for second-order linear elliptic equation. 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 gradient). In 1997, Arbogast et al. [13] 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 [14] and fourth-order elliptic problems [15]. With the development of the expanded mixed finite element method, the method was applied to many evolution equations. In [16], some error estimates of the expanded mixed element for a kind of parabolic equation were given. Woodward and Dawson [17] studied the expanded mixed finite element method for nonlinear parabolic equation. Wu and Chen et al. [18–22] studied the two-grid methods for expanded mixed finite-element solution of semilinear reaction-diffusion equations. Song and Yuan [23] proposed the expanded upwind-mixed multistep method for the miscible displacement problem in three dimensions. Guo and Chen [24] developed and analysed an expanded characteristic-mixed finite element method for a convection-dominated transport problem. In 2010, Chen and Wang [25] proposed an H1-Galerkin expanded mixed method for a nonlinear parabolic equation in porous medium flow, and Liu and Li [26] studied the H1-Galerkin expanded mixed method for pseudo-hyperbolic equation. Liu [27], studied the H1-Galerkin expanded mixed method for RLW-Burgers equation and proved semidiscrete and fully discrete optimal error estimates. Jiang and Li [28] studied an expanded mixed semidiscrete scheme for the problem of purely longitudinal motion of a homogeneous bar. In [29, 30], the expanded mixed covolume method was studied for the linear integrodifferential equation of parabolic type and elliptic problems, respectively. In [31], a posteriori error estimator for expanded mixed hybrid methods was studied and analysed.In 2001, Yang [32] proposed a new mixed finite element method called the splitting positive definite mixed finite element procedure to treat the pressure equation of parabolic type in a nonlinear parabolic system describing a model for compressible flow displacement in a porous medium. Compared to standard mixed methods whose numerical solutions have been quite difficult because of losing positive definite properties, the proposed one does not lead to some saddle point problems. From then on, the method was applied to the hyperbolic equations [33] and pseudo-hyperbolic equations [34].In this paper, our purpose is to propose and analyse a new expanded mixed method based on the positive definite system [32–34] for parabolic integrodifferential equations. Compared to expanded mixed methods, the proposed mixed element system is symmetric positive definite and avoids some saddle point problems. What is more, both the gradient equation and the flux equation are separated from its scalar unknown equation. The existence and uniqueness for semidiscrete scheme are proved and error estimates are derived for both semidiscrete and fully discrete schemes.The layout of the paper is as follows. In Section2, the positive definite expanded mixed weak formulation and semidiscrete mixed scheme are formulated, and the proof of the existence and uniqueness of the discrete solutions is given. Error estimates are derived for both semidiscrete and fully discrete schemes for problems, respectively, in Sections 3 and 4. In Section 5, some numerical results are provided to illustrate the efficiency of our method. Finally in Section 6, we will give some concluding remarks about the positive definite expanded mixed finite element method.Throughout this paper,C will denote a generic positive constant which does not depend on the spatial mesh parameters hu,hσ and time-discretization parameter δ and may be different at their occurrences. Usual definitions, notations, and norms of the Sobolev spaces as in [35–37] are used. We denote the natural inner product in L2(Ω) or [L2(Ω)]d by (·,·) with norm ∥·∥L2(Ω) or ∥·∥L2(Ω) and introduce the function space W=H(div;Ω)={ω∈[L2(Ω)]d;∇·ω∈L2(Ω)}.


## 2. A New Expanded Mixed Variational Formulation

Introducing the auxiliary variables:(2.1)λ=∇u,σ=a(x,t)∇u+b(x,t)∫0t‍∇uds=aλ+b∫0t‍λds,
Then we obtain the equivalent system of parabolic partial integrodifferential equations for the problem (1.1):
(2.2)(a)ut-∇⋅σ=f(x,t),(x,t)∈Ω×J,(b)λ-∇u=0,(x,t)∈Ω×J,(c)σ-aλ-b∫0t‍λds=0,(x,t)∈Ω×J,
with the initial values λ(x,0)=∇u0(x),σ(x,0)=a∇u0(x), and u(x,0)=u0(x).Then, the following expanded mixed weak formulation of (2.2) can be given by:
(2.3)(a)(ut,v)-(∇⋅σ,v)=(f(x,t),v),∀v∈L2(Ω),(b)(λ,w)+(u,∇⋅w)=0,∀w∈W,(c)(σ,z)-(aλ,z)-(b∫0t‍λds,z)=0,∀z∈W.From (2.3)(b) we derive
(2.4)(λt,w)+(ut,∇⋅w)=0.
Taking v=∇·w in (2.3)(a) for w∈W and then substituting it into (2.4), we derive a new equivalent expanded mixed weak formulation of the system (2.3):
(2.5)(a)(λt,w)+(∇⋅σ,∇⋅w)=-(f(x,t),∇⋅w),∀w∈W,(b)(σ,z)-(aλ,z)-(b∫0t‍λds,z)=0,∀z∈W,(c)(ut,v)-(∇⋅σ,v)=(f(x,t),v),∀v∈L2(Ω).Let𝒯hu and 𝒯hσ be two families of quasi-regular partitions of the domain Ω, which may be the same one or not, such that the elements in the partitions have the diameters bounded by hu and hσ, respectively. Let Xhu⊂L2(Ω) and Vhσ⊂W be finite element spaces defined on the partitions 𝒯hu and 𝒯hσ.Now the semidiscrete positive definite expanded mixed finite element method for (2.5) consists in determining (uh,λh,σh)∈Xhu×Vhσ×Vhσ such that
(2.6)(a)(λht,wh)+(∇⋅σh,∇⋅wh)=-(f(x,t),∇⋅wh),∀wh∈Vhσ,(b)(σh,zh)-(aλh,zh)-(b∫0t‍λhds,zh)=0,∀zh∈Vhσ,(c)(uht,vh)-(∇⋅σh,vh)=(f(x,t),vh),∀vh∈Xhu,
with given an initial approximation (uh0,λh0,σh0)∈Xhu×Vhσ×Vhσ.Remark 2.1.
Compared to expanded mixed weak formulation (2.3), the new expanded mixed element system (2.6) is symmetric positive definite, that is to say the gradient function and the flux function system (2.6)(a,b) is symmetric positive definite. And both the gradient equation and the flux equation are separated from its scalar unknown equation (2.6)(c).Theorem 2.2.
There exists a unique discrete solution to the system (2.6).Proof.
Let{ψi(x)}i=1N1 and {φj(x)}j=1N2 be bases of Xhu and Vhσ, respectively. Then, we have the following expressions:
(2.7)uh=∑i=1N1ui(t)ψi(x),λh=∑j=1N2‍λj(t)φj(x),σh=∑j=1N2‍σj(t)φj(x).
Substituting these expressions into (2.6) and choosing vh=ψm,wh=zh=φl, then the problems (2.6) can be written in vector matrix form as: find {U(t),Λ(t),Σ(t)} such that, for all t∈(0,T](2.8)(a)AΛ′(t)+BΣ(t)=F(t),(b)AΣ(t)-CΛ(t)-H∫0t‍Λ(s)ds=0,(c)DU′(t)-EΣ(t)=G(t),
where
(2.9)A=((φj,φl))N2×N2,B=((∇⋅φj,∇⋅φl))N2×N2,C=((aφj,φl))N2×N2,D=((ψi,ψm))N1×N1,H=((bφj,φl))N2×N2,E=((∇⋅φj,ψm))N1×N2,U(t)=(u1(t),u2(t),…,uN1(t))T,Λ(t)=(λ1(t),λ2(t),…,λN2(t))T,Σ(t)=(σ1(t),σ2(t),…,σN2(t))T,F(t)=((-f(t),∇⋅φl))1×N2T,G(t)=((f(t),ψm))1×N1T.
It is easy to see that both A and D are symmetric positive definite. From (2.8), the problems can be written as follows:
(2.10)(a)Λ′(t)=-A-1BA-1(CΛ(t)+H∫0t‍Λ(s)ds)+A-1F(t),(b)Σ(t)=A-1CΛ(t)+A-1H∫0t‍Λ(s)ds,(c)U′(t)=D-1EΣ(t)+D-1G(t).
Thus, by the theory of differential equations [38, 39], (2.10) has a unique solution, and equivalently (2.6) has a unique solution.Remark 2.3.
It is easy to see that the coefficient matrixesA,B,C, and H of system (2.8) are symmetric positive definite. In view of this, the new expanded mixed element system (2.6) is symmetric positive definite.


## 3. Semidiscrete Error Estimates

LetXhu and Vhσ be finite dimensional subspaces of L2(Ω) and W, respectively, with the inverse property (see [36]) and the following approximation properties (see [40–44]): for 0≤p≤+∞ and r,r*,k positive integers
(3.1)inf⁡wh∈Vhσ‖w-wh‖Lp(Ω)≤Chσr+1‖w‖Wr+1,p(Ω),∀w∈H(div⁡;Ω)∩[Wr+1,p(Ω)]d,inf⁡wh∈Vhσ‖∇⋅(w-wh)‖Lp(Ω)≤Chσr*‖∇⋅w‖Wr*,p(Ω),∀w∈H(div⁡;Ω)∩[Wr+1,p(Ω)]d,inf⁡vh∈Xhu‖v-vh‖Lp(Ω)≤Chuk+1‖v‖Wk+1,p(Ω),∀v∈L2(Ω)∩Wk+1,p(Ω),
where r*=r+1 for the Brezzi-Douglas-Fortin-Marini spaces [43] and the Raviart-Thomas spaces [42] and r*=r for the Brezzi-Douglas-Marini spaces [40, 43].For our subsequent error analysis, we introduce two operators. It is well known that, in any one of the classical mixed finite element spaces, there exists an operatorRh from H(div;Ω) onto Vhσ, see [40–44], such that, for 1≤p≤+∞,
(3.2)(∇⋅(σ-Rhσ),ϕh)=0,∀ϕh∈∇⋅Vhσ={ϕh=∇⋅wh,wh∈Vhσ};‖σ-Rhσ‖Lp(Ω)≤Chσr+1‖σ‖Wr+1,p(Ω);‖∇⋅(σ-Rhσ)‖Lp(Ω)≤Chσr*‖∇⋅σ‖Wr*,p(Ω).
We also define the L2-project operator Ph from L2(Ω) onto Xhu such that
(3.3)(v-Phv,vh)=0,∀v∈L2(Ω),vh∈Xhu;‖v-Phv‖Lp(Ω)≤Chuk+1‖v‖Hk+1,p(Ω),∀v∈Hk+1(Ω).
Using the definitions of the operators Rh and Ph, we can easily obtain the following lemma.Lemma 3.1.
Assume that the solution of system (2.5) has the regular properties that ut∈L2(Hk+1(Ω)),    λt,λtt,σt∈L2(Hr+1(Ω)), then one has the following estimates:
(3.4)‖(λ-Rhλ)t‖Lp(Ω)≤Chσr+1‖λt‖Wr+1,p(Ω),‖(λ-Rhλ)tt‖Lp(Ω)≤Chσr+1‖λtt‖Wr+1,p(Ω),‖(σ-Rhσ)t‖Lp(Ω)≤Chσr+1‖σt‖Wr+1,p(Ω),‖(u-Phu)t‖Lp(Ω)≤Chuk+1‖ut‖Hk+1,p(Ω).
Let(3.5)u-uh=u-Phu+Phu-uh=η+ς,λ-λh=λ-Rhλ+Rhλ-λh=ρ+ξ,σ-σh=σ-Rhσ+Rhσ-σh=γ+θ.
Subtracting (2.6) from (2.5) and using projections (3.2) and (3.3), one obtains
(3.6)(ξt,wh)+(∇⋅θ,∇⋅wh)=-(ρt,wh),∀wh∈Vhσ,(θ,zh)-(aξ,zh)-(b∫0t‍ξds,zh)=-(γ,zh)+(aρ,zh)+(b∫0t‍ρds,zh),∀zh∈Vhσ,(ςt,vh)-(∇⋅(σ-σh),vh)=0,∀vh∈Xhu.Theorem 3.2.
Assume that the approximate properties (3.1) hold, and the solution of system (2.5) has regular properties that ut,utt∈L2(Hk+1(Ω)),λt,λtt,σt,σtt∈L2(Hr+1(Ω)). Then one has the error estimates
(3.7)‖λ-λh‖L∞(L2(Ω))+‖σ-σh‖L∞(L2(Ω))≤Chσr+1,‖∇⋅(σ-σh)‖L∞(L2(Ω))≤Chσr*,‖u-uh‖L∞(L2(Ω))≤C(hσr*+huk+1).Proof .
Choosewh=θ in (3.6)(a) and zh=-ξt in (3.6)(b), and add the two equations to obtain
(3.8)12ddt‖a1/2ξ‖L2(Ω)2+‖∇⋅θ‖L2(Ω)2=12(atξ,ξ)+(γ,ξt)-(aρ,ξt)-(ρt,θ)-(b∫0t‍ξds,ξt)-(b∫0t‍ρds,ξt)=12(atξ,ξ)+ddt(γ,ξ)-(γt,ξ)-ddt(aρ,ξ)+(aρt+atρ,ξ)-(ρt,θ)-ddt(b∫0t‍ξds,ξ)+(bt∫0t‍ξds,ξ)+(bξ,ξ)-ddt(b∫0t‍ρds,ξ)+(bt∫0t‍ρds,ξ)+(bρ,ξ).
Integrate with respect to time from 0 to t and apply the Cauchy-Schwarz's inequality and the Young's inequality to obtain
(3.9)‖a1/2ξ‖L2(Ω)2+2∫0t‍‖∇⋅θ‖L2(Ω)2ds≤‖a1/2ξ(0)‖L2(Ω)2+‖ρ‖L2(Ω)2+‖γ‖L2(Ω)2+ε‖ξ‖L2(Ω)2+C∫0t‍(‖ρ‖L2(Ω)2+‖ρs‖L2(Ω)2+‖γs‖L2(Ω)2+‖ξ‖L2(Ω)2+‖θ‖L2(Ω)2)ds.
Choosezh=θ in (3.6)(b) to get
(3.10)‖θ‖L2(Ω)2≤C(‖a1/2ξ‖L2(Ω)2+‖γ‖L2(Ω)2+‖ρ‖L2(Ω)2+∫0t‍(‖ξ‖L2(Ω)2+‖ρ‖L2(Ω)2)ds).
Combining (3.9) and (3.10), we obtain
(3.11)‖ξ‖L2(Ω)2+‖θ‖L2(Ω)2+2∫0t‍‖∇⋅θ‖L2(Ω)2ds≤‖a1/2ξ(0)‖L2(Ω)2+‖ρ‖L2(Ω)2+‖γ‖L2(Ω)2+ε‖ξ‖L2(Ω)2+C∫0t‍(‖ρ‖L2(Ω)2+‖ρs‖L2(Ω)2+‖γs‖L2(Ω)2+‖ξ‖L2(Ω)2+‖θ‖L2(Ω)2)ds.
Using the fact that ξ(0)=0 and Gronwall's lemma, we obtain
(3.12)‖ξ‖L2(Ω)2+‖θ‖L2(Ω)2+2∫0t‍‖∇⋅θ‖L2(Ω)2ds≤‖ρ‖L2(Ω)2+‖γ‖L2(Ω)2+C∫0t‍(‖ρ‖L2(Ω)2+‖ρs‖L2(Ω)2+‖γs‖L2(Ω)2)ds.
Differentiating (3.6)(b) and taking zh=ξt, we obtain
(3.13)(θt,ξt)=‖aξt‖L2(Ω)2+(atξ,ξt)-(γt,ξt)+(atρ+aρt,ξt)+(bt∫0t‍ξds,ξt)+(bξ,ξt)+(bt∫0t‍ρds,ξt)+(bρ,ξt).
Choosing wh=θt in (3.6)(a), we obtain
(3.14)(ξt,θt)+12ddt‖∇⋅θ‖L2(Ω)2=-ddt(ρt,θ)+(ρtt,θ).
Add (3.13) and (3.14) to obtain
(3.15)‖aξt‖L2(Ω)2+12ddt‖∇⋅θ‖L2(Ω)2=-ddt(ρt,θ)+(ρtt,θ)-(atξ,ξt)+(γt,ξt)-(atρ+aρt,ξt)-(bt∫0t‍ξds,ξt)-(bξ,ξt)-(bt∫0t‍ρds,ξt)-(bρ,ξt).
Integrate (3.15) with respect to time from 0 to t to obtain
(3.16)∫0t‍‖ξs‖L2(Ω)2ds+‖∇⋅θ‖L2(Ω)2≤‖ρ‖L2(Ω)2+‖γ‖L2(Ω)2+C(‖ξ‖L2(Ω)2+∫0t‍‖ξ‖L2(Ω)2ds)+C∫0t‍(‖ρ‖L2(Ω)2+‖ρs‖L2(Ω)2+‖ρss‖L2(Ω)2+‖γs‖L2(Ω)2)ds.
Substitute (3.12) into (3.16) to have
(3.17)∫0t‍‖ξs‖L2(Ω)2ds+‖∇⋅θ‖L2(Ω)2≤‖ρ‖L2(Ω)2+‖γ‖L2(Ω)2+C∫0t‍(‖ρ‖L2(Ω)2+‖γ‖L2(Ω)2+‖ρs‖L2(Ω)2+‖ρss‖L2(Ω)2+‖γs‖L2(Ω)2)ds.
Choosing vh=ς in (3.6)(b) and applying Cauchy-Schwarz's inequality, we obtain
(3.18)‖ς‖L2(Ω)≤‖∇⋅(σ-σh)‖L2(Ω).
Using Lemma 3.1, (3.17), and Gronwall's lemma, we get
(3.19)‖ς‖L∞(L2(Ω))≤C(hσr*+huk+1).
Using (3.12), (3.17), (3.19), (3.2), (3.3), and Lemma 3.1, we apply the triangle inequality to complete the proof.


## 4. Fully Discrete Error Estimates

In this section, we get the error estimates of fully discrete schemes. For the backward Euler procedure, let0=t0<t1<t2<⋯<tM=T be a given partition of the time interval [0,T] with step length δ=T/M, for some positive integer M. For a smooth function ϕ on [0,T], define ϕn=ϕ(tn) and ∂tϕn=(ϕn-ϕn-1)/δ. For approximating the integrals, we use the composite left rectangle rule
(4.1)δ∑j=0n-1ϕj≈∫0tn‍ϕ(s)ds.
Note that ϕ∈C1[0,T], the quadrature error satisfies
(4.2)|δ∑j=0n-1ϕj-∫0tn‍ϕ(s)ds|≤Cδ∫0tn‍|ϕs(s)|ds.Equation (2.5) has the following equivalent formulation:
(4.3)(a)(∂tλn,w)+(∇⋅σn,∇⋅w)=(-fn,∇⋅w)+(R1n,w),∀w∈W,(b)(σn,z)-(anλn,z)-(bnδ∑j=0n-1‍λj,z)=-(R3n,z),∀z∈W,(c)(∂tun,v)-(∇⋅σn,v)=(fn,v)+(R2n,v),∀v∈L2(Ω),
where
(4.4)R1n=∂tλn-λt=O(δ∂2λ∂t2),R2n=∂tun-ut=O(δ∂2u∂t2),R3n=δ∑j=0n-1‍λj-∫0tn‍λ(s)ds=O(δ∂λ∂t).
Now we can formulate a fully discrete scheme based on (4.3).Fully discrete scheme: find(uhn,λhn,σhn)∈Xhu×Vhσ×Vhσ,(n=1,2,…,M-1) such that
(4.5)(a)(∂tλhn,wh)+(∇⋅σhn,∇⋅wh)=(-fn,∇⋅wh),∀wh∈Vhσ,(b)(σhn,zh)-(anλhn,zh)-(bnδ∑j=0n-1‍λhj,zh)=0,∀zh∈Vhσ,(c)(∂tuhn,vh)-(∇⋅σhn,vh)=(fn,vh),∀vh∈Xhu,
with given an initial approximation (uh0,λh0,σh0)∈Xhu×Vhσ×Vhσ.For fully discrete error estimates, we now split the errors(4.6)un-uhn=un-Phun+Phun-uhn=ηn+ςn,λn-λhn=λn-Rhλn+Rhλn-λhn=ρn+ξn,σn-σhn=σn-Rhσn+Rhσn-σhn=γn+θn.
From (4.3) to (4.5), we then obtain
(4.7)(a)(∂tξn,wh)+(∇⋅θn,∇⋅wh)=-(∂tρn,wh)+(R1n,wh),∀wh∈Vhσ,(b)(θn,zh)-(anξn,zh)-(bnδ∑j=0n-1‍ξj,zh)=-(γn,zh)+(anρn,zh)+(bnδ∑j=0n-1‍ρj,zh)-(R3n,zh),∀zh∈Vhσ,(c)(∂tςn,vh)-(∇⋅(θn+γn),vh)=(R2n,vh),∀vh∈Xhu.Lemma 4.1.
Assume that the solution of system (2.5) has regular properties that ut∈L2(Hk+1(Ω)),λt, σt∈L2(Hr+1(Ω)). Then one has the estimates
(4.8)max⁡0≤n≤M‖∂t(λ-Rhλ)n‖L2(Ω)+max⁡0≤n≤M‖∂t(σ-Rhσ)n‖L2(Ω)≤Chσr+1,max⁡0≤n≤M‖∂t(u-Phu)n‖L2(Ω)≤Chuk+1.Theorem 4.2.
Assume that∂2u/∂t2,∂u/∂t∈L2(Hk+1(Ω)),∂λ/∂t,∂2λ/∂t2,∂σ/∂t,∂2σ/∂t2∈L2(Hr+1(Ω)), u∈L∞(Hk+1(Ω)), and λ,σ∈L∞(Hr+1(Ω)), then there exists a constant C such that
(4.9)max⁡0≤n≤M‖(λ-λh)n‖L2(Ω)+max⁡0≤n≤M‖(σ-σh)n‖L2(Ω)≤C(hσr+1+δ),max⁡0≤n≤M‖∇⋅(σ-σh)n‖L2(Ω)≤C(hσr*+δ),max⁡0≤n≤M‖(u-uh)n‖L2(Ω)≤C(huk+1+hσr*+δ).Proof .
Setwh=θn in (4.7)(a) and zh=-∂tξn in (4.7)(b) and add the two equations to obtain
(4.10)‖∇⋅θn‖L2(Ω)2+(an∂tξn,ξn)=-(∂tρn,θn)+(R1n,θn)+(γn,∂tξn)-(anρn,∂tξn)-(bnδ∑j=0n-1‍ξj,∂tξn)-(bnδ∑j=0n-1‍ρj,∂tξn)+(R3n,∂tξn).
Note that
(4.11)∂t‖anξn‖L2(Ω)2=(anξn,ξn)-(an-1ξn-1,ξn-1)δ=(anξn,ξn)-(anξn,ξn-1)+(anξn,ξn-1)-(an-1ξn-1,ξn-1)δ=(an∂tξn,ξn)+(an∂tξn,ξn-1)+(∂tanξn-1,ξn-1)=2(an∂tξn,ξn)-‖(an)1/2(ξn-ξn-1)‖L2(Ω)2δ+(∂tanξn-1,ξn-1).
So, we get
(4.12)(an∂tξn,ξn)=12∂t‖anξn‖L2(Ω)2+‖(an)1/2(ξn-ξn-1)‖L2(Ω)22δ-12(∂tanξn-1,ξn-1).
Note that
(4.13)(γn,∂tξn)=(ξn,γn)-(ξn-1,γn-1)δ-(∂tγn,ξn-1),(anρn,∂tξn)=(ξn,anρn)-(ξn-1,an-1ρn-1)δ-(∂tanρn,ξn-1),(R3n,∂tξn)=(ξn,R3n)-(ξn-1,R3n-1)δ-(∂tR3n,ξn-1).(bnδ∑j=0n-1ξj,∂tξn)=(bnδ∑j=0n-1‍ξj,ξn)-(bn-1δ∑j=0n-2‍ξj,ξn-1)δ-(∂tbnδ∑j=0n-1ξj,ξn-1)-(bn-1ξn-1,ξn-1),(bnδ∑j=0n-1ρj,∂tξn)=(bnδ∑j=0n-1‍ρj,ξn)-(bn-1δ∑j=0n-2‍ρj,ξn-1)δ-(∂tbnδ∑j=0n-1ρj,ξn-1)-(bn-1ρn-1,ξn-1).
Substitute (4.12)-(4.13) into (4.10) to get
(4.14)‖∇⋅θn‖L2(Ω)2+12∂t‖anξn‖L2(Ω)2+‖(an)1/2(ξn-ξn-1)‖L2(Ω)22δ=-(∂tρn,θn)+(R1n,θn)+12(∂tanξn-1,ξn-1)+(ξn,γn)-(ξn-1,γn-1)δ-(∂tγn,ξn-1)-(ξn,anρn)-(ξn-1,an-1ρn-1)δ+(∂tanρn,ξn-1)-(bnδ∑j=0n-1‍(ξj+ρj),ξn)-(bn-1δ∑j=0n-2‍(ξj+ρj),ξn-1)δ+(∂tbnδ∑j=0n-1‍(ξj+ρj),ξn-1)+(bn-1(ξn-1+ρn-1),ξn-1)+(ξn,R3n)-(ξn-1,R3n-1)δ-(∂tR3n,ξn-1).
Summing from 1 to n, we find that
(4.15)‖(an)1/2ξn‖L2(Ω)2+2δ∑j=1n‍‖∇⋅θj‖L2(Ω)2+∑j=1n‍‖(an)1/2(ξn-ξn-1)‖L2(Ω)2≤‖(a0)1/2ξ0‖L2(Ω)2-(ξ0,γ0)+(ξ0,a0ρ0)+ε‖ξn‖L2(Ω)2+C‖ρn‖L2(Ω)2+C‖R3n‖L2(Ω)2+C‖γn‖L2(Ω)2+Cδ∑j=1n-1‍‖ξj‖L2(Ω)2+Cδ∑j=1n‍‖θj‖L2(Ω)2+Cδ∑j=1n‍[‖∂tρj‖L2(Ω)2+‖∂tγj‖L2(Ω)2+‖∂tajρj‖L2(Ω)2+‖ρj‖L2(Ω)2+‖R1j‖L2(Ω)2+‖∂tR3j‖L2(Ω)2].
Choose zh=θn in (3.6)(b) to get
(4.16)‖θn‖L2(Ω)2≤C(‖(an)1/2ξn‖L2(Ω)2+‖γn‖L2(Ω)2+‖ρn‖L2(Ω)2)+Cδ∑j=1n-1‍(‖ξj‖L2(Ω)2+‖ρj‖L2(Ω)2).
Substitute (4.16) into (4.15) and note that ξ0=0 to get
(4.17)‖ξn‖L2(Ω)2+‖θn‖L2(Ω)2+2δ∑j=1n‍‖∇⋅θj‖L2(Ω)2≤C‖ρn‖L2(Ω)2+C‖γn‖L2(Ω)2+C‖R3n‖L2(Ω)2+Cδ∑j=1n-1‍‖ξj‖L2(Ω)2+Cδ∑j=1n‍‖θj‖L2(Ω)2+Cδ∑j=1n‍[‖∂tρj‖L2(Ω)2+‖∂tγj‖L2(Ω)2+‖∂tajρj‖L2(Ω)2+‖ρj‖L2(Ω)2+‖R1j‖L2(Ω)2+‖∂tR3j‖L2(Ω)2].
Using Gronwall's lemma, we obtain
(4.18)‖ξn‖L2(Ω)2+(1-Cδ)‖θn‖L2(Ω)2+2δ∑j=1n‍‖∇⋅θj‖L2(Ω)2≤C‖ρn‖L2(Ω)2+C‖γn‖L2(Ω)2+C‖R3n‖L2(Ω)2+Cδ∑j=1n‍[‖∂tρj‖L2(Ω)2+‖∂tγj‖L2(Ω)2+‖∂tajρj‖L2(Ω)2+‖ρj‖L2(Ω)2+‖R1j‖L2(Ω)2+‖∂tR3j‖L2(Ω)2].
Note that(4.19)(a)δ∑j=1n‍‖R1j‖L2(Ω)2≤Cδ3∑j=1n‍‖∂2λ∂t2‖L2(L2(Ω))2≤Cδ2‖∂2λ∂t2‖L∞(L2(Ω))2,(b)δ∑j=1n‍‖∂tρj‖L2(Ω)2≤Cδ∑j=1n‍‖∂ρj∂t‖L2(Ω)2≤CTM⋅nhσ2r+2≤Chσ2r+2,(c)δ∑j=1n‍‖∂tγj‖L2(Ω)2≤Cδ∑j=1n‍‖∂γj∂t‖L2(Ω)2≤CTM⋅nhσ2r+2≤Chσ2r+2,(d)δ∑j=1n‍‖∂tR3j‖L2(Ω)2=Cδ∑j=1n‍‖∫tj-1tj‍(λj-1-λ)dsδ‖L2(Ω)2≤Cδ2‖∂λ∂t‖L∞(L2(Ω))2.
Therefore, substituting the above estimates into (4.18) and choosing δ0 in such a way that for 0<δ≤δ0,(1-Cδ)>0, we obtain
(4.20)‖θn‖L2(Ω)2+‖ξn‖L2(Ω)2+δ∑j=1n‍‖∇⋅θj‖L2(Ω)2≤C(hσ2r+2+δ2).
By (4.7)(b), we obtain
(4.21)(∂tθn,zh)=(anξn-an-1ξn-1δ,zh)-(∂tγn,zh)+(anρn-an-1ρn-1δ,zh)+(bnδ∑j=0n-1‍(ξj+ρj),zh)-(bn-1δ∑j=0n-2‍(ξj+ρj),zh)δ-(∂tR3n,zh)=(anξn-an-1ξn-1δ,zh)-(∂tγn,zh)+(anρn-an-1ρn-1δ,zh)+(bn(ξn-1+ρn-1),zh)+(∂tbnδ∑j=0n-2‍(ξj+ρj),zh)-(∂tR3n,zh).
Set zh=∂tξn in (4.21) to obtain
(4.22)(∂tθn,∂tξn)=(anξn-an-1ξn-1δ,∂tξn)-(∂tγn,∂tξn)+(anρn-an-1ρn-1δ,∂tξn)+(bn(ξn-1+ρn-1),∂tξn)+(∂tbnδ∑j=0n-2‍(ξj+ρj),∂tξn)-(∂tR3n,∂tξn)=‖(an)1/2∂tξn‖L2(Ω)2+(∂tanξn-1,∂tξn)+(an∂tρn,∂tξn)+(∂tanρn-1,∂tξn)-(∂tγn,∂tξn)+(bn(ξn-1+ρn-1),∂tξn)+(∂tbnδ∑j=0n-2‍(ξj+ρj),∂tξn)-(∂tR3n,∂tξn).
Set wh=∂tθn in (4.7)(a) to obtain
(4.23)(∂tξn,∂tθn)+12∂t‖∇⋅θn‖L2(Ω)2+‖∇⋅(θn-θn-1)‖L2(Ω)22δ=-(∂tρn,∂tθn)+(R1n,∂tθn).
Substitute (4.22) into (4.23) to get
(4.24)‖(an)1/2∂tξn‖L2(Ω)2+12∂t‖∇⋅θn‖L2(Ω)2+‖∇⋅(θn-θn-1)‖L2(Ω)22δ=-(∂tρn,∂tθn)-(∂tanξn-1,∂tξn)-(an∂tρn,∂tξn)-(∂tanρn-1,∂tξn)+(∂tγn,∂tξn)-(bn(ξn-1+ρn-1),∂tξn)+(R1n,∂tθn)-(∂tbnδ∑j=0n-2‍(ξj+ρj),∂tξn)+(∂tR3n,∂tξn).
Take zh=∂tθn in (4.21) to obtain
(4.25)‖∂tθn‖L2(Ω)2=(anξn-an-1ξn-1δ,∂tθn)-(∂tγn,∂tθn)+(anρn-an-1ρn-1δ,∂tθn)+(bn(ξn-1+ρn-1),∂tθn)+(∂tbnδ∑j=0n-2‍(ξj+ρj),∂tθn)-(∂tR3n,∂tθn)=(an∂tξn,∂tθn)+(∂tanξn-1,∂tθn)-(∂tγn,∂tθn)+(an∂tρn,∂tθn)+(∂tanρn-1,∂tθn)+(bn(ξn-1+ρn-1),∂tθn)+(∂tbnδ∑j=0n-2‍(ξj+ρj),∂tθn)-(∂tR3n,∂tθn).
Add (4.24) and (4.25) to get
(4.26)‖(an)1/2∂tξn‖L2(Ω)2+‖∂tθn‖L2(Ω)2+12∂t‖∇⋅θn‖L2(Ω)2+‖∇⋅(θn-θn-1)‖L2(Ω)22δ=-(∂tρn,∂tθn)-(∂tanξn-1,∂tξn)-(an∂tρn,∂tξn)-(∂tanρn-1,∂tξn)+(∂tγn,∂tξn)+(an∂tξn,∂tθn)+(∂tanξn-1,∂tθn)-(∂tγn,∂tθn)+(an∂tρn,∂tθn)+(∂tanρn-1,∂tθn)+(R1n,∂tθn)-(∂tR3n,∂tθn+∂tξn)+(bn(ξn-1+ρn-1),∂tθn+∂tξn)+(∂tbnδ∑j=0n-2‍(ξj+ρj),∂tθn+∂tξn).
Apply Cauchy-Schwarz's inequality and Young's inequality to obtain
(4.27)amin⁡‖∂tξn‖L2(Ω)2+‖∂tθn‖L2(Ω)2+12∂t‖∇⋅θn‖L2(Ω)2≤ε(‖∂tξn‖L2(Ω)2+‖∂tθn‖L2(Ω)2)+C[‖ξn-1‖L2(Ω)2+‖∂tρn‖L2(Ω)2+‖∂tγn‖L2(Ω)2+‖ρn-1‖L2(Ω)2+‖R1n‖L2(Ω)2+‖∂tR3n‖L2(Ω)2]+Cδ2∑j=0n-2‍(‖ξj‖L2(Ω)2+‖ρj‖L2(Ω)2).
Using (4.19) and (4.20) and summing from 1 to n, we obtain
(4.28)‖∇⋅θn‖L2(Ω)2+δ∑j=1n‍(‖∂tξn‖L2(Ω)2+‖∂tθn‖L2(Ω)2)≤C(hσ2r+2+δ2).
Choosing vh=ςn in (4.7)(c) and applying Cauchy-Schwarz's inequality, Young's inequality, and (4.28), we have
(4.29)‖ςn‖L2(Ω)2-‖ςn-1‖L2(Ω)2=2δ(∇⋅(σ-σh)n,ςn)+2δ(R2n,ςn)≤Cδ[‖∇⋅(σ-σh)n‖L2(Ω)2+‖R2n‖L2(Ω)2+‖ςn‖L2(Ω)2].
Summing from 1 to n and using the Gronwall lemma, we obtain
(4.30)‖ςn‖L2(Ω)2≤‖ς0‖L2(Ω)2+Cδ∑j=1n‍[‖∇⋅(σ-σh)j‖L2(Ω)2+‖R2j‖L2(Ω)2].
Note that
(4.31)δ∑j=1n‍‖R2j‖L2(Ω)2≤Cδ3∑j=1n‍‖∂2u∂t2‖L2(L2(Ω))2≤Cδ2‖∂2u∂t2‖L∞(L2(Ω))2.
Substituting (4.31) into (4.30) and using (3.2), (4.20), and the triangle inequality, we get
(4.32)‖ςn‖L2(Ω)2≤C(hu2k+2+hσ2r*+δ2).
Combining (3.2), (3.3), (4.20), (4.28), (4.32), and Lemma 4.1, we apply the triangle inequality to complete the proof.


## 5. Numerical Example

In this section, we analyse some numerical results to illustrate the efficiency of the proposed method. We consider the following 2D parabolic partial integrodifferential equations with initial-boundary value condition:(5.1)ut-∇⋅(a(x,t)∇u+b(x,t)∫0t‍∇uds)=f(x,t),(x,t)∈Ω×J,u(x,t)=0,(x,t)∈∂Ω×J-,u(x,0)=sin⁡(πx1)sin⁡(πx2),x∈Ω,
where Ω=[0,1]×[0,1], J=(0,1], a(x,t)=1+x12+2x22, b(x,t)=1+2x12+x22, x=(x1,x2), and f(x,t) is chosen so that the exact solution for the scalar unknown function is
(5.2)u(x,t)=e-tsin⁡(πx1)sin⁡(πx2).
The corresponding exact gradient is
(5.3)λ(x,t)=(λ1,λ2)=(πe-tcos⁡(πx1)sin⁡(πx2),πe-tsin⁡(πx1)cos⁡(πx2)),
and its exact flux is
(5.4)σ(x,t)=(σ1,σ2)=(π(1+(2-e-t)x12+(1+e-t)x22)cos⁡⁡(πx1)sin⁡⁡(πx2),π(1+(2-e-t)x12+(1+e-t)x22)sin⁡(πx1)cos⁡(πx2)).Dividing the domainΩ into the triangulations of mesh size hu=hσ=h uniformly, considering the piecewise constant space Xhu with index k=0 for the scalar unknown function u and the lowest-order Raviart-Thomas triangular space Vhσ [42, 45] for the gradient λ and the flux σ and using the backward Euler procedure with uniform time step length δ=1/M, we obtain some convergence results for ∥u-uh∥L∞(L2(Ω)), ∥λ-λh∥L∞(L2(Ω)), ∥σ-σh∥L∞(L2(Ω)) and ∥σ-σh∥L∞(H(div;Ω)) with h=2δ=2/8,2/16,2/32 in Table 1. With time t=1, h=2δ=2/32, the exact solution u, λ, σ is shown in Figures 1, 3, and 5, respectively, and the corresponding numerical solution uh, λh, σh is shown in Figures 2, 4, and 6, respectively.Table 1
The errors and convergence order.
(
h
,
δ
)
∥
u
-
u
h
∥
L
∞
(
L
2
(
Ω
)
)
Order
∥
λ
-
λ
h
∥
L
∞
(
L
2
(
Ω
)
)
Order
(
2
8
,
1
8
)
1.4527
e
-
002
1.3532
e
-
001
(
2
16
,
1
16
)
6.5250
e
-
003
1.15
6.9202
e
-
002
(
2
32
,
1
32
)
3.0841
e
-
003
1.08
3.5363
e
-
002
(
h
,
δ
)
∥
σ
-
σ
h
∥
L
∞
(
L
2
(
Ω
)
)
Order
∥
σ
-
σ
h
∥
L
∞
(
H
(
div
;
Ω
)
)
(
2
8
,
1
8
)
3.3512
e
-
001
3.8892
e
-
001
(
2
16
,
1
16
)
1.6734
e
-
001
1.00
1.7291
e
-
001
(
2
32
,
1
32
)
8.3746
e
-
002
1.00
8.4431
e
-
002Figure 1
The exact solutionu.Figure 2
The numerical solutionuh.The exact gradientλ1(a) and λ2(b).
(a)
(b)The numerical gradientλ1h (a) and λ2h (b).
(a)
(b)The exact fluxσ1 (a) and σ2 (b).
(a)
(b)The numerical fluxσ1h(a) and σ2h(b).
(a)
(b)We can see from Table1 that the convergence rate is order 1 which confirms the theoretical results of Theorem 4.2 for the above chosen spaces Xhu and Vhσ. The numerical results in Table 1 and Figures 1–6 show that new positive definite expanded mixed scheme is efficient.


## 6.  Concluding Remarks

In the paper, a new expanded mixed finite element method based on a positive definite system is proposed for parabolic partial integrodifferential equation. Compared to expanded mixed method and standard mixed methods, the new expanded mixed element system is symmetric positive definite and both the gradient equation and the flux equation are separated from its scalar unknown equation. The existence and uniqueness for semidiscrete scheme are proved and error estimates are derived for both semidiscrete and fully discrete schemes. Finally, some numerical results are provided to confirm our theoretical analysis. In the near further, we will study the others evolution equations such as hyperbolic wave equation, and miscible displacement of compressible flow in porous media.



---

*Source: 391372-2012-06-28.xml*