ERROR ANALYSIS OF THE GENERALIZED MAC SCHEME
YIN-LIANG HUANG ∗, JIAN-GUO LIU †, AND WEI-CHENG WANG ‡
Abstract.
We present a rigorous convergence analysis for the generalized MAC (gMAC) scheme on curvilinear domains proposed earlier by the authors [HLW]. The error estimate for the velocity field is established by energy estimate utilizing the stream function and discrete identities associated with the spatially compatible discretization. The spatially compatible discretization also induces subtle stabilizing effect that renders the scheme uniform LBB bound even though gMAC is staggered and supported the same way as the Q1− P0 finite element method, which is well known to be unstable under divergence constraint. As a result, full second order error estimate is achieved for both velocity and pressure with minimal regularity requirement.
1. Introduction. A generalized MAC scheme (gMAC) for the Navier-Stokes equation in rotational form
(1.1)
ut+ ω × u + ∇p = ν∇2u+ f on Ω
∇ · u = 0 on Ω
u = 0 on Γ
on curvilinear domains was introduced in [HLW]. With partially staggered grids (velocity components collocated on cell centers, pressure placed on grid points) and centered difference in a locally ’skewed’ coordinate, the discretization preserves crucial identities such as
(1.2) curlh◦ gradh≡ 0, divh◦ curlh≡ 0
and their converse in the discrete setting. The resulting scheme is simple and efficient with full second order accuracy on curvilinear domains. A key ingredient of the scheme is the proper treatment at the boundary, which not only enforces the pressure as discrete Lagrangian multiplier without introduction artificial bound-ary conditions, but also leads to an exact Hodge decomposition for the velocity field, which plays a key role in both stability and efficiency of the scheme.
In this paper, we will show that the spatial compatibility (1.2) (recast as Lemma 2.1 in section 2) also leads to a simple error estimate for the velocity field. Optimal error analysis for the classical MAC scheme was first obtained in [HW], and in [W1] for MAC-like schemes on Cartesian grids. The proof in [HW, W1] is based on high order Strang’s expansion. Here we present an alternative approach that relies mainly on the special structure of the spatial discretization, making use of both the stream function and the discrete differential identities in Lemma 2.1. As a result, we obtain optimal O(h2) error estimate for the velocity
field provided the exact solution satisfies ue∈ L2(0, T ; C4( ¯Ω)) ∩ H1(0, T ; C2( ¯Ω)) and pe∈ L2(0, T ; C3( ¯Ω))
(Theorem 1, section 3.1). This may be the minimal regularity requirement in finite difference setting. On the other hand, the error analysis for the pressure is much more subtle due to lack of evolutionary equation. Our approach is based on establishing the Ladyzhenskaya-Babuˇska-Brezzi condition (also known as the LBB, div-stability or inf-sup condition) for the generalized MAC scheme. The LBB condition provides direct access to pressure error estimate for the dynamic problem (1.1), and is essential to the solvability and uniform estimate for the static Stokes problem. It is worth noting that gMAC is staggered the same way as the Q1−P0finite element method, which is known to violate the uniform LBB condition. The main difference
between the two schemes is the discretization of the viscous term. The spatially compatible discretization associated with gMAC induces subtle stabilizing effect and is key to the uniform LBB bound. As a result, we
∗Department of Mathematics, National Taiwan University, Taipei, 106, Taiwan. (yin1108@gmail.com)
†Department of Physics and Department of Mathematics, Duke University, Durham, NC 27708, USA (jian-guo.liu@duke.edu)
‡Department of Mathematics, National Tsing Hua University, HsinChu, 300, Taiwan. (wangwc@math.nthu.edu.tw) 1
ξ2= constant ξ2= constant ξ2 = constant ξ2 = constant ξ1= constant ξ1= constant ξ1= constant ξ1= constant
Fig. 2.1. The computational domain (left) and physical domain (right)
obtain second order pressure estimate with ue∈ L2(0, T ; C4( ¯Ω)) ∩ H1(0, T ; C3( ¯Ω)) and pe∈ L2(0, T ; C3( ¯Ω))
(Theorem 3, section 3.2).
The rest of the paper is organized as follows. In section 2 we review the generalized MAC scheme, the boundary treatment and relevant identities associated with the spatial discretizations. In section 3.1, we proceed with the basic energy estimate and give a second order error estimate for the velocity and a first order error estimate for the pressure gradient. The full second order pressure estimate utilizing the LBB condition is given in section 3.2. Finally, the technical proof of the LBB condition is given in the Appendix. 2. Generalized MAC Scheme in Curvilinear Coordinate. We briefly summarize the spatial dis-cretization of the generalized MAC scheme in the 2D setting. The details can be found in [HLW]. Denote by x = (x, y) the position vector in the physical domain and (ξ1, ξ2) the coordinate in the computational domain with mesh size ∆ξ1= h1and ∆ξ2= h2. Instead of the default coordinate (ξ1, ξ2), a local coordinate
in the skewed direction (Figure 2.1)
(2.1) ξ 1 ,h2ξ1+ h1ξ2 ph2 1+ h22 , ξ 2 , −h2ξ1+ h1ξ2 ph2 1+ h22 .
is selected to realize the discretization of (1.1).
With the skewed local coordinate, the differential operators can be determined following standard pro-cedure. The superscript ‘
’ is used to denote quantities computed in the skewed variables ξ
α
. The basis vectors and metric tensors with respect to the skewed coordinate (ξ
1
,ξ
2
) are thus denoted by
(2.2) e1= ∂x ∂ξ 1, e2= ∂x ∂ξ 2, g µν =eµ·eν, (2.3) e 1 = ∇ξ 1 , e 2 = ∇ξ 2 , g µν =e µ ·e ν . The Jacobian of the transformation between x and
ξis given by (2.4) pg = det ∂x ∂ξ =qdet(g µν).
The covariant and contra-variant components of a vector field u are defined through
(2.5) u=u 1 e1+u 2 e2=u1e 1 +u2e 2 . 2
The transformation between the covariant and contra-variant components for the metric tensor and for a vector field are given by
(2.6) 2 X γ=1 gµγg γν= δµν and (2.7) u µ = 2 X γ=1 gµγuγ, uν = 2 X γ=1 gγνu γ , µ, ν = 1, 2.
The discretization of (1.1) is based on centered difference approximation of the intrinsic differential operators given below in (2.19)-(2.23). The metric tensors involved there can be calculated either analytically, given explicit form of the mapping (ξ1, ξ2) 7→ (x, y), or numerically from centered difference approximation
of (2.2) for the covariant component, and then from (2.4) and (2.6) for the numerical Jacobian and contra-variant components.
Here for simplicity of presentation, we assume the physical domain Ω is diffeomorphic to a ring, so that a single coordinate chart (ξ1, ξ2) ∈ (0, 1) × S1 is sufficient to represent the computational domain. The generalized MAC scheme can be applied to a generic domain by decomposing it into non-overlapping quadrilateral sub-domains. The details of the discretizations on coordinate interfaces and junctions, as well as the 3D case can be found in [HLW].
We further assume equal spacing in ξ1 and ξ2:
(2.8) h1= h2= h = h √ 2 = 1 N, whereh = ∆ξ 1 = ∆ξ 2 = √2h1h2 h2 1+h22
is the natural grid spacing in the skewed variables. When h1= h2, we also
have √g =√g. The relevant domains for spatial discretizations are summarized as follows:
Ωc, x(ξ1i−1 2, ξ 2 j−1 2) | 1 ≤ i ≤ N; 1 ≤ j ≤ N (2.9) ˚ Ωg, x(ξ1i, ξj2) | 1 ≤ i ≤ N − 1; 1 ≤ j ≤ N (2.10) ¯ Ωg, x(ξ1i, ξj2) | 0 ≤ i ≤ N; 1 ≤ j ≤ N (2.11) Γg, ¯Ωg\ ˚Ωg (2.12) ¯ Ωge, x(ξ 1 i, ξj2) ∈ ¯Ωg| i + j is even (2.13) ¯ Ωgo , x(ξ 1 i, ξj2) ∈ ¯Ωg| i + j is odd (2.14)
By assumption, all scalar and vector fields under consideration, including L2( ¯Ωg, R),ω : ¯Ωg→ R (2.15) L2(Ω c, R2), u : Ωc→ R2 (2.16) and L2( ¯Ωg, R)/R2,p ∈ L2( ¯Ωg, R) X ¯ Ωge ′(p ghp)i,j= 0 = X ¯ Ωgo ′(p ghp)i,j , (2.17) L2 c( ¯Ωg, R),ψ ∈ L2( ¯Ωg, R)
ψi,j= constant on i = 0 and i = N, respectively
(2.18)
are all periodic in ξ2. Here in (2.17), the primed sums denote summing with half weight on boundary grids Γg∩ ¯Ωge and Γg∩ ¯Ωgo respectively. For ω, p ∈ L2( ¯Ω g, R) and u =u 1 e1+u 2 e2=u1e 1 +u2e 2 ∈ L2(Ω c, R2), we define (2.19) ∇h: L2( ¯Ωg, R) 7→ L2(Ωc, R2), ∇hp, (D 1p)e 1 + (D 2p)e 2 (2.20) ∇ ⊥ h : L2( ¯Ωg, R) 7→ L2(Ωc, R2), ∇ ⊥ hω,− D2ω p gh e1+ D1ω p gh e2 (2.21) ∇′h· : L2(Ωc, R2) 7→ L2( ¯Ωg, R), ∇′h· u = 1 p gh D′1( p ghu 1) +D ′ 2( p ghu 2) (2.22) ∇⊥′h · : L2(Ωc, R2) 7→ L2( ¯Ωg, R), ∇⊥′h · u = 1 p gh (D ′ 1u2− D′2u1) and (2.23) △′h: L2( ¯Ωg, R) 7→ L2( ¯Ωg, R), △′hp = 1 p gh 2 X µ,ν=1 D′µ( p ghg µν h Dνp).
The primed operators in (2.21)-(2.23) denote the ‘reduced’ operators following Anderson [AN]. The reduction only takes place near boundary, where all quantities involving the metric tensors located outside the computational domain are set to zero, followed by proper normalization. For example, at interior grids 1 < i < N , (2.23) gives the full Laplacian
( △′hp)i,j= ( △hp)i,j= 1 p ghi,j (q 11 h D 1p +q 12 h D 2p)i+1 2,j+ 1 2 h + (q 21 h D 1p +q 22 h D 2p)i−1 2,j+ 1 2 h −( q11hD 1p +q 12 hD 2p)i−1 2,j−12 h − (q 21 h D 1p +q 22 h D 2p)i+1 2,j−12 h ! (2.24) whereq αβ h = p ghg αβ
h . At a boundary grid, say i = 0, the discrete Laplacian reduces to
(2.25) ( △′hp)0,j = 2 p gh0,j (q 11 hD 1p +q 12 h D 2p)1 2,j+ 1 2 h − (q 21 hD 1p +q 22 hD 2p)1 2,j− 1 2 h ! . The detailed formula of (2.19)-(2.23) can be found in [HLW]. It can be shown that
(2.26) ker(
△′h) = ker(
∇h) = span{1Ω¯
ge, 1Ω¯go}
In case N is odd, ¯Ωge and ¯Ωgo coincide due to periodicity in ξ
2. To be definite, we assume without loss of
generality that N is even and therefore dim ker(
△′h) = dim ker(
∇h) = 2.
The significance of the reduced operator can be seen from the role it plays in the adjointness with respect to the natural inner products:
h u , v iΩc= h 2 N X i=1 N X j=1 (u · v)√gh i−1 2,j−12 = h 2 N X i=1 N X j=1 (u1v 1 +u2v 2 )√gh i−1 2,j−12 = h2 N X i=1 N X j=1 (u 1 v1+u 2 v2)√gh i−1 2,j− 1 2 , u, v ∈ L2(Ωc, R2), (2.27) 4
(2.28) h a , b iΩ¯g = h 2 M X i=0 ′ N X j=1 a b√gh i,j, a, b ∈ L 2( ¯Ω g, R).
More precisely, we have the following Lemma from [HLW] which plays an essential role in the error analysis to be presented below: Lemma 2.1. Let u ∈ L2(Ω c, R2) and a ∈ L2( ¯Ωg, R), we have 1. (2.29) h u , ∇ha iΩc = −h ∇′h· u , a iΩ¯g 2. (2.30) h u , ∇ ⊥ ha iΩc= −h ∇⊥′h · u , a iΩ¯g 3. (2.31) ∇′h· ∇ha = ∇⊥′h · ∇⊥ha = △′ha on ¯Ωg; 4. If a ∈ L2( ¯Ω g, R), then (2.32) ∇′h· ∇⊥ha = ∇⊥′h · ∇ha = 0 on ˚Ωg. In addition, if a ∈ L2 c( ¯Ωg, R), then (2.33) ∇ ′ h·∇ ⊥ ha =∇ ⊥′ h ·∇ ha = 0 on ¯Ωg.
In addition to Lemma 2.1, the reduced operators also provide a way of incorporating the slip, no-penetration conditions at the physical boundary. The resulting scheme for (1.1) is given by
(2.34) ut+ ¯ωu⊥+ ∇hp = ν ∇⊥hω + f on Ωc ω = ∇⊥′h · u on ¯Ωg ∇′h· u = 0 on ¯Ωg
The reduced divergence operator in the third equation of (2.34) has implicitly incorporated the no-penetration condition u · n = 0 in a natural way. On the other hand, the reduced curl operator in the second equation of (2.34) has implicitly incorporated the no-slip condition u × n = 0 on Γg. This can be interpreted as an
implicit form of local vorticity boundary condition.
3. Error Estimate for the Generalized MAC Scheme. We now proceed with our main result, the second order error estimate for the generalized MAC scheme. Rigorous 2nd order error estimate for the classical MAC scheme and some variants were first obtained in [HW] and [W1]. The method used in [HW, W1] is a combination of energy estimate and high order Strang’s expansion. Here we propose an alternative proof and apply it to our scheme. In addition to extending the analysis to curvilinear domains, our method differs from [HW, W1] in several aspects. The first new component in our analysis is to utilize the stream function and combine it with the discrete identity (2.33) in our analysis. As a result, the regularity requirement on the exact solution becomes transparent and less stringent. Secondly, our pressure error analysis is established via uniform inf-sup (LBB) estimate, which is of independent importance and has potential applications in other areas such as computational elasticity and computational electromagnetics. The verification of the inf-sup condition is quite technical and is left in the Appendix.
3.1. Basic Error Estimate. Our first main result, 2nd order error estimate for the velocity field, is obtained from basic energy estimate.
Theorem 1. Assume the mapping x : (ξ1, ξ2) 7→ (x, y) is a C4 bijection from [0, 1] × S1to ¯Ω ⊂ R2. Let
ue∈ L2(0, T ; C4( ¯Ω)) ∩ H1(0, T ; C2( ¯Ω)), pe∈ L2(0, T ; C3( ¯Ω)) be an exact solution of (1.1), and uh, ωh, ph
the numerical solution of (2.34) with initial velocity uh
0 satisfying ∇′h· uh0 = 0. Then we have max [0,T ]ku h − uek2Ωc+ ν Z T 0 kω h − ωek2Ω¯g ≤ K1(kuh0− ue0k2Ωc+ h4), (3.1) Z T 0 k ∇h(ph− pe)k2Ωc ≤ K2 h−2kuh0− ue0k2Ωc+ h2. (3.2)
where K1, K2 are constants depending on T , ν, kuekL2(0,T ;C4( ¯Ω)), kuekH1(0,T ;C2( ¯Ω)), kpekL2(0,T ;C3( ¯Ω)), but
not on h.
A crucial part of Theorem 1 is to utilize the exact stream function ψe to construct a divergence free
approximate solution of the form
(3.3) ua(t, ·) =
∇⊥h(ψe+h
2
ϕ) ∈ L2(Ωc, R2)
such that both ua − ue and uh− ua are O(h2). The latter requires in addition that
∇⊥′h · ua − ωe =
△′h(ψe+h
2
ϕ) − ωe= O(h2). Here ϕ is a correction to be determined from the following Lemma.
Lemma 3.1. If the mapping x : (ξ1, ξ2) 7→ (x, y) is a C4 bijection from [0, 1] × S1 to ¯Ω ⊂ R2 and
ψ ∈ C4( ¯Ω), then △′hψ x(0, ξ2) = 2√2 h g 11∂ 1ψ + g12∂2ψ + △ψ x(0, ξ2) + P L(ψ)(ξ2)h + Qh L(ψ)(ξ2)h 2 , △′hψ x(1, ξ2) = −2√2 h g 11∂ 1ψ + g12∂2ψ + △ψ x(1, ξ2) + PR(ψ)(ξ2)h + Q h R(ψ)(ξ2)h 2 , (3.4) where (3.5) kPL,R(ψ)kC1(S1)+ kQhL,R(ψ)kC0(S1)≤ CkxkC4([0,1]×S1), kx−1kC4( ¯Ω)kψkC4( ¯Ω).
Here in (3.5) and the rest of the paper, we denote by C · · · a positive constant that depends on the arguments inside the bracket.
Proof. We apply the reduced Laplacian (2.25) to ψ on x(0, ξ2) and split it into two parts:
(3.6) △′hψ x(0, ξ2) = 2 √ g(0, ξ2) IL,+(0, ξ 2) + I L,−(0, ξ2)
where IL,±(0, ξ2) are the terms associated with q
αβ (h2, ξ2± h 2) = √ gg αβ (h2, ξ2± h 2), respectively. More precisely, IL,+(0, ξ2), 1 h q 11 D1ψ +˜ q 12 D2ψ˜ (h 2, ξ 2+h 2) = 1 h q11(h 2, ξ 2+h 2) ∂1ψ(˜ h2, ξ2+h2) + h2 24 ∂31ψ(˜ h2, ξ 2+h 2) + K h L,1,+(ψ)(ξ2)h 3 + 1 h q12(h2, ξ2+h2)∂ 2ψ(˜ h2, ξ2+h2) + h2 24 ∂32ψ(˜ h2, ξ 2+h 2) + K h L,2,+(ψ)(ξ2)h 3 (3.7) 6
where ˜ψ = ψ ◦ x and the remainder terms satisfy |Kh
L,ℓ,+(ψ)(ξ2)| ≤ CkxkC4([0,1]×S1)kψkC4( ¯Ω), and we have
assumed that analytic metric tensors have been adopted in the discretization. The analysis below applies to numerical metric tensors without difficulty.
Expand (3.7) around (0, ξ2) and apply (2.2), we get
IL,+(0, ξ2) =1 h q11∂ 1ψ +˜ q 12∂ 2ψ(0, ξ˜ 2) +1 2 ∂1 q 11∂ 1ψ +˜ q 12∂ 2ψ(0, ξ˜ 2) + ˜PL,+(ψ)(ξ2)h + ˜ Qh L,+(ψ)(ξ2)h 2 (3.8) where | ˜Qh L,+(ψ)(ξ2)| ≤ CkxkC4([0,1]×S1), kx−1kC4( ¯Ω)kψkC4( ¯Ω) and (3.9) P˜L,+(ψ)(ξ2), 1 8 ∂21 q 11 ∂1ψ +˜ q 12 ∂2ψ(0, ξ˜ 2) + 1 24 q 11 ∂31ψ +˜ q 12 ∂32ψ(0, ξ˜ 2). Similarly, IL,−(0, ξ2),−1 h q21D 1ψ +˜ q 22D 2ψ˜ (h2, ξ2−h2) = −1 h q21∂ 1ψ +˜ q 22∂ 2ψ(0, ξ˜ 2) + 1 2 ∂2 q 21∂ 1ψ +˜ q 22∂ 2ψ(0, ξ˜ 2) − ˜PL,−(ψ)(ξ2)h + ˜Qh L,−(ψ)(ξ2)h 2 (3.10) where | ˜Qh L,−(ψ)(ξ2)| ≤ CkxkC4([0,1]×S1), kx−1kC4( ¯Ω)kψkC4( ¯Ω) and (3.11) P˜L,−(ψ)(ξ2), 18∂ 2 2 q 21∂ 1ψ +˜ q 22∂ 2ψ(0, ξ˜ 2) + 1 24 q 21∂ 3 1ψ +˜ q 22∂ 3 2ψ(0, ξ˜ 2).
From (2.1), (2.8) and the relation
(3.12) g 11 − 2g 21+ g22= 2g11, g 11 −g 22= 2g12, we have (3.13) 1 h q11∂ 1ψ +˜ q 12∂ 2ψ −˜ q 21∂ 1ψ −˜ q 22∂ 2ψ˜ = √ 2√g h g 11∂ 1ψ + g12∂2ψ.
As a consequence, (3.4) follows from (3.6), (3.8) and (3.10) with PL(ψ)(ξ2),√ 2 g(0, ξ2) P˜L,+(ψ)(ξ 2 ) − ˜PL,−(ψ)(ξ2), QhL(ψ)(ξ2), 2 √ g(0, ξ2) Q˜ h L,+(ψ)(ξ2) + ˜QhL,−(ψ)(ξ2).
The estimate for PR and QhRare similar.
From Lemma 3.1, it is clear that |ωe(t, ·) − ωa(t, ·)| = O(h2) provided ϕ satisfies
ϕ t, x(0, ξ2) = 0, ∂ξ1ϕ t, x(0, ξ2) = − PL(ψe)(t, ξ2) 2√2g11(0, ξ2) , ϕ t, x(1, ξ2) = 0, ∂ξ1ϕ t, x(1, ξ2) = PR(ψe)(t, ξ2) 2√2g11(1, ξ2). (3.14) 7
with PL,R(ψe) obtained by applying Lemma 3.1 to the exact stream function ψe(t, ·).
Such a correction ϕ could be constructed by combining direct tensor products of fL,Rwith proper cutoff
functions in ξ1. However, for regularity consideration, we will elaborate further by mollifying f
L and fR in
the ξ2 direction, with support of the mollifier proportional to the distance to Γ. The singularity induced
by vanishing support of the mollifier near the boundary can be compensated by the condition ϕ = 0 on Γ. More precisely, we have the following Lemma:
Lemma 3.2. Given fL and fR∈ C1(S1), there exists a function ˜ϕ ∈ C2([0, 1] × S1) such that
(3.15) ϕ(0, ξ˜ 2) = 0, ϕ(1, ξ˜ 2) = 0,
(3.16) ∂1ϕ(0, ξ˜ 2) = fL(ξ2), ∂1ϕ(1, ξ˜ 2) = fR(ξ2),
and
(3.17) k ˜ϕkC2([0,1]×S1)≤ C kfLkC1(S1)+ kfRkC1(S1).
Proof. Let η(·) : S17→ R be a standard mollifier with compact support [−δ, δ] ⊂ S1and total mass 1,
Z S1 η(ξ2) dξ2= 1, ηǫ(ξ2), 1 ǫη( ξ2 ǫ ). We define ˜ ϕL(ξ1, ξ2) = ξ1 ηξ 1 ∗ fL(ξ2) = Z S1 η(ξ 2− λ ξ1 )fL(λ) dλ, (3.18) ˜ ϕR(ξ1, ξ2) =(1 − ξ1) η1−ξ 1 ∗ fR(ξ2) = Z S1 η(ξ 2− λ 1 − ξ1)fR(λ) dλ, (3.19) and (3.20) ϕ(ξ˜ 1, ξ2) = ˜ϕL(ξ1, ξ2)Θ(ξ1) − ˜ϕR(ξ1, ξ2)Θ(1 − ξ1),
where Θ is a smooth cutoff function satisfying
(3.21) Θ(ξ1) = 1, 0 ≤ ξ1≤ 1 3 0, 23 ≤ ξ1≤ 1 smoothly connected on 1 3 ≤ ξ 1≤2 3
It is easy to see that ˜ϕ satisfies (3.15). To show that ˜ϕ satisfies (3.16), we note that
∂ξ1(ξ1ηξ 1 (ξ2)) =∂ξ1η ξ2 ξ1 = − ξ2 (ξ1)2η′ ξ2 ξ1 = 1 ξ1ζ ξ2 ξ1 = ζ ξ1 (ξ2), (3.22)
where ζ(ξ2), −ξ2η′(ξ2), is another mollifier with total mass 1 and ζǫ(ξ2),1 ǫζ( ξ2 ǫ). Therefore lim ξ1→0+∂ξ 1ϕ˜L(ξ1, ξ2) = lim ξ1→0+ Z S1 ∂ξ1(ξ1ηξ 1 )(ξ2− λ)fL(λ) dλ = lim ξ1→0+ Z S1 ζξ1(ξ2− λ)fL(λ) dλ = fL(ξ2). (3.23) 8
Similarly,
(3.24) lim
ξ1→1−∂ξ1ϕ˜R(ξ
1, ξ2
) = −fR(ξ2).
To estimate the C2-norm of ˜ϕ(ξ1, ξ2), we first observe that
∂ξ22ϕ˜L(ξ1, ξ2) = Z S1 ∂ξ22η( ξ2− λ ξ1 )fL(λ) dλ = Z S1 ∂λ∂ξ2η(ξ 2 − λ ξ1 )fL(λ) dλ = Z S1 ∂ξ2η( ξ2− λ ξ1 )∂λfL(λ) dλ ≤ kfLkC1(S1)kη′kL1(S1). (3.25) Secondly ∂ξ1ζξ 1 (ξ2− λ) = ∂ ξ1 1 ξ1ζ( ξ2− λ ξ1 ) = − 1 ξ1 1 ξ1ζ( ξ2− λ ξ1 ) + ξ2− λ (ξ1)2 ζ′( ξ2− λ ξ1 ) =1 ξ1∂λ ξ2− λ ξ1 ζ( ξ2− λ ξ1 ) = ∂λZ ξ1 (ξ2− λ) (3.26) where (3.27) Z(ξ2), ξ2ζ(ξ2), Zǫ(ξ2),1 ǫZ ξ2 ǫ From (3.22) and (3.26), we have
|∂ξ21ϕ˜L(ξ1, ξ2)| = Z S1 ∂ξ1ζξ 1 (ξ2− λ)fL(λ) dλ = Z S1 ∂λZξ 1 (ξ2− λ)fL(λ) dλ = Z S1 Zξ1(ξ2− λ)∂λfL(λ) dλ ≤ kfLkC1(S1)kZkL1(S1). (3.28)
Similar calculation leads to
|∂ξ1ϕ˜L(ξ1, ξ2)| ≤ kfLkC0(S1), (3.29) |∂ξ2ϕ˜L(ξ1, ξ2)| ≤ kfLkC0(S1)kη′kL1(S1), (3.30) |∂ξ1∂ξ2ϕ˜L(ξ1, ξ2)| ≤ kfLkC1(S1). (3.31)
The estimate of ˜ϕR(ξ1, ξ2) is also similar. Therefore (3.16) follows.
Proof of Theorem 1. Let ˜ϕ be given by Lemma 3.2 with (3.32) fL(t, ξ2), −PL(ψ
e)(t, ξ2)
2√2g11(0, ξ2) , fR(t, ξ
2), PR(ψe)(t, ξ2)
2√2g11(1, ξ2).
and ϕ = ˜ϕ ◦ x−1. It follows from (3.5) that
kϕ(t, ·)kC2( ¯Ω)≤CkxkC4([0,1]×S1), kx−1kC4( ¯Ω)
kPL(ψe)kC1(S1)+ kPR(ψe)kC1(S1)
≤CkxkC4([0,1]×S1), kx−1kC4( ¯Ω)kψe(t, ·)kC4( ¯Ω).
(3.33)
From the construction of ϕ, it is also clear that
k∂tϕ(t, ·)kC2( ¯Ω)≤CkxkC4([0,1]×S1), kx−1kC4( ¯Ω) k∂tPL(ψe)kC1(S1)+ k∂tPR(ψe)kC1(S1) ≤CkxkC4([0,1]×S1), kx−1kC4( ¯Ω)k∂tψe(t, ·)kC4( ¯Ω). (3.34) 9
Now we define (3.35) ua(t, ·) ,∇ ⊥ h(ψe+h 2 ϕ) ∈ L2(Ω c, R2)
The corresponding approximate vorticity is given by
(3.36) ωa(t, ·) , ∇⊥′h · ua= △′h(ψe+h 2 ϕ) ∈ L2( ¯Ωg, R). It follows that (3.37) ωa = ωe+P L,R(ψe) ± 2 √ 2g11∂ 1ϕ h +Rh 4(ψe) + Rh2(ϕ) h2 on Γg ωe+ ˚Rh 4(ψe) + ˚Rh2(ϕ) h2 on ˚Ωg
where the remainder terms satisfy kRh
k(·)kC0( ¯Ω), k˚Rhk(·)kC0( ¯Ω)≤ CkxkCk([0,1]×S1), kx−1kCk( ¯Ω)k · kCk( ¯Ω).
For brevity, we shall omit the dependence of C on the mapping x from now on. From (3.35), (3.36), (3.37), (3.14) and (2.33), we conclude that
(3.38) |ua(t, ·) − ue(t, ·)| ≤ C1kψe(t, ·)kC3( ¯Ω), kϕ(t, ·)kC1( ¯Ω) h2= C1kue(t, ·)kC2( ¯Ω) h2, (3.39) |∂tua(t, ·) − ∂tue(t, ·)| ≤ C2k∂tψe(t, ·)kC3( ¯Ω), k∂tϕ(t, ·)kC1( ¯Ω) h2= C2k∂tue(t, ·)kC2( ¯Ω) h2, (3.40) |ωa(t, ·) − ωe(t, ·)| ≤ C3kψe(t, ·)kC4( ¯Ω), kϕ(t, ·)kC2( ¯Ω) h2= C3kue(t, ·)kC3( ¯Ω) h2 and (3.41) ∇′h· uh(t, ·) = ∇′h· ua(t, ·) = 0.
We can now write
(3.42) ∂tua+ ¯ωeua⊥+
∇hpe= ν
∇⊥hωe+ E + f on Ωc,
where E is the local truncation error:
(3.43) E = ∂t(ua− ue) + (¯ωeua⊥− ωeue⊥) + (∇ h− ∇)pe+ ν(∇⊥−∇ ⊥ h)ωe. From (3.38, 3.39) and (3.44) ω¯eua⊥− ωe(ue)⊥= ¯ωe(ua − ue)⊥+ (¯ωe − ωe)(ue)⊥,
it is easy to see that
(3.45) Z T 0 kE k2 Ωc≤ C4ν, ku ek H1(0,T ;C2( ¯Ω)), kuekL2(0,T ;C4( ¯Ω)), kpekL2(0,T ;C3( ¯Ω)) h4.
We now proceed to derive an error equation. Define ε(t, ·) = uh(t, ·) − ua(t, ·), we have
(3.46) ∂tε+ (¯ωhuh⊥− ¯ωeua⊥) +
∇h(ph− pe) = ν
∇⊥h(ωh− ωe) − E .
From (2.29) and (3.41), we have h ε , ∇h(ph− pe) iΩc= −h ∇′h· ε , ph− peiΩ¯g = 0, and thus h ε , ∂tεiΩc+ h ε , ¯ω h uh⊥− ¯ωeua⊥iΩc = νh ε , ∇⊥h(ωh− ωe) iΩc− h ε , E iΩc. From (2.30), νh ε ,∇ ⊥ h(ωh− ωe) iΩc = − νh ∇⊥′h · (uh− ua) , ωh− ωeiΩ¯g = − νh ωh− ωa, ωh− ωeiΩ¯g = − νh ωh− ωe, ωh− ωeiΩ¯g− νh ωe− ωa, ωh− ωeiΩ¯g ≤ −ν2kωh− ωek2Ω¯g + ν 2kω e − ωak2Ω¯g.
Since ε = uh− ua is pointwise perpendicular to ¯ωh(uh− ua)⊥, we have
h ε , ¯ωhuh⊥− ¯ωeua⊥iΩc= h ε , ¯ω h(uh − ua)⊥+ (¯ωh − ¯ωe)ua⊥iΩc= h ε , (¯ω h − ¯ωe)ua⊥iΩc.
Thus from Jensen’s inequality, h ε , ¯ωhuh⊥− ¯ωeua⊥iΩc ≤ 1 νku a k2L∞(Ω c)kεk 2 Ωc+ ν 4kω h − ωek2Ω¯g ≤ 1 νku a k2L∞(Ω c)kεk 2 Ωc+ ν 4k¯ω h − ¯ωek2Ωc
In summary, we have shown that
(3.47) 1 2∂tkεk 2 Ωc+ ν 4kω h− ωek2 ¯ Ωg ≤ 1 2 + 1 νku ak2 L∞(Ω c) kεk2 Ωc+ ν 2kω e− ωak2 ¯ Ωg + 1 2kE k 2 Ωc.
In view of (3.40), (3.45), and Gronwall’s inequality, we have
(3.48) max 0≤t≤Tku h − uak2Ωc+ ν Z T 0 kω h − ωek2Ω¯g≤ K1(kuh0− ue0k2Ωc+ h4). The estimate (3.1) follows in view of (3.38).
Next, we proceed to show the preliminary pressure error estimate (3.2). Take the inner product with
∇h(ph− pe) on both sides of (3.46), then apply (2.29) and (3.41), we obtain
h ∇h(ph− pe) , ¯ωhuh⊥− ¯ωeua⊥iΩc+ k ∇h(ph− pe)k2Ωc =νh ∇h(ph− pe) , ∇⊥h(ωh− ωe) iΩc− h ∇h(ph− pe) , E iΩc. (3.49) Thus (3.50) 1 4k ∇h(ph− pe)k2Ωc≤ ν 2 k ∇⊥h(ωh− ωe)k2Ωc+ k¯ω h uh⊥− ¯ωeua⊥k2Ωc+ kE k 2 Ωc. Since k∇ ⊥ h(ωh− ωe)k2Ωc≤ 8 h2kω h− ωek2 ¯ Ωg and (3.51) k¯ωhuh⊥− ¯ωeua⊥k2 Ωc≤ 2ku hk2 L∞(Ω c)kω h− ωek2 ¯ Ωg+ 2kω ek2 C0( ¯Ω)ku h− uak2 Ωc, 11
it follows that 1 4 Z T 0 k ∇h(ph− pe)k2Ωc≤ 8ν2 h2 + 2 max [0,T ]ku h k2L∞(Ω c) Z T 0 kω h − ωek2Ω¯g + 2 max [0,T ]kω e k2C0( ¯Ω) Z T 0 ku h − uak2Ωc+ Z T 0 kE k 2 Ωc. Note that kuhk2L∞(Ωc)≤2kuh− uek2L∞(Ωc)+ 2kuek2C0( ¯Ω) ≤ 2 h2min Ωc √g h kuh− uek2Ωc+ 2ku e k2C0( ¯Ω) ≤min2K1 Ωc √g h (h −2 kuh0− ue0k2Ωc+ h2) + 2kuek2C0( ¯Ω). (3.52) Hence (3.53) max [0,T ]ku h k2L∞(Ω c)≤ C5(1 + h−2kuh0− ue0k2Ωc), where C5= C5T, ν, kuekL2(0,T ;C4( ¯Ω)), kuekH1(0,T ;C2( ¯Ω)), kpekL2(0,T ;C3( ¯Ω)). Consequently, (3.54) Z T 0 k ∇h(ph− pe)k2Ωc≤ K ′ 2 h −2 kuh0− ue0k2Ωc(1 + ku h 0− ue0k2Ωc) + h2 ≤ K2 h −2 kuh0− ue0k2Ωc+ h2.
Remark 1. The error estimate (3.1), (3.2) is subject to appropriate approximation of the initial velocity field. This can be achieved, for example, by taking uh
0 = ua0= ua(0, ·) with ua defined by (3.35). In view of
(3.38), uh 0 = ua0 gives max [0,T ]ku h − uek2Ωc+ ν Z T 0 kω h − ωek2Ω¯g ≤ K1h 4 , (3.55) Z T 0 k ∇h(ph− pe)k2Ωc≤ K2 h2. (3.56)
Similarly, the refined estimates (3.58), (3.59) in Theorem 3 also depend on the approximations of initial velocity uh
0 − ue0 and initial vorticity ωh0 − ω0e = ∇
⊥
h(uh0 − ue0). The choice uh0 = ua0 and consequently
ωh 0 =
∇⊥hua0 gives full second order accuracy in (3.58), (3.59) in view of (3.38) and (3.40).
3.2. The LBB Condition and Refined Pressure Estimate. Since the pressure lacks an evolu-tionary equation, it is in general difficult to obtain optimal estimate from basic energy estimate alone. The pressure estimate (3.2) is only first order. To get full second order accuracy, we resort to the well known inf-sup condition (also known as div-stability condition or Ladyzhenskaya-Babuˇska-Brezzi (LBB) condition). Indeed, we have
Theorem 2 (LBB). If the mapping x : (ξ1, ξ2) 7→ (x, y) is a C1,1 bijection between [0, 1] × S1 and
¯
Ω ⊂ R2, then there exists a constant β > 0 such that
(3.57) inf p∈L2( ¯Ωg,R)/R2 u sup ∈L2(Ωc,R2) h p ,∇ ′ h· u iΩ¯g kpkΩ¯gkuk H1h ≥ β uniformly in h. 12
Here kukH 1 h , (kuk 2 Ωc+ k ∇⊥′h · uk2Ω¯g+ k ∇′h· uk2Ω¯g) 1
2 is the natural norm associated with the discrete
vector Laplacian.
The LBB condition arises naturally as a compatibility condition between the discrete velocity and pres-sure spaces in mixed finite element formulations. It is a fundamental research topic in finite element analysis for steady state computation, yet rarely discussed in finite difference setting or dynamical problems. The uniform estimate (3.57) is not only vital to the pressure error estimate presented here, but also directly affects the condition number and uniform bound of the solution operator for the Stokes’s problem.
The verification of the LBB condition is quite complicated, especially for low order finite element meth-ods. Denote by DOF(u) and DOF(p) the degree of freedom for u and p, respectively. It is clear that, the larger the ratio DOF(u)/DOF(p) is, the more likely (3.57) is to hold and be verified.
To verify the LBB condition, a common approach is to reduce the problem to a patch of elements and conduct local analysis [BN1]. See also [St] for a related approach. This local argument has proved successful for higher order finite element methods with large DOF(u)/DOF(p) ratio. When the ratio is low, the spatial compatibility plays a crucial role in establishing (3.57).
For example, the spaces (Ql, Pl−1) are well known and widely used in mixed finite element formulation.
Using the local argument, it can be shown that LBB condition holds for (Ql, Pl−1) with l ≥ 2 [GR]. However,
this local argument does not apply to the lowest order scheme Q1− P0 (bilinear in velocity components and
piecewise constant in pressure), but only to some of its variants equipped with extra degrees of freedom in the velocity space. In fact, the Q1− P0element is known to violate the LBB condition with β = O(h) [BN2].
In our case, gMAC is staggered and supported the same way as the Q1− P0element but only differs in
the discretization of the vector Laplacian. Both of them have the minimal ratio DOF(u)/DOF(p) = 2 (3, if in 3D) among quadrilateral meshes. The similarity between gMAC and the Q1− P0 element demonstrates
the subtlety and difficulty of Theorem 2. It is also a common belief that staggered grids in some sense implies the inf-sup condition. The contrast between gMAC scheme and the Q1− P0 element indicates the
inf-sup estimate is subtler than plain staggeredness. Our analysis shows that it is in fact encoded in the compatibility of the spatial discretization (2.29)-(2.32). See the Appendix for details.
With the inf-sup condition (3.57) established, the pressure error estimate can be improved to second order with only minor extra regularity requirement on the exact solution:
Theorem 3. Assume the mapping x : (ξ1, ξ2) 7→ (x, y) is a C4 bijection from [0, 1] × S1 to ¯Ω ⊂ R2,
and ue∈ L2(0, T ; C4( ¯Ω)) ∩ H1(0, T ; C3( ¯Ω)), pe∈ L2(0, T ; C3( ¯Ω)). Then ν max [0,T ]kω h − ωek2Ω¯g+ Z T 0 k∂ t(uh− ue)k2Ωc≤ K3( h−2kuh0− ue0k4Ωc+ ku h 0− ue0k2Ωc+ kω h 0− ωe0k2Ω¯g+ h4), (3.58) Z T 0 kp h − pek2Ω¯g≤ K4( h−2kuh0− ue0k4Ωc+ ku h 0− ue0k2Ωc+ kω h 0− ωe0k2Ω¯g+ h4) (3.59)
where K3, K4are constants that depend only on T , ν, kuekH1(0,T ;C3( ¯Ω)),kuekL2(0,T ;C4( ¯Ω)), kpekL2(0,T ;C3( ¯Ω))
and independent of h.
Proof. We start with the estimate (3.58). Take the inner product with ∂tεon both sides of (3.46) and
apply (3.41), we get k∂tεk2Ωc+ h ∂tε, ¯ω huh⊥− ¯ωeua⊥i Ωc= −νh ∂t(ω h − ωa) , ωh− ωeiΩ¯g− h ∂tε, E iΩc. 13
Note that, from (3.51), ν 2∂tkω h − ωek2Ω¯g+ k∂tεk2Ωc ≤|h ∂tε, ¯ωhuh⊥− ¯ωeua⊥iΩc| + ν|h ∂t(ω e − ωa) , ωh− ωeiΩ¯g| +1 4k∂tεk 2 Ωc+ kE k 2 Ωc ≤12k∂tεk2Ωc+ k¯ω h uh⊥− ¯ωeua⊥k2Ωc+ ν 2k∂t(ω e − ωa)k2Ω¯g+ν 2kω h − ωek2Ω¯g+ kE k2Ωc ≤12k∂tεk2Ωc+ 2ku h k2L∞(Ω c)+ ν 2kω h − ωek2Ω¯g+ 2kωek2C0( ¯Ω)ku h − uak2Ωc+ ν 2k∂t(ω e − ωa)k2Ω¯g + kE k2Ωc. From (3.36), we have (3.60) k∂tωe(t, ·) − ∂tωa(t, ·)k2Ω¯g≤ C6k∂tψ e (t, ·)k2C4( ¯Ω), k∂tϕ(t, ·)k2C2( ¯Ω) h4= C6k∂tue(t, ·)k2C3( ¯Ω) h4. Thus, in view of (3.48), we obtain
ν∂tk(ωh− ωe)(t, ·)kΩ2¯g + k∂tε(t, ·)k2Ωc ≤ 4kuhk2L∞(Ω c)+ νkω h − ωek2Ω¯g+ C7ν, K1, kω e kC0( ¯Ω), C6 kuh0− ue0k2Ωc+ h4 + kE k2Ωc. (3.61)
Integrate (3.61) over [0, T ] and apply (3.1), (3.45) and (3.53), we have
(3.62) Z T 0 k∂ tεk2Ωc+ ν max[0,T ]kω h − ωek2Ω¯g ≤ K3( h−2kuh0− ue0k4Ωc+ ku h 0− ue0k2Ωc+ kω h 0− ω0ek2Ω¯g+ h4).
Thus (3.58) follows in view of (3.39).
We now proceed with the pressure error estimate (3.59). Firstly, both phand peare unique up to proper
normalizations. The pressure norm in (3.59) should be understood as kph− pekΩ¯g, min
ph−˜p∈ker(
∇⊥h)
k˜p − pekΩ¯g
Without loss of generality, we assume that Z
Ω
pedx = 0 and ph− pe∈ L2( ¯Ω
g, R)/R2.
Since the assumption of Theorem 2 is weaker than that of Theorem 3, it follows that (3.57) holds and there exists u ∈ L2(Ω c, R2) such that (3.63) −h ∇h(ph− pe) , u iΩc= h p h − pe, ∇′h· u iΩc ≥ βk(p h − pe)kΩ¯gkuk H1h.
Now take the inner product with −u on both sides of (3.46), we get (3.64) h −u , ∂tεiΩc+ h −u , ¯ω h uh⊥− ¯ωeua⊥iΩc− h ∇h(ph− pe) , u iΩc = νh −u , ∇⊥h(ωh− ωe) iΩc+ h u , E iΩc. Therefore βk(ph− pe)k ¯ ΩgkukH 1 h ≤|h u , ∂tεiΩc| + |h u , ¯ω huh⊥− ¯ωeua⊥i Ωc| + |h u , E iΩc| + ν|h ∇⊥′h · u , ωh− ωeiΩ¯g| ≤kukH 1 h k∂tεkΩc+ k¯ω h uh⊥− ¯ωeua⊥kΩc+ kE kΩc+ νkω h − ωekΩ¯g (3.65) 14
Consequently, in view of (3.51), we obtain (3.66) kph− pek2Ω¯g≤ 4 β2 k∂tεk2Ωc+ 2ku h k2L∞(Ω c)+ ν 2 kωh − ωek2Ω¯g+ 2kωek2C0( ¯Ω)ku h − uak2Ωc+ kE k 2 Ωc
Now we integrate from 0 to T on both sides of (3.66) and apply (3.1, 3.39, 3.45, 3.53, 3.58). The refined pressure estimate (3.59) then follows.
Similarly, one can give an L∞(0, T ; L2( ¯Ω
g)) estimate for the pressure error with higher regularity
re-quirement on the exact solution. We state the following Theorem without proof.
Theorem 4. Assume the mapping x : (ξ1, ξ2) 7→ (x, y) is a C4 bijection from [0, 1] × S1 to ¯Ω ⊂ R2,
and ue∈ H1(0, T ; C4( ¯Ω)) ∩ H2(0, T ; C2( ¯Ω)), pe∈ H1(0, T ; C3( ¯Ω)). Then (3.67) max [0,T ]k∂t(u h − ue)k2Ωc+ ν Z T 0 k∂ t(ωh− ωe)kΩ2¯g ≤ K5(χ1+ k∂t(u h − ue)0k2Ωc+ h4) (3.68) max [0,T ]kp h− pek2 ¯ Ωg ≤ K6(χ1+ k∂t(u h− ue) 0k2Ωc+ h4) where (3.69) χ1,h −4 kuh0− ue0k6Ωc+ ku h 0− ue0k2Ωc+ h−1kωh0− ωe0k3Ω¯g+ kωh0− ω0ek2Ω¯g
and K5, K6 are constants that depend on T , ν, kuekH1(0,T ;C4( ¯Ω)), kuekH2(0,T ;C2( ¯Ω)) and kpekH1(0,T ;C3( ¯Ω)),
but not on h.
The k∂t(uh− ue)0k2Ωc term in (3.67) and (3.68) result from Ladyzhenskaya type higher order energy
estimate. An alternative expression in terms of k
∇⊥h(ω0h− ω0e)k2Ωc can be derived as follows.
Rewrite (3.46) at t = 0 as (3.70) ν ∇⊥h(ωh0− ωe0) − (¯ωh0uh0⊥− ¯ωe0ua0⊥) − E0= ∂t(uh− ua)0+ ∇h(ph0− pe0). Since (3.71) h ∂t(uh− ua)0, ∇h(ph0− pe0) iΩc= 0,
it follows that (3.70) is an orthogonal decomposition for ν
∇⊥h(ω0h− ωe0) − (¯ωh0uh0⊥− ¯ω0eua0⊥) − E0. Therefore (3.72) k∂t(uh− ua)0kΩc≤ νk ∇⊥h(ω0h− ω0e)kΩc+ k¯ω h 0uh0⊥− ¯ω0eu0a⊥kΩc+ kE0kΩc.
Moreover, from (3.43), (3.44), (3.38) and (3.39), we have (3.73) kE0k2Ωc ≤ C8ku e 0kC4( ¯Ω), k∂tue0kC2( ¯Ω), kpe0kC3( ¯Ω) h4, and (3.74) k¯ωh0uh0⊥− ¯ω0eua0⊥k2Ωc≤ 2ku h 0k2L∞(Ω c)kω h 0 − ωe0k2Ω¯g+ 2kω0ek2C0( ¯Ω)ku h 0− ua0k2Ωc, χ2, with (3.75) χ2≤ Ckue0kC1( ¯Ω)(kω0h− ω0ek2Ω¯g+ ku0h− ue0k2Ωc+ h−2kuh 0− ue0k2Ωckω h 0− ωe0k2Ω¯g). 15
Therefore (3.76) k∂t(uh− ua)0k2Ωc ≤ νk ∇⊥h(ωh0− ω0e)k2Ωc+ χ2+ C8 h4≤ C9ν, C2, C8] k ∇⊥h(ωh0 − ωe0)k2Ωc+ χ2+ h4. As a result, we can replace the k∂t(uh− ue)0k2Ωc term in (3.67) and (3.68) by k
∇⊥h(ω0h− ωe0)k2Ωc+ χ2, which
is easier to analyze in practice.
For example, if the mapping x is a C5 bijection and ue
0∈ C4( ¯Ω) (and consequently ψe0∈ C5( ¯Ω)), then
following the derivation in the proof of Lemma 3.1, it is easy to see that
(3.77) kPL,R(ψ0e)kC2(S1)≤ CkxkC5([0,1]×S1), kx−1kC5( ¯Ω)kψe0kC5( ¯Ω)
and the construction (3.18)-(3.21) gives
(3.78) k ˜ϕkC3([0,1]×S1)≤ CkxkC5([0,1]×S1), kx−1kC5( ¯Ω)(kfLkC2(S1)+ kfRkC2(S1)).
Together with (3.32), we conclude that
(3.79) kϕ0kC3( ¯Ω) ≤ CkxkC5([0,1]×S1), kx−1kC5( ¯Ω)kψ0ekC5( ¯Ω).
where ϕ0= ˜ϕ0◦ x−1.
Overall, we have the following refined estimate
(3.80) ωa0 = ωe 0+ R4(ψe0) + △ϕ0 h2+Rh 5(ψe0) + Rh3(ϕ0) h3 on Γg ωe 0+ R4(ψe0) + △ϕ0 h2+ ˚Rh 5(ψe0) + ˚Rh3(ϕ0) h3 on ˚Ωg, where (3.81) R4(ψ0e) = 1 24√g 2 X γ=1 ∂3γ( p gg 1γ ∂1ψe+ p gg 2γ ∂2ψe) +∂ γ( p gg 1γ ∂31ψe+ p gg 2γ ∂32ψe)
and the remainder terms satisfy
(3.82) kRhk(·)kC0( ¯Ω), k˚Rhk(·)kC0( ¯Ω)≤ CkxkCk([0,1]×S1), kx−1kCk( ¯Ω)k · kCk( ¯Ω).
It is easy to see from (3.80)-(3.82) that
(3.83) k ∇⊥h(ωa0− ω0e)kΩc≤ C10[kxkC5([0,1]×S1), kx −1k C5( ¯Ω), kψe0kC5( ¯Ω)]h 2 .
In view of Remark 1 and the analysis above, the estimates (3.67) and (3.68) result in full second order accuracy provided the initial data is chosen properly. In particular, we have the following
Corollary 1. Assume the mapping x : (ξ1, ξ2) 7→ (x, y) is a C5 bijection from [0, 1] × S1 to ¯Ω ⊂ R2
and uh 0 = ua0. Then (3.84) max [0,T ]k∂t(u h − ue)k2Ωc+ ν Z T 0 k∂ t(ωh− ωe)k2Ω¯g+ max [0,T ]kp h − pek2Ω¯g ≤ K7h 4 , where K7= K7T, ν, kuekH1(0,T ;C4( ¯Ω)), kuekH2(0,T ;C2( ¯Ω))kpekH1(0,T ;C3( ¯Ω)), kxkC5([0,1]×S1), kx−1kC5( ¯Ω).
It is worth noting that the extra regularity requirement on x is only needed in constructing a good initial data. It will be interesting to see if there exists initial data satisfying (3.83) under a C4bijection.
Acknowledgment. This work is partially supported by National Science Foundation, Center of Math-ematical Sciences at NTHU, National Science Council and the National Center for Theoretical Sciences in Taiwan.
REFERENCES
[AN] C. R. Anderson, Derivation and solution of the discrete pressure equations for the incompressible Navier-Stokes equa-tions, preprint LBL-26353, Lawrence Berkeley Laboratory, Berkeley, CA, 1988.
[BD] F. Bertagnolio and O. Daube, Solution of the div-curl problem in generalized curvilinear coordinates, J. Comp. Phys. 138(1997), pp. 121–152.
[BN1] J. M. Boland and R. A. Nicolaides, Stability of finite elements under divergence constraints, SIAM J. Numer. Anal., 20(1983), pp. 722–731.
[BN2] J. M. Boland and R. A. Nicolaides, On the Stability of Bilinear-Constant Velocity-Pressure Finite Elements, Numer. Math., 44 (1984), pp. 219–222.
[BN3] J. M. Boland and R. A. Nicolaides, Stable and Semistable Low Order Finite Elements for Viscous Flows, SIAM J. Numer. Anal., 22 (1985), pp. 474–492.
[FPT] M. Fortin, R. Peyret and R. Temam, R´esolution num´erique des ´equations de Navier-Stokes pour un fluide incompressible. J. M´ecanique, 10 (1971), pp. 357–390.
[Gu] M. D. Gunzburger, Finite Element Methods for Viscous Incompressible Flows: A Guide to Theory, Practice and Algorithms, Academic Press, New York, 1989.
[GR] V. Girault and P. A. Raviart, Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms, Springer-Verlag, New York, 1986.
[HW] T. Y. Hou and B. R. Wetton, Second-order convergence of a projection scheme for the incompressible Navier-Stokes equations with boundaries, SIAM. J. Numer. Anal., 30 (1993), pp. 609–629.
[HLW] Y.-L. Huang, J.-G. Liu and W.-C. Wang, A generalized MAC scheme on curvilinear domains, preprint.
[MAC] F. H. Harlow and J. E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids, 8 (1965), pp. 2182–2189.
[LLP] J.-G. Liu, J. Liu and R. Pego, Stability and convergence of efficient Navier-Stokes solvers via a commutator estimate, Comm. Pure Appl. Math., 60 (2007), pp. 1443–1487.
[Ne] J.-C. N´ed´elec, A new family of mixed finite elements in R3
, Numer. Math., 50 (1986), pp. 57–81.
[RT] P. A. Raviart and J. M. Thomas, A mixed finite element method for 2nd order elliptic problems, in Mathematical Aspects of Finite Element Methods, Lecture Notes in Mathematics, Springer, Berlin, 606 1975, pp. 292–315. [N1] R. A. Nicolaides, Analysis and Convergence of the MAC Scheme I. The Linear Problem SIAM J. Numer. Anal., 29
(1992), pp. 1579–1591.
[N2] R. A. Nicolaides and X. Wu, Analysis and convergence of the MAC scheme II, Navier-Stokes equations, Math. Comp., 65(1996), pp. 29–44.
[HW1] H. Huang and B. R. Wetton, Discrete compatibility in finite difference methods for viscous incompressible fluid flow, J. Comp. Phys., 126 (1996), pp. 468–478.
[St] R. Stenberg, Analysis of mixed finite element methods for the Stokes problem: a unified approach, Math. Comp., 42 (1984), pp. 9–23.
[Wa] W.-C. Wang, A Jump Condition Capturing Finite Difference Scheme for Elliptic Interface Problems, SIAM J. Sci. Comp., 25 (2004), pp. 1479–1496.
[W1] B. R. Wetton. Analysis of the spatial error for a class of finite difference methods for viscous incompressible flow, SIAM J. Numer. Anal., 34 (1997), pp. 723–755.
Appendix A. Proof of Theorem 2, the LBB Condition.
In this appendix, we proceed with the proof of Theorem 2, the LBB condition for gMAC.
Our proof consists of two parts. In section A.1, Theorem 5, we will prove a (formally) stronger version of the LBB condition for the special case of conformal metrics. That is, mappings with corresponding metric tensors of the form
(A.1) g 11 = g 22 = pg (ξ1, ξ2), g 12 = g 21 ≡ 0. For example, the mapping x
: (ξ1, ξ2) 7→ (x , y ) with x = e2πξ1cos(2πξ2), y = e2πξ1sin(2πξ2), satisfies (A.1) with pg
(ξ1, ξ2) = 4π2e4πξ1
.
To make distinction between conformal and general metrics, we will add the subscript ’+’ for all quantities derived from the conformal coordinate mapping, including the basis vectors and the difference operators and various norms in the rest of the paper.
This strong LBB condition, together with a crucial estimate (Lemma A.5), are then used to give a constructive proof for the general case in section A.2. Note that we do not require the general case to be small perturbation of the conformal one.
A.1. Strong Form of LBB Condition for Conformal Metrics. We start with the following strong version of the LBB condition for conformal metrics satisfying (A.1).
Theorem 5 (Strong form of LBB). If the mapping x
: (ξ1, ξ2) 7→ (x , y ) is a conformal bijection between [0, 1] × S1 and ¯Ω ⊂ R 2
, then given any q ∈ L2( ¯Ωg, R)/R
2
, there exists a vector field v ∈ L2(Ωc, R
2) such that (A.2) ∇ ′ h· v = q and (A.3) γkvkH 1 h ≤ kqkΩ¯ g
where γ > 0 is a constant independent of q, v orh.
It is easy to see that (A.2, A.3) implies the standard LBB condition
(A.4) inf q∈L2( ¯Ω g,R)/R 2 v sup ∈L2(Ω c,R 2) h q , ∇′h· v iΩ¯ g kqkΩ¯ g kvkH 1 h ≥ γ uniformly in h.
In fact, it can be shown that the strong form of LBB condition (A.2, A.3) is equivalent to the standard LBB condition (A.4).
Our approach for Theorem 5 is based on a global construction procedure using Fourier series. We start with a list of notations. Define
Cmi , (√ 2 cos(mπξ1 i), 1 ≤ m ≤ N − 1; cos(mπξi1), m = 0, N, (A.5) Smi ,√2 sin(mπξ1i), 1 ≤ m ≤ N − 1, (A.6) Enj , exp(2nπ√−1ξj2), 0 ≤ n ≤ N − 1. (A.7) 18
It is easy to verify that
(A.8) {Cm⊗ En| 0 ≤ m ≤ N, 0 ≤ n ≤ N − 1}
is an orthonormal basis for L2( ¯Ω
g
, R) with respect to the standard inner product
(A.9) ha, bi0, ¯Ω g = h2 N X i=0 ′ N X j=1 (ab)ij.
where the primed sum denotes half weight at i = 0 and i = N . In addition, (A.8) is an orthonormal eigen-basis for pg △ ′
h. Indeed, since the coordinate mapping is conformal and h1= h2 = h, it follows that
g 11= g 22= 1/q g = 1/pg andg 12= g 21 ≡ 0. Thus (A.10) △ ′ h= 1 pg D12′+ D22′ = 1 pg D21′+ D22+ h2 2 D 2 1′D22. where (A.11) (D21′f )i= 2 h2(f1− f0) i = 0 1 h2(fi+1− 2fi+ fi−1) 1 ≤ i ≤ N − 1 2 h2(fN −1− fN) i = N , (D22g)j= gj+1− 2gj+ gj−1 h2 , 1 ≤ j ≤ N.
The reduction is not needed in the ξ2 direction due to periodicity.
It is east to see that
(A.12) D2 1′Cmi = −λ2mCmi , 0 ≤ i ≤ N, D2 2Enj = −λ22nEnj, 1 ≤ j ≤ N, where (A.13) λm,2 sin(mπh/2) h . Therefore (A.14) pg △ ′ h(C m ⊗ En)ij= −κ 2 mn(Cm⊗ En)ij, where (A.15) κ 2 mn, λ2m+ λ22n− h2 2 λ 2 mλ22n= 4 h2 sin 2(mπh 2 ) cos 2(nπh) + cos2(mπh 2 ) sin 2(nπh).
On the other hand,
(A.16) {Sm⊗ En| 1 ≤ m ≤ N − 1, 0 ≤ n ≤ N − 1} is an orthonormal basis for
(A.17) L20( ¯Ωg, R), ψ ∈ L 2( ¯Ω g, R) ψ0,j = 0 = ψN,j, 1 ≤ j ≤ N 19
with (A.18) D21Smi = −λ2mSmi , 1 ≤ i ≤ N − 1, and (A.19) pg △ h (Sm⊗ En)ij= −κ 2 mn(Sm⊗ En)ij, 1 ≤ i ≤ N − 1.
Note that, however, (A.16) is not an eigen-basis for pg
△ ′ h since pg △ ′ h(S m⊗ En) 6= 0 on Γ g.
Proof of Theorem 5. Given q ∈ L2( ¯Ω
g
, R)/R2, we will construct explicitly a vector field v that satisfies
(A.2) and (A.3). This is done in the following steps: Step 1: Solve for φ ∈ L2( ¯Ω g , R)/R2 from (A.20) △ ′ hφ = q on ¯Ωg .
From (A.10-A.15) above, it is easy to see that we can solve (A.20) by expanding φ and Q = pg
q with respect to the eigen-basis (A.8),
(A.21) φij = N −1 X n=0 N X m=0 ˆ φmnCmi Enj, Qij = (pg q)ij = N −1 X n=0 N X m=0 ˆ QmnCmi Enj.
and compare the coefficients mode by mode to get
(A.22) φˆmn= − 1 κ2mnQˆmn, (m, n) 6= (0, 0), (N, N 2) 0, otherwise. Note that from (A.15), κ
2
mn = 0 if and only if (m, n) = (0, 0) or (N,N2), while ˆQ0,0 = ˆQN,N
2 = 0 as
q ∈ L2( ¯Ω
g
, R)/R2. Thus (A.22) indeed gives the unique solution φ in L2( ¯Ω
g
, R)/R2. It is worth mentioning
that, from the analysis above, we have for φ ∈ L2( ¯Ω
g , R)/R2, h2 N X i=0 ′ N X j=1 φ2ij= N X m=0 N −1 X n=0 | ˆφmn|2≤ 1 κ4min N X m=0 N −1 X n=0 |κ 2 mnφˆmn|2 = h 2 κ4min N X i=0 ′ N X j=1 (pg △ ′ hφ) 2 ij (A.23) where (A.24) κ 2 min, min (m,n)6=(0,0),(N,N 2) κ2mn
Moreover, with straightforward calculation (see also (A.76, A.80) below), it is easy to see that
(A.25) κ 2 min=κ 2 1,0= 4 h2sin 2πh 2 = O(1), 8 ≤κ 2 min< π2.
Thus we have the following estimate for the solution φ ∈ L2( ¯Ω
g , R)/R2: (A.26) kφk2Ω¯ g ≤ g max κ4min k △ ′ hφk 2 ¯ Ω g . 20
Following similar calculations, one can also show that (A.27) h2 N X i=0 ′ N X j=1 φ2ij = N X m=0 N −1 X n=0 | ˆφmn|2≤ 1 κ2min N X m=0 N −1 X n=0 κ2mn| ˆφmn|2= 1 κ2min k ∇hφk 2 Ω c and (A.28) k ∇hφk 2 Ω c = N X m=0 N −1 X n=0 κ2mn| ˆφmn|2≤ 1 κ2min N X m=0 N −1 X n=0 |κ 2 mnφˆmn|2= h2 κ2min N X i=0 ′ N X j=1 (pg △ ′ hφ) 2 ij. Similarly, for ψ ∈ L2 0( ¯Ωg , R), we have (A.29) h2 N −1 X i=1 N X j=1 ψij2 = N −1 X m=1 N −1 X n=0 | ˆψmn|2≤ 1 κ2min N −1 X m=1 N −1 X n=0 κ2mn| ˆψmn|2= 1 κ2min k(pg )−1/2 ∇ ⊥ hψk 2 Ωc (A.30) k(pg )−1/2 ∇ ⊥ hψk 2 Ω c = N −1 X m=1 N −1 X n=0 κ2mn| ˆψmn|2≤ 1 κ2min N −1 X m=1 N −1 X n=0 |κ 2 mnψˆmn|2= h 2 κ2min N −1 X i=1 N X j=1 (pg △ h ψ)2ij.
Here in (A.29, A.30), ˆψmn is the coefficient with respect to the basis (A.16):
(A.31) ψij = N −1 X m=1 N −1 X n=0 ˆ ψmnSmi Enj, ψ ∈ L20( ¯Ωg , R).
As a result, we have the following Poincar´e type inequalities for φ ∈ L2( ¯Ω
g , R)/R2 and ψ ∈ L2 0( ¯Ωg , R) respectively: Lemma A.1. 1. Let φ ∈ L2( ¯Ω g , R)/R2, then (A.32) kφk2Ω¯ g ≤ pg max κ2min k ∇hφk 2 Ω c , k ∇hφk 2 Ω c ≤ pg max κ2min k △ ′ hφk 2 ¯ Ω g . 2. Let ψ ∈ L2 0( ¯Ωg , R), then (A.33) kψk˚2Ω g ≤ pg max κ2minpg min k ∇ ⊥ hψk 2 Ω c , k ∇ ⊥ hψk 2 Ω c ≤ g max κ2min k △ h ψk˚2Ω g .
Note that in (A.33),
(A.34) k △ h ψk2˚Ω g , h2N −1X i=1 N X j=1 pg ( △ h ψ)2 ij,
the integrand is summed over interior grids only. In general,
△ ′ hψ|Γg 6= 0 even if ψ ∈ L 2 0( ¯Ωg , R). Step 2:
With φ obtained in step 1, construct ψ ∈ L2 0( ¯Ωg , R) such that (A.35) D− nψ Γg= ˜Dτφ ˚ Γg 21
and (A.36) k △ h ψk˚Ω g ≤ Ck △ ′ hφkΩ¯g ,
where ˚Γg= {(ξi1, ξj2) | i = 1 or N − 1, 1 ≤ j ≤ N }, C is a constant that only depends on pg
max
and pg
min
, D−
n is the one-sided backward difference with respect to the unit outer normal n,
(A.37) D− nψi,j = ψ0,j− ψ1,j h , i = 0, ψN,j− ψN −1,j h , i = N ;
and ˜Dτ the long-stencil centered difference with respect to the counter-clockwise unit tangent τ ,
(A.38) D˜τφi,j= φ1,j−1− φ1,j+1 2h , i = 1, φN −1,j+1− φN −1,j−1 2h , i = N − 1.
This step is the most technical part of the Theorem. We will detail it in section A.1.1. Step 3:
With ψ given by Step 2, construct explicitly the vector field v as
(A.39) v , ∇hφ − ∇ ⊥ hψ. Since ψ ∈ L2 0( ¯Ωg
, R), it follows from (2.33) that
∇ ′ h· ∇ ⊥ hψ ≡ 0 on ¯Ωg
. Therefore from (2.31), we have
(A.40) ∇ ′ h· v = ∇ ′ h· ∇hφ − ∇ ′ h· ∇ ⊥ hψ = △ ′ hφ = q on ¯Ωg . This gives (A.2).
To see that v indeed satisfies (A.3), we first note from (A.32), (A.33), (A.36) and (A.20) that
(A.41) kvk2Ω c = k ∇hφk 2 Ω c + k ∇ ⊥ hψk 2 Ω c ≤ 1 κ2min (pg max k △ ′ hφk 2 ¯ Ω g + g max k △ h ψk˚2Ω g ) ≤ β1kqk2Ω¯ g
for some constant β1> 0 that only depends on pg
max
. On the other hand, from (2.31) and (2.32), it is easy to see that (A.42) pg ∇ ⊥′ h · v = pg ∇ ⊥′ h · ( ∇hφ − ∇ ⊥ hψ) = pg △ ′ hψ on ˚Ωg .
Note that in general, pg
∇ ⊥′ h · ∇hφ 6= 0 on Γ
g. The net contribution of pg
∇ ⊥′ h · v on (ξ 1 0, ξ2j) for example,
can be calculated by applying the boundary condition (A.35), that is,
(A.43) ψ1,j =
φ1,j+1− φ1,j−1
2 ,
to get (pg ∇ ⊥′ h · v)0,j =(pg ∇ ⊥′ h · ∇hφ − pg ∇ ⊥′ h · ∇ ⊥ hψ)0,j =2φ0,j+1− φ0,j−1− ψ1,j+1− ψ1,j−1 h2 =2φ0,j+1− 2φ0,j−1− (φ1,j+2− φ1,j) − (φ1,j− φ1,j−2) h2 =−(φ1,j+2− 2φ0,j+1+ φ1,j) + (φ1,j− 2φ0,j−1+ φ1,j−2) h2 =(pg △ ′ hφ)0,j−1− (pg △ ′ hφ)0,j+1 2 = 1 2 (pg q)0,j−1− (pg q)0,j+1 . (A.44) The calculation of (pg ∇ ⊥′ h · v)N,j is similar. In summary, we have (A.45) (pg ∇ ⊥′ h · v)i,j= (pg △ h ψ)i,j 1 ≤ i ≤ N − 1; 1 2 (pg q)0,j−1− (pg q)0,j+1 i = 0; 1 2 (pg q)N,j+1− (pg q)N,j−1 i = N. It follows from (A.45), (A.36) and (A.20) that
(A.46) k∇ ⊥′ h · vk 2 ¯ Ω g ≤ β2kqk2Ω¯ g
where β2 is a constant that only depends on pg
. In view of (A.40), (A.41) and (A.46), the estimate (A.3) follows with γ−1=√β1+ β2+ 1. This completes the proof of Theorem 5.
A.1.1. Construction and Estimate ofψ. We now proceed with detailed construction and estimate for the potential ψ asserted in Step 2.
We first define theH 1 2
h norm for periodic grid functions which is essential to our analysis. Let f be a
periodic grid function on {ξ0, ξ1, · · · , ξN ≡ ξ0}, we can expand f(ξj) =PN −1n=0 fˆnEnj and define
(A.47) kfk2 H 1 2 h , N −1 X n=0 | ˆfn|2µ n and for φ ∈ L2( ¯Ω g , R)/R2, (A.48) k ˜Dτφk2 H 1 2 h(˚Γg) , k ˜D2φ(ξ11, ·)k2 H 1 2 h + k ˜D2φ(ξ1N −1, ·)k2 H 1 2 h , where (A.49) µ n, X m∈ Mn cos2(mπh) κ2mn −1 and (A.50) M n, { 1, 2, · · · , N } n = 0; { 0, 1, · · · , N − 1 } n = N/2; { 0, 1, · · · , N } otherwise. 23
is whereκmn6= 0.
We have the following discrete analogue of the trace inequality: Lemma A.2. Let φ ∈ L2( ¯Ωg, R)/R
2. Then (A.51) k ˜Dτφk H12 h(˚Γg) ≤ 4pg max k △ ′ hφkΩ¯g .
Proof. We first show that
(A.52) k ˜D2φ(ξ11, ·)k2 H 1 2 h ≤ 2pg max k △ ′ hφk 2 ¯ Ω g . Since (A.53) E n j+1− Enj−1 2h = ˜λ2nE n j, ˜λ2n , sin(2nπh) h . It follows that (A.54) D˜2φ(ξ11, ξ2j) = N −1 X n=0 XN m=0 √ −1˜λ2nφˆmnCm1 Enj. Hence k ˜D2φ(ξ11, ·)k2 H 1 2 h = N −1 X n=0 N X m=0 ˜ λ2nφˆmnCm1 2 µn = N −1 X n=0 X m∈ M ˜ λ2nφˆmnCm1κmn κmn 2 µn ≤ 2 N −1 X n=0 XN m=0 κ2mn˜λ22n| ˆφmn|2 X m∈ M cos2(mπh) κ2mn µn = 2 N −1 X n=0 N X m=0 κ2mnλ˜22n| ˆφmn|2 (A.55)
where we have used (A.49) in the last equality. Since ˜ λ22n= 4 h2sin 2(nπh) cos2(nπh) ≤ h42min{sin 2(nπh), cos2 (nπh)} ≤ 4 h2 cos2 mπh 2 sin 2(nπh) + sin2 mπh 2 cos 2(nπh)= κ2mn, (A.56)
it follows from (A.55) and (A.56) that
(A.57) k ˜D2φ(ξ11, ·)k2 H 1 2 h ≤ 2 N −1 X n=0 N X m=0 κ4mn| ˆφmn|2≤ 2pg max k △ ′ hφk 2 ¯ Ω g 24
The estimate for k ˜D2φ(ξN −11 , ·)k2
H12 h
is similar and the proof is complete. We proceed with the construction of the potential ψ(ξ1, ξ2). Let
ae(ξj2), ˜ D2φ1,j+ ˜D2φN −1,j 2 = N −1 X n=0 ˆ aenEnj, ao(ξj2), ˜ D2φ1,j− ˜D2φN −1,j 2 = N −1 X n=0 ˆ aonEnj, (A.58) and define (A.59) σ e n, N −1X m=1 m is even ˜ λ2 m κ4mn −1 , σ o n, N −1X m=1 m is odd ˜ λ2 m κ4mn −1 The proposed ψ ∈ L2 0( ¯Ωg , R) is given by (A.60) ψ(ξ1 i, ξj2) = N −1 X m=1 N −1 X n=0 ˆ ψmnSmi Enj where (A.61) ψˆmn, ˆ ae n ˜ λmσ e n √ 2κ 4 mn , if m is even; ˆ ao n ˜ λmσ o n √ 2κ 4 mn , if m is odd.
It is easy to verify that the constructed ψ satisfies (A.35):
−D− nψ(0, ξj2) =D1+ψ(0, ξj2) = ψ(ξ1 1, ξj2) − 0 h = N −1 X n=0 N −1X m=1 m is even ˆ aen ˜ λ2 mσ e n κ4mn + N −1 X m=1 m is odd ˆ aon ˜ λ2 mσ o n κ4mn Enj = N −1 X n=0 ˆ aenσ e n N −1 X m=1 m is even ˜ λ2 m κ4mn + ˆao nσ o n N −1 X m=1 m is odd ˜ λ2 m κ4mn Enj = N −1 X n=0 (ˆaen+ ˆaon)Enj = ˜D2φ1,j = − ˜Dτφ1,j. (A.62)
where we have used Sm 1 =
√
2h˜λmin the second equality above. Similarly,
(A.63) Dn−ψ(1, ξ2j) = D1−ψ(1, ξj2) = 0 − ψ(ξ1 N −1, ξj2) h = N −1 X n=0 (ˆaen− ˆaon)Enj = ˜D2φN −1,j= ˜DτφN −1,j.
In addition, the constructed potential ψ decays at designed rate in Fourier modes. This enables us to give an inverse trace estimate:
Lemma A.3. Let ψ ∈ L2 0( ¯Ωg
, R) be given by (A.60, A.61) and
(A.64) kD− nψk2 H12 h(Γg) , kD+ 1ψ(0, ·)k2 H12 h + kD− 1ψ(1, ·)k2 H12 h .
Then there is a constant C that only depends on pg
such that (A.65) k △ h ψk˚Ω g ≤ CkD−nψk H12 h(Γg) .
Proof. We first expand pg
△
h
ψ with respect to the basis (A.16) on ˚Ω
g : (A.66) pg △ h ψi,j= N −1 X n=0 N −1 X m=1 −κ 2 mnψˆmnSmi Enj, 1 ≤ i ≤ N − 1, 1 ≤ j ≤ N.
From (A.61), we have
(A.67) pg △ h ψi,j= − N −1 X n=0 N −1X m=1 m is even κ2mnˆaen ˜ λmσ e n √ 2κ 4 mn Smi + N −1 X m=1 m is odd κ2mnaˆon ˜ λmσ o n √ 2κ 4 mn Smi Enj,
In view of (A.59), we therefore have
k △ h ψk2˚Ω g ≤ 2pg1 min N −1 X n=0 N −1X m=1 m is even |ˆaen| ˜ λmσ e n κ2mn 2 + N −1 X m=1 m is odd |ˆaon| ˜ λmσ o n κ2mn 2 = 1 2pg min N −1 X n=0 |ˆaenσ e n|2 N −1 X m=1 m is even ˜ λ2 m κ4mn + |ˆao nσ o n|2 N −1 X m=1 m is odd ˜ λ2 m κ4mn = 1 2pg min N −1 X n=0 |ˆaen|2σ e n+ |ˆaon|2σ o n ≤ 2pgC∗ min N −1 X n=0 (|ˆaen|2+ |ˆaon|2)µ
n (see (A.69) below)
= C∗ 4pg min (kae+ aok2 H 1 2 h + kae− aok2 H 1 2 h ) = C∗ 4pg min kD− nψk2 H12 h(Γg) . (A.68)
Here in the second inequality, we have used the estimatesσ
e n≤ C∗µ
n andσ
o n≤ C∗µ
n which will be given in
Lemma A.4 below.
Lemma A.4. There is a constant C∗ such that
(A.69) σ e n≤ C∗µ n, σ o n≤ C∗µ n
uniformly for all 0 ≤ n ≤ N − 1 and h small enough.
Proof. We will only show thatσ e n≤ C∗µ n, or equivalently (A.70) X m∈M n cos2(mπh) κ2mn ≤ C∗ N −1 X m=1 m is even ˜ λ2 m κ4mn , uniformly in n.
The proof forσ
o n≤ C∗µ
n is similar. Denote by
(A.71) sm/2, sin(mπh/2), cm/2, cos(mπh/2).
It follows that LHS of (A.70) = X m∈ Mn c2 m κ2mn ≤ X m∈ Mn 1 κ2mn =h 2 4 X m∈ Mn 1 s2 m/2c2n+ s2nc2m/2 , (L)n, (A.72) RHS of (A.70) = N −1 X m=1 m is even ˜ λ2 m κ4mn =h 2 4 N −1 X m=1 m is even s2 m/2c2m/2 (s2 m/2c2n+ s2nc2m/2)2 , (R)n. (A.73)
It suffices to show that
(A.74) (L)n≤ C∗(R)n
for some constant C∗ independent of n and h. We show it separately for the following cases: Case 1: 1 ≤ n ≤ N/2 − 1 or N/2 + 1 ≤ n ≤ N − 1.
In this case,M
n= { 0, 1, · · · , N }. Since the expressions in (A.72) and (A.73) are symmetric with respect
to n = N/2, it suffices to consider the case of 1 ≤ n ≤ N/2 − 1. Denote by
(A.75) xm/2= mπh 2 , yn= nπh, and let fy(x), 1
sin2x cos2y + sin2y cos2x, x ∈ [0,
π 2], y ∈ [πh, π 2 − πh] (A.76) gy(x), sin2x cos2x
(sin2x cos2y + sin2y cos2x)2, x ∈ [0,
π
2], y ∈ [πh, π 2 − πh] (A.77)
We can rewrite (A.72) and (A.73) as
(L)n= h 2 4 X m∈ Mn 1 s2 m/2c2n+ s2nc2m/2 =h 2 4 N X m=0 fyn(xm/2), (A.78) (R)n= h 2 4 N −1 X m=1 m is even s2 m/2c2m/2 (s2 m/2c2n+ s2nc2m/2)2 = h 2 4 N −1 X m=1 m is even gyn(xm/2). (A.79) Since (A.80) d dx 1 fy(x) = sin(2x) cos(2y), 27
it follows that fy(·) is decreasing if y ∈ [0, π/4] and increasing if y ∈ [π/4, π/2]. In either case, the integral
test can be used to estimate the sum in (A.79) to get
(L)n= h2 4 N X m=0 fyn(xm/2) ≤ h2 4 fyn(x0) + h 4 2 π Z π/2 0 fyn(x) dx = h2 4 1 sin2yn + h 4 sin yncos yn , 1 ≤ n ≤ N 4 h2 4 fyn(xN) + h 4 2 π Z π/2 0 fyn(x) dx = h2 4 1 cos2y n + h 4 sin yncos yn , N 4 ≤ n ≤ N 2 − 1 (A.81)
where we have used the identity Z π/2 0 1 a2sin2x + b2cos2xdx = π 2|ab|. On the other hand,
(A.82) d
dx 1
gy(x)
= 2(tan x cos2y + cot x sin2y)(sec2x cos2y − csc2x sin2y),
it follows that gy(·) is increasing when 0 ≤ x ≤ y, deceasing when y ≤ x ≤ π/2 and therefore attains its
maximum 1
4 sin2y cos2y at x = y. Consequently,
(A.83) πh N −1 X m=1 m is even gyn(xm/2) + max x∈[0,π/2]gyn(x) ≥ Z π/2 0 gyn(x) dx = π
4 sin yncos yn(sin yn+ cos yn)2
where we have used the identity Z π/2
0
sin2x cos2x
(a2sin2x + b2cos2x)2dx =
π
4|ab|(|a| + |b|)2. Therefore from (A.79),
we have (R)n= h2 4 N −1 X m=1 m is even gyn(xm/2) ≥ h 4π Z π/2 0 gyn(x) dx − πh max x∈[0,π/2]gyn(x) = h
16 sin yncos yn(sin yn+ cos yn)2−
h2
16 sin2yncos2yn
. (A.84)
Combining (A.81) and (A.84), we obtain
(A.85) (L)n (R)n ≤
4 cos yn(sin yn+ cos yn)2(h cos yn+ sin yn)
sin yncos yn− h(sin yn+ cos yn)2
, 1 ≤ n ≤ N 4,
4 sin yn(sin yn+ cos yn)2(h sin yn+ cos yn)
sin yncos yn− h(sin yn+ cos yn)2 , N 4 ≤ n ≤ N 2 − 1. That is, (A.86) (L)n (R)n ≤
4(sin yn+ cos yn)2(1 + h cot yn)
1 − 2h(csc(2yn) + 1) , 1 ≤ n ≤
N 4, 4(sin yn+ cos yn)2(1 + h tan yn)
1 − 2h(csc(2yn) + 1) , N 4 ≤ n ≤ N 2 − 1. 28
In either case, we have (A.87) (L)n (R)n ≤ 8(1 + h cot y1) 1 − 2h(csc(2y1) + 1) , 1 ≤ n ≤ N4 8(1 + h tan yN 2−1) 1 − 2h(csc(2yN 2−1) + 1) , N 4 ≤ n ≤ N 2 − 1 = 81 + 1 π 1 − 1π + O(h).
This completes the proof of (A.74) for case 1. Case 2: n = 0.
Using the same argument as in the proof of case 1, we have
(L)0= h 2 4 N X m=1 1 s2 m/2 ≤h 2 4 1 s2 1/2 + N X m=2 1 s2 m/2 ≤h 2 4 csc 2πh 2 + h 2π Z π2 πh 2 csc2θ dθ =h 2 4 csc 2(πh 2 ) + h 2πcot( πh 2 ). (A.88) and (A.89) (R)0= h2 4 N −1 X m=1 m is even c2 m/2 s2 m/2 ≥4πh Z π2 πh cot2θ dθ = h 4π cot(πh) + πh −π2. Hence (A.90) (L)0 (R)0 ≤ h2 4 csc 2(πh 2 ) + h 2πcot( πh 2 ) h 4π cot(πh) + πh −π 2 = 8 + O(h).
This completes the proof of (A.74) for case 2. Case 3: n = N/2. In this case, (A.91) (L)N/2= h2 4 N −1 X m=0 1 c2 m/2 , (R)N/2= h2 4 N −1 X m=1 m is even s2 m/2 c2 m/2 The estimate (A.92) (L)N/2 (R)N/2 ≤ 8 + O(h)
follows from the same argument as in case 2.
In view of (A.87, A.90, A.92), the estimate (A.74) follows and the proof for Lemma A.4 is completed. Denote by
(A.93) kak0= ha, ai
1 2 0, ¯Ω g =h2 N X i=0 ′ N X j=1 a2ij 12 . 29
