CONVERGENCE ANALYSIS OF A COVOLUME ... - Semantic Scholar
CONVERGENCE ANALYSIS OF A COVOLUME SCHEME FOR MAXWELL'S EQUATIONS IN THREE DIMENSIONS R.A. NICOLAIDES AND D-Q WANG Abstract. This paper contains error estimates for covolume discretizations of
Maxwell's equations in three space dimensions. Several estimates are proved. First, an estimate for a semi-discrete scheme is given. Second, the estimate is extended to cover the classical interlaced time marching technique. Third, some of our unstructured mesh results are specialized to rectangular meshes, both uniform and nonuniform. By means of some additional analysis it is shown that the spatial convergence rate is one order higher than for the unstructured case.
1. Introduction Staggered mesh schemes for the numerical solution of Maxwell's equations go back as far as [18]. Over the years this scheme has seen a number of generalizations intended to enhance its usefulness. Mostly, these generalizations are aimed at increasing the geometric complexity that can be handled. Thus, in [3] a tensor formulation was given. This permits the method to be extended from rectangular meshes to meshes de ned by curvilinear coordinates. Two dimensional generalizations using quadrilateral meshes are given in [6] and [8]. For three dimensions, hexahedral mesh formulations were proposed in [7] and [9]. These formulations use interpolation to obtain values for nonbasic eld components{ those which, in the rectangular situation, lie along the primal or dual mesh edges. Another approach is the \control path" method of [4], [5]. The control path method uses the classical nite di erence approach of modifying the nite di erence stencil near the boundary of the domain. A natural generalization of the standard staggered mesh scheme to tetrahedral meshes was given in [10]. A similar technique was used independently in [14] to solve the Navier-Stokes equations. This \covolume" approach does not use any interpolations and is a very natural way to generalize the original rectangular staggered mesh approach. A characteristic feature of covolume schemes is the use of Voronoi-Delaunay mesh pairs to replace the rectangular staggered mesh arrangement. Very little rigorous analysis of staggered mesh schemes for electromagnetics is available. We know only of the rectangular mesh analysis of [11] and [12]. For the incompressible Navier{Stokes equations results may be found in [14], [15] and [16]. Our goal in this paper is to provide a rigorous proof of convergence for the covolume discretization applied to the interior problem for Maxwell's equations. Date : September 15, 1995 and, in revised form, February 20, 1996. 1991 Mathematics Subject Classi cation. Primary 65N30, 65N15, 35L50 . Key words and phrases. Maxwell's equations, covolume schemes, unstructured mesh, three dimensions, convergence . 1
R.A. NICOLAIDES AND D-Q WANG
2
The corresponding analysis for the exterior problem involves the application of radiation conditions at a nite distance from the scatterer and will be given in a subsequent report. The contents of the paper are as follows: in the rst section we state Maxwell's equations for bounded domains with perfect conductor boundary conditions. Following that is a section on Voronoi-Delaunay mesh notations and properties. The main result of section 3 is an error estimate for a semi{discrete covolume approximation to Maxwell's equations. In practice, a kind of leapfrog scheme is used to time march the discrete equations. Section 4 contains an error estimate for the resulting fully discrete approximation to Maxwell's equations. In the last section we specialize our results to the rectangular case. The rate of convergence is shown to be one order higher than for an arbitrary triangulation. 2. Maxwell's Equations Let be a bounded domain in R3 with boundary ? and unit outward normal n. Let the constants and denote, respectively, the electric and magnetic permeabilities of the medium occupying . Then if E(x; t) and H(x; t) denote, respectively, the electric and magnetic elds, Maxwell's equations [2] are: (2.1) Et ? curl H = J in
(0; T ); (2.2)
Ht + curl E = 0
in
(0; T );
(2.3) div( E) = and div( H) = 0 where J = J(x; t) is a known applied current and (x; t) is a charge density. We shall assume a perfect conductor boundary condition so that (2.4) En=0 on ? (0; T ): In addition, initial conditions are prescribed so that (2.5) E(x; 0) = E0(x) and H(x; 0) = H0(x) 8x 2
where E0 and H0 are given functions satisfying div( E0 ) = (x; 0); div( H0 ) = 0 and where it is assumed that @ = div J: (2.6) @t
Let Lp (0; T ; X ) denote the set of all strongly measurable functions u(t; ) from [0; T ] into the Hilbert space X such that
ZT 0
jju(t)jjpX dt < 1
1 p < 1:
are in Lp (0; T ; X ). Also We say u 2 W 1;p (0; T ; X ) if and only if both u and @u @t C m (0; T ; X ) denotes the space of m times continuously di erentiable functions from [0,T] into X . By H (curl; ) and H (div; ) we mean the Hilbert spaces de ned by H (curl; ) : = fv 2 (L2 ( ))3 ; curl v 2 (L2 ( ))3 g; H (div; ) : = fv 2 (L2 ( ))3 ; div v 2 L2 ( )g:
ANALYSIS OF A COVOLUME SCHEME FOR MAXWELL'S EQUATIONS
3
We shall assume the existence of (E; H) to (2.1)-(2.5) such that E; H 2 C 1 (0; T ; (L2( ))3 ) \ C 0(0; T ; H (curl; ) \ H (div; )): For this, it is sucient that J 2 C 0 (0; T ; (L2( ))3 ), 2 C 0 (0; T ; L2( )) and @ 2 L2(0; T ; H ?1( )), see [2]. @t
H H H
n
E
E H
H Figure 1 3. Mesh Notations and Discrete Vector Fields In this section some basic properties of dual mesh systems are introduced. These will lead to a detailed formulation and analysis of our numerical schemes in the following sections. Assume that the polyhedral domain has a primal family of nite element style tetrahedral partitions, parameterized by the maximum side length which is generically denoted by h. We will assume that the ratio of radii of circumscribing spheres and inscribed spheres of all the individual tetrahedra are uniformly bounded above and below as h approaches 0. A dual mesh is formed by connecting adjacent tetrahedral circumcenters and, in the case of tetrahedron with a face on a boundary, by connecting their circumcenters with those of their boundary faces. By elementary geometry these dual edges are perpendicular to the associated tetrahedral faces. These connections also form the edges of a set of polyhedra. It follows from elementary geometry that the edges of tetrahedra are perpendicular to and in one-to-one correspondence with the faces of dual polyhedra or \covolumes". The reciprocal orthogonality between edges and faces is the key to the results which follow. The N nodes of the tetrahedral mesh are assumed to be numbered sequentially in some convenient way, and likewise the T nodes of dual mesh. Similarly the F faces (edges) and M edges (faces) of the primal (dual) mesh are sequentially numbered. The individual tetrahedra, faces, edges and nodes of the primal mesh are denoted by i , j , k and l respectively. Those of the dual mesh are denoted by primed quantities such as j0 . A direction is assigned to each primal edge by the rule that the positive direction is from low to high node number. The dual edges are directed by the corresponding rule. We also denote F1 the number of tetrahedral interior
4
R.A. NICOLAIDES AND D-Q WANG
faces (or dual edges) and M1 the number of tetrahedral interior edges (or dual faces). Let sj denote the area of j and h0 j the length of j0 . In RF1 , where F1 denotes the number of interior primal faces, we will introduce the inner product (; )W de ned by X (u; v)W := (3.1) uj vj sj h0j = (Su; D0v) = (D0 u; Sv) 2 j
and denote the resulting inner product space by U and the associated norm by 1 jjujjW := (u; u)W2 : (3.2) In (3.1), (; ) denotes the standard Euclidean inner product, S := diag(sj ), D0 := diag(h0j ) and W := SD0 are F1 F1 invertible diagonal matrices. The norm de ned by (3.2) is clearly three times a discrete L2 norm. Similarly, we will introduce an inner product in RM , where M denotes the number of primal edges : X (3.3) uj vj s0j hj = (S 0 u; Dv) = (Du; S 0 v) (u; v)W 0 := 0 2
j
and denote the inner product space by U 0 . The associated norm is 1 (3.4) jjujjW 0 := (u; u)W2 0 : The notations in(3.3) and (3.4) correspond to those of (3.1) and (3.2). For example, s0 j denotes the area of dual face 0j . For each primal face i a discrete circulation is de ned by X ~ (3.5) uj hj : (Cu) := i
2@ j
i
Similarly, for each interior covolume face 0i the discrete circulation is X ~0 (3.6) (C 0 u)0 := uj hj : 0 2@0 A tilde on hj or h0j means that the quantity is to be taken with a negative sign if the dual edge is directed against the positive sense of description of @i or @0i i
j
i
respectively, and with a positive sign otherwise. The linear operators C and C 0 map from RM to RF1 and RF1 to RM1 , respectively. For each strictly interior dual edge j0 we can form a vector whose ith component is the sign of the orientation of the edge relative to the orientation of the ith strictly interior dual face. From these vectors we obtain the F1 M1 matrix G de ned as follows: 8 < 1 if j0 is oriented positively along 0i (G)ji := : ?1 if j0 is oriented negatively along 0i 0 if j0 does not meet 0i Let w 2 RM denote the vector whose kth component is the value assigned on the kth primal edge. Denote w1 2 RM1 the restriction of w to the interior primal edges. So wj? = 0 means that the components of w 2 RM related to the boundary edges are zero. A direct computation shows: Cw = GDw1 :
ANALYSIS OF A COVOLUME SCHEME FOR MAXWELL'S EQUATIONS
5
Also it can be veri ed that if v 2 RF1
C 0 v = GT D0 v:
From these two identities we can prove Lemma 3.1. With the above de nitions of w,v and w1 with wj? = 0, (Cw; D0 v) = (C 0 v; Dw1 ): Proof. : This is proved as follows: (C 0 v; Dw1 ) = (GT D0 v; Dw1 ) = (D0 v; GD0 w1 ) = (D0 v; Cw): This lemma provides a discrete analog of the integration formula
Z
curl E H dx =
Z
curl H E dx
which holds when E n = 0 on ?. For each tetrahedron i a discrete ux is de ned by X uj s~j ; 8u 2 RF1 : (Du)i := 2@ j
i
By s~j we mean sj negatively signed if the corresponding velocity component is directed towards the inside of i and positively signed otherwise. In addition to G we will introduce another matrix B1 of dimension F1 T 8 < 1 if j is oriented positively along @ i (B1 )ji := : ?1 if j is oriented negatively along @ i 0 if j does not meet @ i It can be checked directly that D = B1T S: Using B1 we will de ne the di erence operator P by P := D0 ?1 B1 ; 8 2 RT : Now we have
Lemma 3.2.
(3.7) B1T C = 0; and (3.8) (u; P )W = (Du; ); 8 2 RT ; u 2 RF1 : Proof. : For (3.7) see [17] Theorem 1. For (3.8), we have (u; P )W = (SD0 u; D0 ?1 B1 ) = (B1T Su; ) = (Du; ): We remark that (3.7)Zand (3.8) provide the Z discrete analogs of the identity div(curl u) = 0 and Green's formula u grad dx = (div u) dx for 2 H01 ( ); u 2 H (div; ),
respectively. For additional details and relations between other discrete operators in VoronoiDelaunay meshes, see [13] and [17].
R.A. NICOLAIDES AND D-Q WANG
6
4. Semi-Discrete Maxwell's Equations In this section the covolume method is used to obtain a spatial discretization of Maxwell's equations (2.1)-(2.5) from which we will obtain a formulation of a semi-discrete form of Maxwell's equations. A basic error estimate is given in Theorem 4.1. First we introduce, for general eld A, its \face averages" Af 2 RF and its \edge averages " Ae 2 RM as follows (refer to Figure 1): Z A := 1 A n ds; f
si Z A0f := s10 0 A n ds; i Z Ae := h1 A t d; j Z 1 0 Ae := h0 0 A t d; j i
i
j
j
where n and t denote the unit normal vector to the face i (or 0i ) and the unit tangent vector to the edge j (or j0 ), respectively. Error functions for primal edges and faces will be denoted by A := A ? Af ; A := A ? Ae ; A := Af ? Ae : Error functions for dual edges and faces are de ned similarly. As shown in Figure 1, for each tetrahedron i we will use the normal components of the magnetic eld H to its faces and the tangential components of the electric eld E in the directions of its edges. Now integrate both sides of (2.1) over the co-face 0j to obtain Z d ( E ) f j 0 0 0 s (4.1) ? (C H ) = J(x; t) dx: j
e
dt
0
j
j
Here, (Ef )j denotes the average of E nj over the face 0j where the nj is the unit vector in the direction of j , and (C 0 He )0 is the discrete circulation around the face 0j . Similarly, from (2.2), (4.2) s d(Hf )i + (CE ) = 0 j
i
dt
e
i
where (Hf )i denotes the average of H ni where ni is the unit normal to the face i and (CEe ) is the discrete circulation around the face i . Let E and H denote vectors of components in RM1 and RF1 respectively. Then (4.1) and (4.2) (which are exact) suggest the approximations ( where it is implied that components of E associated with boundary edges are zero, i.e. E j? = 0). (4.3) S 0 dE ? C 0 H = J;~ i
(4.4)
dt S dH dt + CE = 0
ANALYSIS OF A COVOLUME SCHEME FOR MAXWELL'S EQUATIONS
7
where J~ 2 RM1 with components given by the right hand of (4.1). Since S and S 0 are invertible, (4.3)-(4.4) is a system of linear ordinary di erential equations, and the existence and uniqueness of a solution follow from well known results. From (4.4) and (3.7) in lemma 3.2 we obtain
dtd (DH ) = dtd (B1T SH ) = ?B1T CE = 0: This shows the sense in which divH = 0 is satis ed at the discrete level in the covolume scheme. By subtracting (4.3) from (4.1) we have
(4.5)
S 0 ddt (E ? Ef0 ) ? C 0 (H ? He0 ) = 0:
Similarly, from (4.2) and (4.4) (4.6)
S ddt (H ? Hf ) + C (E ? Ee ) = 0
where (E; H) denotes the exact solution of (2.1)-(2.5) and where, by the boundary condition (2.4), (4.7) E j? = 0: Note that error in the magnetic eld H satis es the discrete solenoidal condition (4.8) D_H = 0: Multiplying (4.5) by DE and (4.6) by D0 H , and adding we obtain (4.9) (_0E ; E )W 0 + (_H ; H0 )W = (C 0 H0 ; DE ) ? (CE ; D0 H0 ) where the dots denote time di erentiations. The main result in this section is the following theorem: Theorem 4.1. Denote by (E; H ) the solution of (4.3) and (4.4) and by (E; H) 2 W 1;1 (0; T ; (W 1;p( ))3 )2 the solution of (2.1)-(2.5) with p > 2. Then we have the estimate max (jj(E ?Ee )(t)jjW 0 +jj(H ?He )(t)jjW ) Kh(jjE_ jjL1 (0;T ;(W 1 ( ))3 ) +jjH_ jjL1 (0;T ;(W 1 0tT
To prove this theorem we use (4.9) to get
;p
(_E ; E )W 0 + (_H0 ; H0 )W = (E ? 0E ) ; E + (H0 ? H ) ; H0 W W + (C 0 H0 ; DE ) ? (CE ; D0 H0 ): Since E nj? = 0, the components of E restricted to the boundary is 0. So
lemma 3.1 is applied to yield (C 0 H0 ; DE ) ? (CE ; D0 H0 ) = 0: By E ? 0E = Ef ? Ee ; H0 ? H = Hf ? He ;
;p
( ))3 ) ):
R.A. NICOLAIDES AND D-Q WANG
8
we obtain from (4.9) (4.10) 1 d (jj() 12 0 jj2 + jj() 12 jj2 ) = ((H ? H ) ; 0 ) + ((E ? E ) ; ) 0 E W0 f e f e EW H W HW 2 dt where () denotes time di erentiations. Integrating (4.10) from 0 to T and using Cauchy's inequality we obtain (jj() 12 H0 (T )jj2W + jj() 12 E (T )jj2W 0 )
2
ZT 0
(jj() 12 (He ? Hf ) jjW jj() 12 H0 (s)jjW + jj() 12 (Ee ? Ef ) jjW 0 jj() 12 E (s)jjW 0 ) ds:
Let t be such that 1 1 2 H (t)jjW + jj() 2 E (t)jjW 0 ): (jj() 12 H0 (t )jjW + jj() 12 E (t )jjW 0 ) = 0max ( jj ( ) tT Then (jj 12 H0 (t )jjW + jj 12 E (t )jjW 0 )2 2(jj 21 H0 (T )jj2W + jj 21 E (T )jj2W 0 )
4 4 So it follows that (4.11)
Z t Z0 T 0
(jj 12 (He ? Hf ) jjW jj 12 H0 (s)jjW + jj 12 (Ee ? Ef ) jjW 0 jj 21 E (s)jjW 0 ) ds
(jj 21 (He ? Hf ) jjW jj 12 H0 (s)jjW + jj 21 (Ee ? Ef ) jjW 0 jj() 12 E (s)jjW 0 ) ds:
1 1 0 2 (t)jjW + jj() 2 E (t)jjW 0 ) 4 max ( jj ( ) H 0tT
ZT 0
(jj() 12 (He ? Hf ) jjW + jj() 12 (Ee ? Ef ) jjW 0 ) ds:
Theorem 4.1 will follow from the next lemma: Lemma 4.2. Assume that E_ ; H_ 2 (W 1;p( ))3 , p > 2. Then there exists a generic constant K , such that jj(He ? Hf ) jjW KhjH_ j(W 1 ( ))3 : jj(Ee ? Ef ) jjW 0 KhjE_ j(W 1 ( ))3 : where si ); max( s0i ) K K 0 max max ( 0 2 i h2 i ;p
;p
hi
and K 0 is a constant.
i
Proof. For the primal face i ,
Z Z (He ? Hf ) = h10 0 H_ n d l ? s1 H_ n dx: i
i
i
i
By the Sobolev embedding theorem, V : H ! (H_ e ? H_ f ) de nes a bounded linear functional on (W 1;p (j ))3 . For constant H_ this linear functional vanishes in the union of two tetrahedra that share the same face i and consequently j(He ? Hf ) j K ( i [ l )jH_ jW 1 ( [ )3 where K ( i [ l ) is a constant. In order to estimate K ( i [ l ), we use a standard scale change argument as follows : ;p
i
l
ANALYSIS OF A COVOLUME SCHEME FOR MAXWELL'S EQUATIONS
9
Let the primal face i = i \ l be in the xy plane so that i0 is parallel to the z axis , changing the coordinates by
0 = A xy0 ; z = h1 z 0 where A is 2 2 matrix, it results that K ( i [ l ) depends on the quantity 1 0 1 h ? ? 1 ? 1 1 max(jjA jj2 (jdet Ajh1 ) ; 1 ) (det A) where 1p + p10 = 1. Let h0i denote the length of the co-edge and si the area of the primal face. A further calculation then shows that jdet Aj = c1 si , h1 = c2 h0i so that 1?1 p p (AT A) 2 4 T A) 1 det ( A s min ? ? 1 ? 1 i jjA jj2 (jdet Ajh1 ) = c3 c = c 3 4 1 1 1 (si h0i ) (si h0i ) h0 i where ci , i = 1; ; 4 are independent of h0i and si , and min is the least eigenvalue of the positive de nite matrix AT A. Collecting these results, K ( i [ l ) is bounded by s 12 ? 1 0 10 i ;hi : K ( i [ l ) max 1 1 i h 0 i si x y
p
p
p
p
p
p
p
p
p
p
p
p
Now we have
jj(He ? Hf ) jj2W =
F X i=1
K1 K2
si h0i j(He ? Hf ) j2 F X i=1 F X
1? 2 0 20 h i jH_ j2 2 ; 2 i i (W 1 ( [ ))3 h0 i si
si h0 max si
p
p
;p
p
p
i
l
h5? 6 jH_ j2(W 1 ( [ ))3 p
;p
i
l
i=1 where K2 depends only on max( s0i2 ). By Holder's inequality and since h3 F hi 0
K measure( ) we obtain
jj(He ? Hf ) jj2
W
K2
h5? p6
2 X X F F p _ jHj(W 1 ( [ ))3 1 p
;p
i
l
i=1 p?2 6 5 ? p = K3 h F p jH_ j2(W 1;p ( ))3 Kh2 jH_ j2(W 1;p ( ))3 : estimate for jj(Ee ? Ef ) jjW 0 is similar.
?2
p
p
i=1
The proof of the From the estimate above we know the constant K depends on the quantity 0
( si )): max(max ( s0i2 ); max i h2 i
hi
i
10
R.A. NICOLAIDES AND D-Q WANG
Theorem 4.1 now follows from (4.11) and lemma 4.2. 5. The Fully Discrete problem There are many possible time stepping methods that can be applied to (4.3) and (4.4). We will discuss a leapfrog scheme which is very popular in computational electromagnetics (see [18]). In this scheme we approximate E(t) at times tn = n 4 t, 0 n < 1 with a vector fE n g1 n=0 , and H(t) at time tn+ 21 with a vector 1 1 1 n + 2 2 fH gn=0 . The initial value H can be computed using, for example, a Taylor expansion and the given equations (2.1)-(2.2). Given (E n ; Hn+ 21 )n0 , the next approximation (E n+1 ; Hn+ 32 ) is obtained by solving the equations (5.1) S 0 (E n+1 ? E n ) ? 4tC 0 Hn+ 21 = J~n+ 21 (5.2) where
S (Hn+ 32 ? Hn+ 12 ) + 4tC E n+1 = 0 Z (n+1)4t J~n+ 21 = J~ dt: n4t
(5.1), (5.2) is an explicit scheme, so the existence and uniqueness of a solution are apparent. Using the error functions de ned in the last section, we can rewrite (5.1)-(5.2) as 1 (5.3) S 0 (nE+1 ? nE ) = (4t)C 0 Hn+ 2 + Gn ; 3 1 S (nH+ 2 ? Hn+ 2 ) = ?(4t)CEn+1 + Gn : (5.4) By a direct computation, Gn and G n are given by 1 (5.5) Gn = J~n+ 12 ? S 0 (Efn+1 ? Efn ) + 4tC 0 Hen+ 2 ; 1 3 (5.6) G n = S (Hfn+ 2 ? Hfn+ 2 ) ? 4tCEen+1: Our main estimate for the fully discrete scheme is then given in the following theorem: Theorem 5.1. Let (E n ; Hn+ 21 )Nn?01 denote the solution of (5.1)-(5.2), and let (E; H) 2 (H 1 (0; T ; (W 1;p( ))3 ))2 denote the solution of (2.1)-(2.4), p > 2. Under the stability condition (5.7) c 4 t < min(hij )
pM M
3
3 22
where c = ()? 21 is the speed of the light in the medium, M2 is the maximum of the ratios of the maximum to minimum side-lengths over the union of adjacent tetrahedra and M3 is the maximum number of edges over all co-faces, we have the following error estimate for the fully discrete scheme (5.1)-(5.2) (5.8) 1 max (jjE i ? Eei jjW 0 + jjHi+ 12 ? Hei+ 2 jjW ) Kh(jjEjjH 1 (0;T ;(W 1 ( ))3 ) + jjHjjH 1 (0;T ;(W 1 ( ))3 ) ): 0iN ?1 ;p
;p
ANALYSIS OF A COVOLUME SCHEME FOR MAXWELL'S EQUATIONS
11
Proof. Multiplying (5.3) by D(En + En+1 ) and (5.4) by D0 (Hn+ 2 + Hn+ 2 ) respectively and adding all these equations from n = 0; 1; ; N ? 1 gives 1
3
jjEN jj2W 0 +jjHN ? 2 jj2W = 4t(C 0 HN ? 2 ; DEN )+ 1
+
1
NX ?1
NX ?2
i=0
i=0
((Gi ; D(Ei + Ei+1 )) +
NX ?1 i=0
NX ?2
(Ei+1 ?Ei ; Ei )W 0 +
i=0
(Hi+ 2 ?Hi+ 2 ; Hi+ 2 )W 3
(G i ; D0 (Hi+ 2 + Hi+ 2 ))): 3
1
The proof has three steps : (i): First, 1 1 4t(C 0 HN ? 2 ; DEN ) = 4t(C 0 (SD0 )? 21 ((SD0 ) 12 HN ? 2 ); (DS 0 ?1 ) 12 ((DS 0 ) 21 EN )) 1 4tjj(SD0)? 21 C 0 (DS 0 ?1 ) 12 jj2 jjHN ? 2 jjW jjEN jjW 0 : From algebra, jj(SD0 )? 21 C 0 (DS 0 ?1 ) 12 jj2 is the largest singular value of the matrix and by Gershgorin's theorem 3 p 2 jj(SD0 )? 12 C 0 (DS 0 ?1 ) 21 jj2 max( maxij (hij )52 )(2 M3 ) minij (hij ) where maximum above is taken over all of the union of adjacent tetrahedra. Combining the results gives pM M 23 1 3 2 jj N ? 12 jj jj N jj 0 4t(C 0HN ? 2 ; DEN ) 4t 2min( hij ) H W E W
pM M 32 3 2 c(jj N jj2 + jj N ? 21 jj2 ) 4t min( E W0 H W hij ) 1 < (jjEN jj2W 0 + jjHN ? 2 jj2W ):
(ii) Ei+1 ? Ei can be estimated in the following way. We have, by integrating in time jEi+1 ? Ei j jj_E jjL1 (i4tt(i+1)4t): So M1 X s0k hk j(Ei+1 )k ? (Ei )k j2 k=1 Z (i+1)4t M1 X 0 sk hk ( j(_E )k j ds)2 i 4 t k=1 Z (i+1)4t X M1 4t s0k hk j(_E )k j2 ds i4t k=1 Z (i+1)4t = 4t jj_E jj2W 0 ds:
jjEi+1 ? Ei jj2W 0 =
By lemma 4.2
i4t
p
jjEi+1 ? Ei jjW 0 K1h 4tjjE_ jjL2 ((i4t;(i+1)4t);(W 1
;p
( ))3 ) :
1
1
R.A. NICOLAIDES AND D-Q WANG
12
Similarly,
p
jjHi+ 2 ? Hi+ 2 jjW K1 h 4tjjH_ jjL2 ((i4t;(i+1)4t);(W 1 ( ))3 ) : (iii) We have from the de nition of jj jjW 0 (Gi ; D(Ei + Ei+1 )) 2jj(S 0 ?1 )Gi jjW 0 (jjEi jjW 0 + jjEi+1 jjW 0 ): 3
1
From (5.5)
Gil = ?
Z (i+1)4t i4t
;p
((C 0 He )l + 4t(C 0 Hei+ 2 )l ) ds 1
where the subscript corresponds the lth dual edge. By the quadrature rule, Gil vanishes for constant (C 0 He )l in time, so jGil j K2 jj(C 0 H_ e )jjL1 (i4tt(i+1)4t) where, by a scale change, K2 varies as 4t. By Bramble-Hilbert lemma jC 0 H_ e j K3 jH_ j(W 1 ( ))3 ;p
where, by scaling, K3
varies as h2? p . 3
jGil j K4 h2? 3 4 t
So
Z (i+1)4t
p
jGil j2 K4 h4? 6 (4t)3
jH_ j(W 1
;p
i4t
Using Cauchy's inequality
Z (i+1)4t
p
and
l
i4t
( l ))3 ds:
jH_ j2(W 1
( l ))3 ds
;p
Z (i+1)4t X F1 F1 X 6 4 ? 3 i 2 jH_ j2(W 1 ( ))3 ds jGl j K4 h (4t) i 4 t l=1 l=1 2 F F1 1 1? 2 Z (i+1)4t X X 6 p 4 ? 3 _ 1 K4 h (4t) jHj(W 1 ( ))3 ds i4t p
;p
l
p
p
p
l=1
;p
l=1
l
from Holder's inequality. Using (5.7) to estimate 4t we obtain (Gi ; Gi ) K4h3 4 tjH_ j2L2 ((i4t;(i+1)4t);(W 1 ( ))3 ) and from the de nition of jj jjW 0 p jj(S 0 ?1 )Gi jjW 0 K5 h 4tjjHjjH 1 ((i4t;(i+1)4t);(W 1 ( ))3 ) : ;p
;p
By Cauchy's inequality N X
Similarly, and
i=1
N X i=1
N p X
ai N (
i=1
a2i ) 12 and N 4 t = T
jj(S 0 ?1 )Gi jjW 0 K7 hjjHjjH 1 (0;T ;(W 1
;p
p
( ))3 ) :
jjS ?1 G i jjW K8 h 4tjjEjjH 1 ((i4t;(i+1)4t);(W 1 N X i=1
;p
jj(S ?1 )G i jjW K9 hjjEjjH 1 (0;T ;(W 1
;p
( ))3 )
( ))3 ) :
ANALYSIS OF A COVOLUME SCHEME FOR MAXWELL'S EQUATIONS
13
Collecting the terms from (i)-(iii), we obtain nally
max (jjEi jjW 0 +jjHi+ 2 jjW ) Kh(jjEjjH 1 (0;T ;(W 1 0iN ?1 1
;p
( ))3 ) +jjHjjH 1 (0;T ;(W 1;p ( ))3 ) )
and this proves (5.8).
6. On Rectangular Meshes The covolume scheme can be extended to rectangular meshes. In this case both primal and co-face are rectangles and orthogonal to each other. All the duality relations discussed in section 2 are preserved. This method is the standard Yee scheme [18]. Motivated by [11] we will show in this section that on nonuniform but rectangular grids (a.k.a graded grids) with maximum size h the covolume approximation of tangential components of electric eld E and normal components of magnetic eld H are second1 order in3space3 in2 jj jjW 0 and jj jjW respectively. Here we only need (E; H) in (L (0; T ; (H ( )) )) . This improves the norm used in [11], where (E; H) in (L1 (0; T ; (C 3 ( ))3 ))2 was assumed. Using Lemma 3.1 and the error functions introduced in section 4, we can rewrite (4.9) as 1d 2 2 _ 0 (6.1) 2 dt (jjH jjW + jjE jjW 0 ) = (E ; E )W ? (_H ; H )W : We need here two technical lemmas to estimate the terms on the right side of (6.1).
O’
P
4
H
κ
j
E
P
1
i
κ
hx hy Q
P
2
O P
3
z y x
Figure 2
Lemma 6.1. There exists u(t) 2 RF1 such that each component of u is a continuous linear functional of the magnetic eld H and (6.2)
(_H ; H )W = (_H ; u)W
14
with (6.3)
R.A. NICOLAIDES AND D-Q WANG
max(jjujjW ; jjju_ jjW ) Kh2jjHjj(H 3 ( ))3 :
Lemma 6.2. There exists v1 (t) 2 RM1 , v2(t) 2 RF1 such that each component of v1 and v2 is a continuous linear functional of the electric eld E and
(6.4) (_E ; E )W 0 = (v_ 1 ; E )W 0 + (_H ; v_ 2 )W with max(jjv_ 2 jjW ; jjv2 jjW ) Kh2 jjHjj(H 3 ( ))3 : (6.5) jjv_ 1 jjW 0 Kh2 jEj(H 3 ( ))3 ; Assuming these two lemmas we can prove Theorem 6.3. Suppose that (E; H) 2 L1(0; T ; (H 3( ))3 )2 satis es (2.1)-(2.5), and denote by (E; H ) the solution of (4.3)-(4.4) on nonuniform grids with maximum grid size h. Then (6.6) max (jj(E ? Ee )(t)jjW 0 + jj(H ? Hf )(t)jjW ) Kh2 jj(E; H)jj(L1 (0;T ;(H 3 ( ))3 ))2 : 0tT Proof of Theorem 6.3. : Substituting (6.2) and (6.4) into (6.1) and integrating from 0 to t1 , 1 (jj jj2 +jj jj2 )(t ) = ( ; u?v_ 2) (t )+Z t1 (v_ 1 ; ) 0 ( ) d +Z t1 ( ; u_ ?v2 ) ( ) d : H W EW H W 1 2 H W E W0 1 0 0 Applying Lemma 6.1, Lemma 6.2 and Cauchy's inequality proves (6.6).
The method of proving lemma 6.1 is to notice that if h1 6= h2 the quadrature rule
Z h2
?h1
f (x) dx (h1 + h2 )f (0)
is exact for constant functions but after adding the correction term 1 (h2 f (h ) ? h2 f 0(?h )) 1 1 2 2 2 is exact for linear polynomials. In the proof below we will manipulate these rst order correction terms to make them a discrete gradient. Proof of Lemma 6.1. For primal face i , (H )i , the di erence between the average of H ni over the face i and the average of the same quantity along the co-edge i0 = OO0 orthogonal to the face i (see Figure 2), vanishes for the constant eld H. So we write (H )i as (6.7) OO0 (H )i = OO0 ui + u~i where the rst order correction term u~i is u~i := 12 (OP1 2 H3z + h2xH1x + h2y H2y )(O) ? 21 (O0 P1 2 H3z + h2x H1x + h2y H2y )(O0 ) and ui from (6.7) vanishes for piecewise linear polynomial functions. We remark that u~i formed in the way indicated above since there are three co-edges emanating
ANALYSIS OF A COVOLUME SCHEME FOR MAXWELL'S EQUATIONS
15
from the center O. Note also that u~i is a discrete gradient form, i.e. u~ = B1 for some 2 RT . Using (3.8) in lemma 3.2, (_H ; H )W = (_H ; u)W + (_H ; D0 ?1 u~)W = (_H ; u)W + (_H ; P )W = (_H ; u)W + (D_H ; ) = (_H ; u)W : The last step follows from (4.8). Lemma 6.1 follows by a scale change argument for u as in the estimate (6.3). Proof of Lemma 6.2. The proof is similar to that of lemma 6.1. We note that the quadrature rule (see Figure 2)
Z
0
f (y; z ) d f (Q)s0j
j
is exact for constant functions f when the point Q is not the center of the rectangle O0 OBA. By a Taylor expansion we can show that after adding the following correction term 1 1 2 2 2 2 2 [fy (P2 )P2 Q ? fy (P1 )P1 Q ]P3 P4 + 2 [fz (P4 )P4 Q ? fz (P3 )P3 Q ]P1 P2 the quadrature rule is exact for linear polynomials f (y; z ). So on co-face 0j , we split (E )j as (6.8) s0j (E )j = s0j vj1 + vj3 where vj3 := (C 0 vj2 )0 i.e. vj3 is expressed as a discrete curl . Using Pi , i = 1; ; 4 Figure 2 shows how this discrete curl is formed in terms of vj2 : vj2 (P1 ) := 21 (h2x E2x(P1 ) + h2y E1y (P1 )) (6.9) and vj1 is computed from (6.8) and vanishes for piecewise linear polynomial elds E. By summation by parts (_E ; E )W 0 = (v_ 1 ; E )W 0 + (C 0 v_ 2 ; DE ) by (3:3) 0 2 1 0 by Lemma 3:1 = (v_ ; E )W + (D v_ ; CE ) by (4:6) = (v_ 1 ; E )W 0 + (D0 v_ 2 ; S _H) 2 1 by (3:1); = (v_ ; E )W 0 + (v_ ; _H )W and the estimate (6.5) follows from (6.9) and the fact that continuous linear functional vi1 vanishes for linear polynomials E. j
1. 2. 3. 4.
References Philippe G. Ciarlet, The Finite Element Method For Elliptic Problems, North-Holland Publishing Company, 1978 . G.Duvaut and J. Lions, Inequalities in Mechanics and Physics, Springer-Verlag, New York, 1976. R. Holland, Finite-di erence solution of Maxwell's equations in generalized nonorthogonal coordinates, IEEE Trans. Nuclear Science, Vol. NS-30, 1983, pp. 4589-4591. T.G. Jurgens, A. Taflove, K.Umashankar and T.G. Moore, Finite-di erence time-domain modeling of curved surfaces, IEEE Trans. Antennas and Propagation, Vol. 40, 1992, pp. 17031708.
16
R.A. NICOLAIDES AND D-Q WANG
5. T.G. Jurgens and A. Taflove, Three-dimensional contour FDTD modeling of scattering from single and multiple bodies, IEEE Trans. Antennas and Propagation, Vol. 41, 1993, pp. 1429-1438. 6. J.-F. Lee, Numerical solutions of TM scattering using obliquely Cartesian nite di erence time domain algorithm, IEE Proc. H (Microwaves, Antennas and Propagation), Vol. 140, 1992, pp. 23-28. 7. N. Madsen, Divergence preserving discrete surface integral methods for Maxwell's equations using nonorthogonal unstructured grids, J. of Computational Physics, 119, 1995, pp. 34-45. 8. N. Madsen and R.W. Ziolkowski, Numerical solution of maxwell's equations in the time domain using irregular nonorthogonal grids, Wave Motion 10, 1988, pp. 583-596. 9. N. Madsen and R.W. Ziolkowski, A three dimensional modi ed nite volume technique for Maxwell's equations, Electromagnetics, Vol. 10, 1990, pp. 147-161. 10. B.J. McCartin and J.F. Dicello, Three dimensional nite di erence frequency domain scattering computation using the Control Region Approximation, IEEE Trans. Magnetics, Vol. 25, No. 4, 1989, pp. 3092-3094. 11. Peter Monk and Endre Suli, A Convergence Analysis of Yee's Scheme on Nonuniform Grids, SIAM J. of Numerical Anal. Vol 31, 1994, pp. 393-412 . 12. Peter Monk and Endre Suli, Error estimates for Yee's method on nonuniform grids, IEEE Trans. Magnetics, Vol 30, 1994, pp. 3200-3203. 13. R.A. Nicolaides, Direct Discretization of Planar Div-Curl Problems, SIAM J. Numerical Anal. Vol 29, pp. 32-56 (1992) . 14. R.A. Nicolaides, Flow discretization by complementary volume techniques, AIAA paper 89-1978, Proceedings of 9th AIAA CFD Meeting, Bu alo, NY, June 1989. 15. R.A. Nicolaides, Analysis and convergence of the MAC scheme. 1. The linear Problem , SIAM J. Numerical Anal., Vol 29, No. 6, pp. 1579-1591, 1992. 16. R.A. Nicolaides and X.Wu, Analysis and convergence of the MAC scheme. 2. The nonlinear problem in the uniqueness case. To appear in Mathematics of Computation. 17. R.A. Nicolaides and X.Wu, Covolume Solutions of Three Dimensional Div-Curl Equations, submitted to SIAM J. of Numerical Anal. 18. K. Yee, Numerical Solution of Initial Boundary Value Problems Involving Maxwell's Equations in Isotropic Media, IEEE Trans. Antennas and propagation, AP-16 (1966) pp. 302-307 . Department of Mathematics, Carnegie Mellon University, Pittsburgh, PA 15213
Current address : Department of Mathematics, Carnegie Mellon University, Pittsburgh, PA 15213 E-mail address : [email protected] Department of Mathematics, Carnegie Mellon University, Pittsburgh, PA 15213
Current address : Department of Mathematical Sciences, University of Delaware, Newark, DE 19716 E-mail address : [email protected]
Maxwell's equations in three space dimensions. Several estimates are proved. First, an estimate for a semi-discrete scheme is given. Second, the estimate is extended to cover the classical interlaced time marching technique. Third, some of our unstructured mesh results are specialized to rectangular meshes, both uniform and nonuniform. By means of some additional analysis it is shown that the spatial convergence rate is one order higher than for the unstructured case.
1. Introduction Staggered mesh schemes for the numerical solution of Maxwell's equations go back as far as [18]. Over the years this scheme has seen a number of generalizations intended to enhance its usefulness. Mostly, these generalizations are aimed at increasing the geometric complexity that can be handled. Thus, in [3] a tensor formulation was given. This permits the method to be extended from rectangular meshes to meshes de ned by curvilinear coordinates. Two dimensional generalizations using quadrilateral meshes are given in [6] and [8]. For three dimensions, hexahedral mesh formulations were proposed in [7] and [9]. These formulations use interpolation to obtain values for nonbasic eld components{ those which, in the rectangular situation, lie along the primal or dual mesh edges. Another approach is the \control path" method of [4], [5]. The control path method uses the classical nite di erence approach of modifying the nite di erence stencil near the boundary of the domain. A natural generalization of the standard staggered mesh scheme to tetrahedral meshes was given in [10]. A similar technique was used independently in [14] to solve the Navier-Stokes equations. This \covolume" approach does not use any interpolations and is a very natural way to generalize the original rectangular staggered mesh approach. A characteristic feature of covolume schemes is the use of Voronoi-Delaunay mesh pairs to replace the rectangular staggered mesh arrangement. Very little rigorous analysis of staggered mesh schemes for electromagnetics is available. We know only of the rectangular mesh analysis of [11] and [12]. For the incompressible Navier{Stokes equations results may be found in [14], [15] and [16]. Our goal in this paper is to provide a rigorous proof of convergence for the covolume discretization applied to the interior problem for Maxwell's equations. Date : September 15, 1995 and, in revised form, February 20, 1996. 1991 Mathematics Subject Classi cation. Primary 65N30, 65N15, 35L50 . Key words and phrases. Maxwell's equations, covolume schemes, unstructured mesh, three dimensions, convergence . 1
R.A. NICOLAIDES AND D-Q WANG
2
The corresponding analysis for the exterior problem involves the application of radiation conditions at a nite distance from the scatterer and will be given in a subsequent report. The contents of the paper are as follows: in the rst section we state Maxwell's equations for bounded domains with perfect conductor boundary conditions. Following that is a section on Voronoi-Delaunay mesh notations and properties. The main result of section 3 is an error estimate for a semi{discrete covolume approximation to Maxwell's equations. In practice, a kind of leapfrog scheme is used to time march the discrete equations. Section 4 contains an error estimate for the resulting fully discrete approximation to Maxwell's equations. In the last section we specialize our results to the rectangular case. The rate of convergence is shown to be one order higher than for an arbitrary triangulation. 2. Maxwell's Equations Let be a bounded domain in R3 with boundary ? and unit outward normal n. Let the constants and denote, respectively, the electric and magnetic permeabilities of the medium occupying . Then if E(x; t) and H(x; t) denote, respectively, the electric and magnetic elds, Maxwell's equations [2] are: (2.1) Et ? curl H = J in
(0; T ); (2.2)
Ht + curl E = 0
in
(0; T );
(2.3) div( E) = and div( H) = 0 where J = J(x; t) is a known applied current and (x; t) is a charge density. We shall assume a perfect conductor boundary condition so that (2.4) En=0 on ? (0; T ): In addition, initial conditions are prescribed so that (2.5) E(x; 0) = E0(x) and H(x; 0) = H0(x) 8x 2
where E0 and H0 are given functions satisfying div( E0 ) = (x; 0); div( H0 ) = 0 and where it is assumed that @ = div J: (2.6) @t
Let Lp (0; T ; X ) denote the set of all strongly measurable functions u(t; ) from [0; T ] into the Hilbert space X such that
ZT 0
jju(t)jjpX dt < 1
1 p < 1:
are in Lp (0; T ; X ). Also We say u 2 W 1;p (0; T ; X ) if and only if both u and @u @t C m (0; T ; X ) denotes the space of m times continuously di erentiable functions from [0,T] into X . By H (curl; ) and H (div; ) we mean the Hilbert spaces de ned by H (curl; ) : = fv 2 (L2 ( ))3 ; curl v 2 (L2 ( ))3 g; H (div; ) : = fv 2 (L2 ( ))3 ; div v 2 L2 ( )g:
ANALYSIS OF A COVOLUME SCHEME FOR MAXWELL'S EQUATIONS
3
We shall assume the existence of (E; H) to (2.1)-(2.5) such that E; H 2 C 1 (0; T ; (L2( ))3 ) \ C 0(0; T ; H (curl; ) \ H (div; )): For this, it is sucient that J 2 C 0 (0; T ; (L2( ))3 ), 2 C 0 (0; T ; L2( )) and @ 2 L2(0; T ; H ?1( )), see [2]. @t
H H H
n
E
E H
H Figure 1 3. Mesh Notations and Discrete Vector Fields In this section some basic properties of dual mesh systems are introduced. These will lead to a detailed formulation and analysis of our numerical schemes in the following sections. Assume that the polyhedral domain has a primal family of nite element style tetrahedral partitions, parameterized by the maximum side length which is generically denoted by h. We will assume that the ratio of radii of circumscribing spheres and inscribed spheres of all the individual tetrahedra are uniformly bounded above and below as h approaches 0. A dual mesh is formed by connecting adjacent tetrahedral circumcenters and, in the case of tetrahedron with a face on a boundary, by connecting their circumcenters with those of their boundary faces. By elementary geometry these dual edges are perpendicular to the associated tetrahedral faces. These connections also form the edges of a set of polyhedra. It follows from elementary geometry that the edges of tetrahedra are perpendicular to and in one-to-one correspondence with the faces of dual polyhedra or \covolumes". The reciprocal orthogonality between edges and faces is the key to the results which follow. The N nodes of the tetrahedral mesh are assumed to be numbered sequentially in some convenient way, and likewise the T nodes of dual mesh. Similarly the F faces (edges) and M edges (faces) of the primal (dual) mesh are sequentially numbered. The individual tetrahedra, faces, edges and nodes of the primal mesh are denoted by i , j , k and l respectively. Those of the dual mesh are denoted by primed quantities such as j0 . A direction is assigned to each primal edge by the rule that the positive direction is from low to high node number. The dual edges are directed by the corresponding rule. We also denote F1 the number of tetrahedral interior
4
R.A. NICOLAIDES AND D-Q WANG
faces (or dual edges) and M1 the number of tetrahedral interior edges (or dual faces). Let sj denote the area of j and h0 j the length of j0 . In RF1 , where F1 denotes the number of interior primal faces, we will introduce the inner product (; )W de ned by X (u; v)W := (3.1) uj vj sj h0j = (Su; D0v) = (D0 u; Sv) 2 j
and denote the resulting inner product space by U and the associated norm by 1 jjujjW := (u; u)W2 : (3.2) In (3.1), (; ) denotes the standard Euclidean inner product, S := diag(sj ), D0 := diag(h0j ) and W := SD0 are F1 F1 invertible diagonal matrices. The norm de ned by (3.2) is clearly three times a discrete L2 norm. Similarly, we will introduce an inner product in RM , where M denotes the number of primal edges : X (3.3) uj vj s0j hj = (S 0 u; Dv) = (Du; S 0 v) (u; v)W 0 := 0 2
j
and denote the inner product space by U 0 . The associated norm is 1 (3.4) jjujjW 0 := (u; u)W2 0 : The notations in(3.3) and (3.4) correspond to those of (3.1) and (3.2). For example, s0 j denotes the area of dual face 0j . For each primal face i a discrete circulation is de ned by X ~ (3.5) uj hj : (Cu) := i
2@ j
i
Similarly, for each interior covolume face 0i the discrete circulation is X ~0 (3.6) (C 0 u)0 := uj hj : 0 2@0 A tilde on hj or h0j means that the quantity is to be taken with a negative sign if the dual edge is directed against the positive sense of description of @i or @0i i
j
i
respectively, and with a positive sign otherwise. The linear operators C and C 0 map from RM to RF1 and RF1 to RM1 , respectively. For each strictly interior dual edge j0 we can form a vector whose ith component is the sign of the orientation of the edge relative to the orientation of the ith strictly interior dual face. From these vectors we obtain the F1 M1 matrix G de ned as follows: 8 < 1 if j0 is oriented positively along 0i (G)ji := : ?1 if j0 is oriented negatively along 0i 0 if j0 does not meet 0i Let w 2 RM denote the vector whose kth component is the value assigned on the kth primal edge. Denote w1 2 RM1 the restriction of w to the interior primal edges. So wj? = 0 means that the components of w 2 RM related to the boundary edges are zero. A direct computation shows: Cw = GDw1 :
ANALYSIS OF A COVOLUME SCHEME FOR MAXWELL'S EQUATIONS
5
Also it can be veri ed that if v 2 RF1
C 0 v = GT D0 v:
From these two identities we can prove Lemma 3.1. With the above de nitions of w,v and w1 with wj? = 0, (Cw; D0 v) = (C 0 v; Dw1 ): Proof. : This is proved as follows: (C 0 v; Dw1 ) = (GT D0 v; Dw1 ) = (D0 v; GD0 w1 ) = (D0 v; Cw): This lemma provides a discrete analog of the integration formula
Z
curl E H dx =
Z
curl H E dx
which holds when E n = 0 on ?. For each tetrahedron i a discrete ux is de ned by X uj s~j ; 8u 2 RF1 : (Du)i := 2@ j
i
By s~j we mean sj negatively signed if the corresponding velocity component is directed towards the inside of i and positively signed otherwise. In addition to G we will introduce another matrix B1 of dimension F1 T 8 < 1 if j is oriented positively along @ i (B1 )ji := : ?1 if j is oriented negatively along @ i 0 if j does not meet @ i It can be checked directly that D = B1T S: Using B1 we will de ne the di erence operator P by P := D0 ?1 B1 ; 8 2 RT : Now we have
Lemma 3.2.
(3.7) B1T C = 0; and (3.8) (u; P )W = (Du; ); 8 2 RT ; u 2 RF1 : Proof. : For (3.7) see [17] Theorem 1. For (3.8), we have (u; P )W = (SD0 u; D0 ?1 B1 ) = (B1T Su; ) = (Du; ): We remark that (3.7)Zand (3.8) provide the Z discrete analogs of the identity div(curl u) = 0 and Green's formula u grad dx = (div u) dx for 2 H01 ( ); u 2 H (div; ),
respectively. For additional details and relations between other discrete operators in VoronoiDelaunay meshes, see [13] and [17].
R.A. NICOLAIDES AND D-Q WANG
6
4. Semi-Discrete Maxwell's Equations In this section the covolume method is used to obtain a spatial discretization of Maxwell's equations (2.1)-(2.5) from which we will obtain a formulation of a semi-discrete form of Maxwell's equations. A basic error estimate is given in Theorem 4.1. First we introduce, for general eld A, its \face averages" Af 2 RF and its \edge averages " Ae 2 RM as follows (refer to Figure 1): Z A := 1 A n ds; f
si Z A0f := s10 0 A n ds; i Z Ae := h1 A t d; j Z 1 0 Ae := h0 0 A t d; j i
i
j
j
where n and t denote the unit normal vector to the face i (or 0i ) and the unit tangent vector to the edge j (or j0 ), respectively. Error functions for primal edges and faces will be denoted by A := A ? Af ; A := A ? Ae ; A := Af ? Ae : Error functions for dual edges and faces are de ned similarly. As shown in Figure 1, for each tetrahedron i we will use the normal components of the magnetic eld H to its faces and the tangential components of the electric eld E in the directions of its edges. Now integrate both sides of (2.1) over the co-face 0j to obtain Z d ( E ) f j 0 0 0 s (4.1) ? (C H ) = J(x; t) dx: j
e
dt
0
j
j
Here, (Ef )j denotes the average of E nj over the face 0j where the nj is the unit vector in the direction of j , and (C 0 He )0 is the discrete circulation around the face 0j . Similarly, from (2.2), (4.2) s d(Hf )i + (CE ) = 0 j
i
dt
e
i
where (Hf )i denotes the average of H ni where ni is the unit normal to the face i and (CEe ) is the discrete circulation around the face i . Let E and H denote vectors of components in RM1 and RF1 respectively. Then (4.1) and (4.2) (which are exact) suggest the approximations ( where it is implied that components of E associated with boundary edges are zero, i.e. E j? = 0). (4.3) S 0 dE ? C 0 H = J;~ i
(4.4)
dt S dH dt + CE = 0
ANALYSIS OF A COVOLUME SCHEME FOR MAXWELL'S EQUATIONS
7
where J~ 2 RM1 with components given by the right hand of (4.1). Since S and S 0 are invertible, (4.3)-(4.4) is a system of linear ordinary di erential equations, and the existence and uniqueness of a solution follow from well known results. From (4.4) and (3.7) in lemma 3.2 we obtain
dtd (DH ) = dtd (B1T SH ) = ?B1T CE = 0: This shows the sense in which divH = 0 is satis ed at the discrete level in the covolume scheme. By subtracting (4.3) from (4.1) we have
(4.5)
S 0 ddt (E ? Ef0 ) ? C 0 (H ? He0 ) = 0:
Similarly, from (4.2) and (4.4) (4.6)
S ddt (H ? Hf ) + C (E ? Ee ) = 0
where (E; H) denotes the exact solution of (2.1)-(2.5) and where, by the boundary condition (2.4), (4.7) E j? = 0: Note that error in the magnetic eld H satis es the discrete solenoidal condition (4.8) D_H = 0: Multiplying (4.5) by DE and (4.6) by D0 H , and adding we obtain (4.9) (_0E ; E )W 0 + (_H ; H0 )W = (C 0 H0 ; DE ) ? (CE ; D0 H0 ) where the dots denote time di erentiations. The main result in this section is the following theorem: Theorem 4.1. Denote by (E; H ) the solution of (4.3) and (4.4) and by (E; H) 2 W 1;1 (0; T ; (W 1;p( ))3 )2 the solution of (2.1)-(2.5) with p > 2. Then we have the estimate max (jj(E ?Ee )(t)jjW 0 +jj(H ?He )(t)jjW ) Kh(jjE_ jjL1 (0;T ;(W 1 ( ))3 ) +jjH_ jjL1 (0;T ;(W 1 0tT
To prove this theorem we use (4.9) to get
;p
(_E ; E )W 0 + (_H0 ; H0 )W = (E ? 0E ) ; E + (H0 ? H ) ; H0 W W + (C 0 H0 ; DE ) ? (CE ; D0 H0 ): Since E nj? = 0, the components of E restricted to the boundary is 0. So
lemma 3.1 is applied to yield (C 0 H0 ; DE ) ? (CE ; D0 H0 ) = 0: By E ? 0E = Ef ? Ee ; H0 ? H = Hf ? He ;
;p
( ))3 ) ):
R.A. NICOLAIDES AND D-Q WANG
8
we obtain from (4.9) (4.10) 1 d (jj() 12 0 jj2 + jj() 12 jj2 ) = ((H ? H ) ; 0 ) + ((E ? E ) ; ) 0 E W0 f e f e EW H W HW 2 dt where () denotes time di erentiations. Integrating (4.10) from 0 to T and using Cauchy's inequality we obtain (jj() 12 H0 (T )jj2W + jj() 12 E (T )jj2W 0 )
2
ZT 0
(jj() 12 (He ? Hf ) jjW jj() 12 H0 (s)jjW + jj() 12 (Ee ? Ef ) jjW 0 jj() 12 E (s)jjW 0 ) ds:
Let t be such that 1 1 2 H (t)jjW + jj() 2 E (t)jjW 0 ): (jj() 12 H0 (t )jjW + jj() 12 E (t )jjW 0 ) = 0max ( jj ( ) tT Then (jj 12 H0 (t )jjW + jj 12 E (t )jjW 0 )2 2(jj 21 H0 (T )jj2W + jj 21 E (T )jj2W 0 )
4 4 So it follows that (4.11)
Z t Z0 T 0
(jj 12 (He ? Hf ) jjW jj 12 H0 (s)jjW + jj 12 (Ee ? Ef ) jjW 0 jj 21 E (s)jjW 0 ) ds
(jj 21 (He ? Hf ) jjW jj 12 H0 (s)jjW + jj 21 (Ee ? Ef ) jjW 0 jj() 12 E (s)jjW 0 ) ds:
1 1 0 2 (t)jjW + jj() 2 E (t)jjW 0 ) 4 max ( jj ( ) H 0tT
ZT 0
(jj() 12 (He ? Hf ) jjW + jj() 12 (Ee ? Ef ) jjW 0 ) ds:
Theorem 4.1 will follow from the next lemma: Lemma 4.2. Assume that E_ ; H_ 2 (W 1;p( ))3 , p > 2. Then there exists a generic constant K , such that jj(He ? Hf ) jjW KhjH_ j(W 1 ( ))3 : jj(Ee ? Ef ) jjW 0 KhjE_ j(W 1 ( ))3 : where si ); max( s0i ) K K 0 max max ( 0 2 i h2 i ;p
;p
hi
and K 0 is a constant.
i
Proof. For the primal face i ,
Z Z (He ? Hf ) = h10 0 H_ n d l ? s1 H_ n dx: i
i
i
i
By the Sobolev embedding theorem, V : H ! (H_ e ? H_ f ) de nes a bounded linear functional on (W 1;p (j ))3 . For constant H_ this linear functional vanishes in the union of two tetrahedra that share the same face i and consequently j(He ? Hf ) j K ( i [ l )jH_ jW 1 ( [ )3 where K ( i [ l ) is a constant. In order to estimate K ( i [ l ), we use a standard scale change argument as follows : ;p
i
l
ANALYSIS OF A COVOLUME SCHEME FOR MAXWELL'S EQUATIONS
9
Let the primal face i = i \ l be in the xy plane so that i0 is parallel to the z axis , changing the coordinates by
0 = A xy0 ; z = h1 z 0 where A is 2 2 matrix, it results that K ( i [ l ) depends on the quantity 1 0 1 h ? ? 1 ? 1 1 max(jjA jj2 (jdet Ajh1 ) ; 1 ) (det A) where 1p + p10 = 1. Let h0i denote the length of the co-edge and si the area of the primal face. A further calculation then shows that jdet Aj = c1 si , h1 = c2 h0i so that 1?1 p p (AT A) 2 4 T A) 1 det ( A s min ? ? 1 ? 1 i jjA jj2 (jdet Ajh1 ) = c3 c = c 3 4 1 1 1 (si h0i ) (si h0i ) h0 i where ci , i = 1; ; 4 are independent of h0i and si , and min is the least eigenvalue of the positive de nite matrix AT A. Collecting these results, K ( i [ l ) is bounded by s 12 ? 1 0 10 i ;hi : K ( i [ l ) max 1 1 i h 0 i si x y
p
p
p
p
p
p
p
p
p
p
p
p
Now we have
jj(He ? Hf ) jj2W =
F X i=1
K1 K2
si h0i j(He ? Hf ) j2 F X i=1 F X
1? 2 0 20 h i jH_ j2 2 ; 2 i i (W 1 ( [ ))3 h0 i si
si h0 max si
p
p
;p
p
p
i
l
h5? 6 jH_ j2(W 1 ( [ ))3 p
;p
i
l
i=1 where K2 depends only on max( s0i2 ). By Holder's inequality and since h3 F hi 0
K measure( ) we obtain
jj(He ? Hf ) jj2
W
K2
h5? p6
2 X X F F p _ jHj(W 1 ( [ ))3 1 p
;p
i
l
i=1 p?2 6 5 ? p = K3 h F p jH_ j2(W 1;p ( ))3 Kh2 jH_ j2(W 1;p ( ))3 : estimate for jj(Ee ? Ef ) jjW 0 is similar.
?2
p
p
i=1
The proof of the From the estimate above we know the constant K depends on the quantity 0
( si )): max(max ( s0i2 ); max i h2 i
hi
i
10
R.A. NICOLAIDES AND D-Q WANG
Theorem 4.1 now follows from (4.11) and lemma 4.2. 5. The Fully Discrete problem There are many possible time stepping methods that can be applied to (4.3) and (4.4). We will discuss a leapfrog scheme which is very popular in computational electromagnetics (see [18]). In this scheme we approximate E(t) at times tn = n 4 t, 0 n < 1 with a vector fE n g1 n=0 , and H(t) at time tn+ 21 with a vector 1 1 1 n + 2 2 fH gn=0 . The initial value H can be computed using, for example, a Taylor expansion and the given equations (2.1)-(2.2). Given (E n ; Hn+ 21 )n0 , the next approximation (E n+1 ; Hn+ 32 ) is obtained by solving the equations (5.1) S 0 (E n+1 ? E n ) ? 4tC 0 Hn+ 21 = J~n+ 21 (5.2) where
S (Hn+ 32 ? Hn+ 12 ) + 4tC E n+1 = 0 Z (n+1)4t J~n+ 21 = J~ dt: n4t
(5.1), (5.2) is an explicit scheme, so the existence and uniqueness of a solution are apparent. Using the error functions de ned in the last section, we can rewrite (5.1)-(5.2) as 1 (5.3) S 0 (nE+1 ? nE ) = (4t)C 0 Hn+ 2 + Gn ; 3 1 S (nH+ 2 ? Hn+ 2 ) = ?(4t)CEn+1 + Gn : (5.4) By a direct computation, Gn and G n are given by 1 (5.5) Gn = J~n+ 12 ? S 0 (Efn+1 ? Efn ) + 4tC 0 Hen+ 2 ; 1 3 (5.6) G n = S (Hfn+ 2 ? Hfn+ 2 ) ? 4tCEen+1: Our main estimate for the fully discrete scheme is then given in the following theorem: Theorem 5.1. Let (E n ; Hn+ 21 )Nn?01 denote the solution of (5.1)-(5.2), and let (E; H) 2 (H 1 (0; T ; (W 1;p( ))3 ))2 denote the solution of (2.1)-(2.4), p > 2. Under the stability condition (5.7) c 4 t < min(hij )
pM M
3
3 22
where c = ()? 21 is the speed of the light in the medium, M2 is the maximum of the ratios of the maximum to minimum side-lengths over the union of adjacent tetrahedra and M3 is the maximum number of edges over all co-faces, we have the following error estimate for the fully discrete scheme (5.1)-(5.2) (5.8) 1 max (jjE i ? Eei jjW 0 + jjHi+ 12 ? Hei+ 2 jjW ) Kh(jjEjjH 1 (0;T ;(W 1 ( ))3 ) + jjHjjH 1 (0;T ;(W 1 ( ))3 ) ): 0iN ?1 ;p
;p
ANALYSIS OF A COVOLUME SCHEME FOR MAXWELL'S EQUATIONS
11
Proof. Multiplying (5.3) by D(En + En+1 ) and (5.4) by D0 (Hn+ 2 + Hn+ 2 ) respectively and adding all these equations from n = 0; 1; ; N ? 1 gives 1
3
jjEN jj2W 0 +jjHN ? 2 jj2W = 4t(C 0 HN ? 2 ; DEN )+ 1
+
1
NX ?1
NX ?2
i=0
i=0
((Gi ; D(Ei + Ei+1 )) +
NX ?1 i=0
NX ?2
(Ei+1 ?Ei ; Ei )W 0 +
i=0
(Hi+ 2 ?Hi+ 2 ; Hi+ 2 )W 3
(G i ; D0 (Hi+ 2 + Hi+ 2 ))): 3
1
The proof has three steps : (i): First, 1 1 4t(C 0 HN ? 2 ; DEN ) = 4t(C 0 (SD0 )? 21 ((SD0 ) 12 HN ? 2 ); (DS 0 ?1 ) 12 ((DS 0 ) 21 EN )) 1 4tjj(SD0)? 21 C 0 (DS 0 ?1 ) 12 jj2 jjHN ? 2 jjW jjEN jjW 0 : From algebra, jj(SD0 )? 21 C 0 (DS 0 ?1 ) 12 jj2 is the largest singular value of the matrix and by Gershgorin's theorem 3 p 2 jj(SD0 )? 12 C 0 (DS 0 ?1 ) 21 jj2 max( maxij (hij )52 )(2 M3 ) minij (hij ) where maximum above is taken over all of the union of adjacent tetrahedra. Combining the results gives pM M 23 1 3 2 jj N ? 12 jj jj N jj 0 4t(C 0HN ? 2 ; DEN ) 4t 2min( hij ) H W E W
pM M 32 3 2 c(jj N jj2 + jj N ? 21 jj2 ) 4t min( E W0 H W hij ) 1 < (jjEN jj2W 0 + jjHN ? 2 jj2W ):
(ii) Ei+1 ? Ei can be estimated in the following way. We have, by integrating in time jEi+1 ? Ei j jj_E jjL1 (i4tt(i+1)4t): So M1 X s0k hk j(Ei+1 )k ? (Ei )k j2 k=1 Z (i+1)4t M1 X 0 sk hk ( j(_E )k j ds)2 i 4 t k=1 Z (i+1)4t X M1 4t s0k hk j(_E )k j2 ds i4t k=1 Z (i+1)4t = 4t jj_E jj2W 0 ds:
jjEi+1 ? Ei jj2W 0 =
By lemma 4.2
i4t
p
jjEi+1 ? Ei jjW 0 K1h 4tjjE_ jjL2 ((i4t;(i+1)4t);(W 1
;p
( ))3 ) :
1
1
R.A. NICOLAIDES AND D-Q WANG
12
Similarly,
p
jjHi+ 2 ? Hi+ 2 jjW K1 h 4tjjH_ jjL2 ((i4t;(i+1)4t);(W 1 ( ))3 ) : (iii) We have from the de nition of jj jjW 0 (Gi ; D(Ei + Ei+1 )) 2jj(S 0 ?1 )Gi jjW 0 (jjEi jjW 0 + jjEi+1 jjW 0 ): 3
1
From (5.5)
Gil = ?
Z (i+1)4t i4t
;p
((C 0 He )l + 4t(C 0 Hei+ 2 )l ) ds 1
where the subscript corresponds the lth dual edge. By the quadrature rule, Gil vanishes for constant (C 0 He )l in time, so jGil j K2 jj(C 0 H_ e )jjL1 (i4tt(i+1)4t) where, by a scale change, K2 varies as 4t. By Bramble-Hilbert lemma jC 0 H_ e j K3 jH_ j(W 1 ( ))3 ;p
where, by scaling, K3
varies as h2? p . 3
jGil j K4 h2? 3 4 t
So
Z (i+1)4t
p
jGil j2 K4 h4? 6 (4t)3
jH_ j(W 1
;p
i4t
Using Cauchy's inequality
Z (i+1)4t
p
and
l
i4t
( l ))3 ds:
jH_ j2(W 1
( l ))3 ds
;p
Z (i+1)4t X F1 F1 X 6 4 ? 3 i 2 jH_ j2(W 1 ( ))3 ds jGl j K4 h (4t) i 4 t l=1 l=1 2 F F1 1 1? 2 Z (i+1)4t X X 6 p 4 ? 3 _ 1 K4 h (4t) jHj(W 1 ( ))3 ds i4t p
;p
l
p
p
p
l=1
;p
l=1
l
from Holder's inequality. Using (5.7) to estimate 4t we obtain (Gi ; Gi ) K4h3 4 tjH_ j2L2 ((i4t;(i+1)4t);(W 1 ( ))3 ) and from the de nition of jj jjW 0 p jj(S 0 ?1 )Gi jjW 0 K5 h 4tjjHjjH 1 ((i4t;(i+1)4t);(W 1 ( ))3 ) : ;p
;p
By Cauchy's inequality N X
Similarly, and
i=1
N X i=1
N p X
ai N (
i=1
a2i ) 12 and N 4 t = T
jj(S 0 ?1 )Gi jjW 0 K7 hjjHjjH 1 (0;T ;(W 1
;p
p
( ))3 ) :
jjS ?1 G i jjW K8 h 4tjjEjjH 1 ((i4t;(i+1)4t);(W 1 N X i=1
;p
jj(S ?1 )G i jjW K9 hjjEjjH 1 (0;T ;(W 1
;p
( ))3 )
( ))3 ) :
ANALYSIS OF A COVOLUME SCHEME FOR MAXWELL'S EQUATIONS
13
Collecting the terms from (i)-(iii), we obtain nally
max (jjEi jjW 0 +jjHi+ 2 jjW ) Kh(jjEjjH 1 (0;T ;(W 1 0iN ?1 1
;p
( ))3 ) +jjHjjH 1 (0;T ;(W 1;p ( ))3 ) )
and this proves (5.8).
6. On Rectangular Meshes The covolume scheme can be extended to rectangular meshes. In this case both primal and co-face are rectangles and orthogonal to each other. All the duality relations discussed in section 2 are preserved. This method is the standard Yee scheme [18]. Motivated by [11] we will show in this section that on nonuniform but rectangular grids (a.k.a graded grids) with maximum size h the covolume approximation of tangential components of electric eld E and normal components of magnetic eld H are second1 order in3space3 in2 jj jjW 0 and jj jjW respectively. Here we only need (E; H) in (L (0; T ; (H ( )) )) . This improves the norm used in [11], where (E; H) in (L1 (0; T ; (C 3 ( ))3 ))2 was assumed. Using Lemma 3.1 and the error functions introduced in section 4, we can rewrite (4.9) as 1d 2 2 _ 0 (6.1) 2 dt (jjH jjW + jjE jjW 0 ) = (E ; E )W ? (_H ; H )W : We need here two technical lemmas to estimate the terms on the right side of (6.1).
O’
P
4
H
κ
j
E
P
1
i
κ
hx hy Q
P
2
O P
3
z y x
Figure 2
Lemma 6.1. There exists u(t) 2 RF1 such that each component of u is a continuous linear functional of the magnetic eld H and (6.2)
(_H ; H )W = (_H ; u)W
14
with (6.3)
R.A. NICOLAIDES AND D-Q WANG
max(jjujjW ; jjju_ jjW ) Kh2jjHjj(H 3 ( ))3 :
Lemma 6.2. There exists v1 (t) 2 RM1 , v2(t) 2 RF1 such that each component of v1 and v2 is a continuous linear functional of the electric eld E and
(6.4) (_E ; E )W 0 = (v_ 1 ; E )W 0 + (_H ; v_ 2 )W with max(jjv_ 2 jjW ; jjv2 jjW ) Kh2 jjHjj(H 3 ( ))3 : (6.5) jjv_ 1 jjW 0 Kh2 jEj(H 3 ( ))3 ; Assuming these two lemmas we can prove Theorem 6.3. Suppose that (E; H) 2 L1(0; T ; (H 3( ))3 )2 satis es (2.1)-(2.5), and denote by (E; H ) the solution of (4.3)-(4.4) on nonuniform grids with maximum grid size h. Then (6.6) max (jj(E ? Ee )(t)jjW 0 + jj(H ? Hf )(t)jjW ) Kh2 jj(E; H)jj(L1 (0;T ;(H 3 ( ))3 ))2 : 0tT Proof of Theorem 6.3. : Substituting (6.2) and (6.4) into (6.1) and integrating from 0 to t1 , 1 (jj jj2 +jj jj2 )(t ) = ( ; u?v_ 2) (t )+Z t1 (v_ 1 ; ) 0 ( ) d +Z t1 ( ; u_ ?v2 ) ( ) d : H W EW H W 1 2 H W E W0 1 0 0 Applying Lemma 6.1, Lemma 6.2 and Cauchy's inequality proves (6.6).
The method of proving lemma 6.1 is to notice that if h1 6= h2 the quadrature rule
Z h2
?h1
f (x) dx (h1 + h2 )f (0)
is exact for constant functions but after adding the correction term 1 (h2 f (h ) ? h2 f 0(?h )) 1 1 2 2 2 is exact for linear polynomials. In the proof below we will manipulate these rst order correction terms to make them a discrete gradient. Proof of Lemma 6.1. For primal face i , (H )i , the di erence between the average of H ni over the face i and the average of the same quantity along the co-edge i0 = OO0 orthogonal to the face i (see Figure 2), vanishes for the constant eld H. So we write (H )i as (6.7) OO0 (H )i = OO0 ui + u~i where the rst order correction term u~i is u~i := 12 (OP1 2 H3z + h2xH1x + h2y H2y )(O) ? 21 (O0 P1 2 H3z + h2x H1x + h2y H2y )(O0 ) and ui from (6.7) vanishes for piecewise linear polynomial functions. We remark that u~i formed in the way indicated above since there are three co-edges emanating
ANALYSIS OF A COVOLUME SCHEME FOR MAXWELL'S EQUATIONS
15
from the center O. Note also that u~i is a discrete gradient form, i.e. u~ = B1 for some 2 RT . Using (3.8) in lemma 3.2, (_H ; H )W = (_H ; u)W + (_H ; D0 ?1 u~)W = (_H ; u)W + (_H ; P )W = (_H ; u)W + (D_H ; ) = (_H ; u)W : The last step follows from (4.8). Lemma 6.1 follows by a scale change argument for u as in the estimate (6.3). Proof of Lemma 6.2. The proof is similar to that of lemma 6.1. We note that the quadrature rule (see Figure 2)
Z
0
f (y; z ) d f (Q)s0j
j
is exact for constant functions f when the point Q is not the center of the rectangle O0 OBA. By a Taylor expansion we can show that after adding the following correction term 1 1 2 2 2 2 2 [fy (P2 )P2 Q ? fy (P1 )P1 Q ]P3 P4 + 2 [fz (P4 )P4 Q ? fz (P3 )P3 Q ]P1 P2 the quadrature rule is exact for linear polynomials f (y; z ). So on co-face 0j , we split (E )j as (6.8) s0j (E )j = s0j vj1 + vj3 where vj3 := (C 0 vj2 )0 i.e. vj3 is expressed as a discrete curl . Using Pi , i = 1; ; 4 Figure 2 shows how this discrete curl is formed in terms of vj2 : vj2 (P1 ) := 21 (h2x E2x(P1 ) + h2y E1y (P1 )) (6.9) and vj1 is computed from (6.8) and vanishes for piecewise linear polynomial elds E. By summation by parts (_E ; E )W 0 = (v_ 1 ; E )W 0 + (C 0 v_ 2 ; DE ) by (3:3) 0 2 1 0 by Lemma 3:1 = (v_ ; E )W + (D v_ ; CE ) by (4:6) = (v_ 1 ; E )W 0 + (D0 v_ 2 ; S _H) 2 1 by (3:1); = (v_ ; E )W 0 + (v_ ; _H )W and the estimate (6.5) follows from (6.9) and the fact that continuous linear functional vi1 vanishes for linear polynomials E. j
1. 2. 3. 4.
References Philippe G. Ciarlet, The Finite Element Method For Elliptic Problems, North-Holland Publishing Company, 1978 . G.Duvaut and J. Lions, Inequalities in Mechanics and Physics, Springer-Verlag, New York, 1976. R. Holland, Finite-di erence solution of Maxwell's equations in generalized nonorthogonal coordinates, IEEE Trans. Nuclear Science, Vol. NS-30, 1983, pp. 4589-4591. T.G. Jurgens, A. Taflove, K.Umashankar and T.G. Moore, Finite-di erence time-domain modeling of curved surfaces, IEEE Trans. Antennas and Propagation, Vol. 40, 1992, pp. 17031708.
16
R.A. NICOLAIDES AND D-Q WANG
5. T.G. Jurgens and A. Taflove, Three-dimensional contour FDTD modeling of scattering from single and multiple bodies, IEEE Trans. Antennas and Propagation, Vol. 41, 1993, pp. 1429-1438. 6. J.-F. Lee, Numerical solutions of TM scattering using obliquely Cartesian nite di erence time domain algorithm, IEE Proc. H (Microwaves, Antennas and Propagation), Vol. 140, 1992, pp. 23-28. 7. N. Madsen, Divergence preserving discrete surface integral methods for Maxwell's equations using nonorthogonal unstructured grids, J. of Computational Physics, 119, 1995, pp. 34-45. 8. N. Madsen and R.W. Ziolkowski, Numerical solution of maxwell's equations in the time domain using irregular nonorthogonal grids, Wave Motion 10, 1988, pp. 583-596. 9. N. Madsen and R.W. Ziolkowski, A three dimensional modi ed nite volume technique for Maxwell's equations, Electromagnetics, Vol. 10, 1990, pp. 147-161. 10. B.J. McCartin and J.F. Dicello, Three dimensional nite di erence frequency domain scattering computation using the Control Region Approximation, IEEE Trans. Magnetics, Vol. 25, No. 4, 1989, pp. 3092-3094. 11. Peter Monk and Endre Suli, A Convergence Analysis of Yee's Scheme on Nonuniform Grids, SIAM J. of Numerical Anal. Vol 31, 1994, pp. 393-412 . 12. Peter Monk and Endre Suli, Error estimates for Yee's method on nonuniform grids, IEEE Trans. Magnetics, Vol 30, 1994, pp. 3200-3203. 13. R.A. Nicolaides, Direct Discretization of Planar Div-Curl Problems, SIAM J. Numerical Anal. Vol 29, pp. 32-56 (1992) . 14. R.A. Nicolaides, Flow discretization by complementary volume techniques, AIAA paper 89-1978, Proceedings of 9th AIAA CFD Meeting, Bu alo, NY, June 1989. 15. R.A. Nicolaides, Analysis and convergence of the MAC scheme. 1. The linear Problem , SIAM J. Numerical Anal., Vol 29, No. 6, pp. 1579-1591, 1992. 16. R.A. Nicolaides and X.Wu, Analysis and convergence of the MAC scheme. 2. The nonlinear problem in the uniqueness case. To appear in Mathematics of Computation. 17. R.A. Nicolaides and X.Wu, Covolume Solutions of Three Dimensional Div-Curl Equations, submitted to SIAM J. of Numerical Anal. 18. K. Yee, Numerical Solution of Initial Boundary Value Problems Involving Maxwell's Equations in Isotropic Media, IEEE Trans. Antennas and propagation, AP-16 (1966) pp. 302-307 . Department of Mathematics, Carnegie Mellon University, Pittsburgh, PA 15213
Current address : Department of Mathematics, Carnegie Mellon University, Pittsburgh, PA 15213 E-mail address : [email protected] Department of Mathematics, Carnegie Mellon University, Pittsburgh, PA 15213
Current address : Department of Mathematical Sciences, University of Delaware, Newark, DE 19716 E-mail address : [email protected]
