Time domain CFIEs for electromagnetic scattering problems

Applied Numerical Mathematics 79 (2014) 62–78

Contents lists available at SciVerse ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Time domain CFIEs for electromagnetic scattering problems Qiang Chen a,1 , Peter Monk b,∗,1 a b

Department of Computing and Mathematical Sciences, California Institute of Technology, Pasadena, CA 91125, USA Department of Mathematical Sciences, University of Delaware, Newark, DE 19716, USA

a r t i c l e

i n f o

Article history: Available online 26 March 2013 Keywords: Time domain Combined field integral equation Convolution quadrature Error analysis

a b s t r a c t Time domain integral equations complement other methods for solving Maxwell’s equations by handling infinite domains without difficulty and by reducing the computational domain to the surface of the scatterer. In this paper we study the discretization error when convolution quadrature is used to discretize two new regularized combined field integral equation formulations of the problem of computing scattering from a bounded perfectly conducting obstacle. © 2013 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction We shall derive error estimates for approximating a time domain electromagnetic scattering problem via two new regularized time domain combined field integral equations (CFIEs). Time domain integral equations are gaining popularity for solving electromagnetic scattering problems because they handle infinite domains easily and fast solution techniques are available [31,1,19,30]. In addition they reduce the computational domain to the boundary of the scatterer (the advantage of this dimension reduction by solving only on the boundary of the domain is less evident in the time domain because of the need to store a history of the solution there). The standard method for discretizing these equations is a space–time Galerkin method (e.g. [31] for the electric field integral equation case), but then retarded potential integrals need careful evaluation which complicates the use of high order and curvilinear elements. As an alternative we shall analyze the use of convolution quadrature (CQ) to discretize in time [27]. This has already been used to discretize the related time domain electric field integral equation and provides a stable numerical time marching scheme [4,27,34] that can be used for curvilinear and higher order elements [34]. We wish to approximate the electromagnetic field scattered by a bounded metallic object which we model by a perfectly electrically conducting (PEC) boundary condition. The scatter is assumed to occupy an open bounded Lipschitz polyhedral domain having connected complement and denoted by Ω ⊂ R3 . Let Γ := ∂Ω be the boundary of Ω and n be the unit outward normal on the boundary Γ . The unbounded exterior of Ω is denoted by Ω + := R3 \Ω . The domain Ω + outside the scatterer is assumed to be free space and so is homogeneous and isotropic. Then, denoting by T > 0 the final time of the calculation, the time domain electric and magnetic scattered fields (E , H ) satisfy the Maxwell system

∂E − curl H = 0, in Ω + × (0, T ], ∂t ∂H μ0 + curl E = 0, in Ω + × (0, T ], ∂t

ε0

* 1

Corresponding author. E-mail addresses: [email protected] (Q. Chen), [email protected] (P. Monk). Research supported in part by NSF grant DMS-1114889.

0168-9274/$36.00 © 2013 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.apnum.2013.03.005

(1.1a) (1.1b)

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

E ×n=G,

63

on Γ × (0, T ],

E (·, t ) = 0,

in Ω

H (·, t ) = 0,

+

in Ω

(1.1c)

for t  0,

+

(1.1d)

for t  0,

(1.1e)

where the boundary data G is given in terms of the incident field E as follows. The incident field is assumed to be a smooth solution of (1.1a)–(1.1b) in a neighborhood of Ω that vanishes in a neighborhood of Γ for t  0. Then G = −E inc × n on Γ because of the assumed PEC boundary condition on the scatterer. The constants ε0 and μ0 are the permittivity and permeability for free space respectively. After a change of variables (or by choosing appropriate units) [18], we can assume ε0 = 1 and μ0 = 1. By taking the time derivative of the electric equation and using the magnetic equation, we can rewrite Eqs. (1.1) in terms of electric scattered field alone as follows: inc

∂ 2E + curl curl E = 0, in Ω + × (0, T ], ∂t2 E × n = G , on Γ × (0, T ], ∂E E (·, 0) = (·, 0) = 0, in Ω + . ∂t

(1.2a) (1.2b) (1.2c)

The focus of this paper is on solving the time domain Maxwell equations (1.2) by time domain boundary integral equation (TD-BIE) methods. To derive a standard classical TD-BIE, the following time domain single layer ansatz is often assumed for the solution E of problem (1.2):

t  E (x, t ) =

 τ

t  k(x − y, t − τ )ηt (y, τ ) dτ dσ y − ∇

0 Γ

k(x − y, t − τ ) divΓ 0 Γ



η(y, γ ) dγ dτ dσ y ,

(1.3)

0

for x ∈ / Γ , and for some unknown tangential boundary field η and where divΓ denotes the surface divergence. Here the kernel k is the fundamental solution of the wave equation given by

k(x, t ) :=

δ(t − |x|) , 4 π |x |

x = 0,

with δ denoting the delta distribution. Proceeding formally, the field η can then be related to the boundary data G by taking the tangential component of the trace of the above representation (1.3) on Γ . Let Π T denote the tangential projection on Γ , given by

ΠT v = vT := n × (v × n) = v − (v · n)n, and let ∇Γ denote the surface gradient. Then using the continuity of the tangential trace of the single layer for Maxwell’s equations (under suitable conditions on η and Γ ) we obtain the following time domain electric field integral equation (TD-EFIE)

t  V (∂t )η(x, t ) := Π T

 τ

t  k(x − y, t − τ )ηt (y, τ ) dτ dσ y − ∇Γ

0 Γ

= n(x) × G (x),

k(x − y, t − τ ) divΓ 0 Γ



η(y, γ ) dγ dτ dσ y 0

x ∈ Γ.

(1.4)

For a survey of time dependent integral equations and their application to boundary value problems including Maxwell’s equations see [20]. The TD-EFIE is a popular choice when solving Maxwell’s equations, but to avoid possible numerical instability [21] it is often preferred to use a combined field integral equation (CFIE). In particular, the following classical ansatz represents E in terms of a combination of electric and magnetic operators. We suppose there is a tangential field ξ(x, t ) such that the solution E of problem (1.2) has the form, for x ∈ /Γ,

t  E (x, t ) =

k(x − y, t − τ )ξt (y, τ ) dσ y dτ − ∇ 0 Γ

t +



 τ k(x − y, t − τ )

curl 0

 τ

t 

Γ

k(x − y, t − τ ) divΓ 0 Γ



ξ(y, γ ) dγ dσ y dτ , 0

 ξ(y, γ ) dγ dσ y dτ

0

(1.5)

64

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

for some unknown tangential surface vector field ξ . We prefer not to accumulate the charge term as the solution progresses so define

t ψ(x, t ) =

ξ(x, γ ) dγ . 0

Then taking the tangential component trace from the exterior to the boundary Γ , and using the jump conditions for the magnetic vector potential (see Theorem 2.3 later in this paper), we obtain a preliminary version of the combined field integral equation: for x ∈ Γ , t > 0, we seek ψ such that

t  C (∂t )ψ(x, t ) := Π T

t  k(x − y, t − τ )ψtt (y, τ ) dσ y dτ − ∇Γ

0 Γ

+

ψ(x, t ) 2

k(x − y, t − τ ) divΓ ψ(y, τ ) dσ y dτ 0 Γ

t + ΠT

 k(x − y, t − τ )ψ(y, τ ) dσ y dτ

curl 0

Γ

= n(x) × G (x).

(1.6)

As we shall see, this version of the TD-CFIE is not easy to analyze and, motivated by [13], we shall introduce two regularizations that render the TD-CFIE more analytically tractable without destroying the desirable property of avoiding certain interior resonances that appear for the frequency domain EFIE and may result in numerical instability in the time domain. The analysis of CQ applied to a convolution is based on understanding the corresponding convolution in the Fourier– Laplace domain (i.e. after taking the Fourier–Laplace transform in time) [3,27] and we shall describe this technique and recall some results in the next section. After this preliminary discussion, we shall introduce the two regularized TD-CFIEs that are the subject of this paper and establish their mapping properties as required by Lubich’s theory [27]. Then, for one of these regularized systems, we shall proceed to give a complete error analysis. 2. Preliminary results and notation 2.1. Traces and spaces To discuss boundary integral equation methods, some trace operators and corresponding trace spaces are needed, and these are summarized next. For a smooth vector function u defined in Ω , the tangential trace γt− : u → u|Γ × n on Γ is given by





γt− u(x) := lim+ u x − hn(x) × n(x), h→0

for all x ∈ Γ where the normal is well defined. The superscript “−” denotes taking the trace from the interior of Ω onto the boundary Γ . Similarly, if u is a smooth vector field in Ω + , the exterior tangential trace γt+ can also be defined. The tangential component traces γ T± : u → n × (u|Γ × n) and the Neumann traces γ N± : u → curl u|Γ × n on Γ can also be defined with the same meaning of “±” signs as before. For solutions of Maxwell’s equations an appropriate solution space is





H(curl, D ) := u ∈ L2 ( D ), curl u ∈ L2 ( D )

where we shall choose D = Ω or Ω + . It is then common to use a special weighted norm as follows. Let s ∈ C, s > σ0 > 0 for some constant σ0 (s will actually be the Fourier–Laplace transform parameter in later theory). We define the s-dependent norm

u 2Hs (curl, D ) := |s|2 u 2L2 ( D ) + curl u 2L2 ( D ) for u ∈ H(curl, D ). This norm is related to the energy norm in the time domain by noting that multiplication by s in the Fourier–Laplace transform domain corresponds to the time derivative in the time domain. When the boundary Γ is smooth enough, trace spaces for H(curl, D ) have been defined and used for long time [14]. For more general boundaries Γ , definitions of trace spaces are not as clear as in the smooth case. For a Lipschitz (curvilinear) polyhedron, Buffa et al. [6–11] have given a characterization of trace spaces that we adopt here, and we now give a few brief comments about these spaces. 1/ 2 Define the Hilbert space H (Γ ) := γt H1 ( D ), where D = Ω or Ω + . In the rest of this section we shall omit the superscript “±”, which should be clear in the context of D. Then

γt : H1 ( D ) → H1 /2 (Γ ) is continuous and surjective [12]. Let

{Γ1 , . . . , Γn } be the set of open planar faces of Γ , denote ei , j = Γ i ∩ Γ j , i , j = 1, . . . , n, then ei , j is either empty if Γi and Γ j are not adjacent, or a straight line otherwise. For u ∈ C ∞ ( D ), define

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

divΓ

γt u|Γ j ,

divΓ

γt u :=

65

on Γ j ,

((γt u|Γi ) · ni , j + (γt u|Γ j ) · n j ,i )δi , j , on e i , j ,

where δi , j is the delta distribution with support on ei , j , ni , j is the unit outward vector in the plane of Γ j and perpendicular −1/2

to ei , j , and divΓ denotes the usual surface divergence [18,29]. This operator can be extended to H

1/ 2 space of H (Γ ) with the pivot space Lt2 (Γ ) := {ϕˆ ∈ L2 (Γ ), ϕˆ · n = 0 a.e. on Γ }. We now define −1/2

H

(Γ ) which is the dual

  −1/2 (divΓ , Γ ) := ϕˆ ∈ H (Γ ), divΓ ϕˆ ∈ H −1/2 (Γ ) .

Following [7], another kind of trace spaces can be defined 1/2

H⊥ (Γ ) := γ T H1 ( D ), and −1/2

H⊥

  −1/2 (curlΓ , Γ ) := ϕˆ ∈ H⊥ (Γ ), curlΓ ϕˆ ∈ H −1/2 (Γ ) ,

−1/2

1/ 2

where H⊥ (Γ ) is the dual space of H⊥ (Γ ) by pivot space Lt2 (Γ ), and curlΓ denotes the surface curl [18,29]. We now recall the following trace theorem. Theorem 2.1. (See [8, Proposition 4.5, Theorem 5.4].) Given a Lipschitz polyhedral domain Ω , the tangential trace and the tangential component trace operators 1/2 γt : H(curl, Ω) → H− (divΓ , Γ ),

1/2 γT : H(curl, Ω) → H− (curlΓ , Γ ), ⊥

are linear, continuous and surjective. Remark 2.2. Using the above theorem, it is easy to verify the following inequalities, for u ∈ H(curl, D )

γt u H−1/2 (div

Γ ,Γ )



γT u H−1/2 (curl ⊥

−1/2

and given g ∈ H

 c (Γ ) u H(curl, D )  c (Γ, σ0 ) u Hs (curl, D ) ,

Γ ,Γ )

 c (Γ ) u H(curl, D )  c (Γ, σ0 ) u Hs (curl, D ) ,

(divΓ , Γ ), there exists u0 ∈ H(curl, D ), such that γt u0 = g, and

u0 Hs (curl, D )  c (σ0 )|s| u0 H(curl, D )  c (Γ, σ0 )|s| g H−1/2 (div

−1/2

Similar results for g ∈ H⊥

Γ ,Γ )

.

(curlΓ , Γ ) hold.

2.2. Fourier–Laplace domain EFIE and classical CFIE As mentioned in the introduction, the main tool for analyzing CQ methods for time domain integral equations is the analysis of corresponding operators in the Fourier–Laplace transform domain. Of course, for a suitable function f , a formal definition of the Fourier–Laplace transform is

∞ L f (s) :=

e −st f (t ) dt ,

0−

where the Fourier–Laplace parameter s is a complex number and for our application s > σ0 for some constant σ0 > 0. Formally applying the Fourier–Laplace transform to (1.2) we obtain the problem of finding, for each s, E(x) = (L E ) ∈ H(curl, Ω + ) that satisfies

s2 E + curl curl E = 0, E × n = g,

on Γ,

in Ω + ,

(2.1a) (2.1b)

where g(x) := (L G )(x, s) and we usually suppress the s dependence of functions. We can then write down Fourier–Laplace domain boundary integral representations for the solution of the above problem (or equivalently take the Fourier–Laplace transform of TD-BIEs directly). Recall that the fundamental solution of the Laplace domain Maxwell system (2.1a) is

66

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

Φs (x) :=

e −s|x| 4 π |x |

x = 0,

,

which can be obtained through the Laplace transform of k(x, t ) directly. This is analogous to the fundamental solution in frequency domain [18,29]. The following result, which we have already used formally in the introduction, gives the jump properties of the vector single layer potential. Theorem 2.3. (See [28,12].) Let Ω be a Lipschitz domain with boundary Γ and let a ∈ H −1/2 (divΓ , Γ ) denote a tangential vector field. Suppose



A(x) :=

Φs (x, y)a(y) dσ y ,

x∈ / Γ.

Γ 1 Then the potential A ∈ H loc (R3 ). In particular, ∀x ∈ Γ ,

γt+ A(x) = γt− A(x), and

γn+ A(x) = γn− A(x), where γn± A(x) := limh→0+ A(x ± hn(x)) · n(x). Moreover, for x ∈ Γ ,



±

γN A(x) =





curlx Φs (x, y)a(y) × n(x) dσ y ∓

1 2

a(x).

Γ

Remark 2.4. 1. For the simpler case of a smooth domain see [18, Theorems 6.11, 6.12]. 2. Similar continuity and jump conditions hold for time domain potential, for x ∈ / Γ (see [3,24])

t  B(x, t ) :=

k(x − y, t − τ )σ (y, τ ) dτ dσ y . 0 Γ

We now introduce the standard notation



[u ]|Γ = u − Γ − u + Γ which is used to denote the jump of the function u across the boundary Γ . In order to analyze the EFIE (1.4), Terrasse [31] first takes the Fourier–Laplace transform of the representation formula (1.3) to obtain the following representation ansatz for solutions of (2.1a): for x ∈ /Γ,



1



Φs (x, y)ϕˆ (y) dσ y − ∇

E(x) = s

Φs (x, y) divΓ ϕˆ (y) dσ y .

s

Γ

(2.2)

Γ

Taking tangential traces on the boundary using Theorem 2.3 gives the Laplace domain EFIE, for x ∈ Γ ,



V (s)ϕˆ (x) := Π T s

1



Φs (x, y)ϕˆ (y) dσ y − ∇Γ

Φs (x, y) divΓ ϕˆ (y) dσ

s

Γ

Γ

= n(x) × g(x).

(2.3)

This equation can also be obtained by taking the Fourier–Laplace transform of (1.4). For the solution of the partial differential equation (2.1), extending the result of Terrasse [31] to Lipschitz polyhedral domains, we have the following result: −1/2

Theorem 2.5. (See [17,15].) Given a Lipschitz polyhedral domain Ω with the boundary Γ and g ∈ H

a unique solution E ∈ H(curl, Ω + ), and

E Hs (curl,Ω + )  c (Γ, σ0 )|s|2 n × g H−1/2 (curl ⊥

Γ ,Γ )

.

(divΓ , Γ ), problem (2.1) has

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

67

For the EFIE boundary integral problem (2.3), it follows [15,17,32], that the following bounds hold for Lipschitz polyhedral domains (see [31] for the case of a smooth domain): −1/2

Theorem 2.6. The Laplace domain EFIE operator V (s) : H

   V (s)  c (Γ, σ0 )|s|.

(divΓ , Γ ) → H⊥ (curlΓ , Γ ) is bounded, and

Moreover, the operator V (s) is invertible, and

 −1   V (s)  c (Γ, σ0 )|s|2 .

The equivalence of problems (1.2) and (2.1) (based on the Fourier–Laplace transform and inverse Fourier–Laplace transform) implies the existence and uniqueness of the time domain electric field solution E of (1.2), and the invertibility of TD-EFIE operator V (∂t ) in suitable Bochner spaces (see [24,26,27] for the acoustic case). We can now consider a CFIE in the Fourier–Laplace domain. For convenience define



K (s)ϕˆ (x) := Π T





curlx Φs (x, y)ϕˆ (y) dσ y , Γ

I n ϕˆ (x) := ϕˆ (x) × n(x), for any sufficiently smooth tangential vector field ϕˆ and x ∈ Γ . Taking the Fourier–Laplace transform of the representation (1.5) we obtain an ansatz for the solution E of problem (2.1) given by



E(x) = s





Φs (x, y)ϕˆ (y) dσ y − ∇

2

Γ

Φs (x, y)ϕˆ (y) dσ y ,

Φs (x, y) divΓ ϕˆ (y) dσ y + curl Γ

Γ

for x ∈ / Γ and for some unknown tangential vector field ϕˆ . Taking the tangential component trace from the exterior to the boundary Γ , and using the jump condition in Theorem 2.3, we obtain the Laplace domain combined field integral equation (CFIE), for x ∈ Γ ,

1 ΠT E(x) = C s ϕˆ (x) := sV (s)ϕˆ (x) + K (s) + I n ϕˆ (x) = n(x) × g(x),

(2.4)

2

where V (s) is the Fourier–Laplace domain EFIE operator defined in (2.3). This is the Fourier–Laplace version of (1.6). We now consider the existence and uniqueness of solutions of the Fourier–Laplace domain CFIE problem (2.4). If the boundary Γ is smooth, then K (s) : H−1/2 (divΓ , Γ ) → H−1/2 (curlΓ , Γ ) is compact [29, Section 5.5] (the space on the boundary here is defined in the usual way [14]). Then following [12] Gårding’s inequality can be proved for sV (s) + (1/2) I n , which leads to an existence result (because uniqueness for the scattering problem in the Fourier–Laplace domain is easy to verify). But for a Lipschitz boundary Γ , the compactness of K (s) is not available to the best of our knowledge. Another combined field ansatz for E is suggested in the book of Colton and Kress [18]. For x ∈ / Γ , E is represented by







Φs (x, y) V 02 ϕˆ (y) dσ y − s∇

E(x) = s3 Γ

Φs (x, y)ϕˆ (y) dσ y ,

Φs (x, y) divΓ V 02 ϕˆ (y) dσ y + curl Γ

Γ

where V 0 := V (s)|s=0 . Taking the tangential component trace on Γ , we have, for x ∈ Γ the following regularized CFIE



s2 V (s) V 02 ϕˆ (x) +

K (s) +

1 2

In

ϕˆ (x) = n(x) × g(x).

(2.5)

For smooth domains and when s is purely imaginary (the frequency domain case) Colton and Kress use the space 0,α

Td

  (Γ ) := v ∈ T 0,α (Γ ): divΓ v ∈ C 0,α (Γ )

with C 0,α (Γ ) denoting the usual Hölder space and T 0,α (Γ ) denoting the Hölder continuous tangential space. Then they 0,α 0,α prove that s2 V (s) V 02 + K (s) : T d (Γ ) → T d (Γ ) is compact [18, Theorem 6.19]. Then the existence of solutions to (2.5) is proved once uniqueness has been verified. Unfortunately, this result again is only valid for a smooth boundary Γ . We shall address this problem in Section 3 and subsequent sections. These considerations motivated by Buffa and Hiptmair [12] to introduce a new regularized CFIE in the frequency domain suitable for Lipschitz domains. The first of our TD-CFIEs is the extension to the time domain of their method, and we shall describe it in detail in Section 3. 2.3. Convolution quadrature We shall now summarize some aspects of convolution quadrature (CQ). For full details, see the original paper of Lubich [27] and for the time domain EFIE for Maxwell’s equations see [15,17]. The CQ method applies to any convolution

68

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

(although the properties of the kernel need to be understood before an appropriate choice of CQ can be made). As we noted in [15,34,17] this is a standard technique in control theory. In presenting CQ we follow [25,23]. For simplicity let us just consider a pure time convolution and ignore the spatial dependence. In particular let

t K (∂ t ) g :=

k(t − τ ) g (τ ) dτ 0

where k is a possibly singular kernel (in our case, it is a translated delta function) and g is a smooth causal function (g (t ) = 0 if t  0). The first step in discretizing this equation is to formally write it using the Fourier–Laplace transform of

kˆ (s) = L (k)(s), and the inverse Fourier–Laplace transform as

t K (∂ t ) g =



1 2π i

0

=

2π i

kˆ (s)

σ +i R

The function G (t , s) :=

G  = sG + g ,

t 0



σ +i R



1



kˆ (s) exp is(t − τ ) ds g (τ ) dτ

t





exp is(t − τ ) g (τ ) dτ ds. 0

exp(is(t − τ )) g (τ ) dτ satisfies the ordinary differential equation

G (0, s) = 0.

t > 0,

So far we have not discretized the problem, but now G can be discretized using a k-step multistep method with time step t > 0 of the form k 

α j G n+ j−k (s) = t

k 

j =0

  β j sG n+ j −k (s) + gn+ j −k

(2.6)

j =0

where gn+ j −k = g ((n + j − k)t ), and G 0 (s) = 0. Then the discrete convolution is given by



  1 K ∂tt g (nt ) =

2π i

kˆ (s)G n (s) ds.

σ +i R

This would not be very helpful in its current form, however by using generating functions (in effect the z-transform, or discrete Laplace transform) it can be shown [27] that the above equality can be rewritten as a discrete convolution





K ∂tt g (nt ) =

n  j =1

w n−t j g j

where the convolution coefficients w n−t j (which depend on kˆ (s)) are given as follows. Let generating polynomials for the multistep method (2.6)

k

j =0

γ (ζ ) = k

α j ζ k− j

j =0 β j ζ

k− j

.

Then the weights are given by requiring that

kˆ



γ (ζ ) t

=

∞ 

t j w j ζ

j =0

for all ζ in a neighborhood of the 0 ∈ C. Clearly an important case for us is when

k(t ) = δ(t − R ) for different values of R = x  0. In this case

kˆ (s) = exp(−sR ).

γ (ζ ) denote the quotient of

(2.7)

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

69

For the BDF2 second order time stepping scheme, note that [25]

α0 = 3/2,

α1 = −2,

α2 = 1/2,

β0 = 1

and all remaining coefficients are zero in (2.7) and k = 2 so

γ (ζ ) =

1 2

 ζ 2 − 4ζ + 3 .

They then show that

w nt =

1



n!

R

n/2

exp −

2t

3R



2t

Hn

2R

t

where H n is the nth Hermite polynomial. Coefficients for other time stepping schemes can also be computed (for example, for the backward Euler scheme). CQ is not limited to multistep methods. Implicit Runge–Kutta methods are also effective for Maxwell’s equations [33] and have been analyzed for the EFIE by Veit [32]. As might be expected, the Fourier–Laplace transform domain can be used to analyze the error in the CQ approximation in time. Suppose there is a σ0 > 0 such that for all s with (s) > σ0 the transformed kernel kˆ is analytic and satisfies the bound



ˆk(s)  C (σ0 )|s|μ

(2.8)

for some μ > 0, then it turns out that a stable CQ scheme is obtained if the multistep method is A-stable and if no poles on the unit circle |ζ | = 1 [27]. Suppose the exact solution is

γ (ζ ) has

u = K (∂t ) g whereas the approximate CQ solution is





u t = K ∂tt g where the Fourier–Laplace transform of the kernel denoted kˆ satisfies (2.8). Let

en = u (tn ) − u t (tn )

where tn = nt .

Choose, in addition, a final time T = N t > 0. Under these conditions Lubich [27] proves the following general result: Theorem 2.7. Let g ∈ H r ((−∞, ∞)) with g (t ) = 0 for t  0 with r > 1/2 + max(μ, 0). If a pth order A-stable multistep method p with γ having no poles on |ζ | = 1 is used to define K (∂tt ) and if β = min((r − μ) p +1 , r + 1, p ) then

 t

N 

1/2 |en |

2

 C ( T )hβ g H r .

n =0

Remark 2.8. Note that in our application in time.

μ  0 so if the data is such that r = μ + p + 1 we have the full order of convergence

For an application of this theorem to acoustic wave propagation, see Lubich [27] and for electromagnetism see [15,17]. When proving error estimates for a boundary element method including spatial discretization, we must prove the (2.8) for the spatially discrete convolution operator [27]. 3. Laplace domain regularized CFIE-I (RCFIE-I) When Γ is the surface of a bounded Lipschitz polyhedron, we still want to use a CFIE type operator to solve E and avoid possible instabilities caused by resonance frequencies in the frequency domain. Noting in (2.5), the operator V 02 works as a regularizing operator, which makes V (s) V 02 compact. This suggests using regularizing operator in the CFIE operator. In this section, we follow instead the choice of [12] who analyzed a regularized frequency domain CFIE. Following [12], we define M to be a compact operator −1/2

M : H

−1/2

(divΓ , Γ ) → H

−1/2

satisfying, for ξ ∈ H

(divΓ , Γ ),

(divΓ , Γ )

70

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

  M ξ, ξ  :=

M ξ(x) × n(x) · ξ (x) dσx  0, Γ

with the equality if and only if ξ = 0. Here the overbar “ ¯· ” denotes the complex conjugate. Then applying this regularizing operator M within the CFIE operator, we can define a new regularized Fourier–Laplace domain CFIE operator, called RCFIE-I and denoted C s, R , for x ∈ Γ by



C s, R ϕˆ (x) := sV (s)ϕˆ (x) +

1

K (s) +

2

I n M ϕˆ (x).

Correspondingly, the solution E of problem (2.1) has the ansatz, for x ∈ /Γ,





Φs (x, y)ϕˆ (y) dσ y − ∇

E(x) = s2 Γ



Φs (x, y) divΓ ϕˆ (y) dσ y + curl Γ

Φs (x, y) M ϕˆ (y) dσ y , Γ

and, by taking traces on Γ and using the jump conditions, ϕˆ

−1/2 ∈ H (divΓ , Γ ) satisfies, for x ∈ Γ

C s, R ϕˆ (x) = n(x) × g(x).

(3.1)

Two concrete examples of the compact operator M are: −1/2

1. The operator M from [12] which is defined such that if ζ ∈ H

solution of −1/2

w, q + divΓ w, divΓ q = ζ, q for all q ∈ H

−1/2

(divΓ , Γ ) then M ζ = w ∈ H

(divΓ , Γ ) is the unique

(divΓ , Γ ).

Buffa and Hiptmair [12, p. 630] show that this operator is coercive and positive. 2. Another possible choice that keeps entirely in a boundary integral framework is the operator M = Π T S 02 used by [18] (see p. 193) where S 0 is the diagonal operator using the single layer potential for the Helmholtz equation on each component S 0 : H −1/2 (Γ ) → H 1/2 (Γ ) defined by



S 0 ξ(x) = Γ

ξ(y) d A (y). 4π x − y

From [28, discussion on p. 209] we see that S 02 : H −1/2 (Γ ) → H 1 (Γ ) and note that L 2 (Ω) is compactly embedded in H −1/2 (Γ ) so Π T S 02 is indeed compact. We now analyze the regularized scheme. −1/2

Lemma 3.1. Given f ∈ H⊥

−1/2

(curlΓ , Γ ), if there exists ϕˆ ∈ H

−1/2

Proof. Suppose f = 0, and let ϕˆ ∈ H



C s, R ϕˆ = sV (s)ϕˆ +

K (s) +

1 2

(divΓ , Γ ) satisfy

I n M ϕˆ = 0.

Define, for x ∈ /Γ,





 Φs (x, y)ϕˆ (y) dσ y − ∇

u(x) = s2

(divΓ , Γ ) such that C s, R ϕˆ = f, then ϕˆ is unique.

Γ

Φs (x, y) divΓ ϕˆ (y) dσ y + curl Γ

Φs (x, y) M ϕˆ (y) dσ y .

(3.2)

Γ

Using the jump and continuity results in Theorem 2.3, we derive that u satisfies

in Ω ∪ Ω + ,

curl curl u + s2 u = 0,

γT+ u = 0, on Γ. For the exterior part, taking the dot product of the Maxwell equation above by the test function u and using integration by parts, we get



Ω+



curl u(x) 2 dx + s2



Ω+



u(x) 2 dx =



Γ

γN+ u(x) · γT+ u(x) dσx = 0.

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

71



When (s) = 0, the imaginary part of the left hand side gives Ω + |u(x)|2 dx = 0, i.e. u = 0 in Ω + . When (s) = 0, the real part gives the same conclusion. By the definition of u and jump conditions for the vector potentials, we have, on the boundary Γ ,

[γt u] = M ϕˆ ,

[γ N u] = s2 ϕˆ ,

then, using the fact that u = 0 on Ω + we have

γt− u = M ϕˆ ,

γN− u = s2 ϕˆ .

Now, using s2 u as the test function for the above Maxwell system, and integrating by parts, we get



s2



curl u(x) 2 dx + |s|4

Ω







u(x) 2 dx = −|s|4

Ω

ϕˆ (x) · M ϕˆ (x) × n(x) dσx . Γ

For equality, comparing the imaginary part when (s) = 0 and the real part when (s) = 0, we obtain  the above 2 | curl u ( x )| dx = 0. Then comparing the real part, by the definition of the regularizing operator M, we conclude that Ω u = 0 in Ω and ϕˆ = 0 on Γ , which finish the proof. 2 −1/2

Lemma 3.2. The RCFIE-I operator C s, R : H

−1/2

(divΓ , Γ ) → H⊥

(curlΓ , Γ ) is invertible.

Proof. We can derive Gårding’s inequality in which the principle part is given by sV (0) [12] and the compact perturbation is given by s( V (s) − V (0)) and ( K (s) + (1/2) I n ) M (it is here that the compactness of M is used), which implies the invertibility of C s, R by the injectivity given in the proceeding lemma. 2 Lemma 3.2 gives the existence and uniqueness of solutions of the boundary integral equation (3.1). This is not enough for our purpose of analyzing time domain solution error using corresponding Laplace domain operator estimates. We need explicit bounds with respect to s of the inverse operator. We now prove some auxiliary lemmas. −1/2

Lemma 3.3. For a Lipschitz polyhedron Ω with boundary Γ , given g ∈ H

curl curl u + s2 u = 0, u × n = g,

(divΓ , Γ ), the following problem

in Ω ∪ Ω + ,

(3.3a)

on Γ,

(3.3b)

has a unique solution u ∈ H(curl, Ω ∪ Ω + ), and

u Hs (curl,Ω∪Ω + )  c (Γ, σ0 )|s|2 g H−1/2 (div

Γ ,Γ )



.

Remark 3.4. In some cases, we may like to have the boundary data in the form n ×(u × n) = n × g. With a slight modification of the proof and using Theorem 2.1 for γ T± , we have

u Hs (curl,Ω∪Ω + )  c (Γ, σ0 )|s|2 n × g H−1/2 (curl

Γ ,Γ )

⊥

.

Proof. The proof follows the ideas of Terrasse [31] in the case of a smooth domain, and is extended to a Lipschitz domain in [15,17]. For completeness, we shall give a sketch of the proof here. First, consider the interior part,

curl curl u− + s2 u− = 0, −

u × n = g,

in Ω,

on Γ.

From Theorem 2.1, there exists u0 ∈ H(curl, Ω) such that

γt− u0 = g, and u0 Hs (curl,Ω)  c (Γ, σ0 )|s| g H−1/2 (div

consider the inhomogeneous boundary value problem

˜ + s2 u˜ = −curl curl u0 − s2 u0 , curl curl u ˜ × n = 0, u

in Ω,

on Γ.

Define the sesquilinear form



aΩ (w, v) :=

 s2 w(x) · v(x) dx.

curl w(x) · curl v(x) dx + Ω

Ω



Γ ,Γ )

. Then

72

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

Then it is easy to verify the continuity and coercivity (by multiplying by s)



aΩ (w, v)  w H (curl,Ω) v H (curl,Ω) , s



aΩ (w, w)  s w H (curl,Ω) . s |s|

s

˜ and hence Using the continuity and coercivity of aΩ , the Lax–Milgram lemma proves the existence and uniqueness of u, ˜ Furthermore u− = u.

 − u 

Hs (curl,Ω)

 c (Γ, σ0 )|s|2 g H−1/2 (div

Γ ,Γ )



.

A similar argument gives the corresponding estimate for the exterior part. These complete the proof.

2

In the following lemma, we shall give an inequality related to a bound on the map from Neumann data to Dirichlet data. −1/2

Lemma 3.5. Given a Lipschitz polyhedron Ω with boundary Γ and g ∈ H

of (3.3). Then

  [γ N u]

−1/2

H

(divΓ ,Γ )

 c (Γ, σ0 )|s|3 g H−1/2 (div

Γ ,Γ )



(divΓ , Γ ), let u ∈ H(curl, Ω ∪ Ω + ) be the solution

.

Proof. Let u− := u|Ω , then curl u− ∈ L2 (Ω) and curl curl u− = −s2 u− ∈ L2 (Ω), so curl u ∈ H(curl, Ω). Similarly, let u+ := u|Ω + , then curl u+ ∈ H(curl, Ω + ). Applying γt− and γt+ to u− and u+ respectively, and using Theorem 2.1, we have γ N± u± = 1/ 2 γt± curl u± ∈ H− (divΓ , Γ ), and

 ± ±   γ u  −1/2 = γt± curl u± H−1/2 (div ,Γ )  c (Γ, σ0 ) curl u Hs (curl,Ω ± ) N H

(divΓ ,Γ ) Γ

 ±    c (Γ, σ )|s| u H (curl,Ω ± )  c (Γ, σ0 )|s| u Hs (curl,Ω∪Ω + ) . s

Combining with the result in Lemma 3.3,

u Hs (curl,Ω∪Ω + )  c (Γ, σ0 )|s|2 g H−1/2 (div

Γ ,Γ )

and use of this inequality above finishes the proof.

,

2

We now can obtain a bound on C s−, R1 as follows: −1/2

Theorem 3.6. For a Lipschitz polyhedron Ω with boundary Γ , the RCFIE-I operator C s, R : H

invertible, and the inverse is bounded as follows

−1/2

(divΓ , Γ ) → H⊥

(curlΓ , Γ ) is

 −1  C   c (Γ, σ0 )|s|2 . s, R

−1/2

Proof. Given g ∈ H

(divΓ , Γ ), define u in Ω ∪ Ω + by the ansatz (3.2), then u satisfies

[γ N u] = s2 ϕˆ , and

curl curl u + s2 u = 0,

in Ω ∪ Ω + ,

γT+ u = C s, R ϕˆ = n × g, on Γ, s2 [γt u] − M [γ N u] = 0,

on Γ.

Here the last equality is from Theorem 2.3. Using the same argument as in [15,16], the equivalence between the above boundary value problem and the boundary integral equation C s, R ϕˆ = n × g can be verified. By Lemmas 3.3 and 3.5, the following inequality holds

 +  γ u N

−1/2

H

(divΓ ,Γ )

 c (Γ, σ0 )|s|3 n × g H−1/2 (curl ⊥

Γ ,Γ )

(3.4)

. −1/2

Define f := s2 γt+ u − M γ N+ u, then n × f = s2 n × g − n × M γ N+ u ∈ H⊥ u ∈ H(curl, Ω) such that

(curlΓ , Γ ). Consider the interior problem of finding

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

curl curl u + s2 u = 0, s

2

−

73

in Ω,

−

γT u − n × M γN u = n × f, on Γ.

(3.5)

From the proof of Lemma 3.5, we have

 −  γ u N

−1/2

H

 c (Γ, σ0 )|s| u Hs (curl,Ω) .

(divΓ ,Γ )

(3.6)

For problem (3.5), we consider a mixed formulation by giving a new function v satisfying

curl u + sv = 0,

in Ω,

curl v − su = 0,

in Ω,

s

2

−

γT u − sn × M γt− v = n × f, on Γ.

Now define bilinear forms



a(v, w) := s

v(x) · w(x) dx + Ω



s

γt− w(x) · n × M γt− v(x) dσx ,

Γ

 b(u, w) :=

1

u(x) · curl w(x) dσx , Ω



c (u, x) := s

u(x) · z(x) dσx . Ω

Then problem (3.5) is equivalent to finding (u, v) ∈ H(curl, Ω) × H(curl, Ω) such that for all (w, z) ∈ H(curl, Ω) × H(curl, Ω)

1

a(v, w) + b(u, w) = −



s2

γt− w(x) · n(x) × f(x) dσx ,

Γ

b(z, v) − c (u, v) = 0. Letting w = v and z = u and considering the real part, we have

σ v 2L2 (Ω) + σ u 2L2 (Ω) 

1  −  γ v −1/2

n × f H−1/2 (curl ,Γ ) . H

(divΓ ,Γ ) Γ ⊥ | s |2 t

Noting that

 −   γ v t

−1/2

H

(divΓ ,Γ )

=

1 −  γ u −1/2  c (Γ, σ0 ) u Hs (curl,Ω) , |s| N H (divΓ ,Γ )

we have

1

v 2L2 (Ω)  c (Γ, σ0 )

| s |2

u Hs (curl,Ω) n × f H−1/2 (curl

Γ ,Γ )

⊥

,

and

u 2L2 (Ω)  c (Γ, σ0 )

1

| s |2

u Hs (curl,Ω) n × f H−1/2 (curl ⊥

Γ ,Γ )

.

By recalling that v = −(1/s) curl u, combining the above two inequalities gives

u Hs (curl,Ω)  c (Γ, σ0 ) n × f H−1/2 (curl

Γ ,Γ )

⊥

.

Recalling the definition of f and inequality (3.4), we have

n × f H−1/2 (curl ⊥

Γ ,Γ )

   s2 n × gH−1/2 (curl ⊥

Γ ,Γ )

  + n × M γ N+ uH−1/2 (curl

 c (Γ, σ0 )|s| n × g H−1/2 (curl ⊥

Γ ,Γ )

⊥

3

Γ ,Γ )

.

So

 −  γ u N

−1/2

H

(divΓ ,Γ )

 c (Γ, σ0 )|s| u Hs (curl,Ω)  c (Γ, σ0 )|s|4 n × g H−1/2 (curl ⊥

Γ ,Γ )

.

74

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

Finally

ϕˆ H−1/2 (div

Γ ,Γ )



=

1  

 [γ N u]H−1/2 (div

| s |2

Γ ,Γ )



 c (Γ, σ0 )|s|2 n × g H−1/2 (curl ⊥

Γ ,Γ )

,

2

which completes the proof.

4. Laplace domain regularized CFIE-II (RCFIE-II) The previous RCFIE is a direct generalization of previous results for the frequency domain, but the proof of the s dependence of the norm bounds for the operator do not extend easily to spatially discrete operators, and we have been unable to prove convergence of the fully discrete scheme in this case (of course the time semi-discrete scheme converges via Lubich’s theory since we have verified the appropriate mapping properties). An alternative approach which does not use the usual frequency domain choice of weights between the electric and magnetic operators will be considered next. We can modify the previous RCFIE to obtain one that can be more easily analyzed, yet still does not suffer from resonances in the frequency domain. As for the acoustic case [16], we consider a new type of regularized CFIE operator by modifying the compact operator to include time dependence

˜ := M

M s4

,

then new RCFIE-II operator C˜ s, R is given by

˜C s, R ϕˆ := sV (s)ϕˆ + K (s) + 1 I n M ˜ ϕˆ . 2

To obtain the invertibility and boundedness of C˜ s, R , we need the following prerequisite lemmas. −1/2

Lemma 4.1. Given a Lipschitz polyhedron Ω with boundary Γ and ϕˆ ∈ H



Φs (x, y)ϕˆ (y) dσ y ,

u(x) := curl

(divΓ , Γ ), let

x∈ / Γ.

Γ

Then

 ±  γ u N

−1/2

H

(divΓ ,Γ )

 c (Γ, σ )|s|3 n × ϕˆ H−1/2 (curl ⊥

Γ ,Γ )

.

Proof. Noting u ∈ H(curl, Ω ∪ Ω + ), γ N± u = γt± curl u on Γ (by Theorem 2.3 γ N− u = γ N+ u, we shall use the proof), curl u ∈ H(curl, Ω ∪ Ω + ) and curl curl u + s2 u = 0. Then we define, for x ∈ Ω ∪ Ω + ,

γN u in the rest of



Φs (x, y)ϕˆ (y) dσ y

curl u(x) = curl curl Γ





Φs (x, y)ϕˆ (y) dσ y + ∇

= −s2 Γ

Φs (x, y) divΓ ϕˆ (y) dσ y Γ

=: w(x).

(4.1)

Using the jump properties of the vector potential we see that w satisfies

curl curl w + s2 w = 0,

[γ N w] = −s2 ϕˆ , [γt w] = 0,

in Ω ∪ Ω + ,

on Γ,

on Γ.

Using the test function sw and integration by parts we obtain



s



curl w(x) 2 dx + |s|2 s

Ω∪Ω +





w(x) 2 dx = −|s|2 s

Ω∪Ω +

 n(x) × ϕˆ (x) · γt w(x) dσx , Γ

which implies

w 2Hs (curl,Ω∪Ω + )  c (σ0 )|s|3 n × ϕˆ H−1/2 (curl ⊥

Γ ,Γ )

γt w H−1/2 (div

Γ ,Γ )

.

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

75

By Theorem 2.1

γt w H−1/2 (div

Γ ,Γ )



 c (Γ, σ0 ) w Hs (curl,Ω∪Ω + ) ,

so

w Hs (curl,Ω∪Ω + )  c (Γ, σ0 )|s|3 n × ϕˆ H−1/2 (curl ⊥

Γ ,Γ )

.

Finally, we have

γ N u H−1/2 (div

Γ ,Γ )



= γt curl u H−1/2 (div

= γt w H−1/2 (div

Γ ,Γ )



Γ ,Γ )



 c (Γ, σ0 )|s|3 n × ϕˆ H−1/2 (div

Γ ,Γ )

⊥

.

2 −1/2

For convenience, we now define a potential operator ΨDL,s , for ϕˆ ∈ H



(divΓ , Γ ) and x ∈ / Γ , by

Φs (x, y)ϕˆ (y) dσ y .

ΨDL,s ϕˆ (x) := curl Γ

−1/2

Lemma 4.2. Let Ω be a Lipschitz polyhedron with boundary Γ . Then the operator ΨDL,s : H

linear continuous and

(divΓ , Γ ) → H(curl, Ω ∪ Ω + ) is

ΨDL,s  c (Γ, σ0 )|s|2 . −1/2

Proof. Given ϕˆ ∈ H

(divΓ , Γ ), let u = ΨDL,s ϕˆ , then using the jump conditions in Theorem 2.3 [γT u] = n × ϕˆ . Using su as a test function and integration by parts on Ω ∪ Ω + , we get



s



curl u(x) 2 dx + |s|2 s

Ω∪Ω +







u(x) 2 dx = s

Ω∪Ω +

γN u(x) · n(x) × ϕˆ (x) dσx , Γ

which gives

u 2Hs (curl,Ω∪Ω + )  c (Γ, σ0 )|s| γ N u H−1/2 (div

Γ ,Γ )



4

n × ϕˆ H−1/2 (curl

Γ ,Γ )

⊥

ˆ 2 −1/2 , H⊥ (curlΓ ,Γ )

 c (Γ, σ0 )|s| n × ϕ

where the last inequality is from the preceding lemma. This proves the desired inequality. −1/2

Corollary 4.3. The operator K (s) : H

   K (s)  c (Γ, σ0 )|s|2 . −1/2

Proof. Given ϕˆ ∈ H



−1/2

(divΓ , Γ ) → H⊥

(curlΓ , Γ ) is bounded,

(divΓ , Γ ), we have, on Γ ,

1 I n + K (s) ϕˆ = γ T− ΨDL,s ϕˆ . 2

Then

   K (s)ϕˆ 

−1/2

H

1

(divΓ ,Γ )

 ϕˆ H−1/2 (div 2



Γ ,Γ )

  + γT−  ΨDL,s

ϕˆ H−1/2 (div

Γ ,Γ )

.

By the remark following Theorem 2.1 and the previous lemma, the inequality holds

   K (s)ϕˆ 

−1/2

H

(divΓ ,Γ )

which finishes the proof.

 c (Γ, σ0 )|s|2 ϕˆ H−1/2 (div



1



Φs (x, y)ϕˆ (y) dσ y − ∇Γ

Φs (x, y) divΓ ϕˆ (y) dσ y ,

s

Γ

,

2

Recall that, for x ∈ Γ ,

V (s)ϕˆ (x) = sΠ T

Γ ,Γ )



Γ

2

76

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

−1/2

and define the bilinear form on H

 

−1/2

(divΓ , Γ ) × H

(divΓ , Γ ),

 

s2 Φs (x, y)ϕˆ (y) · ξ(x) dσ y dσx +

b s (sϕˆ , ξ ) := Γ Γ

Φs (x, y) divΓ ϕˆ (y) divΓ ξ(x) dσ y dσx . Γ Γ

By results in [15,17], for ϕˆ

−1/2 ∈ H (divΓ , Γ ),





sV (s)ϕˆ , ϕˆ = b s (sϕˆ , ϕˆ )  c (Γ, σ0 ) s ϕˆ 2 −1/2 . (divΓ ,Γ ) H

|s| We finally obtain the following bound of the inverse of RCFIE-II operator C˜ s, R . −1/2

Theorem 4.4. For large enough σ0 (depending on Γ ), and s > σ0 , the RCFIE-II operator C˜ s, R : H

is invertible, and

−1/2

(divΓ , Γ ) → H⊥

(curlΓ , Γ )

 −1  C˜   c (Γ, σ0 )|s|. s, R

−1/2

Proof. Given ϕˆ ∈ H

(divΓ , Γ ),

 





 

C˜ s, R ϕˆ , ϕˆ  = sV (s)ϕˆ , ϕˆ + 1 I n + K (s) M ˜ ϕˆ , ϕˆ



2

 



  1 ˜ ϕˆ , ϕˆ

 sV (s)ϕˆ , ϕˆ −

I n + K (s) M

2

s 1  c (Γ, σ0 ) ϕˆ 2 −1/2 − c (Γ, σ0 ) 2 ϕˆ 2 −1/2 H

H

(divΓ ,Γ ) (divΓ ,Γ ) |s| |s| s  c  (Γ, σ0 ) ϕˆ 2 −1/2 , (divΓ ,Γ ) H

|s|

(4.2)

where the second last inequality is from the preceding corollary, and the last inequality holds for σ0 large enough since s > σ0 . This shows the coercivity of C˜ s, R , and the continuity of C˜ s, R follows from the continuity of the operators V (s), K (s) and I n . Then the Lax–Milgram lemma gives the invertibility of C˜ s, R . Moreover, the above inequality implies the inverse bound. 2 5. Discretization of RCFIE-II operator In this section, we shall outline the discretization of the RCFIE-II operator in the Fourier–Laplace domain. For space dis−1/2 (divΓ , Γ ). cretization, we use Raviart–Thomas elements [5] of degree k, denoted RTk (Γh ) ⊂ H(divΓ , Γ ) to approximate H

Here Γh denotes a non-degenerate regular and quasi-uniform triangulation of Γ with mesh size h > 0 [2]. Raviart– Thomas elements are then defined as follows [5]: On each triangular element K ∈ Γh , uh ∈ RTk (Γh ) is such that uh | K ∈ Pk ( K ) + xpk ( K ), where pk ( K ) denotes polynomials of degree k on Γ and Pk ( K ) = ( pk ( K ))3 . These elements have the property that uh · n ∈ R k (∂ K ), where







R k (∂ K ) := p ∈ L 2 (∂ K ) p |e i ∈ pk (e i ), for each edge e i , i = 1, 2, 3 of K . The degrees of freedom are given by

⎧ ⎪ uh (x) · n(x) p (x) dσx , ∀ p ∈ R k (∂ K ), ⎪ ⎪ ⎪ ⎨ ∂K  ⎪ ⎪ ⎪ uh (x) · P(x) dx, ∀P ∈ Pk−1 ( K ). ⎪ ⎩ K

In engineering, a particular version of these elements called GWP elements [22] is used. These are equivalent to RTk (Γh ) with a special choice of degrees of freedom and a certain explicit basis. We now consider the discretized Fourier–Laplace domain problem for the RCFIE-II operator in variational form: find ϕˆ h ∈ RTk (Γh ) satisfying, for any ξh ∈ RTk (Γh ),

C˜ s, R ϕˆ h , ξh  = n × g, ξh .

(5.1)

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

77

−1/2

Lemma 5.1. Given g ∈ H

(divΓ , Γ ), the discrete Galerkin variational problem (5.1) has a unique solution ϕˆ h ∈ RTk (Γh ) with the quasi-optimal error estimate

ϕˆ − ϕˆ h H−1/2 (div

Γ ,Γ )



 c (Γ, σ0 )|s|2

inf

ξh ∈RTk (Γh )

ϕˆ − ξh H−1/2 (div

Γ ,Γ )



(5.2)

,

and it is stable

ϕˆ h H−1/2 (div

 c (Γ, σ0 )|s| n × g H−1/2 (curl

Γ ,Γ )



Γ ,Γ )

⊥

.

Here the generic constant c is independent of s, h, and ϕˆ is the solution of C˜ s, R ϕˆ = n × g. Proof. By the coercivity of C˜ s, R in Theorem 4.4 for

1/ 2 σ0 large enough, we have ∀ϕˆ h ∈ RTk (Γh ) ⊂ H− (divΓ , Γ ),





C˜ s, R ϕˆ h , ϕˆ h   c (Γ, σ0 ) s ϕˆ h 2 −1/2 . H

(divΓ ,Γ ) |s| −1/2

The continuity of C˜ s, R is given as, for any ϕˆ h and ξh in RTk (Γh ) ∈ H

(divΓ , Γ ),

 





 

C˜ s, R ϕˆ h , ξh   sV (s)ϕˆ h , ξh + 1 I n + K (s) M ˜ ϕˆ h , ξh



2  c (Γ, σ0 )|s| ϕˆ h H−1/2 (div

Γ ,Γ )

ξh H−1/2 (div

Γ ,Γ )

+ c (Γ, σ0 )

1

| s |2

ϕˆ h H−1/2 (div

Γ ,Γ )



ξh H−1/2 (div

Γ ,Γ )

.

Here the first part of the last inequality is from the proof of Theorem 2.6, and the second part of the last inequality is from Corollary 4.3. Then Céa’s lemma shows that

ϕˆ − ϕˆ h H−1/2 (div

Γ ,Γ )



 c (Γ, σ0 )|s|2

inf

ξh ∈RTk (Γh )

ϕˆ − ξh H−1/2 (div

Γ ,Γ )

.

To prove the last part of the lemma, stability is given by

  





s 2



˜ ˆ ˆ ˆ ˆ c (Γ, σ0 )  C s, R ϕh , ϕh  

n(x) × g(x) · ϕh (x) dσx

.

ϕh −1/2 H

(divΓ ,Γ ) |s|

2

Γ

We now have a spatial discretization error theorem by using the above theorem and approximation results in [10]. −1/2

Theorem 5.2. Suppose Ω is a Lipschitz polyhedron with boundary Γ . Given g ∈ H

be the solution of

(divΓ , Γ ), let ϕˆ ∈ Hα (divΓ , Γ ) for some α > 0

−1/2 C˜ s, R ϕˆ , ξ  = n × g, ξ  for all ξ ∈ H (divΓ , Γ ).

Let ϕˆ h ∈ RTk (Γh ) solve problem (5.1), then for any ε > 0, and h  h0 for some h0 > 0, the error is bounded by

ϕˆ − ϕˆ h H−1/2 (div

Γ ,Γ )



 c (Γ, σ0 )hβ |s|2 ϕˆ Hα (divΓ ,Γ ) ,

where β = min{ 32 + k − ε , α +

1 2

− ε , 1 + k + m∗ , α + m∗ }, and m∗ is a constant depending on Γ (see [10] for details).

We have focused on the analysis of CFIE operators in the Fourier–Laplace domain, which then forms the basis for the analysis of the convolution quadrature method applied to the time domain integral equations using Raviart–Thomas elements in space. Because the Fourier–Laplace domain RCFIE-II operator is coercive, the spatially discrete Fourier–Laplace domain RCFIE-II operator satisfies (4.2) and hence has the same s-dependent bound as for the continuous case. It can then be analyzed using Theorem 2.7. This analysis is the direct generalization of the analysis for the acoustic case in the original convolution quadrature paper [27] (see also a similar analysis in electromagnetic case [15,17]). Theorem 5.3. Given the conditions of the previous theorem and boundary data G . Let ψ be a sufficiently smooth solution of C˜ R ψ = n × G , where C˜ R is the time domain RCFIE-II operator corresponding to C˜ s, R . Let ψht be the corresponding full (time and space) discretized solution of the CQ-RTk discretized RCFIE-II using an order p multistep method satisfying the conditions of Theorem 2.7 then

  ψ(·, nt ) − ψ nt 2 −1/2 h

H

(Γ )

  2     + divΓ ψ(·, nt ) − ψhnt  H −1/2 = O t 2p + O h2β .

This verifies the optimal convergence in space and time for the RCFIE-II method provided the solution is smooth enough.

78

Q. Chen, P. Monk / Applied Numerical Mathematics 79 (2014) 62–78

6. Conclusion We have suggested a regularized TD-CFIE suitable for discretization by CQ and Raviart–Thomas finite elements (the RCFIE-II boundary integral equation). We have proved convergence of the fully discrete problem in this case. Clearly a useful next step would be to perform numerical tests of both the regularized and the unregularized CFIE formulations to determine if regularization has a positive effect on the computed solution, and to compare the results to those obtained by an electric field integral equation approach. Acknowledgement The author Q.C. would gratefully thank Dr. L. Banjai for the helpful discussion of the bounds of inverse V (s). References [1] T. Abboud, J.-C. Nédélec, J. Volakis, Stable solution of the retarded potential equations, in: Proc. 17th Ann. Rev. Progress in Appl. Comp. Electromagnetics, Monterey, CA, 2001, pp. 146–151. [2] K. Atkinson, W. Han, Theoretical Numerical Analysis: A Functional Analysis Framework, second ed., Springer, 2005. [3] A. Bamberger, T. HaDuong, Formulation variationnelle espace–temps pour le calcul par potentiel retardé de la diffraction d’une onde acoustique (i), Math. Methods Appl. Sci. 8 (1986) 405–435. [4] L. Banjai, Multistep and multistage convolution quadrature for the wave equation: algorithms and experiments, SIAM J. Sci. Comput. 32 (2010) 2964–2994. [5] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, 1991. [6] A. Buffa, S. Christansen, The electric field integral equation on Lipschitz screens: definitions and numerical approximation, Numer. Math. 94 (2003) 229–267. [7] A. Buffa, P. Ciarlet, On traces for functional spaces related to Maxwell’s equations. Part I: An integration by parts formula in Lipschitz polyhedra, Math. Methods Appl. Sci. 24 (2001) 9–30. [8] A. Buffa, P. Ciarlet, On traces for functional spaces related to Maxwell’s equations. Part II: Hodge decompositions on the boundary of Lipschitz polyhedra and applications, Math. Methods Appl. Sci. 24 (2001) 31–48. [9] A. Buffa, M. Costabel, C. Schwab, Boundary element methods for Maxwell’s equations on non-smooth domain, Numer. Math. 92 (2002) 679–710. [10] A. Buffa, M. Costabel, D. Sheen, On traces for H (curl, Ω) in Lipschitz domains, J. Math. Anal. Appl. 276 (2002) 845–867. [11] A. Buffa, R. Hiptmair, Galerkin boundary element methods for electromagnetic scattering, in: T.J. Barth, M. Griebel, D.E. Keyes, R.M. Nieminen, D. Roose, T. Schlick, M. Ainsworth, P. Davies, D. Duncan, B. Rynne, P. Martin (Eds.), Topics in Computational Wave Propagation, in: Lect. Notes Comput. Sci. Eng., vol. 31, Springer, Berlin, Heidelberg, 2003, pp. 83–124. [12] A. Buffa, R. Hiptmair, A coercive combined field integral equation for electromagnetic scattering, SIAM J. Numer. Anal. 42 (2004) 621–640. [13] A. Buffa, R. Hiptmair, Regularized combined field integral equations, Numer. Math. 100 (2005) 1–19. [14] M. Cessenat, Mathematical Methods in Electromagnetism: Linear Theory and Applications, World Scientific, 1996. [15] Q. Chen, Convolution quadrature applied to time domain acoustic and electromagnetic scattering problems, PhD thesis, University of Delaware, 2011. [16] Q. Chen, P. Monk, Discretization of the time domain CFIE for acoustic scattering problems using convolution quadrature, submitted for publication. [17] Q. Chen, P. Monk, X. Wang, D. Weile, Analysis of convolution quadrature applied to the time electric field integral equation, Commun. Comput. Phys. 11 (2) (2012) 383–399. [18] D. Colton, R. Kress, Inverse Acoustic and Electromagnetic Scattering, second ed., Springer, 1998. [19] K. Cools, F. Andriulli, F. Olyslager, E. Michielssen, Time domain Calderon identities and their application to the integral equation analysis of scattering by PEC objects. Part I: Preconditioning, IEEE Trans. Antennas and Propagation 57 (2009) 2352–2364. [20] M. Costabel, Time-dependent problems with the boundary integral equation method, in: Encyclopedia of Computational Mechanics, vol. 1, Fundamentals, John Wiley & Sons, 2004, pp. 703–721. [21] A. Ergin, B. Shanker, E. Michielssen, Analysis of transient wave scattering from rigid bodies using a Burton–Miller approach, J. Acoust. Soc. Am. 106 (5) (1999) 2396–2404. [22] R. Graglia, D. Wilton, F. Peterson, Higher order interpolatory vector bases for computational electromagnetics, IEEE Trans. Antennas and Propagation 45 (3) (1997) 329–342. [23] W. Hackbusch, W. Kress, S. Sauter, Sparse convolution quadrature for time domain boundary integral formulations of the wave equation, IMA J. Numer. Anal. 29 (2009) 158–179. [24] T. HaDuong, On retarded potential boundary integral equations and their discretisation, in: M. Ainsworth, P. Davies, D. Duncan, P. Martin, B. Rynne (Eds.), Topics in Computational Wave Propagation: Direct and Inverse Problem, Springer-Verlag, Berlin, 2003, pp. 301–336. [25] W. Kress, S. Sauter, Numerical treatment of retarded boundary integral equations by sparse panel clustering, IMA J. Numer. Anal. 28 (2008) 162–185. [26] A.R. Laliena, F.J. Sayas, Theoretical aspects of the application of convolution quadrature to scattering of acoustic waves, Numer. Math. 112 (2009) 637–678. [27] C. Lubich, On the multistep time discretization of linear initial–boundary value problems and their boundary integral equations, Numer. Math. 67 (1994) 365–389. [28] W. McLean, Strongly Elliptic Systems and Boundary Integral Equations, Cambridge University Press, Cambridge, 2000. [29] J. Nédélec, Acoustic and Electromagnetic Equations, Springer, 2000. [30] M. Pocock, M. Bluck, S. Walker, Electromagnetic scattering from 3-D curved dielectric bodies using time-domain integral equations, IEEE Trans. Antennas and Propagation 46 (1998) 1212–1219. [31] I. Terrasse, Resolution mathematique et numerique de equations de Maxwell instationnaires par une methode de potentiels retardes, PhD thesis, Ecole Polytechnique, 1993. [32] A. Veit, Numerical methods for time-domain boundary integral equations, PhD thesis, Universitat Zurich, 2011. [33] X. Wang, D. Weile, Implicit Runge–Kutta methods for the discretization of time domain integral equations, IEEE Trans. Antennas and Propagation 29 (2011) 4651–4663. [34] X. Wang, R. Wildman, D. Weile, P. Monk, A finite difference delay modeling approach to the discretization of the time domain integral equations of electromagnetics, IEEE Trans. Antennas and Propagation 56 (2008) 2442–2452.