SUPERCONVERGENCE ANALYSIS FOR ... - Semantic Scholar
MATHEMATICS OF COMPUTATION Volume 77, Number 262, April 2008, Pages 757–771 S 0025-5718(07)02039-X Article electronically published on November 13, 2007
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS IN DISPERSIVE MEDIA QUN LIN AND JICHUN LI
Abstract. In this paper, we consider the time dependent Maxwell’s equations in dispersive media on a bounded three-dimensional domain. Global superconvergence is obtained for semi-discrete mixed finite element methods for three most popular dispersive media models: the isotropic cold plasma, the one-pole Debye medium, and the two-pole Lorentz medium. Global superconvergence for a standard finite element method is also presented. To our best knowledge, this is the first superconvergence analysis obtained for Maxwell’s equations when dispersive media are involved.
1. Introduction Recently there is a growing interest in finite element modeling and analysis of Maxwell’s equations (e.g. [7, 5, 9, 13, 12, 15, 16, 33, 34, 8]). The readers can find more references in some recent books [2, 18, 35] and conference proceedings [1, 3, 10]. However, most work is restricted to simple media such as free space. Very few papers are devoted to dispersive media using the finite element method (FEM), though there are some publications in finite-difference time-domain (FDTD) modeling of dispersive media started since 1990 [38, Ch. 9]. We want to remark that dispersive media are ubiquitous, for example human tissue, soil, snow, ice, plasma, fiber optics 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. Applications of time-domain finite element method (TDFEM) for dispersive media have appeared only very recently [17, 32]. However, there exists no theoretical error analysis except our initial efforts [20, 21, 22]. Superconvergence of FEM is a phenomenon that the convergence rate exceeds what general cases can provide. Since the 1970s, many studies have been conducted for superconvergence, which can be achieved for smoother solutions with structured grids. More details can be found in [27, 4, 30, 19, 39] and the references therein. In 1994, Monk [34] initiated the investigation on superconvergence for Maxwell’s equations in simple media. Recently, Lin and his collaborators [29, 26] systematically developed global superconvergence, but still restricted it to simple media. To our best knowledge, there exists no work in the literature which studies the superconvergence error analysis of TDFEM for Maxwell’s equations in dispersive media. This paper intends to make an initial effort in this direction. Here we develop Received by the editor May 25, 2006 and, in revised form, January 26, 2007. 2000 Mathematics Subject Classification. Primary 65N30, 35L15, 78Mxx. Key words and phrases. Maxwell’s equations, dispersive media, superconvergence analysis. c 2007 American Mathematical Society
757
758
QUN LIN AND JICHUN LI
global superconvergence error analysis for semi-discrete standard and mixed finite element methods for Maxwell’s equations when dispersive media are involved. The proof uses the powerful integral identity technique developed by Lin’s group in the early 1990s [31, 24, 25] and applied later to many areas [23, 40, 14, 29, 30]. The rest of the paper is organized as follows. In Section 2, some notation and preliminary results are introduced. Then in Section 3, a semi-discrete mixed finite element scheme is developed for the isotropic cold plasma model. Superclose estimates are obtained on the cubic N´ed´elec curl conforming element. Similar results are derived for one-pole Debye medium and two-pole Lorentz medium. In Section 4, global superconvergence is proved for all three dispersive media by using the integral identity technique. In Section 5, we extend the global superconvergence analysis to a standard finite element method.
2. Notation and preliminary results In this paper, C (sometimes with a sub-index) denotes a generic constant, which is independent of the finite element mesh size h. We will introduce some notation to be used in this paper. We define H(curl; Ω) = H α (curl; Ω) = H0 (curl; Ω) =
{v ∈ (L2 (Ω))3 ; ∇ × v ∈ (L2 (Ω))3 }, {v ∈ (H α (Ω))3 ; ∇ × v ∈ (H α (Ω))3 }, {v ∈ H(curl; Ω); n × v = 0 on ∂Ω},
where α ≥ 0 is a real number, and (H α (Ω))3 is the standard Sobolev space equipped with the norm · α and semi-norm | · |α in a bounded polyhedral domain Ω of R3 . Specifically · 0 will mean the (L2 (Ω))3 -norm. Also H(curl; Ω) and H α (curl; Ω) are equipped with norms v0,curl
= (v20 + curl v20 )1/2 ,
vα,curl
= (v2α + curl v2α )1/2 .
We assume that the domain Ω is covered with a regular cubic mesh Th of maximum diameter h. Our mixed finite element spaces are [36]: (1) Uh = {φ ∈ (L2 (Ω))3 ; φ|e ∈ Qk,k−1,k−1 × Qk−1,k,k−1 × Qk−1,k−1,k , ∀ e ∈ Th }, (2) Vh = {ψ ∈ H(curl; Ω); ψ|e ∈ Qk−1,k,k × Qk,k−1,k × Qk,k,k−1 , ∀ e ∈ Th }, where Qi,j,k denotes the space of polynomials whose degrees are less than or equal to i, j, k for x, y, z, respectively. Note that we have [33, p. 114]: (3)
∇ × Vh ⊂ Uh .
Let Ph E ∈ Uh be the standard (L2 (Ω))3 projection operator defined as (4)
(Ph E − E, φ) = 0,
∀ φ ∈ Uh .
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS
759
Furthermore, we define the interpolation operator H I ∈ Vh , which satisfies [29, p. 163] (H − H I ) · τ qdl = 0, ∀ q ∈ Pk (li ), i = 1, · · · , 12, li I ((H − H ) × n) · qdσ = 0, ∀ q ∈ Qk−2,k−1 (σi ) × Qk−1,k−2 (σi ), i = 1, · · · , 6, σi (H − H I ) · qde = 0, ∀ q ∈ Qk−1,k−2,k−2 × Qk−2,k−1,k−2 × Qk−2,k−2,k−1 , e
where li , σi are edges and faces of the element e, τ is the unit tangent vector along edge li , and n is the normal vector on face σi . In this paper, we need the following proven results: Lemma 2.1 ([29, Lemma 3.1]). (∇ × (H − H I )) · φdΩ = O(hk+1 )||H||k+2 ||φ||0
∀ φ ∈ Uh .
Ω
Lemma 2.2 ([29, Lemma 3.2]). (H − H I ) · ψdΩ = O(hk+1 )||H||k+1 ||ψ||0
∀ ψ ∈ Vh .
Ω
Lemma 2.3 ([37, p. 13]). Let f ∈ L1 (0, T ) be a non-negative function, and let g and ϕ be continuous functions on [0, T ]. Moreover, g is non-decreasing. If ϕ satisfies t
ϕ(t) ≤ g(t) +
f (τ )ϕ(τ )dτ
∀ t ∈ [0, T ].
0
Then
ϕ(t) ≤ g(t) exp(
t
f (τ ) dτ )
∀ t ∈ [0, T ].
0
3. Superclose estimates In this section, we shall develop the superclose estimates for all three popular dispersive media models. 3.1. Isotropic cold plasma. The governing equations that describe electromagnetic wave propagation in isotropic non-magnetized cold electron plasma are [11, 6] ∂E 0 (5) = ∇ × H − J, ∂t ∂H (6) = −∇ × E, µ0 ∂t ∂J + νJ = 0 ωp2 E, (7) ∂t where E is the electric field, H is the magnetic field, 0 is the permittivity of free space, µ0 is the permeability of free space, J is the polarization current density, ωp is the plasma frequency, and ν is the electron-neutral collision frequency. Solving (7) with the assumption that the initial electron velocity is zero leads to [6, equation (8)] t t 2 −νt νs 2 e E(x, s)ds = 0 ωp e−ν(t−s) E(x, s)ds, (8) J (x, t; E) = 0 ωp e 0
0
760
QUN LIN AND JICHUN LI
where x ∈ Ω. In this paper we let Ω be a bounded polyhedral domain in R3 with boundary ∂Ω and unit outward normal n. Substituting (8) into (5), we obtain the following system for E and H: (9) (10)
0 E t − ∇ × H + J(E) µ0 H t + ∇ × E
=
0
in Ω × (0, T ),
= 0 in Ω × (0, T ),
where for simplicity we denote J(E) = J (x, t; E). To complete the problem, we assume that the boundary of Ω is a perfect conductor [35]: (11)
n × E = 0 on ∂Ω × (0, T ),
and the initial conditions are (12)
E(x, 0) = E 0 (x)
and H(x, 0) = H 0 (x)
for any x ∈ Ω,
where E 0 and H 0 are some given functions. Furthermore, H 0 satisfies (13)
∇ · (µ0 H 0 ) = 0 in Ω,
H 0 · n = 0 on ∂Ω.
Assuming the existence of smooth solutions to (9)-(13), we obtain the weak formulation: find the solution (E, H) ∈ [C 1 (0, T ; (L2 (Ω))3 ) ∩ C 0 (0, T ; H(curl; Ω))]2 of (9)-(13) such that (14) (15)
0 (E t , φ) − (∇ × H, φ) + (J(E), φ) µ0 (H t , ψ) + (E, ∇ × ψ)
= =
0 ∀ φ ∈ (L2 (Ω))3 , 0 ∀ ψ ∈ H(curl; Ω)
for 0 < t ≤ T with the initial conditions (12). Notice that the boundary condition (11) is used in deriving (15) since (∇ × E, ψ) = n × E, ψ ∂Ω + (E, ∇ × ψ). Now we can construct our semi-discrete mixed method for solving (14)-(15): find (E h , H h ) ∈ [C 1 (0, T ; Uh ) ∩ C 1 (0, T ; Vh )]2 such that (16)
0 (E ht , φh ) − (∇ × H h , φh ) + (J(E h ), φh ) = 0 ∀ φh ∈ Uh ,
(17)
µ0 (H ht , ψ h ) + (E h , ∇ × ψ h ) = 0 ∀ ψ h ∈ Vh
for 0 < t ≤ T, subject to the initial conditions (18)
E h (x, 0) = Ph E 0 (x) and
H h (x, 0) = H I0 (x),
where H I0 ∈ Vh is the interpolation of H 0 defined in §2. Note that (16)-(18) is a system of linear ordinary differential equations, which guarantees the existence and uniqueness of solutions. Theorem 3.1. Let (E(t), H(t)) and (E h (t), H h (t)) be the solutions of (14)-(15) and (16)-(18) at time t, respectively. Then there is a constant C = C(0 , µ0 , ωp , ν), independent of the mesh size h, such that
(19)
µ0 ||(H I − H h )(t)||20 + 0 ||(Ph E − E h )(t)||20 t ≤ Ch2(k+1) [||H(t)||2k+2 + ||H t (t)||2k+1 ]dt, 0
where k is the degree of edge elements in the space V h .
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS
761
Proof. Subtracting (16)-(17) from (14)-(15) with φ = φh and ψ = ψ h , respectively, we have the error equations: (20)
0 ((E − E h )t , φh ) − (∇ × (H − H h ), φh ) + (J(E − E h ), φh ) = 0 ∀ φh ∈ Uh ,
(21)
µ0 ((H − H h )t , ψ h ) + (E − E h , ∇ × ψ h ) = 0 ∀ ψ h ∈ Vh .
Denote ξ(t) = (Ph E − E h )(t), η(t) = (H I − H h )(t). Choosing φh = ξ, ψ h = η in (20)-(21), and rearranging terms lead to 0 (ξt , ξ) − (∇ × η, ξ) = 0 ((Ph E − E)t , ξ) − (∇ × (H I − H), ξ) + (J(E h − E), ξ), µ0 (ηt , η) + (ξ, ∇ × η) = µ0 ((H I − H)t , η) + (Ph E − E, ∇ × η). Adding the two equations above, we obtain
(22)
1 d (µ0 ||η(t)||20 + 0 ||ξ(t)||20 ) 2 dt = −(∇ × (H I − H), ξ) + (J(E h − E), ξ) + µ0 ((H I − H)t , η) =
3
(Err)i
i=1
where we used the definition of operator Ph and the fact that ∇ × Vh ⊂ Uh . Below we will constantly use the basic arithmetic-geometric mean inequality |ab| ≤ δa2 +
(23)
1 2 b , 4δ
for any constant δ > 0. Now we will estimate (Err)i one by one for i = 1, 2, 3. Using Lemma 2.1 and (23), we have (Err)1
= −(∇ × (H I − H), ξ) ≤ Chk+1 ||H(t)||k+2 ||ξ(t)||0 ≤ δ2 0 ||ξ(t)||20 +
C1 h2(k+1) ||H(t)||2k+2 . 4δ2 0
By the linearity of J and the definition of Ph , we obtain (Err)2 (24)
= (J(E h − Ph E) + J(Ph E − E), ξ) = −(J(ξ), ξ) 1 ||J(ξ)||20 . ≤ δ3 0 ||ξ(t)||20 + 4δ3 0
Using the definition of J and the Cauchy-Schwarz inequality, we have t |0 ωp2 e−ν(t−s) ξ(x, s)ds|2 dΩ ||J(ξ)||20 = Ω 0 t t 2 4 |e−ν(t−s) |2 ds)( |ξ(x, s)|2 ds)dΩ ≤ 0 ωp ( Ω 0 0 t 1 2 4 −2νt (1 − e )( |ξ(x, s)|2 ds)dΩ = 0 ωp 2ν Ω 0 20 ωp4 t ≤ (25) ||ξ(s)||20 ds, 2ν 0
762
QUN LIN AND JICHUN LI
which together with (24) gives δ3 0 ||ξ(t)||20
(Err)2 ≤
0 ωp4 + 8δ3 ν
t
||ξ(s)||20 ds. 0
Using Lemma 2.2 and (23), we obtain (Err)3
= µ0 ((H I − H)t , η) ≤ Cµ0 hk+1 ||H t (t)||k+1 ||η(t)||0 ≤ δ4 µ0 ||η(t)||20 +
C2 µ0 h2(k+1) ||H t (t)||2k+1 . 4δ4
Combining all the estimates obtained for (Err)i , i = 1, 2, 3, we shall have
(26)
d µ0 0 ( ||η(t)||20 + ||ξ(t)||20 ) ≤ (δ2 + δ3 )0 ||ξ(t)||20 + δ4 µ0 ||η(t)||20 dt 2 2 0 ωp4 t C1 h2(k+1) C2 µ0 h2(k+1) 2 + ||H(t)||k+2 + ||ξ(s)||20 ds + ||H t (t)||2k+1 . 4δ2 0 8δ3 ν 0 4δ4
Integrating both sides of (26) with respect to t and using the fact that ξ(0) = η(0) = 0, we obtain t µ0 ||η(t)||20 + 0 ||ξ(t)||20 ≤ C3 (µ0 ||η(s)||20 + 0 ||ξ(s)||20 )ds 0 t (27) [||H(t)||2k+2 + ||H t (t)||2k+1 ]dt, +C4 h2(k+1) 0
where we denote C3 = 2 · max{δ2 + δ3 +
tωp4 , δ4 }, 8δ3 ν
C4 = 2 · max{
C1 C2 µ0 , }. 4δ2 0 4δ4
By the Gronwall’s inequality (Lemma 2.3), we obtain t [||H(t)||2k+2 + ||H t (t)||2k+1 ]dt, (28) µ0 ||η(t)||20 + 0 ||ξ(t)||20 ≤ Ch2(k+1) 0
which completes the proof.
3.2. Debye medium. For the single pole model of Debye, the governing equations can be written as [22]: find E and H, which satisfy (29) (30)
0 ∞ E t − ∇ × H +
(s − ∞ )0 ˜ E − J(E) = 0 in Ω × (0, T ), t0 µ0 H t + ∇ × E = 0 in Ω × (0, T ),
with the same boundary and initial conditions as those stated previously for plasma. Here we introduce the pseudo polarization current t (t−s) ˜ ˜ (x, t; E) = (s − ∞ )0 (31) J(E) ≡J e− t0 E(x, s)ds 2 t0 0 in order to carry out our earlier analysis for plasma easily to the Debye medium. Furthermore, ∞ is the permittivity at infinite frequency, s is the permittivity at zero frequency, t0 is the relaxation time, and the rest have the same meaning as those stated previously for the plasma model.
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS
763
From (29)-(30), we can easily obtain the weak formulation: find the solution (E, H) ∈ [C 1 (0, T ; (L2 (Ω))3 ) ∩ C 0 (0, T ; H(curl; Ω))]2 such that (32)
0 ∞ (E t , φ) − (∇ × H, φ) +
(33)
(s − ∞ )0 ˜ (E, φ) − (J(E), φ) = 0 ∀ φ ∈ (L2 (Ω))3 , t0 µ0 (H t , ψ) + (E, ∇ × ψ) = 0 ∀ ψ ∈ H(curl; Ω)
for 0 < t ≤ T with the initial conditions E(x, 0) = E 0 (x)
(34)
and
H(x, 0) = H 0 (x).
The semi-discrete mixed finite element scheme for our Debye model can be formulated as follows: (E h , H h ) ∈ C 1 (0, T ; Uh ) × C 1 (0, T ; Vh ) such that (35)
0 ∞ (E ht , φh ) − (∇ × H h , φh ) +
(36)
(s − ∞ )0 h ˜ h ), φh ) = 0 ∀ φh ∈ Uh , (E , φh ) − (J(E t0 µ0 (H ht , ψ h ) + (E h , ∇ × ψ h ) = 0 ∀ ψ h ∈ Vh
for 0 < t ≤ T, subject to the initial conditions (37)
E h (x, 0) = Ph E 0 (x) and
H h (x, 0) = H I0 (x).
Theorem 3.2. Let (E(t), H(t)) and (E h (t), H h (t)) be the solutions of (32)-(33) and (35)-(36) at time t, respectively. Then there is a constant C = C(0 , µ0 , ∞ , s , t0 ), independent of the mesh size h, such that
(38)
µ0 ||(H I − H h )(t)||20 + 0 ∞ ||(Ph E − E h )(t)||20 t 2(k+1) ≤ Ch [||H(t)||2k+2 + ||H t (t)||2k+1 ]dt, 0
where k is the degree of edge elements in the space V h . Proof. Subtracting (35)-(36) from (32)-(33) gives the error equations: (39)
0 ∞ ((E − E h )t , φh ) − (∇ × (H − H h ), φh ) +
(s − ∞ )0 (E − E h , φh ) t0
˜ −(J(E − E h ), φh ) = 0 ∀ φh ∈ Uh , (40)
µ0 ((H − H h )t , ψ h ) + (E − E h , ∇ × ψ h ) = 0 ∀ ψ h ∈ Vh .
Choosing φh = ξ, ψ h = η in (39)-(40), and rearranging terms lead to (s − ∞ )0 (ξ, ξ) = 0 ∞ ((Ph E − E)t , ξ) t0 (s − ∞ )0 ˜ −(∇ × (H I − H), ξ) + (Ph E − E, ξ) + (J(E − E h ), ξ), t0 µ0 (ηt , η) + (ξ, ∇ × η) = µ0 ((H I − H)t , η) + (Ph E − E, ∇ × η).
0 ∞ (ξt , ξ) − (∇ × η, ξ) +
764
QUN LIN AND JICHUN LI
Adding the above two equations together and using the definition of Ph , we obtain d 0 ∞ µ0 (s − ∞ )0 ( ||ξ(t)||20 + ||η(t)||20 ) + ||ξ(t)||20 dt 2 2 t0 ˜ = −(∇ × (H I − H), ξ) + (J(E − E h ), ξ) + µ0 ((H I − H)t , η) (41) =
3 (Err)i . i=1
The rest of the proof is omitted because it is very similar to the plasma case.
3.3. Lorentz medium. The Lorentzian two pole model can be represented by the governing equations [22]: (42) (43)
ˆ (E) 0 ∞ E t − ∇ × H + J µ0 H t + ∇ × E
= 0 in Ω × (0, T ), = 0 in Ω × (0, T ),
with the same boundary and initial conditions as those stated previously for plasma. ˆ to represent the polarization current for the Lorentz medium: We use J t ˆ ˆ (x, t; E) = β˜ ≡ J (44) J(E) e−δ(t−s) · sin(γ − α(t − s)) · E(x, s)ds, 0 t √ jγ −(δ+jα)(t−s) ˜ (45) e E(x, s)ds) ≡ Im(J(E)), j = −1, = Im(βe 0
2 where β˜ = (s − ∞ )0 ω13 / ω12 − ν4 , and Im(A) denotes the imaginary part of a general complex number A. Furthermore, in addition to the notation defined earlier, ω1 is the resonant frequency, ν is the damping coefficient, δ = ν2 , α = ω12 − δ 2 and ejγ = ωδ1 + j ωα1 . Note that ω1 > δ in real applications [17]. From (42)-(43), we can easily obtain the weak formulation: find the solution (E, H) ∈ [C 1 (0, T ; (L2 (Ω))3 ) ∩ C 0 (0, T ; H(curl; Ω))]2 such that (46) (47)
ˆ φ) 0 ∞ (E t , φ) − (∇ × H, φ) + (J(E), µ0 (H t , ψ) + (E, ∇ × ψ)
= =
0 ∀ φ ∈ (L2 (Ω))3 , 0 ∀ ψ ∈ H(curl; Ω)
for 0 < t ≤ T with the initial conditions (48)
E(x, 0) = E 0 (x)
and
H(x, 0) = H 0 (x).
The semi-discrete mixed finite element scheme for the Lorentz medium is formulated as follows: find E h ∈ Uh , H h ∈ Vh such that (49)
ˆ h ), φh ) 0 ∞ (E ht , φh ) − (∇ × H h , φh ) + (J(E
=
0 ∀ φh ∈ Uh ,
(50)
µ0 (H ht , ψ h )
+ (E , ∇ × ψ h )
=
0 ∀ ψ h ∈ Vh
h
for 0 < t ≤ T, subject to the initial conditions (51)
E h (x, 0) = Ph E 0 (x) and
H h (x, 0) = H I0 (x).
Theorem 3.3. Let (E(t), H(t)) and (E h (t), H h (t)) be the solutions of (46)-(47) and (49)-(50) at time t, respectively. Then there is a constant C = C(0 , µ0 , s , ∞ , ω1 , ν),
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS
765
independent of the mesh size h, such that µ0 ||(H I − H h )(t)||20 + 0 ∞ ||(Ph E − E h )(t)||20 t ≤ Ch2(k+1) [||H(t)||2k+2 + ||H t (t)||2k+1 ]dt,
(52)
0
where k is the degree of edge elements in the space V h . Proof. Subtracting (49)-(50) from (46)-(47) with φ = φh = ξ(t) and ψ = ψ h = η(t), respectively, we can obtain the error equations: 0 ∞ (ξt , ξ) − (∇ × η, ξ) = 0 ∞ ((Ph E − E)t , ξ) ˆ − E h ), ξ), − (∇ × (H I − H), ξ) − (J(E µ0 (ηt , η) + (ξ, ∇ × η) = µ0 ((H I − H)t , η) + (Ph E − E, ∇ × η). Adding the above two equations and using the definition of Ph , we obtain d µ0 0 ∞ ( ||η(t)||20 + ||ξ(t)||20 ) dt 2 2 ˆ = −(∇ × (H I − H), ξ) − (J(E − E h ), ξ) + µ0 ((H I − H)t , η) (53) =
3 (Err)i . i=1
The estimates of (Err)i follow the proof for the plasma case. The only different term is (Err)2 . ˆ and the definition of Ph , we obtain By the linearity of J ˆ ˆ (Err)2 = −(J(E − Ph E) + J(Ph E − E h ), ξ) = −(J(ξ), ξ) ≤
(54)
It is easy to see that
δ3 0 ∞ ||ξ(t)||20 +
1 ˆ (ξ)||2 . ||J 0 4δ3 0 ∞
t |ejγ e−(δ+jα)(t−s) ξ(x, s)ds|2 dΩ) Ω 0 t 1 2 ||ξ(s)||20 ds, ≤ β˜ 2δ 0
ˆ (ξ)||20 ≤ ||J (ξ)||20 = β˜2 ||J
which together with (54) gives
t β˜2 + ||ξ(s)||20 ds. (Err)2 ≤ 8δ3 0 ∞ δ 0 The rest of the proof follows exactly the same as the plasma case. δ3 0 ∞ ||ξ(t)||20
4. Global superconvergence To prove global superconvergence, we need some postprocessing operators introduced by Lin and Yan [29]. e) For each component wj , j = 1, 2, 3 of w ∈ Uh , we define Π2h wj |eˆ ∈ Qk,2k−1,2k−1 (ˆ such that (Π2h wj − wj )q = 0, ∀ q ∈ Qk,k−1,k−1 (ei ), i = 1, 2, 3, 4, ei
where eˆ =
4
i=1 ei .
766
QUN LIN AND JICHUN LI
ˆ 2h is defined for the function v ∈ Vh . For the Another postprocessing operator Π ˆ first component v1 of v, we define Π2h v1 |eˇ ∈ Q2k−1,k,k (ˇ e) such that ˆ 2h v1 − v1 )qdx = 0, ∀ q ∈ Pk−1 (li ), i = 1, · · · , 8, (Π li ˆ 2h v1 − v1 )qdxdy = 0, ∀ q ∈ Qk−1,k−2 (τi ), i = 1, · · · , 4, (Π τi ˆ 2h v1 − v1 )qdxdz = 0, ∀ q ∈ Qk−1,k−2 (τj ), j = 1, · · · , 4, (Π τj
ˆ 2h v1 − v1 )qdxdydz = 0, (Π
∀ q ∈ Qk−1,k−2,k−2 (ei ), i = 1, 2,
ei
where eˇ = e1 ∪e2 , li are edges parallel to the x-axis, τi , τj are surfaces perpendicular ˆ 2h can be defined similarly for the second to the z-axis or y-axis, respectively. Π and third components of v ∈ Vh . Lin and Yan [29] proved the following properties: Lemma 4.1. ˆ 2h v − v||0 ≤ Chk+1 ||v||k+1 , (i) ||Π2h w − w||0 ≤ Chk+1 ||w||k+1 , ||Π k+1 3 (Ω)) , ∀ w, v ∈ (H ˆ 2h v||0 ≤ C||v||0 , ∀ w ∈ Uh , v ∈ Vh , (ii) ||Π2h w||0 ≤ C||w||0 , ||Π ˆ ˆ 2h vI , ∀ w ∈ Uh , v ∈ Vh , (iii) Π2h w = Π2h Ph w, Π2h v = Π where Ph w ∈ Uh and vI ∈ Vh are the interpolations of w and v defined in Section 2. Using these postprocessing operators, we can achieve the following global superconvergence for all three dispersive media: Theorem 4.1. ˆ 2h H h − H||0 ||Π2h E h − E||0 + ||Π
≤ Chk+1 [||E||k+1 + ||H||k+1 + (
t
(||H||2k+2 + ||H t ||2k+1 )ds)1/2 ], 0
where k is the degree of edge elements in the space V h . Proof. By Lemma 4.1 and Theorems 3.1–3.3, we have ||Π2h E h − E||0
(55)
= ||Π2h (E h − Ph E) + (Π2h E − E)||0 ≤ C||E h − Ph E||0 + Chk+1 ||E||k+1 t ≤ Chk+1 [( (||H||2k+2 + ||H t ||2k+1 )ds)1/2 + ||E||k+1 ]. 0
Similarly, we have ˆ 2h H h − H||0 ||Π
(56)
ˆ 2h (H h − H I ) + (Π ˆ 2h H − H)||0 = ||Π ≤ C||H h − H I ||0 + Chk+1 ||H||k+1 t ≤ Chk+1 [( (||H||2k+2 + ||H t ||2k+1 )ds)1/2 + ||H||k+1 ], 0
which, along with (55), completes the proof.
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS
767
5. Extension to a standard finite element method In this section, we will extend the global superconvergence to a standard finite element method. For simplicity, we will restrict our discussion to the cold plasma model, generalizations to other models are similar. Instead of solving the coupled system (5)-(7) with both the electric and magnetic fields as unknowns, we eliminate H by taking the time derivative of (5) and using (6)-(8), to obtain the second order electric field equation (57)
2 0 E tt + ∇ × (µ−1 0 ∇ × E) + 0 ωp E − νJ(E) = 0,
with boundary condition (11) and initial conditions (58)
E(x, 0) = E 0 (x)
and
E t (x, 0) = E 1 (x),
where E 1 (x) = −1 0 ∇ × H 0 (x), which is obtained from (5), (8), and (12). Multiplying (57) by φ ∈ H0 (curl; Ω), and using integration by parts [35, (3.27)], we can easily obtain the weak formulation: find E(t) ∈ H0 (curl; Ω) such that (59)
2 0 (E tt , φ) + µ−1 0 (∇ × E, ∇ × φ) + 0 ωp (E, φ)
− ν(J(E), φ) = 0 φ ∈ H0 (curl; Ω),
subject to the initial conditions (58). Taking the boundary condition (11) into account, we define V 0h = {v h ∈ V h | n × v h = 0 on ∂Ω}. Then we can formulate a standard finite element scheme for (57) as follows: find E h (t) ∈ V 0h such that h h h 0 2 (60) 0 (E htt , φ)+µ−1 0 (∇×E , ∇×φ)+0 ωp (E , φ)−ν(J(E ), φ) = 0 ∀φ ∈ V h ,
subject to the initial conditions (61)
E h (x, 0) = E I0 (x),
E ht (x, 0) = E I1 (x),
where E I0 and E I1 are the interpolations of E 0 and E 1 defined in §2, respectively. Theorem 5.1. Let E(t) and E h (t) be the solutions of (59) and (60) at time t, respectively. Then there is a constant C = C(0 , µ0 , ωp , ν), independent of the mesh size h, such that I h I h 2 2 2 0 ||(E I − E h )t (t)||20 + µ−1 0 ||∇ × (E − E )(t)||0 + 0 ωp ||(E − E )(t)||0 t ≤ Ch2(k+1) [||E||2k+2 + (||E tt ||2k+1 + ||E t ||2k+2 + ||E||2k+1 )ds]. 0
Proof. Denote ξ(t) = (E − E )(t). Subtracting (60) from (59) with φ = ξt (t), and rearranging terms, we obtain the error equation I
h
2 0 (ξtt , ξt ) + µ−1 0 (∇ × ξ, ∇ × ξt ) + 0 ωp (ξ, ξt )
(62)
I = 0 ((E I − E)tt , ξt ) + µ−1 0 (∇ × (E − E), ∇ × ξt )
+ 0 ωp2 (E I − E, ξt ) + ν(J(E − E I ), ξt ) + ν(J(ξ), ξt ).
768
QUN LIN AND JICHUN LI
Integrating (62) with respect to t and using integration by parts and the fact that ξ(0) = ξt (0) = 0, we obtain
(63)
1 2 2 2 (0 ||ξt ||20 + µ−1 0 ||∇ × ξ||0 + 0 ωp ||ξ||0 ) 2 t I ≤ 0 ((E I − E)tt , ξt )ds + µ−1 0 (∇ × (E − E), ∇ × ξ) 0 t (∇ × (E I − E)t , ∇ × ξ)ds − µ−1 0 0 t 2 ((E I − E)t , ξ)ds + 0 ωp (E I − E, ξ) − 0 ωp2 0
t
t
(J(E − E I ), ξt )ds + ν
+ν 0
(J(ξ), ξt )ds =
7
0
(Err)i .
i=1
Now we have to estimate all (Err)i , i = 1, · · · , 7. Using Lemma 2.2 and the inequality (23), we have
t 0
δ1 0 ||ξt ||20 ds 0
Chk+1 ||E tt ||k+1 ||ξt ||0 ds 0
t
≤
t
((E I − E)tt , ξt )ds ≤ 0
= 0
(Err)1
t
+ 0
C0 h2(k+1) ||E tt ||2k+1 ds. 4δ1
Using Lemma 2.1 and the inequality (23) leads to I −1 k+1 (Err)2 = µ−1 ||E||k+2 ||∇ × ξ||0 0 (∇ × (E − E), ∇ × ξ) ≤ µ0 Ch 2 ≤ δ2 µ−1 0 ||∇ × ξ||0 +
Ch2(k+1) ||E||2k+2 , 4δ2 µ0
and (Err)3 = −µ−1 0
t
(∇×(E I −E)t , ∇×ξ)ds ≤ µ−1 0
0
t
≤
δ3 µ−1 0 ||∇
×
ξ||20 ds
0
Ch2(k+1) + 4δ3 µ0
t
Chk+1 ||E t ||k+2 ||∇×ξ||0 ds 0
t
||E t ||2k+2 ds. 0
Similarly, by Lemma 2.2 and the inequality (23), we obtain = 0 ωp2 (E I − E, ξ) ≤ 0 ωp2 Chk+1 ||E||k+1 ||ξ||0
(Err)4
≤ δ4 0 ωp2 ||ξ||20 +
0 ωp2 Ch2(k+1) ||E||2k+1 , 4δ4
and (Err)5
= −0 ωp2
t
0
t
≤ 0
t
((E I − E)t , ξ)ds ≤ 0 ωp2
0 ωp2 Ch2(k+1) δ5 0 ωp2 ||ξ||20 ds + 4δ5
Chk+1 ||E t ||k+1 ||ξ||0 ds
0 t
||E t ||2k+1 ds. 0
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS
769
Note that
t e−ν(t−s) (E − E I )(x, s)ds, ξt (x, t)) (J(E − E I ), ξt ) = 0 ωp2 ( 0 t Chk+1 ||E||k+1 ||ξt ||0 ds, ≤ 0 ωp2 0
from which we have (Err)6
t
(J(E − E ), ξt )ds ≤ I
= ν ≤
0 t
t
Chk+1 ||E||k+1 ||ξt ||0 ds
ν0 ωp2 t 0
δ6 0 ||ξt ||20 ds + 0
0 ωp4 t2 Ch2(k+1) 4δ6
t
||E||2k+1 ds. 0
Finally, using the inequality (23), we have t t t ν2 2 (Err)7 = ν (J(ξ), ξt )ds ≤ δ7 0 ||ξt ||0 ds + ||J (ξ)||20 ds 4δ7 0 0 0 0 t ν0 ωp4 t t ≤ δ7 0 ||ξt ||20 ds + ||ξ||20 ds, 8δ7 0 0 where in the last step we used the estimate (25). Substituting the above estimates for (Err)i into (63), choosing constants δ2 , δ4 < 1 and absorbing the first terms in (Err)2 and (Err)4 by the corresponding terms 2 on the left hand side of (63), we obtain
(64)
1 2 2 2 (0 ||ξt ||20 + µ−1 0 ||∇ × ξ||0 + 0 ωp ||ξ||0 ) 2 t 2 2 2 ≤ C1 (0 ||ξt ||20 + µ−1 0 ||∇ × ξ||0 + 0 ωp ||ξ||0 )ds 0 t (||E tt ||2k+1 + ||E t ||2k+2 + ||E||2k+1 )ds + C2 h2(k+1) 0
+ C3 h
2(k+1)
||E||2k+2 ,
where we absorbed the explicit dependence of physical parameters into the generic constants C1 , C2 and C3 . Our proof concludes by using the Gronwall inequality (Lemma 2.3) to (64). ˆ 2h defined in the last section, we can easUsing the postprocessing operator Π ily achieve the following global superconvergence for the standard finite element method: Theorem 5.2. ˆ 2h E h − E||0 + ||(Π ˆ 2h E h − E)t ||0 ||Π
≤ Chk+1 [||E||k+2 + ||E t ||k+1 + (
t
(||E tt ||2k+1 + ||E t ||2k+2 + ||E||2k+1 )ds)1/2 ], 0
where k is the degree of edge elements in the space V h . Remark 5.1. When Ω is not a cubic domain, global superconvergence of order 1 O(hk+ 2 ) can be achieved on the almost cubic meshes by easily extending the results of Lin and Yan [29, p. 175].
770
QUN LIN AND JICHUN LI
Acknowledgements Jichun Li thanks Dr. Miguel Visbal at the U.S. Air Force Research Laboratory (AFRL) for inspiring him to study the interesting dispersive media problems during his summer research at AFRL. The authors wish to thank the referees for many constructive comments that improved the paper. References [1] M. Ainsworth, P. Davies, D. Duncan, P. Martin and B. Rynne, Topics in computational wave propagation, Lecture Notes in Computational Science and Engineering 31, Springer-Verlag, Berlin, 2003. MR2032865 (2004h:65008) [2] A. Bossavit, Computational Electromagnetism, Academic Press, San Diego, 1998. MR1488417 (99m:78001) [3] C. Carstensen, S. Funken, W. Hackbusch, Ronald H.W. Hoppe and P. Monk (eds.), Computational Electromagnetics, Lecture Notes in Computational Science and Engineering 28, Springer-Verlag, Berlin, 2003. MR1995092 (2004b:78001) [4] C.M. Chen and Y.Q. Huang, High Accuracy Theory of Finite Element Methods (in Chinese), Hunan Science Press, China, 1995. [5] M.-H. Chen, B. Cockburn and F. Reitich, High-order RKDG methods for computational electromagnetics, J. Sci. Comp. 22/23 (2005) 205-226. MR2142195 (2005m:65208) [6] Q. Chen, M. Katsurai and P.H. Aoyagi, An FDTD formulation for dispersive media using a current density, IEEE Trans. Antennas Propagat. 46 (1998) 1739-1746. [7] Z. Chen, Q. Du and J. Zou, Finite element methods with matching and nonmatching meshes for Maxwell equations with discontinuous coefficients, SIAM J. Numer. Anal. 37 (2000) 15421570. MR1759906 (2001h:78044) [8] P. Ciarlet, Jr. and J. Zou, Fully discrete finite element approaches for time-dependent Maxwell’s equations, Numer. Math. 82 (1999) 193-219. MR1685459 (2000c:65083) [9] B. Cockburn, F. Li and C.-W. Shu, Locally divergence-free discontinuous Galerkin methods for the Maxwell equations, J. Comp. Phys. 194 (2004) 588-610. MR2034859 (2004j:78024) [10] G.C. Cohen, E. Heikkola, P. Joly and P. Neittaanm¨ aki (eds.), Mathematical and numerical aspects of wave propagation — WAVES 2003, Springer-Verlag, Berlin, 2003. MR2077970 (2005b:00020) [11] 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 Propagat. 45 (1997) 392-400. [12] L. Demkowicz, Fully automatic hp-adaptivity for Maxwell’s equations, Comput. Methods Appl. Mech. Engrg. 194 (2005) 605-624. MR2105184 (2005m:65296) [13] L. Demkowicz, P. Monk, Ch. Schwab and L. Vardapetyan, Maxwell eigenvalues and discrete compactness in two dimensions, Comput. Math. Appl. 40 (2000) 589-605. MR1772657 (2002e:65166) [14] R.E. Ewing, Y. Lin, T. Sun, J. Wang, and S. Zhang, Sharp L2 -error estimates and superconvergence of mixed finite element methods for non-Fickian flows in porous media, SIAM J. Numer. Anal. 40 (2002) 1538-1560. MR1951906 (2003m:65177) [15] R. Hiptmair, Finite elements in computational electromagnetism, Acta Numerica 11 (2002) 237-339. MR2009375 (2004k:78028) [16] P. Houston, I. Perugia and D. Sch¨ otzau, Energy norm a posteriori error estimation for mixed discontinuous Galerkin approximations of the Maxwell operator, Comput. Methods Appl. Mech. Engrg. 194 (2005) 499-510. MR2105178 (2005h:78027) [17] D. Jiao and J.-M. Jin, Time-domain finite-element modeling of dispersive media, IEEE Microwave and Wireless Components Letters 11 (2001) 220-222. [18] J. Jin, The Finite Element Method in Electromagnetics, 2nd ed., John Wiley & Sons, Inc., New York, 2002. MR1903357 (2004b:78019) [19] M. Krizek and P. Neittaanmaki, Bibliography on Superconvergence, in “Finite Element Methods”, Lecture Notes in Pure and Applied Math., Vol. 196, 1998, pp. 315-348. MR1602730 [20] J. Li, Error analysis of finite element methods for 3-D Maxwell’s equations in dispersive media, Journal of Computational and Applied Mathematics, 188 (2006) 107-120. MR2192630 (2006m:65279)
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS
771
[21] J. Li and Y. Chen, 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. MR2229840 [22] J. Li, Error analysis of fully discrete mixed finite element schemes for 3-D Maxwell’s equations in dispersive media, preprint, Sept. 2005. [23] J. Li and M.F. Wheeler, Uniform convergence and superconvergence of mixed finite element methods on anisotropically refined grids, SIAM J. Numer. Anal. 38 (2000) 770-798. MR1781203 (2001f:65137) [24] Q. Lin, J. Li and A. Zhou, rectangle test for Ciarlet-Raviart scheme, Prof. of Sys. Sci. & Sys. Engrg., Great Wall (H.K.) Culture Publish Co., 1991, pp. 230-231. [25] Q. Lin, J. Li and A. Zhou, A rectangle test for the Stokes equations, Prof. of Sys. Sci. & Sys. Engrg., Great Wall (H.K.) Culture Publish Co., 1991, pp. 240-241. [26] Q. Lin and J.F. Lin, High accuracy approximation of mixed finite element for 2-D Maxwell equations. (Chinese) Acta Math. Sci. Ser. A Chin. Ed. 23 (2003), no. 4, 499-503. MR2090553 (2005i:65145) [27] Q. Lin and J. Xu, Linear finite elements with high accuracy, J. Comp. Math. 3 (1985) 115-133. MR854355 (87k:65141) [28] Q. Lin and N. Yan, Superconvergence of mixed element methods for Maxwell’s equations (in Chinese), Gongcheng Shuxue Xuebao 13 (1996), suppl., 1-10. MR1437477 (98c:65186) [29] Q. Lin and N. Yan, Global superconvergence for Maxwell’s equations, Math. Comp. 69 (1999) 159-176. MR1654029 (2000i:65131) [30] Q. Lin and N. Yan, The Construction and Analysis of High Accurate Finite Element Methods (in Chinese), Hebei University Press, Hebei, China, 1996. [31] Q. Lin, N. Yan and A. Zhou, A rectangle test for interpolated finite elements, Prof. of Sys. Sci. & Sys. Engrg., Great Wall (H.K.) Culture Publish Co., 1991, pp. 217-229. [32] T. Lu, P. Zhang and W. Cai, Discontinuous Galerkin methods for dispersive and lossy Maxwell’s equations and PML boundary conditions, J. Comp. Phys. 200 (2004) 549-580. MR2095277 (2005e:78025) [33] P. Monk, An analysis of N´ed´ elec’s method for the special discretization of Maxwell’s equations, J. Comp. Appl. Math. 47 (1993) 101-121. MR1226366 (94g:65105) [34] P. Monk, Superconvergence of finite element approximations to Maxwell’s equations, Numer. Methods Partial Differential Equations 10 (1994) 793-812. MR1298123 (95h:65090) [35] P. Monk, Finite element methods for Maxwell’s equations, Oxford University Press, 2003. MR2059447 (2005d:65003) [36] J.-C. N´ ed´ elec, Mixed finite elements in R3 , Numer. Math. 35 (1980) 315-341. MR592160 (81k:65125) [37] A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equations, Springer-Verlag, Berlin, 1994. MR1299729 (95i:65005) [38] A. Taflove and C. Hagness, Computational Electrodynamics: The Finite-Difference TimeDomain Method, 2nd ed., Artech House, Norwood, MA, 2000. MR1338377 (96f:78001) [39] L.B. Wahlbin, Superconvergence in Galerkin Finite Element Methods, Springer, Berlin, 1995. MR1439050 (98j:65083) [40] Z. Zhang, N. Yan and T. Sun, Superconvergent derivative recovery for the intermediate finite element family of the second order, IMA Journal of Numerical Analysis 21 (2001) 643-665. MR1844362 (2002f:65171) LSEC, ICMSEC, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100080, China E-mail address: [email protected] Department of Mathematical Sciences, University of Nevada, Las Vegas, 4505 Maryland Parkway, Box 454020, Las Vegas, Nevada 89154-4020 E-mail address: [email protected]
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS IN DISPERSIVE MEDIA QUN LIN AND JICHUN LI
Abstract. In this paper, we consider the time dependent Maxwell’s equations in dispersive media on a bounded three-dimensional domain. Global superconvergence is obtained for semi-discrete mixed finite element methods for three most popular dispersive media models: the isotropic cold plasma, the one-pole Debye medium, and the two-pole Lorentz medium. Global superconvergence for a standard finite element method is also presented. To our best knowledge, this is the first superconvergence analysis obtained for Maxwell’s equations when dispersive media are involved.
1. Introduction Recently there is a growing interest in finite element modeling and analysis of Maxwell’s equations (e.g. [7, 5, 9, 13, 12, 15, 16, 33, 34, 8]). The readers can find more references in some recent books [2, 18, 35] and conference proceedings [1, 3, 10]. However, most work is restricted to simple media such as free space. Very few papers are devoted to dispersive media using the finite element method (FEM), though there are some publications in finite-difference time-domain (FDTD) modeling of dispersive media started since 1990 [38, Ch. 9]. We want to remark that dispersive media are ubiquitous, for example human tissue, soil, snow, ice, plasma, fiber optics 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. Applications of time-domain finite element method (TDFEM) for dispersive media have appeared only very recently [17, 32]. However, there exists no theoretical error analysis except our initial efforts [20, 21, 22]. Superconvergence of FEM is a phenomenon that the convergence rate exceeds what general cases can provide. Since the 1970s, many studies have been conducted for superconvergence, which can be achieved for smoother solutions with structured grids. More details can be found in [27, 4, 30, 19, 39] and the references therein. In 1994, Monk [34] initiated the investigation on superconvergence for Maxwell’s equations in simple media. Recently, Lin and his collaborators [29, 26] systematically developed global superconvergence, but still restricted it to simple media. To our best knowledge, there exists no work in the literature which studies the superconvergence error analysis of TDFEM for Maxwell’s equations in dispersive media. This paper intends to make an initial effort in this direction. Here we develop Received by the editor May 25, 2006 and, in revised form, January 26, 2007. 2000 Mathematics Subject Classification. Primary 65N30, 35L15, 78Mxx. Key words and phrases. Maxwell’s equations, dispersive media, superconvergence analysis. c 2007 American Mathematical Society
757
758
QUN LIN AND JICHUN LI
global superconvergence error analysis for semi-discrete standard and mixed finite element methods for Maxwell’s equations when dispersive media are involved. The proof uses the powerful integral identity technique developed by Lin’s group in the early 1990s [31, 24, 25] and applied later to many areas [23, 40, 14, 29, 30]. The rest of the paper is organized as follows. In Section 2, some notation and preliminary results are introduced. Then in Section 3, a semi-discrete mixed finite element scheme is developed for the isotropic cold plasma model. Superclose estimates are obtained on the cubic N´ed´elec curl conforming element. Similar results are derived for one-pole Debye medium and two-pole Lorentz medium. In Section 4, global superconvergence is proved for all three dispersive media by using the integral identity technique. In Section 5, we extend the global superconvergence analysis to a standard finite element method.
2. Notation and preliminary results In this paper, C (sometimes with a sub-index) denotes a generic constant, which is independent of the finite element mesh size h. We will introduce some notation to be used in this paper. We define H(curl; Ω) = H α (curl; Ω) = H0 (curl; Ω) =
{v ∈ (L2 (Ω))3 ; ∇ × v ∈ (L2 (Ω))3 }, {v ∈ (H α (Ω))3 ; ∇ × v ∈ (H α (Ω))3 }, {v ∈ H(curl; Ω); n × v = 0 on ∂Ω},
where α ≥ 0 is a real number, and (H α (Ω))3 is the standard Sobolev space equipped with the norm · α and semi-norm | · |α in a bounded polyhedral domain Ω of R3 . Specifically · 0 will mean the (L2 (Ω))3 -norm. Also H(curl; Ω) and H α (curl; Ω) are equipped with norms v0,curl
= (v20 + curl v20 )1/2 ,
vα,curl
= (v2α + curl v2α )1/2 .
We assume that the domain Ω is covered with a regular cubic mesh Th of maximum diameter h. Our mixed finite element spaces are [36]: (1) Uh = {φ ∈ (L2 (Ω))3 ; φ|e ∈ Qk,k−1,k−1 × Qk−1,k,k−1 × Qk−1,k−1,k , ∀ e ∈ Th }, (2) Vh = {ψ ∈ H(curl; Ω); ψ|e ∈ Qk−1,k,k × Qk,k−1,k × Qk,k,k−1 , ∀ e ∈ Th }, where Qi,j,k denotes the space of polynomials whose degrees are less than or equal to i, j, k for x, y, z, respectively. Note that we have [33, p. 114]: (3)
∇ × Vh ⊂ Uh .
Let Ph E ∈ Uh be the standard (L2 (Ω))3 projection operator defined as (4)
(Ph E − E, φ) = 0,
∀ φ ∈ Uh .
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS
759
Furthermore, we define the interpolation operator H I ∈ Vh , which satisfies [29, p. 163] (H − H I ) · τ qdl = 0, ∀ q ∈ Pk (li ), i = 1, · · · , 12, li I ((H − H ) × n) · qdσ = 0, ∀ q ∈ Qk−2,k−1 (σi ) × Qk−1,k−2 (σi ), i = 1, · · · , 6, σi (H − H I ) · qde = 0, ∀ q ∈ Qk−1,k−2,k−2 × Qk−2,k−1,k−2 × Qk−2,k−2,k−1 , e
where li , σi are edges and faces of the element e, τ is the unit tangent vector along edge li , and n is the normal vector on face σi . In this paper, we need the following proven results: Lemma 2.1 ([29, Lemma 3.1]). (∇ × (H − H I )) · φdΩ = O(hk+1 )||H||k+2 ||φ||0
∀ φ ∈ Uh .
Ω
Lemma 2.2 ([29, Lemma 3.2]). (H − H I ) · ψdΩ = O(hk+1 )||H||k+1 ||ψ||0
∀ ψ ∈ Vh .
Ω
Lemma 2.3 ([37, p. 13]). Let f ∈ L1 (0, T ) be a non-negative function, and let g and ϕ be continuous functions on [0, T ]. Moreover, g is non-decreasing. If ϕ satisfies t
ϕ(t) ≤ g(t) +
f (τ )ϕ(τ )dτ
∀ t ∈ [0, T ].
0
Then
ϕ(t) ≤ g(t) exp(
t
f (τ ) dτ )
∀ t ∈ [0, T ].
0
3. Superclose estimates In this section, we shall develop the superclose estimates for all three popular dispersive media models. 3.1. Isotropic cold plasma. The governing equations that describe electromagnetic wave propagation in isotropic non-magnetized cold electron plasma are [11, 6] ∂E 0 (5) = ∇ × H − J, ∂t ∂H (6) = −∇ × E, µ0 ∂t ∂J + νJ = 0 ωp2 E, (7) ∂t where E is the electric field, H is the magnetic field, 0 is the permittivity of free space, µ0 is the permeability of free space, J is the polarization current density, ωp is the plasma frequency, and ν is the electron-neutral collision frequency. Solving (7) with the assumption that the initial electron velocity is zero leads to [6, equation (8)] t t 2 −νt νs 2 e E(x, s)ds = 0 ωp e−ν(t−s) E(x, s)ds, (8) J (x, t; E) = 0 ωp e 0
0
760
QUN LIN AND JICHUN LI
where x ∈ Ω. In this paper we let Ω be a bounded polyhedral domain in R3 with boundary ∂Ω and unit outward normal n. Substituting (8) into (5), we obtain the following system for E and H: (9) (10)
0 E t − ∇ × H + J(E) µ0 H t + ∇ × E
=
0
in Ω × (0, T ),
= 0 in Ω × (0, T ),
where for simplicity we denote J(E) = J (x, t; E). To complete the problem, we assume that the boundary of Ω is a perfect conductor [35]: (11)
n × E = 0 on ∂Ω × (0, T ),
and the initial conditions are (12)
E(x, 0) = E 0 (x)
and H(x, 0) = H 0 (x)
for any x ∈ Ω,
where E 0 and H 0 are some given functions. Furthermore, H 0 satisfies (13)
∇ · (µ0 H 0 ) = 0 in Ω,
H 0 · n = 0 on ∂Ω.
Assuming the existence of smooth solutions to (9)-(13), we obtain the weak formulation: find the solution (E, H) ∈ [C 1 (0, T ; (L2 (Ω))3 ) ∩ C 0 (0, T ; H(curl; Ω))]2 of (9)-(13) such that (14) (15)
0 (E t , φ) − (∇ × H, φ) + (J(E), φ) µ0 (H t , ψ) + (E, ∇ × ψ)
= =
0 ∀ φ ∈ (L2 (Ω))3 , 0 ∀ ψ ∈ H(curl; Ω)
for 0 < t ≤ T with the initial conditions (12). Notice that the boundary condition (11) is used in deriving (15) since (∇ × E, ψ) = n × E, ψ ∂Ω + (E, ∇ × ψ). Now we can construct our semi-discrete mixed method for solving (14)-(15): find (E h , H h ) ∈ [C 1 (0, T ; Uh ) ∩ C 1 (0, T ; Vh )]2 such that (16)
0 (E ht , φh ) − (∇ × H h , φh ) + (J(E h ), φh ) = 0 ∀ φh ∈ Uh ,
(17)
µ0 (H ht , ψ h ) + (E h , ∇ × ψ h ) = 0 ∀ ψ h ∈ Vh
for 0 < t ≤ T, subject to the initial conditions (18)
E h (x, 0) = Ph E 0 (x) and
H h (x, 0) = H I0 (x),
where H I0 ∈ Vh is the interpolation of H 0 defined in §2. Note that (16)-(18) is a system of linear ordinary differential equations, which guarantees the existence and uniqueness of solutions. Theorem 3.1. Let (E(t), H(t)) and (E h (t), H h (t)) be the solutions of (14)-(15) and (16)-(18) at time t, respectively. Then there is a constant C = C(0 , µ0 , ωp , ν), independent of the mesh size h, such that
(19)
µ0 ||(H I − H h )(t)||20 + 0 ||(Ph E − E h )(t)||20 t ≤ Ch2(k+1) [||H(t)||2k+2 + ||H t (t)||2k+1 ]dt, 0
where k is the degree of edge elements in the space V h .
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS
761
Proof. Subtracting (16)-(17) from (14)-(15) with φ = φh and ψ = ψ h , respectively, we have the error equations: (20)
0 ((E − E h )t , φh ) − (∇ × (H − H h ), φh ) + (J(E − E h ), φh ) = 0 ∀ φh ∈ Uh ,
(21)
µ0 ((H − H h )t , ψ h ) + (E − E h , ∇ × ψ h ) = 0 ∀ ψ h ∈ Vh .
Denote ξ(t) = (Ph E − E h )(t), η(t) = (H I − H h )(t). Choosing φh = ξ, ψ h = η in (20)-(21), and rearranging terms lead to 0 (ξt , ξ) − (∇ × η, ξ) = 0 ((Ph E − E)t , ξ) − (∇ × (H I − H), ξ) + (J(E h − E), ξ), µ0 (ηt , η) + (ξ, ∇ × η) = µ0 ((H I − H)t , η) + (Ph E − E, ∇ × η). Adding the two equations above, we obtain
(22)
1 d (µ0 ||η(t)||20 + 0 ||ξ(t)||20 ) 2 dt = −(∇ × (H I − H), ξ) + (J(E h − E), ξ) + µ0 ((H I − H)t , η) =
3
(Err)i
i=1
where we used the definition of operator Ph and the fact that ∇ × Vh ⊂ Uh . Below we will constantly use the basic arithmetic-geometric mean inequality |ab| ≤ δa2 +
(23)
1 2 b , 4δ
for any constant δ > 0. Now we will estimate (Err)i one by one for i = 1, 2, 3. Using Lemma 2.1 and (23), we have (Err)1
= −(∇ × (H I − H), ξ) ≤ Chk+1 ||H(t)||k+2 ||ξ(t)||0 ≤ δ2 0 ||ξ(t)||20 +
C1 h2(k+1) ||H(t)||2k+2 . 4δ2 0
By the linearity of J and the definition of Ph , we obtain (Err)2 (24)
= (J(E h − Ph E) + J(Ph E − E), ξ) = −(J(ξ), ξ) 1 ||J(ξ)||20 . ≤ δ3 0 ||ξ(t)||20 + 4δ3 0
Using the definition of J and the Cauchy-Schwarz inequality, we have t |0 ωp2 e−ν(t−s) ξ(x, s)ds|2 dΩ ||J(ξ)||20 = Ω 0 t t 2 4 |e−ν(t−s) |2 ds)( |ξ(x, s)|2 ds)dΩ ≤ 0 ωp ( Ω 0 0 t 1 2 4 −2νt (1 − e )( |ξ(x, s)|2 ds)dΩ = 0 ωp 2ν Ω 0 20 ωp4 t ≤ (25) ||ξ(s)||20 ds, 2ν 0
762
QUN LIN AND JICHUN LI
which together with (24) gives δ3 0 ||ξ(t)||20
(Err)2 ≤
0 ωp4 + 8δ3 ν
t
||ξ(s)||20 ds. 0
Using Lemma 2.2 and (23), we obtain (Err)3
= µ0 ((H I − H)t , η) ≤ Cµ0 hk+1 ||H t (t)||k+1 ||η(t)||0 ≤ δ4 µ0 ||η(t)||20 +
C2 µ0 h2(k+1) ||H t (t)||2k+1 . 4δ4
Combining all the estimates obtained for (Err)i , i = 1, 2, 3, we shall have
(26)
d µ0 0 ( ||η(t)||20 + ||ξ(t)||20 ) ≤ (δ2 + δ3 )0 ||ξ(t)||20 + δ4 µ0 ||η(t)||20 dt 2 2 0 ωp4 t C1 h2(k+1) C2 µ0 h2(k+1) 2 + ||H(t)||k+2 + ||ξ(s)||20 ds + ||H t (t)||2k+1 . 4δ2 0 8δ3 ν 0 4δ4
Integrating both sides of (26) with respect to t and using the fact that ξ(0) = η(0) = 0, we obtain t µ0 ||η(t)||20 + 0 ||ξ(t)||20 ≤ C3 (µ0 ||η(s)||20 + 0 ||ξ(s)||20 )ds 0 t (27) [||H(t)||2k+2 + ||H t (t)||2k+1 ]dt, +C4 h2(k+1) 0
where we denote C3 = 2 · max{δ2 + δ3 +
tωp4 , δ4 }, 8δ3 ν
C4 = 2 · max{
C1 C2 µ0 , }. 4δ2 0 4δ4
By the Gronwall’s inequality (Lemma 2.3), we obtain t [||H(t)||2k+2 + ||H t (t)||2k+1 ]dt, (28) µ0 ||η(t)||20 + 0 ||ξ(t)||20 ≤ Ch2(k+1) 0
which completes the proof.
3.2. Debye medium. For the single pole model of Debye, the governing equations can be written as [22]: find E and H, which satisfy (29) (30)
0 ∞ E t − ∇ × H +
(s − ∞ )0 ˜ E − J(E) = 0 in Ω × (0, T ), t0 µ0 H t + ∇ × E = 0 in Ω × (0, T ),
with the same boundary and initial conditions as those stated previously for plasma. Here we introduce the pseudo polarization current t (t−s) ˜ ˜ (x, t; E) = (s − ∞ )0 (31) J(E) ≡J e− t0 E(x, s)ds 2 t0 0 in order to carry out our earlier analysis for plasma easily to the Debye medium. Furthermore, ∞ is the permittivity at infinite frequency, s is the permittivity at zero frequency, t0 is the relaxation time, and the rest have the same meaning as those stated previously for the plasma model.
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS
763
From (29)-(30), we can easily obtain the weak formulation: find the solution (E, H) ∈ [C 1 (0, T ; (L2 (Ω))3 ) ∩ C 0 (0, T ; H(curl; Ω))]2 such that (32)
0 ∞ (E t , φ) − (∇ × H, φ) +
(33)
(s − ∞ )0 ˜ (E, φ) − (J(E), φ) = 0 ∀ φ ∈ (L2 (Ω))3 , t0 µ0 (H t , ψ) + (E, ∇ × ψ) = 0 ∀ ψ ∈ H(curl; Ω)
for 0 < t ≤ T with the initial conditions E(x, 0) = E 0 (x)
(34)
and
H(x, 0) = H 0 (x).
The semi-discrete mixed finite element scheme for our Debye model can be formulated as follows: (E h , H h ) ∈ C 1 (0, T ; Uh ) × C 1 (0, T ; Vh ) such that (35)
0 ∞ (E ht , φh ) − (∇ × H h , φh ) +
(36)
(s − ∞ )0 h ˜ h ), φh ) = 0 ∀ φh ∈ Uh , (E , φh ) − (J(E t0 µ0 (H ht , ψ h ) + (E h , ∇ × ψ h ) = 0 ∀ ψ h ∈ Vh
for 0 < t ≤ T, subject to the initial conditions (37)
E h (x, 0) = Ph E 0 (x) and
H h (x, 0) = H I0 (x).
Theorem 3.2. Let (E(t), H(t)) and (E h (t), H h (t)) be the solutions of (32)-(33) and (35)-(36) at time t, respectively. Then there is a constant C = C(0 , µ0 , ∞ , s , t0 ), independent of the mesh size h, such that
(38)
µ0 ||(H I − H h )(t)||20 + 0 ∞ ||(Ph E − E h )(t)||20 t 2(k+1) ≤ Ch [||H(t)||2k+2 + ||H t (t)||2k+1 ]dt, 0
where k is the degree of edge elements in the space V h . Proof. Subtracting (35)-(36) from (32)-(33) gives the error equations: (39)
0 ∞ ((E − E h )t , φh ) − (∇ × (H − H h ), φh ) +
(s − ∞ )0 (E − E h , φh ) t0
˜ −(J(E − E h ), φh ) = 0 ∀ φh ∈ Uh , (40)
µ0 ((H − H h )t , ψ h ) + (E − E h , ∇ × ψ h ) = 0 ∀ ψ h ∈ Vh .
Choosing φh = ξ, ψ h = η in (39)-(40), and rearranging terms lead to (s − ∞ )0 (ξ, ξ) = 0 ∞ ((Ph E − E)t , ξ) t0 (s − ∞ )0 ˜ −(∇ × (H I − H), ξ) + (Ph E − E, ξ) + (J(E − E h ), ξ), t0 µ0 (ηt , η) + (ξ, ∇ × η) = µ0 ((H I − H)t , η) + (Ph E − E, ∇ × η).
0 ∞ (ξt , ξ) − (∇ × η, ξ) +
764
QUN LIN AND JICHUN LI
Adding the above two equations together and using the definition of Ph , we obtain d 0 ∞ µ0 (s − ∞ )0 ( ||ξ(t)||20 + ||η(t)||20 ) + ||ξ(t)||20 dt 2 2 t0 ˜ = −(∇ × (H I − H), ξ) + (J(E − E h ), ξ) + µ0 ((H I − H)t , η) (41) =
3 (Err)i . i=1
The rest of the proof is omitted because it is very similar to the plasma case.
3.3. Lorentz medium. The Lorentzian two pole model can be represented by the governing equations [22]: (42) (43)
ˆ (E) 0 ∞ E t − ∇ × H + J µ0 H t + ∇ × E
= 0 in Ω × (0, T ), = 0 in Ω × (0, T ),
with the same boundary and initial conditions as those stated previously for plasma. ˆ to represent the polarization current for the Lorentz medium: We use J t ˆ ˆ (x, t; E) = β˜ ≡ J (44) J(E) e−δ(t−s) · sin(γ − α(t − s)) · E(x, s)ds, 0 t √ jγ −(δ+jα)(t−s) ˜ (45) e E(x, s)ds) ≡ Im(J(E)), j = −1, = Im(βe 0
2 where β˜ = (s − ∞ )0 ω13 / ω12 − ν4 , and Im(A) denotes the imaginary part of a general complex number A. Furthermore, in addition to the notation defined earlier, ω1 is the resonant frequency, ν is the damping coefficient, δ = ν2 , α = ω12 − δ 2 and ejγ = ωδ1 + j ωα1 . Note that ω1 > δ in real applications [17]. From (42)-(43), we can easily obtain the weak formulation: find the solution (E, H) ∈ [C 1 (0, T ; (L2 (Ω))3 ) ∩ C 0 (0, T ; H(curl; Ω))]2 such that (46) (47)
ˆ φ) 0 ∞ (E t , φ) − (∇ × H, φ) + (J(E), µ0 (H t , ψ) + (E, ∇ × ψ)
= =
0 ∀ φ ∈ (L2 (Ω))3 , 0 ∀ ψ ∈ H(curl; Ω)
for 0 < t ≤ T with the initial conditions (48)
E(x, 0) = E 0 (x)
and
H(x, 0) = H 0 (x).
The semi-discrete mixed finite element scheme for the Lorentz medium is formulated as follows: find E h ∈ Uh , H h ∈ Vh such that (49)
ˆ h ), φh ) 0 ∞ (E ht , φh ) − (∇ × H h , φh ) + (J(E
=
0 ∀ φh ∈ Uh ,
(50)
µ0 (H ht , ψ h )
+ (E , ∇ × ψ h )
=
0 ∀ ψ h ∈ Vh
h
for 0 < t ≤ T, subject to the initial conditions (51)
E h (x, 0) = Ph E 0 (x) and
H h (x, 0) = H I0 (x).
Theorem 3.3. Let (E(t), H(t)) and (E h (t), H h (t)) be the solutions of (46)-(47) and (49)-(50) at time t, respectively. Then there is a constant C = C(0 , µ0 , s , ∞ , ω1 , ν),
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS
765
independent of the mesh size h, such that µ0 ||(H I − H h )(t)||20 + 0 ∞ ||(Ph E − E h )(t)||20 t ≤ Ch2(k+1) [||H(t)||2k+2 + ||H t (t)||2k+1 ]dt,
(52)
0
where k is the degree of edge elements in the space V h . Proof. Subtracting (49)-(50) from (46)-(47) with φ = φh = ξ(t) and ψ = ψ h = η(t), respectively, we can obtain the error equations: 0 ∞ (ξt , ξ) − (∇ × η, ξ) = 0 ∞ ((Ph E − E)t , ξ) ˆ − E h ), ξ), − (∇ × (H I − H), ξ) − (J(E µ0 (ηt , η) + (ξ, ∇ × η) = µ0 ((H I − H)t , η) + (Ph E − E, ∇ × η). Adding the above two equations and using the definition of Ph , we obtain d µ0 0 ∞ ( ||η(t)||20 + ||ξ(t)||20 ) dt 2 2 ˆ = −(∇ × (H I − H), ξ) − (J(E − E h ), ξ) + µ0 ((H I − H)t , η) (53) =
3 (Err)i . i=1
The estimates of (Err)i follow the proof for the plasma case. The only different term is (Err)2 . ˆ and the definition of Ph , we obtain By the linearity of J ˆ ˆ (Err)2 = −(J(E − Ph E) + J(Ph E − E h ), ξ) = −(J(ξ), ξ) ≤
(54)
It is easy to see that
δ3 0 ∞ ||ξ(t)||20 +
1 ˆ (ξ)||2 . ||J 0 4δ3 0 ∞
t |ejγ e−(δ+jα)(t−s) ξ(x, s)ds|2 dΩ) Ω 0 t 1 2 ||ξ(s)||20 ds, ≤ β˜ 2δ 0
ˆ (ξ)||20 ≤ ||J (ξ)||20 = β˜2 ||J
which together with (54) gives
t β˜2 + ||ξ(s)||20 ds. (Err)2 ≤ 8δ3 0 ∞ δ 0 The rest of the proof follows exactly the same as the plasma case. δ3 0 ∞ ||ξ(t)||20
4. Global superconvergence To prove global superconvergence, we need some postprocessing operators introduced by Lin and Yan [29]. e) For each component wj , j = 1, 2, 3 of w ∈ Uh , we define Π2h wj |eˆ ∈ Qk,2k−1,2k−1 (ˆ such that (Π2h wj − wj )q = 0, ∀ q ∈ Qk,k−1,k−1 (ei ), i = 1, 2, 3, 4, ei
where eˆ =
4
i=1 ei .
766
QUN LIN AND JICHUN LI
ˆ 2h is defined for the function v ∈ Vh . For the Another postprocessing operator Π ˆ first component v1 of v, we define Π2h v1 |eˇ ∈ Q2k−1,k,k (ˇ e) such that ˆ 2h v1 − v1 )qdx = 0, ∀ q ∈ Pk−1 (li ), i = 1, · · · , 8, (Π li ˆ 2h v1 − v1 )qdxdy = 0, ∀ q ∈ Qk−1,k−2 (τi ), i = 1, · · · , 4, (Π τi ˆ 2h v1 − v1 )qdxdz = 0, ∀ q ∈ Qk−1,k−2 (τj ), j = 1, · · · , 4, (Π τj
ˆ 2h v1 − v1 )qdxdydz = 0, (Π
∀ q ∈ Qk−1,k−2,k−2 (ei ), i = 1, 2,
ei
where eˇ = e1 ∪e2 , li are edges parallel to the x-axis, τi , τj are surfaces perpendicular ˆ 2h can be defined similarly for the second to the z-axis or y-axis, respectively. Π and third components of v ∈ Vh . Lin and Yan [29] proved the following properties: Lemma 4.1. ˆ 2h v − v||0 ≤ Chk+1 ||v||k+1 , (i) ||Π2h w − w||0 ≤ Chk+1 ||w||k+1 , ||Π k+1 3 (Ω)) , ∀ w, v ∈ (H ˆ 2h v||0 ≤ C||v||0 , ∀ w ∈ Uh , v ∈ Vh , (ii) ||Π2h w||0 ≤ C||w||0 , ||Π ˆ ˆ 2h vI , ∀ w ∈ Uh , v ∈ Vh , (iii) Π2h w = Π2h Ph w, Π2h v = Π where Ph w ∈ Uh and vI ∈ Vh are the interpolations of w and v defined in Section 2. Using these postprocessing operators, we can achieve the following global superconvergence for all three dispersive media: Theorem 4.1. ˆ 2h H h − H||0 ||Π2h E h − E||0 + ||Π
≤ Chk+1 [||E||k+1 + ||H||k+1 + (
t
(||H||2k+2 + ||H t ||2k+1 )ds)1/2 ], 0
where k is the degree of edge elements in the space V h . Proof. By Lemma 4.1 and Theorems 3.1–3.3, we have ||Π2h E h − E||0
(55)
= ||Π2h (E h − Ph E) + (Π2h E − E)||0 ≤ C||E h − Ph E||0 + Chk+1 ||E||k+1 t ≤ Chk+1 [( (||H||2k+2 + ||H t ||2k+1 )ds)1/2 + ||E||k+1 ]. 0
Similarly, we have ˆ 2h H h − H||0 ||Π
(56)
ˆ 2h (H h − H I ) + (Π ˆ 2h H − H)||0 = ||Π ≤ C||H h − H I ||0 + Chk+1 ||H||k+1 t ≤ Chk+1 [( (||H||2k+2 + ||H t ||2k+1 )ds)1/2 + ||H||k+1 ], 0
which, along with (55), completes the proof.
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS
767
5. Extension to a standard finite element method In this section, we will extend the global superconvergence to a standard finite element method. For simplicity, we will restrict our discussion to the cold plasma model, generalizations to other models are similar. Instead of solving the coupled system (5)-(7) with both the electric and magnetic fields as unknowns, we eliminate H by taking the time derivative of (5) and using (6)-(8), to obtain the second order electric field equation (57)
2 0 E tt + ∇ × (µ−1 0 ∇ × E) + 0 ωp E − νJ(E) = 0,
with boundary condition (11) and initial conditions (58)
E(x, 0) = E 0 (x)
and
E t (x, 0) = E 1 (x),
where E 1 (x) = −1 0 ∇ × H 0 (x), which is obtained from (5), (8), and (12). Multiplying (57) by φ ∈ H0 (curl; Ω), and using integration by parts [35, (3.27)], we can easily obtain the weak formulation: find E(t) ∈ H0 (curl; Ω) such that (59)
2 0 (E tt , φ) + µ−1 0 (∇ × E, ∇ × φ) + 0 ωp (E, φ)
− ν(J(E), φ) = 0 φ ∈ H0 (curl; Ω),
subject to the initial conditions (58). Taking the boundary condition (11) into account, we define V 0h = {v h ∈ V h | n × v h = 0 on ∂Ω}. Then we can formulate a standard finite element scheme for (57) as follows: find E h (t) ∈ V 0h such that h h h 0 2 (60) 0 (E htt , φ)+µ−1 0 (∇×E , ∇×φ)+0 ωp (E , φ)−ν(J(E ), φ) = 0 ∀φ ∈ V h ,
subject to the initial conditions (61)
E h (x, 0) = E I0 (x),
E ht (x, 0) = E I1 (x),
where E I0 and E I1 are the interpolations of E 0 and E 1 defined in §2, respectively. Theorem 5.1. Let E(t) and E h (t) be the solutions of (59) and (60) at time t, respectively. Then there is a constant C = C(0 , µ0 , ωp , ν), independent of the mesh size h, such that I h I h 2 2 2 0 ||(E I − E h )t (t)||20 + µ−1 0 ||∇ × (E − E )(t)||0 + 0 ωp ||(E − E )(t)||0 t ≤ Ch2(k+1) [||E||2k+2 + (||E tt ||2k+1 + ||E t ||2k+2 + ||E||2k+1 )ds]. 0
Proof. Denote ξ(t) = (E − E )(t). Subtracting (60) from (59) with φ = ξt (t), and rearranging terms, we obtain the error equation I
h
2 0 (ξtt , ξt ) + µ−1 0 (∇ × ξ, ∇ × ξt ) + 0 ωp (ξ, ξt )
(62)
I = 0 ((E I − E)tt , ξt ) + µ−1 0 (∇ × (E − E), ∇ × ξt )
+ 0 ωp2 (E I − E, ξt ) + ν(J(E − E I ), ξt ) + ν(J(ξ), ξt ).
768
QUN LIN AND JICHUN LI
Integrating (62) with respect to t and using integration by parts and the fact that ξ(0) = ξt (0) = 0, we obtain
(63)
1 2 2 2 (0 ||ξt ||20 + µ−1 0 ||∇ × ξ||0 + 0 ωp ||ξ||0 ) 2 t I ≤ 0 ((E I − E)tt , ξt )ds + µ−1 0 (∇ × (E − E), ∇ × ξ) 0 t (∇ × (E I − E)t , ∇ × ξ)ds − µ−1 0 0 t 2 ((E I − E)t , ξ)ds + 0 ωp (E I − E, ξ) − 0 ωp2 0
t
t
(J(E − E I ), ξt )ds + ν
+ν 0
(J(ξ), ξt )ds =
7
0
(Err)i .
i=1
Now we have to estimate all (Err)i , i = 1, · · · , 7. Using Lemma 2.2 and the inequality (23), we have
t 0
δ1 0 ||ξt ||20 ds 0
Chk+1 ||E tt ||k+1 ||ξt ||0 ds 0
t
≤
t
((E I − E)tt , ξt )ds ≤ 0
= 0
(Err)1
t
+ 0
C0 h2(k+1) ||E tt ||2k+1 ds. 4δ1
Using Lemma 2.1 and the inequality (23) leads to I −1 k+1 (Err)2 = µ−1 ||E||k+2 ||∇ × ξ||0 0 (∇ × (E − E), ∇ × ξ) ≤ µ0 Ch 2 ≤ δ2 µ−1 0 ||∇ × ξ||0 +
Ch2(k+1) ||E||2k+2 , 4δ2 µ0
and (Err)3 = −µ−1 0
t
(∇×(E I −E)t , ∇×ξ)ds ≤ µ−1 0
0
t
≤
δ3 µ−1 0 ||∇
×
ξ||20 ds
0
Ch2(k+1) + 4δ3 µ0
t
Chk+1 ||E t ||k+2 ||∇×ξ||0 ds 0
t
||E t ||2k+2 ds. 0
Similarly, by Lemma 2.2 and the inequality (23), we obtain = 0 ωp2 (E I − E, ξ) ≤ 0 ωp2 Chk+1 ||E||k+1 ||ξ||0
(Err)4
≤ δ4 0 ωp2 ||ξ||20 +
0 ωp2 Ch2(k+1) ||E||2k+1 , 4δ4
and (Err)5
= −0 ωp2
t
0
t
≤ 0
t
((E I − E)t , ξ)ds ≤ 0 ωp2
0 ωp2 Ch2(k+1) δ5 0 ωp2 ||ξ||20 ds + 4δ5
Chk+1 ||E t ||k+1 ||ξ||0 ds
0 t
||E t ||2k+1 ds. 0
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS
769
Note that
t e−ν(t−s) (E − E I )(x, s)ds, ξt (x, t)) (J(E − E I ), ξt ) = 0 ωp2 ( 0 t Chk+1 ||E||k+1 ||ξt ||0 ds, ≤ 0 ωp2 0
from which we have (Err)6
t
(J(E − E ), ξt )ds ≤ I
= ν ≤
0 t
t
Chk+1 ||E||k+1 ||ξt ||0 ds
ν0 ωp2 t 0
δ6 0 ||ξt ||20 ds + 0
0 ωp4 t2 Ch2(k+1) 4δ6
t
||E||2k+1 ds. 0
Finally, using the inequality (23), we have t t t ν2 2 (Err)7 = ν (J(ξ), ξt )ds ≤ δ7 0 ||ξt ||0 ds + ||J (ξ)||20 ds 4δ7 0 0 0 0 t ν0 ωp4 t t ≤ δ7 0 ||ξt ||20 ds + ||ξ||20 ds, 8δ7 0 0 where in the last step we used the estimate (25). Substituting the above estimates for (Err)i into (63), choosing constants δ2 , δ4 < 1 and absorbing the first terms in (Err)2 and (Err)4 by the corresponding terms 2 on the left hand side of (63), we obtain
(64)
1 2 2 2 (0 ||ξt ||20 + µ−1 0 ||∇ × ξ||0 + 0 ωp ||ξ||0 ) 2 t 2 2 2 ≤ C1 (0 ||ξt ||20 + µ−1 0 ||∇ × ξ||0 + 0 ωp ||ξ||0 )ds 0 t (||E tt ||2k+1 + ||E t ||2k+2 + ||E||2k+1 )ds + C2 h2(k+1) 0
+ C3 h
2(k+1)
||E||2k+2 ,
where we absorbed the explicit dependence of physical parameters into the generic constants C1 , C2 and C3 . Our proof concludes by using the Gronwall inequality (Lemma 2.3) to (64). ˆ 2h defined in the last section, we can easUsing the postprocessing operator Π ily achieve the following global superconvergence for the standard finite element method: Theorem 5.2. ˆ 2h E h − E||0 + ||(Π ˆ 2h E h − E)t ||0 ||Π
≤ Chk+1 [||E||k+2 + ||E t ||k+1 + (
t
(||E tt ||2k+1 + ||E t ||2k+2 + ||E||2k+1 )ds)1/2 ], 0
where k is the degree of edge elements in the space V h . Remark 5.1. When Ω is not a cubic domain, global superconvergence of order 1 O(hk+ 2 ) can be achieved on the almost cubic meshes by easily extending the results of Lin and Yan [29, p. 175].
770
QUN LIN AND JICHUN LI
Acknowledgements Jichun Li thanks Dr. Miguel Visbal at the U.S. Air Force Research Laboratory (AFRL) for inspiring him to study the interesting dispersive media problems during his summer research at AFRL. The authors wish to thank the referees for many constructive comments that improved the paper. References [1] M. Ainsworth, P. Davies, D. Duncan, P. Martin and B. Rynne, Topics in computational wave propagation, Lecture Notes in Computational Science and Engineering 31, Springer-Verlag, Berlin, 2003. MR2032865 (2004h:65008) [2] A. Bossavit, Computational Electromagnetism, Academic Press, San Diego, 1998. MR1488417 (99m:78001) [3] C. Carstensen, S. Funken, W. Hackbusch, Ronald H.W. Hoppe and P. Monk (eds.), Computational Electromagnetics, Lecture Notes in Computational Science and Engineering 28, Springer-Verlag, Berlin, 2003. MR1995092 (2004b:78001) [4] C.M. Chen and Y.Q. Huang, High Accuracy Theory of Finite Element Methods (in Chinese), Hunan Science Press, China, 1995. [5] M.-H. Chen, B. Cockburn and F. Reitich, High-order RKDG methods for computational electromagnetics, J. Sci. Comp. 22/23 (2005) 205-226. MR2142195 (2005m:65208) [6] Q. Chen, M. Katsurai and P.H. Aoyagi, An FDTD formulation for dispersive media using a current density, IEEE Trans. Antennas Propagat. 46 (1998) 1739-1746. [7] Z. Chen, Q. Du and J. Zou, Finite element methods with matching and nonmatching meshes for Maxwell equations with discontinuous coefficients, SIAM J. Numer. Anal. 37 (2000) 15421570. MR1759906 (2001h:78044) [8] P. Ciarlet, Jr. and J. Zou, Fully discrete finite element approaches for time-dependent Maxwell’s equations, Numer. Math. 82 (1999) 193-219. MR1685459 (2000c:65083) [9] B. Cockburn, F. Li and C.-W. Shu, Locally divergence-free discontinuous Galerkin methods for the Maxwell equations, J. Comp. Phys. 194 (2004) 588-610. MR2034859 (2004j:78024) [10] G.C. Cohen, E. Heikkola, P. Joly and P. Neittaanm¨ aki (eds.), Mathematical and numerical aspects of wave propagation — WAVES 2003, Springer-Verlag, Berlin, 2003. MR2077970 (2005b:00020) [11] 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 Propagat. 45 (1997) 392-400. [12] L. Demkowicz, Fully automatic hp-adaptivity for Maxwell’s equations, Comput. Methods Appl. Mech. Engrg. 194 (2005) 605-624. MR2105184 (2005m:65296) [13] L. Demkowicz, P. Monk, Ch. Schwab and L. Vardapetyan, Maxwell eigenvalues and discrete compactness in two dimensions, Comput. Math. Appl. 40 (2000) 589-605. MR1772657 (2002e:65166) [14] R.E. Ewing, Y. Lin, T. Sun, J. Wang, and S. Zhang, Sharp L2 -error estimates and superconvergence of mixed finite element methods for non-Fickian flows in porous media, SIAM J. Numer. Anal. 40 (2002) 1538-1560. MR1951906 (2003m:65177) [15] R. Hiptmair, Finite elements in computational electromagnetism, Acta Numerica 11 (2002) 237-339. MR2009375 (2004k:78028) [16] P. Houston, I. Perugia and D. Sch¨ otzau, Energy norm a posteriori error estimation for mixed discontinuous Galerkin approximations of the Maxwell operator, Comput. Methods Appl. Mech. Engrg. 194 (2005) 499-510. MR2105178 (2005h:78027) [17] D. Jiao and J.-M. Jin, Time-domain finite-element modeling of dispersive media, IEEE Microwave and Wireless Components Letters 11 (2001) 220-222. [18] J. Jin, The Finite Element Method in Electromagnetics, 2nd ed., John Wiley & Sons, Inc., New York, 2002. MR1903357 (2004b:78019) [19] M. Krizek and P. Neittaanmaki, Bibliography on Superconvergence, in “Finite Element Methods”, Lecture Notes in Pure and Applied Math., Vol. 196, 1998, pp. 315-348. MR1602730 [20] J. Li, Error analysis of finite element methods for 3-D Maxwell’s equations in dispersive media, Journal of Computational and Applied Mathematics, 188 (2006) 107-120. MR2192630 (2006m:65279)
SUPERCONVERGENCE ANALYSIS FOR MAXWELL’S EQUATIONS
771
[21] J. Li and Y. Chen, 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. MR2229840 [22] J. Li, Error analysis of fully discrete mixed finite element schemes for 3-D Maxwell’s equations in dispersive media, preprint, Sept. 2005. [23] J. Li and M.F. Wheeler, Uniform convergence and superconvergence of mixed finite element methods on anisotropically refined grids, SIAM J. Numer. Anal. 38 (2000) 770-798. MR1781203 (2001f:65137) [24] Q. Lin, J. Li and A. Zhou, rectangle test for Ciarlet-Raviart scheme, Prof. of Sys. Sci. & Sys. Engrg., Great Wall (H.K.) Culture Publish Co., 1991, pp. 230-231. [25] Q. Lin, J. Li and A. Zhou, A rectangle test for the Stokes equations, Prof. of Sys. Sci. & Sys. Engrg., Great Wall (H.K.) Culture Publish Co., 1991, pp. 240-241. [26] Q. Lin and J.F. Lin, High accuracy approximation of mixed finite element for 2-D Maxwell equations. (Chinese) Acta Math. Sci. Ser. A Chin. Ed. 23 (2003), no. 4, 499-503. MR2090553 (2005i:65145) [27] Q. Lin and J. Xu, Linear finite elements with high accuracy, J. Comp. Math. 3 (1985) 115-133. MR854355 (87k:65141) [28] Q. Lin and N. Yan, Superconvergence of mixed element methods for Maxwell’s equations (in Chinese), Gongcheng Shuxue Xuebao 13 (1996), suppl., 1-10. MR1437477 (98c:65186) [29] Q. Lin and N. Yan, Global superconvergence for Maxwell’s equations, Math. Comp. 69 (1999) 159-176. MR1654029 (2000i:65131) [30] Q. Lin and N. Yan, The Construction and Analysis of High Accurate Finite Element Methods (in Chinese), Hebei University Press, Hebei, China, 1996. [31] Q. Lin, N. Yan and A. Zhou, A rectangle test for interpolated finite elements, Prof. of Sys. Sci. & Sys. Engrg., Great Wall (H.K.) Culture Publish Co., 1991, pp. 217-229. [32] T. Lu, P. Zhang and W. Cai, Discontinuous Galerkin methods for dispersive and lossy Maxwell’s equations and PML boundary conditions, J. Comp. Phys. 200 (2004) 549-580. MR2095277 (2005e:78025) [33] P. Monk, An analysis of N´ed´ elec’s method for the special discretization of Maxwell’s equations, J. Comp. Appl. Math. 47 (1993) 101-121. MR1226366 (94g:65105) [34] P. Monk, Superconvergence of finite element approximations to Maxwell’s equations, Numer. Methods Partial Differential Equations 10 (1994) 793-812. MR1298123 (95h:65090) [35] P. Monk, Finite element methods for Maxwell’s equations, Oxford University Press, 2003. MR2059447 (2005d:65003) [36] J.-C. N´ ed´ elec, Mixed finite elements in R3 , Numer. Math. 35 (1980) 315-341. MR592160 (81k:65125) [37] A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equations, Springer-Verlag, Berlin, 1994. MR1299729 (95i:65005) [38] A. Taflove and C. Hagness, Computational Electrodynamics: The Finite-Difference TimeDomain Method, 2nd ed., Artech House, Norwood, MA, 2000. MR1338377 (96f:78001) [39] L.B. Wahlbin, Superconvergence in Galerkin Finite Element Methods, Springer, Berlin, 1995. MR1439050 (98j:65083) [40] Z. Zhang, N. Yan and T. Sun, Superconvergent derivative recovery for the intermediate finite element family of the second order, IMA Journal of Numerical Analysis 21 (2001) 643-665. MR1844362 (2002f:65171) LSEC, ICMSEC, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100080, China E-mail address: [email protected] Department of Mathematical Sciences, University of Nevada, Las Vegas, 4505 Maryland Parkway, Box 454020, Las Vegas, Nevada 89154-4020 E-mail address: [email protected]
