- Email: [email protected]
A stabilized mixed finite element method for steady and unsteady reaction–diffusion equations
Available online at www.sciencedirect.com
ScienceDirect Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117 www.elsevier.com/locate/cma
A stabilized mixed finite element method for steady and unsteady reaction–diffusion equations Hongfei Fu a,∗ , Hui Guo a , Jian Hou b , Junlong Zhao a a College of Science, China University of Petroleum, Qingdao 266580, China b College of Petroleum Engineering, China University of Petroleum, Qingdao 266580, China
Received 18 June 2015; received in revised form 1 December 2015; accepted 15 January 2016 Available online 12 February 2016
Abstract In this paper, we propose a new mixed finite element method, called stabilized mixed finite element method, for the approximation of steady reaction–diffusion partial differential equations (PDEs). The method is obtained by translating the primal second-order PDEs into a first-order mixed system, and then adding some suitable elementwise residual terms multiplied by a stabilization parameter to the weak formulation. The new method is compatible, i.e., the added terms equal to zero in the continuous case. Furthermore, it is mesh-independent, i.e., the stabilization parameter is independent of the mesh size. We prove both coercive and continuous properties in a weighted norm for the corresponding new mixed bilinear formulation. These assure that the finite element function spaces do not require to satisfy the classical Ladyzhenkaya–Babuska–Brezzi (LBB) consistency condition. Therefore, the widely used Lagrange finite element can be adopted. A simple proof of a priori error estimate with lower order regularity requirement is discussed, and numerical experiments confirm the efficiency and reliability of the new stabilized mixed method. Finally, the method is applied to solving unsteady reaction–diffusion equations. Error estimates are also given, and numerical examples still support the theoretical analysis very well. c 2016 Elsevier B.V. All rights reserved. ⃝
MSC: 65N15; 65N30 Keywords: Reaction–diffusion equation; Stabilized mixed finite element; Classical mixed finite element; A priori error estimate; Numerical experiments
1. Introduction Let Ω be a bounded convex domain in Rd (d = 1, 2, 3) with Lipschitz boundary ∂Ω . In this paper, we are interested in proposing a stabilized mixed finite element method for the following linear reaction–diffusion equations: −div(A∇ y) + cy = f, in Ω , (1.1) y = 0, on ∂Ω . ∗ Corresponding author.
E-mail addresses: [email protected] (H. Fu), [email protected] (H. Guo), [email protected] (J. Hou). http://dx.doi.org/10.1016/j.cma.2016.01.010 c 2016 Elsevier B.V. All rights reserved. 0045-7825/⃝
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
103
Depending on a simplified model for particular real problems, the variable y can represent a temperature, a pressure, or a concentration of a chemical species. The matrix A = A(x) describes the diffusion feature, and c = c(x) stands for a reaction rate. Detailed description is in the following section. As in the real problems, the gradient of the primal variable is also meaningful. For example, in the flow problem, it stands for a velocity; or, in the temperature control problem, large temperature gradients during cooling or heating may lead to its destruction. Therefore, in these cases people pay their special attention on the flux of the primal variable. Let us introduce the flux σ = −A∇ y, then (1.1) is equivalent to the following mixed form in Ω , divσ + cy = f, (1.2) A−1 σ + ∇ y = 0, in Ω , y = 0, on ∂Ω . It is well known that (1.2) is widely studied by many authors using the classical mixed finite element method; see, for example, [1–5]. The most important merit of the mixed finite element method is that the flux σ can be approximated simultaneously to a same order of accuracy as the primal variable y, without any accuracy lost compared to the standard finite element method. However, the finite element function spaces for y and σ have to be specially chosen such that they strictly satisfy the LBB consistency condition. In this sense, some of the best known and widely used finite element spaces are thereby excluded. In particular, in high-dimensional problems, it is not a straightforward thing to construct such kinds of mixed finite element function spaces. Meanwhile, the practical computation for classical mixed finite element methods is complicated, people need to solve a saddle-point problem which is symmetric but indefinite. To overcome these difficulties, some stabilization techniques have been studied, for example, least-squares methods [6,7], H 1 -Galerkin methods [8], as well as splitting methods [9,10]. In this paper, different from the previous approaches, we shall construct a stabilized mixed finite element method for the problem under consideration. The method can be easily obtained by adding some suitable elementwise residual terms to a dual-mixed weak formulation. We prove both coercive and continuous properties in a weighted norm for the corresponding new mixed bilinear formulation. Thus, the new method can avoid the drawbacks in selecting finite element function spaces. It is also competitive in practical computation, as the coefficient matrix of the leading system of linear algebraic equation is positive definite, and there are fewer degree of freedoms (Dofs) compared to the corresponding classical mixed finite element method. The rest of the paper is organized as follows: In Section 2, starting from the classical mixed formulation, we first introduce the stabilized mixed weak formulation, and then prove the existence and uniqueness of its solution. Finally, discretization using standard Lagrange elements is discussed. In Section 3, a simple error analysis for the stabilized mixed finite element method is given. In Section 4, we conduct some numerical experiments to verify the theoretical analysis. We also compare the results of the new proposed stabilized mixed method with those obtained by the classical mixed method. In Section 5, we extent the stabilized method to solving unsteady reaction–diffusion problems. Error estimate and numerical test are also presented. Finally, some concluding remarks are given in the last section. 2. Stabilized mixed finite element method Let ω ⊂ Ω . Throughout this paper, we adopt the standard notations L p (ω) (1 ≤ p ≤ ∞) for Lebesgue space of real-valued functions with norm ∥ · ∥0, p,ω , and W m, p (ω) (1 ≤ p ≤ ∞) for Sobolev spaces endowed with norm ∥ · ∥m, p,ω and semi-norm | · |m, p,ω . For p = 2, we denote W m,2 (ω) = H m (ω) and we drop the subscript p = 2 in the corresponding norms and semi-norms. Besides, the subscript ω will also be omitted in the norms if ω = Ω . To give a detailed description of the reaction–diffusion problem under consideration in a mixed weak formulation, we shall take the function spaces V = L 2 (Ω ),
W = H (div; Ω ) , {τ ∈ L 2 (Ω )d : divτ ∈ L 2 (Ω )},
Y = H01 (Ω ),
Σ = L 2 (Ω )d .
and
Besides, in this paper we need the following assumptions. (H-1) The solution (y, σ ) of (1.2) satisfies the following regularity: y ∈ H k+1 (Ω ),
σ ∈ H r (Ω )d .
104
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
(H-2) A = (ai, j (x))d×d is symmetric, and there exist positive constants α and β such that α|X |2 ≤ X T AX ≤ β|X |2 ,
∀X ∈ Rd .
(H-3) c = c(x) is bounded below and up by two constants c∗ and c∗ such that 0 ≤ c∗ ≤ c ≤ c∗ . Now, let us first recall the classical mixed weak formulation of (1.2). It aims to find (y, σ ) ∈ V × W such that −1 (A σ, τ ) − (y, divτ ) = 0, ∀τ ∈ W, (2.1) (divσ, v) + (cy, v) = ( f, v), ∀v ∈ V or equivalently (A−1 σ, τ ) − (y, divτ ) + (divσ, v) + (cy, v) = ( f, v),
∀(v, τ ) ∈ V × W.
(2.2)
This idea is firstly proposed by Raviart and Thomas to second-order elliptic problems, see [1]. From then on, it is widely discussed and used to solve various equations. It is also extended to combine with other numerical methods, such as discontinuous Galerkin method, finite volume method, and so on. We cannot even briefly summarize it here. An alternative approach uses the concept of least-squares residual stabilization as proposed by Hughes et al. [11] for solving advection–diffusion equations, and theoretically analyzed by Franca and Stenberg [12] for the elasticity equations. They call the new method as Galerkin least-squares method. After then, an unusual stabilized finite element method, see [13–15] is discussed and analyzed. The difference between these two stabilized methods appears in the sign of second-order term in the test function, see [13] for details. Inspired by the previous work, we hope to adopt the structure of these stabilized finite element methods, and propose a stabilized mixed finite element method to conquer the stability problem in classical mixed finite element method. To do this, let us first rewrite Eq. (2.2) as follows: (A−1 σ, τ ) + (∇ y, τ ) − (σ, ∇v) + (cy, v) = ( f, v),
∀(v, τ ) ∈ Y × Σ .
(2.3) (A−1 σ
Then we augment the dual-mixed formulation (2.3) by a non-symmetric residual term + ∇ y, τ − A∇v)e on each element e to propose a stabilized mixed method. In this sense, the new mixed weak formulation of (1.2) reads as finding (y, σ ) ∈ Y × Σ such that aδ (y, σ ; v, τ ) = ( f, v),
∀(v, τ ) ∈ Y × Σ .
(2.4)
Here the new bilinear form aδ (y, σ ; v, τ ) , (A−1 σ, τ ) + (∇ y, τ ) − (σ, ∇v) + (cy, v) −
e∈T
δ(A−1 σ + ∇ y, τ − A∇v)e ,
(2.5)
h
where δ > 0 is an elementwise constant parameter, (·, ·)e denotes the inner product in L 2 (e)d , and the element e belongs to a regular triangulation T h of the given domain Ω , which is defined below. To prove the existence and uniqueness of the solution pair to the new stabilized mixed weak formulation (2.4), we first define a weighted δ- stabilization norm corresponding to the bilinear form aδ (, ; , ) as follows: |||{y, σ }|||2stab = (A−1 σ, σ ) + δ(A∇ y, ∇ y)e + (cy, y). (2.6) e∈T h
Then, we can derive the following coercive and bounded results which are proved similarly in Ref. [16], but with only one stabilization parameter here. For readers’ convenience, we still include the proof here. Proposition 2.1. For 0 < δ0 ≤ δ ≤ δ1 < 1, we have aδ (y, σ ; y, σ ) ≥ cδ |||{y, σ }|||2stab , Moreover, for all (y, σ ; v, τ ) ∈ (Y
∀(y, σ ) ∈ Y × Σ .
(2.7)
× Σ )2
aδ (y, σ ; v, τ ) ≤ Cδ |||{y, σ }|||stab |||{v, τ }|||stab . Here δi (i = 0, 1) are two positive constants, and cδ = 1 − δ1 and Cδ = max{4, 2 + 1/δ0 }.
(2.8)
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
Proof. On the one hand, it follows from the definition that aδ (y, σ ; y, σ ) = (A−1 σ, σ ) + (cy, y) − δ(A−1 σ, σ )e + δ(A∇ y, ∇ y)e . e∈T h
105
(2.9)
e∈T h
Hence, for δ ≤ δ1 < 1, the bilinear form a(, ; , ) is coercive on Y × Σ , and satisfy (2.7). On the other hand, it follows from the definition that aδ (y, σ ; v, τ ) = (A−1 σ, τ ) + (∇ y, τ ) − (σ, ∇v) + (cy, v) − δ(A−1 σ, τ )e e∈T h
+
δ(σ, ∇v)e −
e∈T h
δ(∇ y, τ )e +
e∈T h
δ(A∇ y, ∇v)e .
(2.10)
e∈T h
Note that i
|ai bi | ≤
ai2
1 2
i
1
2
.
bi2
i
Then for 0 < δ0 ≤ δ ≤ δ1 < 1 we have 1 2 aδ (y, σ ; v, τ ) ≤ 4(A−1 σ, σ ) + (A∇ y, ∇ y) + 2 δ(A∇ y, ∇ y)e + (cy, y) e∈T h
1 2 4(A−1 τ, τ ) + (A∇v, ∇v) + 2 δ(A∇v, ∇v)e + (cv, v) e∈T h
−1
≤ 4(A
1 σ, σ ) + 2 + δ0
δ(A∇ y, ∇ y)e + (cy, y)
21
e∈T h
21 1 4(A−1 τ, τ ) + 2 + δ(A∇v, ∇v)e + (cv, v) δ0 h e∈T
≤ Cδ |||{y, σ }|||δ |||{v, τ }|||δ , which proves the bounded result (2.8).
(2.11)
Proposition 2.2. Let f ∈ L 2 (Ω ). Under assumptions (H-2)–(H-3) , Problem (2.4) admits a unique solution pair (y, σ ) ∈ Y × Σ . Moreover, the solution exhibits the improved regularity (y, σ ) ∈ H 2 (Ω ) × H 1 (Ω )d , and satisfies the following stability estimate ∥y∥2 + ∥σ ∥1 ≤ C∥ f ∥.
(2.12)
Proof. Proposition 2.1 implies that the bilinear form a(, ; , ) is coercive and bounded. Thus the Lax–Milgram lemma implies the existence and uniqueness of the solution pair (y, σ ) ∈ Y × Σ to problem (2.4). Furthermore, stability result (2.12) follows from the elliptic regularity depending on the shape of the domain. Now we consider the finite element discretization. Let T h be a family of regular partition of Ω , such that Ω = ∪e∈T h e. Let h e denote the diameter of the element e ∈ T h , and h = maxe∈T h h e . As the bilinear form a(, ; , ) is coercive and bounded, therefore we can choose the standard Lagrange element as the flux function spaces. That is to say, the LBB consistency condition which is required by classical mixed element spaces, for example, RT elements [1], now is not necessary. For two integers k ≥ 1, r ≥ 0, we construct the finite element subspaces Y h ⊂ Y and Σ h ⊂ Σ in a standard way: Y h = {vh ∈ C 0 (Ω ) : vh |e ∈ Pk (e), ∀e ∈ T h , vh = 0 on ∂Ω }, Σ h = {τh ∈ L 2 (Ω )d : (τh )i |e ∈ Pr (e), i = 1, 2, . . . , d, ∀e ∈ T h }. Here Pm denotes polynomials of total degree at most m. We can also adopt r th order RT elements as Σ h .
(2.13)
106
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
Then the discretization of the state equation (2.4) reads as follows: Find (yh , σh ) ∈ Y h × Σ h such that aδ (yh , σh ; vh , τh ) = ( f, vh ),
∀(vh , τh ) ∈ Y h × Σ h .
(2.14)
Remark 2.3. Note that the new stabilized mixed finite element method is compatible. Besides, it satisfies the following orthogonality relation: aδ (y − yh , σ − σh ; vh , τh ) = 0,
∀(vh , τh ) ∈ Y h × Σ h .
(2.15)
Remark 2.4. It is well known that the discretization of (2.1) using the classical mixed finite element spaces, see, for example, RT element [1], BDDM element [4], BDFM element [5], BDM element [3], will lead to a typical saddle-point problem. Therefore, the popular conjugate gradient or algebraic multigrid solvers cannot be adopted. However, the new stabilized mixed finite element method used in (2.14) not only has a free choice of finite element function spaces (see (2.13)), but also produces a positive definite linear algebraic equation (see Proposition 2.1). Therefore, it can be solved rapidly and efficiently with less degree of freedoms (Dofs) compared to the classical mixed method. We will illustrate an iteration algorithm for decoupling the finite element solutions yh and σh of (2.14) in Section 4. 3. A priori error estimates In this section, we aim to derive optimal a priori error estimate for the stabilized mixed finite element approximation of steady reaction–diffusion equations. In the following, let us first recall two important results without proof. Lemma 3.1. Let Ih be the standard Lagrange interpolation operator defined in Ref. [17]. Then there is a constant C > 0 such that ∥v − Ih v∥0,e + h e |v − Ih v|1,e ≤ Ch m e |v|m,e , for any v ∈
H m (e)
and e ∈ T
(3.1)
h.
Theorem 3.2. Let (y, σ ) and (yh , σh ) be the solutions of (2.4) and (2.14), respectively. Assume that the solutions y ∈ H01 (Ω ) ∩ H k+1 (Ω ) and σ ∈ H r (Ω )d . Then the following estimates hold (3.2) |||{y − yh , σ − σh }|||stab ≤ Ch min(k,r ) ∥y∥ H k+1 (Ω ) + ∥σ ∥ H r (Ω )d . Proof. Following from the orthogonality property (2.15), the theorem can be easily proved by the interpolation estimation theory in Lemma 3.1. Let Ih y ∈ Y h and Ih σ ∈ Σ h be defined as the kth and r th order Lagrange interpolations of y and σ , respectively. Then we conclude from (2.15) that aδ (Ih y − yh , Ih σ − σh ; vh , τh ) = −aδ (y − Ih y, σ − Ih σ ; vh , τh ).
(3.3)
Let vh = Ih y − yh and τh = Ih σ − σh in Eq. (3.3). It then follows from Proposition 2.1 and Cauchy–Schwarz inequality that cδ |||{Ih y − yh , Ih σ − σh }|||2δ ≤ Cδ |||{y − Ih y, σ − Ih σ }|||δ |||{Ih y − yh , Ih σ − σh }|||δ .
(3.4)
This implies that |||{Ih y − yh , Ih σ − σh }|||stab ≤ C|||{y − Ih y, σ − Ih σ }|||stab .
(3.5)
Thus we deduce by the triangle inequality and Lemma 3.1 that |||{y − yh , σ − σh }|||stab ≤ |||{y − Ih y, σ − Ih σ }|||stab + |||{Ih y − yh , Ih σ − σh }|||stab ≤ C|||{Ih y − y, Ih σ − σ }|||stab ≤ Ch min(k,r ) ∥y∥ H k+1 (Ω ) + ∥σ ∥ H r (Ω )d , which proves Theorem 3.2.
(3.6)
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
107
4. Numerical experiments In this section, we are going to carry out some numerical experiments to numerically support the a priori error estimates derived in the last section. To solve the problem numerically, we use a C++ software package: AFEPack, see Ref. [18] for the details. Let Ω = [0, 1] × [0, 1] and the stabilized parameter δ = 0.75 in the whole domain Ω . We consider piecewise linear C 0 elements for both the approximations of the primal variable y and its flux σ for the new proposed stabilized mixed finite element method (2.14). To show more convincingly the benefits of the stabilized mixed method, we also compare the results with those obtained using the classical mixed finite element approximation of (2.1). As the classical mixed method is only stable for special choice of matched mixed finite element spaces, we consider here the lowest-order RT elements (see [19,1]) for the approximation of the flux σ , while matched piecewise constant elements are employed for the approximation of the primal variable y. In this part, we consider the following iteration algorithm for solving the coupled mixed system. Algorithm. Step 1. Give an initial value σh0 and tolerance TOL = 1.0E−05. Step 2. Solve yhm ∈ V h such that δ(A∇ yhm , ∇vh ) + (cyhm , vh ) = ( f, vh ) + (1 − δ)(A−1 σhm−1 , ∇vh ). Step 3. Solve σhm ∈ W h such that (A−1 σhm , τh ) = −(∇ yhm , τh ). Step 4. Stop if ∥σhm − σhm−1 ∥ L 2 (Ω )2 < TOL is satisfied, and let yh = yhm , σh = σhm . Otherwise set m := m + 1 and go to Step 2 again.
Example 4.1. For the first example, we consider a diffusion-dominated reaction–diffusion problem with coefficients A = I and c = 1.0E−03. Let the solutions be as follows: y(x) = sin(πx1 ) sin(π x2 ), cos(π x1 ) sin(π x2 ) σ (x) = −π . sin(π x1 ) cos(π x2 ) The source function f can be computed using the above functions. Table 1 displays the weighted errors with respect to the decreasing mesh size. We can observe that the convergence orders are consistent with Theorem 3.2. We also show the main CPU time for the iteration algorithm in Table 1. It can be found that 52.26 s of CPU time is employed when h = 1/128 for the stabilized mixed finite element method using an AMG (algebraic multigrid method) solver. However, to obtain a same stop criterion TOL = 1.0E−05, we can observe from Table 2 that the classical mixed finite element method using a preconditioned GMRES (generalized minimal residual method) solver costs 279.08 s of CPU time, and at the meantime the error |||{y − yh , σ − σh }||| RT := ∥y − yh ∥+∥σ −σh ∥ H (div;Ω ) = 7.18E−02. We should notice that here a preconditioner is adopted for fast computation of the nonsymmetric algebraic equation systems result from the classical mixed method. Otherwise, the classical mixed method shall spend more CPU time for a fine mesh size. Figs. 1–2 show the profiles of the numerical primal variable and its flux with h = 1/128 using the stabilized mixed method. It can be seen that the numerical figures are also agreement with the analytical solutions very well. Example 4.2. For the second example, we consider a reaction–diffusion problem with coefficients A = I and c = 1, where the reaction and diffusion are of the same importance. Let the solutions be as follows: y(x) = x1 (1 − x1 )x2 (1 − x2 ) exp(x1 − x2 ), (2x1 − 1)x2 (1 − x2 ) −1 σ (x) = exp(x1 − x2 ) +y . (2x2 − 1)x1 (1 − x1 ) 1 The source function f can be derived using the above functions.
108
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
Fig. 1. The numerical primal variable with h = 1/128.
Table 1 Results for the stabilized mixed finite element method with δ = 0.75. Mesh
||{y − yh , σ − σh }|| δ
Order
CPU time (s)
16*16 32*32 64*64 128*128
1.34E−01 6.70E−02 3.34E−02 1.67E−02
– 1.00 1.00 1.00
0.73 3.14 13.06 52.26
Table 2 Results for the classical mixed finite element method. Mesh
||{y − yh , σ − σh }|| RT
Order
CPU time (s)
16*16 32*32 64*64 128*128
5.73E−01 2.87E−01 1.43E−01 7.18E−02
– 1.00 1.01 0.99
0.19 1.51 14.41 279.08
Table 3 Results for the stabilized mixed finite element method with δ = 0.75. Mesh
||{y − yh , σ − σh }|| δ
Order
CPU time (s)
16*16 32*32 64*64 128*128
1.34E−02 6.66E−03 3.29E−03 1.64E−03
— 1.01 1.02 1.00
0.6 2.45 9.79 40.12
We present the corresponding errors and CPU times for the stabilized/classical mixed finite element methods with respect to different mesh sizes in Tables 3 and 4. Numerical figures of the primal and flux variables using the stabilized mixed method are shown in Figs. 3 and 4. Upon the same stop criterion TOL = 1.0E−05, the classical mixed finite element method spends 80.49 s of CPU time, which is about twice of the stabilized mixed finite element method for the mesh size h = 1/128. It is clear that the new stabilized mixed method is much more efficient for large computation.
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
109
Fig. 2. The numerical flux variable with h = 1/128.
0
Fig. 3. The numerical primal variable with h = 1/128.
Table 4 Results for the classical mixed finite element method. Mesh
||{y − yh , σ − σh }|| RT
Order
CPU time (s)
16*16 32*32 64*64 128*128
3.64E−02 1.83E−02 9.29E−03 5.72E−03
– 0.99 0.98 0.70
0.18 1.31 9.95 80.49
Example 4.3. For the third example, we consider a reaction-dominated reaction–diffusion problem with coefficients A = ε I and c = 1. Let ε = 4.0E−04 and the solution pair as follows: (x1 − 0.5)2 + (x2 − 0.5)2 y(x) = exp − , √ ε √ x − 0.5 σ (x) = 2 εy 1 . x2 − 0.5
110
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
0
0
Fig. 4. The numerical flux variable with h = 1/128.
Table 5 Results for the stabilized mixed finite element method with δ = 0.75. Mesh
||{y − yh , σ − σh }|| δ
Order
CPU time (s)
16*16 32*32 64*64 128*128
7.60E−03 3.49E−03 1.71E−03 8.49E−04
– 1.12 1.03 1.01
0.27 1.2 4.72 18.82
Table 6 Results for the classical mixed finite element method. Mesh
||{y − yh , σ − σh }|| RT
Order
CPU time (s)
16*16 32*32 64*64 128*128
2.25E−02 1.13E−02 5.67E−03 2.84E−03
– 0.99 0.99 1.00
0.11 0.28 1.35 7.65
The source function f can be computed using the above functions. We present weighted errors and CPU times with respect to different mesh sizes in Table 5. Numerical figures of the primal and flux variables are shown in Figs. 5 and 6. This example provides a very sharp solution around the point (0.5, 0.5). However, the test results using the stabilized mixed finite element method also support the theoretical analysis in Theorem 3.2 very well. As this example presents a nearly zero solution in most area of the domain, the preconditioned GMRES solver embodies a high efficiency because of the sparsity of the leading algebraic matrix, and from Table 6 we can see that only 7.65 s is employed by the classical mixed finite element method with error 2.84E−03 when h = 1/128. However, we could still anticipate that when the mesh size is decreasing, the stabilized mixed method shall be more competitive compared to the classical one. As when h is getting smaller, there will be a more increasing degree of freedoms (Dofs) for the classical mixed method than that for the stabilized mixed method. In summary, numerical experiments in this section show the high efficiency and reliability of the stabilized mixed finite element method, which coincides with the theoretical analysis. In most cases, few CPU time and small errors are observed compared to the classical mixed finite element method even using an efficient preconditioned GMRES solver.
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
111
Fig. 5. The numerical primal variable with h = 1/128.
Fig. 6. The numerical flux variable with h = 1/128.
5. Application to unsteady reaction–diffusion equations 5.1. Preliminary and discretization In this section, we consider an extension of stabilized mixed finite element method for the following unsteady reaction–diffusion equations: yt − div(A∇ y) + cy = f, in Ω × (0, T ], y = 0, on ∂Ω × (0, T ], (5.1) y(x, 0) = y0 (x), in Ω where T < +∞ denotes the final time, and y0 (x) is a prescribed initial data. To define a fully-discrete scheme, we introduce a time partition of [0, T ] by t i = i∆t, i = 0, 1, 2, . . . , N with ∆t = T /N where N is a positive integer. Let f i = f (·, t i ). We denote by L s (0, T ; W m, p (Ω )) the Banach space of
112
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
all L s integrable functions from (0, T ] into W m, p (Ω ) with norm ∥v∥ L s (0,T ;W m, p (Ω )) =
1s
T
for s ∈ [1, ∞)
∥v∥m, p dt 0
and the standard modification for s = ∞. Similarly, let σ = −A∇ y, then the first equation in (5.1) is equivalent to the following mixed form yt + divσ + cy = f, in Ω × (0, T ], A−1 σ + ∇ y = 0, in Ω × (0, T ]. Therefore, a stabilized mixed weak formulation of (5.2) reads as finding (y, σ ) : (0, T ] → Y × Σ such that (yt , v) + aδ (y, σ ; v, τ ) = ( f, v), ∀(v, τ ) ∈ Y × Σ , y(0) = y0 .
(5.2)
(5.3)
Let Y h and Σ h be the finite element spaces as defined in (2.13). Then the discretization of Eq. (5.3) reads as follows: Find (yhi , σhi ) ∈ Y h × Σ h for i = 1, 2, . . . , N such that i i−1 yh − yh , v + a (y i , σ i ; v , τ ) = ( f i , v ), ∀(v , τ ) ∈ Y h × Σ h , h δ h h h h h h h ∆t (5.4) 0 yh = y0,h , where y0,h ∈ Y h is an approximation of y0 which will be specified later on. 5.2. A priori error estimate In this section, based on the error analysis for the steady problem in Section 3, we use an elliptic projection to give a priori error estimate for the unsteady reaction–diffusion problem. To this aim, we define ( yh , σh ) : (0, T ] → Y h × Σ h such that aδ (y − yh , σ − σh ; vh , τh ) = 0,
∀(vh , τh ) ∈ Y h × Σ h .
(5.5)
Then it follows from (5.5) and Theorem 3.2, we have the following auxiliary result. Lemma 5.1. Let (y, σ ) and ( yh , σh ) be the solutions of (5.3) and (5.5), respectively. For any t ∈ (0, T ], assume that the solutions y(t) ∈ H01 (Ω ) ∩ H k+1 (Ω ) and σ (t) ∈ H r (Ω )d . Then the following estimates hold |||{y − yh , σ − σh }|||stab ≤ Ch min(k,r ) ∥y∥ H k+1 (Ω ) + ∥σ ∥ H r (Ω )d . (5.6) In the following, we just need to estimate the bounds for yh − yh and σh − σh . It is summarized as the following lemma. Lemma 5.2. Let (yh , σh ) and ( yh , σh ) be the solutions of (5.4) and (5.5), respectively. Assume that yh (0) = yh (0) and the solutions y ∈ H 1 (0, T ; H k+1 (Ω )) ∩ H 2 (0, T ; L 2 (Ω )) and σ ∈ H 1 (0, T ; H r (Ω )d ). Then we have ∥ yh − yh ∥ L ∞ (0,T ;L 2 (Ω )) + |||{ yh − yh , σh − σh }||| L 2 (0,T ;stab) ≤ Ch min(k,r ) ∥yt ∥ L 2 (0,T ;H k+1 (Ω )) + ∥σt ∥ L 2 (0,T ;H r (Ω )d ) + C∆t∥ytt ∥2L 2 (0,T ;L 2 (Ω )) , where |||{ξ, ζ }|||2L 2 (0,T ;stab) :=
N
2
∆t|||{ξ i , ζ i }|||stab
i=1
is a time-dependent discrete norm.
(5.7)
113
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
Proof. For simplicity, we split the errors at time t = t i as follows: y i − yhi = (y i − yhi ) + ( yhi − yhi ) = ρ i + ξ i , σ i − σhi = (σ i − σhi ) + ( σhi − σhi ) = ηi + ζ i . It then follows from (5.3) to (5.4) and the auxiliary projection (5.5) that for any (vh , τh ) ∈ Y h × Σ h i i ξ − ξ i−1 ρ − ρ i−1 i i , vh + aδ (ξ , ζ ; vh , τh ) = − , vh − (χ i , vh ), ∆t ∆t
(5.8)
where χ i := yti −
y i − y i−1 . ∆t
Let the test functions vh = ξ i , τh = ζ i in Eq. (5.8). Note that i ξ − ξ i−1 i 1 i 2 ,ξ ≥ ∥ξ ∥ − ∥ξ i−1 ∥2 . ∆t 2∆t Then we conclude from Proposition 2.1 that the left-hand side of (5.8) is bounded below by 1 i 2 2 ∥ξ ∥ − ∥ξ i−1 ∥2 + cδ |||{ξ i , ζ i }|||stab . 2∆t
(5.9)
For the right-hand side of (5.8), it is a matter of calculation that i i t t i i−1 |ρ − ρ | = ρt dt ≤ |ρt | dt, t i−1 t i−1 i 1 ti t i (t − ti−1 )ytt dt ≤ |χ | = |ytt | dt. ∆t t i−1 t i−1 Then, following Cauchy–Schwarz inequality we have
ρ i − ρ i−1 i ,ξ ∆t
1 ≤ ∆t
ti
t i−1
∥ρt ∥dt · ∥ξ i ∥ ≤
1 ∥ρt ∥2L 2 (t i−1 ,t i ;L 2 (Ω )) + C∥ξ i ∥2 , ∆t
(5.10)
(χ i , ξ i ) ≤ ∆t∥ytt ∥2L 2 (t i−1 ,t i ;L 2 (Ω )) + C∥ξ i ∥2 .
(5.11)
The inequalities (5.10)–(5.11) can be combined with (5.9) to give the following recursion relation 1 i 2 2 ∥ξ ∥ − ∥ξ i−1 ∥2 + cδ |||{ξ i , ζ i }|||stab 2∆t 1 ∥ρt ∥2L 2 (t i−1 ,t i ;L 2 (Ω )) + ∆t∥ytt ∥2L 2 (t i−1 ,t i ;L 2 (Ω )) + C∥ξ i ∥2 . ≤ ∆t
(5.12)
Let the initial value y0,h = yh (0), i.e., ξ 0 = 0. Then if (5.12) is multiplied by 2∆t and summed over i from 1 to N , we have ∥ξ N ∥2 + C
N i=1
2
∆t|||{ξ i , ζ i }|||stab ≤ C∥ρt ∥2L 2 (0,T ;L 2 (Ω )) + C(∆t)2 ∥ytt ∥2L 2 (0,T ;L 2 (Ω )) + C
Then Lemma 5.2 follows from Lemma 5.1 and the discrete Gronwall’s lemma.
N
∆t∥ξ i ∥2 .
i=1
Finally, collecting the bounds in Lemmas 5.1, 5.2, and using the triangle inequality we can easily derive the following main theorem.
114
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
Theorem 5.3. Let (y, σ ) and (yh , σh ) be the solutions of (5.3) and (5.4), respectively. Assume that all conditions in Lemmas 5.1 and 5.2 are valid. Furthermore, suppose that y ∈ L ∞ (0, T ; H01 (Ω ) ∩ H k+1 (Ω )) and σ ∈ L ∞ (0, T ; H r (Ω )d ). Then the following estimate holds ∥y − yh ∥ L ∞ (0,T ;L 2 (Ω )) + |||{y − yh , σ − σh }||| L 2 (0,T ;stab) ≤ C(h min(k,r ) + ∆t).
(5.13)
Remark 5.4. To improve the convergence order in the temporal direction, one can adopt the following Crank–Nicolson stabilized mixed finite element discretization of (5.3) as yhi − yhi−1 1 i− 1 i− 1 , vh + aδ (yh 2 , σh 2 ; vh , τh ) = ( f i− 2 , vh ), ∀(vh , τh ) ∈ Y h × Σ h , (5.14) ∆t where yhi + yhi−1 , 2 Then one can expect i− 12
yh
=
i− 12
σh
=
σhi + σhi−1 . 2
∥y − yh ∥ L ∞ (0,T ;L 2 (Ω )) + |||{y − yh , σ − σh }||| L 2 (0,T ;stab) ≤ C(h min(k,r ) + (∆t)2 ).
(5.15)
5.3. Numerical tests In this subsection, we consider piecewise linear elements for both the approximations of y and σ , i.e., k = r = 1. Besides, an iteration algorithm as the steady case for solving the coupled mixed system is also employed, while the stop criterion is replaced by ∥σhN ,m − σhN ,m−1 ∥ L 2 (Ω )2 < TOL. Example 5.1. In this example, we numerically test a reaction-dominated reaction–diffusion problem on Ω = [−1, 1] × [−1, 1] and T = 0.5. Let the small diffusion coefficient A = ε I with ε = 4.0E−03, the reaction rate c = 1, and the sharp solutions be as follows: 1.0 − cos(t − x1 ) 1.0 − cos(t − x2 ) y(x, t) = exp − exp − ; ε ε sin(t − x1 ) σ (x, t) = −y . sin(t − x2 ) The source function f can be computed using the above functions. We list the time-dependent |||{y − yh , σ − σh }||| L 2 (0,T ;stab) errors and main CPU times with respect to different spatial mesh sizes and temporal sizes in Table 7. A first-order spatial and temporal convergence can be observed for the stabilized mixed method, which is in accordance with the theoretical analysis in Theorem 5.3 very well. It is well known that the solutions have very steep moving fronts, and when the time is marching the position of moving fronts is shifted along the straight line (x1 , x2 ) = (t, t). These phenomena can be observed clearly from Figs. 7–9, where the contour plots of numerical primal and flux variables are shown at t = 0, 0.25, and 0.5 respectively for h = 1/64. Furthermore, it is apparent that the maximum numerical value of the primal variable is indeed centered at (t, t) for different t. This is in accordance with the analytical solution very well. Example 5.2. For the second example, we consider a diffusion-dominated reaction–diffusion problem with coefficients A = I and c = 1.0E−03 on Ω = [0, 1] × [0, 1] and T = 1. Let the solutions be as follows: y(x, t) = 16t x1 (1 − x1 )x2 (1 − x2 ), (2x1 − 1)x2 (1 − x2 ) σ (x, t) = 16t . (2x2 − 1)x1 (1 − x1 ) The source function f can be derived using the above functions.
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
Fig. 7. Contour plots of the numerical primal variable at t = 0, 0.25, and 0.5.
Fig. 8. Contour plots of the first component of numerical flux variable at t = 0, 0.25, and 0.5.
Fig. 9. Contour plots of the second component of numerical flux variable at t = 0, 0.25, and 0.5. Table 7 Errors, convergence orders, and CPU times with δ = 0.75. ∆t
Mesh
||{y − yh , σ − σh }|| L 2 (0,T ;stab)
Order
CPU time (s)
1/8 1/16 1/32 1/64
16*16 32*32 64*64 128*128
1.72E−01 9.42E−02 4.71E−02 2.37E−02
– 0.83 1.00 0.99
1.18 10.7 83.48 673.16
115
116
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
Fig. 10. Contour plots of the numerical primal variable at t = 0, 0.5, and 1.0.
Fig. 11. Contour plots of the first component of numerical flux variable at t = 0, 0.5, and 1.0. Table 8 Errors, convergence orders, and CPU times with δ = 0.75. ∆t
Mesh
||{y − yh , σ − σh }|| L 2 (0,T ;stab)
Order
CPU time (s)
1/8 1/16 1/32 1/64
8*8 16*16 32*32 64*64
1.92E−01 9.22E−02 4.51E−02 2.21E−02
– 1.06 1.03 1.03
1.63 12.93 101.06 811.12
In Table 8 we present the errors and main CPU time for the stabilized mixed finite element method. Besides, in Figs. 10–12 the contours of numerical primal and flux variables are also plotted at t = 0, 0.5, and 1.0 respectively for h = 1/64. From these figures we can see the changing of height for the numerical solutions as the time marches. Both the errors and figures show that the stabilized mixed method is efficient and it matches very well with the theoretical analysis. 6. Concluding remarks In this paper, we propose a stabilized mixed finite element discretization of both steady and unsteady reaction–diffusion equations. Some a priori error estimates are discussed. Compared to the classical mixed finite element method, the new method has a free choice of mixed element function spaces, and it leads to a positive definite linear algebraic equation which is easily be solved by the popular solvers. Finally, numerical tests also show the efficiency and reliability of the proposed stabilized mixed finite element method. We shall consider the application of the new method in other equations, such as Stokes equation. Also a posteriori error estimates shall be discussed in the future.
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
117
Fig. 12. Contour plots of the second component of numerical flux variable at t = 0, 0.5, and 1.0.
Acknowledgments The authors would like to thank the editor and the anonymous referee for their valuable comments and suggestions on an earlier version of this paper. This work was supported by the National Natural Science Foundation of China (Nos. 11201485, 11571367), the Promotive Research Fund for Excellent Young and Middle-aged Scientists of Shandong Province (No. BS2013NJ001), and the Fundamental Research Funds for the Central Universities (Nos. 14CX02217A, 15CX08004A). The first author was also supported by the China Scholarship Council (No. 201506455014). References [1] P.A. Raviart, J.M. Thomas, A mixed finite element method for 2nd order elliptic problems, in: Lecture Notes in Mathematics, vol. 606, Springer, Berlin, 1977, pp. 292–315. [2] J. Douglas, J.E. Roberts, Global estimates for mixed methods for second order elliptic equations, Math. Comp. 44 (1985) 39–52. [3] F. Brezzi, J. Douglas Jr., L.D. Marini, Two family of mixed finite elements for second order elliptic problems, Numer. Math. 47 (1985) 217–235. [4] F. Brezzi, J. Douglas Jr., R. Duran, L.D. Marini, Mixed finite elements for second order elliptic problems in three space variables, Numer. Math. 51 (1987) 237–250. [5] F. Brezzi, J. Douglas Jr., M. Fortin, L.D. Marini, Efficient rectangular mixed finite elements in two and three space variables, RAIRO Model. Math. Numer. Anal. 21 (1987) 581–603. [6] A.I. Pehlivanov, G.F. Carey, D. Lazarov, Least-squares mixed finite elements for second-order elliptic problems, SIAM J. Numer. Anal. 31 (1994) 1368–1377. [7] Z. Cai, R. Lazarov, T.A. Manteuffel, S.F. Mccormick, First-Order system least squares for second-order partial differential equations: Part I, SIAM J. Numer. Anal. 31 (1994) 1785–1799. [8] A.K. Pani, An H 1 -Galerkin mixed finite element method for parabolic partial differential equations, SIAM J. Numer. Anal. 35 (1998) 712–727. [9] D.P. Yang, A splitting positive defnite mixed element method for miscible displacement of compressible flow in porous media, Numer. Methods Partial Differential Equations 17 (2001) 229–249. [10] H. Rui, S. Kim, S.D. Kim, A remark on least-squares mixed element methods for reaction–diffusion problems, J. Comput. Appl. Math. 202 (2007) 230–236. [11] T.J.R. Hughes, L.P. Franca, G. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations, Comput. Methods Appl. Mech. Engrg. 73 (1989) 173–189. [12] L.P. Franca, R. Stenberg, Error analysis of Galerkin least squares methods for the elasticity equations, SIAM J. Numer. Anal. 28 (1991) 1680–1697. [13] L.P. Franca, S.L. Frey, T.J.R. Hughes, Stabilized finite element methods: I. Application to the advective-diffusive model, Comput. Methods Appl. Mech. Engrg. 95 (1992) 253–276. [14] L.P. Franca, C. Farhat, Bubble functions prompt unusual stabilized finite element methods, Comput. Methods Appl. Mech. Engrg. 123 (1995) 299–308. [15] L.P. Franca, F. Valentin, On an improved unusual stabilized finite element method for the advective-reactive-diffusive equation, Comput. Methods Appl. Mech. Engrg. 190 (2000) 1785–1800. [16] H. Fu, H. Rui, J. Hou, H. Li, A stabilized mixed finite element method for elliptic optimal control problems, J. Sci. Comput. (2015) http://dx.doi.org/10.1007/s10915-015-0050-3. [17] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2002. [18] R. Li, W.B. Liu, http://circus.math.pku.edu.cn/AFEPack. [19] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, 1991.
ScienceDirect Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117 www.elsevier.com/locate/cma
A stabilized mixed finite element method for steady and unsteady reaction–diffusion equations Hongfei Fu a,∗ , Hui Guo a , Jian Hou b , Junlong Zhao a a College of Science, China University of Petroleum, Qingdao 266580, China b College of Petroleum Engineering, China University of Petroleum, Qingdao 266580, China
Received 18 June 2015; received in revised form 1 December 2015; accepted 15 January 2016 Available online 12 February 2016
Abstract In this paper, we propose a new mixed finite element method, called stabilized mixed finite element method, for the approximation of steady reaction–diffusion partial differential equations (PDEs). The method is obtained by translating the primal second-order PDEs into a first-order mixed system, and then adding some suitable elementwise residual terms multiplied by a stabilization parameter to the weak formulation. The new method is compatible, i.e., the added terms equal to zero in the continuous case. Furthermore, it is mesh-independent, i.e., the stabilization parameter is independent of the mesh size. We prove both coercive and continuous properties in a weighted norm for the corresponding new mixed bilinear formulation. These assure that the finite element function spaces do not require to satisfy the classical Ladyzhenkaya–Babuska–Brezzi (LBB) consistency condition. Therefore, the widely used Lagrange finite element can be adopted. A simple proof of a priori error estimate with lower order regularity requirement is discussed, and numerical experiments confirm the efficiency and reliability of the new stabilized mixed method. Finally, the method is applied to solving unsteady reaction–diffusion equations. Error estimates are also given, and numerical examples still support the theoretical analysis very well. c 2016 Elsevier B.V. All rights reserved. ⃝
MSC: 65N15; 65N30 Keywords: Reaction–diffusion equation; Stabilized mixed finite element; Classical mixed finite element; A priori error estimate; Numerical experiments
1. Introduction Let Ω be a bounded convex domain in Rd (d = 1, 2, 3) with Lipschitz boundary ∂Ω . In this paper, we are interested in proposing a stabilized mixed finite element method for the following linear reaction–diffusion equations: −div(A∇ y) + cy = f, in Ω , (1.1) y = 0, on ∂Ω . ∗ Corresponding author.
E-mail addresses: [email protected] (H. Fu), [email protected] (H. Guo), [email protected] (J. Hou). http://dx.doi.org/10.1016/j.cma.2016.01.010 c 2016 Elsevier B.V. All rights reserved. 0045-7825/⃝
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
103
Depending on a simplified model for particular real problems, the variable y can represent a temperature, a pressure, or a concentration of a chemical species. The matrix A = A(x) describes the diffusion feature, and c = c(x) stands for a reaction rate. Detailed description is in the following section. As in the real problems, the gradient of the primal variable is also meaningful. For example, in the flow problem, it stands for a velocity; or, in the temperature control problem, large temperature gradients during cooling or heating may lead to its destruction. Therefore, in these cases people pay their special attention on the flux of the primal variable. Let us introduce the flux σ = −A∇ y, then (1.1) is equivalent to the following mixed form in Ω , divσ + cy = f, (1.2) A−1 σ + ∇ y = 0, in Ω , y = 0, on ∂Ω . It is well known that (1.2) is widely studied by many authors using the classical mixed finite element method; see, for example, [1–5]. The most important merit of the mixed finite element method is that the flux σ can be approximated simultaneously to a same order of accuracy as the primal variable y, without any accuracy lost compared to the standard finite element method. However, the finite element function spaces for y and σ have to be specially chosen such that they strictly satisfy the LBB consistency condition. In this sense, some of the best known and widely used finite element spaces are thereby excluded. In particular, in high-dimensional problems, it is not a straightforward thing to construct such kinds of mixed finite element function spaces. Meanwhile, the practical computation for classical mixed finite element methods is complicated, people need to solve a saddle-point problem which is symmetric but indefinite. To overcome these difficulties, some stabilization techniques have been studied, for example, least-squares methods [6,7], H 1 -Galerkin methods [8], as well as splitting methods [9,10]. In this paper, different from the previous approaches, we shall construct a stabilized mixed finite element method for the problem under consideration. The method can be easily obtained by adding some suitable elementwise residual terms to a dual-mixed weak formulation. We prove both coercive and continuous properties in a weighted norm for the corresponding new mixed bilinear formulation. Thus, the new method can avoid the drawbacks in selecting finite element function spaces. It is also competitive in practical computation, as the coefficient matrix of the leading system of linear algebraic equation is positive definite, and there are fewer degree of freedoms (Dofs) compared to the corresponding classical mixed finite element method. The rest of the paper is organized as follows: In Section 2, starting from the classical mixed formulation, we first introduce the stabilized mixed weak formulation, and then prove the existence and uniqueness of its solution. Finally, discretization using standard Lagrange elements is discussed. In Section 3, a simple error analysis for the stabilized mixed finite element method is given. In Section 4, we conduct some numerical experiments to verify the theoretical analysis. We also compare the results of the new proposed stabilized mixed method with those obtained by the classical mixed method. In Section 5, we extent the stabilized method to solving unsteady reaction–diffusion problems. Error estimate and numerical test are also presented. Finally, some concluding remarks are given in the last section. 2. Stabilized mixed finite element method Let ω ⊂ Ω . Throughout this paper, we adopt the standard notations L p (ω) (1 ≤ p ≤ ∞) for Lebesgue space of real-valued functions with norm ∥ · ∥0, p,ω , and W m, p (ω) (1 ≤ p ≤ ∞) for Sobolev spaces endowed with norm ∥ · ∥m, p,ω and semi-norm | · |m, p,ω . For p = 2, we denote W m,2 (ω) = H m (ω) and we drop the subscript p = 2 in the corresponding norms and semi-norms. Besides, the subscript ω will also be omitted in the norms if ω = Ω . To give a detailed description of the reaction–diffusion problem under consideration in a mixed weak formulation, we shall take the function spaces V = L 2 (Ω ),
W = H (div; Ω ) , {τ ∈ L 2 (Ω )d : divτ ∈ L 2 (Ω )},
Y = H01 (Ω ),
Σ = L 2 (Ω )d .
and
Besides, in this paper we need the following assumptions. (H-1) The solution (y, σ ) of (1.2) satisfies the following regularity: y ∈ H k+1 (Ω ),
σ ∈ H r (Ω )d .
104
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
(H-2) A = (ai, j (x))d×d is symmetric, and there exist positive constants α and β such that α|X |2 ≤ X T AX ≤ β|X |2 ,
∀X ∈ Rd .
(H-3) c = c(x) is bounded below and up by two constants c∗ and c∗ such that 0 ≤ c∗ ≤ c ≤ c∗ . Now, let us first recall the classical mixed weak formulation of (1.2). It aims to find (y, σ ) ∈ V × W such that −1 (A σ, τ ) − (y, divτ ) = 0, ∀τ ∈ W, (2.1) (divσ, v) + (cy, v) = ( f, v), ∀v ∈ V or equivalently (A−1 σ, τ ) − (y, divτ ) + (divσ, v) + (cy, v) = ( f, v),
∀(v, τ ) ∈ V × W.
(2.2)
This idea is firstly proposed by Raviart and Thomas to second-order elliptic problems, see [1]. From then on, it is widely discussed and used to solve various equations. It is also extended to combine with other numerical methods, such as discontinuous Galerkin method, finite volume method, and so on. We cannot even briefly summarize it here. An alternative approach uses the concept of least-squares residual stabilization as proposed by Hughes et al. [11] for solving advection–diffusion equations, and theoretically analyzed by Franca and Stenberg [12] for the elasticity equations. They call the new method as Galerkin least-squares method. After then, an unusual stabilized finite element method, see [13–15] is discussed and analyzed. The difference between these two stabilized methods appears in the sign of second-order term in the test function, see [13] for details. Inspired by the previous work, we hope to adopt the structure of these stabilized finite element methods, and propose a stabilized mixed finite element method to conquer the stability problem in classical mixed finite element method. To do this, let us first rewrite Eq. (2.2) as follows: (A−1 σ, τ ) + (∇ y, τ ) − (σ, ∇v) + (cy, v) = ( f, v),
∀(v, τ ) ∈ Y × Σ .
(2.3) (A−1 σ
Then we augment the dual-mixed formulation (2.3) by a non-symmetric residual term + ∇ y, τ − A∇v)e on each element e to propose a stabilized mixed method. In this sense, the new mixed weak formulation of (1.2) reads as finding (y, σ ) ∈ Y × Σ such that aδ (y, σ ; v, τ ) = ( f, v),
∀(v, τ ) ∈ Y × Σ .
(2.4)
Here the new bilinear form aδ (y, σ ; v, τ ) , (A−1 σ, τ ) + (∇ y, τ ) − (σ, ∇v) + (cy, v) −
e∈T
δ(A−1 σ + ∇ y, τ − A∇v)e ,
(2.5)
h
where δ > 0 is an elementwise constant parameter, (·, ·)e denotes the inner product in L 2 (e)d , and the element e belongs to a regular triangulation T h of the given domain Ω , which is defined below. To prove the existence and uniqueness of the solution pair to the new stabilized mixed weak formulation (2.4), we first define a weighted δ- stabilization norm corresponding to the bilinear form aδ (, ; , ) as follows: |||{y, σ }|||2stab = (A−1 σ, σ ) + δ(A∇ y, ∇ y)e + (cy, y). (2.6) e∈T h
Then, we can derive the following coercive and bounded results which are proved similarly in Ref. [16], but with only one stabilization parameter here. For readers’ convenience, we still include the proof here. Proposition 2.1. For 0 < δ0 ≤ δ ≤ δ1 < 1, we have aδ (y, σ ; y, σ ) ≥ cδ |||{y, σ }|||2stab , Moreover, for all (y, σ ; v, τ ) ∈ (Y
∀(y, σ ) ∈ Y × Σ .
(2.7)
× Σ )2
aδ (y, σ ; v, τ ) ≤ Cδ |||{y, σ }|||stab |||{v, τ }|||stab . Here δi (i = 0, 1) are two positive constants, and cδ = 1 − δ1 and Cδ = max{4, 2 + 1/δ0 }.
(2.8)
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
Proof. On the one hand, it follows from the definition that aδ (y, σ ; y, σ ) = (A−1 σ, σ ) + (cy, y) − δ(A−1 σ, σ )e + δ(A∇ y, ∇ y)e . e∈T h
105
(2.9)
e∈T h
Hence, for δ ≤ δ1 < 1, the bilinear form a(, ; , ) is coercive on Y × Σ , and satisfy (2.7). On the other hand, it follows from the definition that aδ (y, σ ; v, τ ) = (A−1 σ, τ ) + (∇ y, τ ) − (σ, ∇v) + (cy, v) − δ(A−1 σ, τ )e e∈T h
+
δ(σ, ∇v)e −
e∈T h
δ(∇ y, τ )e +
e∈T h
δ(A∇ y, ∇v)e .
(2.10)
e∈T h
Note that i
|ai bi | ≤
ai2
1 2
i
1
2
.
bi2
i
Then for 0 < δ0 ≤ δ ≤ δ1 < 1 we have 1 2 aδ (y, σ ; v, τ ) ≤ 4(A−1 σ, σ ) + (A∇ y, ∇ y) + 2 δ(A∇ y, ∇ y)e + (cy, y) e∈T h
1 2 4(A−1 τ, τ ) + (A∇v, ∇v) + 2 δ(A∇v, ∇v)e + (cv, v) e∈T h
−1
≤ 4(A
1 σ, σ ) + 2 + δ0
δ(A∇ y, ∇ y)e + (cy, y)
21
e∈T h
21 1 4(A−1 τ, τ ) + 2 + δ(A∇v, ∇v)e + (cv, v) δ0 h e∈T
≤ Cδ |||{y, σ }|||δ |||{v, τ }|||δ , which proves the bounded result (2.8).
(2.11)
Proposition 2.2. Let f ∈ L 2 (Ω ). Under assumptions (H-2)–(H-3) , Problem (2.4) admits a unique solution pair (y, σ ) ∈ Y × Σ . Moreover, the solution exhibits the improved regularity (y, σ ) ∈ H 2 (Ω ) × H 1 (Ω )d , and satisfies the following stability estimate ∥y∥2 + ∥σ ∥1 ≤ C∥ f ∥.
(2.12)
Proof. Proposition 2.1 implies that the bilinear form a(, ; , ) is coercive and bounded. Thus the Lax–Milgram lemma implies the existence and uniqueness of the solution pair (y, σ ) ∈ Y × Σ to problem (2.4). Furthermore, stability result (2.12) follows from the elliptic regularity depending on the shape of the domain. Now we consider the finite element discretization. Let T h be a family of regular partition of Ω , such that Ω = ∪e∈T h e. Let h e denote the diameter of the element e ∈ T h , and h = maxe∈T h h e . As the bilinear form a(, ; , ) is coercive and bounded, therefore we can choose the standard Lagrange element as the flux function spaces. That is to say, the LBB consistency condition which is required by classical mixed element spaces, for example, RT elements [1], now is not necessary. For two integers k ≥ 1, r ≥ 0, we construct the finite element subspaces Y h ⊂ Y and Σ h ⊂ Σ in a standard way: Y h = {vh ∈ C 0 (Ω ) : vh |e ∈ Pk (e), ∀e ∈ T h , vh = 0 on ∂Ω }, Σ h = {τh ∈ L 2 (Ω )d : (τh )i |e ∈ Pr (e), i = 1, 2, . . . , d, ∀e ∈ T h }. Here Pm denotes polynomials of total degree at most m. We can also adopt r th order RT elements as Σ h .
(2.13)
106
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
Then the discretization of the state equation (2.4) reads as follows: Find (yh , σh ) ∈ Y h × Σ h such that aδ (yh , σh ; vh , τh ) = ( f, vh ),
∀(vh , τh ) ∈ Y h × Σ h .
(2.14)
Remark 2.3. Note that the new stabilized mixed finite element method is compatible. Besides, it satisfies the following orthogonality relation: aδ (y − yh , σ − σh ; vh , τh ) = 0,
∀(vh , τh ) ∈ Y h × Σ h .
(2.15)
Remark 2.4. It is well known that the discretization of (2.1) using the classical mixed finite element spaces, see, for example, RT element [1], BDDM element [4], BDFM element [5], BDM element [3], will lead to a typical saddle-point problem. Therefore, the popular conjugate gradient or algebraic multigrid solvers cannot be adopted. However, the new stabilized mixed finite element method used in (2.14) not only has a free choice of finite element function spaces (see (2.13)), but also produces a positive definite linear algebraic equation (see Proposition 2.1). Therefore, it can be solved rapidly and efficiently with less degree of freedoms (Dofs) compared to the classical mixed method. We will illustrate an iteration algorithm for decoupling the finite element solutions yh and σh of (2.14) in Section 4. 3. A priori error estimates In this section, we aim to derive optimal a priori error estimate for the stabilized mixed finite element approximation of steady reaction–diffusion equations. In the following, let us first recall two important results without proof. Lemma 3.1. Let Ih be the standard Lagrange interpolation operator defined in Ref. [17]. Then there is a constant C > 0 such that ∥v − Ih v∥0,e + h e |v − Ih v|1,e ≤ Ch m e |v|m,e , for any v ∈
H m (e)
and e ∈ T
(3.1)
h.
Theorem 3.2. Let (y, σ ) and (yh , σh ) be the solutions of (2.4) and (2.14), respectively. Assume that the solutions y ∈ H01 (Ω ) ∩ H k+1 (Ω ) and σ ∈ H r (Ω )d . Then the following estimates hold (3.2) |||{y − yh , σ − σh }|||stab ≤ Ch min(k,r ) ∥y∥ H k+1 (Ω ) + ∥σ ∥ H r (Ω )d . Proof. Following from the orthogonality property (2.15), the theorem can be easily proved by the interpolation estimation theory in Lemma 3.1. Let Ih y ∈ Y h and Ih σ ∈ Σ h be defined as the kth and r th order Lagrange interpolations of y and σ , respectively. Then we conclude from (2.15) that aδ (Ih y − yh , Ih σ − σh ; vh , τh ) = −aδ (y − Ih y, σ − Ih σ ; vh , τh ).
(3.3)
Let vh = Ih y − yh and τh = Ih σ − σh in Eq. (3.3). It then follows from Proposition 2.1 and Cauchy–Schwarz inequality that cδ |||{Ih y − yh , Ih σ − σh }|||2δ ≤ Cδ |||{y − Ih y, σ − Ih σ }|||δ |||{Ih y − yh , Ih σ − σh }|||δ .
(3.4)
This implies that |||{Ih y − yh , Ih σ − σh }|||stab ≤ C|||{y − Ih y, σ − Ih σ }|||stab .
(3.5)
Thus we deduce by the triangle inequality and Lemma 3.1 that |||{y − yh , σ − σh }|||stab ≤ |||{y − Ih y, σ − Ih σ }|||stab + |||{Ih y − yh , Ih σ − σh }|||stab ≤ C|||{Ih y − y, Ih σ − σ }|||stab ≤ Ch min(k,r ) ∥y∥ H k+1 (Ω ) + ∥σ ∥ H r (Ω )d , which proves Theorem 3.2.
(3.6)
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
107
4. Numerical experiments In this section, we are going to carry out some numerical experiments to numerically support the a priori error estimates derived in the last section. To solve the problem numerically, we use a C++ software package: AFEPack, see Ref. [18] for the details. Let Ω = [0, 1] × [0, 1] and the stabilized parameter δ = 0.75 in the whole domain Ω . We consider piecewise linear C 0 elements for both the approximations of the primal variable y and its flux σ for the new proposed stabilized mixed finite element method (2.14). To show more convincingly the benefits of the stabilized mixed method, we also compare the results with those obtained using the classical mixed finite element approximation of (2.1). As the classical mixed method is only stable for special choice of matched mixed finite element spaces, we consider here the lowest-order RT elements (see [19,1]) for the approximation of the flux σ , while matched piecewise constant elements are employed for the approximation of the primal variable y. In this part, we consider the following iteration algorithm for solving the coupled mixed system. Algorithm. Step 1. Give an initial value σh0 and tolerance TOL = 1.0E−05. Step 2. Solve yhm ∈ V h such that δ(A∇ yhm , ∇vh ) + (cyhm , vh ) = ( f, vh ) + (1 − δ)(A−1 σhm−1 , ∇vh ). Step 3. Solve σhm ∈ W h such that (A−1 σhm , τh ) = −(∇ yhm , τh ). Step 4. Stop if ∥σhm − σhm−1 ∥ L 2 (Ω )2 < TOL is satisfied, and let yh = yhm , σh = σhm . Otherwise set m := m + 1 and go to Step 2 again.
Example 4.1. For the first example, we consider a diffusion-dominated reaction–diffusion problem with coefficients A = I and c = 1.0E−03. Let the solutions be as follows: y(x) = sin(πx1 ) sin(π x2 ), cos(π x1 ) sin(π x2 ) σ (x) = −π . sin(π x1 ) cos(π x2 ) The source function f can be computed using the above functions. Table 1 displays the weighted errors with respect to the decreasing mesh size. We can observe that the convergence orders are consistent with Theorem 3.2. We also show the main CPU time for the iteration algorithm in Table 1. It can be found that 52.26 s of CPU time is employed when h = 1/128 for the stabilized mixed finite element method using an AMG (algebraic multigrid method) solver. However, to obtain a same stop criterion TOL = 1.0E−05, we can observe from Table 2 that the classical mixed finite element method using a preconditioned GMRES (generalized minimal residual method) solver costs 279.08 s of CPU time, and at the meantime the error |||{y − yh , σ − σh }||| RT := ∥y − yh ∥+∥σ −σh ∥ H (div;Ω ) = 7.18E−02. We should notice that here a preconditioner is adopted for fast computation of the nonsymmetric algebraic equation systems result from the classical mixed method. Otherwise, the classical mixed method shall spend more CPU time for a fine mesh size. Figs. 1–2 show the profiles of the numerical primal variable and its flux with h = 1/128 using the stabilized mixed method. It can be seen that the numerical figures are also agreement with the analytical solutions very well. Example 4.2. For the second example, we consider a reaction–diffusion problem with coefficients A = I and c = 1, where the reaction and diffusion are of the same importance. Let the solutions be as follows: y(x) = x1 (1 − x1 )x2 (1 − x2 ) exp(x1 − x2 ), (2x1 − 1)x2 (1 − x2 ) −1 σ (x) = exp(x1 − x2 ) +y . (2x2 − 1)x1 (1 − x1 ) 1 The source function f can be derived using the above functions.
108
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
Fig. 1. The numerical primal variable with h = 1/128.
Table 1 Results for the stabilized mixed finite element method with δ = 0.75. Mesh
||{y − yh , σ − σh }|| δ
Order
CPU time (s)
16*16 32*32 64*64 128*128
1.34E−01 6.70E−02 3.34E−02 1.67E−02
– 1.00 1.00 1.00
0.73 3.14 13.06 52.26
Table 2 Results for the classical mixed finite element method. Mesh
||{y − yh , σ − σh }|| RT
Order
CPU time (s)
16*16 32*32 64*64 128*128
5.73E−01 2.87E−01 1.43E−01 7.18E−02
– 1.00 1.01 0.99
0.19 1.51 14.41 279.08
Table 3 Results for the stabilized mixed finite element method with δ = 0.75. Mesh
||{y − yh , σ − σh }|| δ
Order
CPU time (s)
16*16 32*32 64*64 128*128
1.34E−02 6.66E−03 3.29E−03 1.64E−03
— 1.01 1.02 1.00
0.6 2.45 9.79 40.12
We present the corresponding errors and CPU times for the stabilized/classical mixed finite element methods with respect to different mesh sizes in Tables 3 and 4. Numerical figures of the primal and flux variables using the stabilized mixed method are shown in Figs. 3 and 4. Upon the same stop criterion TOL = 1.0E−05, the classical mixed finite element method spends 80.49 s of CPU time, which is about twice of the stabilized mixed finite element method for the mesh size h = 1/128. It is clear that the new stabilized mixed method is much more efficient for large computation.
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
109
Fig. 2. The numerical flux variable with h = 1/128.
0
Fig. 3. The numerical primal variable with h = 1/128.
Table 4 Results for the classical mixed finite element method. Mesh
||{y − yh , σ − σh }|| RT
Order
CPU time (s)
16*16 32*32 64*64 128*128
3.64E−02 1.83E−02 9.29E−03 5.72E−03
– 0.99 0.98 0.70
0.18 1.31 9.95 80.49
Example 4.3. For the third example, we consider a reaction-dominated reaction–diffusion problem with coefficients A = ε I and c = 1. Let ε = 4.0E−04 and the solution pair as follows: (x1 − 0.5)2 + (x2 − 0.5)2 y(x) = exp − , √ ε √ x − 0.5 σ (x) = 2 εy 1 . x2 − 0.5
110
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
0
0
Fig. 4. The numerical flux variable with h = 1/128.
Table 5 Results for the stabilized mixed finite element method with δ = 0.75. Mesh
||{y − yh , σ − σh }|| δ
Order
CPU time (s)
16*16 32*32 64*64 128*128
7.60E−03 3.49E−03 1.71E−03 8.49E−04
– 1.12 1.03 1.01
0.27 1.2 4.72 18.82
Table 6 Results for the classical mixed finite element method. Mesh
||{y − yh , σ − σh }|| RT
Order
CPU time (s)
16*16 32*32 64*64 128*128
2.25E−02 1.13E−02 5.67E−03 2.84E−03
– 0.99 0.99 1.00
0.11 0.28 1.35 7.65
The source function f can be computed using the above functions. We present weighted errors and CPU times with respect to different mesh sizes in Table 5. Numerical figures of the primal and flux variables are shown in Figs. 5 and 6. This example provides a very sharp solution around the point (0.5, 0.5). However, the test results using the stabilized mixed finite element method also support the theoretical analysis in Theorem 3.2 very well. As this example presents a nearly zero solution in most area of the domain, the preconditioned GMRES solver embodies a high efficiency because of the sparsity of the leading algebraic matrix, and from Table 6 we can see that only 7.65 s is employed by the classical mixed finite element method with error 2.84E−03 when h = 1/128. However, we could still anticipate that when the mesh size is decreasing, the stabilized mixed method shall be more competitive compared to the classical one. As when h is getting smaller, there will be a more increasing degree of freedoms (Dofs) for the classical mixed method than that for the stabilized mixed method. In summary, numerical experiments in this section show the high efficiency and reliability of the stabilized mixed finite element method, which coincides with the theoretical analysis. In most cases, few CPU time and small errors are observed compared to the classical mixed finite element method even using an efficient preconditioned GMRES solver.
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
111
Fig. 5. The numerical primal variable with h = 1/128.
Fig. 6. The numerical flux variable with h = 1/128.
5. Application to unsteady reaction–diffusion equations 5.1. Preliminary and discretization In this section, we consider an extension of stabilized mixed finite element method for the following unsteady reaction–diffusion equations: yt − div(A∇ y) + cy = f, in Ω × (0, T ], y = 0, on ∂Ω × (0, T ], (5.1) y(x, 0) = y0 (x), in Ω where T < +∞ denotes the final time, and y0 (x) is a prescribed initial data. To define a fully-discrete scheme, we introduce a time partition of [0, T ] by t i = i∆t, i = 0, 1, 2, . . . , N with ∆t = T /N where N is a positive integer. Let f i = f (·, t i ). We denote by L s (0, T ; W m, p (Ω )) the Banach space of
112
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
all L s integrable functions from (0, T ] into W m, p (Ω ) with norm ∥v∥ L s (0,T ;W m, p (Ω )) =
1s
T
for s ∈ [1, ∞)
∥v∥m, p dt 0
and the standard modification for s = ∞. Similarly, let σ = −A∇ y, then the first equation in (5.1) is equivalent to the following mixed form yt + divσ + cy = f, in Ω × (0, T ], A−1 σ + ∇ y = 0, in Ω × (0, T ]. Therefore, a stabilized mixed weak formulation of (5.2) reads as finding (y, σ ) : (0, T ] → Y × Σ such that (yt , v) + aδ (y, σ ; v, τ ) = ( f, v), ∀(v, τ ) ∈ Y × Σ , y(0) = y0 .
(5.2)
(5.3)
Let Y h and Σ h be the finite element spaces as defined in (2.13). Then the discretization of Eq. (5.3) reads as follows: Find (yhi , σhi ) ∈ Y h × Σ h for i = 1, 2, . . . , N such that i i−1 yh − yh , v + a (y i , σ i ; v , τ ) = ( f i , v ), ∀(v , τ ) ∈ Y h × Σ h , h δ h h h h h h h ∆t (5.4) 0 yh = y0,h , where y0,h ∈ Y h is an approximation of y0 which will be specified later on. 5.2. A priori error estimate In this section, based on the error analysis for the steady problem in Section 3, we use an elliptic projection to give a priori error estimate for the unsteady reaction–diffusion problem. To this aim, we define ( yh , σh ) : (0, T ] → Y h × Σ h such that aδ (y − yh , σ − σh ; vh , τh ) = 0,
∀(vh , τh ) ∈ Y h × Σ h .
(5.5)
Then it follows from (5.5) and Theorem 3.2, we have the following auxiliary result. Lemma 5.1. Let (y, σ ) and ( yh , σh ) be the solutions of (5.3) and (5.5), respectively. For any t ∈ (0, T ], assume that the solutions y(t) ∈ H01 (Ω ) ∩ H k+1 (Ω ) and σ (t) ∈ H r (Ω )d . Then the following estimates hold |||{y − yh , σ − σh }|||stab ≤ Ch min(k,r ) ∥y∥ H k+1 (Ω ) + ∥σ ∥ H r (Ω )d . (5.6) In the following, we just need to estimate the bounds for yh − yh and σh − σh . It is summarized as the following lemma. Lemma 5.2. Let (yh , σh ) and ( yh , σh ) be the solutions of (5.4) and (5.5), respectively. Assume that yh (0) = yh (0) and the solutions y ∈ H 1 (0, T ; H k+1 (Ω )) ∩ H 2 (0, T ; L 2 (Ω )) and σ ∈ H 1 (0, T ; H r (Ω )d ). Then we have ∥ yh − yh ∥ L ∞ (0,T ;L 2 (Ω )) + |||{ yh − yh , σh − σh }||| L 2 (0,T ;stab) ≤ Ch min(k,r ) ∥yt ∥ L 2 (0,T ;H k+1 (Ω )) + ∥σt ∥ L 2 (0,T ;H r (Ω )d ) + C∆t∥ytt ∥2L 2 (0,T ;L 2 (Ω )) , where |||{ξ, ζ }|||2L 2 (0,T ;stab) :=
N
2
∆t|||{ξ i , ζ i }|||stab
i=1
is a time-dependent discrete norm.
(5.7)
113
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
Proof. For simplicity, we split the errors at time t = t i as follows: y i − yhi = (y i − yhi ) + ( yhi − yhi ) = ρ i + ξ i , σ i − σhi = (σ i − σhi ) + ( σhi − σhi ) = ηi + ζ i . It then follows from (5.3) to (5.4) and the auxiliary projection (5.5) that for any (vh , τh ) ∈ Y h × Σ h i i ξ − ξ i−1 ρ − ρ i−1 i i , vh + aδ (ξ , ζ ; vh , τh ) = − , vh − (χ i , vh ), ∆t ∆t
(5.8)
where χ i := yti −
y i − y i−1 . ∆t
Let the test functions vh = ξ i , τh = ζ i in Eq. (5.8). Note that i ξ − ξ i−1 i 1 i 2 ,ξ ≥ ∥ξ ∥ − ∥ξ i−1 ∥2 . ∆t 2∆t Then we conclude from Proposition 2.1 that the left-hand side of (5.8) is bounded below by 1 i 2 2 ∥ξ ∥ − ∥ξ i−1 ∥2 + cδ |||{ξ i , ζ i }|||stab . 2∆t
(5.9)
For the right-hand side of (5.8), it is a matter of calculation that i i t t i i−1 |ρ − ρ | = ρt dt ≤ |ρt | dt, t i−1 t i−1 i 1 ti t i (t − ti−1 )ytt dt ≤ |χ | = |ytt | dt. ∆t t i−1 t i−1 Then, following Cauchy–Schwarz inequality we have
ρ i − ρ i−1 i ,ξ ∆t
1 ≤ ∆t
ti
t i−1
∥ρt ∥dt · ∥ξ i ∥ ≤
1 ∥ρt ∥2L 2 (t i−1 ,t i ;L 2 (Ω )) + C∥ξ i ∥2 , ∆t
(5.10)
(χ i , ξ i ) ≤ ∆t∥ytt ∥2L 2 (t i−1 ,t i ;L 2 (Ω )) + C∥ξ i ∥2 .
(5.11)
The inequalities (5.10)–(5.11) can be combined with (5.9) to give the following recursion relation 1 i 2 2 ∥ξ ∥ − ∥ξ i−1 ∥2 + cδ |||{ξ i , ζ i }|||stab 2∆t 1 ∥ρt ∥2L 2 (t i−1 ,t i ;L 2 (Ω )) + ∆t∥ytt ∥2L 2 (t i−1 ,t i ;L 2 (Ω )) + C∥ξ i ∥2 . ≤ ∆t
(5.12)
Let the initial value y0,h = yh (0), i.e., ξ 0 = 0. Then if (5.12) is multiplied by 2∆t and summed over i from 1 to N , we have ∥ξ N ∥2 + C
N i=1
2
∆t|||{ξ i , ζ i }|||stab ≤ C∥ρt ∥2L 2 (0,T ;L 2 (Ω )) + C(∆t)2 ∥ytt ∥2L 2 (0,T ;L 2 (Ω )) + C
Then Lemma 5.2 follows from Lemma 5.1 and the discrete Gronwall’s lemma.
N
∆t∥ξ i ∥2 .
i=1
Finally, collecting the bounds in Lemmas 5.1, 5.2, and using the triangle inequality we can easily derive the following main theorem.
114
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
Theorem 5.3. Let (y, σ ) and (yh , σh ) be the solutions of (5.3) and (5.4), respectively. Assume that all conditions in Lemmas 5.1 and 5.2 are valid. Furthermore, suppose that y ∈ L ∞ (0, T ; H01 (Ω ) ∩ H k+1 (Ω )) and σ ∈ L ∞ (0, T ; H r (Ω )d ). Then the following estimate holds ∥y − yh ∥ L ∞ (0,T ;L 2 (Ω )) + |||{y − yh , σ − σh }||| L 2 (0,T ;stab) ≤ C(h min(k,r ) + ∆t).
(5.13)
Remark 5.4. To improve the convergence order in the temporal direction, one can adopt the following Crank–Nicolson stabilized mixed finite element discretization of (5.3) as yhi − yhi−1 1 i− 1 i− 1 , vh + aδ (yh 2 , σh 2 ; vh , τh ) = ( f i− 2 , vh ), ∀(vh , τh ) ∈ Y h × Σ h , (5.14) ∆t where yhi + yhi−1 , 2 Then one can expect i− 12
yh
=
i− 12
σh
=
σhi + σhi−1 . 2
∥y − yh ∥ L ∞ (0,T ;L 2 (Ω )) + |||{y − yh , σ − σh }||| L 2 (0,T ;stab) ≤ C(h min(k,r ) + (∆t)2 ).
(5.15)
5.3. Numerical tests In this subsection, we consider piecewise linear elements for both the approximations of y and σ , i.e., k = r = 1. Besides, an iteration algorithm as the steady case for solving the coupled mixed system is also employed, while the stop criterion is replaced by ∥σhN ,m − σhN ,m−1 ∥ L 2 (Ω )2 < TOL. Example 5.1. In this example, we numerically test a reaction-dominated reaction–diffusion problem on Ω = [−1, 1] × [−1, 1] and T = 0.5. Let the small diffusion coefficient A = ε I with ε = 4.0E−03, the reaction rate c = 1, and the sharp solutions be as follows: 1.0 − cos(t − x1 ) 1.0 − cos(t − x2 ) y(x, t) = exp − exp − ; ε ε sin(t − x1 ) σ (x, t) = −y . sin(t − x2 ) The source function f can be computed using the above functions. We list the time-dependent |||{y − yh , σ − σh }||| L 2 (0,T ;stab) errors and main CPU times with respect to different spatial mesh sizes and temporal sizes in Table 7. A first-order spatial and temporal convergence can be observed for the stabilized mixed method, which is in accordance with the theoretical analysis in Theorem 5.3 very well. It is well known that the solutions have very steep moving fronts, and when the time is marching the position of moving fronts is shifted along the straight line (x1 , x2 ) = (t, t). These phenomena can be observed clearly from Figs. 7–9, where the contour plots of numerical primal and flux variables are shown at t = 0, 0.25, and 0.5 respectively for h = 1/64. Furthermore, it is apparent that the maximum numerical value of the primal variable is indeed centered at (t, t) for different t. This is in accordance with the analytical solution very well. Example 5.2. For the second example, we consider a diffusion-dominated reaction–diffusion problem with coefficients A = I and c = 1.0E−03 on Ω = [0, 1] × [0, 1] and T = 1. Let the solutions be as follows: y(x, t) = 16t x1 (1 − x1 )x2 (1 − x2 ), (2x1 − 1)x2 (1 − x2 ) σ (x, t) = 16t . (2x2 − 1)x1 (1 − x1 ) The source function f can be derived using the above functions.
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
Fig. 7. Contour plots of the numerical primal variable at t = 0, 0.25, and 0.5.
Fig. 8. Contour plots of the first component of numerical flux variable at t = 0, 0.25, and 0.5.
Fig. 9. Contour plots of the second component of numerical flux variable at t = 0, 0.25, and 0.5. Table 7 Errors, convergence orders, and CPU times with δ = 0.75. ∆t
Mesh
||{y − yh , σ − σh }|| L 2 (0,T ;stab)
Order
CPU time (s)
1/8 1/16 1/32 1/64
16*16 32*32 64*64 128*128
1.72E−01 9.42E−02 4.71E−02 2.37E−02
– 0.83 1.00 0.99
1.18 10.7 83.48 673.16
115
116
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
Fig. 10. Contour plots of the numerical primal variable at t = 0, 0.5, and 1.0.
Fig. 11. Contour plots of the first component of numerical flux variable at t = 0, 0.5, and 1.0. Table 8 Errors, convergence orders, and CPU times with δ = 0.75. ∆t
Mesh
||{y − yh , σ − σh }|| L 2 (0,T ;stab)
Order
CPU time (s)
1/8 1/16 1/32 1/64
8*8 16*16 32*32 64*64
1.92E−01 9.22E−02 4.51E−02 2.21E−02
– 1.06 1.03 1.03
1.63 12.93 101.06 811.12
In Table 8 we present the errors and main CPU time for the stabilized mixed finite element method. Besides, in Figs. 10–12 the contours of numerical primal and flux variables are also plotted at t = 0, 0.5, and 1.0 respectively for h = 1/64. From these figures we can see the changing of height for the numerical solutions as the time marches. Both the errors and figures show that the stabilized mixed method is efficient and it matches very well with the theoretical analysis. 6. Concluding remarks In this paper, we propose a stabilized mixed finite element discretization of both steady and unsteady reaction–diffusion equations. Some a priori error estimates are discussed. Compared to the classical mixed finite element method, the new method has a free choice of mixed element function spaces, and it leads to a positive definite linear algebraic equation which is easily be solved by the popular solvers. Finally, numerical tests also show the efficiency and reliability of the proposed stabilized mixed finite element method. We shall consider the application of the new method in other equations, such as Stokes equation. Also a posteriori error estimates shall be discussed in the future.
H. Fu et al. / Comput. Methods Appl. Mech. Engrg. 304 (2016) 102–117
117
Fig. 12. Contour plots of the second component of numerical flux variable at t = 0, 0.5, and 1.0.
Acknowledgments The authors would like to thank the editor and the anonymous referee for their valuable comments and suggestions on an earlier version of this paper. This work was supported by the National Natural Science Foundation of China (Nos. 11201485, 11571367), the Promotive Research Fund for Excellent Young and Middle-aged Scientists of Shandong Province (No. BS2013NJ001), and the Fundamental Research Funds for the Central Universities (Nos. 14CX02217A, 15CX08004A). The first author was also supported by the China Scholarship Council (No. 201506455014). References [1] P.A. Raviart, J.M. Thomas, A mixed finite element method for 2nd order elliptic problems, in: Lecture Notes in Mathematics, vol. 606, Springer, Berlin, 1977, pp. 292–315. [2] J. Douglas, J.E. Roberts, Global estimates for mixed methods for second order elliptic equations, Math. Comp. 44 (1985) 39–52. [3] F. Brezzi, J. Douglas Jr., L.D. Marini, Two family of mixed finite elements for second order elliptic problems, Numer. Math. 47 (1985) 217–235. [4] F. Brezzi, J. Douglas Jr., R. Duran, L.D. Marini, Mixed finite elements for second order elliptic problems in three space variables, Numer. Math. 51 (1987) 237–250. [5] F. Brezzi, J. Douglas Jr., M. Fortin, L.D. Marini, Efficient rectangular mixed finite elements in two and three space variables, RAIRO Model. Math. Numer. Anal. 21 (1987) 581–603. [6] A.I. Pehlivanov, G.F. Carey, D. Lazarov, Least-squares mixed finite elements for second-order elliptic problems, SIAM J. Numer. Anal. 31 (1994) 1368–1377. [7] Z. Cai, R. Lazarov, T.A. Manteuffel, S.F. Mccormick, First-Order system least squares for second-order partial differential equations: Part I, SIAM J. Numer. Anal. 31 (1994) 1785–1799. [8] A.K. Pani, An H 1 -Galerkin mixed finite element method for parabolic partial differential equations, SIAM J. Numer. Anal. 35 (1998) 712–727. [9] D.P. Yang, A splitting positive defnite mixed element method for miscible displacement of compressible flow in porous media, Numer. Methods Partial Differential Equations 17 (2001) 229–249. [10] H. Rui, S. Kim, S.D. Kim, A remark on least-squares mixed element methods for reaction–diffusion problems, J. Comput. Appl. Math. 202 (2007) 230–236. [11] T.J.R. Hughes, L.P. Franca, G. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations, Comput. Methods Appl. Mech. Engrg. 73 (1989) 173–189. [12] L.P. Franca, R. Stenberg, Error analysis of Galerkin least squares methods for the elasticity equations, SIAM J. Numer. Anal. 28 (1991) 1680–1697. [13] L.P. Franca, S.L. Frey, T.J.R. Hughes, Stabilized finite element methods: I. Application to the advective-diffusive model, Comput. Methods Appl. Mech. Engrg. 95 (1992) 253–276. [14] L.P. Franca, C. Farhat, Bubble functions prompt unusual stabilized finite element methods, Comput. Methods Appl. Mech. Engrg. 123 (1995) 299–308. [15] L.P. Franca, F. Valentin, On an improved unusual stabilized finite element method for the advective-reactive-diffusive equation, Comput. Methods Appl. Mech. Engrg. 190 (2000) 1785–1800. [16] H. Fu, H. Rui, J. Hou, H. Li, A stabilized mixed finite element method for elliptic optimal control problems, J. Sci. Comput. (2015) http://dx.doi.org/10.1007/s10915-015-0050-3. [17] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2002. [18] R. Li, W.B. Liu, http://circus.math.pku.edu.cn/AFEPack. [19] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, 1991.