Analysis of a time-domain finite element method for 3-D Maxwell’s equations in dispersive media

Comput. Methods Appl. Mech. Engrg. 195 (2006) 4220–4229 www.elsevier.com/locate/cma

Analysis of a time-domain finite element method for 3-D Maxwells equations in dispersive media Jichun Li a b

a,*

, Yitung Chen

b

Department of Mathematical Sciences, University of Nevada, Las Vegas 4505 Maryland PKWY, Box 454020, Las Vegas, NV 89154-4020, USA Department of Mechanical Engineering, University of Nevada, Las Vegas 4505 Maryland PKWY, Box 454020, Las Vegas, NV 89154-4020, USA Received 9 February 2005; received in revised form 18 July 2005; accepted 9 August 2005

Abstract We consider the time dependent Maxwells equations in dispersive media on a bounded domain in three-dimensional space. A fully discrete finite element scheme is developed to approximate the electric field equation derived from the Maxwells equations. Optimal energy-norm error estimates are proved for Ne´de´lec curl-conforming edge elements. This is the first finite element error analysis for Maxwells equations in dispersive media.  2005 Elsevier B.V. All rights reserved. MSC: 65N30; 35L15; 78-08. Keywords: Maxwells equations; Dispersive media; Finite element method

1. Introduction Recently there is a growing interest in finite element modeling and analysis of Maxwells equations (e.g. [6,10–13,17,18]). The readers can find more references in the books [3,15,19,22] and some recent conference proceedings [1,4,7]. However, most work is restricted to non-dispersive media such as air in the free space. Very little work is devoted to dispersive media using finite element method (FEM), though there is some work in finite-difference time-domain (FDTD) modeling of dispersive media started since 1990 [23, Chapter 9]. The applications of time-domain finite element method (TDFEM) [16] for the dispersive media were seen only very recently [14]. We want to remark that dispersive media are ubiquitous, for example human tissue, soil, snow, ice, plasma, optical fibers and radar-absorbing materials. In order to accurately perform wideband electromagnetic simulations, we have to consider the effect of medium dispersion in the modeling equations. To our best knowledge, there exists no work in the literature studying the error analysis of TDFEM for the Maxwells equations in dispersive media. This paper intends to make an initial effort in this direction. Therefore, for simplicity, we only consider cold plasma (such as ionosphere and magnetosphere) in this paper. Analysis of other dispersive media will be discussed in our forthcoming paper. The governing equations that describe electromagnetic wave propagation in isotropic nonmagnetized cold electron plasma are [9,5]

*

Corresponding author. Tel.: +1 702 895 0365; fax: +1 702 895 4343. E-mail address: [email protected] (J. Li).

0045-7825/$ - see front matter  2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2005.08.002

J. Li, Y. Chen / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4220–4229

oE ¼ r  H  J; ot oH ¼ r  E; l0 ot oJ þ mJ ¼ 0 x2p E; ot

4221

ð1Þ

0

ð2Þ ð3Þ

where E is the electric field, H is the magnetic field, 0 is the permittivity of free space, l0 is the permeability of free space, J is the polarization current density, xp is the plasma frequency, and m is the electron-neutral collision frequency. Note that m = 0 reduces to the collisionless case discussed by Young [24]. Solving (3) with the assumption that the initial electron velocity is zero leads to [5, Eq. (8)] Z t Z t 2 mt ms 2 e Eðx; sÞ ds ¼ 0 xp emðtsÞ Eðx; sÞ ds; ð4Þ Jðx; tÞ ¼ 0 xp e 0

0

where x 2 X. Here we let X be a bounded polyhedral domain in R3 with boundary oX and unit outward normal n. Instead of solving the coupled system (1)–(3) with both the electric and magnetic fields as unknowns, we eliminate H, by taking the time derivative of (1) and using (2)–(4), to obtain the second order electric field equation 2 2 0 E tt þ r  ðl1 0 r  EÞ þ 0 xp E  0 xp mwðEÞ ¼ 0;

where we denote Z t wðEÞ ¼ emðtsÞ Eðx; sÞ ds.

ð5Þ

ð6Þ

0

Moreover, we assume that the boundary of X is a perfect conductor so that nE ¼0

on oX  ð0; T Þ.

ð7Þ

In addition, we assume the same initial conditions as [18,6] Eðx; 0Þ ¼ E 0 ðxÞ

and

Hðx; 0Þ ¼ H 0 ðxÞ x 2 X;

where E0 and H0 are given functions. Using (1) and (4), we obtain the initial conditions for (5) as follows: Eðx; 0Þ ¼ E 0 ðxÞ

and

E t ðx; 0Þ ¼ E 1 ðxÞ;

ð8Þ

where E 1 ðxÞ ¼ 1 0 r  H 0 ðxÞ. The purpose of the paper is to analyze the convergence of a fully discrete finite element scheme for the electric field equation (5). We prove an optimal error estimate in H(curl; X) for Ne´de´lecs H(curl; X) conforming tetrahedral elements [20] under the assumption that the solution is smooth enough. Our proofs are based on the techniques developed by Ciarlet and Zou [6] for non-dispersive media. The rest of the paper is organized as follows. In Section 2, we present our fully discrete finite element scheme using the Ne´de´lec second type curl conforming elements on tetrahedra. Section 3 is devoted to the proofs of the error estimates for the scheme. 2. Fully discrete finite element schemes We start this section by introducing some notation to be used in this paper. We define H ðcurl; XÞ ¼ fv 2 ðL2 ðXÞÞ3 ; r  v 2 ðL2 ðXÞÞ3 g; H a ðcurl; XÞ ¼ fv 2 ðH a ðXÞÞ3 ; r  v 2 ðH a ðXÞÞ3 g; H 0 ðcurl; XÞ ¼ fv 2 H ðcurl; XÞ; n  v ¼ 0 on oXg; where a P 0 is a real number, and (Ha(X))3 is the standard Sobolev space [21] equipped with the norm k Æ ka and semi-norm j Æ ja. Specifically, k Æ k0 denotes the (L2(X))3-norm. Also, H(curl; X) and Ha(curl; X) are equipped with the norm 2

2 1=2

kvk0;curl ¼ ðkvk0 þ kcurl vk0 Þ kvka;curl ¼

ðkvk2a

þ

;

kcurl vk2a Þ1=2 .

In this paper, C denotes a generic constant which is independent of both the time step s and the finite element mesh size h.

4222

J. Li, Y. Chen / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4220–4229

Now let us define the second Ne´de´lec spaces based on tetrahedra [20] to be used in this paper. Hence we assume that Th is a shape regular triangulation of X with a mesh size h made of tetrahedra. Following [20,18], we denote Pk the space of polynomials of total degree at most k P 1, P~ k the subspace of Pk consisting of homogeneous polynomials of degree exactly k, and n o Dk ¼ ðP k1 Þ3  pðxÞxjp 2 Pe k1 . Let K 2 Th be a non-degenerated tetrahedron with face f and edge e. Let se be a vector in the direction of e and u 2 (W1,l(K))3 for some l > 2. Furthermore, we define the following three sets of degrees of freedom: Z  M e ðuÞ ¼ ðu  se Þq dsjq 2 P k ðeÞ on the six edges of K ; ð9Þ Ze  u  q dAjq 2 Dk1 ðf Þ tangent to the face f for the four faces f of K ; ð10Þ M f ðuÞ ¼ f Z  M K ðuÞ ¼ u  q dV jq 2 Dk2 ðKÞ . ð11Þ K

Ne´de´lec [20] proves that these degrees of freedom are (Pk)3 unisolvent and H(curl; X) conforming. We like to mention that in (9) we can select any arbitrary base in Pk(e) and define degrees of freedom by using that base. If we switch to some other base, the degrees of freedom will be different, but overall error estimate will be the same. The same comment applies for definition (10) and (11). Define n o 3 V h ¼ vh 2 H ðcurl; XÞjvh jK 2 ðP k Þ ; 8K 2 T h ; then any function in Vh can be uniquely defined by the degrees of freedom (9)–(11) on each K 2 Th. Hence for any u 2 (W1, l(X))3, l > 2, we can define the interpolation operator Phu 2 Vh such that PhujK 2 VhjK has the same moments (9)–(11) as u on K for each K 2 Th. Taking the boundary condition (7) into account, we define V 0h ¼ fvh 2 V h jn  vh ¼ 0 on oXg. The following interpolation properties for the Ne´de´lec spaces are given in [20] and re-stated clearly in [18, p. 718]. Lemma 2.1 (1) If u 2 (Hl+1(X))3, 1 6 l 6 k, then ku  Ph uk0 þ hku  Ph uk0;curl 6 Chlþ1 kuklþ1 . (2) If u 2 H0(curl; X) is such that u 2 (H1(X))3 and $ · u 2 (Hl(X))3, 1 6 l 6 k, then kr  u  r  Ph uk0 6 Chl kr  ukl . To define our fully discrete scheme, we divide the time interval (0, T) into M uniform subintervals by points 0 = t0 < t1 <    < tM = T, where tk = ks, and denote the kth subinterval by Ik = (tk1, tk]. Moreover, we define uk = u(Æ, ks) for 0 6 k 6 M, and denote the first and second order backward finite differences uk  uk1 os uk  os uk1 ; o2s uk ¼ . s s Now we can formulate our fully discrete finite element scheme for the electric field equations (5)–(8) as follows: os u k ¼

E 0h ¼ Ph E 0 ;

E 0h  E 1 h ¼ sPh E 1 ;

and for k = 1, 2, . . . , M, find 0 ðo2s E kh ; vÞ

þ l1 0 ðr 

E kh

2

E kh ; r

V 0h

ð12Þ

such that

 vÞ þ 0 x2p ðE kh ; vÞ  0 x2p mðwk1 8v 2 V 0h ; h ; vÞ ¼ 0

ð13Þ

where 1 wkh ¼ ems wk1 þ sðems E k1 þ E kh Þ 8k P 1. ð14Þ h h 2 Here E kh denotes the finite element solution of E at time t = tk. Obviously, for each k = 1, 2, . . . , M, the system (13) has a unique solution E kh by Lax–Milgram lemma [19, p. 20], since the coefficient matrix for E kh is symmetric positive definite. Notice that our scheme is constructed as an extension of [6, p. 197] to the dispersive media [14]. w0h ¼ 0;

J. Li, Y. Chen / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4220–4229

4223

3. Finite element error estimates 3.1. Preliminary estimates This section is devoted to proving two lemmas that will be used in the proof of the error estimates for the fully discrete finite element scheme (12)–(14). Since our proofs are based on the techniques of [6], we shall use many notations of [6] as well, in order for readers to follow our proofs easily. Following [6, p. 203], we need to assume that the electric field E is well defined in terms of time variable t on the interval [2s, T]. We also need the projection operator P h : H 0 ðcurl; XÞ ! V 0h defined by [19, p. 171] 8v 2 V 0h .

ðu  P h u; vÞ þ ðr  ðu  P h uÞ; r  vÞ ¼ 0

ð15Þ

Obviously, Ph is well-defined in H0(curl; X). Furthermore, for a > 1/2, Ph has the property ku  P h uk0;curl 6 ku  Ph uk0;curl

8u 2 H 0 ðcurl; XÞ \ H a ðcurl; XÞ.

Lemma 3.1 [6, p. 203]. For B = H1(curl; X) or B = (Ha(X))3 with a P 0, we have the following estimates: Z k 1 t 2 k 2 kos u kB 6 kut ðtÞkB dt 8u 2 H 1 ð0; T ; BÞ; s tk1 Z k 1 t 2 k 2 kutt ðtÞk2B dt 8u 2 H 2 ð0; T ; BÞ; kos u kB 6 s tk2 Z tk 2 2 k 2 k kuttt ðtÞkB dt 8u 2 H 3 ð0; T ; BÞ. kos ut  os u kB 6 Cs tk2

Lemma 3.2. For B = H1(curl; X) or B = (Ha(X))3 with a P 0, we have the following estimates:  2 Z Z  k 1  u   6s uðtÞ dt kut ðtÞk2B dt   s Ik Ik B

8u 2 H 1 ð0; T ; BÞ

and  2 Z Z  k1 1  u   uðtÞ dt 6 s kut ðtÞk2B dt  k s k I

8u 2 H 1 ð0; T ; BÞ.

I

B

Proof. Squaring both sides of the identity 1 u  s k

Z

1 uðtÞ dt ¼ s Ik

Z

Z Ik

!

tk

us ðsÞ ds dt

t

and using the Cauchy–Schwarz inequality, we have 0 2 1 2  2 Z  Z Z tk Z Z Z tk    k 1  1 1     2 u   6 @  A¼ uðtÞ dt 1 dt u ðsÞ ds dt u ðsÞ ds   dt  s s       k k k s Ik s2 s I I t I t ! Z k ! Z Z tk t 1 2 2 6 1 ds jus ðsÞj ds dt s Ik t t ! Z k ! Z Z tk Z t 1 2 2 2 1 ds jus ðsÞj ds dt ¼ s jus ðsÞj ds; 6 k s Ik k1 k1 t t I which completes the proof of the first part by integrating both sides in X. The proof of the second part follows exactly the same way except using the following identity: 1 uk1  s

Z

1 uðtÞ dt ¼ k s I

Z

Z Ik

t

tk1

! us ðsÞ ds dt.



ð16Þ

ð17Þ ð18Þ ð19Þ

4224

J. Li, Y. Chen / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4220–4229

Lemma 3.3. Let wk  w(Ek) defined by (6), wkh defined by (14). Then for any 1 6 n 6 M, we have " jwnh

n 2

 w j 6 CT s

n X

jE kh

k 2

E j þs

4

Z

#

tn 2

2

2

ðjEj þ jE t j þ jE tt j Þ dt .

0

k¼0

Proof. Letting t = tk in the definition of w (see (6)), we have k

k

w  wðE Þ ¼ e

mtk

Z

tk ms

e Eðx; sÞ ds ¼ e

mtk

"Z

0

which gives k

w ¼e

ms

w

k1

þe

mtk

tk1 ms

e Eðx; sÞ ds þ

Z

ms

e Eðx; sÞ ds ; tk1

0

Z

#

tk

ems Eðx; sÞ ds.

ð20Þ

Ik

Subtracting (20) from (14), we obtain Z  1 mtk k k k k1 k1 ms mtk mtk1 k1 ms ðe E h þ e E h Þ  e Eðx; sÞ ds wh  w ¼ e ðwh  w Þ þ e Ik 2 Z i 1 h mtk k k1 k k1 k1 mtk mtk1 e ¼ ems ðwk1  w Þ þ e ðE  E Þ þ e ðE  E Þ ds h h h Ik 2 Z  3 X 1 mtk k k k1 ðe E þ emt E k1 Þ  ems Eðx; sÞ ds ¼ þ emt ðERRÞi . Ik 2 i¼1

ð21Þ

Next we will estimate (ERR)i for i = 1, 2, 3 one by one. For (ERR)1, we have jðERRÞ1 j 6 ems jwk1  wk1 j 6 jwk1  wk1 j. h h

ð22Þ

For (ERR)2, we have s s  E k1 Þj 6 ðjE kh  E k j þ jE k1  E k1 jÞ. jðERRÞ2 j ¼ jðE kh  E k Þ þ ems ðE k1 h h 2 2

ð23Þ

To estimate (ERR)3, we need the following identity: Z



tk

tk1

Z k 1 1 t f ðsÞ  ðf ðtk1 Þ þ f ðtk ÞÞ ds ¼ fss ðsÞ  ðs  tk1 Þðs  tk Þ ds; 2 2 tk1

ð24Þ

which can be proved easily by integration by parts. Applying the Cauchy–Schwarz inequality to (24) with f(s) = emsE(x,s), we obtain 1 k jðERRÞ3 j 6 emt 2

Z

jfss ðsÞ  ðs  t

Ik

k1

s2 k Þðs  t Þj ds 6 emt 2 k

Z

jfss ðsÞj ds Ik

Z s2 mtk e ems jm2 Eðx; sÞ þ 2mE s ðx; sÞ þ E ss ðx; sÞj ds k 2 I Z 6 Cs2 ðjEðx; sÞj þ jE s ðx; sÞj þ jE ss ðx; sÞjÞ ds.

¼

ð25Þ

Ik

Combining (21)–(25), we have jwkh

k

w j

jwk1 h

w

k1

s j 6 ðjE kh  E k j þ jE k1  E k1 jÞ þ Cs2 h 2

Z

ðjEðx; sÞj þ jE s ðx; sÞj þ jE ss ðx; sÞjÞ ds.

ð26Þ

Ik

Summing both sides of (26) over k = 1,2, . . . , n, and using the fact w0h ¼ w0 ¼ 0 lead to jwnh  wn j 6

Z n n X sX k1 2 ðjE kh  E k j þ jE k1  E jÞ þ Cs ðjEðx; sÞj þ jE s ðx; sÞj þ jE ss ðx; sÞjÞ ds. h 2 k¼1 Ik k¼1

ð27Þ

J. Li, Y. Chen / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4220–4229

4225

Squaring both sides of (27) and using Cauchy–Schwarz inequality, we have Z tn n X 2 2 2 n n 2 k k 2 k1 k1 2 2 4 n jwh  w j 6 Cs n ðjE h  E j þ jE h  E j Þ þ Cs t ðjEðx; sÞj þ jE s ðx; sÞj þ jE ss ðx; sÞj Þ ds 0

k¼1

6 Cstn "

n X

2

jE kh  E k j þ Cs4 tn

Z

6 CT s

2

2

2

ðjEj þ jE s j þ jE ss j Þ ds

0

k¼0 n X

tn

jE kh

k 2

E j þs

4

Z

#

tn 2

2

2

ðjEj þ jE s j þ jE ss j Þ ds ;

ð28Þ

0

k¼0

where we use the fact that tn 6 T for any 1 6 n 6 M.

h

Our main result is the following optimal error estimate in the energy norm. Theorem 3.1. Let E and E nh be the solutions of the electric field equations (5)–(8) and the finite element scheme (12)–(14) at time t and tn respectively. Furthermore, we assume that 3

r  r  E t ðtÞ 2 ðL2 ðXÞÞ ;

3

EðtÞ and E t ðtÞ 2 ðH lþ1 ðXÞÞ ;

0 6 t 6 T; 1 6 l 6 k

and 3

3

E ttt ðtÞ 2 ðL2 ðXÞÞ ;

E tt ðtÞ 2 ðH lþ1 ðXÞÞ ; s 6 t 6 T ; 1 6 l 6 k.

Then there is a constant C = C(T, 0, l0, xp, m, E), independent of both the time step s and mesh size h, such that 2

2

max ðkos E nh  E nt k0 þ kE nh  E n k0;curl Þ 6 Cðs2 þ h2l Þ;

16n6M

where 1 6 l 6 k, and k is the degree of edge elements in the space Vh. Remark 3.1. Our error estimates are obtained provided that the continuous solution E is smooth enough. Rigorous regularity analysis for the solutions of Maxwells equations in dispersive media is very challenging and still an open problem, though some interesting work has been done for Maxwells equations in free space (see [2,8] and references therein). 3.2. Proofs of the error estimates To obtain the error estimate for E k  E kh , we multiply (5) by s1 v 2 V 0h and integrate the resultant over X in space and over Ik in time to obtain Z   Z   Z  1 1 1 0 ðos E kt ; vÞ þ r  E dt; r  v þ 0 x2p E dt; v  0 x2p m wðEÞdt; v ¼ 0 8v 2 V 0h . ð29Þ sl0 s Ik s Ik Ik Let gkh ¼ E kh  P h E k . Subtracting (29) from (13) and rearranging the terms using the following four identities and Galerkin orthogonal properties: ðiÞ 0 ðo2s E kh ; vÞ  0 ðos E kt ; vÞ ¼ 0 ðo2s ðE kh  P h E k Þ; vÞ þ 0 ðo2s P h E k  o2s E k ; vÞ þ 0 ðo2s E k  os E kt ; vÞ; Z  1 k ðiiÞ l1 ðr  E ; r  vÞ  r  Edt; r  v h 0 sl0 Ik   Z 1 k k k k k 1 1 1 ¼ l0 ðr  ðE h  P h E Þ; r  vÞ þ l0 ðr  ðP h E  E Þ; r  vÞ þ l0 r  E  r  E dt; r  v ; s Ik  Z   Z  1 1 ðiiiÞ 0 x2p ðE kh ; vÞ  0 x2p E dt; v ¼ 0 x2p ðE kh  P h E k ; vÞ þ 0 x2p ðP h E k  E k ; vÞ  0 x2p Edt  E k ; v ; s Ik s Ik  Z    Z 1 1 k1 k1 k1 2 2 2 ðivÞ 0 x2p mðwk1 ; vÞ   x m wðEÞdt; v ¼  x mðw  w ; vÞ þ  x m w  wðEÞ dt; v ; 0 0 0 p p p h h s Ik s Ik

4226

J. Li, Y. Chen / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4220–4229

we have 0 ðo2s gkh ; vÞ þ

1 ðr  gkh ; r  vÞ þ 0 x2p ðgkh ; vÞ l0

k k ¼ 0 ðos E kt  o2s E k ; vÞ þ 0 ðo2s E k  o2s P h E k ; vÞ þ l1 0 ðr  ðE  P h E Þ; r  vÞ     Z Z 1 1 k k k k 2 2 r  E  r  Edt; r  v þ  x ðE  P E ; vÞ   x E  E dt; v  l1 0 p h 0 p 0 s Ik s Ik   Z 1 þ 0 x2p mðwk1  wk1 ; vÞ þ 0 x2p m wk1  wðEÞ dt; v . h s Ik

As in [6], taking v ¼ sos gkh ¼ gkh  gk1 above and using the inequality h 1 aða  bÞ P ða2  b2 Þ 2 lead to 0 x2p 0 1 2 2 k 2 k1 2 ðkðos gkh k20  kos gk1 ðkgkh k20  kgk1 k Þ þ ðkr  g k  kr  g k Þ þ h h 0 h h k0 Þ 0 0 2l0 2 2 k k k 6 0 sðos E kt  o2s E k ; os gkh Þ þ 0 sðo2s E k  o2s P h E k ; os gkh Þ þ l1 0 sðr  ðE  P h E Þ; r  os gh Þ     Z Z 1 1 k k k k 1 k 2 k 2 k r  E dt; r  os gh þ 0 xp sðE  P h E ; os gh Þ  0 xp s E  E dt; os gh  l0 s r  E  s Ik s Ik   X Z 8 1 k1 k1 k1 2 k 2 k þ 0 xp msðwh  w ; os gh Þ þ 0 xp ms w  wðEÞ dt; os gh ¼ ðIÞi . s Ik i¼1

ð30Þ

Then we shall estimate (I)i one by one for i = 1, 2, . . . , 8. For (I)1, using the Cauchy–Schwarz inequality and (19), we obtain 1 1 1 2 2 2 ðIÞ1 6 s0 kos E kt  o2s E k k0 þ s0 kos gkh k0 6 s0 kos gkh k0 þ Cs2 0 2 2 2

Z

tk 2

tk2

kE ttt k0 dt.

ð31Þ

For (I)2, using the Cauchy–Schwarz inequality, (16) and Lemma 2.1, we have 1 1 1 2 2 2 2 ðIÞ2 6 s0 kos gkh k0 þ s0 ko2s E k  P h o2s E k k0 6 s0 kos gkh k0 þ Cs0 h2ðlþ1Þ ko2s E k klþ1 2 2 2 Z tk 1 2 2ðlþ1Þ k 2 6 s0 kos gh k0 þ C0 h kE tt klþ1 dt ðby (18)Þ. 2 k2 t

ð32Þ

Using the definition of Ph and the same technique as the one used for (I)2, we have the estimate for (I)3: 1 s 1 s 2ðlþ1Þ k 2 k k k k 2 ðIÞ3 ¼ l1 kE k  P h E k k20 6 s0 kos gkh k20 þ C h kE klþ1 . 0 sðE  P h E ; os gh Þ 6 s0 kos gh k0 þ 2 2 20 l0 2 0 l20

ð33Þ

To analyze (I)4, we use Greens formula and the boundary condition (7) to derive [6, p. 205] ðIÞ4 ¼ l1 0 ¼ l1 0

Z Ik

Z

r  ðE  E k Þ dt; r  os gkh rr

Z

Ik



t tk



¼ l1 0

E s ds dt; os gkh



Z Ik

r  r  ðE  E k Þ dt; os gkh

 pffiffiffiffiffiffi k ¼ s0 os gh ;

Z Z t  2  1 1  2  rr 6 s0 kos gkh k0 þ E ds dt s ; 2  k 2 2 l s k 0 0

I

t

2

1 pffiffiffiffiffiffi l0 s0

0

2

where we used the inequality ða; bÞ 6 12 ðjjajj0 þ jjbjj0 Þ in the last step.

Z Ik



rr

Z

t

  E s ds dt

tk

ð34Þ

J. Li, Y. Chen / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4220–4229

Rb Rb Rb Using the inequality j a f ðxÞgðxÞ dxj2 6 ð a jf ðxÞj2 dxÞð a jgðxÞj2 dxÞ twice, we have Z 2 Z t  2 Z Z t      rr    E ds dt ¼ 1  ðr  r  E Þ ds dt s s  k    I Ik tk tk 2 ! Z  Z Z t Z   2   6 1 dt  r  ðr  E s Þ ds dt ¼ s Ik

Ik

Z

Z

tk 2

! Z

Ik

3

t

Z

!

tk

I

ð35Þ

t

2

jr  ðr  E s Þj ds dt

1 ds

6s 6s

tk

Z t 2    r  ðr  E s Þ ds dt  k  k

4227

ð36Þ

t 2

jr  ðr  E s Þj ds;

ð37Þ

Ik

which gives Z Z t  2 Z   2 3  rr  E ds dt 6 s kr  ðr  E t Þk0 dt. s   Ik

tk

Ik

0

Hence 1 s2 2 ðIÞ4 6 s0 kos gkh k0 þ 2 20 l20

Z

2

Ik

kr  ðr  E t Þk0 dt.

Using Lemma 2.1, (16) and the Cauchy–Schwarz inequality, we obtain ðIÞ5 6

0 x2p s k 0 x2p s 0 x2p s 2 2 2 2 kE  P h E k k0 þ kos gkh k0 6 C0 x2p sh2ðlþ1Þ kE k klþ1 þ kos gkh k0 . 2 2 2

ð38Þ

Using Lemma 3.2 and the Cauchy–Schwarz inequality, we shall have  2 Z Z  0 x2p s  k 1 0 x2p s 1 1 2 2 k 2 2 2   E  kos gh k0 6 0 xp s E dt þ kE t k0 dt þ 0 x2p skos gkh k0 . ðIÞ6 6 k 2  2 s Ik 2 2 I 0 Using Lemma 3.3 and Cauchy–Schwarz inequality, we obtain " # Z T k X 1 1 2 2 2 j j 2 2 k 2 2 4 kE h  E k0 þ CT s ðkEk0 þ kE t k0 þ kE tt k0 Þ dt . ðIÞ7 6 0 xp mskos gh k0 þ 0 xp ms CT s 2 2 0 j¼0 Moreover, by (16) and Lemma 2.1, we have 2

2

2

2

2

kE jh  E j k0 6 2ðkE jh  P h E j k0 þ kP h E j  E j k0 Þ 6 2½kgjh k0 þ Ch2ðlþ1Þ kE j klþ1 . Therefore, we have

" k X 1 2 k 2 2 2 ðIÞ7 6 0 xp mskos gh k0 þ C0 xp mT s kg0h k20 þ h2ðlþ1Þ kE 0 k2lþ1 þ ðkgjh k20 þ h2ðlþ1Þ kE j k2lþ1 Þ 2 j¼1 # Z T 2 2 2 þ s3 ðkEk0 þ kE t k0 þ kE tt k0 Þ dt . 0

By Lemma 3.2 and the Cauchy–Schwarz inequality, we easily come to Z 1 1 ðIÞ8 6 0 x2p ms2 kwt ðEÞk20 dt þ 0 x2p mskos gkh k20 . 2 2 Ik By the definition of w (see (6)), we easily obtain Z t wt ¼ EðtÞ  m emðtsÞ Eðx; sÞ ds; 0

which leads to

Z t 2 Z t Z t    2 2 2 2 2 jwt j 6 2jEðx; tÞj þ 2m2  emðtsÞ Eðx; sÞ ds 6 2jEðx; tÞj þ 2m2 jemðtsÞ j ds jEðx; sÞj ds 0

2

¼ 2jEðx; tÞj þ mð1  e2mt Þ

Z

0

t 2

2

jEðx; sÞj ds 6 2jEðx; tÞj þ m

0

Z

0

t 2

jEðx; sÞj ds. 0

ð39Þ

ð40Þ

4228

J. Li, Y. Chen / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4220–4229

Therefore we have kwt k20

6

2kEk20

þm

Z

t

kEk20 ds. 0

Following [6, p. 206], now summing both sides of (30) over k = 1, 2, . . . , n and using those estimates obtained above for (I)i, 1 6 i 6 8, we have 0 x2p n 2 0 1 kos gnh k20 þ kgh k0 kr  gnh k20 þ 2l0 2 2 0 x2p 0 2 0 1 2 2 kgh k0 6 kos g0h k0 þ kr  g0h k0 þ 2l0 2 2 Z Z n T X 2 2ðlþ1Þ k 2 2 þ Cs kos gh k0 þ Cs 0 kE ttt k0 dt þ C0 h k¼1 2ðlþ1Þ

þC þ

sh 0 l20

s

n X

2

2

kE k klþ1 þ

k¼1

C0 x2p sh2ðlþ1Þ

n X

s 20 l20

2 kE k klþ1

k¼1

T

s

Z

2

kE tt klþ1 dt

T 2

kr  r  E t k0 dt

0

1 þ 0 x2p s2 2

Z

T 2

kE t k0 dt þ ðIIÞ7 þ ðIIÞ8 ;

ð41Þ

0

where the constant C currently depends on 0, l0, xp, and m. Now we shall estimate the terms (II)7 and (II)8. For (II)7, we have # " Z T n k X X 2 2 2 j 2 2ðlþ1Þ 0 2 2ðlþ1Þ j 2 2 0 2 3 ðIIÞ7 ¼ Cs kgh k0 þ h kE klþ1 þ ðkgh k0 þ h kE klþ1 Þþ s ðkEk0 þ kE t k0 þ kE tt k0 Þ dt k¼1

"

6 CsT

0

j¼1 2 kg0h k0

þ

2 h2ðlþ1Þ kE 0 klþ1

þ

n X

2 ðkgjh k0

þ

2 h2ðlþ1Þ kE j klþ1 Þ

þ s3

Z

where we use the fact ns = T. (II)8 can be estimated as follows:  Z t Z n Z  n Z  X X 2 2 2 2 2 ðIIÞ8 ¼ Cs kEk0 þ kEk0 ds dt 6 Cs kEk0 þ " ¼ Cs

2

Ik

n Z X k¼1

Ik

0 2 kEk0 dt

þ ns

Z

#

T 2 kEk0 ds

#

T 2 ðkEk0

þ

2 kE t k0

þ

2 kE tt k0 Þ dt

;

ð43Þ

0

j¼1

k¼1

ð42Þ

k¼1

Ik

¼ Cs2 ð1 þ T Þ

0

T 2 kEk0 ds

 dt

0

Z

T 2

kEk0 dt.

ð44Þ

0

Putting (43) and (44) into (41) and using s

n X

kE k k2lþ1 6 T max kEðtÞk2lþ1 06t6T

k¼1

in (41) lead to 0 x2p n 2 0 1 kos gnh k20 þ kgh k0 kr  gnh k20 þ 2l0 2 2 Z T n X 2ðlþ1Þ 2 4 6 CQ0 þ Q1 ðEÞ  ðs þ h Þ þ Cs ðkEk20 þ kE t k20 þ kE tt k20 Þ dt þ Cs ðkos gkh k20 þ kr  gkh k20 þ kgkh k20 Þ; 0

ð45Þ

k¼1

where 2

2

2

Q0 ¼ kos g0h k0 þ kr  g0h k0 þ kg0h k0 ; Z T Z T 2 2 2 2 2 2 Q1 ðEÞ ¼ ðkE ttt k0 þ kE tt klþ1 Þ dt þ ðkr  r  E t k0 þ kE t k0 þ kEk0 Þ dt þ max kEðtÞklþ1 . s

0

06t6T

From now on, the rest of the proofs are simple modifications of [6]. First let us consider the initial error. Using the definition of projection Ph, and E 0h ¼ Ph E 0 , we have kr  g0h k0 ¼ kr  ðE 0h  P h E 0 Þk0 ¼ kr  ðPh E 0  P h E 0 Þk0 ¼ kr  P h ðPh E 0  E 0 Þk0 6 kPh E 0  E 0 k0;curl 6 Chl kE 0 klþ1 .

ð46Þ

J. Li, Y. Chen / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4220–4229

4229

As for the estimate of g0h , we have kg0h k0 ¼ kE 0h  P h E 0 k0 ¼ kP h ðPh E 0  E 0 Þk0 6 kPh E 0  E 0 k0;curl 6 Chl kE 0 klþ1 .

ð47Þ

The estimate of os g0h follows easily from [6, p. 206] kos g0h k0 6 s1 kEðsÞ  Eð0Þ þ sE t ð0Þk0;curl þ kPh E t ð0Þ  E t ð0Þk0;curl 6 Cs sup kE tt k1 þ Chl kE t ð0Þklþ1 .

ð48Þ

ðs;0Þ

Therefore 2

2

2

Q0 6 Cs2 sup kE tt k1 þ Ch2l ðkE 0 klþ1 þ kE t ð0Þklþ1 Þ.

ð49Þ

ðs;0Þ

Substituting (49) into (45) and applying the discrete Gronwalls inequality, we obtain 2

2

2

max ðkos gnh k0 þ kr  gnh k0 þ kgnh k0 Þ 6 Cðs2 þ h2l Þ.

16n6M

ð50Þ

Finally, using the triangle inequality to os E nh  E nt ¼ ðos E nh  P h os E n Þ þ ðP h os E n  os E n Þ þ ðos E n  E nt Þ and E nh  E n ¼ ðE nh  P h E n Þ þ ðP h E n  E n Þ; we have completed the proof of the error estimate 2

2

max ðkos E nh  E nt k0 þ kE nh  E n k0;curl Þ 6 Cðs2 þ h2l Þ.

16n6M

Acknowledgements The authors wish to thank the three anonymous referees for many constructive comments for improving the paper. The author JL wants to thank Dr. Miguel Visbal at Air Vehicles Directorate of US Air Force Research Laboratory (AFRL) for pointing out this research topic to him. The work was partially supported by the Summer Faculty Program of AFRL. References [1] M. Ainsworth, P. Davies, D. Duncan, P. Martin, B. Rynne (Eds.), Topics in Computational Wave Propagation, Lecture Notes in Computational Science and Engineering 31, Springer-Verlag, Berlin, 2003. [2] F. Assous, P. Ciarlet Jr., P.-A. Raviart, E. Sonnendrcker, Characterization of the singular part of the solution of Maxwells equations in a polyhedral domain, Math. Method. Appl. Sci. 22 (1999) 485–499. [3] A. Bossavit, Computational Electromagnetism, Academic Press, San Diego, 1998. [4] C. Carstensen, S. Funken, W. Hackbusch, Ronald H.W. Hoppe, P. Monk (Eds.), Computational Electromagnetics, Lecture Notes in Computational Science and Engineering 28, Springer-Verlag, Berlin, 2003. [5] Q. Chen, M. Katsurai, P.H. Aoyagi, An FDTD formulation for dispersive media using a current density, IEEE Trans. Antennas Propag. 46 (1998) 1739–1746. [6] P. Ciarlet Jr., J. Zou, Fully discrete finite element approaches for time-dependent Maxwells equations, Numer. Math. 82 (1999) 193–219. [7] G.C. Cohen, E. Heikkola, P. Joly, P. Neittaanma¨ki (Eds.), Mathematical and Numerical Aspects of Wave Propagation—WAVES 2003, SpringerVerlag, Berlin, 2003. [8] M. Costabel, M. Dauge, Singularities of electromagnetic fields in polyhedral domains, Arch. Ration. Mech. Anal. 151 (2000) 221–276. [9] S.A. Cummer, An analysis of new and existing FDTD methods for isotropic cold plasma and a method for improving their accuracy, IEEE Trans. Antennas Propag. 45 (1997) 392–400. [10] L. Demkowicz, Fully automatic hp-adaptivity for Maxwells equations, Comput. Methods Appl. Mech. Engrg. 194 (2005) 605–624. [11] L. Demkowicz, P. Monk, L. Vardapetyan, W. Rachowicz, de Rham diagram for hp finite element spaces, Comput. Math. Appl. 39 (2000) 29–38. [12] J. Gopalakrishnan, J.E. Pasciak, L. Demkowicz, Analysis of a multigrid algorithm for time harmonic Maxwell equations, SIAM J. Numer. Anal. 42 (2004) 90–108. [13] R. Hiptmair, Finite elements in computational electromagnetism, Acta Numer. 11 (2002) 237–339. [14] D. Jiao, J.-M. Jin, Time-domain finite-element modeling of dispersive media, IEEE Microw. Wirel. Compon. Lett. 11 (2001) 220–222. [15] J. Jin, The Finite Element Method in Electromagnetics, 2nd ed., John Wiley & Sons, Inc., New York, 2002. [16] J.-F. Lee, R. Lee, A. Cangellaris, Time-domain finite-element methods, IEEE Trans. Antennas Propag. 45 (1997) 430–442. [17] Q. Lin, N. Yan, Global superconvergence for Maxwells equations, Math. Comput. 69 (1999) 159–176. [18] P. Monk, Analysis of a finite element method for Maxwells equations, SIAM J. Numer. Anal. 29 (1992) 714–729. [19] P. Monk, Finite Element Methods for Maxwells Equations, Oxford University Press, 2003. [20] J.-C. Ne´de´lec, A new family of mixed finite elements in R3 , Numer. Math. 50 (1986) 57–81. [21] J.T. Oden, G.F. Carey, Finite Elements: Mathematical Aspects, Vol. IV, Prentice-Hall, 1983. [22] P.P. Silvester, R.L. Ferrari, Finite Elements for Electrical Engineers, 3rd ed., Cambridge University Press, 1996. [23] A. Taflove, C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed., Artech House, Norwood MA, 2000. [24] J.L. Young, A higher order FDTD method for EM propagation in a collisionless cold plasma, IEEE Trans. Antennas Propag. 44 (1996) 1283–1289.