# A Coupling Method of New EMFE and FE for Fourth-Order Partial Differential Equation of Parabolic Type

Authors: Yang Liu; Hong Li; Zhichao Fang; Siriguleng He; Jinfeng Wang
Journal: Advances in Mathematical Physics
 (2013)

---

## Abstract

We propose and analyze a new numerical method, called a coupling method based on a new expanded mixed finite element (EMFE) and finite element (FE), for fourth-order partial differential equation of parabolic type. We first reduce the fourth-order parabolic equation to a coupled system of second-order equations and then solve a second-order equation by FE method and approximate the other one by a new EMFE method. We find that the new EMFE method’s gradient belongs to the simple square integrable(L2(Ω))2 space, which avoids the use of the classical H(div; Ω) space and reduces the regularity requirement on the gradient solution λ=∇u. For a priori error estimates based on both semidiscrete and fully
discrete schemes, we introduce a new expanded mixed projection and some important lemmas.
We derive the optimal a priori error estimates in L2 and H1-norm for both the scalar unknown u and the diffusion term γ and a priori error estimates in (L2)2-norm for its gradient λ
and its flux σ (the coefficients times the negative gradient). Finally, we provide some numerical results to illustrate the efficiency of our method.


---

## Body

## 1. Introduction

In recent years, many researchers have studied some numerical methods for fourth-order elliptic equations [1–6], fourth-order parabolic equations [5–9], fourth-order wave equations [10, 11], and so on. Chen [1] proposed and analyzed an expanded mixed finite element method for fourth-order elliptic problems. In [2], Chen et al. studied an anisotropic nonconforming element for fourth-order elliptic singular perturbation problem. In [3–6], some mixed finite element (MFE) methods were studied for fourth-order linear/nonlinear elliptic equations. In [7], the FE method was studied for nonlinear Cahn-Hilliard equation. The optimal-order error estimates were obtained in L2-norm by means of an FE biharmonic projection approximation. In [12], the MFE methods were studied for solving a fourth-order nonlinear reaction diffusion equation. In [8], an H1-Galerkin MFE method was studied for solving the fourth-order parabolic partial differential equations. Optimal error estimates were derived for both semidiscrete and fully discrete schemes for problems in one space dimension, and error estimates were derived for semidiscrete scheme for several space dimensions, and the stability for fully discrete scheme was proved by the iteration method. In [9], the FE method was studied for fourth-order nonlinear parabolic equation. He et al. [10] studied and analyzed the explicit/implicit MFE methods for a class of fourth-order wave equations. Shi and Peng [13] studied the finite element methods for fourth-order eigenvalue problems on anisotropic meshes. Liu et al. [11] studied a C1-conforming FE method for nonlinear fourth-order wave equation. In this paper, we consider the following fourth-order partial differential equations of parabolic type:
(1)ut+∇·(b(x,t)∇(∇·(a(t)∇u)))=f(x,t),(x,t)∈Ω×J,u(x,t)=Δu(x,t)=0,(x,t)∈∂Ω×J-,u(x,0)=u0(x),x∈Ω,
where Ω⊂R2 is a bounded convex polygonal domain with the Lipschitz continuous boundary ∂Ω and J=(0,T] is the time interval with 0<T<∞. The initial value u0(x) and f(x,t) are given functions; coefficients a=a(t) and b=b(x,t) satisfy
(2)A1:0<a0≤a(t)≤a1<+∞,|at(t)|≤c0,A2:0<b0≤b(x,t)≤b1<+∞,|bt(x,t)|≤c0,
for some positive constants a0,  a1,  b0,  b1, and c0.Chen proposed the EMFE method for second-order linear/quasilinear elliptic equation [14–16] and fourth-order linear elliptic equation [1]. Then the method was applied to other partial differential equations [17, 18]. Compared to standard MFE methods, the EMFE method is expanded in the sense that three variables are explicitly approximated, namely, the scalar unknown, its gradient, and its flux. In recent years, many numerical method based on the expanded mixed scheme, such as two-grid EMFE methods [19–22], H1-Galerkin EMFE methods [23–25], expanded mixed covolume method [26], expanded mixed hybrid method [27], and positive definite EMFE method [28], have been proposed and discussed.In [29], we developed and analyzed a new EMFE method for second-order elliptic problems based on the mixed weak formulation [30–32]. The new EMFE method [29], whose gradient belongs to the square integrable (L2(Ω))2 space instead of the classical H(div;Ω) space, is different from Chen’s EMFE method [14–16].In this paper, we present a new coupling method of new EMFE scheme [29] and FE scheme for fourth-order partial differential equation of parabolic type. We first introduce an auxiliary variable γ=-∇·(a(t)∇u) to reduce the fourth-order parabolic equation to a coupled system of second-order equations, then solve a second-order equation by FE method. For the other one, we introduce the two auxiliary variables, λ=∇u and σ=-a(t)∇u=-aλ, to split the equation into a first-order system of equations, then approximate that by a new EMFE method. Because the new EMFE method’s gradient belongs to the simple square integrable (L2(Ω))2 space, the regularity requirement on the gradient solution λ=∇u is reduced. We apply the new coupling method of new EMFE and FE to derive a priori error estimates based on both semidiscrete and fully discrete schemes; that is to say, we obtain the optimal a priori error estimates in L2 and H1-norm for both the scalar unknown u and the diffusion term γ and a priori error estimates in (L2)2-norm for its gradient λ and its flux σ. Finally, we provide some numerical results to verify our theoretical analysis.Throughout this paper,C will denote a generic positive constant which does not depend on the spatial mesh parameter h or the time step Δt. At the same time, we denote the natural inner product in L2(Ω) or (L2(Ω))2 by (·,·) with the corresponding norm ∥·∥. The other notations and definitions for the Sobolev spaces as in [33, 34] are used.


## 2. New Expanded Mixed Scheme

We first introduce an auxiliary variableγ=-∇·(a(t)∇u) to obtain the following coupled system of second-order equations: (3a)ut-∇·(b(x,t)∇γ)=f(x,t),(3b)γ+∇·(a(t)∇u)=0,(3c)u(x,t)=γ(x,t)=0.For (3b), we introduce the two auxiliary variables, λ=∇u and σ=-a(t)∇u=-aλ, to obtain the following first-order system of equations: (4a)ut-∇·(b∇γ)=f(x,t),(4b)γ-∇·σ=0,(4c)σ+aλ=0,(4d)λ-∇u=0,(4e)u(x,t)=γ(x,t)=0.Using Green’s formula, the new mixed weak formulation of (4a)–(4e) is to determine {u,γ,λ,σ}:[0,T]↦H01(Ω)×H01(Ω)×(L2(Ω))2×(L2(Ω))2 such that (5a)(ut,v)+(b∇γ,∇v)=(f,v),∀v∈H01(Ω),(5b)(γ,g)+(σ,∇g)=0,∀g∈H01(Ω),(5c)(σ,z)+(aλ,z)=0,∀z∈(L2(Ω))2,(5d)(λ,w)-(∇u,w)=0,∀w∈(L2(Ω))2.Then, the semidiscrete coupled EMFE scheme-FE scheme for (5a)–(5d) is to find {uh,γh,λh,σh}:[0,T]↦Vh×Vh×Wh×Wh such that (6a)(uht,vh)+(b∇γh,∇vh)=(f,vh),∀vh∈Vh,(6b)(γh,gh)+(σh,∇gh)=0,∀gh∈Vh,(6c)(σh,zh)+(aλh,zh)=0,∀zh∈Wh,(6d)(λh,wh)-(∇uh,wh)=0,∀wh∈Wh.By the theory of differential equations [35], we can prove the existence and uniqueness of solution for semidiscrete scheme (6a)–(6d).Remark 1.
From the new coupling weak formulation (5a)–(5d), we can find that the gradient belongs to the weaker square integrable (L2(Ω))2 space taking the place of the classical H(div;Ω) space. It is easy to see that our method reduces the regularity requirement on the gradient solution λ=∇u.Remark 2.
In the semidiscrete coupling scheme (6a)–(6d), the MFE space (Vh,Wh) based on the FE pair P1-P02 is chosen as follows:
(7)Vh={vh∈C0(Ω)∩H01(Ω)∣vh∈P1(K),∀K∈𝒦h},Wh={ttwh=(w1h,w2h)∈(L2(Ω))2∣w1h,w2h∈P0(K),∀K∈𝒦htt}.
From [30, 31], we know that MFE space (Vh,Wh) satisfies the discrete LBB condition.Remark 3.
From (6a)–(6d), we know that (6a) is an FE scheme for uh,γh,vh∈Vh and that (6b)–(6d) is a new EMFE scheme for (γh,σh,λh), (gh,zh,wh)∈Vh×Wh×Wh. So the scheme (6a)–(6d) is called as the coupling system of new EMFE scheme and FE scheme.Remark 4.
The finite element spacesVh and Wh for the new EMFE method are chosen as finite dimensional subspaces of H01(Ω) and (L2(Ω))2, respectively. However, Vh and Wh for H1-Galerkin MFE method in two space variables in [8] are finite dimensional subspaces of H01(Ω) and H(div;Ω), respectively. Compared to the complex H(div;Ω) space used in the H1-Galerkin MFE method, we use the weaker (L2(Ω))2 space for our EMFE method.


## 3.  Semidiscrete Error Estimates

For deriving some a priori error estimates of our method, we first introduce the projection operator𝒫h and the new expanded mixed elliptic projection operator ℛh associated with the coupled equations.Lemma 5.
One now defines a linear projection operator𝒫h:H01(Ω)→Vh as
(8)(b∇(γ-𝒫hγ),∇vh)=0,∀vh∈Vh,
and then
(9)∥γ-𝒫hγ∥L2(Ω)+h∥γ-𝒫hγ∥H1≤Ch2∥γ∥H2,∥γt-𝒫hγt∥L2(Ω)≤Ch2(∥γ∥H2+∥γt∥H2),∥γtt-𝒫hγtt∥L2(Ω)≤Ch2(∥γ∥H2+∥γt∥H2+∥γtt∥H2).Let(ℛhu,ℛhλ,ℛhσ):[0,T]→Vh×Wh×Wh be given by the following new expanded mixed relations [29] (10a)(σ-ℛhσ,∇gh)=0,∀gh∈Vh,(10b)(σ-ℛhσ,zh)+(a(λ-ℛhλ),zh)=0,∀zh∈Wh,(10c)(λ-ℛhλ,wh)-(∇(u-ℛhu),wh)=0,∀wh∈Wh.Then the following two important lemmas based on the new expanded mixed projection (10a)–(10c) hold.Lemma 6.
There is a constantC>0 independent of the spatial mesh parameter h such that
(11)∥λ-ℛhλ∥≤Ch(∥λ∥(H1)2+∥u∥H2),∥σ-ℛhσ∥≤Ch(∥u∥H2+∥σ∥(H1)2+∥λ∥(H1)2),∥∇(u-ℛhu)∥≤Ch(∥u∥H2+∥λ∥(H1)2).Lemma 7.
There is a constantC>0 independent of the spatial mesh parameter h such that
(12)∥u-ℛhu∥≤Ch2∥u∥H2,∥ut-ℛhut∥≤Ch2∥ut∥H2,∥utt-ℛhutt∥≤Ch2∥utt∥H2,∥u-ℛhu∥1≤Ch(∥u∥H2+∥λ∥(H1)2).In [29], we can obtain the detailed proof for Lemmas 6-7.For a priori error estimates, we write the errors as(13)γ-γh=γ-𝒫hγ+𝒫hγ-γh=ψ+ϖ;u-uh=u-ℛhu+ℛhu-uh=η+ς;λ-λh=λ-ℛhλ+ℛhλ-λh=δ+θ;σ-σh=σ-ℛhσ+ℛhσ-σh=ρ+ξ.Using (5a), (5b), (5c), (5d), (6a), (6b), (6c), (6d), (8), (10a), (10b), and (10c) we can obtain the error equations (14a)(ςt,vh)+(b∇ϖ,∇vh)=-(ηt,vh),∀vh∈Vh,(14b)(ϖ,gh)+(ξ,∇gh)=-(ψ,gh),∀gh∈Vh,(14c)(ξ,zh)+(aθ,zh)=0,∀zh∈Wh,(14d)(θ,wh)-(∇ς,wh)=0,∀wh∈Wh.We will prove the error estimates for semidiscrete scheme.Theorem 8.
Suppose thatγ,γt,ut∈L2(H2(Ω)), λ,σ∈L∞((H1(Ω))2), and u,γ∈L∞(H2(Ω)); then there exists a constant C>0 independent of the spatial mesh parameter h such that
(15)∥u-uh∥≤Ch2[(∫0t‍∥|+|∥2ds)1/2∥u∥L∞(H2)+∥γ∥L∞(H2)+(∫0t‍|||■|||2ds)1/2],∥u-uh∥1≤Ch[(∫0t‍∥|+|∥2ds)1/2∥λ∥L∞((H1)2)+∥u∥L∞(H2)+∥γ∥L∞(H2)+(∫0t‍|||■|||2ds)1/2],∥λ-λh∥≤Ch[(∫0t‍∥|+|∥2ds)1/2∥λ∥L∞((H1)2)+∥γ∥L∞(H2)+(∫0t‍|||■|||2ds)1/2],∥σ-σh∥≤Ch[(∫0t‍∥|+|∥2ds)1/2∥u∥L∞(H2)+∥σ∥L∞((H1)2)+∥λ∥L∞((H1)2)+∥γ∥L∞(H2)+(∫0t‍|||■|||2ds)1/2],
where |||■|||2≜∥ut∥H22+∥γ∥H22+∥γt∥H22.Proof.
Choosevh=ϖ, gh=ςt, zh=∇ςt,  and  wh=∇ςt in (14a)–(14d) to obtain (16a)(ςt,ϖ)+∥b1/2∇ϖ∥2=-(ηt,ϖ),(16b)-(ϖ,ςt)-(ξ,∇ςt)=(ψ,ςt),(16c)(ξ,∇ςt)+(aθ,∇ςt)=0,(16d)-(aθ,∇ςt)+12ddt∥a1/2∇ς∥2=12(at∇ς,∇ς).Adding the above four equations and using Cauchy-Schwarz inequality and Young’s inequality, we have
(17)12ddt∥a1/2∇ς∥2+∥b1/2∇ϖ∥2=-(ηt,ϖ)+(ψ,ςt)+12(at∇ς,∇ς)=-(ηt,ϖ)+ddt(ψ,ς)-(ψt,ς)+12(at∇ς,∇ς)≤ddt(ψ,ς)+∥ηt∥∥ϖ∥+C(∥ς∥12+∥ψt∥2).
Integrate with respect to time from 0 to t to obtain
(18)a0∥∇ς∥2+2b0∫0t‍∥∇ϖ∥2ds≤2(ψ,ς)+2∫0t‍∥ηt∥∥ϖ∥ds+C∫0t‍(∥ς∥12+∥ψt∥2)ds.
Noting that ς,ϖ∈H01(Ω), we use Poincaré inequality to have
(19)∥ς∥1≤C1∥∇ς∥,∥ϖ∥≤C2∥∇ϖ∥.
Using (18), (19), Cauchy-Schwarz inequality and Young inequality, we have
(20)a0∥∇ς∥2+2b0∫0t‍∥∇ϖ∥2ds≤2∥ψ∥∥ς∥+2∫0t‍∥ηt∥∥ϖ∥ds+C∫0t‍(∥ς∥12+∥ψt∥2)ds≤2C1∥ψ∥∥∇ς∥+2C2∫0t‍∥ηt∥∥∇ϖ∥ds+C∫0t‍(C12∥∇ς∥2+∥ψt∥2)ds≤8C12a0∥ψ∥2+a02∥∇ς∥2+4C22b0∫0t‍∥ηt∥2ds+b0∫0t‍∥∇ϖ∥2ds+C∫0t‍(C12∥∇ς∥2+∥ψt∥2)ds.
Using Gronwall’s lemma, we get
(21)∥∇ς∥2+∫0t‍∥∇ϖ∥2ds≤C(∥ψ∥2+∫0t‍(∥ηt∥2+∥ψ∥2+∥ψt∥2)ds).
Substitute (19) into (21) to have
(22)∥ς∥2+∫0t‍∥ϖ∥2ds≤C(∥ψ∥2+∫0t‍(∥ηt∥2+∥ψ∥2+∥ψt∥2)ds).
Taking wh=θ in (14d) and using (22), we get
(23)∥θ∥2≤C(∥ψ∥2+∫0t‍(∥ηt∥2+∥ψ∥2+∥ψt∥2)ds).
Choose zh=ξ in (14c) and use (22) to obtain
(24)∥ξ∥2≤C(∥ψ∥2+∫0t‍(∥ηt∥2+∥ψ∥2+∥ψt∥2)ds).
Combining Lemmas 5–7, (21)–(24) and using the triangle inequality, we get the error estimates for Theorem 8.Remark 9.
From Theorem8, we can see that
(25)∥u-uh∥L∞(L2)+h∥u-uh∥L∞(H1)=O(h2);
that is to say, a priori error estimates in L2 and H1-norms for the scalar unknown u are optimal.Theorem 10.
Suppose that∂ju/∂tj,∂jγ/∂tj∈L2(H2(Ω)) and γ,γt,ut∈L∞(H2(Ω)); then there exists a constant C>0 independent of the spatial mesh parameter h such that
(26)∥γ-γh∥+h∥γ-γh∥1≤Ch2[∑j=02‍∥∂ju∂tj∥H22∥ut∥L∞(H2)+∥γ∥L∞(H2)+∥γt∥L∞(H2)+(∫0t‍(∑j=02‍∥∂ju∂tj∥H22+∥∂jγ∂tj∥H22)ds)1/2].Proof.
Differentiating (14b), (14c), and (14d) with respect to time t, respectively, we get (27a)(ςt,vh)+(b∇ϖ,∇vh)=-(ηt,vh),(27b)(ϖt,gh)+(ξt,∇gh)=-(ψt,gh),(27c)(ξt,zh)+(atθ,zh)+(aθt,zh)=0,(27d)(θt,wh)-(∇ςt,wh)=0.Choose vh=ϖt, gh=ςt, zh=∇ςt,  and  wh=∇ςt in (27a)–(27d) to get (28a)(ςt,ϖt)+12ddt∥b1/2∇ϖ∥2=-(ηt,ϖt)+12(bt∇ϖ,∇ϖ),(28b)-(ϖt,ςt)-(ξt,∇ςt)=(ψt,ςt),(28c)(ξt,∇ςt)+(atθ,∇ςt)+(aθt,∇ςt)=0,(28d)-(aθt,∇ςt)+∥a1/2∇ςt∥2=0.Adding the four equations for (28a)–(28d), we have
(29)12ddt∥b1/2∇ϖ∥2+∥a1/2∇ςt∥2=-(ηt,ϖt)+(ψt,ςt)-(atθ,∇ςt)+12(bt∇ϖ,∇ϖ)=-ddt(ηt,ϖ)+(ηtt,ϖ)+ddt(ψt,ς)-(ψtt,ς)-(atθ,∇ςt)+12(bt∇ϖ,∇ϖ).
Integrate with respect to time from 0 to t and use (19), Cauchy-Schwarz inequality, and Young inequality to get
(30)b0∥∇ϖ∥2+2a0∫0t‍∥∇ςt∥2ds≤-2(ηt,ϖ)+2C2∫0t‍∥ηtt∥∥∇ϖ∥ds+2(ψt,ς)+C∫0t‍(∥ψtt∥2+∥ς∥2+∥θ∥2)ds+a0∫0t‍∥∇ςt∥2ds≤8C22b0∥ηt∥2+b02∥∇ϖ∥2+∥ψt∥2+∥ς∥2+a0∫0t‍∥∇ςt∥2ds+C∫0t‍(∥ψtt∥2+∥ς∥2+∥θ∥2+∥ηtt∥2+∥∇ϖ∥2)ds.
Substituting (22) and (23) into (30) and using Gronwall lemma, we get
(31)∥∇ϖ∥2+∫0t‍∥∇ςt∥2ds≤C∥ηt∥2+∥ψ∥2+∥ψt∥2+C∫0t‍(∥η∥2+∥ηt∥2+∥ηtt∥2+∥ψ∥2+∥ψt∥2+∥ψtt∥2)ds.
Use (19) to obtain
(32)∥ϖ∥2≤C∥ηt∥2+∥ψ∥2+∥ψt∥2+C∫0t‍(∥η∥2+∥ηt∥2+∥ηtt∥2+∥ψ∥2+∥ψt∥2+∥ψtt∥2)ds.
Combining Lemmas 5–7, (31), and (32) and using the triangle inequality, we obtain the error estimates for Theorem 10.Remark 11.
From Theorem10, we can see that
(33)∥γ-γh∥L∞(L2)+h∥γ-γh∥L∞(H1)=O(h2);
that is to say, a priori error estimates in L2 and H1-norm for the diffusion term γ are optimal.


## 4. Fully Discrete Error Estimates

In the following analysis, we will derive a priori error estimates based on fully discrete backward Euler scheme. Let0=t0<t1<t2<⋯<tM=T be a given partition of the time interval [0,T] with step length Δt=T/M and nodes tn=nΔt, for some positive integer M. For a smooth function ϕ on [0,T], define ϕn=ϕ(tn).For (5a), (5b), (5c), and (5d), we have the following equivalent formulation:(34a)(un-un-1Δt,v)+(bn∇γn,∇v)=(fn+R1n,v),∀v∈H01(Ω),(34b)(γn,g)+(σn,∇g)=0,∀g∈H01(Ω),(34c)(σn,z)+(anλn,z)=0,∀z∈(L2(Ω))2,(34d)(λn,w)-(∇un,w)=0,∀w∈(L2(Ω))2,where
(35)R1n=un-un-1Δt-ut(tn)=1Δt∫tn-1tn‍(tn-1-s)uttds.
Now a fully discrete procedure is to find (uhn,γhn,λhn,σhn)∈Vh×Vh×Wh×Wh(n=0,1,…,M) such that (36a)(uhn-uhn-1Δt,vh)+(bn∇γhn,∇vh)=(fn,vh),∀vh∈Vh,(36b)(γhn,gh)+(σhn,∇gh)=0,∀gh∈Vh,(36c)(σhn,zh)+(anλhn,zh)=0,∀zh∈Wh,(36d)(λhn,wh)-(∇uhn,wh)=0,∀wh∈Wh.For deriving the fully discrete error estimates, we now decompose the errors as
(37)γn-γhn=γn-𝒫hγn+𝒫hγn-γhn=ψn+ϖn;un-uhn=un-ℛhun+ℛhun-uhn=ηn+ςn;λn-λhn=λn-ℛhλn+ℛhλn-λhn=δn+θn;σn-σhn=σn-ℛhσn+ℛhσn-σhn=ρn+ξn.
We will prove the theorem for the fully discrete error estimates.Theorem 12.
Suppose thatγ,γt,ut∈L2(H2(Ω)), λ,σ∈L∞((H1(Ω))2), γ∈L∞(H2(Ω)), and utt∈L2(L2(Ω)); then there exists a constant C>0 independent of the spatial mesh parameter h and the time step Δt such that
(38)∥uJ-uhJ∥≤Ch2(∥u∥L∞(H2)+|||◆|||)+CΔt∥utt∥L2(L2),∥uJ-uhJ∥1≤Ch(∥u∥L∞(H2)+∥λ∥L∞((H1)2)+|||◆|||)+CΔt∥utt∥L2(L2),∥λJ-λhJ∥≤Ch(∥λ∥L∞((H1)2)+|||◆|||)+CΔt∥utt∥L2(L2),∥σJ-σhJ∥≤Ch(∥σ∥L∞((H1)2)+|||◆|||)+CΔt∥utt∥L2(L2),
where |||◆|||≜∥ut∥L2(H2)+∥γ∥L∞(H2)+∥γ∥L2(H2)+∥γt∥L2(H2).Proof.
Using (8), (10a), (10b), (10c), (34a), (34b), (34c), (34d), (36a), (36b), (36c), and (36d) at t=tn, we obtain the error equations(39a)(ςn-ςn-1Δt,vh)+(bn∇ϖn,∇vh)=-(ηn-ηn-1Δt,vh)+(R1n,vh),∀vh∈Vh,(39b)(ϖn,gh)+(ξn,∇gh)=-(ψn,gh),∀gh∈Vh,(39c)(ξn,zh)+(anθn,zh)=0,∀zh∈Wh,(39d)(θn,wh)-(∇ςn,wh)=0,∀wh∈Wh.Take vh=ϖn in (39a), gh=(ςn-ςn-1)/Δt in (39b), zh=(∇ςn-∇ςn-1)/Δt in (39c), and wh=(∇ςn-∇ςn-1)/Δt in (39d) to get
(40)∥(bn)1/2∇ϖn∥2=-(ςn-ςn-1Δt,ϖn)-(ηn-ηn-1Δt,ϖn)+(R1n,ϖn)=(ξn,∇ςn-∇ςn-1Δt)-(ηn-ηn-1Δt,ϖn)+(R1n,ϖn)+(ψn,ςn-ςn-1Δt)=-(anθn,∇ςn-∇ςn-1Δt)-(ηn-ηn-1Δt,ϖn)+(R1n,ϖn)+(ψn,ςn-ςn-1Δt)=-(an∇ςn-∇ςn-1Δt,∇ςn)-(ηn-ηn-1Δt,ϖn)+(R1n,ϖn)+(ψn,ςn-ςn-1Δt)=-12Δt(∥(an)1/2∇ςn∥2-∥(an-1)1/2∇ςn-1∥2)-12Δt∥an(∇ςn-∇ςn-1)∥2+12(an-an-1Δt∇ςn-1,∇ςn-1)-(ηn-ηn-1Δt,ϖn)+(R1n,ϖn)+(ψn,ςn-ςn-1Δt).
By (40), we have
(41)∥(bn)1/2∇ϖn∥2+12Δt∥an(∇ςn-∇ςn-1)∥2+12Δt(∥(an)1/2∇ςn∥2-∥(an-1)1/2∇ςn-1∥2)=12(an-an-1Δt∇ςn-1,∇ςn-1)-(ηn-ηn-1Δt,ϖn)+(R1n,ϖn)+(ψn,ςn-ςn-1Δt)=12(1Δt∫tn-1tn‍at(s)ds∇ςn-1,∇ςn-1)-(ηn-ηn-1Δt,ϖn)+(R1n,ϖn)+(ψn,ςn)-(ψn-1,ςn-1)Δt-(ςn-1,ψn-ψn-1Δt).
Multiplying by 2Δt, summing (41) from n=1 to J, and using (19), the resulting equation becomes
(42)a0∥∇ςJ∥2+2b0Δt∑n=1J‍∥∇ϖn∥2≤2Δt∑n=1J‍(∥ηn-ηn-1Δt∥+∥R1n∥)∥ϖn∥+supt∈[tn-1,tn]|at(t)|Δt∑n=1J-1‍∥∇ςn∥2+Δt∑n=1J‍(∥ςn-1∥2+∥ψn-ψn-1Δt∥2)+∥ψJ∥∥ςJ∥≤CΔt∑n=1J‍(∥ηn-ηn-1Δt∥2+∥R1n∥2)+b0Δt∑n=1J‍∥∇ϖn∥2+CΔt∑n=1J‍(∥∇ςn-1∥2+∥ψn-ψn-1Δt∥2)+C∥ψJ∥2+a02∥∇ςJ∥2.
Use (42) and Gronwall lemma to get
(43)∥∇ςJ∥2+Δt∑n=1J‍∥∇ϖn∥2≤C∥ψJ∥2+CΔt∑n=1J‍(∥ηn-ηn-1Δt∥2+∥ψn-ψn-1Δt∥2+∥R1n∥2).
Note that
(44)∥R1n∥2≤CΔt∫tn-1tn‍∥utt∥2ds,∥ηn-ηn-1Δt∥2≤1Δt∫tn-1tn‍∥ηt(s)∥2ds,∥ψn-ψn-1Δt∥2≤1Δt∫tn-1tn‍∥ψt(s)∥2ds.
Substitute (44) into (43) to get
(45)∥∇ςJ∥2+Δt∑n=1J‍∥∇ϖn∥2≤C(∥ψJ∥2+∫t0tJ‍(∥ηt(s)∥2+∥ψt(s)∥2)ds+(Δt)2∫t0tJ‍∥utt∥2ds).
Combining (19) and (45), we obtain
(46)∥ςJ∥2+Δt∑n=1J‍∥ϖn∥2≤C(∥ψJ∥2+∫t0tJ‍(∥ηt(s)∥2+∥ψt(s)∥2)ds+(Δt)2∫t0tJ‍∥utt∥2ds).
Taking wh=θn in (39d) and zh=ξn in (39c) and using (45), we have
(47)∥θJ∥2+∥ξJ∥2≤C(∥ψJ∥2+∫t0tJ‍(∥ηt(s)∥2+∥ψt(s)∥2)ds+(Δt)2∫t0tJ‍∥utt∥2ds).
Combining Lemmas 5–7, (45), (46), and (47) with the triangle inequality, we accomplish the proof for Theorem 12.In the following analysis, we will derive the optimal a priori error estimates inL2 and H1-norm for γ.Theorem 13.
Suppose thatγ,γt,ut∈L2(H2(Ω)),  ∂2u/∂t2,  ∂4u/∂t4∈L2(L2(Ω)),  γ∈L∞(H2(Ω)), and ∂3u/∂t3∈L∞(L2(Ω)); then there exists a constant C>0 independent of the spatial mesh parameter h and time discretization parameter Δt such that
(48)∥γJ-γhJ∥j≤Ch2-j(∥γ∥L∞(H2)+∥γ∥L2(H2)+∥ut∥L2(H2)+∥γt∥L2(H2))+Δt(∥∂2u∂t2∥L2(L2)+∥∂3u∂t3∥L∞(L2)+∥∂4u∂t4∥L2(L2)),j=0,1.Proof.
Use (39a)–(39d) to get(49a)(ςn-ςn-1Δt,vh)+(bn∇ϖn,∇vh)=-(ηn-ηn-1Δt,vh)+(R1n,vh),∀vh∈Vh,(49b)(ϖn-ϖn-1Δt,gh)+(ξn-ξn-1Δt,∇gh)=-(ψn-ψn-1Δt,gh),∀gh∈Vh,(49c)(ξn-ξn-1Δt,zh)+(anθn-an-1θn-1Δt,zh)=0,∀zh∈Wh,(49d)(θn-θn-1Δt,wh)-(∇ςn-∇ςn-1Δt,wh)=0,∀wh∈Wh.Taking vh=(ϖn-ϖn-1)/Δt,  gh=(ςn-ςn-1)/Δt,  zh=(∇ςn-∇ςn-1)/Δt, and wh=(∇ςn-∇ςn-1)/Δt, respectively, in (49a)–(49d), we get
(50)12Δt(∥(bn)1/2∇ϖn∥2-∥(bn-1)1/2∇ϖn-1∥2)+12Δt∥bn(∇ϖn-∇ϖn-1)∥2=(bn∇ϖn,∇ϖn-ϖn-1Δt)+12(bn-bn-1Δt∇ϖn-1,∇ϖn-1)=-(ςn-ςn-1Δt,ϖn-ϖn-1Δt)-(ηn-ηn-1Δt,ϖn-ϖn-1Δt)+(R1n,ϖn-ϖn-1Δt)+12(bn-bn-1Δt∇ϖn-1,∇ϖn-1)=(ξn-ξn-1Δt,∇ςn-∇ςn-1Δt)+(ψn-ψn-1Δt,ςn-ςn-1Δt)+(R1n,ϖn-ϖn-1Δt)-(ηn-ηn-1Δt,ϖn-ϖn-1Δt)+12(bn-bn-1Δt∇ϖn-1,∇ϖn-1)=-(an(θn-θn-1)Δt,∇ςn-∇ςn-1Δt)-(an-an-1Δtθn-1,∇ςn-∇ςn-1Δt)+(ψn-ψn-1Δt,ςn-ςn-1Δt)+(R1n,ϖn-ϖn-1Δt)-(ηn-ηn-1Δt,ϖn-ϖn-1Δt)+12(bn-bn-1Δt∇ϖn-1,∇ϖn-1)=-∥(an)1/2∇ςn-∇ςn-1Δt∥2-(an-an-1Δtθn-1,∇ςn-∇ςn-1Δt)+(ψn-ψn-1Δt,ςn-ςn-1Δt)+(R1n,ϖn-ϖn-1Δt)-(ηn-ηn-1Δt,ϖn-ϖn-1Δt)+12(bn-bn-1Δt∇ϖn-1,∇ϖn-1).
By using (50), we obtain
(51)12Δt(∥(bn)1/2∇ϖn∥2-∥(bn-1)1/2∇ϖn-1∥2)+12Δt∥bn(∇ϖn-∇ϖn-1)∥2+∥(an)1/2∇ςn-∇ςn-1Δt∥2=-(an-an-1Δtθn-1,∇ςn-ςn-1Δt)+(ψn-ψn-1Δt,ςn-ςn-1Δt)+(ϖn,R1n)-(ϖn-1,R1n-1)Δt-(R1n-R1n-1Δt,ϖn-1)-((ηn-ηn-1)/Δt,ϖn)-((ηn-1-ηn-2)/Δt,ϖn-1)Δt+(ηn-2ηn-1+ηn-2(Δt)2,ϖn-1)+12(bn-bn-1Δt∇ϖn-1,∇ϖn-1).
Multiply (51) by 2Δt, sum from n=1 to J, and use (19) to get
(52)b0∥∇ϖJ∥2+2a0Δt∑n=1J‍∥∇ςn-∇ςn-1Δt∥2≤2Δt∑n=1J‍(∥an-an-1Δtθn-1∥∥∇ςn-∇ςn-1Δt∥+∥ψn-ψn-1Δt∥∥ςn-ςn-1Δt∥)+2∥ϖJ∥∥R1J∥+2Δt∑n=1J‍∥R1n-R1n-1Δt∥∥ϖn-1∥+2∥ηJ-ηJ-1Δtc∥∥ϖJ∥+2Δt∑n=1J‍∥ηn-2ηn-1+ηn-2(Δt)2∥∥ϖn-1∥+CΔt∑n=1J‍∥∇ϖn-1∥2≤CΔt∑n=1J‍(∥θn-1∥+∥ψn-ψn-1Δt∥2+∥ηn-2ηn-1+ηn-2(Δt)2∥2+∥R1n-R1n-1Δt∥2)+a0Δt∑n=1J‍∥∇ςn-∇ςn-1Δt∥2+C(∥R1J∥2+∥ηJ-ηJ-1Δt∥2)+b02∥∇ϖJ∥2+CΔt∑n=1J‍∥∇ϖn-1∥2.
Use Gronwall lemma to obtain
(53)∥∇ϖJ∥2≤CΔt∑n=1J‍(∥θn-1∥+∥ψn-ψn-1Δt∥2+∥R1n-R1n-1Δt∥2+∥ηn-2ηn-1+ηn-2(Δt)2∥2)+C(∥R1J∥2+∥ηJ-ηJ-1Δt∥2).
Note that
(54)∥R1n-R1n-1Δt∥2=∥un-2un-1+un-2(Δt)2-ut(tn)-ut(tn-1)Δt∥2∥R1n-R1n-1Δt∥2=∥un-2un-1+un-2(Δt)2-utt(tΔ*)∥2∥R1n-R1n-1Δt∥2≤C[Δt∫tn-2tn‍∥∂4u∂t4∥ds+Δt∥∂3u∂t3∥L∞(L2)]2∥R1n-R1n-1Δt∥2≤C(Δt)2∫tn-2tn‍∥∂4u∂t4∥2ds+(Δt)2∥∂3u∂t3∥L∞(L2)2,wheretΔ*∈(tn-1,tn).
Substitute (44), (47) and (54) into (53) to obtain
(55)∥∇ϖJ∥2≤C(∥∂4u∂t4∥L2(L2)2∥ψ∥L∞(L2)2+∫t0tJ‍(∥ηt(s)∥2+∥ψt(s)∥2)ds+(Δt)2[∥∂3u∂t3∥L∞(L2)2+∥∂2u∂t2∥L2(L2)2+∥∂4u∂t4∥L2(L2)2]).
Use (55) and Poincaré inequality (19) to get
(56)∥ϖJ∥2≤C(∥∂4u∂t4∥L2(L2)2∥ψ∥L∞(L2)2+∫t0tJ‍(∥ηt(s)∥2+∥ψt(s)∥2)ds+(Δt)2[∥∂3u∂t3∥L∞(L2)2+∥∂2u∂t2∥L2(L2)2+∥∂4u∂t4∥L2(L2)2]).
Combining Lemmas 6-7, (55), (56) and using the triangle inequality, we complete the proof.


## 5. Numerical Results

In this section, we would like to give some numerical results for the coupling method of EMFE method and FE method proposed and analyzed in this paper. We consider the following initial-boundary value problem of fourth-order parabolic system:(57)ut+∇·(b∇(∇·(a(t)∇u)))=f(x,t),(x,t)∈Ω×J,u(x,t)=Δu(x,t)=0,(x,t)∈∂Ω×J¯,u(x,0)=sin(πx1)sin(πx2),x=(x1,x2)∈Ω,
where Ω=[0,1]×[0,1], J=(0,1], a(t)=1+t2,  b=1, and f(x,t) is chosen so that the exact solution for the scalar unknown function is
(58)u(x,t)=e-2tsin(πx1)sin(πx2),
its exact gradient is
(59)λ=∇u=e-2tπ[cos(πx1)sin(πx2)sin(πx1)cos(πx2)],
its exact flux function is
(60)σ=-aλ=e-2t(1+t2)π[cos(πx1)sin(πx2)sin(πx1)cos(πx2)],
and the diffusion term is
(61)γ=∇·σ=2π2e-2t(1+t2)sin(πx1)sin(πx2).Dividing the domainΩ into the triangulations of mesh size h uniformly and choosing the piecewise linear space Vh with index k=1 and h=2Δt=2/8,2/16,2/32,2/64, we obtain the optimal a priori error estimates in L2 and H1-norm for the scalar unknown u in Table 1. Similarly, from Table 2, we can find that both ∥γ-γh∥L∞(H1(Ω)) and ∥γ-γh∥L∞(L2(Ω)) for the diffusion term γ are optimal, too. At the same time, we take the corresponding piecewise constant space Wh with index k=0 for the gradient λ and the flux σ and obtain some convergence results in (L2)2-norm with h=2Δt=2/8,2/16,2/32,2/64 in Table 3.Table 1
L
2 and H1-norms and convergence order for u.
(
h
,
Δ
t
)
∥
u
-
u
h
∥
L
∞
(
L
2
(
Ω
)
)
Order
∥
u
-
u
h
∥
L
∞
(
H
1
(
Ω
)
)
Order
(
2
/
8
,
1
/
8
)
2.8933
e
-
002
2.3766
e
-
001
(
2
/
16
,
1
/
16
)
8.2158
e
-
003
1.8162
1.2986
e
-
001
0.8719
(
2
/
32
,
1
/
32
)
2.1012
e
-
003
1.9672
6.8466
e
-
002
0.9235
(
2
/
64
,
1
/
64
)
5.3511
e
-
004
1.9733
3.5245
e
-
002
0.9580Table 2
L
2 and H1-norms and convergence order for γ.
(
h
,
Δ
t
)
∥
γ
-
γ
h
∥
L
∞
(
L
2
(
Ω
)
)
Order
∥
γ
-
γ
h
∥
L
∞
(
H
1
(
Ω
)
)
Order
(
2
/
8
,
1
/
8
)
3.4238
e
-
001
4.5457
e
+
000
(
2
/
16
,
1
/
16
)
9.2701
e
-
002
1.8849
2.5429
e
+
000
0.8381
(
2
/
32
,
1
/
32
)
2.3840
e
-
002
1.9592
1.3493
e
+
000
0.9143
(
2
/
64
,
1
/
64
)
5.8649
e
-
003
2.0232
6.9558
e
-
001
0.9559Table 3
(
L
2
)
2-norm and convergence order for λ and σ.
(
h
,
Δ
t
)
∥
σ
-
σ
h
∥
L
∞
(
(
L
2
(
Ω
)
)
2
)
Order
∥
λ
-
λ
h
∥
L
∞
(
(
L
2
(
Ω
)
)
2
)
Order
(
2
/
8
,
1
/
8
)
2.3958
e
-
001
2.3589
e
-
001
(
2
/
16
,
1
/
16
)
1.3011
e
-
001
0.8808
1.2960
e
-
001
0.8641
(
2
/
32
,
1
/
32
)
6.8501
e
-
002
0.9255
6.8434
e
-
002
0.9213
(
2
/
64
,
1
/
64
)
3.5250
e
-
002
0.9585
3.5242
e
-
002
0.9574Takingt=1 with h=2Δt=2/32, we show the surfaces of the exact solution u, σ, λ, and γ in Figures 1, 3, 5, and 7, respectively. Figures 2, 4, 6, and 8 show the surfaces of the numerical solution uh, σh, λh and γh at t=1 with h=2Δt=2/16, respectively.Figure 1
Surface for exact solutionu.Figure 2
Surface for numerical solutionuh.Figure 3
Surface for exact solutionσ=(σ1,σ2).Figure 4
Surface for numerical solutionσh=(σ1h,σ2h).Figure 5
Surface for exact solutionλ=(λ1,λ2).Figure 6
Surface for numerical solutionλh=(λ1h,λ2h).Figure 7
Surface for exact solutionγ.Figure 8
Surface for numerical solutionγh.From the convergence results in Tables1–3 and Figures 1–8, we can find that the numerical results verify our theoretical analysis.


## 6. Concluding Remarks

In this article, we study a new coupling method of new EMFE scheme [29] and FE scheme for fourth-order partial differential equation of parabolic type. The new EMFE method’s gradient belongs to the simple square integrable (L2(Ω))2 space; therefore, the regularity requirement on the gradient solution λ=∇u is reduced. We obtain the optimal priori error estimates in L2 and H1-norm for both the scalar unknown u and the diffusion term γ and the priori error estimates in (L2)2-norm for its gradient λ and its flux σ. Finally, a numerical example is provided to verify the efficiency of our methods. At the same time, we apply the new coupling method based on new EMFE method and FE method to solving the fourth-order eigenvalue problems, which will be shown in another article.In the near future, we will study the new expanded mixed finite element method for other partial differential equations, such as fourth-order wave equations. Moreover, we will apply some new techniques based on the combination of two-grid methods [36, 37] and new (expanded) mixed methods [29–32, 38] for solving fourth-order nonlinear elliptic equation and fourth-order nonlinear parabolic equations.



---

*Source: 787891-2013-04-18.xml*