# Finite element methods with matching and nonmatching meshes for Maxwell equations with disc

SIAM J. NUMER. ANAL. Vol. 37, No. 5, pp. 1542–1570

c 2000 Society for Industrial and Applied Mathematics

FINITE ELEMENT METHODS WITH MATCHING AND NONMATCHING MESHES FOR MAXWELL EQUATIONS WITH DISCONTINUOUS COEFFICIENTS?

ZHIMING CHEN? , QIANG DU? , AND JUN ZOU§ Abstract. We investigate the ?nite element methods for solving time-dependent Maxwell equations with discontinuous coe?cients in general three-dimensional Lipschitz polyhedral domains. Both matching and nonmatching ?nite element meshes on the interfaces are considered, and optimal error estimates for both cases are obtained. The analysis of the latter case is based on an abstract framework for nested saddle point problems, along with a characterization of the trace space for H (curl ; D), a new extension theorem for H (curl ; D) functions in any Lipschitz domain D, and a novel compactness argument for deriving discrete inf-sup conditions. Key words. ?nite element method, Maxwell equations, interface problems, nonmatching meshes, error estimates, saddle point formulation, extension theorem AMS subject classi?cations. 65M60, 65N30, 35Q60, 46E35 PII. S0036142998349977

1. Introduction. This paper is concerned with the following Maxwell equations in a dielectric medium: (1.1) (1.2) (1.3) (1.4) ε ?t E ? curl H =J in ? × (0, T ), ? ?t H + curl E =0 in ? × (0, T ), div(ε E) = ρ in ? × (0, T ), div(? H) =0 in ? × (0, T ).

Here ? ? R3 is a simply-connected Lipschitz polyhedral domain with connected boundary which is occupied by the dielectric material. E and H are the electric and magnetic ?elds, and J and ρ are the current density and charge density. We assume that the permeability parameter ? and the permittivity parameter ε of the medium are discontinuous across an interface Γ ? ?, where Γ is the boundary of a simply? 1 ? ? and ?2 = ? \ ? ? 1 . ?2 is also connected Lipschitz polyhedral domain ?1 with ? assumed to be simply-connected which, in turn, implies that Γ is connected. Without loss of generality we consider only the case with ε and ? being two piecewise constant functions in the domain ?, namely, ε= ε1 ε2 in in ?1 , ?2 , ?= ?1 ?2 in in ?1 , ?2 ,

? Received by the editors January 18, 1999; accepted for publication (in revised form) September 30, 1999; published electronically May 4, 2000. http://www.siam.org/journals/sinum/37-5/34997.html ? Institute of Mathematics, Academia Sinica, Beijing 100080, People’s Republic of China (zmchen@math03.math.ac.cn). The work of this author was partially supported by the China NSF under grant 19771080, China National Key Project “Large Scale Scienti?c and Engineering Computing.” ? Department of Mathematics, Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong, and the Department of Mathematics, Iowa State University, Ames, IA 50011 (madu@ust.hk). The work of this author was partially supported by the USNSF under grant DMS9796208 and by a grant from HKRGC. § Department of Mathematics, The Chinese University of Hong Kong, Shatin, N.T., Hong Kong (zou@math.cuhk.edu.hk). The work of this author was partially supported by Hong Kong RGC grants CUHK 338/96E and CUHK 4004/98P.

1542

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1543

and εi , ?i (i = 1, 2) are positive constants. It is known that the electric and magnetic ?elds E and H must satisfy the following jump conditions across the interface Γ : (1.5) (1.6) [E × n] = 0, [H × n] = JΓ , [ε E · n] = ρΓ , [? H · n] = 0 ,

where n is the unit outward normal to ? ?1 and ρΓ and JΓ are the surface charge and current density. Throughout the paper, the jump of any function A across the interface Γ is de?ned as [A] := A2 |Γ ? A1 |Γ with Ai = A|?i , i = 1, 2. We supplement Maxwell equations (1.1)–(1.6) with the initial conditions E(x, 0) = E0 (x), H(x, 0) = H0 (x), x∈?.

Instead of solving the fully coupled Maxwell system (1.1)–(1.6), we are interested in ?nding the electric and magnetic ?elds separately. To do so, we ?rst eliminate the magnetic ?eld H in (1.1)–(1.2) to obtain the electric ?eld equations, (1.7) (1.8) ε ?tt E + curl (??1 curl E) =?t J div (ε E) =ρ in in ? × (0, T ), ? × (0, T ),

with the following interface and boundary conditions: (1.9) (1.10) [E × n] = 0, [εE · n] = ρΓ , E × n = 0 on ? ? × (0, T ). [??1 curl E × n] = ??t JΓ ,

Similarly, we obtain the magnetic ?eld equations from (1.1)–(1.6), (1.11) (1.12) ? ?tt H + curl (ε?1 curl H) = ? curl (ε?1 J) in ? × (0, T ) , div (? H) = 0 in ? × (0, T ) ,

with the following interface and boundary conditions: (1.13) [H × n] = JΓ , [? H · n] = 0, [ε?1 curl H × n] = ?[ε?1 J × n] , (1.14) curl H × n = ?J × n on ? ? × (0, T ). Due to the practical interests, there has been a great deal of work on the numerical approximation to time-dependent Maxwell equations and also on the convergence analysis of numerical schemes for stationary Maxwell equations and related models; see, for example, [1, 3, 10, 13, 18, 20, 25] and the references therein. On the contrary, not too much work was available in the literature concerning the convergence analysis or error estimates for the fully discrete numerical methods for the time-dependent Maxwell systems. For some recent work in this aspect, we refer readers to [12, 19] and [24] for time-dependent Maxwell systems with continuous coe?cients. In this paper we will study a fully discrete ?nite element method for the timedependent Maxwell equations with discontinuous coe?cients. There have been numerous studies on the use of ?nite element methods for elliptic equations having discontinuous coe?cients (see, for example, [5, 11, 14, 7]). For the study presented

1544

ZHIMING CHEN, QIANG DU, AND JUN ZOU

here, we will concentrate on the electric ?eld equations (1.7)–(1.10). There is no essential di?culty for the extension of our numerical analysis here to the case for magnetic ?eld equation (1.11)–(1.14). We will discretize the system (1.7)–(1.10) using an implicit scheme in time and a N? ed? elec’s edge element scheme in space (cf. [21, 22]). In particular, we will investigate ?nite element methods with both matching and nonmatching grids on the interface Γ. As we will see later, the convergence analysis on the nonmatching case seems to be much more technical than the matching case. For example, a major obstacle is how to ?nd some appropriate ?nite element spaces which satisfy suitable inf-sup conditions. The paper is organized as follows. In section 2, we study the stationary model elliptic problem, (1.15) (1.16) curl (α curl A) + γβ A = f div (β A) = g in ?, in ?,

with the following interface conditions and boundary conditions: (1.17) (1.18) [A × n] = 0, A×n=0 [β A · n] = gΓ , on ? ?. [α curl A × n] = fΓ on Γ,

Here α and β are piecewise constants in ?; i.e., α = αi , β = βi in ?i for i = 1, 2, and αi , βi are positive constants. For technical reasons, γ is a chosen constant which may be zero in some cases (as in section 3) or strictly positive in some other cases (as in section 4). Such a lower-order term is required to ensure that the corresponding bilinear form induces a full H (curl ; ?)-norm in the latter cases. Finite element methods with matching and nonmatching meshes on the interface will be considered in sections 3 and 4. The results for the above stationary system (1.15)–(1.16) are needed in the convergence analysis of ?nite element methods for the time-dependent equations (1.7)–(1.10), which will be discussed in section 5. A number of interesting theoretical results concerning the weak formulation and the ?nite element approximation will be stated there while some technical proofs for an extension theorem and an abstract framework for the convergence analysis will be provided in an appendix at the end of the paper. A few conclusion remarks will be furnished in section 6. 2. Stationary model problem. For the convenience of presentation, we ?rst give some notation that will be used throughout the paper: H (curl ; ?) = {v ∈ L2 (?)3 ; curl v ∈ L2 (?)3 } , H α (curl ; ?) = {v ∈ H α (?)3 ; curl v ∈ H α (?)3 } (α>0), H0 (curl ; ?) = {v ∈ H (curl ; ?); v × n = 0 on ? ?} . The spaces H (curl ; ?) and H α (curl ; ?) are equipped with the norms

2 ||v||curl ,? = ||v||2 0,? + ||curl v||0,? 1 /2

, .

2 ||v||α,curl ,? = ||v||2 α,? + ||curl v||α,?

1 /2

Here and in what follows, · 0,? denotes the L2 (?)3 -norm (or the L2 (?)-norm for scalar functions) and · α,? and | · |α,? denote the norm and the seminorm of the Sobolev space H α (?)3 (or H α (?) for scalar functions). Similar de?nitions are adopted

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1545

for ?1 and ?2 . If no confusion is caused, we will often omit ?, ?1 , and ?2 in the subscripts of the norms. The constant C will always represent a generic constant independent of the time step τ and the mesh size h, unless otherwise speci?ed. We now discuss the variational formulation and the ?nite element approximation for the stationary model equations (1.15)–(1.18). As usual, we introduce a Lagrangian multiplier p to relax the divergence condition (1.16); hence, we consider the following system: (2.1) (2.2) (2.3) (2.4) (2.5) curl (α curl A) + γβ A ? β ?p = f in ?, div (β A) = g in ?, [A × n] = 0, [β A · n] = gΓ on Γ, [p] = 0, [α curl A × n] = fΓ on Γ, p = 0, A × n = 0 on ? ?.

Here γ is a constant which may be taken to be zero unless otherwise stated. Note that if the source terms f and g satisfy div f = γg, then we know the solution p of the above system is equal to zero by taking the divergence of (2.1) and using the boundary condition p = 0 on ? ?. This justi?es the introduction of the Lagrangian multiplier p. To establish an appropriate variational formulation for the system (2.1)–(2.5), we need a new trace space of H0 (curl ; ?) on the interface Γ. We next introduce this trace space and present some of its properties for the later use. We know that any v ∈ H0 (curl ; ?) has a tangential trace v × n in H ?1/2 (Γ)3 , de?ned by (2.6) or v × n, ?

Γ

v × n, ?

Γ

=

?1

v · curl ? dx ?

?1

curl v · ?dx ? ? ∈ H 1 (?1 )3

=

?2

curl v · ?dx ?

?2

v · curl ?dx ?? ∈ H 1 (?2 )3 ∩ H0 (curl ; ?2 ),

where ·, · Γ denotes the duality pairing between H ?1/2 (Γ)3 and H 1/2 (Γ)3 (or the duality pairing between H ?1/2 (Γ) and H 1/2 (Γ) for scalar functions) and H0 (curl ; ?2 ) = {v ∈ H (curl ; ?2 ) ; v × n = 0 on ? ?}. However, this characterization of the trace v × n is rather inconvenient in applications since this trace mapping from H0 (curl ; ?) to H ?1/2 (Γ)3 is not onto. To overcome this di?culty, we introduce the following trace space on Γ for functions in H0 (curl ; ?): T (Γ) = {s ∈ H ?1/2 (Γ)3 ; ? v ∈ H0 (curl ; ?) such that v × n = s on Γ}. It is not di?cult to see that T (Γ) is a Banach space under the norm: (2.7) s

T (Γ)

= inf { v

curl ,? ; v

∈ H0 (curl ; ?) and v × n = s on Γ}.

We note that for domains with smooth boundary, T (Γ) is well known, and it is often denoted by H ?1/2 (div , Γ) (cf. [8]). For domains with a Lipschitz boundary, we next provide a new characterization of this space.

1546

ZHIMING CHEN, QIANG DU, AND JUN ZOU

For convenience, we denote X1 = H (curl ; ?1 ) and X2 = H0 (curl ; ?2 ). Note that the left-hand side of (2.6) does not make sense if ? belongs to X1 , but the right-hand side of (2.6) is well-de?ned for all v and ? in X1 ; thus, for any s ∈ T (Γ), instead of (2.6) we can de?ne (2.8) s, ?

1 ,Γ

=

?1

v · curl ? dx ?

?1

curl v · ? dx

? ? ∈ X1 ,

where v ∈ H0 (curl , ?) such that v × n = s on Γ. Using the Green’s formula (2.6) and the density of H 1 (?1 ) in X1 , we know that (2.8) is independent of the choice of v ∈ H (curl ; ?) such that v × n = s on Γ. Thus (2.8) is well-de?ned. Similarly, we can de?ne (2.9) s, ?

2 ,Γ

=

?2

curl v · ?dx ?

?2

v · curl ?dx

?? ∈ X2 ,

where v ∈ H0 (curl ; ?) such that s = v × n on Γ. It is clear from the de?nitions (2.7)–(2.9) that (2.10) s, ?

1 ,Γ

? s, ?

2 ,Γ

≤ s

T (Γ)

?

X1 ×X2

? ? ∈ X1 × X2 ,

which implies that s, · 1,Γ ? s, · 2,Γ de?nes a continuous linear functional in X1 ×X2 . Hence by Riesz representation theorem, there exists a Q ∈ X1 × X2 satisfying (2.11) (curl Q · curl ? + Q · ?) dx = s, ?

?1 ∪?2 1 ,Γ

? s, ?

2 ,Γ

? ? ∈ X1 × X2

and the following identity holds: (2.12) Q

X1 ×X2

=

?∈X1 ×X2

sup

s, ?

1 ,Γ

? s, ?

X1 ×X2

2 ,Γ

?

.

1 (?i )3 Now applying the Green’s formula (2.6) to (2.11) with test functions ? ∈ H0 and ? ∈ H 1 (?i )3 , respectively, we obtain

curl curl Q + Q = 0 in ?1 ∪ ?2 , curl Q × n = 0 on ? ?, curl Q × n = s on Γ where the last two relations are understood in the sense of H ?1/2 (? ?)3 and H ?1/2 (Γ)3 , respectively. This yields curl Q ∈ H0 (curl ; ?) and (2.13) s

T (Γ)

≤ curl Q

curl ,?

= Q

X1 ×X2 .

Combining (2.10)–(2.13) we have proved the following lemma. Lemma 2.1. For any s ∈ T (Γ), we have the equality s

T (Γ)

=

?∈X1 ×X2

sup

s, ?

1 ,Γ

? s, ?

X1 ×X2

2 ,Γ

?

.

A direct consequence of this lemma is that T (Γ) is indeed a Hilbert space. In fact, if

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1547

si ∈ T (Γ) and Qi ∈ X1 × X2 (i = 1, 2) is the corresponding solution of (2.11), then we can de?ne an inner product ((·, ·))Γ as follows: ((s1 , s2 ))Γ =

?1 ∪?2 1 /2

(curl Q1 · curl Q2 + Q1 · Q2 )dx .

Clearly, s T (Γ) = ((s, s))Γ . In practice, Lemma 2.1 is rather inconvenient as it uses information from both ?1 and ?2 to de?ne the norm on T (Γ). To overcome the inconvenience, we note that (2.14) |s |1,Γ = sup ?∈X1 s, ? 1,Γ , ? X1 |s |2,Γ = sup ?∈X2 s, ? 2,Γ ? X2

are also norms of T (Γ). It is clear that |s |1,Γ ≤ s

T (Γ) ,

|s |2,Γ ≤ s

T (Γ)

?s ∈ T (Γ).

Both | · |1,Γ and | · |2,Γ are, in fact, equivalent norms on T (Γ) to · T (Γ) . To show this, we need the following extension theorem which may be of independent interest. Lemma 2.2 (extension theorem). Assume U is a bounded domain in R3 with a Lipschitz boundary ?U . Let U ?? D. Then there exists a bounded linear operator E : H (curl ; U ) → H (curl ; R3 ) such that E v = v on U and supp(E v) ? D ? v ∈ H (curl ; U ). Because the proof of Lemma 2.2 has not been seen in the literature and the full proof of the result is somewhat involved, we will describe it in more detail in the appendix, section A1. Lemma 2.3. The norms | · |1,Γ and | · |2,Γ are equivalent to · T (Γ) . Proof. We need only to prove that |s |1,Γ ≤ C |s |2,Γ ?s ∈ T (Γ). For any ?1 ∈ X1 , we use Lemma 2.2 to get a function ? ∈ H0 (curl ; ?) such that ? = ?1 on ?1 and ? Note that s, ? Thus, we have

1 ,Γ X2

≤ ?

curl ;?

≤ C ? ?1

X1 .

=

s, ?

2 ,Γ

for any ? ∈ H0 (curl ; ?) by the Green’s formula. s, ? 2,Γ = sup ?1 ∈X1 ?1 X1 s, ? 2,Γ s, ? 2,Γ = |s |2,Γ . ≤ C ? sup ? curl ;? ?∈X2 ? X2

|s |1,Γ = sup ?1 ∈X1 ≤ C?

s, ?1 ?1

1 ,Γ X1

sup ?∈H0 (curl ;?)

This completes the proof. 3. Finite element methods with a matching grid. In this section we discuss the ?nite element method based on a matching ?nite element grid on the interface Γ; i.e., the restrictions from two ?nite element triangulations in ?1 and ?2 match with each other on Γ. So both triangulations from ?1 and ?2 are combined into a

1548

ZHIMING CHEN, QIANG DU, AND JUN ZOU

standard triangulation of the whole domain ?. In this case, we set γ = 0 in the system 1 (?) (conforming across (2.1)–(2.5) and we require that A ∈ H (curl ; ?) and p ∈ H0 the interface). This leads to the following weak formulation for the system (2.1)–(2.5) with γ = 0. 1 Find (A, p) ∈ H0 (curl ; ?) × H0 (?) such that (3.1) (3.2) a(A, B) + b(B, p) =

?

f · Bdx ? fΓ , B g qdx + gΓ , q

? Γ

2 ,Γ

? B ∈ H0 (curl ; ?),

b(A, q ) =

2 3

? q ∈ H1 0 (?),

where we assume f ∈ L (?) , g ∈ L2 (?), fΓ ∈ T (Γ) and gΓ ∈ H ?1/2 (Γ). The two bilinear forms a(·, ·) and b(·, ·) are de?ned by a(B, D) =

?

α curl B · curl Ddx ? B, D ∈ X, β B · ?qdx

?

b(B, q ) = ?

? B ∈ X, q ∈ Q,

1 (?). with X = H0 (curl ; ?) and Q = H0 The following compactness result is due to Weber [28]. Lemma 3.1. Let β be a piecewise constant function in ? and {En } be a sequence in L2 (?)3 satisfying

En

0 ,?

≤ C,

curl En

0 ,?

≤ C,

div β En

0 ,?

≤ C,

En × n|? ? = 0

2

for some constant C ; then {En } has a convergent subsequence in L (?)3 . With the above lemma, we may get the following theorem. Theorem 3.2. There exists a unique solution (A, p) ∈ X × Q to (3.1)–(3.2). The proof is omitted as it following from a standard argument (see, for example, [21]) based on the Babuska–Brezzi theory. We next discuss the ?nite element discretization of (3.1)–(3.2). We note that for the matching grid case, studies of ?nite element discretization for the Maxwell systems have been given in [2, 20]. The main error estimates given in this section are, in fact, similar to the existing ones in the literature. We include the discussion here for the sake of completeness, and some of the technical results are also useful for the later discussion on nonmatching grids. Let T h be a shape regular triangulation of ? which matches with the interface Γ, i.e., K ∩Γ = ?

? ?

? K ∈ T h.

Here K is the interior of K . Let Sh be the standard continuous piecewise linear ed? elec H (curl , ?)-conforming edge element space ?nite element space and Vh be the N? de?ned by (3.3) Vh = {vh ∈ H (curl ; ?); vh = aK + bK × x on K ? K ∈ T h} ,

where aK and bK are two constant vectors. It is known (cf. [21]) that any function vh ∈ Vh is uniquely determined by the degrees of freedom in the set ME (v) of the moments on each element K ∈ T h , which is given by ME (v) =

e

v · τ ds ;

e is an edge of K

.

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1549

Here τ is the unit vector along the edge e. For any K ∈ T h and p > 2, we denote Yp (K ) = {v ∈ Lp (K )3 ; curl v ∈ Lp (K )3 ; and v × n ∈ Lp (?K )3 }. We know from [2, Lemma 4.7] that the integrals required in the de?nition of ME (v) make sense for any v ∈ Yp (K ). Thus for any v ∈ H s (?)3 with curl v ∈ Lp (?)3 , where s > 1/2 and p > 2, we can de?ne an interpolation πh v such that πh v ∈ Vh , and πh v has the same degrees of freedom as v on each K in T h . Following the same argument as the one used in [12], we can prove the following interpolation error estimate. Lemma 3.3. For 1/2 < s ≤ 1, we have v ? πh v

0,K

+ curl (v ? πh v)

0,K

≤ C hs K v

s,curl ,K

? v ∈ H s (curl ; K ),

where hK is the diameter of K ∈ T h . Furthermore, under the assumptions on the domain ? given in section 1, the interpolation operator πh has the following property [21]. 1 Lemma 3.4. Let v = ?p with p ∈ H0 (?). Then, if v is regular enough to ensure 1 the existence of πh v, we have πh v = ?ph for some ph ∈ Sh ∩ H0 (?). Now de?ne two ?nite element subspaces of Vh and Sh : Xh = Vh ∩ H0 (curl ; ?),

1 Qh = Sh ∩ H0 (?).

Then the ?nite element approximation of (3.1)–(3.2) can be formulated as follows. Find (Ah , ph ) ∈ Xh × Qh such that (3.4) (3.5) a(Ah , Bh ) + b(Bh , ph ) =

?

f · Bh dx ? fΓ , Bh g qh dx + gΓ , qh

? Γ

2 ,Γ

? Bh ∈ X h ,

b(Ah , qh ) =

? qh ∈ Qh .

The following theorem is the main result of this section. Theorem 3.5. Assume the solution (A, p) of problems (3.1)–(3.2) has the following regularity: for some 1/2 < s ≤ 1, A ∈ H s (curl ; ?i ), p ∈ H 1+s (?i ), i = 1, 2.

Then there exists a unique solution (Ah , ph ) ∈ Xh × Qh to the discrete problem (3.4)– (3.5), and we have the error estimate

2

(3.6)

A ? Ah

curl ,? + p ? ph

1 ,?

≤ C hs

i=1

{ A

s,curl ,?i

+ p

1+s,?i }.

To show the theorem, we need the following embedding result proved in [2] for any polyhedral domains. Lemma 3.6. Let XN (?) and XT (?) be the following subspaces of L2 (?)3 : XN (?) = {v ∈ L2 (?)3 ; curl v ∈ L2 (?)3 , div v ∈ L2 (?), v × n = 0 on ? ?}, XT (?) = {v ∈ L2 (?)3 ; curl v ∈ L2 (?)3 , div v ∈ L2 (?), v · n = 0 on ? ?}. If the domain ? is a Lipschitz polyhedral domain in R3 , then there exists a real number r > 1/2 such that XT (?) and XN (?) are continuously embedded into H r (?)3 .

1550

ZHIMING CHEN, QIANG DU, AND JUN ZOU

Proof of Theorem 3.5. Let Nh = {Bh ∈ Xh ; b(Bh , qh ) = 0 ? qh ∈ Qh }. Then Theorem 3.5 follows immediately from Lemma 3.3, from the coercivity, and from the inf-sup condition below by applying the standard Babuska–Brezzi theory (for example, cf. [6]): (3.7) (3.8) sup

Bh ∈Xh

a(Bh , Bh ) ≥ C0 Bh b(Bh , qh ) ≥ C1 qh Bh curl ,?

2 curl ; ? 1 ,?

? Bh ∈ N h ,

? qh ∈ Qh .

The inf-sup condition (3.8) is a consequence of the Poinc? are inequality by taking Bh = ?qh ∈ Xh for any ?xed qh ∈ Qh . The coercivity (3.7) was proved in [21] under the assumption that the ?nite element triangulation T h is quasi-uniform. We now 1 give a proof assuming only that T h is shape regular. For any Bh ∈ Nh , let ? ∈ H0 (?) be the solution of the following problem: ?? · ?vdx =

? ? 1 Bh · ?vdx ? v ∈ H0 (?) .

Set w = Bh ? ??; then it is obvious to see (3.9) curl w = curl Bh , div w = 0 in ?; w × n = 0 on ? ? .

Hence w ∈ XN (?) and so w ∈ H r (?) for some r > 1/2 by Lemma 3.6. Now since curl w = curl Bh ∈ L∞ (?)3 , we know that πh w is well-de?ned. Note that πh Bh = Bh ; we obtain by Lemma 3.4 that (3.10) Bh = πh w + ? ? h for some ?h ∈ Qh .

On the other hand, by using Lemma 3.3, (3.9), and the local inverse estimate, we have w ? πh w

0,K

≤ Chr K ( w r,K + curl w r,K ) = Chr K ( w r,K + curl Bh r,K ) ≤ C w r,K + C curl Bh 0,K ,

which yields, by using Lemma 3.6 and (3.9), πh w

0 ,?

≤C( w ≤C( w

0 ,?

+ w r,? + curl Bh 0,? ) 0,? + curl w 0,? + curl Bh

0 ,? )

0 ,? )

(3.11)

≤ C ( w 0,? + curl Bh ≤ C curl Bh 0,? .

Here we have used Lemma 3.1 to conclude that w

0 ,?

≤ C curl w

0 ,?

= C curl Bh

0 ,? .

Now the theorem follows by multiplying (3.10) by β Bh , then using (3.11) and the fact that Bh ∈ Nh . 4. Finite element method with a nonmatching grid. The ?nite element method discussed in section 3 for solving the system (2.1)–(2.5) is based on a matching ?nite element mesh on the interface Γ. This imposes a serious restriction, especially

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1551

in three dimensions, on the triangulations in ?1 and ?2 : both must match with each other on Γ. We are now going to relax this restriction and consider a nonmatching mesh on the interface that allows the two triangulations in ?1 and ?2 to be generated independently. This advantage, however, brings some di?culty to the convergence analysis since the resulting ?nite element spaces will be nonconforming for both the unknown A and Lagrangian multiplier p. Similarly to the treatment of the divergence condition, we will also deal with the constraints [A × n] = 0 and [p] = 0 on Γ by a Lagrangian multiplier approach. As we shall see, the introduction of these new multipliers leads to a nested saddle point problem, for which we need to ?rst generalize the standard saddle point theory (for example, cf. [6, 17]). 4.1. Abstract framework. Let X, Q, M be three real Hilbert spaces and a : X × X → R, b : X × Q → R, c : Q × M → R be three continuous bilinear functionals. Given f ∈ X ′ , g ∈ Q′ , and χ ∈ M ′ , we consider the following problem: Find (u, p, λ) ∈ X × Q × M such that (4.1) (4.2) (4.3) a(u, v ) + b(v, p) = f, v b(u, q ) + c(q, λ) = g, q c(p, ?) = χ, ? ? v ∈ X, ? q ∈ Q, ? ? ∈ M.

In the above we have used the same notation ·, · to denote all three duality pairings. This saddle point problem and its discretization were considered earlier in [29]. The results to be shown below are more general than the ones in [29] in the sense that our statement on the coercivity and the inf-sup condition for both continuous and discrete cases is more precise and much less restrictive. We will give all the proofs in the appendix, section A2. We ?rst introduce two subspaces N1 ? Q and N2 ? X as follows: (4.4) (4.5) N1 = {q ∈ Q; c(q, ?) = 0 ? ? ∈ M }, N2 = {v ∈ X ; b(v, q ) = 0 ? q ∈ N1 }.

Note that in the de?nition of N2 , the test function q is required only to be in N1 , which naturally provides a relaxed condition on the spaces than that given in [29]. We have the following result on the existence and uniqueness of the solution for the system (4.1)–(4.3). Lemma 4.1. Assume that a(·, ·) is N2 -coercive, i.e., (4.6) a(v, v ) ≥ a0 v

2 X

? v ∈ N2 ,

and the following inf-sup conditions hold: (4.7) (4.8) inf sup b(v, q ) ≥ b0 , v X q Q c(q, ?) ≥ c0 q Q ? M

q ∈N1 v ∈X

?∈M q ∈Q

inf sup

for some positive constants a0 , b0 , c0 . Then the problem (4.1)–(4.3) has a unique solution (u, p, λ) in X × Q × M . Next let Xh ? X, Qh ? Q, and Mh ? M be three ?nite dimensional subspaces. We introduce the corresponding approximation of (4.1)–(4.3) as follows.

1552

ZHIMING CHEN, QIANG DU, AND JUN ZOU

Find (uh , ph , λh ) ∈ Xh × Qh × Mh such that (4.9) (4.10) (4.11) a(uh , vh ) + b(vh , ph ) = f, vh b(uh , qh ) + c(qh , λh ) = g, qh c(ph , ?h ) = χ, ?h ? vh ∈ Xh , ? qh ∈ Qh , ? ?h ∈ Mh .

Corresponding to (4.4)–(4.5) we set (4.12) (4.13) N1h = {qh ∈ Qh ; c(qh , ?h ) = 0 N2h = {vh ∈ Xh ; b(vh , qh ) = 0 ? ?h ∈ Mh }, ? qh ∈ N1h }.

We have the following result concerning the discrete problem (4.9)–(4.11). Lemma 4.2. Assume the following conditions hold. (i) a(·, ·) is N2h -coercive; i.e., there exists a constant a? > 0 such that (4.14) a(vh , vh ) ≥ a? vh

2 X

? vh ∈ N2h .

(ii) There exists a constant b? > 0 such that (4.15)

qh ∈N1h vh ∈Xh

inf

sup

b(vh , qh ) ≥ b? . vh X qh Q

(iii) There exists a constant c? > 0 such that (4.16)

?h ∈Mh qh ∈Qh

inf

sup

c(qh , ?h ) ≥ c? . qh Q ?h M

Then the discrete problem (4.9)–(4.11) has a unique solution (uh , ph , λh ) in the space Xh × Qh × Mh with the following error estimate, u ? uh ≤C

X

+ p ? ph u ? vh

Q X

+ λ ? λh + inf

qh ∈Qh

M

vh ∈Xh

inf

p ? qh

Q

+ inf

?h ∈Mh

λ ? ?h

M

,

where the constant C depends only on a? , b? , c? , and on the operator norms a , b , and c of the bilinear functional a(·, ·), b(·, ·), and c(·, ·), respectively. 4.2. The variational formulation. In this subsection we introduce a new weak formulation for the system (2.1)–(2.5) with the constant coe?cient γ for the lowerorder term being strictly positive. Note that the results of section 3 for matching ?nite element grids are applicable to both stationary and time-dependent Maxwell equations even in the absence of the lower-order term. For nonmatching ?nite element grids, however, without this additional term, we have some technical di?culty in the veri?cation of the coercivity of the bilinear form corresponding to the term curl (α curl A). Though such veri?cation has been made in the matching grid case; see Lemma 4.7. Let us now introduce the following spaces: X1 = H (curl , ?1 ), Q1 = H 1 (?1 ), Then set X = X1 × X2 , Q = Q1 × Q2 × T (Γ), M = H ?1/2 (Γ) X2 = {v ∈ H (curl , ?2 ); v × n = 0 on ? ?}, Q2 = {v ∈ H 1 (?2 ); v = 0 on ? ?}.

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1553

and de?ne three bilinear forms a : X × X → R, as follows:

2

a(A, B) =

?

γβ A · B dx +

i=1 2 ?i

αi curl Ai · curl Bi dx ? A, B ∈ X, ? s, B1 ? B ∈ X, (q, s) ∈ Q,

b(B, (q, s)) = ?

i=1

βi ?qi · Bi dx + s, B2

?i Γ

2 ,Γ

1 ,Γ

c((q, s), ?) = q1 ? q2 , ?

? (q, s) ∈ Q, ? ∈ M.

Now applying the standard technique of integration by parts leads immediately to the following weak formulation of problem (2.1)–(2.5). Problem (P). Given (f , g ) ∈ L2 (?)4 , ?nd (A, (p, t), λ) ∈ X × Q × M such that

2

(4.17) (4.18) (4.19)

a(A, B) + b(B, (p, t)) =

i=1 2 ?i

f · Bi dx ? fΓ , B2 g qi dx + gΓ , q2

i=1 ?i

2 ,Γ

? B ∈ X,

b(A, (q, s)) + c((q, s), λ) = c((p, t), ?) = 0 ? ? ∈ M.

Γ

? (q, s) ∈ Q,

The above weak formulation is di?erent from the one used in section 3; it is more consistent with the ?nite element discretization on a nonmatching grid on the interface Γ. Note that the above system can also be derived based on an optimal control formulation of the interface problem (1.15)–(1.18). Similar approaches have been studied for the Poisson-type equations in [14]. The above formulation can be used as a basis for the further development of nonoverlapping domain decomposition methods. First, we have the following result on the existence and uniqueness of solutions to the Problem (P). Theorem 4.3. There exists a unique solution (A, (p, t), λ) ∈ X × Q × M to the 1 system (4.17)–(4.19). Moreover, (A, p) ∈ H0 (curl , ?) × H0 (?) satis?es (2.1)–(2.2) in the sense of distribution, and the Lagrangian multiplier (t, λ) ∈ T (Γ) × H ?1/2 (Γ) satis?es the following relations: t = α1 curl A1 × n = α2 curl A2 × n ? fΓ λ = β1 A1 · n = β2 A2 · n ? gΓ in T (Γ), in H ?1/2 (Γ).

Proof. We only prove that the system (4.17)–(4.19) has a unique solution; the rest of the theorem can be easily obtained by an appropriate application of the Green’s formula. For the proof, we need to verify the three conditions in Lemma 4.1. First, it is obvious that the bilinear form a : X × X → R is N2 -elliptic; i.e., (4.20) a(B, B) ≥ a0 B

2 X

? B ∈ N2

for some constant a0 > 0. Second, we show that there exists a constant b0 > 0 such that (4.21) sup

B∈X

b(B, (q, s)) ≥ b0 (q, s) B X

Q

? (q, s) ∈ N1 × T (Γ).

1 (?). For any given s ∈ T (Γ), let Q ∈ H (curl , ?1 ) be the We note that N1 = H0 solution of the following problem:

(4.22)

?1

(curl Q · curl B + Q · B)dx = s, B

1 ,Γ

? B ∈ H (curl ; ?1 ).

1554 It is easy to see that (4.23)

ZHIMING CHEN, QIANG DU, AND JUN ZOU

|s |1,Γ = Q

curl ,?1 .

1 (?) be the solution of the following elliptic interface problem: Now let u ∈ H0

(4.24)

?

β ?u · ?v dx =

?1

1 β1 Q · ?v dx ? v ∈ H0 (?).

1 (?) which satis?es It is clear that (4.24) has a unique solution u ∈ H0

(4.25)

?u

0 ,?

≤C Q

0 ,? 1 .

1 With these preparations, for any q ∈ H0 (?) and s ∈ T (Γ), we de?ne

? = B

??q1 + ?u1 ? Q ??q2 + ?u2

in ?1 , in ?2 ,

? ∈ X and where qi = q |?i , ui = u|?i , i = 1, 2. It is obvious that B ? B (4.26)

X

= ? ?q1 + ?u1 ? Q 0,?1 + curl Q ≤ C ( ?q 0,? + Q curl ,?1 )

0 ,? 1

+ ? ?q 2 + ? u 2

0 ,? 2

where we have used (4.25). Note that since q1 = q2 , u1 = u2 on Γ, we have ?(q1 ? q2 ) × n = 0, ?(u1 ? u2 ) × n = 0 on Γ.

Thus, by using (4.24) and (4.22), we are able to obtain ? , (q, s)) = ? b(B

?

β ?q · (??q + ?u)dx +

?1 1 ,Γ

β1 ?q1 · Q dx + s, Q

1 ,Γ

=

?

β |?q |2 dx + s, Q β |?q |2 dx + Q

?

=

2 curl ,?1 ,

which yields, together with (4.26), (4.23), and Lemma 2.3, ? , (q, s)) b(B ≥ C ( ?q ? X B

0 ,?

+ Q

curl ,?1 )

≥C( q

1 ,?

+ s

T (Γ) )

.

This completes the proof of (4.21). Finally, we verify that there exists a constant c0 > 0 such that (4.27) c((q, s), ?) ≥ c0 ? (q, s) Q (q,s)∈Q sup

M

? ? ∈ M.

This follows immediately from the trace theorem c((q, s), ?) q1 , ? Γ = ? ≥ sup (q, s) Q (q,s)∈Q q1 ∈H 1 (?1 ) q1 1,?1 sup

M.

From (4.20), (4.21), (4.27), and Lemma 4.1 we conclude that the system (4.17)–(4.19) has a unique solution (A, (p, t), λ) ∈ X × Q × M .

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1555

4.3. Finite element discretization. In this subsection we propose a ?nite element method for solving the problem (4.17)–(4.19), which allows a nonmatching ?nite element grid on the interface Γ. Let T h1 and T h2 be a shape regular triangulation of ?1 and ?2 , respectively. They induce naturally two ?nite element triangulations Γh1 and Γh2 on the interface Γ. Let Γh0 be an another shape regular triangulation over Γ. Note that Γhi , i = 0, 1, 2, are allowed to be di?erent from each other. However, we make the following reasonable assumption: (H1) Each triangle in Γh1 and Γh2 must be contained in some triangle of Γh0 . Let Xhi ? Xi , i = 1, 2, be the N? ed? elec edge element space de?ned over T hi (cf. (3.3)), let Qhi ? Qi , i = 1, 2, be the standard piecewise linear ?nite element space over T hi , let Mh0 be the standard piecewise constant ?nite element space over Γh0 . Moreover, we de?ne Th0 (Γ) = {sh ? T (Γ); sh = (ατ + β τ × x) × n on any τ ∈ Γh0 , ατ , β τ ∈ R3 }. Note that on planar portions of the interface, the space Th0 (Γ) coincides with the lowest order Raviart–Thomas element (see [26]) which has been used often for electromagnetic integral equation calculations (cf. [4]). Now set Xh = Xh1 × Xh2 ? X, Qh = Qh1 × Qh2 × Th0 (Γ) ? Q, Mh = Mh0 ? M.

We will assume the following two inf-sup conditions: (H2) There exists a constant β ? > 0 independent of h0 , h1 , h2 such that sup

Bhi ∈Xhi

sh , Bhi i,Γ ≥ β ? sh Bhi curl ,?i

T (Γ)

? sh ∈ Th0 (Γ) ,

i = 1 or 2.

(H3) There exists a constant γ ? > 0 independent of h0 , h1 , h2 such that sup

qhi ∈Qhi

qhi , ?h Γ ≥ γ ? ?h qhi 1,?i

?1/2,Γ

? ?h ∈ Mh0 ,

i = 1 or 2 .

The assumptions (H2)–(H3) indicate that the mesh Γh0 should be coarse enough compared with the meshes T h1 or T h2 in order to stabilize the e?ect of the introduced Lagrangian multipliers. We will verify these two assumptions in the next subsection using a general compactness argument. Now we are in a position to introduce the discrete version of Problem (P). Problem (Ph ). Find (Ah , (ph , th ), λh ) ∈ Xh × Qh × Mh such that

2

(4.28) a(Ah , Bh ) + b(Bh , (ph , th )) =

i=1 ?i 2

f · Bhi dx ? fΓ , Bh2 g qhi dx + gΓ , qh2

i=1 ?i

2 ,Γ

? Bh ∈ X h ,

(4.29) b(Ah , (qh , sh )) + c((qh , sh ), λh ) = (4.30) c((ph , th ), ?h ) = 0 ? ?h ∈ Mh .

Γ

? (qh , sh ) ∈ Qh ,

The following theorem summarizes our main results of this section.

1556

ZHIMING CHEN, QIANG DU, AND JUN ZOU

Theorem 4.4. Under the assumptions (H1)–(H3) the discrete Problem (Ph ) has a unique solution (Ah , (ph , th ), λh ) ∈ Xh × Qh × Mh . It satis?es the error estimate, for some generic constant C > 0,

2

( Ai ? Ahi

i=1

curl ,?i 2

+ pi ? phi

1 ,? i )

+ t ? th inf

T (Γ) 2

+ λ ? λh pi ? qhi

?1/2,Γ

≤C

Bh ∈Xh

inf

Ai ? Bhi

i=1

curl ,?i

+

(qh ,sh )∈Qh

1 ,? i

+ t ? sh

T (Γ)

i=1

+ inf

?h ∈Mh

λ ? ?h

?1/2,Γ

.

Moreover, if f ∈ H s (?1 )3 × H s (?2 )3 , g ∈ H s (?1 ) × H s (?2 ), gΓ ∈ H s?1/2 (Γ), fΓ = ψ ×n for some ψ ∈ H s (curl ; ?1 ) and the solution (A, p) of problem (4.17)–(4.19) has the regularity, A ∈ H s (curl , ?1 ) × H s (curl , ?2 ), p ∈ H 1+s (?1 ) × H 1+s (?2 ), where s ∈ (1/2, 1]; then we have the following error estimate:

2

{ Ai ? Ahi

i=1 2

curl ,?i

+ pi ? phi + pi

1 ,? i }

+ t ? th + f

T (Γ)

+ λ ? λh

?1/2,Γ

(4.31)

≤C

i=1

hs i ( Ai

s,curl ,?i

1+s,?i

s,?i )

+ Chs 0 λ

?1/2+s,Γ .

We are going to apply the abstract framework in Lemma 4.2 to prove Theorem 4.4. We ?rst note that N1h = {(qh , sh ) ∈ Qh ; c((qh , sh ), ?h ) = 0 ? ?h ∈ Mh } = {(qh , sh ) ∈ Qh1 × Qh2 × Th0 (Γ); qh1 ? qh2 , ?h Γ = 0 ??h ∈ Mh0 }. For those functions in N1h , we have the following equivalence. Lemma 4.5. For any qh ∈ Qh1 × Qh2 , the following two relations are equivalent: (4.32)

Γ

(qh1 ? qh2 ) ?h dσ = 0 ?(qh1 ? qh2 ) · sh dσ = 0

Γ

?h , ? ?h ∈ M 0 ? sh ∈ Th0 (Γ),

(4.33)

? h = {?h ∈ Mh ; ?h , 1 Γ = 0} . where M 0 0 Proof. We denote by T h0 any shape regular triangulation of ?1 whose restriction ed? elec H (curl , ?1 )-conforming edge on Γ coincides with Γh0 , and let Xh0 be the N? element space over T h0 . Then from the de?nition of Th0 (Γ) we know that Th0 (Γ) = {?h × n; ?h ∈ Xh0 }. Thus (4.33) is equivalent to ?(qh1 ? qh2 ) · (?h × n)dσ = 0 ? ?h ∈ Xh0 .

Γ

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1557

Then, by integration by parts, we obtain (cf. [15, Lemma 2.1]) ?(qh1 ? qh2 ) · (?h × n)dσ = ? (qh1 ? qh2 )(curl ?h · n)dσ ? ?h ∈ Xh0 .

Γ

Γ

Now set Nh0 = {curl ?h · n ; ?h ∈ Xh0 }; then for our purpose it su?ces to show ? h = Nh . It is clear that Nh ? M ? h . Let nf , nv , ne stand for the number of faces, M 0 0 0 0 vertices, edges of the triangulation Γh0 , respectively. Using the equivalence between a two-dimensional tangential vector ?eld (in this case the tangential components of the elements of Xh0 on the interface Γ), de?ned on the edges of the triangulation having zero circulation on each triangle, and that being induced by a scalar twodimensional potential ?eld de?ned on the vertices of the triangulation [23], we have ? h due to the dim Nh0 = ne ? (nv ? 1) which is equal to nf ? 1, the dimension of M 0 Euler formula. This completes the proof. Furthermore, we verify the following inf-sup condition. Lemma 4.6. Under the assumptions (H2)–(H3) we have (4.34) sup

Bh ∈Xh

b(Bh , (qh , sh )) ≥ b? (qh , sh ) Bh X

Qh

? (qh , sh ) ∈ N1h ,

where b? > 0 is a constant independent of h0 , h1 , h2 . Proof. Without loss of generality we assume that (H2) is valid for i = 1. Thus for ? h ∈ Xh such that any sh ∈ Th0 (Γ), there exists a Q 1 1 (4.35) β ? sh

T (Γ)

≤

?h Γ sh , Q 1 . ? Qh1 curl ,?1

Then we de?ne Qh1 ∈ Xh1 to be the solution of the following problem: (curl Qh1 · curl Bh1 + Qh1 · Bh1 )dx =

?1

sh , Bh1

1 ,Γ

? Bh1 ∈ Xh1 .

? h as test functions, respectively, and using (4.35), we Taking Bh = Qh1 and Bh = Q 1 obtain (4.36) β ? sh

T (Γ)

≤ Qh1

curl ;?1

= |s |1,Γ ≤ sh

T (Γ) .

Now we use Qh1 to de?ne (uh , θh ) ∈ (Qh1 × Qh2 ) × Mh to be the solution of the following discrete elliptic interface problem:

2

βi ?uhi · ?vhi dx + vh1 ? vh2 , θh

i=1 ?i

Γ

(4.37) (4.38) uh1 ? uh2 , ?h

Γ

=

?1

β1 Qh1 · ?vh1 dx ? vh ∈ Qh1 × Qh2 , ? ?h ∈ Mh .

=0

Under the assumption (H3), this problem has a unique solution (uh , θh ). Choosing vh to be uh in (4.37) and using (4.38) we obtain

2

(4.39)

i=1

?uhi

0 ,? i

≤ C Qh1

0 ,? 1 .

1558 Now set

ZHIMING CHEN, QIANG DU, AND JUN ZOU

?h = B

??qh1 + ?uh1 ? Qh1 ??qh2 + ?uh2

in ?1 , in ?2 .

? h ∈ Xh , and we have from (4.39) that It is clear that B ?h B (4.40)

X

= ?uh1 ? ?qh1 ? Qh1

2

0 ,?

+ curl Qh1

curl ,?1 .

0 ,? 1

+ ?uh2 ? ?qh2

0 ,? 2

≤C

i=1

?qhi

0 ,?

+ Qh1

But from Lemma 4.5 we know that ?(qh1 ? qh2 ), sh

Γ

= 0,

?(uh1 ? uh2 ) , sh

Γ

= 0 ? sh ∈ Th0 (Γ),

which yields, together with taking vh = qh in (4.37) and using (4.36),

2

? h , (qh , sh )) = b(B

i=1 ?i 2

βi |?qhi |2 dx + sh , Qh1 ?qhi

i=1 2 0 ,? i

1 ,Γ

≥ C

+ Qh1

2 curl ,?1

? (qh , sh ) ∈ N1h .

Hence using (4.36) again we easily derive ? h , (qh , sh )) b(B ≥C ?h X B

2

?qhi

i=1

0 ,? i

+ sh

T (Γ)

? (qh , sh ) ∈ N1h .

Finally note the fact that qh2 = 0 on ? ? and Γ qh1 dσ = Γ qh2 dσ due to (qh , sh ) ∈ N1h , and we have by means of the standard argument for Poincar? e inequality (cf. [9]) that

2 2

qhi

i=1

1 ,? i

≤C

i=1

?qhi

0 ,? i .

This completes the proof of Lemma 4.6. Proof of Theorem 4.4. We know from (H3) that c : Q × M → R satis?es the inf-sup condition sup

(qh ,sh )∈Qh

c((qh , sh ), ?h ) ≥ γ ? ?h (qh , sh ) Q

M.

Thus, with this, (4.34), and the obvious coercivity of a(·, ·) we are able to apply Lemma 4.2 to conclude that

2

( Ai ? Ahi

i=1

curl ,?i

+ pi ? phi

1 ,? i ) 2

+ t ? th inf

T (Γ)

+ λ ? λh

1 /2 ,Γ

≤C (4.41) +

Bh ∈Xh

inf

Ai ? Bhi t ? sh

curl ,?i

+

i=1

qhi ∈Qhi

pi ? qhi .

1 ,? i

sh ∈Th0 (Γ)

inf

T (Γ)

+ inf

?h ∈Mh

λ ? ?h

?1/2,Γ

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1559

Now using the standard interpolation estimates (cf. [9]) and Lemma 3.3, we get

Bhi ∈Xhi

inf

Ai ? Bhi

curl ,?i

+

qhi ∈Qhi

inf

pi ? qhi + pi

1 ,? i 1+s,?i

(4.42) (4.43)

?h ∈Mh

≤ Chs i inf λ ? ?h

?1/2,Γ

Ai

s,curl ,?i

,

≤ Chs 0 λ

?1/2+s,Γ .

Next we introduce a triangulation T h0 in ?1 whose restriction on Γ coincides with Γh0 and let Xh0 be the N? ed? elec H (curl , ?1 )-conforming edge element over the mesh T h0 . Then from the de?nition of Th0 (Γ) we know that Th0 (Γ) = {?h × n; ?h ∈ Xh0 } . Now using the fact that t = α1 curl A1 × n := Q × n, we can easily show that, by (2.1), Q ∈ H s (curl , ?1 ) and (4.44) Q

s,curl ,?1

≤ C ( A1

s,curl ,?1

+ p1

1+s,?1

+ f

s,?1 ).

Thus we obtain by Lemma 2.3 and Lemma 3.3 that (4.45)

sh ∈Th0 (Γ)

inf

t ? sh

T (Γ)

≤C

inf Q ? ?h ?h ∈Xh0

curl ,?1

≤ Chs 0 Q

s,curl ,?1 .

Now Theorem 4.4 follows from (4.41)–(4.45). To conclude this subsection, we give some remarks on the technical di?culty encountered if the lower-order term γβ A in (2.1) is not included. First we show the following. Lemma 4.7. There exists a constant a? (h) > 0 such that

2

(4.46)

i=1

curl Bhi

2 0 ,? i

≥ a? (h) Bh

2 0 ,?

? Bh ∈ N 2 h .

Proof. First we note that N 2h = Bh ∈ X h ; sh , Bh1

2 1 ,Γ

? sh , Bh2

2 ,Γ

=0

and

βi ?qhi · Bhi dx = 0 ? (qh , sh ) ∈ N1h .

i=1 ?i

To see (4.46) we need only to show that for any Bh ∈ N2h , curl Bhi = 0 in ?i , (i = 1, 2) implies Bh = 0. Since curl Bhi = 0 in ?i , there exists ?hi ∈ Qhi such that Bhi = ??hi in ?i . The choice of ?h1 is unique up to a constant, and we thus may let (4.47) Now since Bh ∈ N2h , we get (4.48) and

2

?h1 , 1

Γ

= ?h2 , 1

Γ

.

sh , ??h2

2 ,Γ

? sh , ??h1

1 ,Γ

= 0 ? sh ∈ Th0 (Γ)

(4.49)

i=1 ?i

βi ?qhi · ??hi dx = 0

1560

ZHIMING CHEN, QIANG DU, AND JUN ZOU

for any qh ∈ Qh such that (4.50)

Γ

(qh1 ? qh2 ) ?h ds = 0 ? ?h ∈ Mh .

Note that sh , ??hi i,Γ = sh , ??hi Γ since sh , ??hi ∈ L2 (Γ)3 , i = 1, 2. Thus by Lemma 4.5 we know that the above ?h ∈ Qh1 × Qh2 satis?es ? h. (?h1 ? ?h2 ) ?h ds = 0 ? ?h ∈ M

Γ

Coupled with (4.47), we have that ?h ∈ Qh1 × Qh2 satis?es (4.50). Thus we can take qhi = ?hi in (4.49) and obtain ??hi = 0, or ?hi = 0 (i = 1, 2) since we have (4.47) and ?h2 = 0 on ? ?. This completes the proof. The inequality (4.46) does not imply that the bilinear form corresponding to the term curl (αcurl A) ful?lls the N2h -coercivity condition uniformly in h, thus the abstract results in section 4.1 cannot be applied without the lower-order term. Whether the dependence of a? (h) in (4.46) on h can be removed remains an interesting open question. 4.4. Veri?cation of inf -sup conditions. In this subsection we show that the assumptions (H2)–(H3) are valid at least when the mesh size h1 or h2 is suitably small compared with h0 . Let us ?rst introduce a projection operator Rh1 from H (curl , ?1 ) to Xh1 . For any Q ∈ H (curl , ?1 ), Rh1 Q ∈ Xh1 is the unique solution of the equation (4.51)

?1

curl (Rh1 Q ? Q) · curl Bh1 + (Rh1 Q ? Q) · Bh1 dx = 0 ? Bh1 ∈ Xh1 .

It is obvious that Rh1 Q ? Q

curl ,?1

≤

Bh1 ∈Xh1

inf

Q ? Bh1

curl ,?1 ,

which, together with Lemma 3.3 and the standard density argument, yields (4.52) Rh1 Q ? Q

curl ,?1

→0

as

h1 → 0.

Similar to the proof of Lemma 2 in [27], we can show the following lemma which indicates that the convergence in (4.52) is uniform in any compact subset of H (curl , ?1 ). Lemma 4.8. Let W be a ?xed compact subset of H (curl , ?1 ). Then given any ?1 = h ? 1 (ε, W ) > 0 such that for any Q ∈ W and any 0 < h1 < ε > 0, there exists a h ? 1, h Rh1 Q ? Q

curl ,?1

≤ε

? Q ∈ W.

Now we can state the following theorem for (H2). Theorem 4.9. For any given triangulation Γh0 on the interface Γ, there exists ? ? a constant h? 1 = h1 (h0 ) such that for any h1 < h1 , we have (4.53) sup

Bh1 ∈Xh1

sh , Bh1 1,Γ ≥ β ? sh Bh1 curl ,?1

T (Γ)

? sh ∈ Th0 (Γ),

where β ? > 0 is a constant independent of h0 , h1 , h2 .

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1561

Proof. Let Dh0 = {sh ∈ Th0 (Γ) ; |s |1,Γ = 1} be the unit sphere in Th0 (Γ). It is clear by Lemma 2.3 that Dh0 is compact in T (Γ) since Th0 (Γ) is a ?nite dimensional subspace of T (Γ). For any sh ∈ Th0 (Γ), let Q = Q(sh ) ∈ H (curl , ?1 ) be the unique solution of the problem (4.54)

?1

(curl Q · curl ? + Q · ?)dx = sh , ?

1 ,Γ

? ? ∈ H (curl , ?1 ).

From the de?nition (2.14) we have (4.55) Q

curl ,?1

=

|s |1,Γ .

Let W := {Q = Q(sh ); sh ∈ Dh0 } be the subspace of H (curl , ?1 ). It is easy to check ? that W is compact in H (curl , ?1 ). Thus by Lemma 4.8, there exists an h? 1 = h1 (h0 ) ? such that for h1 < h1 , Q ? Rh1 Q

curl ,?1

≤ 1/2 ? Q ∈ W.

Thus, for any sh ∈ Th0 (Γ) satisfying |s |1,Γ = 1, we obtain sh , Rh1 Q 1,Γ = Rh1 Q Rh1 Q curl ,?1

curl ,?1

≥ Q

curl ,?1

? Q ? Rh1 Q

curl ,?1

≥

1 , 2

where we have used (4.51), (4.54), and (4.55). This completes the proof of Theorem 4.9 by Lemma 2.3. We remark that similar results are valid for the inf-sup condition (H3) using the same arguments as that in the proof of Theorem 4.9. Moreover, by switching the de?nition of the projection to that on the domain ?2 , we can easily restate the theorem with h1 being replaced by h2 , with the corresponding spaces and norms de?ned on ?2 . The requirement on either h1 or h2 to be suitable small compared with h0 in Theorem 4.9 has also been used in [14] to verify the discrete inf-sup conditions for the ?nite element approximations of an interface problem for the Poisson equations with nonmatching ?nite element meshes. In the case where we expect that the mesh in one of the subdomains is much ?ner than the other (say, h1 ? h2 ), the condition on h1 being suitably small compared with h0 would appear only to be a very mild restriction on the meshes. How this restriction a?ects the practical performance of the numerical method and whether this condition could be lifted will be issues to be further examined, both theoretically and through numerical testings, in the future. 5. The time-dependent Maxwell equations. Finally we turn our attention to the main aim of the paper, i.e., to investigate ?nite element methods for the time-dependent electric ?eld equations (1.7)–(1.10). Again we introduce a Lagrangian multiplier for the divergence constraint (1.8) and consider the problem, (5.1) (5.2) (5.3) (5.4) (5.5) ε ?tt E + curl (??1 curl E) ? ε ?p = div (ε E) = [E × n] = 0, [ε E · n] = [p] = 0, [??1 curl E × n] = p = 0, E × n = f in ? × (0, T ), ρ ρΓ fΓ 0 in on on on ? × (0, T ), Γ × (0, T ), Γ × (0, T ), ? ? × (0, T ),

1562

ZHIMING CHEN, QIANG DU, AND JUN ZOU

together with the initial conditions E(x, 0) = E0 (x) and Et (x, 0) = E1 (x) in ?,

where f = ?t J and fΓ = ?t JΓ , and E1 (x) = ε?1 (J(x, 0) + curl H0 (x)), which is obtained from (1.1) with t = 0. We are going to approximate the system (5.1)–(5.5) using an implicit ?nite difference scheme in time and the edge element method in space. We shall study only the case with a nonmatching mesh below; the treatment of the case with a matching grid is similar and in fact much simpler. All the notations used in this section are carried over from those in section 4, unless otherwise speci?ed. Let X, Q, and M be the Banach spaces de?ned in section 4.2; and then we introduce three bilinear forms a : X × X → R, b : X × Q → R, and c : Q × M → R as follows:

2

a(A, B) =

i=1 2 ?i

1 ?? i curl Ai · curl Bi dx ? A, B ∈ X,

b(B, (q, s)) =

i=1 ?i

?εi ?qi · Bi dx + s, B2

Γ

2 ,Γ

? s, B1

1 ,Γ

? B ∈ X, (q, s) ∈ Q,

c((q, s), ?) = q1 ? q2 , ?

? (q, s) ∈ Q, ? ∈ M.

Then the weak formulation of (5.1)–(5.5) reads as follows. Find (E, (p, t), λ) in the following spaces E ∈ H 2 (0, T ; L2 (?)3 ) ∩ L2 (0, T ; X ), (p, t) ∈ L2 (0, T ; Q), λ ∈ L2 (0, T ; M ) such that it satis?es the initial conditions (5.6) E(x, 0) = E0 (x), ?t E(x, 0) = E1 (x), x∈?

and the equations (5.7) (5.8) (5.9) (ε ?tt E, B) + a(E, B) + b(B, (p, t)) = (f , B) ? fΓ , B2 b(E, (q, s)) + c((q, s), λ) = (ρ, q ) + ρΓ , q2 Γ c((p, t), ?) = 0 ? ? ∈ M

2 ,Γ

? B ∈ X,

? (q, s) ∈ Q,

for almost everywhere (a.e.) t ∈ (0, T ). To see that the system (5.6)–(5.9) admits a unique solution (E, (p, t), λ), one can apply the standard argument of method of lines combining with the results in section 4.2; we omit the details here. We next introduce a fully discrete scheme for the system (5.6)–(5.9). Let us ?rst divide the time interval (0, T ) into M equally spaced subintervals using the nodal points 0 = t0 < t1 < · · · < tM = T,

2 2 3 with tn = nτ and τ = T /M . For a given sequence {un }M n=0 in L (?) or L (?) , we de?ne the ?rst and second order backward ?nite di?erences:

?τ un =

un ? un?1 , τ

2 n ?τ u =

?τ un ? ?τ un?1 . τ

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1563

For a continuous mapping u : [0, T ] → L2 (?) or L2 (?)3 , we de?ne 1 u (·) = u(·, nτ ) and u ? (·) = τ

n n tn

u(·, s)ds

tn?1

for 1 ≤ n ≤ M . Let Xh ? X, Qh ? Q, and Mh ? M be the ?nite element spaces de?ned in section 4.3; then the fully discrete ?nite element approximation to (5.6)–(5.9) can be formulated as follows. n n n For n = 0, 1, . . . , M, ?nd En h ∈ Xh , (ph , th ) ∈ Qh , and λh ∈ Mh such that (5.10) E0 h = πh E0 ,

?1 E0 h ? Eh = τ πh E1 ,

and for any Bh ∈ Xh , (qh , sh ) ∈ Qh and ?h ∈ Mh the following equations hold:

2 n n n ?n ?n Eh , Bh ) + a(En (5.11) (ε?τ h , Bh ) + b(Bh , (ph , th )) = (f , Bh ) ? fΓ , Bh2 2,Γ , n ρn , q h ) + ρ ?n (5.12) b(En Γ , qh2 Γ , h , (qh , sh )) + c((qh , sh ), λh ) = (?

(5.13)

n c((pn h , th ), ?h ) = 0.

From Theorem 4.4 we know that under the hypotheses (H1)–(H3), the system (5.11)– n n n 2 (5.13) has a unique solution {En h , (ph , th ), λh } at each time step n, with γ = ε/τ here. Moreover, we have the following main results of this section. Theorem 5.1. Assume that for some 1/2 < s ≤ 1, the solution of the continuous problem (5.7)–(5.9) has the regularity E ∈ H 2 (0, T ; H s (curl , ?1 ) × H s (curl , ?2 )), p ∈ H 1 (0, T ; H 1+s (?1 ) × H 1+s (?2 )). Moreover, let ρ ∈ H 1 (0, T ; H s (?1 ) × H s (?2 )), f ∈ H 1 (0, T ; H s (?1 )3 × H s (?2 )3 ), ρΓ ∈ H 1 (0, T ; H s?1/2 (Γ)) and fΓ = ? × n for some ? ∈ H 1 (0, T ; H s (curl , ?1 )). Then we have the following error estimates:

2 1≤n≤M 2 n En h ?E i=1 curl ,?i

E ∈ H 3 (0, T ; L2 (?)),

max

n ?τ En h ? Et

0 ,?

+

≤Cτ +

i=1

s Ci hs i + C0 h0 .

Proof. The proof is standard in the sense that we ?rst estimate the errors between the discrete time-dependent solution and the so-called elliptic projection of the exact solution. The desired error estimates will then follow from the triangle inequality and the error estimates obtained for the elliptic projection (namely, Theorem 4.4). We de?ne the projection operator Ph : X × Q × M → Xh × Qh × Mh to be (Ah , (ph , th ), λh ) = Ph (A, (p, t), λ), for any (A, (p, t), λ) in X × Q × M, with (Ah , (ph , th ), λh ) in Xh × Qh × Mh satisfying (Ah ? A, Bh ) + a(Ah ? A, Bh ) + b(Bh , (ph ? p, th ? t)) = 0 ? Bh ∈ Xh , b(Ah ? A, (qh , sh )) + c((qh , sh ), λh ? λ) = 0 ? (qh , sh ) ∈ Qh , c((ph ? p, th ? t), ?h ) = 0 ? ?h ∈ Mh .

1564

ZHIMING CHEN, QIANG DU, AND JUN ZOU

We remark that the term (Ah ? A, Bh ) above is not necessary for the case with matching ?nite element grids; see section 3. Now, if we let B = τ ?1 Bh ∈ Xh , (q, s) = τ ?1 (qh , sh ) ∈ Qh , and ? = τ ?1 ?h ∈ Mh in (5.7)–(5.9) and integrate over (tn?1 , tn ), we obtain

n ?n pn , ? tn )) = (? f n , Bh ) ? ? fΓ , Bh2 2,Γ , (5.14) (ε ?τ En t , Bh ) + a(E , Bh ) + b(Bh , (? n n n n ? ? ρ , qh ) + ρ ?Γ , qh2 Γ , (5.15) b(E , (qh , sh )) + c((qh , sh ), λ ) = (? n ?n (5.16) c((? p , t ), ?h ) = 0 .

Letting

n n n n n n ?n) ? n pn , ? tn ), λ , ζh , δh ) = (En (ηh h , (ph , th ), λh ) ? Ph (E , (?

and subtracting equations (5.14)–(5.16) from the equations (5.11)–(5.13), together with the de?nition of the projection operator Ph , we derive (5.17) (5.18) (5.19)

n n 2 2 n ?n ηh , Bh ) + a(ηh , Bh ) + b(Bh , ζh ) = ε (?τ En (ε ?τ t ? ?τ Ph E , Bh ) ?n ? E ? n , Bh ) ?Bh ∈ Xh , +(Ph E n n b(ηh , (qh , sh )) + c((qh , sh ), δh ) = 0 ?(qh , sh ) ∈ Qh , n , ?h ) = 0 ??h ∈ Mh . c(ζh

n in (5.17) and using (5.18)–(5.19), we then have Taking Bh = 2τ ?τ ηh n?1 2 n?1 n?1 n 2 n n ε ?τ ηh ? ε ?τ ηh + a(ηh , ηh ) ? a(ηh , ηh ) n n 2 n ? ?n ? E ? n , ?τ η n ). ≤ 2τ ε (?τ Et ? ?τ Ph E , ?τ ηh ) + 2τ (Ph E h

Now, using the discrete Gronwall’s inequality and the estimates we have obtained for n the elliptic projections Ph (cf. Theorem 4.4), we get the estimate on ?τ ηh 0,? and n 2 n curl ηh 0,?i . The remaining L -norm ηh 0,? follows from the identity

n n ηh 2 0 ,?

=

k=0

n n 0 n τ (?τ ηh , ηh ) + (ηh , ηh ).

This implies the ?nal estimate given in the theorem with the help of the triangle inequality. We omit the details. 6. Conclusion. In this paper, we have studied the ?nite element approximations to the stationary and time-dependent Maxwell equations in a polyhedral domain using matching and nonmatching grids. We focus on the particular case where the coe?cients are allowed to display discontinuous behavior as they vary from subdomain to subdomain. In many practical electromagnetic applications, such scenarios arise frequently due to the spatial inhomogeneities. Aside from the technical results we have proved in the paper, it is worthwhile to point out that the freedom in choosing nonmatching meshes for di?erent subdomains will be a nice feature when developing e?ective numerical methods to simulate the complicated spatial structures. The abstract framework for the nonmatching grids outlined here will also be useful to the development of domain decomposition methods for the resulting linear (or even nonlinear) algebraic systems. We will pursue this and other issues as well as actual numerical testings in the future. Appendix. The proofs of some technical results quoted earlier in the paper are provided in this appendix.

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1565

A1. Proof of Lemma 2.2. The argument is similar to the proof for extensions of W 1,p functions (1 ≤ p < ∞); see [16], for example. But our construction here contains essential di?erences. Step 1. Since ?U is Lipschitz, for any point x on ?U, there exist, upon rotating and relabeling the coordinate axes if necessary, a system of orthogonal coordinates (y1 , y2 , y3 ), a cube Cx containing x, Cx = Π3 i=1 (?ai , ai ), and a Lipschitz continuous function Φ : (?a1 , a1 ) × (?a2 , a2 ) → (?a3 , a3 ) such that U ∩ Cx = {y ∈ Cx ; y3 > Φ(y1 , y2 )}, ?U ∩ Cx = {y ∈ Cx ; y3 = Φ(y1 , y2 )}.

′ ′ = (?a1 , a1 ) × (?a2 , a2 ) × (?a3 /2, a3 /2). We write in be the reduced cube Cx Let Cx what follows ′ , U + = U ∩ Cx ′ ?. U ? = Cx ?U

It is clear that if we de?ne n(y ) = ?Φ ?Φ (y1 , y2 ), (y1 , y2 ), ?1 ?y1 ?y2 ?y = (y1 , y2 , y3 ) ∈ Cx ,

then n is normal to ?U for y ∈ ?U ∩ Cx . ? )3 and suppose for the moment that supp(v) ? Step 2. Let v = (v1 , v2 , v3 ) ∈ C 1 (U ′ ? . Note that if y = Qx + b is the coordinate transformation with Q ∈ R3×3 being Cx ∩ U orthogonal matrix, then by the same technique as used by N? ed? elec [21], we know that ? (y ) = Qv(QT y ? QT b) v is in H (curl , U + ) in the coordinates y if and only if v ∈ H (curl ; U + ) in the coordinates x. Thus, without lost of generality, we can take Q as the identity matrix and b = 0 in the following. Let e3 = (0, 0, 1) and z = y + 2(Φ(y1 , y2 ) ? y3 )e3 . Notice that ? ? , and we set ? + for y ∈ U we have z ∈ U ?+ v+ (y ) = v(y ) if y ∈ U ? ?. v? (y ) = v(z ) + 2v3 (z )n(y ) if y ∈ U

′ Note that, on ?U ∩ Cx , we obviously have z = y and

v+ (y ) × n = v? (y ) × n . Step 3. We have (6.1) v?

curl ;U ?

≤C v

curl ;U

.

To prove this, let {Φk } be a sequence of C ∞ functions such that (cf. [16, p. 136]) Φk ≥ Φ, sup DΦk

k L∞

< +∞, Φk → Φ, DΦk → DΦ uniformly a.e.,

and let z k = y + 2(Φk (y1 , y2 ) ? y3 )e3 with nk (y ) = ? Φk ? Φk (y1 , y2 ), (y1 , y2 ), ?1 ?y1 ?y2 ?y = (y1 , y2 , y3 ) ∈ Cx .

1566 ? ? , denote For y ∈ U

ZHIMING CHEN, QIANG DU, AND JUN ZOU

vk (y ) = v(z k ) + 2v3 (z k )nk (y ). ? ? , simple calculations yield Then, for y ∈ U ?y × v(z k ) = ?z × v(z k ) + 2nk (y ) × ?y × v3 (z k )nk (y ) = ? v(z k ) , ?z3 ?v3 k ?z v3 (z k ) + 2 (z )nk (y ) × nk (y ) ?z3

+v3 (z k )?y × nk (y ) = ?z v3 (z k ) × nk (y ) . Thus we have, as k → ∞, curl vk (y ) → curl z v(z ) + 2n(y ) × ?v3 ?v2 ?v3 ?v1 ? , ? , 0 (z ) ?z3 ?z1 ?z3 ?z2

∞ a.e. for y in U ? . Now we obtain, for any ? ∈ C0 (U ? )3 ,

v? · curl ? dy = lim

U?

k→∞

vk · curl ? dy

U?

= lim =

k→∞

curl vk · ? dy

U?

curl z v(z ) + 2n(y ) ×

U? L∞

?v1 ?v3 ?v2 ?v3 ? , ? , 0 (z ) · ? dy . ?z3 ?z1 ?z3 ?z2 n

L∞

Recall that DΦ

< +∞, we get v

?

< +∞, and thus

curl ;U

curl ;U ?

≤C v

by change of variable formula. Step 4. De?ne ? + ? v v? Ev = ? 0 Ev

? +, in U in U ? , ? + ∪ U ? ). in R3 ? (U

′ ′ ? D. Now it is easy and supp(E v) ? Cx Note that E v × n is continuous on ?U ∩ Cx 3 to see by using (6.1) that E v ∈ H (curl ; R ) and curl ;R3

≤C v

curl ;U .

′ ?. ∩U This completes the proof in the case that v is C 1 , with support in Cx 1 ? Step 5. The rest of the argument is standard. We ?rst assume v ∈ C (U ) but then drop the restriction on its support. By using the compactness of ?U and the partition of unity, we can show that there exists a E v ∈ H (curl ; R3 ) with support in D such that

Ev

curl ;R3

≤C v

curl ;U .

? ) due to the Finally, if v ∈ H (curl ; U ), we approximate v by functions vk ∈ C 1 (U 1 ? density of C (U ) in H (curl ; U ) (see, for example, [17]), and set E v = lim E vk .

k→∞

This concludes the proof of the extension theorem.

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1567

A2. Proofs of the abstract results in section 4.1. Here, the proofs for the abstract results stated in section 4.1 are provided. Let X, Q, and M be three real Hilbert spaces with norms · X , · Q , and · M , respectively, and let X ′ , Q′ , and M ′ be their corresponding dual spaces with the norms · X ′ , · Q′ , and · M ′ . We will use the same notation ·, · to denote the duality pairings between X and X ′ , Q and Q′ , M and M ′ . Suppose there are three given continuous bilinear forms a : X × X → R, b : X × Q → R, c : Q × M → R,

and their operator norms are denoted as a , b , and c . With the bilinear forms b(·, ·) and c(·, ·), we associate two linear operators B ∈ L(X, Q′ ) and C ∈ L(Q, M ′ ) with their dual operators B ′ ∈ L(Q, X ′ ) and C ′ ∈ L(M, Q′ ) de?ned as follows: B ′ q, v = q, Bv = b(v, q ) ? v ∈ X, q ∈ Q, C ′ ?, q = ?, Cq = c(q, ?) ? q ∈ Q, ? ∈ M. ? ∈ X ′; f ?, v = 0, v ∈ V }. Later we will use the notation V = ker(B ) and V 0 = {f To study the existence and uniqueness of the problem (4.1)–(4.3) we ?rst recall the classical results of the Babuska–Brezzi theory (cf., for example, [6], [17]). Lemma A.1. The following three properties are equivalent: (i) There exists a constant β > 0 such that (6.2) inf sup b(v, q ) ≥ β. v X q Q

q ∈Q v ∈X

(ii) The operator B ′ is an isomorphism from Q onto V 0 and B′q

X′

≥β q

Q

? q ∈ Q.

(iii) The operator B is an isomorphism from V ⊥ onto Q′ and Bv

Q′

≥β v

X

? v ∈ V ⊥.

Lemma A.2. Assume that the bilinear form a(·, ·) is V -elliptic, i.e., there exists a constant α > 0 such that a(v, v ) ≥ α v

2 X

? v ∈ V,

and the bilinear form b(·, ·) satis?es the inf-sup condition (6.2). Then there exists a unique solution (u, p) ∈ X × Q to the following problem: a(u, v ) + b(v, p) = f, v b(u, q ) = g, q ? v ∈ X, ? q ∈ Q.

We now come to the proofs of Lemmas 4.1–4.2 in section 4.1. Proof of Lemma 4.1. For any given χ ∈ M ′ , by Lemma A.1 and the inf-sup ⊥ condition (4.8) there exists a unique p⊥ ∈ N1 such that C p⊥ = χ. Now by Lemma A.2 and (4.6)–(4.7), we know that there is a unique (u, p0 ) ∈ N2 × N1 such that (6.3) (6.4) Au + B ′ p0 = f ? B ′ p⊥ in X ′ , ′ Bu = g in N1 .

1568

ZHIMING CHEN, QIANG DU, AND JUN ZOU

0 From (6.4) we know that g ? Bu ∈ N1 , where 0 ′ N1 = {? ∈ N1 ; l, q = 0 ? q ∈ N1 }.

Thus, again by (4.8) and Lemma A.1, there exists a unique λ ∈ M such that C ′ λ = g ? Bu in Q′ . This proves Lemma 4.1 by letting p = p0 + p⊥ ∈ Q. For the proof of Lemma 4.2 we need the following lemma, for which we introduce Qh (χ) = {qh ∈ Xh ; c(qh , ?h ) = χ, ?h ? ?h ∈ Mh }, Xh (g ) = {vh ∈ Xh ; b(vh , qh ) = g, qh ? qh ∈ N1h }. Clearly we have Qh (0) = N1h and Xh (0) = N2h . Lemma A.3. With the inf-sup conditions (4.15)–(4.16), we have the following estimates:

qh ∈Qh (χ)

inf

p ? qh u ? wh

Q

≤ ≤

wh ∈Xh (g )

inf

X

c c? b 1+ ? b 1+

qh ∈Qh

inf

p ? qh u ? vh

Q,

vh ∈Xh

inf

X

+

c b?

?h ∈Mh

inf

λ ? ?h

M.

The proof of Lemma A.3 follows as a minor modi?cation of the proof for the classical Babuska–Brezzi theory (cf., for example, [17, p. 114]). Proof of Lemma 4.2. We argue the proof in the following three steps. Step 1. Since uh ∈ Xh (g ), then for any wh ∈ Xh (g ), we have vh = uh ? wh ∈ N2h . Thus, using (4.2) and (4.10), we obtain a(vh , vh ) = a(uh ? wh , vh ) = a(uh ? u, vh ) + a(u ? wh , vh ) = b(vh , p ? ph ) + a(u ? wh , vh ). Now for any qh ∈ Qh (χ) we have ph ? qh ∈ Qh (0) = N1h , which yields b(vh , ph ? qh ) = 0 ? qh ∈ Qh (χ). Note that vh ∈ N2h ; hence, for any qh ∈ Qh (χ) and wh ∈ Xh (g ), a? vh

2 X

≤ a(vh , vh ) = b(vh , p ? qh ) + a(u ? wh , vh ) ? wh ∈ Xh (g ),

qh ∈ Qh (χ),

which, together with the triangle inequality and Lemma A.3, implies u ? uh

X

≤C

vh ∈Xh

inf

u ? vh

X

+ inf

qh ∈Qh

p ? qh

Q

+ inf

?h ∈Mh

λ ? ?h

M

.

Step 2. For any ?h ∈ Mh , it follows from (4.16), (4.2), and (4.10) that λh ? ?h

M

≤ =

c(qh , λh ? ?h ) 1 sup c? qh ∈Qh qh Q

1 b(u ? uh , qh ) + c(qh , λ ? ?h ) sup c? qh ∈Qh qh Q 1 ≤ ? ( b u ? uh X + c λ ? ?h M ), c

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1569

which, along with the result from Step 1, yields λ ? λh

M

≤C

vh ∈Xh

inf

u ? vh

X

+ inf

qh ∈Qh

p ? qh

Q

+ inf

?h ∈Mh

λ ? ?h

M

.

Step 3. For any qh ∈ Qh (χ), since ph ? qh ∈ N1h , by means of (4.15), (4.1), and (4.9) we have ph ? q h

Q

≤ =

b(vh , ph ? qh ) 1 sup b? vh ∈Xh vh X

a(u ? uh , vh ) + b(vh , p ? qh ) 1 sup b? vh ∈Xh vh X 1 ? qh ∈ Qh (χ), ≤ ? a u ? uh X + b p ? qh Q b

which, combined with Lemma A.3 and the result from Step 1, implies p ? ph

Q

≤C

vh ∈Xh

inf

u ? vh

X

+ inf

qh ∈Qh

p ? qh

Q

+ inf

?h ∈Mh

λ ? ?h

M

.

This completes the proof of Lemma 4.2. Acknowledgments. The authors wish to thank an anonymous referee and the associate editor for many constructive comments that improved the paper.

REFERENCES ¨cker, A ?nite element formulation [1] J. J. Ambrosiano, S. T. Brandon, and E. Sonnendru of the Darwin PIC model for use on unstructured grids, J. Comput. Phys., 121 (1995), pp. 281–297. [2] C. Amrouche, C. Bernardi, M. Dauge, and V. Girault, Vector potentials in threedimensional nonsmooth domains, Math. Methods Appl. Sci., 21 (1998), pp. 823–864. ?, P.-A. Raviart, and J. Segre ?, On a ?nite-element [3] F. Assous, P. Degond, E. Heintze method for solving the three-dimensional Maxwell equations, J. Comput. Phys., 109 (1993), pp. 222–237. [4] A. Bendali, Numerical analysis of the exterior boundary value problem for the time harmonic Maxwell equations by a boundary ?nite element method. I and II, Math. Comp., 43 (1984), pp. 29–86. [5] J. Bramble and T. King, A ?nite element method for interface problems in domains with smooth boundaries and interfaces, Adv. Comput. Math., 6 (1996), pp. 109–138. [6] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991. [7] Y. Cao and M. D. Gunzburger, Least-squares ?nite element approximations to solutions of interface problems, SIAM J. Numer. Anal., 35 (1998), pp. 393–405. [8] M. Cessenat, Mathematical Methods in Electromagnetism, World Scienti?c, River Edge, NJ, 1998. [9] P. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [10] P. Ciarlet, Jr. and J. Zou, Finite element convergence for the Darwin model to Maxwell’s equations, RAIRO Mod? el Math. Anal. Num? er., 31 (1997), pp. 213–250. [11] Z. Chen and J. Zou, Finite element methods and their convergence for elliptic and parabolic interface problems, Numer. Math., 79 (1998), pp. 175–202. [12] P. Ciarlet, Jr. and J. Zou, Fully discrete ?nite element approaches for time-dependent Maxwell equations, Numer. Math., 82 (1999), pp. 193–219. [13] P. Degond and P.-A. Raviart, An analysis of the Darwin model of approximation to Maxwell’s equations, Forum Math., 4 (1992), pp. 13–44. [14] Q. Du and M. Gunzburger, A gradient method approach to optimization based multidisciplinary simulations and nonoverlapping domain decomposition methods, SIAM J. Numer. Anal., 37 (2000), pp. 1513–1541.

1570

ZHIMING CHEN, QIANG DU, AND JUN ZOU

[15] F. Dubois, Discrete vector potential representation of a divergence-free vector ?eld in threedimensional domains: Numerical analysis of a model problem, SIAM J. Numer. Anal., 27 (1990), pp. 1103–1141. [16] L.C. Evans and R.F. Gariepy, Measure Theory and Fine Properties of Functions, CRC Press, Boca Raton, FL, 1992. [17] V. Girault and P.-A. Raviart, Finite Element Methods for Navier–Stokes Equations, Springer-Verlag, Berlin, 1986. [18] D. W. Hewett and C. Nielson, A multidimensional quasineutral plasma simulation model, J. Comput. Phys., 29 (1978), pp. 219–236. [19] C. Makridakis and P. Monk, Time-discrete ?nite element schemes for Maxwell’s equations, RAIRO Mod? el Math. Anal. Num? er., 29 (1995), pp. 171–197. [20] P. Monk, Analysis of a ?nite element method for Maxwell’s equations, SIAM J. Numer. Anal., 29 (1992), pp. 714–729. ? de ?lec, Mixed ?nite elements in R3 , Numer. Math., 35 (1980), pp 315–341. [21] J. Ne ?lec, A new family of mixed ?nite elements in R3 , Numer. Math., 50 (1986), pp. 57–81. ? de [22] J. Ne [23] R. A. Nicolaides, Direct discretization of planar div-curl problems, SIAM J. Numer. Anal., 29 (1992), pp. 32–56. [24] R. Nicolaides and D.-Q. Wang, Convergence analysis of a covolume scheme for Maxwell’s equations in three dimensions, Math. Comp., 67 (1998), pp. 947–963. [25] P.-A. Raviart, Finite Element Approximation of the Time-Dependent Maxwell Equations, Tech. report GdR SPARCH 6, Ecole Polytechnique, France, 1993. [26] P.-A. Raviart and J. Thomas, A mixed ?nite element method for 2nd order elliptic problems, in Mathematical Aspects of the Finite Element Method, A. Dold and B. Eckman, eds., Lecture Notes in Math. 606, Springer, New York, 1977. [27] A. Schatz and J. Wang, Some new error estimates for Ritz–Galerkin methods with minimal regularity assumptions, Math. Comp., 65 (1996), pp. 19–27. [28] C. Weber, A local compactness theorem for Maxwell’s equations, Math. Methods Appl. Sci., 2 (1980), pp. 12–25. [29] L.-A. Ying and S.N. Atluri, A hybrid ?nite element method for Stokes ?ow: Part II–Stability and convergence studies, Comput. Methods Appl. Mech. Engrg., 36 (1983), pp. 36–60.