SIAM J. NUMER. ANAL. Vol. 35, No. 2, pp. 712–727, April 1998
c 1998 Society for Industrial and Applied Mathematics
017
AN H 1 -GALERKIN MIXED FINITE ELEMENT METHOD FOR PARABOLIC PARTIAL DIFFERENTIAL EQUATIONS∗ AMIYA K. PANI† Abstract. In this paper, an H 1 -Galerkin mixed finite element method is proposed and analyzed for parabolic partial differential equations with nonselfadjoint elliptic parts. Compared to the standard H 1 -Galerkin procedure, C 1 -continuity for the approximating finite dimensional subspaces can be relaxed for the proposed method. Moreover, it is shown that the finite element approximations have the same rates of convergence as in the classical mixed method, but without LBB consistency condition and quasiuniformity requirement on the finite element mesh. Finally, a better rate of convergence for the flux in L2 -norm is derived using a modified H 1 -Galerkin mixed method in two and three space dimensions, which confirms the findings in a single space variable and also improves upon the order of convergence of the classical mixed method under extra regularity assumptions on the exact solution. Key words. parabolic initial and boundary value problem, mixed finite element method, H 1 Galerkin, LBB condition, elliptic projection, semidiscrete scheme, backward Euler’s method, error estimates, Gronwall’s lemma AMS subject classifications. 65M60, 65M15, 65M12 PII. S0036142995280808
1. Introduction. In [9], [22], and [30], an H 1 -Galerkin method which requires C - continuity on the finite element spaces is discussed for parabolic problems. Apart from difficulties in constructing such spaces, computationally more attractive and widely used piecewise linear elements are excluded from this class of finite dimensional spaces. In order to relax C 1 -smoothness and use C 0 -elements, we first split the problem into a first-order system and then propose a nonsymmetric version of a least square method that is an H 1 -Galerkin procedure for the solution and its flux. The present method which we call an H 1 -Galerkin mixed finite element procedure even yields a better rate of convergence for the flux than the conventional use of linear elements. Earlier mixed finite element methods are developed and analyzed in [1], [2], [11] for elliptic equations, [7], [17] for parabolic problems, and [8], [15] for wave equations. Moreover, there is a growing list of literature on mixed methods in oil reservoir engineering problems; see [10], [12], [13], and [28] for references therein. However, this procedure has to satisfy LBB consistency condition on the approximating subspaces and restricts the choice of finite element spaces. For example, the Raviart–Thomas spaces of index r ≥ 1 are usually used for standard mixed methods. In this paper, an H 1 -Galerkin method is applied to a mixed system in p and its flux u. The approximating finite element spaces Vh and Wh are allowed to be of differing polynomial degrees. Hence, estimations have been obtained which distinguish the better approximation properties of Vh and Wh . Compared to standard mixed methods, the proposed one is not subject to LBB consistency condition. Although we require extra regularity on the solution, a better order of convergence for the flux 1
∗ Received
by the editors February 1, 1995; accepted for publication (in revised form) February 3, 1997. http://www.siam.org/journals/sinum/35-2/28080.html † Department of Mathematics, Indian Institute of Technology, Bombay, Powai, Mumbai- 400 076, India (
[email protected]). Present address: Department of Mathematics, Universidade Federal Do Parana, Centro Polit´ecnico, Cx.P. 19081, CEP: 81531-990, Curitiba, PR Brazil. 712
H 1 -MIXED METHOD
713
in L2 -norm is obtained (see Remarks 2.1, 4.1 and section 6). Moreover, we note that for L2 and H 1 -error estimates we do not impose quasiuniformity condition on the finite element mesh. The numerical experiments conforming the theoretical rate of convergence and comparision with classical mixed methods will be reported in a separate paper. Recently, an H 1 -Galerkin mixed finite element method was introduced and a brief discussion on the error analysis is presented for the unidimensional single phase Stefan problem in [21]. Subsequently, Pani and Anderssen [23] have extended this method to a parameter estimation problem in parabolic equations without rigorous error analysis. Related studies on least square mixed methods for elliptic equations can be found in [4], [5], [6], [14], [16], [18], [19], [20], [24], [25], and references therein. In section 2, a parabolic equation in a single space variable is considered and optimal error estimates are discussed for the semidiscrete case. In two and three space dimensions, a similar analysis is carried out in section 3. Moreover, the rates of convergence obtained coincide with those in Johnson and Thom´ee [17] using classical mixed method and the Raviart–Thomas finite element spaces, but compared to a one-dimensional case, the results in this section are not sharp. Therefore, in section 4, based on Neittaanm¨ aki and Saranen [18], [19], [20] a modification of H 1 Galerkin mixed finite element method is analyzed and semidiscrete error estimates are established. In section 5, a backward Euler scheme is briefly described along with a priori error bounds for the modified H 1 -Galerkin mixed method. Finally in section 6, extensions with concluding remarks are discussed. Throughout this paper, C will denote a generic positive constant which does not depend on h. 2. Parabolic problem in one space dimension. We consider the following one-dimensional parabolic partial differential equation: (2.1)
pt − (apx )x + bpx + cp = f (x, t) (x, t) ∈ (0, 1) × J,
with Dirichlet boundary condition p(0, t) = p(1, t) = 0, t ∈ J,
(2.2) and initial condition (2.3)
p(x, 0) = p0 , x ∈ I = (0, 1),
where pt = ∂p/∂t, px = ∂p/∂x and J = (0, T ] with T < ∞. The coefficients a, b, c are smooth functions of x and a is bounded below by a positive constant say a0 . For H 1 -Galerkin mixed finite element procedure, we first split the parabolic equation into a first-order system. Now introduce apx = u, and set α(x) = 1/a(x), β(x) = α(x)b(x). Then (2.1) can be written in the form of a first-order system as (2.4)
px = α(x)u,
pt − ux + βu + cp = f.
Denoting (·, ·) by natural inner product in L2 (I), let H01 = {v ∈ H 1 (I) : v(0) = v(1) = 0}. Further, we shall use the classical Sobolev spaces W m,p (I), 1 ≤ p ≤ ∞ and call it as W m,p with norm k · km,p . For the purpose of H 1 -Galerkin method, we now consider the following weak formulation: ‘Find a pair {p, u} : [0, T ] → H01 × H 1 such that (2.5) (2.6)
(px , vx ) = (α(x)u, vx ) , v ∈ H01 , (αut , w) + (ux , wx ) = (βu, wx ) + (cp, wx ) − (f, wx ) ,
w ∈ H 1 .’
714
AMIYA K. PANI
For the first term in (2.6), we have used integration by parts and the Dirichlet boundary conditions pt (0) = pt (1) = 0. Let Vh and Wh be finite dimensional subspaces of H01 and H 1 , respectively, with the following approximation properties: For 1 ≤ p ≤ ∞, and k > 0, r > 0 integers, inf kv − vh kLp (I) + hkv − vh kW 1,p (I) ≤ Chk+1 kvkW k+1,p , v ∈ H01 ∩ W k+1,p (I), vh ∈Vh
and inf
wh ∈Wh
kw − wh kLp (I) + hkw − wh kW 1,p (I) ≤ Chr+1 kwkW r+1,p , w ∈ W r+1,p (I).
The semidiscrete H 1 -Galerkin mixed finite element method for (2.4) is determined by a pair {ph , uh } ∈ Vh × Wh such that (2.7) (2.8)
(phx , vhx ) = (α(x)uh , vhx ) , vh ∈ Vh , (αuht , wh ) + (uhx , whx ) = (βuh , whx ) + (cph , whx ) − (f, whx ) ,
w ∈ Wh
with given (ph (0), uh (0)). The above equations yield a system of differential algebraic equations. This is of index one, since the stiffness matrix associated with (phx , vhx ) is positive definite. Therefore, the system (2.7)–(2.8) is uniquely solvable for a consistent initial condition, see [3]. Following Wheeler [31], we now define elliptic projections {˜ uh , p˜h } : [0, T ] → Wh × Vh as (2.9)
A(u − u ˜h , wh ) = 0,
∀wh ∈ Wh ,
(2.10)
(px − p˜hx , vhx ) = 0,
∀vh ∈ Vh ,
where A(u, v) = (ux , vx ) − (βu, vx ) + λ (u, v). Here λ is chosen approximately so that A is H 1 coercieve, i.e, A(v, v) ≥ α0 kvk21 , v ∈ H 1 . Moreover, it is not hard to check that A(., .) is bounded. Let ρ = u− u ˜h and η = p− p˜h . It is now quite standard to obtain the following estimate for ρ and η, (2.11)
kρkj + kρt kj ≤ Chr+1−j (kukr+1 + kut kr+1 ) , j = 0, 1,
and (2.12)
kηkj + kηt kj ≤ Chk+1−j (kpkk+1 + kpt kk+1 ) , j = 0, 1.
Further for j = 0, 1 and 1 ≤ p ≤ ∞, we have (2.13)
kρkW j,p (I) ≤ Chr+1−j kukW r+1,p (I) ,
and (2.14)
kηkW j,p (I) ≤ Chk+1−j kpkW k+1,p (I) .
Note that for p = ∞, we need the quasiuniformity condition on the finite element mesh. Semidiscrete error estimates. In the later part of this section, we shall discuss a priori error estimates for the semidiscrete Galerkin procedure. Let u − uh = (u − uh − uh ) = ρ + ξ and p − ph = (p − p˜h ) + (˜ ph − ph ) = η + ζ . Using (2.5)–(2.6), u ˜h ) + (˜
H 1 -MIXED METHOD
715
(2.7)–(2.8) and auxiliary projections (2.9)–(2.10), we now obtain error equations in ξ and ζ as (ζx , vhx ) = (α(x)ρ, vhx ) + (α(x)ξ, vhx ) , vh ∈ Vh ,
(2.15) and (2.16)
(αξt , wh ) + A (ξ, wh ) = − (α(x)ρt , wh ) + λ (ρ, wh ) + λ (ξ, wh ) + (cζ, whx ) + (cη, whx ) .
˜h (0). Then the error THEOREM 2.1. With u0 = ap0x , assume that uh (0) = u u − uh and p − ph can be estimated by the inequalities k(p − ph )(t)k + k(u − uh )(t)k + hk(p − ph )(t)k1 ≤ Chmin(k+1,r+1) [kpkL∞ (H k+1 ) + kukL∞ (H r+1 ) + kut kL2 (H r+1 ) ]. Proof. Since estimates of ρ and η can be found out from (2.11)–(2.12), it suffices to estimate ξ and ζ. Choose vh = ζ in (2.15) to have kζx k ≤ C(kρk + kξk).
(2.17)
Further set wh = ξ in (2.16) to obtain for small , (2.18)
1 d kα 2 ξk2 + (2α0 − )kξk21 ≤ C kρt k2 + kηk2 + kζk2 + kρk2 + Ckξk2 . dt
Since kζk ≤ kζx k, for ζ ∈ H01 , use (2.17) in (2.18) and then integrate with respect to time. An application of Gronwall’s lemma now yields Z t Z t kξ(s)k21 ds ≤ Ch2min(k+1,r+1) kp(s)k2k+1 + ku(s)k2r+1 + kut k2r+1 ds. kξk2 + 0
0
Using (2.11) and the above estimate in (2.17), we have from the Poincar´e inequality kζ(t)k ≤ kζ(t)x k ≤ Chmin(k+1,r+1) kpkL2 (H k+1 ) + kukL∞ (H r+1 ) + kut kL2 (H r+1 ) . Apply the triangle inequality with (2.11) and (2.12) to complete the rest of the proof. Note that it is possible to choose uh (0) as L2 projection of u0 onto Wh . ˜h (0). Then there exists a constant C > 0 THEOREM 2.2. Assume that uh (0) = u independent of h, such that k(u − uh )(t)k1 ≤Chmin(k+1,r) [kpkL∞ (H k+1 ) + kpt kL2 (H k+1 ) +kukL∞ (H r ) + kut kL2 (H r ) ], and for 1 < p ≤ ∞ k(p − ph )(t)kLp + k(u − uh )(t)kLp ≤ Chmin(k+1,r+1) [kpkL∞ (W k+1,p ) + kpt kL2 (H k+1 ) + kukL∞ (W r+1,p ) + kut kL2 (H r+1 ) ]. Proof. Since from Theorem 2.1, we have a superconvergence result for ζx , it is sufficient to obtain a superconvergence estimate for ξ in H 1 -norm. Choose wh = ξt in (2.16) to have 1
kα 2 ξt k2 +
1 d A (ξ, ξ) = − (α(x)ρt , ξt ) + λ (ρ, ξt ) + λ (ξ, ξt ) 2 dt d (cη, ξx ) − (cηt , ξx ) − ((cζ)x , ξt ) . + dt
716
AMIYA K. PANI
On integrating with respect to time from 0 to t and applying the standard Young’s inequality, we obtain Z t Z t kξk21 + kξ(s)t k2 ds ≤ C[kξ(t)k2 + kη(t)k2 + (kρ(s)k2 + kρt (s)k2 0 0 Z t 2 2 + kηt (s)k + kζx (s)k ) ds] + C kξ(s)k21 ds. 0
Using Theorem 2.1, an application of Gronwall’s Lemma now yields " Z t 2 2min(k+1,r+1) ku(s)k2r+1 kp(t)k2k+1 + ku(t)k2r+1 + kξk1 ≤ Ch 0
+
kpt (s)k2k+1
+
#
kut k2r+1 )ds
.
For H 1 -estimate, it is possible to reduce the regularity of u and ut . Since for 1 ≤ p ≤ ∞, we have kξ(t)kLp ≤ kξ(t)k1 , ξ ∈ H 1 and kζ(t)kLp ≤ kζx (t)k, ζ ∈ H01 , apply the superconvergence result with estimates (2.13)–(2.14) and then the triangle inequality completes the rest of the proof. Remark 2.1. (i) From Theorem 2.1 with k = r kp − ph kL∞ (L2 ) + ku − uh kL∞ (L2 ) ≤ Chr+1 [kpkL∞ (H r+2 ) + kpkL2 (H r+2 ) ]. Further, from the superconvergence results for k˜ ph − ph k1 and the estimate for kp − p˜h k1 , we have k(p − ph )(t)k1 ≤ Chmin(k,r+1) kpkL∞ (H k ) + kukL2 (H r+1 ) + kut kL2 (H r+1 ) . For k = r + 1, we obtain kp − ph kL∞ (H 1 ) = O(hk ), whereas ku − uh kL∞ (H 1 ) = O(hr ). Hence, the order of convergence corresponds to the degree of the polynomials used in the corresponding finite element spaces. For the one-dimensional selfadjoint elliptic equation, similar estimates are obtained by Pehlivanov et al. [24]. (ii) For C 0 -Lagrange elements with k = 3 and r = 1, the classical mixed finite element method fails, whereas the present method converges with order of convergence O(h2 ) for kp − ph k1 and ku − uh k-norms. (iii) Note that the coupling between p and u is mainly through c-term. If c = 0, we obtain ku − uh k ≤ Chr+1 . Further from Theorem 2.2, we have the superconvergence result kξ(t)k1 ≤ Chr+1 . Then the use of Sobolev imbedding theorem yields kξ(t)kL∞ ≤ Chr+1 , and hence, k(u − uh )(t)kL∞ ≤ Chr+1 . In this case the degree k for Vh does not influence the L2 and L∞ -estimates of the error u − uh . (iv) Let us compare the above results with Johnson and Thom´ee [17]. In one space dimension (see [17]) Vh consists of discontinuous piecewise linear polynomials
H 1 -MIXED METHOD
717
and the space Wh contains C 0 -piecewise quadratic polynomials. The results obtained there are as Z t 2 kpt (s)k2 ds , k(p − ph )(t)k ≤ Ch kp(t)k2 + 0
and
" k(u − uh )(t)k ≤ Ch
2
Z
1/2
t
kpt (s)k22
kp(t)k3 +
# ds .
0
Of course with C 0 -piecewise linear polynomial spaces for both Vh and Wh , using the proposed method we obtain the following error estimates: " Z t 1/2 # 2 2 kpt (s)k3 ds . k(p − ph )(t)k + k(u − uh )(t)k ≤ Ch kp(t)k3 + 0
Moreover, if c = 0 from (iii), we have for the present case with C 0 -piecewise quadratic polynomial space Wh " Z t 1/2 # 3 2 k(u − uh )(t)k ≤ Ch kp(t)k4 + kpt (s)k4 ds . 0
Therefore, with higher regularity we obtain better estimates. 3. Parabolic equation in several space variables. Let Ω be a bounded domain in Rd , d = 2, 3 with smooth boundary ∂Ω. For simplicity, we assume that the domain Ω is a two-dimensional polygon or a three-dimensional polytope. We now consider the following initial boundary value problem: (3.1)
pt − ∇ · (a(x)∇p) + b(x) · ∇p + c(x)p = f (x, t), (x, t) ∈ Ω × J, p = 0, (x, t) ∈ ∂Ω × J,
and p(x, 0) = p0 , x ∈ Ω, where pt = ∂p/∂t and J = (0, T ] with T < ∞. Introducing αu = ∇p with α = 1/a, (3.1) can be rewritten as a system αu = ∇p,
(3.2) and
pt − ∇·u + β · u + cp = f. Let W = q ∈ (L (Ω)) : ∇ · q ∈ L2 (Ω) with norm kqkH(div,Ω) = (kqk2 + k∇ · qk2 )1/2 . Then weak formulation is now defined to be a pair {p, u} : [0, T ] → H01 × W satisfying (3.3)
(3.4)
2
d
(a∇p, ∇v) = (u, ∇v) , v ∈ H01 ,
and (3.5)
(αut , w) + (∇·u, ∇·w) − (β · u, ∇·w) = (cp, ∇·w) − (f, ∇·w) , w ∈ W.
For our subsequent use, we shall employ the classical Hilbert Sobolev spaces H m (Ω) and shall call this H m with norm k · km . Let (H m )d = Hm denote the corresponding product space with usual product norm.
718
AMIYA K. PANI
3.1. Semidiscrete H 1 -Galerkin mixed finite element procedure. Let Th be a partition of Ω into a finite number of elements called simplexes, i.e, Ω = ∪K∈Th K with h = max {dim(K) : K ∈ Th }. Let Vh and Wh , respectively, be finite dimensional subspaces of H01 and W satisfying the following approximation properties: For k > 0, r > 0 integers inf {kv − vh k + hkv − vh k1 } ≤ Chk+1 kvkk+1 , v ∈ H k+1 ∩ H01 ,
vh ∈Vh
and inf
qh ∈Wh
{kq − qh k + hkq − qh k1 } ≤ Chr+1 kqkr+1 , q ∈ Hr+1 .
Standard examples of such spaces are as follows: Vh = vh ∈ C 0 (Ω) : vh |K ∈ Pk (K), ∀K ∈ Th , vh = 0 on ∂Ω , and Wh = {qh ∈ W : (qh )i |K ∈ Pr (K), i = 1, 2, . . . , d, ∀K ∈ Th } , where Ps (K) is the space of polynnomials of degree ≤ s on K. Other examples of approximating spaces can be found in Raviart and Thomas [27]. Now semidiscrete H 1 -Galerkin finite element procedure for the system is determined as a pair {ph , uh } : [0, T ] → Vh × Wh satisfying (3.6)
(a∇ph , ∇vh ) = (uh , ∇vh ) , vh ∈ Vh ,
and (3.7)
(αuht , wh ) + (∇·uh , ∇·wh ) − (β · uh , ∇·wh ) = (cph , ∇·wh ) − (f, ∇·wh ) , wh ∈ Wh
with the appropriately chosen initial pair {ph (0), uh (0)} to be defined later. Again following Wheeler [31], let us define Ritz projection p˜h ∈ Vh of p as (∇(p − p˜h ), ∇vh ) = 0, ∀vh ∈ Vh . ˜ I ∈ Wh be the standard finite element interpolant of u. Let ρ = u − u ˜I Further, let u and η = p − p˜h . Then for nonnegative integers k and r (3.8)
kηk + hk∇ηk ≤ Chk+1 kpkk+1 ,
and (3.9)
kρk + hkρkH(div,Ω) ≤ Chr+1 kukr+1 .
˜I + u ˜ I − uh = ρ + ξ and p − ph = p − p˜h + p˜h − ph = η + ζ. Let u − uh = u − u For semidiscrete error analysis, we have from (3.4)–(3.5), (3.6)–(3.7) and auxiliary projections the error equations in ξ and ζ as (3.10)
(∇ζ, ∇vh ) = (α(u − uh ), ∇vh ) , vh ∈ Vh ,
H 1 -MIXED METHOD
719
and (3.11)
(αξt , wh ) + (∇ · ξ, ∇ · wh ) = − (αρt , wh ) − (∇ · ρ, ∇ · wh ) + (c(η + ζ), ∇ · wh ) + (β · (ξ + ρ), ∇ · wh ) , wh ∈ Wh .
In the following we shall obtain semidiscrete error estimates for u − uh and p − ph . THEOREM 3.1. With u0 = a∇p0 , assume that ku0 − u0h kH(div,Ω) ≤ Chr ku0 kr+1 . Then there is a constant C independent of h such that " (3.12)
k∇(p − ph )(t)k + k∇ · (u − uh )(t)k ≤ Ch 2
2
2min(r,k)
Z + kp(t)k2k+1 + ku(t)k2r+1 +
ku0 k2r+1 + kp0 k2k+1
t
ku(s)k2r+1 0
#
+ kp(s)k2k+1 + kut (s)k2r+1 + kpt (s)k2k+1 ds and
"
(3.13)
k(p − ph )(t)k + k(u − uh )(t)k ≤ Ch 2
2
Z
2min(r,k+1)
t
ku0 k2r+1 + kp(t)k2k+1
ku(s)k2r+1 + kut (s)k2r + kp(s)k2k+1 ds
+ ku(t)k2r +
# .
0
Proof. Choose vh = ζ in (3.10) to have k∇ζk ≤ C(kρk + kξk).
(3.14)
Further, set wh = ξ in (3.11) and use Cauchy–Schwartz inequality with kζk ≤ k∇ζk to obtain i h d 1/2 2 kα ξk + k∇ · ξk2 ≤ C kρt k2 + kξk2 + kρk2H(div,Ω) + kηk2 + k∇ζk2 . dt Now substituting the estimate of k∇ζk, integrate the resulting inequality with respect to time from 0 to t and apply Young’s inequality to have Z t Z t kρt (s)k2 + kρ(s)k2H(div,Ω) + kη(s)k2 ds] k∇ · ξk2 ds ≤ C[kξ(0)k2 + kξ(t)k2 + 0 0 Z t kξ(s)k2 ds. + 0
Apply Gronwall’s lemma to obtain Z t Z t 2 2 2min(r,k+1) k∇ · ξk ds ≤ Ch kξ(0)k2 + ku(s)k2r+1 kξ(t)k + 0 0 2 2 + kut (s)kr + kp(s)kk+1 ds .
720
AMIYA K. PANI
To obtain an estimate for ∇ζ, apply the above estimate with (3.9) in (3.14) to have " k∇ζk ≤ Ch 2
2min(r,k+1)
Z
t
ku(0)k2r+1 + ku(t)k2r +
ku(s)k2r+1 0
+
kut (s)k2r
+
kp(s)k2k+1
# ds
,
and hence, using (3.8)–(3.9) we obtain the estimates of p−ph and u−uh . To estimate ∇ · ξ, choose wh = ξt in (3.11) to have 1 d d k∇ · ξk2 = −(αρt , ξt ) + (β · ξ + ∇ · ρ + β · ρ + cη, ∇ · ξ) 2 dt dt − (β · ξt + ∇ · ρt + β · ρt + cηt , ∇ · ξ) − (∇(cζ), ξt ) .
kα1/2 ξt k2 +
For the last term we have used the Gauss divergence theorem and ζ ∈ H01 . On integration with respect to time and using the standard Young’s inequality, we obtain Z
t 0
i h kξt (s)k2 ds + k∇ · ξ(t)k2 ≤ C kξ(0)k2H(div,Ω) + kρ(0)k2H(div,Ω) + kη(0)k2 Z t 2 2 2 kρt (s)k2H(div,Ω) + C kρ(t)kH(div,Ω) + kξ(t)k + kη(t)k + 0 Z t + k∇ζ(s)k2 + kηt (s)k2 ds + C k∇ · ξ(s)k2 ds. 0
Using the estimates of k∇ζk, kξk and (3.9) apply Gronwall’s lemma to have Z
t 0
kξt (s)k2 ds + k∇ · ξ(t)k2 ≤ Ch2min(r,k+1) ku0 k2r+1 + kp0 k2k+1 + ku(t)k2r+1 Z t 2 2 2 2 + kp(s)kk+1 + ku(s)kr+1 + kut (s)kr+1 + kpt (s)kk+1 ds . 0
The triangle inequality with (3.8) and (3.9) now completes the rest of the proof. Remark 3.1. (i) Estimates in (3.12) indicate that for k = r, the error estimate k∇ · (u − uh )(t)k is optimal in the stated norm. However, for k 6= r, this estimate distinguishes the better approximation properties of Vh or Wh . For k + 1 = r in (3.13) we can decrease the influence of Vh on the rate of convergence for uh . (ii) We now compare (3.13) with the order of convergence of the mixed finite element method presented in Johnson and Thom´ee [17]. Assume that(Vh , Wh )is a pair of Raviart–Thomas spaces. Then components of Wh contain incomplete polynomials of degree r = k + 1 on each finite element K ∈ Th . From Theorem 2.1 in Johnson and Thom´ee [17], we obtain the following estimates: k(p − ph )(t)k + k(u − uh )(t)k ≤ C(u, p)hr . Hence, the rate of convergence in both the cases coincide, but for the present method LBB condition has been avoided. Moreover, the quasiuniformity assumption is not required in our analysis.
H 1 -MIXED METHOD
721
4. Modified H 1 -Galerkin mixed finite element procedure. The results obtained in one-dimensional situations are quite sharp compared to the results in section 3. The reason for this is that in section 2, we have used elliptic projection (Wheeler’s technique) for the flux (see (2.9)), whereas in case of several spatial vari˜ I of the flux u is used as an intermediate projection. Therefore, ables interpolant u we have only derived optimal estimates for k∇ · (u − uh )k and suboptimal estimates in L2 -norm. But as we have remarked earlier we obtain in section 3 the same order of convergence as in Johnson and Thom´ee [17]. In this section, we shall modify the Galerkin method to obtain sharper estimates as in one-dimensional cases. With αu = ∇p, write (3.1) as pt − ∇·u + β · u + cp = f, (x, t) ∈ Ω × J, ∇ × (αu) = 0, (x, t) ∈ Ω × J, (n ∧ αu) = 0, (x, t) ∈ ∂Ω × J, p(0) = p0 , x ∈ Ω, where n is the outward normal and ∧ denotes the exterior product. For the weak formulation, let W = {w ∈ (H 1 )d : n ∧ αw = 0 on ∂Ω, d = 2, 3}. Using the Gauss divergence theorem we now have a pair {p, u} : [1, T ] → H01 × W such that (4.1)
(∇p, ∇v) = (αu, ∇v), v ∈ H01 ,
(4.2)
(αut , w) + A(u, w) − (β · u, ∇ · w) = (−f + cp, ∇ · w), w ∈ H1 ,
where A(φ, w) = (∇ · φ, ∇ · w) + (∇ × αφ, ∇ × αw). We note that A(·, ·) satisfies the following coercivity condition: A(φ, φ) ≥ µ0 kφk21 , for some positive constant µ0 ; see [19]. For the Galerkin procedure, we first consider a finite dimensional space Vh as in section 3 and then define ¯ d : (wh )i |K ∈ Pr (K), i = 1, . . . , d, ∀K ∈ Th Wh = {wh ∈ C(Ω) (n ∧ αwh ) = 0, at the nodes on ∂Ω}. Since (n ∧ αwh ) = 0 only at the boundary nodes, the finite element space Wh is not a subspace of W and hence, we have nonconforming methods. Note that the above finite dimensional spaces satisfy the same approximation properties as in section 3; see [19] for r = 1. A modified H 1 -Galerkin mixed finite element approximation is determined as a pair {ph , uh } : [0, T ] → Vh × Wh such that (4.3) (4.4)
(∇ph , ∇vh ) = (αuh , ∇vh ), v ∈ Vh , (αuht , wh ) + A(uh , wh ) − (β · uh , ∇ · wh ) = (−f + cph , ∇ · wh ), wh ∈ Wh ,
with uh (0) to be defined later.
722
AMIYA K. PANI
˜ h } ∈ Vh × Wh as Now define auxiliary projections {˜ ph , u (4.5)
(∇(p − p˜h ), ∇vh ) = 0, vh ∈ Vh ,
(4.6)
˜ h , wh ) = 0, wh ∈ Wh . A(u − u
˜ h = ρ. Following Neittaanm¨ aki and Saranen [18]–[20], we Let p − p˜h = η and u − u have kρkj + kρt kj ≤ Chr+1−j [kukr+1 + kut kr+1 ], j = 0, 1. Note that the related difficulties with the nonconforming finite element method will mainly show up in the error estimates for ρ. When the paper was under revision the author came across a paper by A.I. Pehlivanov and G. F. Carey [26], in which a similar analysis was discussed in the context of elliptic equation. With an appropriate modification of the analysis of Pehlivanov and Carey, it is possible to obtain the above error estimates for ρ and ρt . For semidiscrete error estimates, we now split the errors p − ph = (p − p˜h ) + (˜ ph − ˜ h ) + (˜ uh − uh ) = ρ + ξ. Below, we shall state and ph ) = η + ζ and u − uh = (u − u prove our main theorem in this section. ˜ h (0) with u(0) = a∇p0 . Then there exists THEOREM 4.1. Assume that uh (0) = u a positive constant C independent of h such that k(p − ph )(t)k + k(u − uh )(t)k + hk(p − p)(t)k1 ≤ Chmin(k+1,r+1) [kpkL∞ (H k+1 ) + kukL∞ (Hr+1 ) + kut kL2 (Hr+1 ) ]. Further, we have k(u − uh )(t)k1 ≤ C hmin(k+1,r) [kpkL∞ (H k+1 ) + kukL∞ (Hr ) + kpt kL2 (H k+1 ) + kut kL2 (Hr ) ]. Proof. From (4.1)–(4.6) we have (∇ζ, ∇vh ) = (α(ρ + ξ), ∇vh ), vh ∈ Vh ,
(4.7) and (4.8)
(αξt , wh ) + A(ξ, wh ) = (β · (ρ + ξ), ∇ · wh ) − (αρt , wh ) + (c(ζ + η), ∇ · wh ), wh ∈ Wh .
Choose vh = ζ in (4.7) and wh = ξ in (4.8) to obtain k∇ζk ≤ C(kρk + kξk),
(4.9) and
1 1 d kα 2 ξk2 + (µ0 − C()) kξk21 ≤ C[kρk2 + kρt k2 + kξk2 + kηk2 + k∇ζk2 ]. 2 dt Here we have used Cauchy–Schwartz and Young’s inequalities with > 0. Now choose appropriately so that (µ0 −C()) > 0. Integrating the resulting equation with respect to time, use (4.9) and then apply Gronwall’s lemma to have Z t Z t kξ(s)k21 ds ≤ C [kρ(s)k2 + kρt (s)k2 + kη(s)k2 ] ds. kξ(t)k2 +
0
0
The triangle inequality with estimates ρ and η now completes the proof of the first estimate.
H 1 -MIXED METHOD
723
For a superconvergence result in kξk1 , choose wh = ξt in (4.8) and obtain 1 d d A(ξ, ξ) = −(αρt , ξt ) + (β · (ρ + ξ) + cη, ∇ · ξ) 2 dt dt − (β · (ρt + ξt ) + cηt , ∇ · ξ) + (∇(cζ), ξt ).
1
kα 2 ξt (t)k2 +
On integrating with respect to time, apply Young’s inequality and estimate (4.9) to have Z t Z t 2 2 kξ(t)k1 + kξt (s)k ds ≤ C (kρ(s)k2 + kρt (s)k2 + kηt (s)k2 + kξ(s)k2 ) ds 0 0 Z t + kη(t)k2 + kρ(t)k2 + kξ(t)k2 + C kξ(s)k21 ds. 0
An application of Gronwall’s lemma yields kξ(t)k21 ≤ C hmin(k+1,r+1) ku(t)k2r+1 + kp(t)k2k+1 Z +
t
(kpt (s)k2k+1 + ku(s)k2r+1 + kut (s)k2r+1 ) ds .
0
It is an easy excercise to reduce the regularity requirement on u for the estimate kξk1 , if only the optimal estimate of u − uh in H 1 -norm is required. The triangle inequality now completes the rest of the proof. Remark 4.1. (i) With k = r we have kp−ph kL∞ (L2 ) +ku−uh kL∞ ((L2 )d ) = O(hr+1 ). (ii) Compared to [17] the present analysis yields better results for L2 -estimate of u − uh with a higher regularity assumption on the exact solution. When c = 0 in (3.1), we have again ku − uh k = O(hr+1 ) even if k 6= r with k < r. ˜ h (0), then using Sobolev’s (iii) When d = 2 that is Ω ⊂ R2 and uh (0) = u imbedding theorem and superconvergence estimates for the ∇ζ and ∇ξ, we have 1 1 2 [kpkL∞ (H k+1 ) kζ(t)kL∞ ≤ Chmin(k+1,r+1) log h + kukL∞ (Hr+1 ) + kut kL2 (Hr+1 ) ]. Similarly, 1/2 1 log ≤ Ch [kpkL∞ (H k+1 ) h + kukL∞ (Hr+1 ) + kpt kL2 (H k+1 ) + kut kL2 (Hr+1 ) ].
kξ(t)kL∞
min(k+1,r+1)
From Thom´ee [29], it follows that the error η in the Ritz projection satisfies (with k = 1 say) 1 hk+1 kpkL∞ (W k+1,∞ ) . kη(t)kL∞ ≤ C log h In this case, we have to assume the quasiuniformity condition on the finite element space for the maximum norm estimate. Therefore, for d = 2, we obtain the following maximum norm estimates, 1 k(p − ph )(t)kL∞ ≤ C log hmin(k+1,r+1) [kpkL∞ (H k+1 ) + kukL∞ (W r+1,∞ ) h + kut kL2 (Hr+1 ) ].
724
AMIYA K. PANI
Finally, it is possible to obtain a maximum norm estimate for the flux, provided the maximum norm estimate for ρ is available. 5. Fully discrete schemes. In this section, we shall briefly describe the backward Euler method for approximating {p, u} in (4.1)–(4.2) and discuss the related error estimates. Note that the error analysis for higher-order schemes like Crank– Nicolson and two-step backward difference methods will follow similarly with appropriate changes. Therefore, we shall not pursue these methods further in this paper. For the backward Euler procedure, let 0 = t0 < t1 < · · · < tM = T be a given partition of the time interval [0,T] with step length ∆t = T /M , for some ¯ n= positive integer M . For a smooth function φ on [0,T], define φn = φ(tn ) and ∂φ n n−1 n n /∆t. Let P and U , respectively, be the approximations of p and u at φ −φ t = tn which we shall define through the following scheme. Given {P n−1 , Un−1 } in Vh × Wh , we now determine a pair {P n , Un } in Vh × Wh satisfying (5.1)
(∇P n , ∇vh ) = (αUn , ∇vh ), v ∈ Vh ,
(5.2)
¯ n , wh ) + A(Un , wh ) − (β · Un , ∇ · wh ) (α∂U = (−f n + cP n , ∇ · wh ), wh ∈ Wh ,
with P 0 = p0,h to be defined later. For fully discrete error estimates, we now split the errors p(tn ) − P n = (p(tn ) − ˜ h (tn ))+(˜ p˜h (tn ))+(˜ ph (tn )−P n ) = η n +ζ n and u(tn )−Un = (u(tn )− u uh (tn )−Un ) = n n n n ρ + ξ . Since the estimates of η and ρ can be easily found out it is enough to estimate ζ n and ξ n . Note that the error equations in ζ n and ξ n may be written as (∇ζ n , ∇vh ) = (α(ρn + ξ n ), ∇vh ), vh ∈ Vh ,
(5.3) and (5.4)
¯ n , wh ) + A(ξ n , wh ) = (β · (ρn + ξ n ), ∇ · wh ) − (α∂ρ ¯ n + ασ n , wh ) (α∂ξ + (c(ζ n + η n ), ∇ · wh ), wh ∈ Wh ,
¯ n ). where σ n = ut (tn ) − ∂u(t ˜ h (0) with u(0) = a∇p0 . Then there exists THEOREM 5.1. Assume that U0 = u a positive constant C independent of h and ∆t such that for 0 < ∆t ≤ ∆t0 and J = 0, 1, . . . , M kp(tJ ) − P J k + ku(tJ ) − UJ k + hkp(tJ ) − P J k1 ≤ Chmin(k+1,r+1) [kpkL∞ (H k+1 ) + kukL∞ (Hr+1 ) + kut kL2 (Hr+1 ) ] + C∆tkutt kL2 ((L2 )d ) . Further, we have ku(tJ ) − UJ k1 ≤ Chmin(k,r) [kpkL∞ (H k+1 ) + kukL∞ (Hr+1 ) + kpt kL2 (H k+1 ) + kut kL2 (Hr+1 ) ] + C∆tkutt kL2 ((L2 )d ) . Proof. Choose vh = ζ n in (5.3) to have for n = 0, 1, . . . , M (5.5)
k∇ζ n k ≤ C(kρn k + kξ n k).
Set wh = ξ in (5.4) and use Cauchy–Schwarz as well as Young’s inequality to obtain 1¯ 1 n 2 ¯ n k2 + kξ n k2 + kη n k2 + kσ n k2 ]. ∂kα 2 ξ k + (µ0 − C()) kξ n k21 ≤ C[kρn k2 + k∂ρ 2
H 1 -MIXED METHOD
Here we have also used (5.5). Note that Z kσ n k2 ≤ ∆t
tn
725
kutt (s)k2 ds,
tn−1
and ¯ n k2 ≤ 1 k∂ρ ∆t
Z
tn
kρt (s)k2 ds.
tn−1
On substitution and summing from n = 1, . . . , J with = µ0 /4, the resulting equation becomes " J J X X J 2 n 2 kξ k1 ≤ C ∆t (kρn k2 + kη n k2 ) (1 − C∆t)kξ k + µ0 n=0
Z
tJ
n=0
Z kρt (s)k ds + (∆t) 2
+ 0
2
tJ
# kutt (s)k ds 2
+ C∆t
0
J−1 X
kξ n k2 .
n=0
Choose ∆t0 in such a way that for 0 < ∆t ≤ ∆t0 , (1 − C∆t) > 0. Then an application of Gronwall’s lemma with triangle inequality completes the L2 -estimate. Further use (5.5) at n = J and the triangle inequality to complete the error estimate for kp(tn ) − P n k1 . ¯ n in (5.4) to obtain For the superconvergence result in kξ n k1 , set wh = ∂ξ ¯ n ) = −(α∂ρ ¯ n + ασ n , ∂ξ ¯ n k2 + A(ξ n , ∂ξ ¯ n ) + ∂(β ¯ · (ρn + ξ n ) + cη n , ∇ · ξ n ) kα1/2 ∂ξ ¯ n ) + c∂η ¯ n , ∇ · ξ n−1 ) − (∇ · (cζ n ), ∂ξ ¯ n ). ¯ n + ∂ξ − (β · (∂ρ Note that ¯ n) = A(ξ n , ∂ξ
1¯ ∆t ¯ n ¯ n ∂A(ξ n , ξ n ) + A(∂ξ , ∂ξ ). 2 2
Now summing from n = 1 to J, we use standard arguments to obtain ∆t
J X
kξ n k2 + 2(µ0 − C∆t)kξ J k21 ≤ Ch2min(k+1,r+1) [kpk2L∞ (H k+1 )
n=1
Z +
kuk2L∞ (Hr+1 )
+
tJ
(kpt (s)k2k+1 + kut (s)k2r+1 ) ds]
0
Z + (∆t)2 0
tJ
kutt (s)k2 ds + C∆t
J−1 X
kξ n k21 .
n=0
Again apply Gronwall’s lemma to obtain a superconvergence result for kξ J k1 . A use of triangle inequality now completes the rest of the proof. 6. Extensions and concluding remarks. In this paper a priori error estimates are presented for the H 1 -Galerkin mixed finite element method. The Galerkin approximation in section 3 has the same rate of convergence as in the case of classical mixed finite element method, but without LBB consistency condition and also without quasiuniform assumption on the finite element mesh. Since the estimate for ku − uh k is not optimal in two- and three-space dimensions, a modification is proposed and
726
AMIYA K. PANI
analyzed in section 4 to derive a sharper estimate in L2 -norm. This confirms the findings in the one-dimensional case. Compared to classical mixed method, the present method allows us to use two different finite element spaces for approximating p and its flux u. In particular, use of piecewise linear polynomial spaces yields O(h2 ) order of convergence in both p − ph and u − uh in L2 -norm. Although higher regularity on the solution is assumed, better convergence results are proved specifically for the flux. In the future, analysis of this method will be extended to second-order wave equations and miscible displacement problems in porous media. Acknowledgments. The author is grateful to the referee for his/her valuable suggestions and remarks. REFERENCES ´ , AND M. FORTIN, Mixed finite elements for second [1] F. BREZZI, J. DOUGLAS, JR., R. DURAN order elliptic problems in three variables, Numer. Math., 51 (1987), pp. 237–250. [2] F. BREZZI, J. DOUGLAS, JR., AND L. D. MARINI, Two families of mixed finite elements for second order elliptic problems, Numer. Math., 47 (1985), pp. 217–235. [3] S. L. CAMBELL, K. E. BRENAN, AND L. R. PETZOLD, Numerical Solution of Initial Value Problems in Differential-Algebraic Equations, American Elsevier Science, New York, 1989. [4] G. F. CAREY AND Y. SHEN, Convergence studies of least-square finite elements for first order systems, Comm. Appl. Numer. Methods, 5 (1989), pp. 427–434. [5] C. L. CHANG, A least-squares finite element method for the Helmholtz equation, Comput. Methods Appl. Mech. Engrg., 83 (1990), pp. 1–7. [6] T. F. CHEN, On least-square approximations to compressible flow problems, Numer. Methods Partial Differential Equations, 2 (1986), pp. 207–228. [7] S. H. CHOU AND Q. LI, Error estimates for mixed finite element methods for nonlinear parabolic problems, Numer. Methods Partial Differential Equations, 8 (1992), pp. 395–404. [8] L. C. COWSAR, T. F. DUPONT, AND M. F. WHEELER, A priori estimates for mixed finite element methods for the wave equation, Comput. Methods Appl. Mech. Engrg., 82 (1990), pp. 205–222. [9] J. DOUGLAS, JR., T. F. DUPONT, AND M. F. WHEELER, H 1 -Galerkin methods for the Laplace and heat equations, Mathematical Aspect of Finite Elements in Partial Differential Equations, C. de Boor, ed., Academic Press, New York, 1975, pp. 383–415. [10] J. DOUGLAS, JR., R. E. EWING, AND M. F. WHEELER, Approximation of the pressure by a mixed method in the simulation of miscible displacement, RAIRO Anal. Numer., 17 (1983), pp. 17–33. [11] J. DOUGLAS, JR., AND J. E. ROBERTS, Global estimates for mixed methods for the second order elliptic equations, Math. Comp., 44 (1985), pp. 39–52. ´ , On the approximation of miscible displacement in porous media by a method of [12] R. G. DURAN characteristics combined with a mixed method, SIAM J. Numer. Anal., 25 (1988), pp. 989– 1001. [13] R. E. EWING, T. F. RUSSEL, AND M. F. WHEELER, Convergence analysis of an approximation of miscible displacement in porous media by mixed finite elements and a modified method of characteristics, Comput. Methods Appl. Mech. Engrg., 47 (1984), pp. 73–92. [14] G. J. FIX, M. D. GUNZBURGER, AND R. A. NICOLAIDES, On the mixed finite element methods for first order elliptic systems, Numer. Math., 37 (1981), pp. 29–48. [15] T. GEVECI, On the application of mixed finite element methods to the wave equation, Math. Model. Numer. Anal, 22 (1988), pp. 243–250. ¨ [16] J. HASLINGER AND P. NEITTAANMAKI , On different finite element methods for approximating the gradient of the solution to the Helmholtz equation, Comput. Methods Appl. Mech. Engrg., 42 (1984), pp. 131–148. ´ , Error estimates for some mixed finite element methods for [17] C. JOHNSON AND V. THOMEE parabolic type problems, RAIRO Numer. Anal., 15 (1981), pp. 41–78. ¨ AND J. SARANEN, On finite element approximation of the gradient for [18] P. NEITTAANMAKI solution of Poisson equations, Numer. Math., 37 (1981), pp. 333–337. ¨ [19] P. NEITTAANMAKI AND J. SARANEN, On the finite element approximation of vector fields by curl and divergence, Math. Methods Appl. Sci., 3 (1981), pp. 328–335. ¨ AND J. SARANEN, A mixed finite element method for the heat flow problem, [20] P. NEITTAANMAKI BIT, 21 (1981), pp. 342–346.
H 1 -MIXED METHOD
727
[21] A. K. PANI, Mixed Finite Element Method for a Stefan Problem, Tech. report, Department of Mathematics, IIT, Bombay, 1988. [22] A. K. PANI AND P. C. DAS, An H 1 -Galerkin method for quasilinear parabolic partial differential equations, Methods of Functional Analysis in Approximation Theory, C. A. Micchelli, D. V. Pai, and B. V. Limaye, eds., ISNM 76, Birkh¨ auser-Verlag, Basel, 1986, pp. 357–370. [23] A. K. PANI AND R. S. ANDERSSEN, Finite element methods for identification of parameters in parabolic problems, Inverse Problems in Partial Differential Equations, Proc. Centre for Mathematics and Its Applications, Vol. 31, A. K. Pani and R. S. Anderssen, eds., Australian National University, Canberra, 1992, pp. 208–221. [24] A. I. PEHLIVANOV, G. F. CAREY, R. D. LAZAROV, AND Y. SHEN, Convergence analysis of least-square mixed finite elements, Computing, 51 (1993), pp. 111–123. [25] A. I. PEHLIVANOV, G. F. CAREY, AND R. D. LAZAROV, Least-square mixed finite elements for second order elliptic problems, SIAM J. Numer. Anal., 31 (1994), pp. 1368–1377. [26] A. I. PEHLIVANOV AND G. F. CAREY, Error estimates for least-squares mixed finite elements, M 2 AN, Math. Model. Numer. Anal., 28 (1994), pp. 499–516. [27] P. A. RAVIART AND J. M. THOMAS, A mixed finite element method for second order elliptic problems, Lecture Notes in Mathematics 606, Springer-Verlag, New York, 1977, pp. 293– 315. [28] T. F. RUSSEL AND M. F. WHEELER, Finite element and finite difference methods for continuous flows in porous media, The Mathematics of Reservoir Simulation, R.E. Ewing, ed., SIAM, Phildelphia, 1983, pp. 35–106. ´ , Galerkin finite element methods for parabolic problems, Lecture Notes in Mathe[29] V. THOMEE matics 1054, Springer-Verlag, New York, 1984. ´ AND L. B. WAHLBIN, On Galerkin methods in semilinear parabolic problems, [30] V. THOMEE SIAM J. Numer. Anal., 12 (1989), pp. 376–389. [31] M. F. WHEELER, A priori L2 -error estimates for Galerkin approximations to parabolic differential equations, SIAM J. Numer. Anal., 10 (1973), pp. 723–749.
