A new coupling of mixed finite element and ... - Semantic Scholar
A new coupling of mixed finite element and boundary element methods for an exterior Helmholtz problem in the plane∗ ´rquez‡, Salim Meddahi§ Gabriel N. Gatica†, Antonio Ma
Abstract This paper deals with the scattering of time harmonic electromagnetic waves by an infinitely long cylinder containing a non-homogeneous conducting medium. More precisely, we study the transverse magnetic field that solves an interface problem holding between the cross section of the cylinder and the exterior two-dimensional free space. We apply a dual-mixed variational formulation in the obstacle coupled with a boundary integral equation method in the unbounded homogeneous space. A Fredholm alternative is utilized to prove that this continuous formulation is well posed. We define the corresponding discrete scheme by using the lowest order rotated Raviart-Thomas finite elements for the magnetic field and spectral elements for the boundary unknown. Then, we show that the resulting Galerkin scheme is uniquely solvable and convergent, and prove optimal error estimates. Finally, we illustrate our analysis with some results from computational experiments. Key words. mixed-FEM, BEM, coupling, Helmholtz equation. Mathematics subject classifications (1991). 65N30, 65N38, 65N15, 35Q60
1
Introduction
In this paper we consider the exterior Helmholtz equation that arises from a time harmonic electromagnetic scattering problem posed in the plane. The scatterer is a heterogeneous conducting body embedded in an infinite, homogeneous, isotropic, and non-conducting medium. Since the underlying partial differential equation is coercive, we observe that the treatment of the usual primal variational formulation is a rather straightforward application of the Fredholm alternative (see, e.g. [4, 9] and the references therein for more details). On the other hand, it is well known (see, e.g. [2, 11]) that a dual-mixed variational formulation is more adequate to obtain a good approximation of the vectorial unknown, namely the magnetic field in the conductor. A procedure in this direction is developed recently in [6], where the coupling of mixed finite element (mixed-FEM) and boundary element methods (BEM) is employed for the corresponding discrete scheme. In fact, the approach in [6] uses a dual-mixed variational formulation in a bounded domain delimited by a C ∞ closed curve, which is provided with the non-local boundary conditions ∗
This research was partially supported by CONICYT-Chile through the FONDAP Program in Applied Mathematics, by the Direcci´ on de Investigaci´ on of the Universidad de Concepci´ on through the Advanced Research Groups Program, and by the Ministerio de Educaci´ on y Ciencia of Spain, through the project No. MTM2004-05417. † GI2 MA, Departamento de Ingenier´ıa Matem´ atica, Facultad de Ciencias F´ısicas y Matem´ aticas, Universidad de Concepci´ on, Casilla 160-C, Concepci´ on, Chile, e-mail: [email protected] ‡ Departamento de Construcci´ on e Ingenier´ıa de Fabricaci´ on, Universidad de Oviedo, Oviedo, Espa˜ na, e-mail: [email protected] § Departamento de Matem´ aticas, Facultad de Ciencias, Universidad de Oviedo, Calvo Sotelo s/n, Oviedo, Espa˜ na, e-mail: [email protected]
1
arising from the application of the boundary integral equation method in the unbounded region. The corresponding analysis considers first the auxiliary problem obtained from the original formulation after substituting all the acoustic boundary integral operators by the corresponding potential boundary integral operators. Then, the properties of these operators and suitable duality arguments allow to prove a global inf-sup condition that implies the well posedness of the auxiliary problem. Afterwards, the Fredholm alternative is applied to the original problem to conclude that it is also well posed. Nevertheless, a drawback of the resulting coupled method is the fact that it is not well-posed if the square of the wave number is a Neumann eigenvalue of the Laplace operator in the interior domain. Thus, the corresponding numerical method may exhibit an unstable behavior in the neighborhood of these singular coefficients. Furthermore, the anonymous referee of this paper suggested us that eliminating the scalar variable, which constitutes the Lagrange multiplier in [6], would lead to a variational formulation whose structure agrees with the abstract framework provided by Buffa in [3] for a certain type of non-coercive problems. Consequently, the purpose of the present paper is to follow the above suggestion and develop a new coupling of mixed finite element and boundary element methods for the exterior Helmholtz equation described before. More precisely, we show that the powerful tool from [3] simplifies considerably the analysis presented in [6] for both the continuous and the discrete problem. In addition, the coupling of mixed-FEM and BEM is modified here in such a way that, on the contrary to [6], the discrete scheme now becomes free from spurious solutions. The rest of the paper is organized as follows. In Sections 2, 3 and 4 we introduce the model problem, derive our coupled formulation and prove its well-posedness. In Section 5 we introduce the Galerkin mixed-FEM/BEM discretization and provide the associated convergence result. Finally, Section 6 is devoted to the numerical results. √ In the sequel, we deal with complex valued functions and the symbol ı is used for −1. We let α and |α| be the conjugate and the modulus, respectively, of a complex number α ∈ C. We use boldface letters to denote vectors, vector-valued functions or vector-valued function spaces. The symbol |·| is 2 X qi qi for all q := (q1 , q2 ) ∈ C2 . used to denote the 2-norm for vectors, that is |q| = q · q := i=1
2 In addition, Z given a bounded domain Ω ⊆ R with a sufficiently smooth boundaryqΓ, we denote f g for all f, g ∈ L2 (Ω). Thus, the norm in L2 (Ω) is given by kvk0,Ω := f, g := v, v . Also, Ω
for any non-negative integer m, H m (Ω) is the usual Hilbertian Sobolev space endowed with the norm k·km,Ω . We denote the corresponding seminorm |·|m,Ω ∞ be On the other hand, we also need to introduce periodic Sobolev spaces. To this end, we let C2π the space of 2π-periodic infinitely differentiable complex valued functions of a single variable, and for ∞ define its Fourier coefficients by each g ∈ C2π 1 gb(k) := 2π
Z
2π
g(t) e−kıt dt
0
∀ k ∈ Z.
∞ with respect to the Then, given r ∈ R, the Hilbertian Sobolev space H r is the completion of C2π !1/2 X norm kgkr := . Moreover, the H 0 -bilinear form defined as h λ, η i := (1 + k2 )r |b g (k)|2
Z
2π
k∈Z
λ(t) η(t) dt can be extended to represent the duality between H −r and H r for all r > 0. We will
0
keep the same notation for this duality bracket. Throughout this paper, C denotes a generic constant independent of the discretization parameters, which usually takes different values at different places.
2
2
Statement of the problem
We consider an incident electromagnetic wave upon a body containing a non-homogeneous conducting medium and we are interested in determining the scattered electromagnetic wave. The obstacle is assumed to be an infinitely long cylinder parallel to the x3 -axis whose cross section is a bounded open domain Ωc . We denote by Ω′c its complement in R2 . Next, we let ǫ, µ and σ be the electric permittivity, the magnetic permeability, and the conductivity, respectively. They are assumed to be constant along each line parallel to the x3 -axis. These coefficients are bounded real valued scalar functions that satisfy almost everywhere in Ωc ǫ0 ≤ ǫ(x) ≤ ǫ1 ,
µ0 ≤ µ(x) ≤ µ1 ,
and
0 ≤ σ(x) ≤ σ1 ,
(2.1)
where ǫ0 , ǫ1 , µ0 , µ1 , and σ1 are positive constants, and ǫ0 and µ0 are the electric permittivity and magnetic permeability, respectively, of the free space. Moreover, since the medium is dielectric and homogeneous, outside the obstacle there holds ǫ(x) = ǫ0 ,
µ(x) = µ0 ,
and
in Ω′c .
σ(x) = 0
The incident electromagnetic field is supposed to exhibit a time-harmonic behavior with frequency ω and to be polarized perpendicular to the x3 -axis. Under these conditions, the electric and magnetic fields will also have a time harmonic behavior with frequency ω. Their amplitudes E and H, which are independent of x3 , are given by 0 h1 (x1 , x2 ) 1 and H = √1 h2 (x1 , x2 ) , 0 E= √ ǫ0 µ0 u(x1 , x2 ) 0
where the complex-valued functions u and h = (h1 , h2 )t are now the unknowns of the problem. Then, from the time dependent form of Maxwell’s equations, we know that u and h satisfy the equations curl h + ı k a u = 0 in R2 , (2.2) curl u − ı k b h = 0 in R2 ,
where
∂h2 ∂h1 − curl h := ∂x1 ∂x2
Here, k = ω a(x) =
√
and
curl u :=
ǫ0 µ0 is the wave number, b(x) =
∂u ∂u ,− ∂x2 ∂x1
t
.
µ(x) , and the index of refraction is given by µ0
σ(x) ǫ(x) +ı for all x ∈ R2 . Eliminating h from (2.2) yields ǫ0 ω ǫ0 curl (b−1 curl u) − k2 a u = 0
in R2 ,
u = ui + us in R2 ,
(2.3)
where ui and us denote the incident and scattered waves, respectively, and ui is given, for example, by the plane wave ui (x) = eı k x·d in the direction determined by the unit vector d. Furthermore, the behavior of the scattered electromagnetic field at infinity (known as the Silver-Muller condition) implies the following Sommerfeld radiation condition: s √ ∂u s r − ıku → 0 when r = |x| → ∞ . (2.4) ∂r 3
Note that our hypothesis on the physical coefficients of the problem imply that 1 ≤ b(x) ≤
µ1 µ0
for a.e. x ∈ R2
(2.5)
and a(x) := aR (x) + ı aI (x) is a bounded complex-valued function with 1 ≤ aR (x) ≤
ǫ1 ǫ0
and
0 ≤ aI (x) ≤
σ1 ω ǫ0
for a.e. x ∈ R2 .
(2.6)
Moreover, in the unbounded homogeneous region Ω′c we have b(x) = 1 and
a(x) = 1.
We recall that, in practice, aI = 0 means that no loss of energy occurs in the medium. Now, we follow [5] and make assumptions on the coefficients in order to guarantee uniqueness of solution for problem (2.3). For this purpose, we first introduce a sufficiently large C ∞ closed curve Γ such that its interior region, denoted from now on by Ω, is convex and contains the heterogeneous medium Ωc . The outward unit vector normal to Γ is represented by ν. Then, we assume that there exists a set {Ω0 , Ω1 , · · · , ΩN } of open and disjoint subdomains with Lipschitz-continuous boundaries such that Ω = Ω1 ∪ Ω2 ∪ · · · ∪ ΩN , Ω0 = R2 − Ω , and in each subdomain Ωi , i ∈ {0, 1, · · · , M }, the coefficients a and b satisfy one of the following two conditions: (C1) b and aR are functions in H 3 (Ωi ) and aI vanishes identically in Ωi ; or (C2) aI (x) ≥ a0 > 0 almost everywhere in Ωi . Theorem 2.1. Assume that ui = 0. Then problem (2.3) only admits the trivial solution. Proof. We refer to Theorem 3.1 of [5] for a proof of this result.
3
The coupling of mixed-FEM and BEM
Let x : R → R2 be a regular 2π-periodic parametric representation of the curve Γ, that is |x′ (s)| > 0 ∀s ∈ R and
x(s) = x(t)
⇔
t − s ∈ 2 π Z.
We define the parameterized trace on Γ as the linear continuous extension of γ : C ∞ (Ω) → L2 (0, 2π) u 7→ γ(u) := u|Γ ◦ x
(3.1)
to H 1 (Ω), which yields a bounded and onto linear application γ : H 1 (Ω) → H 1/2 (see Theorem 8.15 in [7]). Similarly, the parameterized tangential trace is defined as the linear continuous extension of γτ : [C ∞ (Ω)]2 → H −1/2 q 7→ γτ (q) := (q|Γ ◦ x) · x′
(3.2)
to H(curl; Ω), which yields a bounded and onto linear application γτ : H(curl; Ω) → H −1/2 . We recall here that H(curl; Ω) := q ∈ L2 (Ω) : curl q ∈ L2 (Ω) , 4
and its Hilbertian norm is given by kqk2H(curl;Ω) := kqk20,Ω + kcurl qk20,Ω . In order to deduce a mixed variational formulation of (2.3) we introduce the vectorial unknown p := b−1 curl u in Ω, so that b p := curl u in Ω . (3.3) Note that p is related to the magnetic field by p = ı k h. With this notation, the restriction of the first equation of (2.3) to Ω becomes curl p − k2 a u = 0
in
Ω.
(3.4)
Then, testing (3.3) with the vector-valued function q ∈ H(curl; Ω) and integrating by parts yields b p, q − u, curl q − h γτ (q), γ(u) i = 0. (3.5)
Now, we proceed differently from [6] and replace u = k−2 a−1 curl p from (3.4) into the second term of (3.5) to obtain the following variational formulation in the bounded domain Ω: Find p ∈ H(curl; Ω) such that b p, q − k−2 a−1 curl p, curl q − h γτ (q), γ(u) i = 0 ∀ q ∈ H(curl; Ω). (3.6)
Next, we employ a boundary integral equation method to complete (3.6) with the corresponding variational formulation in the unbounded exterior domain Ω0 , thus providing, in particular, a suitable expression for γ(u) (see (3.11) below). More precisely, given an arbitrary real number η > 0, we consider the representation formula Z 2π Z 2π ∂E ′ s E(x, x(t)) ψ(t) dt ∀ x ∈ Ω0 , (3.7) (x, x(t)) |x (t)| ψ(t) dt + ı η u (x) = ∂ν y 0 0 where E(x, y) :=
ı (1) H (k|x − y|) 4 0 (1)
is the radial outgoing fundamental solution of the Helmholtz equation, H0 is the Hankel function of the first kind and order zero, and ψ is an unknown density. Expression (3.7) is the parameterized version of the combined single-double layer potential sometimes referred to as a Brakhage-Werner potential (see [4, 13]). Note that the function |x′ (t)| appears only in one of the integrands, which aims to obtain parameterized integral equations without cumbersome weights (see [10] for more details). We remark that a function us given by (3.7) satisfies the Helmholtz equation ∆ us + k2 us = 0 in Ω0 and the sommerfeld radiation condition (2.4) for any 2 π-periodic density ψ : [0, 2π] → C. Thus, if we want to obtain a global solution for problem (2.3), we have just to impose the following transmission conditions on Γ Z 2π Z 2π ∂E E(x, x(t)) ψ(t) dt + , (x, x(t)) |x′ (t)| ψ(t) dt + ı η us − = x 7→ ∂ν Γ Γ y 0 0 (3.8) Z 2π Z 2π s ∂u ∂E ∂ ′ x 7→ E(x, x(t)) ψ(t) dt + , (x, x(t)) |x (t)| ψ(t) dt + ı η = ∂ν Γ− ∂ν ∂ν y Γ 0 0 where Γ+ and Γ− denote respectively exterior and interior limits onto Γ. In order to establish (3.8) in terms of boundary integral operators, we now introduce the parameterized versions of the single and double layer acoustic potentials Z 2π Z 2π K(s, t) g(t) dt , V (s, t) g(t) dt and Dg(s) := Sg(s) := 0
0
5
where V (s, t) := and K(s, t) := − (1)
ı (1) H (k |x(s) − x(t)|) 4 0
k ı (1) x′ (t) (x1 (t) − x1 (s)) − x′1 (t) (x2 (t) − x2 (s)) H1 (k |x(t) − x(s)|) 2 , 4 |x(t) − x(s)|
H1 being the Hankel function of the first kind and order one. We also introduce the hypersingular operator H which is related to the single layer operator, via tangential derivatives, as follows (see [8, page 295]): h H ψ, ϕ i = h ϕ′ , Sψ ′ i − k2 h Se ψ, ϕ i ∀ ψ, ϕ ∈ H 1/2 , (3.9) where Se is the integral operator whose kernel is given by Ve (t, s) := x′ (t) · x′ (s) V (t, s). We now deduce from the first equation of (2.3) that the tangential components of b−1 curl u must be continuous through Γ. Thus, recalling that b(x) = 1 for all x ∈ Ω0 , noting that the unit tangential x′ vector is given by |x ′ | , and using (3.2), we easily find that ′
γτ (p) = |x |
∂u ◦x ∂ν
i
′
= ρ + |x |
∂us ◦x ∂ν
,
(3.10)
∂ui ◦ x . Hence, applying the jump properties of the layer potentials and paramwhere ρ := |x | ∂ν eterizing the resulting integrals on Γ, the equations in (3.8) become i
′
1 γ(u) − γ(ui ) = ( I + D) ψ + ı η S ψ , 2
(3.11)
1 γτ (p) − ρi = − H ψ − ı η ( I − D t ) ψ , (3.12) 2 where I denotes, from now on, a generic identity operator and D t is the transposed of D. In this way, combining (3.6), (3.11), and (3.12), we obtain the following mixed-FEM/BEM formulation of (2.3): Find (p, ψ) ∈ H(curl; Ω) × H 1/2 such that 1 b p, q − k−2 a−1 curl p, curl q − h γτ (q), ( I + D)ψ + ı η Sψ i = h γτ (q), γ(ui ) i, 2 1 1 1 1 1 i h γτ (p) + Hψ + ı η ( I − D t )ψ, ϕ i = hρ , ϕi 2 2 2 2 2
(3.13)
for all (q, ϕ) ∈ H(curl; Ω) × H 1/2 . In fact, note that replacing γ(u) from (3.11) into (3.6) leads to the first equation of (3.13), whereas the second one arises after multiplying (3.12) by 12 and then testing it with ϕ ∈ H 1/2 . Next, defining the space M := H(curl; Ω) × H 1/2 provided with its natural inner product, we find that problem (3.13) can be rewritten, equivalently, as: Find (p, ψ) ∈ M such that A((p, ψ), (q, ϕ)) = F((q, ϕ))
∀ (q, ϕ) ∈ M ,
(3.14)
where A : M × M → C is the bilinear form defined by 1 1 A((p, ψ), (q, ϕ)) := b p, q − k−2 a−1 curl p, curl q − h γτ (q), ψ i + h γτ (p), ϕ i 2 2 1 1 1 + h Hψ, ϕ i + ı η h ψ, ϕ i − h γτ (q), Dψ i − ı η h γτ (q), Sψ i − ı η h ψ, Dϕ i , 2 4 2
6
(3.15)
and F : M → C is the linear functional given by F((q, ϕ)) := h γτ (q), γ(ui ) i +
1 i h ρ ,ϕ i. 2
In the following lemma, we collect some useful properties of the operators S, D and H, which will be used below in Section 4 to prove the well-posedness of (3.14). Lemma 3.1. The mappings S : H p → H p+1 , D : H p → H p+2 and S − Λ0 : H p → H p+3 are bounded for any p ∈ R, where Z 2π 4 1 t − s ϕ(s) ds . Λ0 (ϕ)(t) := − sin2 log 2π 0 e 2 Moreover, the operator H0 : H 1/2 → H −1/2
ϕ 7→ H0 ϕ := − (Λ0 ϕ′ )′ +
is H 1/2 −elliptic:
h H0 ϕ, ϕ i ≥
√
2 π kϕk21/2
1 2π
Z
2π
ϕ dt
0
∀ ϕ ∈ H 1/2 .
(3.16)
Proof. See Lemma 4.1 of [10].
4
Existence and uniqueness of solution
In this section we follow [3] and employ a suitable decomposition of H(curl; Ω) to show that (3.14) becomes a compact perturbation of a well-posed problem. For this purpose, we introduce the linear operator R : H(curl; Ω) → H(curl; Ω) q 7→ R(q) := curl v ,
where v ∈ H 1 (Ω) is the unique (up to an additive constant) solution of the auxiliary problem R curl q ∂v in Ω , = 0 on Γ , −∆ v = curl q − Ω |Ω| ∂ν with |Ω| denoting the mesure of Ω. It follows that R curl q curl R(q) = curl q − Ω in Ω |Ω|
and
γτ (R(q)) = 0 .
(4.1)
(4.2)
Then, using the continuous dependence result for (4.1) we find that kR(q)kH(curl;Ω) ≤ C kcurl qk0,Ω
∀q ∈ H(curl; Ω),
which shows that R is bounded. Moreover, it is easy to see that R is actually a projector and hence there holds H(curl; Ω) = R(H(curl; Ω)) ⊕ (I − R)(H(curl; Ω)) . (4.3) 7
On the other hand, owing to the regularity result of the Poisson equation with Neumann boundary conditions, we know that v ∈ H 2 (Ω) and kvk2,Ω ≤ CR kcurl qk0,Ω , which yields R(q) ∈ H1 (Ω) and kR(q)k1,Ω ≤ CR kcurl qk0,Ω
∀ q ∈ H(curl; Ω) .
(4.4)
We now employ the stable decomposition (4.3) to reformulate (3.14) in a more suitable form. Indeed, writing p = R(p) + (I − R)(p) and q = R(q) + (I − R)(q) in (3.15), and using that γτ (R(q)) = 0 (cf. (4.2)), we find that A can be decomposed as A((p, ψ), (q, ϕ)) = A0 ((p, ψ), (q, ϕ)) + K((p, ψ), (q, ϕ)) , where A0 : M × M → C and K : M × M → C are the bounded bilinear forms given by A0 ((p, ψ), (q, ϕ)) := − b R(p), R(q) − k−2 a−1 curl R(p), curl R(q) + b (I − R)(p), (I − R)(q) + k−2 a−1 curl (I − R)(p), curl (I − R)(q) −
(4.5)
1 1 1 1 h γτ (I − R)(q), ψ i + h γτ (I − R)(p), ϕ i + h H0 ψ, ϕ i + ı η h ψ, ϕ i , 2 2 2 4
and K((p, ψ), (q, ϕ)) := 2 (b R(p), R(q)) + (b (I − R)(p), R(q)) + (b R(p), (I − R)(q)) − 2 k−2 (a−1 curl (I − R)(p), curl (I − R)(q)) − k−2 (a−1 curl R(p), curl (I − R)(q)) 1 (H − H0 )ψ, ϕ 2 1 − h γτ (q), Dψ i − ı η h γτ (q), Sψ i − ı η h ψ, Dϕ i . 2
− k−2 (a−1 curl (I − R)(p), curl R(q)) +
(4.6)
Next, we let M′ be the dual of M pivotal to L2 (Ω) × L2 (0, 2π), whence M ⊂ L2 (Ω) × L2 (0, 2π) ⊂ M′ with dense inclusions, and denote by [·, ·] the duality bracket between M′ and M. Then, we introduce the bounded linear operators A : M → M′ , A0 : M → M′ , and K : M → M′ , induced by the bilinear forms A, A0 , and K, respectively, that is [A(u), v] := A(u, v) ,
[A0 (u), v] := A0 (u, v) ,
and [K(u), v] := K(u, v)
∀ u, v ∈ M ,
so that A = A0 + K and (3.14) can be reformulated as: Find (p, ψ) ∈ M such that (A0 + K)((p, ψ)) = F .
(4.7)
In the following two lemmas we prove that A0 is an isomorphism and that K is compact, respectively. Lemma 4.1. There exists a constant C > 0 such that kA0 ((p, ψ))kM′ :=
sup (q,ϕ)∈M
(q,ϕ)6=0
|A0 ((p, ψ), (q, ϕ))| ≥ C k(p, ψ)kM k(q, ϕ)kM
∀ (p, ψ) ∈ M .
(4.8)
In addition, there holds sup |A0 ((p, ψ), (q, ϕ))| > 0
∀ (q, ϕ) ∈ M ,
(p,ψ) ∈ M
8
(q, ϕ) 6= 0 .
(4.9)
Proof. Let us introduce the linear and bounded operator Ξ:M → M
(4.10)
(p, ψ) → Ξ(p, ψ) := ((I − 2R)(p), ψ) .
It is easy to see that Ξ is bijective and Ξ−1 = Ξ. Moreover, since R is a projector, we find that R (I − 2R)(p) = − R(p)
and (I − R) (I − 2R)(p) = (I − R)(p)
∀ p ∈ H(curl; Ω) , (4.11)
which yields A0 ((p, ψ), Ξ(p, ψ)) = b R(p), R(p) + k−2 a−1 curl R(p), curl R(p) + b (I − R)(p), (I − R)(p) + k−2 a−1 curl (I − R)(p), curl (I − R)(p) −
(4.12)
1 1 1 1 h γτ (I − R)(p), ψ i + h γτ (I − R)(p), ψ i + h H0 ψ, ψ i + ı η h ψ, ψ i . 2 2 2 4
It follows that o n aR R(p) Re A0 ((p, ψ), Ξ(p, ψ)) = b R(p), R(p) + k−2 curl R(p), curl |a|2 aR 1 + b (I − R)(p), (I − R)(p) + k−2 (I − R)(p) + h H0 ψ, ψ i , curl (I − R)(p), curl |a|2 2 and hence, according to (2.5), (2.6), and (3.16), we deduce that o n o n Re A0 ((p, ψ), Ξ(p, ψ)) ≥ β kR(p)k2H(curl;Ω) + k(I − R)(p)k2H(curl;Ω) + kψk21/2 ,
(4.13)
with β > 0 depending on k, ǫ0 , ǫ1 , σ1 , and ω. Then, using the triangle inequality we easily obtain kpk2H(curl;Ω) ≤ 2 kR(p)k2H(curl;Ω) + 2 k(I − R)(p)k2H(curl;Ω) , which, together with (4.13), gives n o Re A0 ((p, ψ), Ξ(p, ψ)) ≥ C k(p, ψ)k2M
∀(p, ψ) ∈ M .
Moreover, (4.14) and the boundedness of Ξ imply that n o Re A0 ((p, ψ), Ξ(p, ψ)) ≥ C k(p, ψ)kM kΞ(p, ψ)kM
(4.14)
∀(p, ψ) ∈ M ,
from which the inf-sup condition (4.8) follows straightforwardly. On the other hand, applying the estimate (4.14) with (p, ψ) = Ξ(q, ϕ) for a given (q, ϕ) ∈ M, (q, ϕ) 6= 0, and recalling that Ξ−1 = Ξ, we get o n (4.15) Re A0 (Ξ(q, ϕ), (q, ϕ)) ≥ C kΞ(q, ϕ)k2M > 0 , which provides (4.9) and completes the proof.
Lemma 4.2. The operator K is compact. 9
Proof. Let us decompose K = K1 + K2 + K3 , with Kj : M → M′ , j ∈ {1, 2, 3}, given by [K1 (p, ψ), (q, ϕ)] := 2 (b R(p), R(q)) + (b (I − R)(p), R(q)) + (b R(p), (I − R)(q)) , [K2 (p, ψ), (q, ϕ)] := − 2 k−2 (a−1 curl (I − R)(p), curl (I − R)(q)) − k−2 (a−1 curl R(p), curl (I − R)(q)) − k−2 (a−1 curl (I − R)(p), curl R(q)) , and [K3 (p, ψ), (q, ϕ)] :=
1 1 (H − H0 )ψ, ϕ − h γτ (q), Dψ i − ı η h γτ (q), Sψ i − ı η h ψ, Dϕ i . 2 2
The boundedness of R : H(curl; Ω) → H1 (Ω) (cf. (4.4)) and the compact embedding of H1 (Ω) into L2 (Ω), imply that K1 is compact. In addition, it is clear from (4.2) that curl (I − R) : H(curl; Ω) → L2 (Ω)
q → curl (I − R)(q) =
R
Ω
curl q |Ω|
is a finite rank operator, whence K2 also becomes compact. Finally, the compactness of K3 follows directly from Lemma 3.1. We are now in a position to show the unique solvability of (3.13) (equivalently (3.14) or (4.7)). Theorem 4.3. There exists a unique (p, ψ) ∈ M solution of problem (3.13). Proof. It follows from Lemmas 4.1 and 4.2 that A is a Fredholm operator of index zero and then, the present proof simply reduces to show uniqueness for (3.13). In fact, let (p0 , ψ0 ) be a solution of (3.13) with ui = 0 and define the function u0 (x) if x ∈ Ω, u e(x) := z(x) if x ∈ Ω0 ,
where
u0 (x) := k−2 a−1 curl p0 (x) and z(x) :=
Z
0
2π
∂E (x, x(t)) |x′ (t)| ψ0 (t) dt + ı η ∂ν y
Z
2π
E(x, x(t)) ψ0 (t) dt .
0
From the first equation in (3.13) we easily see that p0 = b−1 curl u0 in Ω, whence curl (b−1 curl u0 ) − k2 a u0 = 0
in Ω
1 and γ(u0 ) = ( I + D) ψ0 + ı η S ψ0 . 2
(4.16)
Further, the second equation in (3.13) gives directly 1 γτ (p0 ) = − H ψ0 − ı η ( I − D t ) ψ0 . 2
(4.17)
On the other hand, it is clear that z solves the Helmholtz equation ∆z + k2 z = 0 in Ω0 10
(4.18)
and satisfies the Sommerfeld radiation condition. In addition, the jump properties of the double layer potential and the normal derivative of the single layer potential through Γ imply (cf. [13]): γ(z) ∂z |x (t)| ◦x ∂ν ′
= =
1 ( I + D) ψ0 + ı η S ψ0 , 2 1 − H ψ0 − ı η ( I − D t ) ψ0 . 2
(4.19)
As a consequence of (4.16), (4.17), (4.18), and (4.19), we deduce that u e solves (2.3) with datum ui = 0 and satisfies the Sommerfeld condition (2.4). In this way, Theorem 2.1 ensures that u e vanishes identically, which yields u0 = 0, and therefore p0 = 0 and 1 ( I + D) ψ0 + ı η S ψ0 = 0 . 2
1 Finally, since ( I + D) + ı η S is one to one (see [13, Theorem 3.3.9]), we obtain that ψ0 = 0. 2
5
The discrete problem
We first construct a triangulation covering exactly the domain Ω. To this end, given N ∈ N, we let 0 = t0 < t1 < · · · < tN = 2π be a uniform partition of [0, 2π] with tj+1 −tj = h = 2π N for j ∈ {0, 1, ..., N −1}, and denote by Ωh the polygonal domain whose vertices are {x(t1 ), x(t2 ), ..., x(tN )}. Then we let Th0 be a regular triangulation of Ωh by triangles T of diameter hT such that h := sup{hT : T ∈ Th0 }. Further, we assume that any vertex of a triangle touching ∂Ωh belongs to {x(t1 ), x(t2 ), ..., x(tN )}. Thus, we replace each triangle T ∈ Th0 with one side along ∂Ωh , by the corresponding curved triangle with one side along Γ. In this way, we obtain from Th0 a triangulation Th of Ω made up of straight and curved triangles. Let T be a curved triangle of Th . We denote its vertices by aT1 , aT2 and aT3 , numbered in such a way that aT2 and aT3 are the endpoints of the curved side of T . Let i ∈ {0, 1, ..., N − 1} be such that x(ti ) = aT2 and x(ti+1 ) = aT3 . Then, ϕ(t) := x (ti + t h), t ∈ [0, 1], is a parameterization of the curved side of T . Next, we let Tb be the reference triangle with vertices b a1 := (0, 0), b a2 := (1, 0) b3 := (0, 1), and consider the affine map GT defined by GT (b and a ai ) = aTi for i ∈ {1, 2, 3}. Also, we introduce the function ΘT : Tb → R2 given by ΘT (b x) :=
x ˆ1 ϕ(ˆ x2 ) − (1 − x ˆ2 )aT2 − x ˆ2 aT3 1−x ˆ2
b ∈ Tb , ∀x
where the limiting value has to be taken as x ˆ2 goes to 1, and define FT := GT + ΘT : Tb → R2 . Now, if T is a straight (interior) triangle, we just take the curving perturbation ΘT ≡ 0 and thus FT is the usual affine map GT from the reference triangle. Lemma 5.1. Let T ∈ Th be a curved triangle. Then, there exists h0 > 0 such that for all h ∈ (0, h0 ) FT is a diffeomorphism of class C ∞ that maps one-to-one Tˆ onto the curved triangle T in such a way ai ) = aTi for all i ∈ {1, 2, 3}. In addition, the Jacobian JT of this mapping does not vanish in that FT (b a neighborhood of Tˆ, and there exist positive constants Ci , i ∈ {1, 2, 3}, independent of T , such that C1 h2T ≤ |JT (b x)| ≤ C2 h2T and max |∂ α ΘT (b x)| ≤ C3 hk+1 T b∈Tb x
∀ α ∈ N2 , 11
b ∈ Tb , ∀x |α| = k ,
(5.1) ∀k ≥ 1.
(5.2)
Proof. See [12] or Theorem 22.4 of [14]. Further estimates concerning curved triangles are collected in the following lemma from [1]. Lemma 5.2. Assume that T is a curved triangle and let BT := D GT . A function v belongs to H 1 (T ) if and only if the function vb = v ◦ FT belongs to H 1 (Tb). In this case, there exist C3 , C4 > 0, independent of T , such that o n 2 (5.3) k kB k kb v k + |b v | |v|1,T ≤ C3 |det BT |1/2 kB−1 T b b T 0,T 1,T , and
|b v |1,Tb ≤ C4 |det BT |−1/2 kBT k kvk1,T .
(5.4)
Proof. See Lemma 2.3 in [1]. We now introduce the first order rotated Raviart-Thomas space W(Tb) ⊆ H(rot; Tb) (cf. [2, 11]), which is given by: a1 x b2 b W(T ) := +b : a1 , a2 , b ∈ C , a2 −b x1
and define
W(T ) :=
n
b ◦ F−1 q = (D FT )−t q T ,
q:
b ∈ W(Tb) q
o
,
where (D FT )−t is the inverse of the transpose of D FT . Then, we seek an approximation of the unknown p in the finite dimensional space Wh := { q ∈ H(curl; Ω) :
q|T ∈ W(T )
∀ T ∈ Th } .
Next, let us denote by Eh : H1 (Ω) → Wh the equilibrium interpolation operator (see [2, 11]). Then, applying (5.1), (5.3), (5.4) and the Bramble-Hilbert Lemma, we deduce that k(I − Eh )(q)kH(curl;Ω) ≤ C h
(
kqk1,Ω +
X
T ∈Th
kcurl qk21,T
!1/2 )
(5.5)
for all q ∈ H1 (Ω) such that (curl q)|T ∈ H 1 (T ), ∀ T ∈ Th . In particular, using Lemmas 5.1 and 5.2, and the estimate (5.5), we can prove the following result. Lemma 5.3. There exists a constant C > 0, independent of h, such that o n k(I − Eh )(q)kH(curl;Ω) ≤ C h kqk1,Ω + kcurl qk0,Ω for all q ∈ H1 (Ω) such that curl q −
R
Ω
curl q |Ω| R
Proof. Let q ∈ H1 (Ω) such that curl q − (5.5), that k(I − Eh )(q)kH(curl;Ω) ≤ C h
(
Ω
kqk1,Ω +
(5.6)
∈ curl (Wh ). curl q |Ω|
= curl r, with r ∈ Wh . It follows, according to
X
T ∈Th
12
kcurl
qk20,T
+ |curl
r|21,T
!1/2 )
,
(5.7)
and then, the rest of the proof reduces to show that |curl r|21,T can be estimated in terms of kcurl qk20,T . Indeed, we first observe that if T is a straight triangle then curl (r|T ) is constant, whence |curl r|1,T vanishes. We now consider a curved triangle T . In this case, (5.3) implies that o n 2 [ [ k kB k k curl rk + | curl r| |curl r|1,T ≤ C3 |det BT |1/2 kB−1 T T 0,Tb 1,Tb , \ and hence, the identity curl (r|T ) =
1 JT
curl b r, together with the fact that curl b r is constant, give
o n −1 2 \rk0,T (b x )| max |J (b x )| kcurl k kB k + max |∇J |curl r|1,T ≤ C |det BT |1/2 kB−1 T T T T ≤ C
|det BT
|1/2
min |JT (b x)|1/2
kB−1 T k
b∈Tb x
n
b∈Tb x
b∈Tb x
x)| max |JT (b x)| kBT k2 + max |∇JT−1 (b b∈Tb x
b∈Tb x
o
kcurl rk0,T .
Next, we deduce from (5.1) and (5.2) that x)| = max max |∇JT−1 (b b∈Tb x
b ∈Tb x
|∇ JT (b x)| ≤ C h−1 T , 2 JT (b x)
2 which, together with the usual estimates for kB−1 T k and kBT k , yields R
Ω curl q
|curl r|1,T ≤ C kcurl rk0,T = C curl q − ≤ C kcurl qk0,T , |Ω| 0,T
(5.8)
where C > 0 is independent of T . Finally, (5.7) and (5.8) provide (5.6) and finish the proof.
It remains to define the finite element subspace for the boundary unknown ψ. In fact, given n ∈ N, we let Tn be the subspace of H 1/2 spanned by the 2π-periodic trigonometric functions of degree not greater than n, that is n n X X bj sin(jt) ∀ t ∈ [0, 2π] , aj , bj ∈ C . aj cos(jt) + Tn := ϕ ∈ H 1/2 : ϕ(t) = j=0
j=1
We note here that the orthogonal projection operator Pn : L2 (0, 2π) → Tn satisfies the following error estimate (cf. Theorem 8.2.1 of [13]) r−s 2 kϕkr ∀ ϕ ∈ Hr , ∀ r ≥ s . (5.9) kϕ − Pn (ϕ)ks ≤ n Now, given h > 0 and n ∈ N, we let Mh,n = Wh × Tn and introduce the following discrete version of (3.14): Find (ph , ψn ) ∈ Mh,n such that A((ph , ψn ), (q, ϕ)) = F((q, ϕ))
∀ (q, ϕ) ∈ Mh,n .
(5.10)
We now aim to prove the well-posedness of (5.10). For this purpose, as established by a classical result on projection methods for Fredholm operators of index zero (see, e.g. [7, Theorem 13.7]), we need to show that the Galerkin scheme associated to the isomorphism A0 is well-posed. Precisely, the following lemma establishes the corresponding discrete inf-sup condition for the bilinear form A0 on the finite element subspace Mh,n . 13
Lemma 5.4. There exist constants C, h0 > 0, independent of h and n, such that for each h ≤ h0 there holds |A0 ((p, ψ), (q, ϕ))| sup ≥ C k(p, ψ)kM ∀ (p, ψ) ∈ Mh,n . (5.11) k(q, ϕ)kM (q,ϕ)∈Mh,n (q,ϕ)6=0
Proof. Let us introduce the linear operator Ξh,n : Mh,n → Mh,n
(5.12)
(p, ψ) 7→ Ξh,n (p, ψ) := ((I − 2 Eh R)(p), ψ) ,
which is the discrete version of Ξ (cf. (4.10)). It follows, applying (5.6), (4.4), and (4.2), that kΞh,n (p, ψ) − Ξ(p, ψ)kM = 2 k(I − Eh ) (R(p))kH(curl;Ω) o n ≤ C h kcurl pk0,Ω , ≤ C h kR(p)k1,Ω + kcurl (R(p))k0,Ω
(5.13)
which, combined with the boundedness of Ξ, implies that Ξh is also bounded, with a constant independent of h and n. Hence, using (4.14), the boundedness of A0 , and (5.13), we find that for each (p, ψ) ∈ Mh,n there holds n o n o n o Re A0 ((p, ψ), Ξh,n (p, ψ)) = Re A0 ((p, ψ), Ξ(p, ψ)) + Re A0 ((p, ψ), (Ξh,n − Ξ)(p, ψ)) ≥ C k(p, ψ)k2M − C0 k(p, ψ)kM kΞh,n (p, ψ) − Ξ(p, ψ)kM ≥ (C − C0 h) k(p, ψ)k2M ,
which yields, for each h ≤ h0 , with h0 sufficiently small, n o Re A0 ((p, ψ), Ξh,n (p, ψ)) ≥ C k(p, ψ)k2M ≥ C k(p, ψ)kM kΞh,n (p, ψ)kM
∀ (p, ψ) ∈ Mh,n .
In this way, the discrete inf-sup condition (5.11) follows straightforwardly from the above estimate. The well-posedness and convergence of the discrete scheme (5.10) can now be established. ˜ 0 ∈ ]0, h0 ] and Theorem 5.5. Let h0 be the constant provided by Lemma 5.4. Then, there exist h ˜ n0 ∈ N such that for each h ≤ h0 and n ≥ n0 , the Galerkin scheme (5.10) has a unique solution (ph , ψn ) ∈ Mh,n . In addition, there exist c1 , c2 > 0, independent of h and n, such that k(ph , ψn )kM ≤ c1
sup (qh ,ϕn )∈Mh,n (qh ,ϕn )6=0
|F((qh , ϕn ))| k(qh , ϕn )kM
and k(p, ψ) − (ph , ψn )kM ≤ c2
inf
(qh ,ϕn )∈Mh,n
k(p, ψ) − (qh , ϕn )kM ,
where (p, ψ) is the unique solution of (3.14). Furthermore, if p ∈ H1 (Ω) and curl p ∈ H 1 (Ω), then there holds σ 2 kψkσ+1/2 k(p, ψ) − (ph , ψn )kM ≤ c3 h kpk1,Ω + h kcurl pk1,Ω + ∀ σ > 0 , (5.14) n where c3 > 0 is also independent of h and n. 14
Proof. We know from Lemmas 4.1 and 4.2, and Theorem 4.3 that the operator A is a compact perturbation of the isomorphism A0 , and that the corresponding problem (3.14) is well posed. In addition, it follows from Lemma 5.4 that the Galerkin scheme associated to A0 is stable and convergent for each h ≤ h0 . Also, it is easy to show that the approximation properties (5.5) and (5.9) and the density of smooth functions in M yield lim inf k(p, ψ) − (qh , ϕn )kM = 0 ∀ (p, ψ) ∈ M . (5.15) h→0
n→∞
(qh ,ϕn )∈Mh,n
Consequently, the proof of the first part is a direct application of [7, Theorem 13.7], whereas the rate of convergence (5.14) follows from the approximation properties (5.5) and (5.9), and the fact that the boundary unknown ψ belongs to H σ for each σ ≥ 0. Indeed, since u solves the homogeneous Helmholtz equation in Ω′c , it becomes analytic in a sufficiently small neighborhood of Γ and hence γ(us ) := γ(u − ui ) ∈ H σ for each σ ≥ 0 (assuming ui sufficiently smooth). Therefore, since 1 1 σ σ −1 (γ(us )) ∈ H σ 2 I + D + ı η S : H → H is an isomorphism, we conclude that ψ = ( 2 I + D + ı η S) for each σ ≥ 0. We remark here that if the solution u of the original problem (2.3) belongs to H 2 (Ω), then it follows easily from (3.3) and (3.4) that p ∈ H1 (Ω), curl p ∈ H 1 (Ω), and kpk1,Ω + kcurl pk1,Ω ≤ C kuk2,Ω , whence the a priori error estimate (5.14) becomes σ 2 k(p, ψ) − (ph , ψn )kM ≤ C h kuk2,Ω + kψkσ+1/2 ∀σ > 0. n
6
Numerical results
In this section we present an example illustrating the performance of our mixed finite element scheme (5.10). We consider the domain Ωc as the annular region determined by the circle centered at the origin with radius 0.7 and the ellipse of minor and major semiaxis 1.3 and 1.5, respectively. The C ∞ closed curve Γ is taken then as this ellipse. In addition, the coefficients a and b are given in Ωc by (1 + |x|2 ) + ı|x|4 if 0.7 < |x| < 1 , 1 + (1 − |x|4 )2 if 0.7 < |x| < 1 , a(x) = and b(x) = 2+ı 1 otherwise , otherwise . (1)
We compute the right hand side in such a way that the Hankel function u(x) := H0 (k|x|) is the exact solution of the problem. In what follows, k2 = 10.50, which is an approximation of the sixth eigenvalue of the Laplacian in Ω with a Neumann boundary condition on Γ. Thus, our method is tested in the case where the numerical scheme given in [6] may exhibit an unstable behavior because of the lack of uniqueness. The numerical results shown below were obtained in a Pentium Xeon computer with dual processors, using a MATLAB code. In Table 1, we take h = 2 π/192 and decrease the spectral parameter n until we obtain the smallest value that preserves the order of accuracy. We can see that the number of degrees of freedom on the boundary is drastically reduced. Hence, for relatively low wave numbers k, it is possible to eliminate the boundary variable at matricial level (by a static condensation procedure) and solve the reduced linear system for the finite element variable. Then, in Table 2 we fix n = 16 and change h in order to 15
check out the asymptotic convergence predicted by the theory for p − ph . We also display there the experimental rate of convergence defined by r(p) :=
logkp − ph kH(curl;Ω) − logkp − ph˜ kH(curl;Ω) , ˜ log(h) − log(h)
˜ denote two consecutive mesh sizes. The computed values of r(p) show that the error where h and h kp − ph kH(curl;Ω) decreases linearly with respect to the mesh parameter h. Moreover, if we denote by Qh the L2 (Ω)-orthogonal projection onto the space of piecewise constant function associated to the triangulation Th , we can use the postprocess computation uh = k−2 Qh (a−1 curl ph )
in
Ω.
to obtain an approximation uh of the scalar variable u. Note from (3.4) that u = k−2 a−1 curl p in Ω. The behavior of the error ku − uh k0,Ω and the corresponding experimental rate of convergence r(u) :=
logku − uh k0,Ω − logku − uh˜ k0,Ω , ˜ log(h) − log(h)
are also depicted in Table 2. 2n 128 64 32 16 8
kp − ph kH(curl;Ω) 2.64 × 10−1 2.64 × 10−1 2.64 × 10−1 2.68 × 10−1 3.86 × 10+0
Table 1: Convergence history for different values of n when k2 = 10.50 and h = 2 π/192.
h 2π/48 2π/64 2π/96 2π/128 2π/192 2π/256
kp − ph kH(curl;Ω) 1.00 × 10+0 8.21 × 10−1 5.31 × 10−1 4.02 × 10−1 2.64 × 10−1 1.99 × 10−1
r(p) — 0.70 1.07 0.96 1.03 0.98
ku − uh k0,Ω 6.04 × 10−3 3.67 × 10−3 1.92 × 10−3 9.91 × 10−4 4.21 × 10−4 2.55 × 10−4
r(u) — 1.72 1.59 2.31 2.10 1.74
Table 2: Convergence history for different values of h when k2 = 10.50 and n = 16.
References [1] C. Bernardi, Optimal finite-element interpolation on curved domains. SIAM Journal on Numerical Analysis, vol. 26, pp. 1212-1240, (1989). [2] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods. Springer, Berlin Heidelberg New York, 1991. 16
[3] A. Buffa, Remarks on the discretization of some noncoercive operator with applications to heterogeneous Maxwell equations. SIAM Journal on Numerical Analysis, vol. 43, pp. 1-18, (2005). [4] D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Second Edition. Springer-Verlag, Berlin, 1998. [5] J. Coyle and P. Monk, Scattering of time-harmonic electromagnetic waves by anisotropic inhomogeneous scatterers or impenetrable obstacles. SIAM Journal on Numerical Analysis, vol. 37, pp. 1590-1617, (2000). [6] G.N. Gatica and S. Meddahi, On the coupling of MIXED-FEM and BEM for an exterior Helmholtz problem in the plane. Numerische Mathematik, vol. 100, pp. 663-695, (2005). [7] R. Kress, Linear Integral Equations, Second Edition. Springer-Verlag, New York, 1999. [8] W. McLean, Strongly Eliptic Systems and Boundary Integral Equations. Cambridge University Press, 2000. ´rquez and V. Selgas, Computing acoustic waves in an inhomogeneous [9] S. Meddahi, A. Ma medium of the plane by a coupling of spectral and finite elements. SIAM Journal on Numerical Analysis, vol. 41, pp. 1729-1750, (2003). [10] S. Meddahi and F. J. Sayas, Analysis of a new BEM-FEM coupling for two-dimensional fluidsolid interaction. Numerical Methods for Partial Differential Equations, vol. 21, 6, pp. 1017-1042, (2005). [11] P. A. Raviart and J. M. Thomas, A mixed finite element method for 2nd order elliptic problems. In: Mathematical Aspects of the Finite Element Method, Lecture Notes in Mathematics, vol. 606, Springer-Verlag, Berlin, 1977. [12] R. Scott, Interpolated boundary conditions in the finite element method, SIAM Journal on Numerical Analysis, vol. 12 , pp. 404-427, (1975). [13] J. Sarannen and G. Vainikko, Periodic Integral and Pseudodifferential Equations with Numerical Approximation. Springer Verlag, 2002. ˇ ´ıˇ [14] A. Zen sek, Nonlinear Elliptic and Evolution Problems and their Finite Element Approximations. Academic Press, London, 1990.
17
Abstract This paper deals with the scattering of time harmonic electromagnetic waves by an infinitely long cylinder containing a non-homogeneous conducting medium. More precisely, we study the transverse magnetic field that solves an interface problem holding between the cross section of the cylinder and the exterior two-dimensional free space. We apply a dual-mixed variational formulation in the obstacle coupled with a boundary integral equation method in the unbounded homogeneous space. A Fredholm alternative is utilized to prove that this continuous formulation is well posed. We define the corresponding discrete scheme by using the lowest order rotated Raviart-Thomas finite elements for the magnetic field and spectral elements for the boundary unknown. Then, we show that the resulting Galerkin scheme is uniquely solvable and convergent, and prove optimal error estimates. Finally, we illustrate our analysis with some results from computational experiments. Key words. mixed-FEM, BEM, coupling, Helmholtz equation. Mathematics subject classifications (1991). 65N30, 65N38, 65N15, 35Q60
1
Introduction
In this paper we consider the exterior Helmholtz equation that arises from a time harmonic electromagnetic scattering problem posed in the plane. The scatterer is a heterogeneous conducting body embedded in an infinite, homogeneous, isotropic, and non-conducting medium. Since the underlying partial differential equation is coercive, we observe that the treatment of the usual primal variational formulation is a rather straightforward application of the Fredholm alternative (see, e.g. [4, 9] and the references therein for more details). On the other hand, it is well known (see, e.g. [2, 11]) that a dual-mixed variational formulation is more adequate to obtain a good approximation of the vectorial unknown, namely the magnetic field in the conductor. A procedure in this direction is developed recently in [6], where the coupling of mixed finite element (mixed-FEM) and boundary element methods (BEM) is employed for the corresponding discrete scheme. In fact, the approach in [6] uses a dual-mixed variational formulation in a bounded domain delimited by a C ∞ closed curve, which is provided with the non-local boundary conditions ∗
This research was partially supported by CONICYT-Chile through the FONDAP Program in Applied Mathematics, by the Direcci´ on de Investigaci´ on of the Universidad de Concepci´ on through the Advanced Research Groups Program, and by the Ministerio de Educaci´ on y Ciencia of Spain, through the project No. MTM2004-05417. † GI2 MA, Departamento de Ingenier´ıa Matem´ atica, Facultad de Ciencias F´ısicas y Matem´ aticas, Universidad de Concepci´ on, Casilla 160-C, Concepci´ on, Chile, e-mail: [email protected] ‡ Departamento de Construcci´ on e Ingenier´ıa de Fabricaci´ on, Universidad de Oviedo, Oviedo, Espa˜ na, e-mail: [email protected] § Departamento de Matem´ aticas, Facultad de Ciencias, Universidad de Oviedo, Calvo Sotelo s/n, Oviedo, Espa˜ na, e-mail: [email protected]
1
arising from the application of the boundary integral equation method in the unbounded region. The corresponding analysis considers first the auxiliary problem obtained from the original formulation after substituting all the acoustic boundary integral operators by the corresponding potential boundary integral operators. Then, the properties of these operators and suitable duality arguments allow to prove a global inf-sup condition that implies the well posedness of the auxiliary problem. Afterwards, the Fredholm alternative is applied to the original problem to conclude that it is also well posed. Nevertheless, a drawback of the resulting coupled method is the fact that it is not well-posed if the square of the wave number is a Neumann eigenvalue of the Laplace operator in the interior domain. Thus, the corresponding numerical method may exhibit an unstable behavior in the neighborhood of these singular coefficients. Furthermore, the anonymous referee of this paper suggested us that eliminating the scalar variable, which constitutes the Lagrange multiplier in [6], would lead to a variational formulation whose structure agrees with the abstract framework provided by Buffa in [3] for a certain type of non-coercive problems. Consequently, the purpose of the present paper is to follow the above suggestion and develop a new coupling of mixed finite element and boundary element methods for the exterior Helmholtz equation described before. More precisely, we show that the powerful tool from [3] simplifies considerably the analysis presented in [6] for both the continuous and the discrete problem. In addition, the coupling of mixed-FEM and BEM is modified here in such a way that, on the contrary to [6], the discrete scheme now becomes free from spurious solutions. The rest of the paper is organized as follows. In Sections 2, 3 and 4 we introduce the model problem, derive our coupled formulation and prove its well-posedness. In Section 5 we introduce the Galerkin mixed-FEM/BEM discretization and provide the associated convergence result. Finally, Section 6 is devoted to the numerical results. √ In the sequel, we deal with complex valued functions and the symbol ı is used for −1. We let α and |α| be the conjugate and the modulus, respectively, of a complex number α ∈ C. We use boldface letters to denote vectors, vector-valued functions or vector-valued function spaces. The symbol |·| is 2 X qi qi for all q := (q1 , q2 ) ∈ C2 . used to denote the 2-norm for vectors, that is |q| = q · q := i=1
2 In addition, Z given a bounded domain Ω ⊆ R with a sufficiently smooth boundaryqΓ, we denote f g for all f, g ∈ L2 (Ω). Thus, the norm in L2 (Ω) is given by kvk0,Ω := f, g := v, v . Also, Ω
for any non-negative integer m, H m (Ω) is the usual Hilbertian Sobolev space endowed with the norm k·km,Ω . We denote the corresponding seminorm |·|m,Ω ∞ be On the other hand, we also need to introduce periodic Sobolev spaces. To this end, we let C2π the space of 2π-periodic infinitely differentiable complex valued functions of a single variable, and for ∞ define its Fourier coefficients by each g ∈ C2π 1 gb(k) := 2π
Z
2π
g(t) e−kıt dt
0
∀ k ∈ Z.
∞ with respect to the Then, given r ∈ R, the Hilbertian Sobolev space H r is the completion of C2π !1/2 X norm kgkr := . Moreover, the H 0 -bilinear form defined as h λ, η i := (1 + k2 )r |b g (k)|2
Z
2π
k∈Z
λ(t) η(t) dt can be extended to represent the duality between H −r and H r for all r > 0. We will
0
keep the same notation for this duality bracket. Throughout this paper, C denotes a generic constant independent of the discretization parameters, which usually takes different values at different places.
2
2
Statement of the problem
We consider an incident electromagnetic wave upon a body containing a non-homogeneous conducting medium and we are interested in determining the scattered electromagnetic wave. The obstacle is assumed to be an infinitely long cylinder parallel to the x3 -axis whose cross section is a bounded open domain Ωc . We denote by Ω′c its complement in R2 . Next, we let ǫ, µ and σ be the electric permittivity, the magnetic permeability, and the conductivity, respectively. They are assumed to be constant along each line parallel to the x3 -axis. These coefficients are bounded real valued scalar functions that satisfy almost everywhere in Ωc ǫ0 ≤ ǫ(x) ≤ ǫ1 ,
µ0 ≤ µ(x) ≤ µ1 ,
and
0 ≤ σ(x) ≤ σ1 ,
(2.1)
where ǫ0 , ǫ1 , µ0 , µ1 , and σ1 are positive constants, and ǫ0 and µ0 are the electric permittivity and magnetic permeability, respectively, of the free space. Moreover, since the medium is dielectric and homogeneous, outside the obstacle there holds ǫ(x) = ǫ0 ,
µ(x) = µ0 ,
and
in Ω′c .
σ(x) = 0
The incident electromagnetic field is supposed to exhibit a time-harmonic behavior with frequency ω and to be polarized perpendicular to the x3 -axis. Under these conditions, the electric and magnetic fields will also have a time harmonic behavior with frequency ω. Their amplitudes E and H, which are independent of x3 , are given by 0 h1 (x1 , x2 ) 1 and H = √1 h2 (x1 , x2 ) , 0 E= √ ǫ0 µ0 u(x1 , x2 ) 0
where the complex-valued functions u and h = (h1 , h2 )t are now the unknowns of the problem. Then, from the time dependent form of Maxwell’s equations, we know that u and h satisfy the equations curl h + ı k a u = 0 in R2 , (2.2) curl u − ı k b h = 0 in R2 ,
where
∂h2 ∂h1 − curl h := ∂x1 ∂x2
Here, k = ω a(x) =
√
and
curl u :=
ǫ0 µ0 is the wave number, b(x) =
∂u ∂u ,− ∂x2 ∂x1
t
.
µ(x) , and the index of refraction is given by µ0
σ(x) ǫ(x) +ı for all x ∈ R2 . Eliminating h from (2.2) yields ǫ0 ω ǫ0 curl (b−1 curl u) − k2 a u = 0
in R2 ,
u = ui + us in R2 ,
(2.3)
where ui and us denote the incident and scattered waves, respectively, and ui is given, for example, by the plane wave ui (x) = eı k x·d in the direction determined by the unit vector d. Furthermore, the behavior of the scattered electromagnetic field at infinity (known as the Silver-Muller condition) implies the following Sommerfeld radiation condition: s √ ∂u s r − ıku → 0 when r = |x| → ∞ . (2.4) ∂r 3
Note that our hypothesis on the physical coefficients of the problem imply that 1 ≤ b(x) ≤
µ1 µ0
for a.e. x ∈ R2
(2.5)
and a(x) := aR (x) + ı aI (x) is a bounded complex-valued function with 1 ≤ aR (x) ≤
ǫ1 ǫ0
and
0 ≤ aI (x) ≤
σ1 ω ǫ0
for a.e. x ∈ R2 .
(2.6)
Moreover, in the unbounded homogeneous region Ω′c we have b(x) = 1 and
a(x) = 1.
We recall that, in practice, aI = 0 means that no loss of energy occurs in the medium. Now, we follow [5] and make assumptions on the coefficients in order to guarantee uniqueness of solution for problem (2.3). For this purpose, we first introduce a sufficiently large C ∞ closed curve Γ such that its interior region, denoted from now on by Ω, is convex and contains the heterogeneous medium Ωc . The outward unit vector normal to Γ is represented by ν. Then, we assume that there exists a set {Ω0 , Ω1 , · · · , ΩN } of open and disjoint subdomains with Lipschitz-continuous boundaries such that Ω = Ω1 ∪ Ω2 ∪ · · · ∪ ΩN , Ω0 = R2 − Ω , and in each subdomain Ωi , i ∈ {0, 1, · · · , M }, the coefficients a and b satisfy one of the following two conditions: (C1) b and aR are functions in H 3 (Ωi ) and aI vanishes identically in Ωi ; or (C2) aI (x) ≥ a0 > 0 almost everywhere in Ωi . Theorem 2.1. Assume that ui = 0. Then problem (2.3) only admits the trivial solution. Proof. We refer to Theorem 3.1 of [5] for a proof of this result.
3
The coupling of mixed-FEM and BEM
Let x : R → R2 be a regular 2π-periodic parametric representation of the curve Γ, that is |x′ (s)| > 0 ∀s ∈ R and
x(s) = x(t)
⇔
t − s ∈ 2 π Z.
We define the parameterized trace on Γ as the linear continuous extension of γ : C ∞ (Ω) → L2 (0, 2π) u 7→ γ(u) := u|Γ ◦ x
(3.1)
to H 1 (Ω), which yields a bounded and onto linear application γ : H 1 (Ω) → H 1/2 (see Theorem 8.15 in [7]). Similarly, the parameterized tangential trace is defined as the linear continuous extension of γτ : [C ∞ (Ω)]2 → H −1/2 q 7→ γτ (q) := (q|Γ ◦ x) · x′
(3.2)
to H(curl; Ω), which yields a bounded and onto linear application γτ : H(curl; Ω) → H −1/2 . We recall here that H(curl; Ω) := q ∈ L2 (Ω) : curl q ∈ L2 (Ω) , 4
and its Hilbertian norm is given by kqk2H(curl;Ω) := kqk20,Ω + kcurl qk20,Ω . In order to deduce a mixed variational formulation of (2.3) we introduce the vectorial unknown p := b−1 curl u in Ω, so that b p := curl u in Ω . (3.3) Note that p is related to the magnetic field by p = ı k h. With this notation, the restriction of the first equation of (2.3) to Ω becomes curl p − k2 a u = 0
in
Ω.
(3.4)
Then, testing (3.3) with the vector-valued function q ∈ H(curl; Ω) and integrating by parts yields b p, q − u, curl q − h γτ (q), γ(u) i = 0. (3.5)
Now, we proceed differently from [6] and replace u = k−2 a−1 curl p from (3.4) into the second term of (3.5) to obtain the following variational formulation in the bounded domain Ω: Find p ∈ H(curl; Ω) such that b p, q − k−2 a−1 curl p, curl q − h γτ (q), γ(u) i = 0 ∀ q ∈ H(curl; Ω). (3.6)
Next, we employ a boundary integral equation method to complete (3.6) with the corresponding variational formulation in the unbounded exterior domain Ω0 , thus providing, in particular, a suitable expression for γ(u) (see (3.11) below). More precisely, given an arbitrary real number η > 0, we consider the representation formula Z 2π Z 2π ∂E ′ s E(x, x(t)) ψ(t) dt ∀ x ∈ Ω0 , (3.7) (x, x(t)) |x (t)| ψ(t) dt + ı η u (x) = ∂ν y 0 0 where E(x, y) :=
ı (1) H (k|x − y|) 4 0 (1)
is the radial outgoing fundamental solution of the Helmholtz equation, H0 is the Hankel function of the first kind and order zero, and ψ is an unknown density. Expression (3.7) is the parameterized version of the combined single-double layer potential sometimes referred to as a Brakhage-Werner potential (see [4, 13]). Note that the function |x′ (t)| appears only in one of the integrands, which aims to obtain parameterized integral equations without cumbersome weights (see [10] for more details). We remark that a function us given by (3.7) satisfies the Helmholtz equation ∆ us + k2 us = 0 in Ω0 and the sommerfeld radiation condition (2.4) for any 2 π-periodic density ψ : [0, 2π] → C. Thus, if we want to obtain a global solution for problem (2.3), we have just to impose the following transmission conditions on Γ Z 2π Z 2π ∂E E(x, x(t)) ψ(t) dt + , (x, x(t)) |x′ (t)| ψ(t) dt + ı η us − = x 7→ ∂ν Γ Γ y 0 0 (3.8) Z 2π Z 2π s ∂u ∂E ∂ ′ x 7→ E(x, x(t)) ψ(t) dt + , (x, x(t)) |x (t)| ψ(t) dt + ı η = ∂ν Γ− ∂ν ∂ν y Γ 0 0 where Γ+ and Γ− denote respectively exterior and interior limits onto Γ. In order to establish (3.8) in terms of boundary integral operators, we now introduce the parameterized versions of the single and double layer acoustic potentials Z 2π Z 2π K(s, t) g(t) dt , V (s, t) g(t) dt and Dg(s) := Sg(s) := 0
0
5
where V (s, t) := and K(s, t) := − (1)
ı (1) H (k |x(s) − x(t)|) 4 0
k ı (1) x′ (t) (x1 (t) − x1 (s)) − x′1 (t) (x2 (t) − x2 (s)) H1 (k |x(t) − x(s)|) 2 , 4 |x(t) − x(s)|
H1 being the Hankel function of the first kind and order one. We also introduce the hypersingular operator H which is related to the single layer operator, via tangential derivatives, as follows (see [8, page 295]): h H ψ, ϕ i = h ϕ′ , Sψ ′ i − k2 h Se ψ, ϕ i ∀ ψ, ϕ ∈ H 1/2 , (3.9) where Se is the integral operator whose kernel is given by Ve (t, s) := x′ (t) · x′ (s) V (t, s). We now deduce from the first equation of (2.3) that the tangential components of b−1 curl u must be continuous through Γ. Thus, recalling that b(x) = 1 for all x ∈ Ω0 , noting that the unit tangential x′ vector is given by |x ′ | , and using (3.2), we easily find that ′
γτ (p) = |x |
∂u ◦x ∂ν
i
′
= ρ + |x |
∂us ◦x ∂ν
,
(3.10)
∂ui ◦ x . Hence, applying the jump properties of the layer potentials and paramwhere ρ := |x | ∂ν eterizing the resulting integrals on Γ, the equations in (3.8) become i
′
1 γ(u) − γ(ui ) = ( I + D) ψ + ı η S ψ , 2
(3.11)
1 γτ (p) − ρi = − H ψ − ı η ( I − D t ) ψ , (3.12) 2 where I denotes, from now on, a generic identity operator and D t is the transposed of D. In this way, combining (3.6), (3.11), and (3.12), we obtain the following mixed-FEM/BEM formulation of (2.3): Find (p, ψ) ∈ H(curl; Ω) × H 1/2 such that 1 b p, q − k−2 a−1 curl p, curl q − h γτ (q), ( I + D)ψ + ı η Sψ i = h γτ (q), γ(ui ) i, 2 1 1 1 1 1 i h γτ (p) + Hψ + ı η ( I − D t )ψ, ϕ i = hρ , ϕi 2 2 2 2 2
(3.13)
for all (q, ϕ) ∈ H(curl; Ω) × H 1/2 . In fact, note that replacing γ(u) from (3.11) into (3.6) leads to the first equation of (3.13), whereas the second one arises after multiplying (3.12) by 12 and then testing it with ϕ ∈ H 1/2 . Next, defining the space M := H(curl; Ω) × H 1/2 provided with its natural inner product, we find that problem (3.13) can be rewritten, equivalently, as: Find (p, ψ) ∈ M such that A((p, ψ), (q, ϕ)) = F((q, ϕ))
∀ (q, ϕ) ∈ M ,
(3.14)
where A : M × M → C is the bilinear form defined by 1 1 A((p, ψ), (q, ϕ)) := b p, q − k−2 a−1 curl p, curl q − h γτ (q), ψ i + h γτ (p), ϕ i 2 2 1 1 1 + h Hψ, ϕ i + ı η h ψ, ϕ i − h γτ (q), Dψ i − ı η h γτ (q), Sψ i − ı η h ψ, Dϕ i , 2 4 2
6
(3.15)
and F : M → C is the linear functional given by F((q, ϕ)) := h γτ (q), γ(ui ) i +
1 i h ρ ,ϕ i. 2
In the following lemma, we collect some useful properties of the operators S, D and H, which will be used below in Section 4 to prove the well-posedness of (3.14). Lemma 3.1. The mappings S : H p → H p+1 , D : H p → H p+2 and S − Λ0 : H p → H p+3 are bounded for any p ∈ R, where Z 2π 4 1 t − s ϕ(s) ds . Λ0 (ϕ)(t) := − sin2 log 2π 0 e 2 Moreover, the operator H0 : H 1/2 → H −1/2
ϕ 7→ H0 ϕ := − (Λ0 ϕ′ )′ +
is H 1/2 −elliptic:
h H0 ϕ, ϕ i ≥
√
2 π kϕk21/2
1 2π
Z
2π
ϕ dt
0
∀ ϕ ∈ H 1/2 .
(3.16)
Proof. See Lemma 4.1 of [10].
4
Existence and uniqueness of solution
In this section we follow [3] and employ a suitable decomposition of H(curl; Ω) to show that (3.14) becomes a compact perturbation of a well-posed problem. For this purpose, we introduce the linear operator R : H(curl; Ω) → H(curl; Ω) q 7→ R(q) := curl v ,
where v ∈ H 1 (Ω) is the unique (up to an additive constant) solution of the auxiliary problem R curl q ∂v in Ω , = 0 on Γ , −∆ v = curl q − Ω |Ω| ∂ν with |Ω| denoting the mesure of Ω. It follows that R curl q curl R(q) = curl q − Ω in Ω |Ω|
and
γτ (R(q)) = 0 .
(4.1)
(4.2)
Then, using the continuous dependence result for (4.1) we find that kR(q)kH(curl;Ω) ≤ C kcurl qk0,Ω
∀q ∈ H(curl; Ω),
which shows that R is bounded. Moreover, it is easy to see that R is actually a projector and hence there holds H(curl; Ω) = R(H(curl; Ω)) ⊕ (I − R)(H(curl; Ω)) . (4.3) 7
On the other hand, owing to the regularity result of the Poisson equation with Neumann boundary conditions, we know that v ∈ H 2 (Ω) and kvk2,Ω ≤ CR kcurl qk0,Ω , which yields R(q) ∈ H1 (Ω) and kR(q)k1,Ω ≤ CR kcurl qk0,Ω
∀ q ∈ H(curl; Ω) .
(4.4)
We now employ the stable decomposition (4.3) to reformulate (3.14) in a more suitable form. Indeed, writing p = R(p) + (I − R)(p) and q = R(q) + (I − R)(q) in (3.15), and using that γτ (R(q)) = 0 (cf. (4.2)), we find that A can be decomposed as A((p, ψ), (q, ϕ)) = A0 ((p, ψ), (q, ϕ)) + K((p, ψ), (q, ϕ)) , where A0 : M × M → C and K : M × M → C are the bounded bilinear forms given by A0 ((p, ψ), (q, ϕ)) := − b R(p), R(q) − k−2 a−1 curl R(p), curl R(q) + b (I − R)(p), (I − R)(q) + k−2 a−1 curl (I − R)(p), curl (I − R)(q) −
(4.5)
1 1 1 1 h γτ (I − R)(q), ψ i + h γτ (I − R)(p), ϕ i + h H0 ψ, ϕ i + ı η h ψ, ϕ i , 2 2 2 4
and K((p, ψ), (q, ϕ)) := 2 (b R(p), R(q)) + (b (I − R)(p), R(q)) + (b R(p), (I − R)(q)) − 2 k−2 (a−1 curl (I − R)(p), curl (I − R)(q)) − k−2 (a−1 curl R(p), curl (I − R)(q)) 1 (H − H0 )ψ, ϕ 2 1 − h γτ (q), Dψ i − ı η h γτ (q), Sψ i − ı η h ψ, Dϕ i . 2
− k−2 (a−1 curl (I − R)(p), curl R(q)) +
(4.6)
Next, we let M′ be the dual of M pivotal to L2 (Ω) × L2 (0, 2π), whence M ⊂ L2 (Ω) × L2 (0, 2π) ⊂ M′ with dense inclusions, and denote by [·, ·] the duality bracket between M′ and M. Then, we introduce the bounded linear operators A : M → M′ , A0 : M → M′ , and K : M → M′ , induced by the bilinear forms A, A0 , and K, respectively, that is [A(u), v] := A(u, v) ,
[A0 (u), v] := A0 (u, v) ,
and [K(u), v] := K(u, v)
∀ u, v ∈ M ,
so that A = A0 + K and (3.14) can be reformulated as: Find (p, ψ) ∈ M such that (A0 + K)((p, ψ)) = F .
(4.7)
In the following two lemmas we prove that A0 is an isomorphism and that K is compact, respectively. Lemma 4.1. There exists a constant C > 0 such that kA0 ((p, ψ))kM′ :=
sup (q,ϕ)∈M
(q,ϕ)6=0
|A0 ((p, ψ), (q, ϕ))| ≥ C k(p, ψ)kM k(q, ϕ)kM
∀ (p, ψ) ∈ M .
(4.8)
In addition, there holds sup |A0 ((p, ψ), (q, ϕ))| > 0
∀ (q, ϕ) ∈ M ,
(p,ψ) ∈ M
8
(q, ϕ) 6= 0 .
(4.9)
Proof. Let us introduce the linear and bounded operator Ξ:M → M
(4.10)
(p, ψ) → Ξ(p, ψ) := ((I − 2R)(p), ψ) .
It is easy to see that Ξ is bijective and Ξ−1 = Ξ. Moreover, since R is a projector, we find that R (I − 2R)(p) = − R(p)
and (I − R) (I − 2R)(p) = (I − R)(p)
∀ p ∈ H(curl; Ω) , (4.11)
which yields A0 ((p, ψ), Ξ(p, ψ)) = b R(p), R(p) + k−2 a−1 curl R(p), curl R(p) + b (I − R)(p), (I − R)(p) + k−2 a−1 curl (I − R)(p), curl (I − R)(p) −
(4.12)
1 1 1 1 h γτ (I − R)(p), ψ i + h γτ (I − R)(p), ψ i + h H0 ψ, ψ i + ı η h ψ, ψ i . 2 2 2 4
It follows that o n aR R(p) Re A0 ((p, ψ), Ξ(p, ψ)) = b R(p), R(p) + k−2 curl R(p), curl |a|2 aR 1 + b (I − R)(p), (I − R)(p) + k−2 (I − R)(p) + h H0 ψ, ψ i , curl (I − R)(p), curl |a|2 2 and hence, according to (2.5), (2.6), and (3.16), we deduce that o n o n Re A0 ((p, ψ), Ξ(p, ψ)) ≥ β kR(p)k2H(curl;Ω) + k(I − R)(p)k2H(curl;Ω) + kψk21/2 ,
(4.13)
with β > 0 depending on k, ǫ0 , ǫ1 , σ1 , and ω. Then, using the triangle inequality we easily obtain kpk2H(curl;Ω) ≤ 2 kR(p)k2H(curl;Ω) + 2 k(I − R)(p)k2H(curl;Ω) , which, together with (4.13), gives n o Re A0 ((p, ψ), Ξ(p, ψ)) ≥ C k(p, ψ)k2M
∀(p, ψ) ∈ M .
Moreover, (4.14) and the boundedness of Ξ imply that n o Re A0 ((p, ψ), Ξ(p, ψ)) ≥ C k(p, ψ)kM kΞ(p, ψ)kM
(4.14)
∀(p, ψ) ∈ M ,
from which the inf-sup condition (4.8) follows straightforwardly. On the other hand, applying the estimate (4.14) with (p, ψ) = Ξ(q, ϕ) for a given (q, ϕ) ∈ M, (q, ϕ) 6= 0, and recalling that Ξ−1 = Ξ, we get o n (4.15) Re A0 (Ξ(q, ϕ), (q, ϕ)) ≥ C kΞ(q, ϕ)k2M > 0 , which provides (4.9) and completes the proof.
Lemma 4.2. The operator K is compact. 9
Proof. Let us decompose K = K1 + K2 + K3 , with Kj : M → M′ , j ∈ {1, 2, 3}, given by [K1 (p, ψ), (q, ϕ)] := 2 (b R(p), R(q)) + (b (I − R)(p), R(q)) + (b R(p), (I − R)(q)) , [K2 (p, ψ), (q, ϕ)] := − 2 k−2 (a−1 curl (I − R)(p), curl (I − R)(q)) − k−2 (a−1 curl R(p), curl (I − R)(q)) − k−2 (a−1 curl (I − R)(p), curl R(q)) , and [K3 (p, ψ), (q, ϕ)] :=
1 1 (H − H0 )ψ, ϕ − h γτ (q), Dψ i − ı η h γτ (q), Sψ i − ı η h ψ, Dϕ i . 2 2
The boundedness of R : H(curl; Ω) → H1 (Ω) (cf. (4.4)) and the compact embedding of H1 (Ω) into L2 (Ω), imply that K1 is compact. In addition, it is clear from (4.2) that curl (I − R) : H(curl; Ω) → L2 (Ω)
q → curl (I − R)(q) =
R
Ω
curl q |Ω|
is a finite rank operator, whence K2 also becomes compact. Finally, the compactness of K3 follows directly from Lemma 3.1. We are now in a position to show the unique solvability of (3.13) (equivalently (3.14) or (4.7)). Theorem 4.3. There exists a unique (p, ψ) ∈ M solution of problem (3.13). Proof. It follows from Lemmas 4.1 and 4.2 that A is a Fredholm operator of index zero and then, the present proof simply reduces to show uniqueness for (3.13). In fact, let (p0 , ψ0 ) be a solution of (3.13) with ui = 0 and define the function u0 (x) if x ∈ Ω, u e(x) := z(x) if x ∈ Ω0 ,
where
u0 (x) := k−2 a−1 curl p0 (x) and z(x) :=
Z
0
2π
∂E (x, x(t)) |x′ (t)| ψ0 (t) dt + ı η ∂ν y
Z
2π
E(x, x(t)) ψ0 (t) dt .
0
From the first equation in (3.13) we easily see that p0 = b−1 curl u0 in Ω, whence curl (b−1 curl u0 ) − k2 a u0 = 0
in Ω
1 and γ(u0 ) = ( I + D) ψ0 + ı η S ψ0 . 2
(4.16)
Further, the second equation in (3.13) gives directly 1 γτ (p0 ) = − H ψ0 − ı η ( I − D t ) ψ0 . 2
(4.17)
On the other hand, it is clear that z solves the Helmholtz equation ∆z + k2 z = 0 in Ω0 10
(4.18)
and satisfies the Sommerfeld radiation condition. In addition, the jump properties of the double layer potential and the normal derivative of the single layer potential through Γ imply (cf. [13]): γ(z) ∂z |x (t)| ◦x ∂ν ′
= =
1 ( I + D) ψ0 + ı η S ψ0 , 2 1 − H ψ0 − ı η ( I − D t ) ψ0 . 2
(4.19)
As a consequence of (4.16), (4.17), (4.18), and (4.19), we deduce that u e solves (2.3) with datum ui = 0 and satisfies the Sommerfeld condition (2.4). In this way, Theorem 2.1 ensures that u e vanishes identically, which yields u0 = 0, and therefore p0 = 0 and 1 ( I + D) ψ0 + ı η S ψ0 = 0 . 2
1 Finally, since ( I + D) + ı η S is one to one (see [13, Theorem 3.3.9]), we obtain that ψ0 = 0. 2
5
The discrete problem
We first construct a triangulation covering exactly the domain Ω. To this end, given N ∈ N, we let 0 = t0 < t1 < · · · < tN = 2π be a uniform partition of [0, 2π] with tj+1 −tj = h = 2π N for j ∈ {0, 1, ..., N −1}, and denote by Ωh the polygonal domain whose vertices are {x(t1 ), x(t2 ), ..., x(tN )}. Then we let Th0 be a regular triangulation of Ωh by triangles T of diameter hT such that h := sup{hT : T ∈ Th0 }. Further, we assume that any vertex of a triangle touching ∂Ωh belongs to {x(t1 ), x(t2 ), ..., x(tN )}. Thus, we replace each triangle T ∈ Th0 with one side along ∂Ωh , by the corresponding curved triangle with one side along Γ. In this way, we obtain from Th0 a triangulation Th of Ω made up of straight and curved triangles. Let T be a curved triangle of Th . We denote its vertices by aT1 , aT2 and aT3 , numbered in such a way that aT2 and aT3 are the endpoints of the curved side of T . Let i ∈ {0, 1, ..., N − 1} be such that x(ti ) = aT2 and x(ti+1 ) = aT3 . Then, ϕ(t) := x (ti + t h), t ∈ [0, 1], is a parameterization of the curved side of T . Next, we let Tb be the reference triangle with vertices b a1 := (0, 0), b a2 := (1, 0) b3 := (0, 1), and consider the affine map GT defined by GT (b and a ai ) = aTi for i ∈ {1, 2, 3}. Also, we introduce the function ΘT : Tb → R2 given by ΘT (b x) :=
x ˆ1 ϕ(ˆ x2 ) − (1 − x ˆ2 )aT2 − x ˆ2 aT3 1−x ˆ2
b ∈ Tb , ∀x
where the limiting value has to be taken as x ˆ2 goes to 1, and define FT := GT + ΘT : Tb → R2 . Now, if T is a straight (interior) triangle, we just take the curving perturbation ΘT ≡ 0 and thus FT is the usual affine map GT from the reference triangle. Lemma 5.1. Let T ∈ Th be a curved triangle. Then, there exists h0 > 0 such that for all h ∈ (0, h0 ) FT is a diffeomorphism of class C ∞ that maps one-to-one Tˆ onto the curved triangle T in such a way ai ) = aTi for all i ∈ {1, 2, 3}. In addition, the Jacobian JT of this mapping does not vanish in that FT (b a neighborhood of Tˆ, and there exist positive constants Ci , i ∈ {1, 2, 3}, independent of T , such that C1 h2T ≤ |JT (b x)| ≤ C2 h2T and max |∂ α ΘT (b x)| ≤ C3 hk+1 T b∈Tb x
∀ α ∈ N2 , 11
b ∈ Tb , ∀x |α| = k ,
(5.1) ∀k ≥ 1.
(5.2)
Proof. See [12] or Theorem 22.4 of [14]. Further estimates concerning curved triangles are collected in the following lemma from [1]. Lemma 5.2. Assume that T is a curved triangle and let BT := D GT . A function v belongs to H 1 (T ) if and only if the function vb = v ◦ FT belongs to H 1 (Tb). In this case, there exist C3 , C4 > 0, independent of T , such that o n 2 (5.3) k kB k kb v k + |b v | |v|1,T ≤ C3 |det BT |1/2 kB−1 T b b T 0,T 1,T , and
|b v |1,Tb ≤ C4 |det BT |−1/2 kBT k kvk1,T .
(5.4)
Proof. See Lemma 2.3 in [1]. We now introduce the first order rotated Raviart-Thomas space W(Tb) ⊆ H(rot; Tb) (cf. [2, 11]), which is given by: a1 x b2 b W(T ) := +b : a1 , a2 , b ∈ C , a2 −b x1
and define
W(T ) :=
n
b ◦ F−1 q = (D FT )−t q T ,
q:
b ∈ W(Tb) q
o
,
where (D FT )−t is the inverse of the transpose of D FT . Then, we seek an approximation of the unknown p in the finite dimensional space Wh := { q ∈ H(curl; Ω) :
q|T ∈ W(T )
∀ T ∈ Th } .
Next, let us denote by Eh : H1 (Ω) → Wh the equilibrium interpolation operator (see [2, 11]). Then, applying (5.1), (5.3), (5.4) and the Bramble-Hilbert Lemma, we deduce that k(I − Eh )(q)kH(curl;Ω) ≤ C h
(
kqk1,Ω +
X
T ∈Th
kcurl qk21,T
!1/2 )
(5.5)
for all q ∈ H1 (Ω) such that (curl q)|T ∈ H 1 (T ), ∀ T ∈ Th . In particular, using Lemmas 5.1 and 5.2, and the estimate (5.5), we can prove the following result. Lemma 5.3. There exists a constant C > 0, independent of h, such that o n k(I − Eh )(q)kH(curl;Ω) ≤ C h kqk1,Ω + kcurl qk0,Ω for all q ∈ H1 (Ω) such that curl q −
R
Ω
curl q |Ω| R
Proof. Let q ∈ H1 (Ω) such that curl q − (5.5), that k(I − Eh )(q)kH(curl;Ω) ≤ C h
(
Ω
kqk1,Ω +
(5.6)
∈ curl (Wh ). curl q |Ω|
= curl r, with r ∈ Wh . It follows, according to
X
T ∈Th
12
kcurl
qk20,T
+ |curl
r|21,T
!1/2 )
,
(5.7)
and then, the rest of the proof reduces to show that |curl r|21,T can be estimated in terms of kcurl qk20,T . Indeed, we first observe that if T is a straight triangle then curl (r|T ) is constant, whence |curl r|1,T vanishes. We now consider a curved triangle T . In this case, (5.3) implies that o n 2 [ [ k kB k k curl rk + | curl r| |curl r|1,T ≤ C3 |det BT |1/2 kB−1 T T 0,Tb 1,Tb , \ and hence, the identity curl (r|T ) =
1 JT
curl b r, together with the fact that curl b r is constant, give
o n −1 2 \rk0,T (b x )| max |J (b x )| kcurl k kB k + max |∇J |curl r|1,T ≤ C |det BT |1/2 kB−1 T T T T ≤ C
|det BT
|1/2
min |JT (b x)|1/2
kB−1 T k
b∈Tb x
n
b∈Tb x
b∈Tb x
x)| max |JT (b x)| kBT k2 + max |∇JT−1 (b b∈Tb x
b∈Tb x
o
kcurl rk0,T .
Next, we deduce from (5.1) and (5.2) that x)| = max max |∇JT−1 (b b∈Tb x
b ∈Tb x
|∇ JT (b x)| ≤ C h−1 T , 2 JT (b x)
2 which, together with the usual estimates for kB−1 T k and kBT k , yields R
Ω curl q
|curl r|1,T ≤ C kcurl rk0,T = C curl q − ≤ C kcurl qk0,T , |Ω| 0,T
(5.8)
where C > 0 is independent of T . Finally, (5.7) and (5.8) provide (5.6) and finish the proof.
It remains to define the finite element subspace for the boundary unknown ψ. In fact, given n ∈ N, we let Tn be the subspace of H 1/2 spanned by the 2π-periodic trigonometric functions of degree not greater than n, that is n n X X bj sin(jt) ∀ t ∈ [0, 2π] , aj , bj ∈ C . aj cos(jt) + Tn := ϕ ∈ H 1/2 : ϕ(t) = j=0
j=1
We note here that the orthogonal projection operator Pn : L2 (0, 2π) → Tn satisfies the following error estimate (cf. Theorem 8.2.1 of [13]) r−s 2 kϕkr ∀ ϕ ∈ Hr , ∀ r ≥ s . (5.9) kϕ − Pn (ϕ)ks ≤ n Now, given h > 0 and n ∈ N, we let Mh,n = Wh × Tn and introduce the following discrete version of (3.14): Find (ph , ψn ) ∈ Mh,n such that A((ph , ψn ), (q, ϕ)) = F((q, ϕ))
∀ (q, ϕ) ∈ Mh,n .
(5.10)
We now aim to prove the well-posedness of (5.10). For this purpose, as established by a classical result on projection methods for Fredholm operators of index zero (see, e.g. [7, Theorem 13.7]), we need to show that the Galerkin scheme associated to the isomorphism A0 is well-posed. Precisely, the following lemma establishes the corresponding discrete inf-sup condition for the bilinear form A0 on the finite element subspace Mh,n . 13
Lemma 5.4. There exist constants C, h0 > 0, independent of h and n, such that for each h ≤ h0 there holds |A0 ((p, ψ), (q, ϕ))| sup ≥ C k(p, ψ)kM ∀ (p, ψ) ∈ Mh,n . (5.11) k(q, ϕ)kM (q,ϕ)∈Mh,n (q,ϕ)6=0
Proof. Let us introduce the linear operator Ξh,n : Mh,n → Mh,n
(5.12)
(p, ψ) 7→ Ξh,n (p, ψ) := ((I − 2 Eh R)(p), ψ) ,
which is the discrete version of Ξ (cf. (4.10)). It follows, applying (5.6), (4.4), and (4.2), that kΞh,n (p, ψ) − Ξ(p, ψ)kM = 2 k(I − Eh ) (R(p))kH(curl;Ω) o n ≤ C h kcurl pk0,Ω , ≤ C h kR(p)k1,Ω + kcurl (R(p))k0,Ω
(5.13)
which, combined with the boundedness of Ξ, implies that Ξh is also bounded, with a constant independent of h and n. Hence, using (4.14), the boundedness of A0 , and (5.13), we find that for each (p, ψ) ∈ Mh,n there holds n o n o n o Re A0 ((p, ψ), Ξh,n (p, ψ)) = Re A0 ((p, ψ), Ξ(p, ψ)) + Re A0 ((p, ψ), (Ξh,n − Ξ)(p, ψ)) ≥ C k(p, ψ)k2M − C0 k(p, ψ)kM kΞh,n (p, ψ) − Ξ(p, ψ)kM ≥ (C − C0 h) k(p, ψ)k2M ,
which yields, for each h ≤ h0 , with h0 sufficiently small, n o Re A0 ((p, ψ), Ξh,n (p, ψ)) ≥ C k(p, ψ)k2M ≥ C k(p, ψ)kM kΞh,n (p, ψ)kM
∀ (p, ψ) ∈ Mh,n .
In this way, the discrete inf-sup condition (5.11) follows straightforwardly from the above estimate. The well-posedness and convergence of the discrete scheme (5.10) can now be established. ˜ 0 ∈ ]0, h0 ] and Theorem 5.5. Let h0 be the constant provided by Lemma 5.4. Then, there exist h ˜ n0 ∈ N such that for each h ≤ h0 and n ≥ n0 , the Galerkin scheme (5.10) has a unique solution (ph , ψn ) ∈ Mh,n . In addition, there exist c1 , c2 > 0, independent of h and n, such that k(ph , ψn )kM ≤ c1
sup (qh ,ϕn )∈Mh,n (qh ,ϕn )6=0
|F((qh , ϕn ))| k(qh , ϕn )kM
and k(p, ψ) − (ph , ψn )kM ≤ c2
inf
(qh ,ϕn )∈Mh,n
k(p, ψ) − (qh , ϕn )kM ,
where (p, ψ) is the unique solution of (3.14). Furthermore, if p ∈ H1 (Ω) and curl p ∈ H 1 (Ω), then there holds σ 2 kψkσ+1/2 k(p, ψ) − (ph , ψn )kM ≤ c3 h kpk1,Ω + h kcurl pk1,Ω + ∀ σ > 0 , (5.14) n where c3 > 0 is also independent of h and n. 14
Proof. We know from Lemmas 4.1 and 4.2, and Theorem 4.3 that the operator A is a compact perturbation of the isomorphism A0 , and that the corresponding problem (3.14) is well posed. In addition, it follows from Lemma 5.4 that the Galerkin scheme associated to A0 is stable and convergent for each h ≤ h0 . Also, it is easy to show that the approximation properties (5.5) and (5.9) and the density of smooth functions in M yield lim inf k(p, ψ) − (qh , ϕn )kM = 0 ∀ (p, ψ) ∈ M . (5.15) h→0
n→∞
(qh ,ϕn )∈Mh,n
Consequently, the proof of the first part is a direct application of [7, Theorem 13.7], whereas the rate of convergence (5.14) follows from the approximation properties (5.5) and (5.9), and the fact that the boundary unknown ψ belongs to H σ for each σ ≥ 0. Indeed, since u solves the homogeneous Helmholtz equation in Ω′c , it becomes analytic in a sufficiently small neighborhood of Γ and hence γ(us ) := γ(u − ui ) ∈ H σ for each σ ≥ 0 (assuming ui sufficiently smooth). Therefore, since 1 1 σ σ −1 (γ(us )) ∈ H σ 2 I + D + ı η S : H → H is an isomorphism, we conclude that ψ = ( 2 I + D + ı η S) for each σ ≥ 0. We remark here that if the solution u of the original problem (2.3) belongs to H 2 (Ω), then it follows easily from (3.3) and (3.4) that p ∈ H1 (Ω), curl p ∈ H 1 (Ω), and kpk1,Ω + kcurl pk1,Ω ≤ C kuk2,Ω , whence the a priori error estimate (5.14) becomes σ 2 k(p, ψ) − (ph , ψn )kM ≤ C h kuk2,Ω + kψkσ+1/2 ∀σ > 0. n
6
Numerical results
In this section we present an example illustrating the performance of our mixed finite element scheme (5.10). We consider the domain Ωc as the annular region determined by the circle centered at the origin with radius 0.7 and the ellipse of minor and major semiaxis 1.3 and 1.5, respectively. The C ∞ closed curve Γ is taken then as this ellipse. In addition, the coefficients a and b are given in Ωc by (1 + |x|2 ) + ı|x|4 if 0.7 < |x| < 1 , 1 + (1 − |x|4 )2 if 0.7 < |x| < 1 , a(x) = and b(x) = 2+ı 1 otherwise , otherwise . (1)
We compute the right hand side in such a way that the Hankel function u(x) := H0 (k|x|) is the exact solution of the problem. In what follows, k2 = 10.50, which is an approximation of the sixth eigenvalue of the Laplacian in Ω with a Neumann boundary condition on Γ. Thus, our method is tested in the case where the numerical scheme given in [6] may exhibit an unstable behavior because of the lack of uniqueness. The numerical results shown below were obtained in a Pentium Xeon computer with dual processors, using a MATLAB code. In Table 1, we take h = 2 π/192 and decrease the spectral parameter n until we obtain the smallest value that preserves the order of accuracy. We can see that the number of degrees of freedom on the boundary is drastically reduced. Hence, for relatively low wave numbers k, it is possible to eliminate the boundary variable at matricial level (by a static condensation procedure) and solve the reduced linear system for the finite element variable. Then, in Table 2 we fix n = 16 and change h in order to 15
check out the asymptotic convergence predicted by the theory for p − ph . We also display there the experimental rate of convergence defined by r(p) :=
logkp − ph kH(curl;Ω) − logkp − ph˜ kH(curl;Ω) , ˜ log(h) − log(h)
˜ denote two consecutive mesh sizes. The computed values of r(p) show that the error where h and h kp − ph kH(curl;Ω) decreases linearly with respect to the mesh parameter h. Moreover, if we denote by Qh the L2 (Ω)-orthogonal projection onto the space of piecewise constant function associated to the triangulation Th , we can use the postprocess computation uh = k−2 Qh (a−1 curl ph )
in
Ω.
to obtain an approximation uh of the scalar variable u. Note from (3.4) that u = k−2 a−1 curl p in Ω. The behavior of the error ku − uh k0,Ω and the corresponding experimental rate of convergence r(u) :=
logku − uh k0,Ω − logku − uh˜ k0,Ω , ˜ log(h) − log(h)
are also depicted in Table 2. 2n 128 64 32 16 8
kp − ph kH(curl;Ω) 2.64 × 10−1 2.64 × 10−1 2.64 × 10−1 2.68 × 10−1 3.86 × 10+0
Table 1: Convergence history for different values of n when k2 = 10.50 and h = 2 π/192.
h 2π/48 2π/64 2π/96 2π/128 2π/192 2π/256
kp − ph kH(curl;Ω) 1.00 × 10+0 8.21 × 10−1 5.31 × 10−1 4.02 × 10−1 2.64 × 10−1 1.99 × 10−1
r(p) — 0.70 1.07 0.96 1.03 0.98
ku − uh k0,Ω 6.04 × 10−3 3.67 × 10−3 1.92 × 10−3 9.91 × 10−4 4.21 × 10−4 2.55 × 10−4
r(u) — 1.72 1.59 2.31 2.10 1.74
Table 2: Convergence history for different values of h when k2 = 10.50 and n = 16.
References [1] C. Bernardi, Optimal finite-element interpolation on curved domains. SIAM Journal on Numerical Analysis, vol. 26, pp. 1212-1240, (1989). [2] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods. Springer, Berlin Heidelberg New York, 1991. 16
[3] A. Buffa, Remarks on the discretization of some noncoercive operator with applications to heterogeneous Maxwell equations. SIAM Journal on Numerical Analysis, vol. 43, pp. 1-18, (2005). [4] D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Second Edition. Springer-Verlag, Berlin, 1998. [5] J. Coyle and P. Monk, Scattering of time-harmonic electromagnetic waves by anisotropic inhomogeneous scatterers or impenetrable obstacles. SIAM Journal on Numerical Analysis, vol. 37, pp. 1590-1617, (2000). [6] G.N. Gatica and S. Meddahi, On the coupling of MIXED-FEM and BEM for an exterior Helmholtz problem in the plane. Numerische Mathematik, vol. 100, pp. 663-695, (2005). [7] R. Kress, Linear Integral Equations, Second Edition. Springer-Verlag, New York, 1999. [8] W. McLean, Strongly Eliptic Systems and Boundary Integral Equations. Cambridge University Press, 2000. ´rquez and V. Selgas, Computing acoustic waves in an inhomogeneous [9] S. Meddahi, A. Ma medium of the plane by a coupling of spectral and finite elements. SIAM Journal on Numerical Analysis, vol. 41, pp. 1729-1750, (2003). [10] S. Meddahi and F. J. Sayas, Analysis of a new BEM-FEM coupling for two-dimensional fluidsolid interaction. Numerical Methods for Partial Differential Equations, vol. 21, 6, pp. 1017-1042, (2005). [11] P. A. Raviart and J. M. Thomas, A mixed finite element method for 2nd order elliptic problems. In: Mathematical Aspects of the Finite Element Method, Lecture Notes in Mathematics, vol. 606, Springer-Verlag, Berlin, 1977. [12] R. Scott, Interpolated boundary conditions in the finite element method, SIAM Journal on Numerical Analysis, vol. 12 , pp. 404-427, (1975). [13] J. Sarannen and G. Vainikko, Periodic Integral and Pseudodifferential Equations with Numerical Approximation. Springer Verlag, 2002. ˇ ´ıˇ [14] A. Zen sek, Nonlinear Elliptic and Evolution Problems and their Finite Element Approximations. Academic Press, London, 1990.
17
