PUBLICATIONS DE L’INSTITUT MATHÉMATIQUE Nouvelle série, tome 96 (110) (2014), 103–123 DOI: 10.2298/PIM1410103F

EXACT NONREFLECTING BOUNDARY CONDITIONS FOR EXTERIOR WAVE EQUATION PROBLEMS

Silvia Falletta and Giovanni Monegato

Abstract. We consider the classical wave equation problem defined on the exterior of a bounded 2D space domain, possibly having far field sources. We consider this problem in the time domain, but also in the frequency domain. For its solution we propose to associate with it a boundary integral equation (BIE) defined on an artificial boundary surrounding the region of interest. This boundary condition is nonreflecting (or transparent) for both outgoing and incoming waves and it does not have to include necessarily the problem datum supports. The problem physical domain can even be a multi-domain, defined by the union of several disjoint domains. These domains can be convex or nonconvex.

This transparent boundary condition is imposed pointwise on the chosen artificial boundary; therefore, its (space collocation) discretization can be coupled with a (space) finite difference or finite element method for the associated PDE problem. In the time-domain case, a classical (explicit or implicit) time integrator is also used. We present a consistency result for the BIE discretization and a sample of the intensive numerical testing we have performed.

1. Introduction

We consider some wave propagation problems defined on unbounded space domains, possibly having far field sources. We will consider these problems in the time domain, but also in the frequency domain. Although we will limit our description to the 2D case, the proposed approaches can be easily extended to 3D problems.

Since usually one has to compute the solution of the given PDE problem only in a bounded area surrounding a physical (bounded) domain, this area is then

2010 Mathematics Subject Classification: 65M60, 65M06, 65M38. Key words and phrases: Wave equation, Helmholtz equation, exterior problems, absorbing

boundary conditions, numerical methods.

This work was supported by the Ministero dell’Istruzione, dell’Università e della Ricerca of Italy, under the research program PRIN09: Boundary element methods for time dependent problems, and by the GNCS-INDAM 2013 research program: Metodi fast per la risoluzione numerica di sistemi di equazioni integro-differenziali.

defined by introducing an artificial outer boundary B, where a nonreflecting (or transparent) condition is imposed. The problem physical domain can even be a multi-domain, defined by the union of several disjoint domains. These domains can be convex or nonconvex. To efficiently deal with these general situations, the shape of the chosen artificial boundary should be the most appropriate one for the given problem. It could be a nonconvex closed curve or even the union of (disjoint) closed curves. It should not necessarily include the problem datum supports, in particular when these are away from the computational domain. Thus an appropriate boundary condition on B should be nonreflecting for both outgoing and incoming waves.

The most commonly used transparent conditions are of local type, hence approximate, to reduce their computational cost. However, they hold only for convex artificial boundaries, having particular shapes; moreover, they do not satisfy most of the above requirements. For a survey of them, see [8, 9].

As already proposed since many years for some elliptic problems, (see, for example, [10, 20]), nonreflecting boundary conditions can be defined by proper boundary integral equations (BIE), defining a relationship between the solution of the differential problem and its normal derivative on B. Such conditions are exact, hence nonlocal, but their computational cost is higher than that of the local ones. Furthermore, till now they have been coupled mainly with finite element methods (FEM) by means of a variational formulation; a coupling with a finite difference (FD) scheme can be found, for example, in [18]. However, these NRBC have all the good properties mentioned above. In the next two sections we will describe a much cheaper approach, which allows the coupling of a BIE also with a finite difference (FD) scheme. Furthermore, as shown in [7], in some cases the computational cost due to this approach can be significantly reduced.

In this paper, which is a continuation of [7], we propose to solve the PDE problem in the chosen spatial region, by taking as boundary conditions the given Dirichlet or Neumann condition on the physical boundary Γ, and the BIE that we impose pointwise on the chosen artificial boundary B. Its (approximate) solution is obtained by applying a finite element or finite difference scheme, associated with a discretization of the chosen (space) computational region, and, in the case of a time-dependent problem, a classical (explicit or implicit) time integrator. At the same time, the BIE condition on B is discretized, with respect to the space variable, by using a classical collocation method. In the case of the time-domain wave equation, since the associated integral equation is of space-time type, we also need to discretize its time integral. This is performed by means of a convolution quadrature (see [5–7, 14]).

In Section 2 we consider the case of the classical wave equation in the time domain and describe the above mentioned approach. In Section 3, the proposed method is applied to the Helmholtz equation. A sample of the intensive numerical testing we have performed is reported in Section 4. This includes also cases of multi-scattering and of artificial boundaries with corners. We discuss all the computational issues that make the proposed approach efficient, and compare this later with some of those of local type.

2. Time-domain problems

We consider the exterior problem for the classical (nonhomogeneous) wave equation with constant speed, that for simplicity we assume to be unitary, with a Dirichlet boundary condition; the Neumann case is very similar. Thus, in the case of the Dirichlet condition this is:

(2.1)

u_{tt}(x, t) − Δu(x, t) = f(x, t) in Ω^{e} × (0, T)
u(x, t) = g(x, t) in Γ × (0, T)
u(x, 0) = u_{0}(x) in Ω^{e}
u_{t}(x, 0) = v_{0}(x) in Ω^{e},

where Ω^{e} is the exterior region defined by a bounded domain Ω^{i}, whose boundary is denoted by Γ and is assumed to be smooth. The problem data are assumed to have local supports and to satisfy the smoothness and compatibility conditions required by the theoretical results we will use and by the numerical methods we will apply. Nevertheless, in one example we will consider also the case of a concentrated source f.

The domain Ω^{i} can also be a multi-domain, that is, a union of a finite number of disjoint bounded domains Ω_{j}^{i}, each one having its own boundary Γ_{j}. In this case we define Ω^{i} = ⋃ _{j} Ω_{j}^{i} and Γ = ⋃ _{j} Γ_{j}. The boundaries Γ and Γ_{j} do not have to be necessarily convex (closed) curves.

To solve the above problem, we introduce an artificial boundary B which divides the original infinite domain Ω^{e} into two domains: the finite computational domain Ω of interest, and an infinite residual domain D. In the multi-domain case, each element Ω_{j}^{i} is surrounded by an artificial boundary B_{j}, with ⋂ _{j} B_{j} = ∅, and defines B = ⋃ _{j} B_{j}.

Then, by analyzing the problem in D, we obtain a nonreflecting relation on B between u and its (outward) normal derivative. This relation is used as a boundary condition on B, to obtain a well posed problem in Ω. Finally, this problem is solved in Ω by coupling a standard numerical method (finite differences or finite elements), for the space approximation, with a time integrator.

In [7] the following space-time boundary integral equation has been used to define a nonreflecting boundary condition for the Dirichlet problem:

\frac{1}{2}u(x, t) = Vλ(x, t) − Ku(x, t) + I_{u_{0}} (x, t) + I_{v_{0}} (x, t) − I_{f} (x, t), x ∈ B

where: λ(x, t) = ∂_{n}u(x, t) (n = n_{x} is the unit outward normal vector with respect to D at the point x ∈ B), and

(2.2)  Vλ(x, t) = ∫_{0}^{t} ∫_{B} G(x − y, t − τ)λ(y, τ) dB_{y}dτ
(2.3)  Ku(x, t) = ∫_{0}^{t} ∫_{B} ∂_{n_{y}} G(x − y, t − τ)u(y, τ) dB_{y}dτ,

I_{u_{0}} (x, t) = \frac{∂}{∂t} ∫_{D} u_{0}(y)G(x − y, t) dy, I_{v_{0}} (x, t) = ∫_{D} v_{0}(y)G(x − y, t) dy

(2.4)  I_{f} (x, t) = ∫_{0}^{t} ∫_{D} f(y, τ)G(x − y, t − τ) dy dτ.

The function

G(x, t) = \frac{1}{2π} \frac{H(t − ∥x∥)}{\sqrt{t^{2} − ∥x∥^{2}}}

is the fundamental solution of the wave equation.

This BIE defines an integral relationship that u and λ = ∂_{n_{x}} u must satisfy if no spurious reflections have to be produced. It represents a (global) boundary condition, which is nonreflecting for both outgoing and incoming waves. Its operators V and K are bounded when acting between the following spaces (see [11]):

V : H_{0}^{r+1} (0, T; H^{−1/2} (B)) → H_{0}^{r}(0, T; H^{1/2} (B)) (2.5)

K : H_{0}^{r+3/2} (0, T; H^{1/2} (B)) → H_{0}^{r}(0, T; H^{1/2} (B))

for r ⩾ 0. Thus, in the domain of interest Ω the problem reformulates

(2.6)

u_{tt}(x, t) − Δu(x, t) = f(x, t) in Ω × (0, T)
u(x, t) = g(x, t) in Γ × (0, T)
\frac{1}{2}u(x, t) − Vλ(x, t) + Ku(x, t) = I_{u_{0}} (x, t) + I_{v_{0}} (x, t) + I_{f} (x, t),
(x, t) ∈ B × (0, T)
u(x, 0) = u_{0}(x) in Ω
u_{t}(x, 0) = v_{0}(x) in Ω

Remark 2.1. Each of the above “volume” integrals I_{u_{0}} , I_{v_{0}} , I_{f} is present only if the (local) support of its datum is not in Ω (for its numerical evaluation see [5, 7]). When this occurs, the datum must be replaced by the trivial one in the corresponding equation of (2.6) defined in Ω.

Remark 2.2. In the case of a Neumann boundary condition on Γ, it is sufficient to replace in (2.6) the Dirichlet condition u(x, t) = g(x, t) on Γ by ∂_{n_{x}} u(x, t) = g(x, t), while the nonreflecting condition on B remains unaltered.

As shown in [7, Remark 2.5], it is straightforward to show that if the original problem (2.1) has a unique solution, for example in C^{1}([0, T]; H^{1}(Ω^{e})), then also (2.6) has a unique solution in C^{1}([0, T]; H^{1}(Ω)).

Thus at this point we solve a new problem (2.6) by applying to it a numerical method. This is obtained by coupling a proper discretization of our NRBC boundary condition with a finite difference or finite element method. Since the construction of a finite difference or finite element method is standard, we briefly recall the main steps which lead to the discretization of the BIE. For more details see [5].

After having defined a uniform partition: t_{n} = nΔ_{t}, n = 0, …, N, Δ_{t} = T/N, of the integration time interval [0, T], we discretize the convolution time integrals in (2.2) and (2.3) by using the Lubich convolution quadrature rule [13]

∫_{0}^{t_{n}} K(t_{n} − τ) φ(τ) dτ ≈ ∑_{j=0}^{n} w_{n−j} (Δ_{t})φ(t_{j}), n = 0, …, N

where

w_{m}(Δ_{t}) = \frac{1}{2πı}∫_{|z|=ρ} K̂ ( \frac{γ(z)}{Δ_{t}}) z^{−(m+1)} dz;

K̂(s) is the Laplace transform of the kernel K(t) and ρ is a (small) properly chosen positive real. By introducing the polar coordinate z = ρe^{i̇φ}, and applying the trapezoidal rule with L = 2N equal steps of length 2π/L to the corresponding integral, we obtain

w_{m}(Δ_{t}) ≈ \frac{ρ^{−m}}{L} ∑_{l=0}^{L−1} K̂ ( \frac{γ(ρ exp(i̇l2π/L))}{Δ_{t}} ) exp(−i̇ml2π/L).

By using the FFT, all the coefficients w_{m}, m = 0, …, N, can be simultaneously computed with O(N logN) flops.

In the following, we will denote by w_{m}^{V} and w_{m}^{K} the quadrature rule coefficients associated with the convolution kernels in (2.2) and (2.3), respectively. In general, the presence of the upper indices V, K means that the corresponding quantities are associated with the operators defined in (2.2) and (2.3).

For the space discretization, we describe only the case of a single domain Ω^{i}, being straightforward its generalization to multiple scattering. Thus, first we introduce a parametrization of the curve B, x = ψ (x) = (ψ_{1}(x), ψ_{2}(x)) and y = ψ (y) = (ψ_{1}(y), ψ_{2}(y)), with x, y ∈ [a, b]. Notice that this requirement is not a restriction. Indeed, since the contour B can be arbitrarily chosen, we can always define a smooth parametric curve having the desired shape.

At every time instant t_{j} we approximate the (unknown) function u(ψ (x), t_{j}) and its normal derivative λ(ψ (x), t_{j}) by continuous piecewise linear interpolants, associated with a uniform partition {x_{k}}_{k=1}^{M+1} of the parametrization interval [a, b], whose size is defined by the quantity Δ_{x} = max_{1⩽k⩽M} (x_{k+1} − x_{k}). These interpolants are written in the form

u(ψ (x), t_{j}) ≈ u_{Δ_{x}}^{ψ} (x, t_{j}) := ∑_{i=1}^{M+1} u_{i}^{j} N_{i}(x), u_{i}^{j} = u(ψ (x_{i}), t_{j})

λ(ψ (x), t_{j}) ≈ λ_{Δ_{x}}^{ψ} (x, t_{j}) := ∑_{i=1}^{M+1} λ_{i}^{j} N_{i}(x), λ_{i}^{j} = λ(ψ (x_{i}), t_{j})

where {N_{i}(x)} are the classical Lagrangian basis functions of local degree 1, u_{1}^{j} = u_{M+1}^{j} and λ_{1}^{j} = λ_{M+1}^{j} . These in turns define the associated interpolants of u(x, t_{j}) and λ(x, t_{j}) on the curve B:

u(x, t_{j}) ≈ u_{Δ_{x}} (x, t_{j}) := ∑_{i=1}^{M+1} u_{i}^{j} b_{i}(x), u_{i}^{j} = u(ψ (x_{i}), t_{j})

λ(x, t_{j}) ≈ λ_{Δ_{x}} (x, t_{j}) := ∑_{i=1}^{M+1} λ_{i}^{j} b_{i}(x), λ_{i}^{j} = λ(ψ (x_{i}), t_{j})

with b_{i}(x) = b_{i}(ψ (x)) = N_{i}(x).

Since the role of the NRBC is to define on B a relationship between the (outgoing/incoming) wave and its normal derivative, which prevents the raising of (spurious) incoming/outgoing waves, the more accurate is the discretized relationship the more transparent this will be. To this end, having chosen a continuous piecewise linear (space) approximant for u(ψ (x), t_{j}), we use an approximant of the same type also for λ(ψ (x), t_{j}). Then, we collocate the BIE at the boundary mesh points x_{m} for all m = 1, …, M, n = 0, …, N. We obtain

(Ku)(x_{m}, t_{n}) ≈ ∑_{j=0}^{n} ∑_{i=1}^{M} u_{i}^{j} ∫_{B} w_{n−j}^{K} (Δ_{t}; ∥x_{m} − y∥)b_{i}(y) dB_{y},

(Vλ)(x_{m}, t_{n}) ≈ ∑_{j=0}^{n} ∑_{i=1}^{M} λ_{i}^{j} ∫_{B} w_{n−j}^{V} (Δ_{t}; ∥x_{m} − y∥)b_{i}(y) dB_{y}.

For the truncation errors associated with the discrete operators

(K_{Δ} u)(x, t_{n}) := ∑_{j=0}^{n} ∫_{B} w_{n−j}^{K} (Δ_{t}; ∥x − y∥)u_{Δ_{x}} (y, t_{j}) dB_{y}

(V_{Δ} λ)(x, t_{n}) := ∑_{j=0}^{n} ∫_{B} w_{n−j}^{V} (Δ_{t}; ∥x − y∥)λ_{Δ_{x}} (y, t_{j}) dB_{y}

we have obtained in [7] the following consistency estimates.

Proposition 2.1. Under the assumption u ∈ C^{4}([0, T]; H^{r}(D)), r ⩾ 7/2, we have

max_{0⩽n⩽N} ∥(K − K_{Δ} )u(·, t_{n})∥_{L^{2} (B)} = O(Δ_{t}^{2}) + O(Δ_{x}^{3/2} ),
max_{0⩽n⩽N} ∥(V − V_{Δ} )λ(·, t_{n})∥_{L^{2} (B)} = O(Δ_{t}^{2}) + O(Δ_{x}^{2}).

Similar bounds can be derived under a milder smoothness assumption on u, with respect to the space variable. In this case, as it will be pointed out in the proof, from the theoretical point of view it is sufficient to approximate λ(ψ (x), t_{j}) by a piecewise constant function. See however the remark made after the proof of the next proposition.

Proposition 2.2. Let u ∈ C^{4}([0, T], H^{1}(D)), and assume that for any given t, u is a piecewise H^{5/2} (D) function, whose normal derivative on B has at most finite jumps at a finite number of points z_{ℓ} = z_{ℓ}(t), ℓ = 1, …, ν(t) ⩽ ν_{T}. Then, we have

max_{0⩽n⩽N} ∥(K − K_{Δ} )u(·, t_{n})∥_{L^{2} (B)} = O(Δ_{t}^{2}) + O(Δ_{x}),
max_{0⩽n⩽N} ∥(V − V_{Δ} )λ(·, t_{n})∥_{L^{2} (B)} = O(Δ_{t}^{2}) + O(Δ_{x}^{1/2} ).

Proof. For any given t = t_{n}, n = 0, …, N, we rewrite the operator discretization errors as

(K − K_{Δ} )u = (Ku − Ku_{Δ_{x}} ) + (Ku_{Δ_{x}} − K_{Δ} u) =: R_{1}^{K} + R_{2}^{K}
(V − V_{Δ} )λ = (Vλ − Vλ_{Δ_{x}} ) + (Vλ_{Δ_{x}} − V_{Δ} λ) =: R_{1}^{V} + R_{2}^{V}

Taking into account the mapping properties of the operators K, V (see [7]) we obtain the bounds

(2.7)  ∥R_{1}^{K}∥_{H^{1/2} (B)} ⩽ C∥u − u_{Δ_{x}} ∥_{H^{1/2} (B)}
(2.8)  ∥R_{1}^{V}∥_{H^{1/2} (B)} ⩽ C∥λ − λ_{Δ_{x}} ∥_{H^{−1/2} (B)} ⩽ C∥λ − λ_{Δ_{x}} ∥_{L^{2} (B)}

where, here and in the following, the constants C, which in general will take different values on different occurrences, depends only on T, not on t_{n}. To obtain the behaviors, in terms of Δ_{x}, of the last bounds in (2.7) and (2.8), we perform some fairly standard steps.

First we introduce in bound (2.7) the (smooth) curve parametrization ψ (x), which reduces the integration over the curve B to an equivalent integral defined on an interval I = [a, b], that above we have subdivided into subintervals I_{k} = (x_{k}, x_{k+1} ), whose maximum length was denoted by Δ_{x}. For notational simplicity, for any given instant t_{n} we set ũ(x) = u(ψ (x), t_{n}) and λ̃(x) = λ(ψ (x), t_{n}); therefore all the functions defined below in the proof will also depend on t_{n}.

Then, we consider the intervals I_{k} that do not contain a z_{ℓ} = z_{ℓ}(t_{n}) point, that is, where ũ ∈ H^{2}(I_{k}). Note that in these intervals we have ũ ∈ C^{1}(\overline{I}_{k}). In each of these intervals we set e_{k} = ũ − Π_{k}^{1} ũ, Π_{k}^{1} ũ being the linear polynomial which interpolates the function ũ at the endpoints of the interval I_{k}, we have e_{k}(x_{k}) = e_{k}(x_{k+1} ) = 0. Therefore, under the smoothness assumption we have made on u, there exists a point ξ_{k} ∈ I_{k} where e′_{k}(ξ_{k}) = 0. Thus,

e′_{k}(y) = ∫_{ξ_{k}}^{y} e′′_{k}(x) dx = ∫_{ξ_{k}}^{y} ũ′′(x) dx, ∀y ∈ [x_{k}, x_{k+1} ]

wherefrom we obtain

|e′_{k}(y)| ⩽ ∫_{x_{k}}^{x_{k+1}} |ũ′′(x)|dx ⩽ Δ_{x}^{1/2} ( ∫_{x_{k}}^{x_{k+1}} |ũ′′(x)|^{2}dx) ^{1/2}

that is,

(2.9)  ∫_{x_{k}}^{x_{k+1}} |e′_{k}(y)|^{2}dy ⩽ Δ_{x}^{2} ∫_{x_{k}}^{x_{k+1}} |ũ′′(x)|^{2}dx.

Since

|e_{k}(x)| = |∫_{x_{k}}^{x} e′_{k}(y) dy| ⩽ |∫_{x_{k}}^{x_{k+1}} e′_{k}(y) dy|

from (2.9) we also get

|e_{k}(x)| ⩽ Δ_{x}^{3/2} ( ∫_{x_{k}}^{x_{k+1}} |ũ′′(x)|^{2}dx) ^{1/2}

wherefrom

∫_{x_{k}}^{x_{k+1}} |e_{k}(x)|^{2}dx ⩽ Δ_{x}^{4} ∫_{x_{k}}^{x_{k+1}} |ũ′′(x)|^{2}dx.

Summing all these terms we obtain

(2.10)  ∑_{k,z_{ℓ} ∉I_{k}} ∫_{x_{k}}^{x_{k+1}} |e_{k}(x)|^{2}dx ⩽ Δ_{x}^{4} ∑_{k,z_{ℓ} ∉I_{k}} ∫_{x_{k}}^{x_{k+1}} |ũ′′(x)|^{2}dx,

(2.11)  ∑_{k,z_{ℓ} ∉I_{k}} ∫_{x_{k}}^{x_{k+1}} |e′_{k}(x)|^{2}dx ⩽ Δ_{x}^{2} ∑_{k,z_{ℓ} ∉I_{k}} ∫_{x_{k}}^{x_{k+1}} |ũ′′(x)|^{2}dx ⩽ CΔ_{x}^{2}.

For at most ν_{T} intervals I_{k} containing a point z_{ℓ} where ũ′(x) has a finite jump, a simple calculation gives

(2.12)  ∑_{k,z_{ℓ}∈I_{k}} ∫_{x_{k}}^{x_{k+1}} |e_{k}(x)|^{2}dx ⩽ CΔ_{x}^{3},

(2.13)  ∑_{k,z_{ℓ}∈I_{k}} ∫_{x_{k}}^{x_{k+1}} |e′_{k}(x)|^{2}dx ⩽ CΔ_{x}.

Adding (2.12) to (2.10), and (2.13) to (2.11), we obtain

∥u − u_{Δ_{x}} ∥_{L^{2}} = O(Δ_{x}^{3/2} ) and ∥u − u_{Δ_{x}} ∥_{H^{1}} = O(Δ_{x}^{1/2} )

uniformly with respect to t_{n}. Recalling the well known interpolation inequality ∥ · ∥_{H^{1/2}} ⩽ \sqrt{∥ · ∥_{L^{2}} ∥ · ∥_{H^{1}}}, hence taking advantage of the inequality ∥ · ∥_{L^{2}} ⩽ ∥ · ∥_{H^{1/2}} , we finally have max_{0⩽n⩽N} ∥R_{1}^{K}∥_{L^{2}} = O(Δ_{x}).

In the case of R_{1}^{V}, we note that λ̃(x) is a piecewise H^{1} function, having finite jumps at the abscissas z_{ℓ}. Thus in this case it is more natural to interpolate λ̃(x), at the mesh points x_{k}, by a piecewise constant function. For the associated error ϵ_{k} = λ̃ − Π_{k}^{0} λ̃ we have

(2.14)  ∑_{k,z_{ℓ} ∉I_{k}} ∫_{x_{k}}^{x_{k+1}} |ϵ_{k}(x)|^{2}dx ⩽ CΔ_{x}^{3},

∑_{k,z_{ℓ}∈I_{k}} ∫_{x_{k}}^{x_{k+1}} |ϵ_{k}(x)|^{2}dx ⩽ CΔ_{x},

from which it follows max_{0⩽n⩽N} ∥R_{1}^{V}∥_{L^{2}} = O(Δ_{x}^{1/2} ).

To bound the error terms R_{2}^{K} we first rewrite it as

R_{2}^{K}(x, t_{n}) = ∫_{B} ∑_{i=1}^{M+1} E_{i}^{K}(x, y, t_{n})N_{i}^{B}(y) dB_{y}

having set

E_{i}^{K}(x, y, t_{n}) = ∫_{0}^{t_{n}} K^{K}(∥x −y∥; t_{n} −τ)λ(y_{i}, τ) dτ −∑_{j=0}^{n} ω_{n−j}^{K} (Δ_{t}; ∥x −y∥)λ(y_{i}, t_{j}),

where y_{i} = ψ (y_{i}).

We remark that at any point x ∈ B, until an incoming or outgoing wave has not reached this point, we have u(x, t) = λ(x, t) ≡ 0 in an interval [0, t_{0}], for some t_{0} > 0. Therefore, also the compatibility conditions required by [14, Theorem 3.1] are satisfied. Thus, recalling the error estimates derived in this theorem, under the smoothness assumption we have made we immediately obtain the (uniform) bound |E_{i}^{K}(x, y, t_{n})| ⩽ C_{T}^{B}Δ_{t}^{2}, where the constant C_{T}^{B} does not depend on i, x, y and n. Since N_{i}^{B}(y) = N_{i}(y) ⩾ 0 and ∑ _{i=1}^{M+1} N_{i}(y) = 1, we finally have

max_{0⩽n⩽N} ∥R_{2}^{K}∥_{L^{2} (B)} ⩽ CΔ_{t}^{2}.

The error term R_{2}^{V} can be bounded by proceeding as we did for R_{2}^{K}. The only difference concerns the term

E_{i}^{V}(x, y, t_{n}) = ∫_{0}^{t_{n}} K^{V}(∥x − y∥; t_{n} − τ)λ(y_{i}, τ) dτ − ∑_{j=0}^{n} ω_{n−j}^{V} (Δ_{t}; ∥x − y∥)λ(y_{i}, t_{j}).

In this case, to apply [14, Theorem 3.1], we set K^{V}(r, t) = r^{−ϵ} [r^{ϵ}K^{V}(r, t)] = r^{−ϵ} K_{1}^{V}(r, t), where ϵ > 0 can be taken arbitrarily small. Note that if K̂^{V}(rs) is the Laplace transform of K^{V}, then r^{ϵ}K̂_{1}(rs) is the Laplace transform of K_{1}. This latter transform now satisfies the requirement of the above mentioned theorem. Furthermore, we also have ω_{n−j}^{V} (Δ_{t}; r) = r^{−ϵ} ω_{n−j}^{V_{1}} (Δ_{t}; r), where V_{1} denotes operator (2.2) with kernel K_{1}. Thus, we can apply the above mention theorem and obtain

|E_{i}^{V}(x, y, t_{n})| ⩽ \frac{C_{T}^{B}}{∥x − y∥^{ϵ}} Δ_{t}^{2},

where the constant C_{T}^{B} does not depend on i, x, y and n. Since N_{i}^{B}(y) = N_{i}(y) ⩾ 0 and ∑ _{i=1}^{M+1} N_{i}(y) = 1, we finally have max_{0⩽n⩽N} ∥R_{2}^{V}∥_{L^{2} (B)} ⩽ CΔ_{t}^{2}. □

Remark 2.3. We note that in Proposition 2.2, the maximum number of intervals I_{k} containing points of B, where the solution u is assumed to be less smooth, is bounded by a number ν_{T} that does not depend on Δ_{x}. In general, it is much smaller than that of the remaining partition intervals. Actually, the number of the latter ones tends to infinity, as Δ_{x} → 0, and the error contribution due to these intervals becomes dominant. Thus, as Δ_{x} decreases, but does not take sufficiently small values, the expected error behavior, with respect to the space discretization, is O(Δ_{x}^{3/2} ) for ∥u − u_{Δ_{x}} ∥_{H^{1/2} (B)} and O(Δ_{x}) for ∥λ − λ_{Δ_{x}} ∥_{L^{2} (B)} .

In the latter case, if we assume that u is a piecewise H^{7/2} (D) function, so that λ̃ is a piecewise H^{2}(B) function, and we interpolate, as we did for ũ, λ̃ by a piecewise linear function, the error bound (2.14) can be replaced by O(Δ_{x}^{4}), and the above O(Δ_{x}) estimate for ∥λ − λ_{Δ} ∥_{L^{2} (B)} by O(Δ_{x}^{3/2} ).

Finally, using the vectorial notation, and setting

V_{m,i}^{j} := ∫_{B} w_{j}^{V}(Δ_{t}; ∥x_{m} − y∥)b_{i}(y) dB_{y},
K_{m,i}^{j} := ∫_{B} w_{j}^{K}(Δ_{t}; ∥x_{m} − y∥)b_{i}(y) dB_{y},

the full discretization of our NRBC takes the form

\frac{1}{2}u^{n} − ∑_{j=0}^{n} V^{n−j} λλλ^{j} + ∑_{j=0}^{n} K^{n−j} u^{j} = I_{u_{0}}^{n} + I_{v_{0}}^{n} + I_{f}^{n}, n = 1, …, N

This is then coupled with the system of equations generated by the FD or FE space discretization of the PDE and by the chosen time ODE integrator (see [7]).

Remark 2.4. To reduce the overall computational cost, an FFT routine is used for the simultaneous computation of all the elements of the matrices V^{n} and K^{n}. A further significant CPU reduction can be obtained when the chosen artificial boundary is a circle and we choose a uniform partition on it. In this case, the matrices V^{j} and K^{j} have a Toeplitz structure (therefore only one row must be constructed and stored for each matrix, thus saving space memory and CPU time). For more complex geometries, V^{j} and K^{j} can be approximated by sparse matrices (see [7]).

3. Frequency domain

When the wave equation problem considered in the previous section is reformulated in the frequency domain, and for simplicity we take f = u_{0} = v_{0} ≡ 0, this is reduced to an (exterior) boundary value problem for the Helmholtz equation of the following type:

(3.1)

Δû(x) + k^{2} û(x), = 0  x ∈ Ω^{e}
û(x) = ĝ_{k}(x), x ∈ Γ

where û satisfies the Sommerfeld radiation condition at infinity.

For the solution of this problem, many papers have been written in the last thirty years. Among them, we recall the variational coupling of FEM and boundary element methods (BEM). For a few examples of this approach, see [12, 15, 20]. In general, the NRBC is imposed in a weak form.

To solve problem (3.1), we consider the NRBC given by the following boundary integral relationship

\frac{1}{2} û(x) = Vλ̂(x) − Kû(x) x ∈ B,

where

Vψ(x) := ∫_{B} G(x − y)ψ(y) dB_{y},
Kφ(x) := ∫_{B} ∂_{n_{y}} G(x − y)φ(y) dB_{y}, G(x) = \frac{1}{2π} K_{0}(ωx),

which is the stationary analogue of the one we have used in the previous section. Contrary to what that has been generally done in the literature, when dealing with NRBC of global type, we will impose it in its strong form, i.e., pointwise. This is to reduce the computational cost of its discretization, that we will perform by applying a classical collocation boundary element method.

The mapping properties of the new operators V and K are well known (see for example, [19]). Actually, those we have reported in (2.5) are obtained from the corresponding ones derived for the Helmholtz case we are considering in this section. The latter ones are V : H^{−1/2} (B) → H^{1/2} (B), K : H^{1/2} (B) → H^{1/2} (B) Therefore, also for their discretization, which coincides with the space discretization of the associated time dependent case described in the previous section, we have the same consistency estimates we have obtained in Propositions 2.1 and 2.2 for the space discretization. Also in this case Remark 2.4 applies.

The reformulation of the above Helmholtz problem, in the chosen computational domain Ω, becomes

Δû(x) + k^{2} û(x) = 0  in Ω
û(x) = ĝ_{k}(x) on Γ
\frac{1}{2} û(x) − Vλ̂(x) + Kû(x) = 0  on B

As in the time dependent case, we solve this problem by coupling a classical finite difference or finite element method with a classical (collocation) discretization of the NRBC. In the next section we present some examples that have been solved by applying this approach.

Unfortunately, as for the time dependent case, no stability and convergence (theoretical) results are known. On the contrary, for the variational coupling of FEM and BEM, stability and convergence results have been derived (see, for example, [15,20]). Nevertheless, as for the time dependent case, the proposed approach has shown to be very efficient both for small and relatively large values of the wave number k. We have had good results even when k is extremely close (15 significant digits) to a value which is critical for the invertibility of the BIE (see [2]).

4. Numerical results

4.1. Examples of time domain problems. In this section, we present some examples of the numerical testing we have performed for solving time domain problems, by using the approach discussed in Section 2. To measure the accuracy of the approximations we construct, we take a reference “exact” solution obtained by applying the Lubich-collocation boundary element method described in [5] with a very fine discretization. Once the density function is retrieved, the solution at any point in the infinite domain Ω^{e} is defined by computing the associated potential (see [5] for details). In the following, this solution will be denoted by the acronym BEM. To confirm its accuracy, we have first taken a sufficiently large artificial boundary and set on it the trivial Dirichlet condition, hence applied on the new computational domain the Crank–Nicolson/FE method (see [7]), with sufficiently fine time step and triangulation. The agreement of the solutions produced by these two approaches has been verified in all the examples which we will present. We have however preferred the former, because, to obtain a sufficiently high relative accuracy (≈ 1E − 6), needed to estimate the accuracy of the approximants, the FE method requires an excessively fine space triangulation.

When the chosen artificial boundary B is a circle, we compare the solution obtained by using the exact NRBC (in the following denoted by the acronym ABC

Figure 1. Example 1. FEM: triangulation of the annulus nt = 214, nB = 34 (left plot) and approximate solutions for Δ_{t} = 40/128 (right plot).

Figure 2. Example 1. FD: approximate solutions for Δ_{t} = 40/128 (left plot) and Δ_{t} = 200/640 (right plot).

Lubich) with the classical (and cheap) first order Engquist–Majda NRBC (see [4]) and with the Sommerfeld NRBC (see [8]). In these examples, we have also applied the second order Engquist–Majda [4] and Bayliss–Turkel [3] methods. However, when using the latter NRBC, because of the discrepancy between the boundary B_{Δ} of the FE computational domain and that (B) of the NRBC, and of the presence of the tangential derivative, the solutions we have obtained are not satisfactory (they are very much underestimated). Note that while on B this derivative is smooth, on B_{Δ} it is only piecewise continuous. Indeed, when we rewrite the original PDE problem in polar coordinates, and apply the finite-difference method for the space discretization, the results produced by the second order Engquist–Majda and Bayliss–Turkel methods are the expected ones.

Example 1. As a first example, we apply our numerical scheme to the homogeneous case of Problem (2.1): the source f and the initial data u_{0} and v_{0} are zero throughout the infinite exterior domain Ω^{e}. The boundary Γ is the circle of radius r = 0.25, where we prescribe the Dirichlet condition g(x, t) = t^{3}e^{−0.05(x_{1}^{2} +x_{2}^{2} −\sqrt{2}t) ^{2}} for all t ⩾ 0. Clearly, the solution of this problem is a radial function.

We first choose a circular artificial boundary with radius R = 0.5, so that Ω is the annulus bounded internally by Γ and externally by B. For the finite element discretization, we perform an (approximate) domain triangulation in the cartesian coordinates. This means that the Ω domain is itself approximated by an inscribed polygon. Instead, for the discretization of our NRBC we use the parametric representation given for the boundary B. The boundary mesh is defined by the boundary points of the above domain triangulation, which in our case turns out to be (slightly) nonuniform. For the space discretization we choose an unstructured triangular mesh of n_{t} triangles, having n_{B} (not equally spaced) points on the boundary B. We apply the Crank–Nicolson scheme in time (see [7]).

Figure 3. Example 1. Solution at P = ( \frac{\sqrt{2}}{4}, \frac{\sqrt{2}}{4} ), Δ_{t} = 40/256

In Figure 1, we plot the triangulation of the annulus (left plot) and the behavior of the solution at a mesh point x ∈ B, in the time interval [0, 40]. In particular, in the right plot we compare the approximants produced by the different NRBC we have considered, taking n_{t} = 214 triangles in Ω and N = 128, with the exact solution.

We also apply the finite-difference method for the space discretization, by using polar coordinates to transform the original domain into the rectangle [r, R]×[0, 2π]. In Figure 2, left plot, we compare the solution obtained by using our NRBC with the one produced by the second order Bayliss–Turkel artificial boundary conditions, both for the discretization of the transformed domain into 16 × 16 subrectangles, and choosing a uniform partition of the time interval [0, 40] into 128 subintervals. With the same space refinement, in the right plot, we show the solution obtained by applying the exact NRBC for the final time T = 200. This last example shows that our NRBC produces an accurate solution even for large time integration intervals.

As already remarked, the proposed NRBC allows the treatment of incoming and outgoing waves and has the property of being suitable for artificial boundaries of general shapes, even of nonconvex type, and having also corners, if necessary. Indeed, for some geometries of the physical domain boundary Γ, or of the domain of interest Ω, the choice of a circular artificial boundary B can be wasteful both from the computational and space memory point of view. For the same problem, we consider here, for example, a square and a nonconvex L-shape artificial boundaries. In both cases corners are present in B. In this situation, the boundary mesh nodes coinciding with corner points are slightly shifted to the left.

In Figure 3 we show the bounded computational domains and the behavior of the solutions, obtained by applying the new NRBC, at the point P = (\sqrt{2}/4, \sqrt{2}/4). Note that P coincides with a corner of the square artificial boundary, and with an inner corner of the L-shape artificial boundary.

Example 2. In this second example we consider nontrivial data: in particular, for simplicity, we choose u_{0} = 0, v_{0} = 0, and f ≠ 0. In this case, the artificial boundary B is chosen in such a way that the source f is locally supported in the residual infinite domain D, while it is zero in Ω. Therefore, the artificial boundary condition reads:

\frac{1}{2}u(x, t) − Vλ(x, t) + Ku(x, t) = I_{f} (x, t) in B × (0, T].

In particular, we consider a source concentrated at a point x_{0}: f(x, t) = δ(x − x_{0}) × sin(5t), having a periodic constant oscillatory behavior. With this choice, the volume integral I_{f} (see (2.4)) has the following simpler form:

I_{f} (x, t) = ∫_{0}^{t} sin(5τ)G(x − x_{0}, t − τ) dτ.

For the computation of the volume integral I_{f} , we apply the Lubich convolution rule.

We have compared the solution obtained by the above mentioned approach and the usual one, which consists of including the source in the finite computational domain Ω. We recall that, if the source is away from the area of interest, the latter approach would require a much larger domain Ω, thus wasting computational time and space memory. In our test we have placed the source f at x_{0} = (10, 0); Γ and B are the circles of radius r_{0} = 1 and R, respectively, both centered at the origin. We have chosen first R = 2 (Ω = Ω_{1}) and then R = 12 (Ω = Ω_{2}). In the first case, f is external to the finite computational domain, while in the second case it is included in the annulus bounded by Γ and B. The solutions produced by the two approaches are very similar. For simplicity, in Figure 4 we have reported the plot of the solution we have obtained at the mesh point x_{P_{1}} = (1.9995, 0.0436) in the case of the external source (R = 2), in the time interval [0, 40] (left plot) and [0, 160] (right plot).

Figure 4. Example 2. The external source case: T = 40, N = 512 (left plot), T = 160, N = 2048 (right plot) for x_{0} = (10, 0).

It is worth noting that it is more efficient to include the source term f in the I_{f} term of our NRBC, than to have to treat it as the right-hand side of the wave equation. In particular, for far away sources, to include the support of a nontrivial data in the computational domain, one has to choose it unnecessarily large. In the above case, for example, when Ω = Ω_{2}, to obtain a reasonable accurate solution we had to choose a very fine domain triangulation. Instead, for the discretization of Ω_{1}, a triangulation with n_{t} = 406 triangles has been more than sufficient to get a similar accuracy.

Example 3. When a scatterer is made up of several disjoint obstacles, and we are interested in the behavior of the solution only in a neighborhood of each obstacle, the use of a single artificial boundary to enclose the entire scattering region could be too expensive. In general, it is preferable to enclose each subscatterer by a separate artificial boundary. Then, the exact NRBC is defined on B = ⋃ _{j} B_{j}, where B_{j} denotes the artificial boundary that surrounds a single computational subdomain Ω_{j}. The global finite computational domain Ω = ⋃ _{j} Ω_{j} is therefore a nonconvex domain and the waves that propagate are both incoming and outgoing in each subdomain (see Figure 5).

In this example, we consider the case of two disjoint scatterers: the first one is a disk, whose boundary Γ_{1} is the circle centered at the origin and with radius r_{1} = 2, while the second one is a disk whose boundary Γ_{2} is the circle centered at C = (14, 0) and with radius r_{2} = 3. We choose the artificial boundary B_{1} as the circle centered at the origin and with radius R_{1} = 8, and the artificial boundary B_{2} as the circle centered at C and with radius R_{2} = 5. We prescribe homogeneous

Figure 5. Multiple scattering: a single artificial boundary (left) and separate artificial boundaries (right)

Figure 6. Example 3. Snapshots of the solution at different times.

boundary conditions on Γ_{1} and on Γ_{2}. We consider null source and initial velocity, while the initial configuration is u_{0}(x) = e^{−5((x_{1} −5)^{2} +x_{2}^{2}} .

Although u_{0} does not have a local support (and thus contradicts one of our assumptions), it decays exponentially fast away from its center x = (5, 0), in such a way that from the computational point of view it can be regarded as supported in a disk with radius smaller than 3 (at distance 2.7 from its center, it assumes approximately values of the order 10^{−16} ). The support of u_{0} is therefore included in Ω_{1} (see Figure 6, first plot). The disks bounded by Γ_{1} and Γ_{2} represent two soft obstacles that act as a reflecting body. With these choices, the data compatibility conditions are satisfied and the Crank–Nicolson scheme can be applied. In Figure 6 we show some snapshots of the solution obtained with our exact NRBC at different time steps, by using nt = 41458 and nt = 10782 triangles in Ω_{1} and in Ω_{2}, respectively, and a uniform partition of the time interval [0, 20] into N = 256 subintervals. No (significant) spurious reflections are produced.

4.2. Examples of frequency domain problems. Example 4. We consider here the Helmholtz problem (3.1), Ω^{e} being the region exterior to the unit disk. On the unit circle we prescribe the Dirichlet condition

(4.1)  û(x) = \frac{ı}{4}H_{0}^{(1)}(k|x − x_{0}|), x ∈ Γ

where H_{0}^{(1)} is the Hankel function of the first kind of order zero. In this case, the exact solution is known and is the field produced by the point source at x_{0}, whose expression is given by (4.1) for x ∈ ℝ^{2}. In the following numerical tests, we study

Figure 7. Example 4. Exact and approx solutions on B: k = 1e − 02, M_{ρ} = 256 and M_{θ} = 256

Figure 8. Example 4. L^{2} error behavior w.r to mesh refinement (left plot) and ratio Cond(ABC Lubich)/Cond(Bayliss–Turkel) (right plot) for k = 1e − 02 Figure 9. Example 4. Exact and approx solutions on B: k = 1, M_{ρ} = 256 and M_{θ} = 256

Figure 10. Example 4. L^{2} error behavior w.r to mesh refinement (left plot) and ratio Cond(ABC Lubich)/Cond(Bayliss–Turkel) (right plot) for k = 1

the behavior of the numerical solution and the error for some values of the wave number k.

We have applied both the finite difference scheme and the finite element method in space. Since the two methods have produced very similar results, we show here only the first one. To fix the ideas, we choose the artificial boundary as the circle of radius R = 2. By using the polar coordinates x = ρ cos(θ), y = ρ sin(θ), ρ ∈ [r, R], θ ∈ [0, 2π], we transform the bounded computational domain into the rectangle [r, R] × [0, 2π]. The Helmholtz equation in polar coordinates reads

( û_{ρρ} + \frac{1}{ρ^{2}} û_{θθ} + \frac{1}{ρ}û_{ρ}) + k^{2} û = 0. (ρ, θ) ∈ [r, R] × [0, 2π].

Figure 11. Example 4. Exact and approx solutions on B: k = 10, M_{ρ} = 256 and M_{θ} = 256

Figure 12. Example 4. L^{2} error behavior w.r to mesh refinement (left plot) and ratio Cond(ABC Lubich)/Cond(Bayliss–Turkel) (right plot) for k = 10

Moreover, we compare our exact NRBC with the second order approximate Bayliss– Turkel NRBC, which is given by

û_{ρρ} + ( \frac{3}{R} − 2ık) û_{ρ} + ( \frac{5}{4R^{2}} − \frac{3ık}{R}− k^{2}) û = 0.

We consider a discretization of the transformed domain into rectangular cells, obtained by a uniform partition of [r, R] and [0, 2π] into M_{ρ} and M_{θ} subintervals, respectively. Then, we apply the classical second order finite difference scheme for the Laplace operator. The final scheme is obtained by coupling the latter one with the collocation discretization of our NRBC. In the following numerical tests we set x_{0} = (0, 0.5), and we consider different values of k: k = 1E − 02, 1, 10 and 1E + 02.

In Figures 7, 9 and 11 we compare the real and imaginary parts of the numerical solution at the artificial boundary B with the exact ones, for k = 1E − 02, 1, 10,

Figure 13. Example 4. L^{2} error behavior w.r to mesh refinement (left plot) and ratio Cond(ABC Lubich)/Cond(Bayliss–Turkel) (right plot) for k = 100

Figure 14. Example 4. Ratio CPU(ABC Lubich)/CPU(Bayliss– Turkel) required to solve the problem by using the exact NRBC and the second order Bayliss–Turkel ABC

respectively. We omit the graph corresponding to k = 1E + 02 because of the extremely high oscillatory behavior of the solution. In Figures 8, 10, 12 and 13, on the left plots, we show the behavior of the L^{2}-norm error with respect to the mesh parameter M = M_{ρ} = M_{θ}, for all the values of the wave numbers k = 1E −02, 1, 10 and 1E + 02. In the same figure, on the right plots, we compare the condition numbers (Cond(·)) of the matrix associated with the complete discretized scheme for the two approaches. In particular, we show the ratio between the condition number when our NRBC is considered, and the condition number when the second order Bayliss–Turkel NRBC is applied. Finally, in Figure 14, we plot the ratio between the CPU times required to solve the problem by using the exact NRBC and the second order Bayliss–Turkel NRBC (note that this quantity is independent of the value of k).

References

1. M. Abramowitz, I. Stegun, (eds.) Handbook of Mathematical Functions, National Bureau of Standards, Washington D.C., 1966.

2. S. Amini, S. M. Kirkup, Solution of Helmholtz equation in the exterior domain by elementary boundary integral methods, J. Comput. Phys. 118 (1995) 208–221.

3. A. Bayliss, E. Turkel, Radiation boundary conditions for wave-like equations, Comm. Pure Appl. Math. 33 (1980) 707–725.

4. B. Engquist, A. Majda, Absorbing boundary conditions for numerical simulation of waves, Math. Comput. 31 (1977) 629–651.

5. S. Falletta, G. Monegato, L. Scuderi, Space-time BIE methods for non homogeneous exterior wave equation problems. The Dirichlet case, IMA J. Num. Anal. 32 (2012) 202–226.

6. , Space-time BIE methods for non homogeneous exterior wave equation problems. The Neumann case, IMA J. Num. Anal. 34 (2014)390–434.

7. S. Falletta, G. Monegato, An exact non reflecting boundary condition for 2D time-dependent wave equation problems, Wave Motion 51 (2014) 168–192.

8. D. Givoli, Numerical Methods for Problems in Infinite Domains, Elsevier, Amsterdam, 1992.

9. , High-order local non-reflecting boundary conditions: a review, Wave Motion 39 (2004) 319–326.

10. C. Johnson, J. C. Nedelec, On the coupling of boundary integral and finite element methods, Math. Comp. 35 (1980) 1063–1079.

11. 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.

12. R. Li, On the coupling of BEM and FEM for exterior problems for the Helmholtz equation, Math. Comp. 68 (1999) 945–953.

13. Ch. Lubich, Convolution quadrature and discretized operational calculus. I, Numer. Math. 52 (1988) 129–145.

14. , On the multistep time discretization of linear initial-boundary value problems and their boundary integral equations, Numer. Math. 67 (1994) 365–389.

15. M. Masmoudi, Numerical solution for exterior problems, Numer. Math. 51 (1987) 87–101.

16. G. Monegato, L. Scuderi, M. P. Stanić, Lubich convolution quadratures and their applications to space-time BIEs, Numer. Algorithms 56 (2010) 405–436.

17. G. Monegato, L. Scuderi, Potential evaluation in space-time BIE formulations of 2D wave equation problems, J. Comp. Appl. Math 243 (2013), 60–79.

18. R. Rangogni, M. Reali, The coupling of the finite difference method and the boundary element method, Appl. Math. Modelling 6 (1992) 233–236.

19. S. Sauter, Ch. Schwab, Boundary Element Methods, Springer, Heidelberg, 2010.

20. J. Song, L. Kai-Tai, The coupling of finite element method and boundary element method for two-dimensional Helmholtz equation in an exterior domain, J. Comput. Math. 5 (1987) 21–37.

Department of Mathematical Sciences Politecnico di Torino Torino Italy

 silvia.falletta@polito.it
 giovanni.monegato@polito.it