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

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 (𝐿²(Ω))² space, which avoids the use of the classical H(div; Ω) space and reduces the regularity requirement on the gradient solution 𝝀=∇𝑢. 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 𝐿² and 𝐻¹-norm for both the scalar unknown 𝑢 and the diffusion term γ and a priori error estimates in (𝐿²)²-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.

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 𝐿²-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 𝐻¹-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 𝐶¹-conforming FE method for nonlinear fourth-order wave equation. In this paper, we consider the following fourth-order partial differential equations of parabolic type: 𝑢_𝑡(𝑏()∇((𝑎(𝑡)(),()𝑢()()()+∇⋅𝐱,𝑡∇⋅∇𝑢)))=𝑓𝐱,𝑡𝐱,𝑡∈Ω×𝐽,𝐱,𝑡=Δ𝑢𝐱,𝑡=0,𝐱,𝑡∈𝜕Ω×𝑢()𝐽,𝐱,0=𝑢₀(𝐱),𝐱∈Ω, (1) where Ω⊂𝑅² is a bounded convex polygonal domain with the Lipschitz continuous boundary 𝜕Ω and 𝐽=(0,𝑇] is the time interval with 0<𝑇<∞. The initial value 𝑢₀(𝐱) and 𝑓(𝐱,𝑡) are given functions; coefficients 𝑎=𝑎(𝑡) and 𝑏=𝑏(𝐱,𝑡) satisfy 𝐀₁∶0<𝑎₀≤𝑎(𝑡)≤𝑎₁<+∞,𝑎_𝑡(𝑡)≤𝑐₀,𝐀₂∶0<𝑏₀≤𝑏(𝐱,𝑡)≤𝑏₁<+∞,𝑏_𝑡(𝐱,𝑡)≤𝑐₀, (2) for some positive constants 𝑎₀,  𝑎₁,  𝑏₀,  𝑏₁, and 𝑐₀.

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], 𝐻¹-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 (𝐿²(Ω))² space instead of the classical 𝐇(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 𝛾=−∇⋅(𝑎(𝑡)∇𝑢) 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, 𝝀=∇𝑢 and 𝝈=−𝑎(𝑡)∇𝑢=−𝑎𝝀, 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 (𝐿²(Ω))² space, the regularity requirement on the gradient solution 𝝀=∇𝑢 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 𝐿² and 𝐻¹-norm for both the scalar unknown 𝑢 and the diffusion term 𝛾 and a priori error estimates in (𝐿²)²-norm for its gradient 𝝀 and its flux 𝝈. Finally, we provide some numerical results to verify our theoretical analysis.

Throughout this paper, 𝐶 will denote a generic positive constant which does not depend on the spatial mesh parameter ℎ or the time step Δ𝑡. At the same time, we denote the natural inner product in 𝐿²(Ω) or (𝐿²(Ω))² 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 𝛾=−∇⋅(𝑎(𝑡)∇𝑢) to obtain the following coupled system of second-order equations: 𝑢_𝑡𝑏()(),(𝑎(𝑡))𝑢()()−∇⋅𝐱,𝑡∇𝛾=𝑓𝐱,𝑡𝛾+∇⋅∇𝑢=0,𝐱,𝑡=𝛾𝐱,𝑡=0. (3a)(3b)(3c)For (3b), we introduce the two auxiliary variables, 𝝀=∇𝑢 and 𝝈=−𝑎(𝑡)∇𝑢=−𝑎𝝀, to obtain the following first-order system of equations: 𝑢_𝑡(),𝑢()()−∇⋅𝑏∇𝛾=𝑓𝐱,𝑡𝛾−∇⋅𝝈=0,𝝈+𝑎𝝀=0,𝝀−∇𝑢=0,𝐱,𝑡=𝛾𝐱,𝑡=0. (4a)(4b)(4c)(4d)(4e)Using Green’s formula, the new mixed weak formulation of (4a)–(4e) is to determine {𝑢,𝛾,𝝀,𝝈}∶[0,𝑇]↦𝐻¹₀(Ω)×𝐻¹₀(Ω)×(𝐿²(Ω))²×(𝐿²(Ω))² such that 𝑢_𝑡,𝑣+𝑏∇𝛾,∇𝑣=𝑓,𝑣,∀𝑣∈𝐻¹₀(Ω),𝛾,𝑔+𝝈,∇𝑔=0,∀𝑔∈𝐻¹₀(Ω),(𝝈,𝐳)+(𝑎𝝀,𝐳)=0,∀𝐳∈𝐿²(Ω)²,(𝝀,𝐰)−(∇𝑢,𝐰)=0,∀𝐰∈𝐿²(Ω)². (5a)(5b)(5c)(5d)Then, the semidiscrete coupled EMFE scheme-FE scheme for (5a)–(5d) is to find {𝑢_ℎ,𝛾_ℎ,𝝀_ℎ,𝝈_ℎ}∶[0,𝑇]↦𝑉_ℎ×𝑉_ℎ×𝐖_ℎ×𝐖_ℎ such that 𝑢_ℎ𝑡,𝑣_ℎ+𝑏∇𝛾_ℎ,∇𝑣_ℎ=𝑓,𝑣_ℎ,∀𝑣_ℎ∈𝑉_ℎ,𝛾_ℎ,𝑔_ℎ+𝝈_ℎ,∇𝑔_ℎ=0,∀𝑔_ℎ∈𝑉_ℎ,𝝈_ℎ,𝐳_ℎ+𝑎𝝀_ℎ,𝐳_ℎ=0,∀𝐳_ℎ∈𝐖_ℎ,𝝀_ℎ,𝐰_ℎ−∇𝑢_ℎ,𝐰_ℎ=0,∀𝐰_ℎ∈𝐖_ℎ. (6a)(6b)(6c)(6d)

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 (𝐿²(Ω))² space taking the place of the classical 𝐇(div;Ω) space. It is easy to see that our method reduces the regularity requirement on the gradient solution 𝝀=∇𝑢.

Remark 2. In the semidiscrete coupling scheme (6a)–(6d), the MFE space (𝑉_ℎ,𝐖_ℎ) based on the FE pair 𝑃₁−𝑃²₀ is chosen as follows: 𝑉_ℎ=𝑣_ℎ∈𝐶⁰(Ω)∩𝐻¹₀(Ω)∣𝑣_ℎ∈𝑃₁(𝐾),∀𝐾∈𝒦_ℎ,𝐖_ℎ=𝐰_ℎ=𝑤_1ℎ,𝑤_2ℎ∈𝐿²(Ω)²∣𝑤_1ℎ,𝑤_2ℎ∈𝑃₀(𝐾),∀𝐾∈𝒦_ℎ. (7) From [30, 31], we know that MFE space (𝑉_ℎ,𝐖_ℎ) satisfies the discrete LBB condition.

Remark 3. From (6a)–(6d), we know that (6a) is an FE scheme for 𝑢_ℎ,𝛾_ℎ,𝑣_ℎ∈𝑉_ℎ and that (6b)–(6d) is a new EMFE scheme for (𝛾_ℎ,𝝈_ℎ,𝝀_ℎ), (𝑔_ℎ,𝐳_ℎ,𝐰_ℎ)∈𝑉_ℎ×𝐖_ℎ×𝐖_ℎ. So the scheme (6a)–(6d) is called as the coupling system of new EMFE scheme and FE scheme.

Remark 4. The finite element spaces 𝑉_ℎ and 𝐖_ℎ for the new EMFE method are chosen as finite dimensional subspaces of 𝐻¹₀(Ω) and (𝐿²(Ω))², respectively. However, 𝑉_ℎ and 𝐖_ℎ for 𝐻¹-Galerkin MFE method in two space variables in [8] are finite dimensional subspaces of 𝐻¹₀(Ω) and 𝐇(div;Ω), respectively. Compared to the complex 𝐇(div;Ω) space used in the 𝐻¹-Galerkin MFE method, we use the weaker (𝐿²(Ω))² 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 𝒫_ℎ and the new expanded mixed elliptic projection operator ℛ_ℎ associated with the coupled equations.

Lemma 5. One now defines a linear projection operator 𝒫_ℎ∶𝐻¹₀(Ω)→𝑉_ℎ as 𝑏∇𝛾−𝒫_ℎ𝛾,∇𝑣_ℎ=0,∀𝑣_ℎ∈𝑉_ℎ, (8) and then 𝛾−𝒫_ℎ𝛾_𝐿2(Ω)+ℎ𝛾−𝒫_ℎ𝛾_𝐻1≤𝐶ℎ²𝛾_𝐻2,𝛾_𝑡−𝒫_ℎ𝛾_{𝑡𝐿2(Ω)}≤𝐶ℎ²𝛾_𝐻2+𝛾_𝑡𝐻2,𝛾_𝑡𝑡−𝒫_ℎ𝛾_{𝑡𝑡𝐿2(Ω)}≤𝐶ℎ²𝛾_𝐻2+𝛾_𝑡𝐻2+𝛾_{𝑡𝑡𝐻2}. (9)

Let (ℛ_ℎ𝑢,ℛ_ℎ𝝀,ℛ_ℎ𝝈)∶[0,𝑇]→𝑉_ℎ×𝐖_ℎ×𝐖_ℎ be given by the following new expanded mixed relations [29] 𝝈−ℛ_ℎ𝝈,∇𝑔_ℎ=0,∀𝑔_ℎ∈𝑉_ℎ,𝝈−ℛ_ℎ𝝈,𝐳_ℎ+𝑎𝝀−ℛ_ℎ𝝀,𝐳_ℎ=0,∀𝐳_ℎ∈𝐖_ℎ,𝝀−ℛ_ℎ𝝀,𝐰_ℎ−∇𝑢−ℛ_ℎ𝑢,𝐰_ℎ=0,∀𝐰_ℎ∈𝐖_ℎ. (10a)(10b)(10c)

Then the following two important lemmas based on the new expanded mixed projection (10a)–(10c) hold.

Lemma 6. There is a constant 𝐶>0 independent of the spatial mesh parameter ℎ such that 𝝀−ℛ_ℎ𝝀≤𝐶ℎ‖𝝀‖_(𝐻¹₎²+‖𝑢‖_𝐻²,𝝈−ℛ_ℎ𝝈≤𝐶ℎ‖𝑢‖_𝐻²+‖𝝈‖_(𝐻¹₎²+‖𝝀‖_(𝐻¹₎²,∇𝑢−ℛ_ℎ𝑢≤𝐶ℎ‖𝑢‖_𝐻²+‖𝝀‖_(𝐻¹₎². (11)

Lemma 7. There is a constant 𝐶>0 independent of the spatial mesh parameter ℎ such that 𝑢−ℛ_ℎ𝑢≤𝐶ℎ²‖𝑢‖_𝐻²,𝑢_𝑡−ℛ_ℎ𝑢_𝑡≤𝐶ℎ²𝑢_𝑡𝐻2,𝑢_𝑡𝑡−ℛ_ℎ𝑢_𝑡𝑡≤𝐶ℎ²𝑢_{𝑡𝑡𝐻2},𝑢−ℛ_ℎ𝑢₁≤𝐶ℎ‖𝑢‖_𝐻²+‖𝝀‖_(𝐻¹₎². (12)

In [29], we can obtain the detailed proof for Lemmas 6-7.

For a priori error estimates, we write the errors as 𝛾−𝛾_ℎ=𝛾−𝒫_ℎ𝛾+𝒫_ℎ𝛾−𝛾_ℎ=𝜓+𝜛;𝑢−𝑢_ℎ=𝑢−ℛ_ℎ𝑢+ℛ_ℎ𝑢−𝑢_ℎ=𝜂+𝜍;𝝀−𝝀_ℎ=𝝀−ℛ_ℎ𝝀+ℛ_ℎ𝝀−𝝀_ℎ=𝜹+𝜽;𝝈−𝝈_ℎ=𝝈−ℛ_ℎ𝝈+ℛ_ℎ𝝈−𝝈_ℎ=𝝆+𝝃. (13)

Using (5a), (5b), (5c), (5d), (6a), (6b), (6c), (6d), (8), (10a), (10b), and (10c) we can obtain the error equations 𝜍_𝑡,𝑣_ℎ+𝑏∇𝜛,∇𝑣_ℎ=−𝜂_𝑡,𝑣_ℎ,∀𝑣_ℎ∈𝑉_ℎ,𝜛,𝑔_ℎ+𝝃,∇𝑔_ℎ=−𝜓,𝑔_ℎ,∀𝑔_ℎ∈𝑉_ℎ,𝝃,𝐳_ℎ+𝑎𝜽,𝐳_ℎ=0,∀𝐳_ℎ∈𝐖_ℎ,𝜽,𝐰_ℎ−∇𝜍,𝐰_ℎ=0,∀𝐰_ℎ∈𝐖_ℎ. (14a)(14b)(14c)(14d)We will prove the error estimates for semidiscrete scheme.

Theorem 8. Suppose that 𝛾,𝛾_𝑡,𝑢_𝑡∈𝐿²(𝐻²(Ω)), 𝝀,𝝈∈𝐿^∞((𝐻¹(Ω))²), and 𝑢,𝛾∈𝐿^∞(𝐻²(Ω)); then there exists a constant 𝐶>0 independent of the spatial mesh parameter ℎ such that 𝑢−𝑢_ℎ≤𝐶ℎ²‖𝑢‖_𝐿^∞_(𝐻²₎+𝛾_{𝐿∞(𝐻2)}+^𝑡₀|||■|||²𝑑𝑠^1/2,𝑢−𝑢_ℎ1≤𝐶ℎ‖𝝀‖_{𝐿∞((𝐻1)}²₎+‖𝑢‖_𝐿^∞_(𝐻²₎+𝛾_{𝐿∞(𝐻2)}+^𝑡₀|||■|||²𝑑𝑠^1/2,𝝀−𝝀_ℎ≤𝐶ℎ‖𝝀‖_𝐿^∞_((𝐻¹₎²₎+𝛾_{𝐿∞(𝐻2)}+^𝑡₀|||■|||²𝑑𝑠^1/2,𝝈−𝝈_ℎ≤𝐶ℎ‖𝑢‖_𝐿^∞_(𝐻²₎+‖𝝈‖_𝐿^∞_((𝐻¹₎²₎+‖𝝀‖_𝐿^∞_((𝐻¹₎²₎+𝛾_{𝐿∞(𝐻2)}+^𝑡₀|||■|||²𝑑𝑠^1/2, (15) where |||■|||²≜‖𝑢_𝑡‖²_𝐻2+‖𝛾‖²_𝐻2+‖𝛾_𝑡‖²_𝐻2.

Proof. Choose 𝑣_ℎ=𝜛, 𝑔_ℎ=𝜍_𝑡, 𝐳_ℎ=∇𝜍_𝑡,  and  𝐰_ℎ=∇𝜍_𝑡 in (14a)–(14d) to obtain 𝜍_𝑡,𝜛+𝑏^1/2∇𝜛²=−𝜂_𝑡,−,𝜛𝜛,𝜍_𝑡−𝝃,∇𝜍_𝑡=𝜓,𝜍_𝑡,𝝃,∇𝜍_𝑡+𝑎𝜽,∇𝜍_𝑡−=0,𝑎𝜽,∇𝜍_𝑡+12𝑑𝑎𝑑𝑡^1/2∇𝜍²=12𝑎_𝑡∇𝜍,∇𝜍. (16a)(16b)(16c)(16d)Adding the above four equations and using Cauchy-Schwarz inequality and Young’s inequality, we have 12𝑑𝑎𝑑𝑡^1/2∇𝜍²+𝑏^1/2∇𝜛²=−𝜂_𝑡,𝜛+𝜓,𝜍_𝑡+12𝑎_𝑡𝜂∇𝜍,∇𝜍=−_𝑡+𝑑,𝜛−𝜓𝑑𝑡𝜓,𝜍_𝑡+1,𝜍2𝑎_𝑡≤𝑑∇𝜍,∇𝜍+𝜂𝑑𝑡𝜓,𝜍_𝑡‖𝜛‖+𝐶‖𝜍‖²₁+𝜓_𝑡². (17) Integrate with respect to time from 0 to 𝑡 to obtain 𝑎₀‖∇𝜍‖²+2𝑏₀^𝑡₀‖∇𝜛‖²𝑑𝑠≤2𝜓,𝜍+2^𝑡₀𝜂_𝑡‖𝜛‖𝑑𝑠+𝐶^𝑡₀‖𝜍‖²₁+𝜓_𝑡²𝑑𝑠. (18) Noting that 𝜍,𝜛∈𝐻¹₀(Ω), we use Poincaré inequality to have ‖𝜍‖₁≤𝐶₁‖∇𝜍‖,‖𝜛‖≤𝐶₂‖∇𝜛‖. (19) Using (18), (19), Cauchy-Schwarz inequality and Young inequality, we have 𝑎₀‖∇𝜍‖²+2𝑏₀^𝑡₀‖∇𝜛‖²𝜓𝜍‖𝑑𝑠≤2‖+2^𝑡₀𝜂_𝑡‖𝜛‖𝑑𝑠+𝐶^𝑡₀‖𝜍‖²₁+𝜓_𝑡²𝑑𝑠≤2𝐶₁𝜓‖∇𝜍‖+2𝐶₂^𝑡₀𝜂_𝑡‖∇𝜛‖𝑑𝑠+𝐶^𝑡₀𝐶²₁‖∇𝜍‖²+𝜓_𝑡²≤𝑑𝑠8𝐶²₁𝑎₀𝜓²+𝑎₀2‖∇𝜍‖²+4𝐶²₂𝑏₀^𝑡₀𝜂_𝑡²𝑑𝑠+𝑏₀^𝑡₀‖∇𝜛‖²𝑑𝑠+𝐶^𝑡₀𝐶²₁‖∇𝜍‖²+𝜓_𝑡²𝑑𝑠. (20) Using Gronwall’s lemma, we get ‖∇𝜍‖²+^𝑡₀‖∇𝜛‖²𝜓𝑑𝑠≤𝐶²+^𝑡₀𝜂_𝑡²+𝜓²+𝜓_𝑡²𝑑𝑠. (21) Substitute (19) into (21) to have ‖𝜍‖²+^𝑡₀‖𝜛‖²𝜓𝑑𝑠≤𝐶²+^𝑡₀𝜂_𝑡²+𝜓²+𝜓_𝑡²𝑑𝑠. (22) Taking 𝐰_ℎ=𝜽 in (14d) and using (22), we get ‖𝜽‖²≤𝐶𝜓²+^𝑡₀𝜂_𝑡²+𝜓²+𝜓_𝑡²𝑑𝑠. (23) Choose 𝐳_ℎ=𝝃 in (14c) and use (22) to obtain ‖𝝃‖²≤𝐶𝜓²+^𝑡₀𝜂_𝑡²+𝜓²+𝜓_𝑡²𝑑𝑠. (24) Combining Lemmas 5–7, (21)–(24) and using the triangle inequality, we get the error estimates for Theorem 8.

Remark 9. From Theorem 8, we can see that 𝑢−𝑢_{ℎ𝐿∞(𝐿2)}+ℎ𝑢−𝑢_{ℎ𝐿∞(𝐻1)}=𝑂ℎ²; (25) that is to say, a priori error estimates in 𝐿² and 𝐻¹-norms for the scalar unknown 𝑢 are optimal.

Theorem 10. Suppose that 𝜕^𝑗𝑢/𝜕𝑡^𝑗,𝜕^𝑗𝛾/𝜕𝑡^𝑗∈𝐿²(𝐻²(Ω)) and 𝛾,𝛾_𝑡,𝑢_𝑡∈𝐿^∞(𝐻²(Ω)); then there exists a constant 𝐶>0 independent of the spatial mesh parameter ℎ such that 𝛾−𝛾_ℎ+ℎ𝛾−𝛾_ℎ1≤𝐶ℎ²𝑢_{𝑡𝐿∞(𝐻2)}+𝛾_{𝐿∞(𝐻2)}+𝛾_{𝑡𝐿∞(𝐻2)}+^𝑡₀²_𝑗=0𝜕^𝑗𝑢𝜕𝑡^𝑗2_𝐻2+𝜕^𝑗𝛾𝜕𝑡^𝑗2_𝐻2𝑑𝑠^1/2. (26)

Proof. Differentiating (14b), (14c), and (14d) with respect to time 𝑡, respectively, we get 𝜍_𝑡,𝑣_ℎ+𝑏∇𝜛,∇𝑣_ℎ=−𝜂_𝑡,𝑣_ℎ,𝜛_𝑡,𝑔_ℎ+𝝃_𝑡,∇𝑔_ℎ=−𝜓_𝑡,𝑔_ℎ,𝝃_𝑡,𝐳_ℎ+𝑎_𝑡𝜽,𝐳_ℎ+𝑎𝜽_𝑡,𝐳_ℎ𝜽=0,_𝑡,𝐰_ℎ−∇𝜍_𝑡,𝐰_ℎ=0. (27a)(27b)(27c)(27d)Choose 𝑣_ℎ=𝜛_𝑡, 𝑔_ℎ=𝜍_𝑡, 𝐳_ℎ=∇𝜍_𝑡,  and  𝐰_ℎ=∇𝜍_𝑡 in (27a)–(27d) to get 𝜍_𝑡,𝜛_𝑡+12𝑑𝑏𝑑𝑡^1/2∇𝜛²=−𝜂_𝑡,𝜛_𝑡+12𝑏_𝑡,−𝜛∇𝜛,∇𝜛_𝑡,𝜍_𝑡−𝝃_𝑡,∇𝜍_𝑡=𝜓_𝑡,𝜍_𝑡,𝝃_𝑡,∇𝜍_𝑡+𝑎_𝑡𝜽,∇𝜍_𝑡+𝑎𝜽_𝑡,∇𝜍_𝑡−=0,𝑎𝜽_𝑡,∇𝜍_𝑡+𝑎^1/2∇𝜍_𝑡²=0. (28a)(28b)(28c)(28d)Adding the four equations for (28a)–(28d), we have 12𝑑𝑏𝑑𝑡^1/2∇𝜛²+𝑎^1/2∇𝜍_𝑡²=−𝜂_𝑡,𝜛_𝑡+𝜓_𝑡,𝜍_𝑡−𝑎_𝑡𝜽,∇𝜍_𝑡+12𝑏_𝑡𝑑∇𝜛,∇𝜛=−𝜂𝑑𝑡^𝑡,𝜛+𝜂_𝑡𝑡+𝑑,𝜛𝜓𝑑𝑡^𝑡−𝜓,𝜍_𝑡𝑡,𝜍−𝑎_𝑡𝜽,∇𝜍_𝑡+12𝑏_𝑡∇𝜛,∇𝜛. (29) Integrate with respect to time from 0 to 𝑡 and use (19), Cauchy-Schwarz inequality, and Young inequality to get 𝑏₀‖∇𝜛‖²+2𝑎₀^𝑡₀∇𝜍_𝑡²𝜂𝑑𝑠≤−2_𝑡,𝜛+2𝐶₂^𝑡₀𝜂_𝑡𝑡‖∇𝜛‖𝑑𝑠+2𝜓_𝑡,𝜍+𝐶^𝑡₀𝜓_𝑡𝑡²+‖𝜍‖²+‖𝜽‖²𝑑𝑠+𝑎₀^𝑡₀∇𝜍_𝑡²≤𝑑𝑠8𝐶²₂𝑏₀𝜂_𝑡²+𝑏₀2‖∇𝜛‖²+𝜓_𝑡²+‖𝜍‖²+𝑎₀^𝑡₀∇𝜍_𝑡²𝑑𝑠+𝐶^𝑡₀𝜓_𝑡𝑡²+‖𝜍‖²+‖𝜽‖²+𝜂_𝑡𝑡²+‖∇𝜛‖²𝑑𝑠. (30) Substituting (22) and (23) into (30) and using Gronwall lemma, we get ‖∇𝜛‖²+^𝑡₀∇𝜍_𝑡²𝜂𝑑𝑠≤𝐶_𝑡²+𝜓²+𝜓_𝑡²+𝐶^𝑡₀𝜂²+𝜂_𝑡²+𝜂_𝑡𝑡²+𝜓²+𝜓_𝑡²+𝜓_𝑡𝑡²𝑑𝑠. (31) Use (19) to obtain ‖𝜛‖²≤𝐶𝜂_𝑡²+𝜓²+𝜓_𝑡²+𝐶^𝑡₀𝜂²+𝜂_𝑡²+𝜂_𝑡𝑡²+𝜓²+𝜓_𝑡²+𝜓_𝑡𝑡²𝑑𝑠. (32) Combining Lemmas 5–7, (31), and (32) and using the triangle inequality, we obtain the error estimates for Theorem 10.

Remark 11. From Theorem 10, we can see that 𝛾−𝛾_{ℎ𝐿∞(𝐿2)}+ℎ𝛾−𝛾_{ℎ𝐿∞(𝐻1)}=𝑂ℎ²; (33) that is to say, a priori error estimates in 𝐿² and 𝐻¹-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. Let 0=𝑡₀<𝑡₁<𝑡₂<⋯<𝑡_𝑀=𝑇 be a given partition of the time interval [0,𝑇] with step length Δ𝑡=𝑇/𝑀 and nodes 𝑡_𝑛=𝑛Δ𝑡, for some positive integer 𝑀. For a smooth function 𝜙 on [0,𝑇], define 𝜙^𝑛=𝜙(𝑡_𝑛).

For (5a), (5b), (5c), and (5d), we have the following equivalent formulation:𝑢^𝑛−𝑢^𝑛−1+𝑏Δ𝑡,𝑣^𝑛∇𝛾^𝑛,∇𝑣=𝑓^𝑛+𝑅^𝑛₁,𝑣,∀𝑣∈𝐻¹₀(Ω),𝛾^𝑛,𝑔+𝝈^𝑛,∇𝑔=0,∀𝑔∈𝐻¹₀(Ω),𝝈^𝑛,𝐳+𝑎^𝑛𝝀^𝑛,𝐳=0,∀𝐳∈𝐿²(Ω)²,𝝀^𝑛,𝐰−∇𝑢^𝑛,𝐰=0,∀𝐰∈𝐿²(Ω)², (34a)(34b)(34c)(34d)where 𝑅^𝑛₁=𝑢^𝑛−𝑢^𝑛−1Δ𝑡−𝑢_𝑡𝑡_𝑛=1Δ𝑡^𝑡𝑛_𝑡𝑛−1𝑡_𝑛−1−𝑠𝑢_𝑡𝑡𝑑𝑠. (35) Now a fully discrete procedure is to find (𝑢^𝑛_ℎ,𝛾^𝑛_ℎ,𝝀^𝑛_ℎ,𝝈^𝑛_ℎ)∈𝑉_ℎ×𝑉_ℎ×𝐖_ℎ×𝐖_ℎ(𝑛=0,1,…,𝑀) such that 𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1Δ𝑡,𝑣_ℎ+𝑏^𝑛∇𝛾^𝑛_ℎ,∇𝑣_ℎ=𝑓^𝑛,𝑣_ℎ,∀𝑣_ℎ∈𝑉_ℎ,𝛾^𝑛_ℎ,𝑔_ℎ+𝝈^𝑛_ℎ,∇𝑔_ℎ=0,∀𝑔_ℎ∈𝑉_ℎ,𝝈^𝑛_ℎ,𝐳_ℎ+𝑎^𝑛𝝀^𝑛_ℎ,𝐳_ℎ=0,∀𝐳_ℎ∈𝐖_ℎ,𝝀^𝑛_ℎ,𝐰_ℎ−∇𝑢^𝑛_ℎ,𝐰_ℎ=0,∀𝐰_ℎ∈𝐖_ℎ. (36a)(36b)(36c)(36d)For deriving the fully discrete error estimates, we now decompose the errors as 𝛾^𝑛−𝛾^𝑛_ℎ=𝛾^𝑛−𝒫_ℎ𝛾^𝑛+𝒫_ℎ𝛾^𝑛−𝛾^𝑛_ℎ=𝜓^𝑛+𝜛^𝑛;𝑢^𝑛−𝑢^𝑛_ℎ=𝑢^𝑛−ℛ_ℎ𝑢^𝑛+ℛ_ℎ𝑢^𝑛−𝑢^𝑛_ℎ=𝜂^𝑛+𝜍^𝑛;𝝀^𝑛−𝝀^𝑛_ℎ=𝝀^𝑛−ℛ_ℎ𝝀^𝑛+ℛ_ℎ𝝀^𝑛−𝝀^𝑛_ℎ=𝜹^𝑛+𝜽^𝑛;𝝈^𝑛−𝝈^𝑛_ℎ=𝝈^𝑛−ℛ_ℎ𝝈^𝑛+ℛ_ℎ𝝈^𝑛−𝝈^𝑛_ℎ=𝝆^𝑛+𝝃^𝑛. (37) We will prove the theorem for the fully discrete error estimates.

Theorem 12. Suppose that 𝛾,𝛾_𝑡,𝑢_𝑡∈𝐿²(𝐻²(Ω)), 𝝀,𝝈∈𝐿^∞((𝐻¹(Ω))²), 𝛾∈𝐿^∞(𝐻²(Ω)), and 𝑢_𝑡𝑡∈𝐿²(𝐿²(Ω)); then there exists a constant 𝐶>0 independent of the spatial mesh parameter ℎ and the time step Δ𝑡 such that 𝑢^𝐽−𝑢^𝐽_ℎ≤𝐶ℎ²‖𝑢‖_𝐿^∞_(𝐻²₎+|||◆|||+𝐶Δ𝑡𝑢_{𝑡𝑡𝐿2(𝐿2)},𝑢^𝐽−𝑢^𝐽_ℎ1≤𝐶ℎ‖𝑢‖_𝐿^∞_(𝐻²₎+‖𝝀‖_𝐿^∞_((𝐻¹₎²₎+◆𝑢||||||+𝐶Δ𝑡_{𝑡𝑡𝐿2(𝐿2)},𝝀^𝐽−𝝀^𝐽_ℎ≤𝐶ℎ‖𝝀‖_𝐿^∞_((𝐻¹₎²₎+|||◆|||+𝐶Δ𝑡𝑢_{𝑡𝑡𝐿2(𝐿2)},𝝈^𝐽−𝝈^𝐽_ℎ≤𝐶ℎ‖𝝈‖_𝐿^∞_((𝐻¹₎²₎+|||◆|||+𝐶Δ𝑡𝑢_{𝑡𝑡𝐿2(𝐿2)}, (38) where |||◆|||≜‖𝑢_𝑡‖_𝐿2(𝐻2)+‖𝛾‖_{𝐿∞(𝐻2)}+‖𝛾‖_𝐿2(𝐻2)+‖𝛾_𝑡‖_𝐿2(𝐻2).

Proof. Using (8), (10a), (10b), (10c), (34a), (34b), (34c), (34d), (36a), (36b), (36c), and (36d) at 𝑡=𝑡_𝑛, we obtain the error equations𝜍^𝑛−𝜍^𝑛−1Δ𝑡,𝑣_ℎ+𝑏^𝑛∇𝜛^𝑛,∇𝑣_ℎ𝜂=−^𝑛−𝜂^𝑛−1Δ𝑡,𝑣_ℎ+𝑅^𝑛₁,𝑣_ℎ,∀𝑣_ℎ∈𝑉_ℎ,𝜛^𝑛,𝑔_ℎ+𝝃^𝑛,∇𝑔_ℎ=−𝜓^𝑛,𝑔_ℎ,∀𝑔_ℎ∈𝑉_ℎ,𝝃^𝑛,𝐳_ℎ+𝑎^𝑛𝜽^𝑛,𝐳_ℎ=0,∀𝐳_ℎ∈𝐖_ℎ,𝜽^𝑛,𝐰_ℎ−∇𝜍^𝑛,𝐰_ℎ=0,∀𝐰_ℎ∈𝐖_ℎ. (39a)(39b)(39c)(39d)Take 𝑣_ℎ=𝜛^𝑛 in (39a), 𝑔_ℎ=(𝜍^𝑛−𝜍^𝑛−1)/Δ𝑡 in (39b), 𝐳_ℎ=(∇𝜍^𝑛−∇𝜍^𝑛−1)/Δ𝑡 in (39c), and 𝐰_ℎ=(∇𝜍^𝑛−∇𝜍^𝑛−1)/Δ𝑡 in (39d) to get 𝑏^𝑛1/2∇𝜛^𝑛2𝜍=−^𝑛−𝜍^𝑛−1Δ𝑡,𝜛^𝑛−𝜂^𝑛−𝜂^𝑛−1Δ𝑡,𝜛^𝑛+𝑅^𝑛₁,𝜛^𝑛=𝝃^𝑛,∇𝜍^𝑛−∇𝜍^𝑛−1−𝜂Δ𝑡^𝑛−𝜂^𝑛−1Δ𝑡,𝜛^𝑛+𝑅^𝑛₁,𝜛^𝑛+𝜓^𝑛,𝜍^𝑛−𝜍^𝑛−1𝑎Δ𝑡=−^𝑛𝜽^𝑛,∇𝜍^𝑛−∇𝜍^𝑛−1−𝜂Δ𝑡^𝑛−𝜂^𝑛−1Δ𝑡,𝜛^𝑛+𝑅^𝑛₁,𝜛^𝑛+𝜓^𝑛,𝜍^𝑛−𝜍^𝑛−1𝑎Δ𝑡=−^𝑛∇𝜍^𝑛−∇𝜍^𝑛−1Δ𝑡,∇𝜍^𝑛−𝜂^𝑛−𝜂^𝑛−1Δ𝑡,𝜛^𝑛+𝑅^𝑛₁,𝜛^𝑛+𝜓^𝑛,𝜍^𝑛−𝜍^𝑛−11Δ𝑡=−𝑎2Δ𝑡^𝑛1/2∇𝜍^𝑛2−𝑎^𝑛−11/2∇𝜍^𝑛−12−1𝑎2Δ𝑡^𝑛∇𝜍^𝑛−∇𝜍^𝑛−12+12𝑎^𝑛−𝑎^𝑛−1Δ𝑡∇𝜍^𝑛−1,∇𝜍^𝑛−1−𝜂^𝑛−𝜂^𝑛−1Δ𝑡,𝜛^𝑛+𝑅^𝑛₁,𝜛^𝑛+𝜓^𝑛,𝜍^𝑛−𝜍^𝑛−1.Δ𝑡 (40) By (40), we have 𝑏^𝑛1/2∇𝜛^𝑛2+1𝑎2Δ𝑡^𝑛∇𝜍^𝑛−∇𝜍^𝑛−12+1𝑎2Δ𝑡^𝑛1/2∇𝜍^𝑛2−𝑎^𝑛−11/2∇𝜍^𝑛−12=12𝑎^𝑛−𝑎^𝑛−1Δ𝑡∇𝜍^𝑛−1,∇𝜍^𝑛−1−𝜂^𝑛−𝜂^𝑛−1Δ𝑡,𝜛^𝑛+𝑅^𝑛₁,𝜛^𝑛+𝜓^𝑛,𝜍^𝑛−𝜍^𝑛−1=1Δ𝑡21Δ𝑡^𝑡𝑛_𝑡𝑛−1𝑎_𝑡(𝑠)𝑑𝑠∇𝜍^𝑛−1,∇𝜍^𝑛−1−𝜂^𝑛−𝜂^𝑛−1Δ𝑡,𝜛^𝑛+𝑅^𝑛₁,𝜛^𝑛+𝜓^𝑛,𝜍^𝑛−𝜓^𝑛−1,𝜍^𝑛−1−𝜍Δ𝑡_𝑛−1,𝜓^𝑛−𝜓^𝑛−1.Δ𝑡 (41) Multiplying by 2Δ𝑡, summing (41) from 𝑛=1 to 𝐽, and using (19), the resulting equation becomes 𝑎₀∇𝜍^𝐽2+2𝑏₀Δ𝑡^𝐽_𝑛=1∇𝜛^𝑛2≤2Δ𝑡^𝐽_𝑛=1𝜂^𝑛−𝜂^𝑛−1+𝑅Δ𝑡^𝑛1𝜛^𝑛+sup_{𝑡𝑡∈𝑛−1,𝑡𝑛}𝑎_𝑡(𝑡)Δ𝑡^𝐽−1_𝑛=1∇𝜍^𝑛2+Δ𝑡^𝐽_𝑛=1𝜍^𝑛−12+𝜓^𝑛−𝜓^𝑛−1Δ𝑡²+𝜓^𝐽𝜍^𝐽≤𝐶Δ𝑡^𝐽_𝑛=1𝜂^𝑛−𝜂^𝑛−1Δ𝑡²+𝑅^𝑛₁²+𝑏₀Δ𝑡^𝐽_𝑛=1∇𝜛^𝑛2+𝐶Δ𝑡^𝐽_𝑛=1∇𝜍^𝑛−12+𝜓^𝑛−𝜓^𝑛−1Δ𝑡²+𝐶𝜓^𝐽2+𝑎₀2∇𝜍^𝐽2. (42) Use (42) and Gronwall lemma to get ∇𝜍^𝐽2+Δ𝑡^𝐽_𝑛=1∇𝜛^𝑛2≤𝐶𝜓^𝐽2+𝐶Δ𝑡^𝐽_𝑛=1𝜂^𝑛−𝜂^𝑛−1Δ𝑡²+𝜓^𝑛−𝜓^𝑛−1Δ𝑡²+𝑅^𝑛₁². (43) Note that 𝑅^𝑛₁²≤𝐶Δ𝑡^𝑡𝑛_𝑡𝑛−1𝑢_𝑡𝑡²𝜂𝑑𝑠,_𝑛−𝜂^𝑛−1Δ𝑡²≤1Δ𝑡^𝑡𝑛_𝑡𝑛−1𝜂_𝑡(𝑠)²𝜓𝑑𝑠,_𝑛−𝜓^𝑛−1Δ𝑡²≤1Δ𝑡^𝑡𝑛_𝑡𝑛−1𝜓_𝑡(𝑠)²𝑑𝑠. (44) Substitute (44) into (43) to get ∇𝜍^𝐽2+Δ𝑡^𝐽_𝑛=1∇𝜛^𝑛2≤𝐶𝜓^𝐽2+^𝑡𝐽_𝑡0𝜂_𝑡(𝑠)²+𝜓_𝑡(𝑠)²+()𝑑𝑠Δ𝑡^2𝑡𝐽_𝑡0𝑢_𝑡𝑡²𝑑𝑠. (45) Combining (19) and (45), we obtain 𝜍^𝐽2+Δ𝑡^𝐽_𝑛=1𝜛^𝑛2≤𝐶𝜓^𝐽2+^𝑡𝐽_𝑡0𝜂_𝑡(𝑠)²+𝜓_𝑡(𝑠)²+()𝑑𝑠Δ𝑡^2𝑡𝐽_𝑡0𝑢_𝑡𝑡²𝑑𝑠. (46) Taking 𝐰_ℎ=𝜽^𝑛 in (39d) and 𝐳_ℎ=𝝃^𝑛 in (39c) and using (45), we have 𝜽^𝐽2+𝝃^𝐽2≤𝐶𝜓^𝐽2+^𝑡𝐽_𝑡0𝜂_𝑡(𝑠)²+𝜓_𝑡(𝑠)²+()𝑑𝑠Δ𝑡^2𝑡𝐽_𝑡0𝑢_𝑡𝑡²𝑑𝑠. (47) 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 in 𝐿² and 𝐻¹-norm for 𝛾.

Theorem 13. Suppose that 𝛾,𝛾_𝑡,𝑢_𝑡∈𝐿²(𝐻²(Ω)),  𝜕²𝑢/𝜕𝑡²,  𝜕⁴𝑢/𝜕𝑡⁴∈𝐿²(𝐿²(Ω)),  𝛾∈𝐿^∞(𝐻²(Ω)), and 𝜕³𝑢/𝜕𝑡³∈𝐿^∞(𝐿²(Ω)); then there exists a constant 𝐶>0 independent of the spatial mesh parameter ℎ and time discretization parameter Δ𝑡 such that 𝛾^𝐽−𝛾^𝐽_ℎ𝑗≤𝐶ℎ^2−𝑗𝛾_{𝐿∞(𝐻2)}+𝛾_𝐿2(𝐻2)+𝑢_{𝑡𝐿2(𝐻2)}+𝛾_{𝑡𝐿2(𝐻2)}𝜕+Δ𝑡²𝑢𝜕𝑡²_𝐿2(𝐿2)+𝜕³𝑢𝜕𝑡³_{𝐿∞(𝐿2)}+𝜕⁴𝑢𝜕𝑡⁴_𝐿2(𝐿2),𝑗=0,1. (48)

Proof. Use (39a)–(39d) to get𝜍^𝑛−𝜍^𝑛−1Δ𝑡,𝑣_ℎ+𝑏^𝑛∇𝜛^𝑛,∇𝑣_ℎ𝜂=−^𝑛−𝜂^𝑛−1Δ𝑡,𝑣_ℎ+𝑅^𝑛₁,𝑣_ℎ,∀𝑣_ℎ∈𝑉_ℎ,𝜛^𝑛−𝜛^𝑛−1Δ𝑡,𝑔_ℎ+𝝃^𝑛−𝝃^𝑛−1Δ𝑡,∇𝑔_ℎ𝜓=−^𝑛−𝜓^𝑛−1Δ𝑡,𝑔_ℎ,∀𝑔_ℎ∈𝑉_ℎ,𝝃^𝑛−𝝃^𝑛−1Δ𝑡,𝐳_ℎ+𝑎^𝑛𝜽^𝑛−𝑎^𝑛−1𝜽^𝑛−1Δ𝑡,𝐳_ℎ=0,∀𝐳_ℎ∈𝐖_ℎ,𝜽^𝑛−𝜽^𝑛−1Δ𝑡,𝐰_ℎ−∇𝜍^𝑛−∇𝜍^𝑛−1Δ𝑡,𝐰_ℎ=0,∀𝐰_ℎ∈𝐖_ℎ. (49a)(49b)(49c)(49d)Taking 𝑣_ℎ=(𝜛^𝑛−𝜛^𝑛−1)/Δ𝑡,  𝑔_ℎ=(𝜍^𝑛−𝜍^𝑛−1)/Δ𝑡,  𝐳_ℎ=(∇𝜍^𝑛−∇𝜍^𝑛−1)/Δ𝑡, and 𝐰_ℎ=(∇𝜍^𝑛−∇𝜍^𝑛−1)/Δ𝑡, respectively, in (49a)–(49d), we get 1𝑏2Δ𝑡^𝑛1/2∇𝜛^𝑛2−𝑏^𝑛−11/2∇𝜛^𝑛−12+1𝑏2Δ𝑡^𝑛(∇𝜛^𝑛−∇𝜛^𝑛−1)²=𝑏^𝑛∇𝜛^𝑛𝜛,∇^𝑛−𝜛^𝑛−1+1Δ𝑡2𝑏^𝑛−𝑏^𝑛−1Δ𝑡∇𝜛^𝑛−1,∇𝜛^𝑛−1𝜍=−^𝑛−𝜍^𝑛−1,𝜛Δ𝑡^𝑛−𝜛^𝑛−1−𝜂Δ𝑡_𝑛−𝜂^𝑛−1,𝜛Δ𝑡^𝑛−𝜛^𝑛−1+𝑅Δ𝑡_𝑛1,𝜛^𝑛−𝜛^𝑛−1+1Δ𝑡2𝑏^𝑛−𝑏^𝑛−1Δ𝑡∇𝜛^𝑛−1,∇𝜛^𝑛−1=𝝃^𝑛−𝝃^𝑛−1,Δ𝑡∇𝜍^𝑛−∇𝜍^𝑛−1+𝜓Δ𝑡^𝑛−𝜓^𝑛−1,𝜍Δ𝑡^𝑛−𝜍^𝑛−1+𝑅Δ𝑡_𝑛1,𝜛^𝑛−𝜛^𝑛−1−𝜂Δ𝑡^𝑛−𝜂^𝑛−1,𝜛Δ𝑡^𝑛−𝜛^𝑛−1+1Δ𝑡2𝑏^𝑛−𝑏^𝑛−1Δ𝑡∇𝜛^𝑛−1,∇𝜛^𝑛−1𝑎=−^𝑛𝜽^𝑛−𝜽^𝑛−1,Δ𝑡∇𝜍^𝑛−∇𝜍^𝑛−1−𝑎Δ𝑡_𝑛−𝑎^𝑛−1𝜽Δ𝑡^𝑛−1,∇𝜍^𝑛−∇𝜍^𝑛−1+𝜓Δ𝑡_𝑛−𝜓^𝑛−1,𝜍Δ𝑡^𝑛−𝜍^𝑛−1+𝑅Δ𝑡_𝑛1,𝜛^𝑛−𝜛^𝑛−1−𝜂Δ𝑡^𝑛−𝜂^𝑛−1,𝜛Δ𝑡^𝑛−𝜛^𝑛−1+1Δ𝑡2𝑏^𝑛−𝑏^𝑛−1Δ𝑡∇𝜛^𝑛−1,∇𝜛^𝑛−1=−𝑎^𝑛1/2∇𝜍^𝑛−∇𝜍^𝑛−1Δ𝑡²−𝑎^𝑛−𝑎^𝑛−1𝜽Δ𝑡^𝑛−1,∇𝜍^𝑛−∇𝜍^𝑛−1+𝜓Δ𝑡_𝑛−𝜓^𝑛−1,𝜍Δ𝑡^𝑛−𝜍^𝑛−1+𝑅Δ𝑡_𝑛1,𝜛^𝑛−𝜛^𝑛−1−𝜂Δ𝑡_𝑛−𝜂^𝑛−1,𝜛Δ𝑡^𝑛−𝜛^𝑛−1+1Δ𝑡2𝑏^𝑛−𝑏^𝑛−1Δ𝑡∇𝜛^𝑛−1,∇𝜛^𝑛−1. (50) By using (50), we obtain 1𝑏2Δ𝑡^𝑛1/2∇𝜛^𝑛2−𝑏^𝑛−11/2∇𝜛^𝑛−12+1𝑏2Δ𝑡^𝑛∇𝜛^𝑛−∇𝜛^𝑛−12+𝑎^𝑛1/2∇𝜍^𝑛−∇𝜍^𝑛−1Δ𝑡²𝑎=−^𝑛−𝑎^𝑛−1𝜽Δ𝑡^𝑛−1𝜍,∇^𝑛−𝜍^𝑛−1+𝜓Δ𝑡_𝑛−𝜓^𝑛−1,𝜍Δ𝑡^𝑛−𝜍^𝑛−1+𝜛Δ𝑡_𝑛,𝑅^𝑛₁−𝜛^𝑛−1,𝑅₁^𝑛−1−𝑅Δ𝑡_𝑛1−𝑅₁^𝑛−1Δ𝑡,𝜛^𝑛−1−𝜂^𝑛−𝜂^𝑛−1/Δ𝑡,𝜛^𝑛−𝜂^𝑛−1−𝜂^𝑛−2/Δ𝑡,𝜛^𝑛−1+𝜂Δ𝑡_𝑛−2𝜂^𝑛−1+𝜂^𝑛−2(Δ𝑡)²,𝜛^𝑛−1+12𝑏^𝑛−𝑏^𝑛−1Δ𝑡∇𝜛^𝑛−1,∇𝜛^𝑛−1. (51) Multiply (51) by 2Δ𝑡, sum from 𝑛=1 to 𝐽, and use (19) to get 𝑏₀∇𝜛^𝐽2+2𝑎₀Δ𝑡^𝐽_𝑛=1∇𝜍^𝑛−∇𝜍^𝑛−1Δ𝑡²≤2Δ𝑡^𝐽_𝑛=1𝑎^𝑛−𝑎^𝑛−1𝜽Δ𝑡^𝑛−1∇𝜍^𝑛−∇𝜍^𝑛−1+𝜓Δ𝑡_𝑛−𝜓^𝑛−1𝜍Δ𝑡^𝑛−𝜍^𝑛−1𝜛Δ𝑡+2^𝐽𝑅^𝐽₁+2Δ𝑡^𝐽_𝑛=1𝑅^𝑛₁−𝑅₁^𝑛−1𝜛Δ𝑡^𝑛−1𝜂+2^𝐽−𝜂^𝐽−1𝑐𝜛Δ𝑡^𝐽+2Δ𝑡^𝐽_𝑛=1𝜂^𝑛−2𝜂^𝑛−1+𝜂^𝑛−2(Δ𝑡)²𝜛^𝑛−1+𝐶Δ𝑡^𝐽_𝑛=1∇𝜛^𝑛−12≤𝐶Δ𝑡^𝐽_𝑛=1𝜽^𝑛−1+𝜓^𝑛−𝜓^𝑛−1Δ𝑡²+𝜂^𝑛−2𝜂^𝑛−1+𝜂^𝑛−2(Δ𝑡)²²+𝑅^𝑛₁−𝑅₁^𝑛−1Δ𝑡²+𝑎₀Δ𝑡^𝐽_𝑛=1∇𝜍^𝑛−∇𝜍^𝑛−1Δ𝑡²+𝐶𝑅^𝐽₁²+𝜂^𝐽−𝜂^𝐽−1Δ𝑡²+𝑏₀2∇𝜛^𝐽2+𝐶Δ𝑡^𝐽_𝑛=1∇𝜛^𝑛−12. (52) Use Gronwall lemma to obtain ∇𝜛^𝐽2≤𝐶Δ𝑡^𝐽_𝑛=1𝜽^𝑛−1+𝜓^𝑛−𝜓^𝑛−1Δ𝑡²+𝑅^𝑛₁−𝑅₁^𝑛−1Δ𝑡²+𝜂^𝑛−2𝜂^𝑛−1+𝜂^𝑛−2(Δ𝑡)²²+𝐶𝑅^𝐽₁²+𝜂^𝐽−𝜂^𝐽−1Δ𝑡². (53) Note that 𝑅^𝑛₁−𝑅₁^𝑛−1Δ𝑡²=𝑢^𝑛−2𝑢^𝑛−1+𝑢^𝑛−2(Δ𝑡)²−𝑢_𝑡(𝑡_𝑛)−𝑢_𝑡(𝑡_𝑛−1)Δ𝑡²=𝑢^𝑛−2𝑢^𝑛−1+𝑢^𝑛−2(Δ𝑡)²−𝑢_𝑡𝑡(𝑡_Δ∗)²≤𝐶Δ𝑡^𝑡𝑛_𝑡𝑛−2𝜕⁴𝑢𝜕𝑡⁴𝜕𝑑𝑠+Δ𝑡³𝑢𝜕𝑡³_{𝐿∞(𝐿2)}²≤𝐶(Δ𝑡)^2𝑡𝑛_𝑡𝑛−2𝜕⁴𝑢𝜕𝑡⁴²𝑑𝑠+(Δ𝑡)²𝜕³𝑢𝜕𝑡³²_{𝐿∞(𝐿2)},where𝑡_Δ∗∈𝑡_𝑛−1,𝑡_𝑛. (54) Substitute (44), (47) and (54) into (53) to obtain ∇𝜛^𝐽2≤𝐶𝜓²_𝐿∞𝐿2+^𝑡𝐽_𝑡0𝜂_𝑡(𝑠)²+𝜓_𝑡(𝑠)²+()𝑑𝑠Δ𝑡²𝜕³𝑢𝜕𝑡³²_𝐿∞𝐿2+𝜕²𝑢𝜕𝑡²²_𝐿2𝐿2+𝜕⁴𝑢𝜕𝑡⁴²_𝐿2𝐿2. (55) Use (55) and Poincaré inequality (19) to get 𝜛^𝐽2≤𝐶𝜓²_𝐿∞𝐿2+^𝑡𝐽_𝑡0𝜂_𝑡(𝑠)²+𝜓_𝑡(𝑠)²+()𝑑𝑠Δ𝑡²𝜕³𝑢𝜕𝑡³²_𝐿∞𝐿2+𝜕²𝑢𝜕𝑡²²_𝐿2𝐿2+𝜕⁴𝑢𝜕𝑡⁴²_𝐿2𝐿2. (56) 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: 𝑢_𝑡(((𝑎(𝑡)(),()𝑢()()()+∇⋅𝑏∇∇⋅∇𝑢)))=𝑓𝐱,𝑡𝐱,𝑡∈Ω×𝐽,𝐱,𝑡=Δ𝑢𝐱,𝑡=0,𝐱,𝑡∈𝜕Ω×𝑢()𝐽,𝐱,0=sin𝜋𝑥₁sin𝜋𝑥₂,𝐱=𝑥₁,𝑥₂∈Ω, (57) where Ω=[0,1]×[0,1], 𝐽=(0,1], 𝑎(𝑡)=1+𝑡²,  𝑏=1, and 𝑓(𝐱,𝑡) is chosen so that the exact solution for the scalar unknown function is 𝑢(𝐱,𝑡)=𝑒^−2𝑡sin𝜋𝑥₁sin𝜋𝑥₂, (58) its exact gradient is 𝝀=∇𝑢=𝑒^−2𝑡𝜋cos𝜋𝑥₁sin𝜋𝑥₂sin𝜋𝑥₁cos𝜋𝑥₂, (59) its exact flux function is 𝝈=−𝑎𝝀=𝑒^−2𝑡1+𝑡²𝜋cos𝜋𝑥₁sin𝜋𝑥₂sin𝜋𝑥₁cos𝜋𝑥₂, (60) and the diffusion term is 𝛾=∇⋅𝝈=2𝜋²𝑒^−2𝑡1+𝑡²sin𝜋𝑥₁sin𝜋𝑥₂. (61)

Dividing the domain Ω into the triangulations of mesh size ℎ uniformly and choosing the piecewise linear space 𝑉_ℎ with index 𝑘=1 and ℎ=√2Δ𝑡=√2/8,√2/16,√2/32,√2/64, we obtain the optimal a priori error estimates in 𝐿² and 𝐻¹-norm for the scalar unknown 𝑢 in Table 1. Similarly, from Table 2, we can find that both ‖𝛾−𝛾_ℎ‖_{𝐿∞(𝐻1(Ω))} and ‖𝛾−𝛾_ℎ‖_{𝐿∞(𝐿2(Ω))} for the diffusion term 𝛾 are optimal, too. At the same time, we take the corresponding piecewise constant space 𝐖_ℎ with index 𝑘=0 for the gradient 𝝀 and the flux 𝝈 and obtain some convergence results in (𝐿²)²-norm with ℎ=√2Δ𝑡=√2/8,√2/16,√2/32,√2/64 in Table 3.

Table 1 

𝐿² and 𝐻¹-norms and convergence order for 𝑢.

Table 2 

𝐿² and 𝐻¹-norms and convergence order for 𝛾.

Table 3 

(𝐿²)²-norm and convergence order for 𝜆 and 𝜎.

Taking 𝑡=1 with ℎ=√2Δ𝑡=√2/32, we show the surfaces of the exact solution 𝑢, 𝝈, 𝝀, and 𝛾 in Figures 1, 3, 5, and 7, respectively. Figures 2, 4, 6, and 8 show the surfaces of the numerical solution 𝑢_ℎ, 𝝈_ℎ, 𝝀_ℎ and 𝛾_ℎ at 𝑡=1 with ℎ=√2Δ𝑡=√2/16, respectively.

Figure 1 

Surface for exact solution 𝑢.

Figure 2 

Surface for numerical solution 𝑢ℎ.

Figure 3 

Surface for exact solution 𝝈=(𝝈1,𝝈2).

Figure 4 

Surface for numerical solution 𝝈ℎ=(𝝈1ℎ,𝝈2ℎ).

Figure 5 

Surface for exact solution 𝝀=(𝝀1,𝝀2).

Figure 6 

Surface for numerical solution 𝝀ℎ=(𝝀1ℎ,𝝀2ℎ).

Figure 7 

Surface for exact solution 𝛾.

Figure 8 

Surface for numerical solution 𝛾ℎ.

From the convergence results in Tables 1–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 (𝐿²(Ω))² space; therefore, the regularity requirement on the gradient solution 𝝀=∇𝑢 is reduced. We obtain the optimal priori error estimates in 𝐿² and 𝐻¹-norm for both the scalar unknown 𝑢 and the diffusion term 𝛾 and the priori error estimates in (𝐿²)²-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.

Acknowledgments

The authors thank the referees and editor for their valuable suggestions and comments, which greatly improve the paper. This work is supported by National Natural Science Fund (11061021), Natural Science Fund of Inner Mongolia Autonomous Region (2012MS0108, 2012MS0106), Scientific Research Projection of Higher Schools of Inner Mongolia (NJZZ12011, NJ10006, NJZY13199), and Program of Higher-level talents of Inner Mongolia University (125119, 30105-125132).