 # Numerical solution of a one-dimensional nonlocal Helmholtz equation by Perfectly Matched Layers

We consider the computation of a nonlocal Helmholtz equation by using Perfectly Matched Layer (PML). We first derive the nonlocal PML equation by extending PML modifications from the local operator to the nonlocal operator of integral form. We then give stability estimates of some weighted average value of the nonlocal Helmholtz solution and prove that (i) the weighted average value of the nonlocal PML solution decays exponentially in PML layers in one case; (ii) in the other case, the weighted average value of the nonlocal Helmholtz solution itself decays exponentially outside some domain. Particularly for a typical kernel function γ_1(s)=1/2 e^-| s|, we obtain the Green's function of the nonlocal Helmholtz equation, and use the Green's function to further prove that (i) the nonlocal PML solution decays exponentially in PML layers in one case; (ii) in the other case, the nonlocal Helmholtz solution itself decays exponentially outside some domain. Based on our theoretical analysis, the truncated nonlocal problems are discussed and an asymptotic compatibility scheme is also introduced to solve the resulting truncated problems. Finally, numerical examples are provided to verify the effectiveness and validation of our nonlocal PML strategy and theoretical findings.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The development of nonlocal models has grown impressively over the last decade because of its huge potential of emerging applications in various research areas, such as the peridynamical theory of continuum mechanics, the nonlocal wave propagation, and the modeling of nonlocal diffusion process [7, 15, 29, 36, 41]. In this paper, we consider the computation of a nonlocal Helmholtz equation on the whole real axis

 Lδu(x)−k2u(x)=f(x),x∈R, (1.1)

where is a constant related to the traditional wavenumber for the local Helmholtz equation, the source is supported on , and the nonlocal operator is defined as

 Lδu(x)=∫R(u(x)−u(y))γδ(y−x)dy. (1.2)

The kernel function in (1.2) is determined by a rescaling of a parent kernel through

 γδ(s)=1δ3γ1(sδ), (1.3)

where the parameter represents the range/radius of nonlocal interaction, and is piecewisely smooth, and satisfies:

• nonnegativeness: ;

• symmetry in : ;

• finite horizon: , such that if ;

• the second moment condition

Recently, much works are carried out for the simulation of nonlocal problems with free or fixed boundary conditions. There are applications in which the simulation of an infinite medium may be useful, such as wave or crack propagation in whole space. The nonlocal Helmholtz equation can be used to describe the nonlocal wave propagation. In fact, it can be derived from the nonlocal wave equation

 (∂2t+Lδ)u(x,t)=f(x,t),x∈R, (1.4)

where represents the displacement field, is the space-time source term with compact support at all time. If we make the ansatz that is a superposition of the time-harmonic sources . Then, for each , the corresponding mode satisfies

 Lδ(u(x)e−ikt)−k2u(x)e−ikt=f(x)e−ikt.

Thus the time domains solution is the sum of the time-harmonic modes over all possible values of .

The aim of the paper is to develop an efficient computation of a nonlocal Helmholtz equation on the whole real axis. Absorbing boundary conditions (ABCs) are a successful approach to simulate the wave behaviors of a physical domain of interest by imposing a suitable boundary condition, to absorb the impinging wave at artificial boundaries. For the construction of tractional ABCs, it is well studied for local problems [2, 20, 19, 18, 30, 21, 3], and there are also much progress for nonlocal problems, see [28, 39, 40, 38]. In this paper, we consider to construct the perfectly matched layer (PML) as ABCs for a 1D nonlocal Helmholtz equation. The PML, originally proposed by Berenger , has the two important features: (i) the wave in a special designed layer decays exponentially, and (ii), if the wave reflects off the truncated boundary, the returning waves after one round trip through the absorbing layer are very tiny [1, 9, 6, 13, 34, 24, 11, 10, 4, 25, 8, 23]. Specifically for peridynamics, Wildman and Gazonas  present a PML by treating the nonlocal kernel as the convolution of the displacement with the second derivative of a nascent Dirac delta distribution. In this paper, we derive a PML different from  and has a simpler structure, and more importantly, give the stability estimate theoretically.

The contributions of this paper are given as follows:

• We extend the strategy of the PML method from the local operator to nonlocal operator, and obtain a nonlocal PML equation directly from the weak form of nonlocal Helmholtz equation. The resulting nonlocal PML equation converges to the corresponding local PML equation while the nonlocal interaction horizon vanishes. Such consistency is useful to demonstrate the validation/verification of our analytical continuation for the nonlocal operator.

• The properties of the nonlocal Helmholtz solution and nonlocal PML solution are analyzed. To do so, we introduce the weighted average value with respect to the nonlocal Helmholtz solution , defined by

 uw(x):=∫Ru(t+x)κ(t)dt,

where the weight is given in (3.6)-(3.7). For general kernel function , the stability estimates of and its analytical continuation are established. Specifically, we prove in Theorem 1: (1) the weighted average value of the nonlocal PML solution exponentially decays in PML layers under a quantitative condition depending on the wavenumber (taken smaller value) and the kernel ; (2) under another quantitative condition, the weighted average value itself of the nonlocal Helmholtz solution exponentially decays outside of a domain large enough. The result (1) shows the PML is efficient to absorb the waves when they impinge the PML layers. The result (2) suggests that one can directly truncate the nonlocal Helmholtz equation by putting forced boundary layers (i.e. homogeneous Dirichlet boundary constrains) at sufficiently large instead of using PML.

• In particular, the refined estimates of the nonlocal Helmholtz solution and its PML solution are further established for a typical kernel . To this end, we first derive the exact formula of the Green’s function for nonlocal Helmholtz equation (1.1), and then present the stability estimates of and instead of their weighted average values. Specifically, we prove in Theorem 3 that: (1) if , the nonlocal PML solution exponentially decays in PML layers; (2) if , the nonlocal Helmholtz solution itself exponentially decays outside . The resulting Green’s function for this typical kernel shows an attractive property that the local and nonlocal Helmholtz solutions have an intimate connection, i.e., they can be expressed by each other.

The behaviors of nonlocal Helmholtz (PML) solutions and Green’s function for the Helmholtz equation with are newly obtained, although the Green’s functions of some nonlocal problems have been previously studied [35, 22, 36, 27, 37]

which are mainly based on Fourier transforms and Laplace transforms, but can not be extended to the nonlocal Helmholtz equation.

On the other hand, the nonlocal operator has an intimate connection with some local differential operator. More specifically, as the nonlocal horizon vanishes, see [15, 41, 14, 17], the nonlocal operator converges to the second-order differential operator

 limδ→0Lδu(x)=−∂2xu(x):=L0u(x) (1.5)

under the above second moment condition. Such consistency is quite useful not only for the modeling, but also for the validation/verification of numerical simulations. For discrete schemes, it is useful to make use of the asymptotic compatibility (AC) scheme, a concept developed in [31, 32] and further extended in [16, 33], to discretize the nonlocal operator to preserve such (analogous) limit (1.5) in a discrete level. Recently, Du et al. [17, 28] proposed an AC scheme to discretize the nonlocal operator in unbounded multi-scale mediums. In this paper we extend the method  to discretize the nonlocal operator with PML modifications.

The outline of this paper is as follows. In Section 2 we propose our PML technique for the nonlocal Helmholtz equation by extending complex coordinate transforms for the local Helmholtz equation. Section 3 is devoted to analyzing the nonlocal Helmholtz solution and give the truncated nonlocal problems based on our theoretical analysis. In particular, for a special kernel we give the Green’s function of the nonlocal problem. In Section 4, we introduce the AC scheme for numerically solving the truncated problems. In Section 5, some numerical tests are provided to verify the theoretical results and the effectiveness of our PML.

## 2 PML for the local and nonlocal Helmholtz equations

In this section we first review the classical PML method for local Helmholtz equations, and recall some properties of local PML solutions. We then employ the weak form of nonlocal equation (1.1) to derive the corresponding nonlocal PML equation, and finally consider the local limits between the local and nonlocal PML equations as the nonlocal horizon vanishes.

### 2.1 PML for the local Helmholtz equation

We now recall the PML method for a local Helmholtz equation

 L0uloc(x)−k2uloc(x)=f(x) (2.1)

with Sommerfeld radiation condition. As discussed in [9, 12, 8, 13, 24], the PML modifications can be viewed as a complex coordinate stretching of the original scattering problem by constructing an analytic continuation to the complex plane, namely,

 ~x(x)=∫x0ω(t)dt=∫x01+iσ(t)dt, (2.2)

where , called the absorption function, is equal to in and positive outside . is the imaginary unit.

Setting , and using the relation

 ∂x→1ω(x)∂x=11+iσ(x)∂x,

the modified modal solution satisfies the governing equation

 ~L0uloc(x)−k2ω(x)uloc(x)=f(x) (2.3)

with the local PML operator given by

 ~L0uloc(x)=−∂x(1ω(x)∂x~uloc(x)).

The PML method for local Helmholtz equations has been well studied, see [25, 26]. Here we recall some properties of solutions as follows:

• The local Helmholtz solution is oscillating and possibly doesn’t vanish as .

• The local PML solution is an analytic continuation of in the complex coordinate. The solution is unique in the sense that over the domain .

• The stability estimate of the solution to the local Helmholtz equation is given as

 ∥uloc∥H1(Ω)+k∥uloc∥L2(Ω)≤C∥f∥L2(Ω).
• The analytical continuation changes oscillating waves into exponentially decaying waves outside the region of interest,

 |~u(x)|≤Ce−k|∫x0σ(t)dt|∥f∥L2(Ω)for |x|>l. (2.4)

### 2.2 PML for the nonlocal Helmholtz equation

The original PML technique is proposed for the standard PDEs, such as the local Helmholtz equation shown above and the Maxwell’s equation. We point out that PML is for the problems with some far field boundary conditions. An typical example is the Sommerfeld radiation condition for the “local” Helmholtz equation in a homogeneous medium. However, it is unknown for the nonlocal Helmholtz equation (1.1). Therefore, in this paper we simply assume that a suitable boundary condition at infinity is imposed to exclude energy incoming from infinity and only to allow energy outgoing to infinity. Based on this assumption, we here manage to extend the PML to solve the nonlocal Helmholtz equation.

To do so, we introduce the PML by considering the weak form of Eq. (1.1)

 (Lδu,v)−k2(u,v)=(f,v),∀v∈C∞0(R), (2.5)

where denotes the inner product in the complex valued -space, and

 (Lδu,v)=12∫R∫R(u(x)−u(y))(¯v(x)−¯v(y))γδ(y−x)dydx. (2.6)

Here represents the complex conjugate of . We now apply the same transform as (2.2) and immediately produce the corresponding differential forms by

 x→~x=∫x0ω(t)dt, y→~y=∫y0ω(t)dt, (2.7) dx→ω(x)dx, dy→ω(y)dy. (2.8)

Setting , , and , we can transform (2.5) into the following nonlocal equation with PML modifications

 12∫R∫R(~u(x)−~u(y))(~¯v(x)−~¯v(y))γδ(~y−~x)ω(x)ω(y)dydx (2.9) −k2∫R~u(x)~¯v(x)ω(x)dx=∫Rf(x)~¯v(x)ω(x)dx.

Using the facts and , we have the strong form of (2.9) as

 ~Lδ~u(x)−k2ω(x)~u(x)=f(x), (2.10)

where the nonlocal operator with PML modifications is given by

 ~Lδ~u(x)=∫R(~u(x)−~u(y))γδ(~y−~x)ω(x)ω(y)dy. (2.11)

Thus, we have that the solution of Eq. (2.10) is an analytic continuation of the solution of Eq. (1.1) in the complex coordinate, and it holds that for . We emphasize that the kernel function in the nonlocal PML operator (2.11) must be the analytic continuation of the original kernel in complex coordinates. From Figure 2 in the section of numerical examples, one can see that such analytic continuation of the kernel plays an important role to get the correct solution in simulations.

On the other hand, it is well known that the nonlocal operator has an intimate connection with the local differential operator shown in Eq. (1.5). It is interesting to ask if the nonlocal PML equation (2.10) converges to its corresponding local PML equation (2.3) as the nonlocal interaction horizon vanishes. Such consistency is useful to demonstrate the validation/verification of our analytical continuation for the nonlocal operator. In fact, by using the second moment condition of the kernel function and Taylor expansion, it is straightforward to verify that

 limδ→0+(~Lδ~u(x),v(x))=(~L0~u(x),v(x)). (2.12)

Thus, the local limit (2.12) implies that the nonlocal PML equation (2.10) will converge to the corresponding local PML equation (2.3) as the nonlocal interaction horizon tends to zero.

## 3 The exponentially decaying waves

In the previous section, the “a priori” estimates of solutions are introduced for the local Helmholtz equation and its PML equation, see [25, 26]. These properties are useful for the further theoretical analysis of the PML method. Here we also want to study the solution properties of the nonlocal Helmholtz equation (1.1) and the nonlocal PML equation (2.9). Generally, the properties of nonlocal solutions are not trivial to be explored since it is hard to exactly express the Green’s function for the general kernel. To do so, we first introduce the following function

 Gx0(x)=⎧⎪⎨⎪⎩C1(x0)e−i~kx,x≤x0withC1(x0)=−ei~kx02i~k,C2(x0)ei~kx,x>x0withC2(x0)=−e−i~kx02i~k, (3.1)

where is the solution of the following equation

 μ(~k) :=1δ2∫R(1−ei~kδs)γ1(s)ds=k2. (3.2)

The identity (3.2) implies that is only dependent of the kernel , and the wavenumber , and converges to as under the assumption of second moment condition of the kernel .

Note that it also holds since the kernel is symmetric, which implies that there exists in the set of solutions to Eq. (3.2) such that or if . Specifically, there exists a positive value which is

 k0=sup~k∈R+√μ(~k),

such that if , and if . In this paper, we only need to take the root of identity (3.2) as a positive number or a complex number with positive imaginary part.

To make clear the relation between and , we here show two examples.

• Example 1. Taking the kernel , we have

 μ(~k)=~k2/(δ2~k2+1).

The direct calculation shows that one typical solution of Eq. (3.2) is

 ~k=k√11−(δk)2. (3.3)

It’s clear that if and if .

• Example 2. Taking the kernel , we have

 μ(~k)=4δ2(1−e−(δ~k)2/4).

One typical solution of Eq. (3.2) is

 ~k=2δ√−log(1−(kδ)2/4). (3.4)

It is clear that if , and if .

To explore the properties of the solution to the nonlocal Helmholtz equation (1.1), we now introduce the definition of a weighted average value with respect to by

 uw(x):=∫Ru(t+x)κ(t)dt, (3.5)

where the weight is given by for and for with

 w1(t)=12δ2i~k∫tδ−∞(ei~k(t−δs)−e−i~k(t−δs))γ1(s)ds=1δ2~k∫tδ−∞sin~k(t−δs)γ1(s)ds, (3.6) w2(t)=12δ2i~k∫+∞tδ(e−i~k(t−δs)−ei~k(t−δs))γ1(s)ds=−1δ2~k∫+∞tδsin~k(t−δs)γ1(s)ds. (3.7)

The weighted average value satisfies the following identity.

###### Lemma 1

The weighted average value defined by (3.5) satisfies

 uw(x)=∫RGx(y)f(y)dy. (3.8)

For brevity, we leave the proof of Lemma 1 in the appendix 1.

###### Theorem 1

Assume that is a constant only depending on , and is the root of relation (3.2). Then the weighted average value of the solution to the nonlocal Helmholtz equation (1.1), and its analytic continuation by complex coordinate transform (2.7) have the following estimates:

• If is positive real and , it holds that

 |uw(x)|H1(Ω)+~k∥uw(x)∥L2(Ω)≤C∥f∥L2(Ω), (3.9) |~uw(x)|≤Ce−~k|∫x0σ(t)dt|∥f∥L2(Ω),for |x|>l. (3.10)
• If is a complex number with and , it holds that

 |uw(x)|H1(Ω)+|~k|∥uw(x)∥L2(Ω)≤C∥f∥L2(Ω), (3.11) |uw(x)|≤C1|~k|e−I(~k)(|x|−l)∥f∥L2(Ω),for |x|>l. (3.12)

Proof. (1) If is positive, (3.1) is actually the Green’s function of the following local Helmholtz equation 

 −∂2x¯u(x)−~k2¯u(x)= f(x) x∈R, ∣∣∂ruw−i~kuw∣∣→ 0 as r=|x|→+∞.

From Eq. (3.8), is the solution of the Helmholtz problem above. By using the results for local problems in [25, 26] (see Eq. (2.4)), we can directly prove the stability estimate (3.9) for and the exponentially decaying estimate (3.10) for .

(2) From Eq. (3.1) and (3.8), we get

 ∥uw∥2L2(Ω)= ∫Ω∣∣∫ΩGx(y)f(y)dy∣∣2dx≤∥f∥2L2(Ω)∫Ω∫Ω∣∣Gx(y)∣∣2dydx ≤ 14|~k|2∥f∥2L2(Ω)(∫Ωe−2I(~k)xdx∫x−le2I(~k)ydy+∫Ωe2I(~k)xdx∫lxe−2I(~k)ydy) = 14|~k|2∥f∥2L2(Ω)⋅(2lI(~k)+12I(~k)2(e−4I(~k)l−1)) = 2l2|~k|2⋅4I(~k)l+e−4I(~k)l−1(4I(~k)l)2∥f∥2L2(Ω) ≤ l2|~k|2∥f∥2L2(Ω),

where we have used the fact that .

By similar arguments, we have

 |uw|2H1(Ω)≤l2∥f∥2L2(Ω).

We now prove that is exponentially decaying as . For , it holds

 |uw(x)|= ∣∣∣∫ΩGx(y)f(y)dy∣∣∣=∣∣∣∫ΩC0ei~k(x−y)f(y)dy∣∣∣ ≤ 12|~k|(∫Ωe−2I(~k)(x−y)dy)12∥f∥L2(Ω) ≤ √2l2|~k|e−I(~k)(x−l)∥f∥L2(Ω).

Similarly, we can get for . This completes the proof.

###### Remark 1

For the main results in Theorem 1, it is interesting to point out that

• For the case of , the estimate (3.10) shows that the weighted average value of the solution to the nonlocal PML equation (2.10) decays exponentially as .

• These results are still valid for the case , i.e., . This means that decays as . In this situation, we can truncate the kernel when for (a finite constant) for a given tolerance . Taking the weighted average value with a truncated kernel by

 ^uwδ(x):=∫δ^lγ−δ^lγu(t+x)κ(t)dt,

we then have

 |~uwδ(x)−^uwδ(~x)|≤2ϵ∥~u(x)∥L1(R)/(~k2δ3).

This implies we can use the truncated kernel to replace the original one with a tolerance error.

### 3.1 A refined estimate for a typical exponential kernel γ1(s)=12e−|s|

In the previous subsection, we considered the general compactly supported kernel function and have proved that the weighted average value of the nonlocal PML solution or the nonlocal Helmholtz solution decays exponentially as (see Theorem 1). Here we consider a special kernel , and directly analyze the behavior of or as instead of their weighted average values. To do so, we first consider the Green’s function for the following equation

 LδGLδx0(x)−k2GLδx0(x)=D(x−x0), (3.13)

where we use to represent the Dirac delta function differing from the nonlocal horizon . The following theorem will present the explicit formula of .

###### Theorem 2

For the kernel function , the Green’s function , the solution of Eq. (3.13), has the following exact expression

 GLδx0(x)=1(1−(δk)2)2Gx0(x)+δ21−(δk)2D(x−x0), (3.14)

where is given by (3.1) with .

Proof. By definition of (3.5) and a simple calculation, at point we get

 uw(x0)= 12δ(1+~k2δ2)∫Ru(t+x0)e−|t|δdt = 1−(δk)22δ∫Ru(t+x0)e−|t|δdt, (3.15)

where is the solution to Eq. (1.1), and in the last identity we use the fact , which can be directly calculated from the relation . On the other hand, at point also satisfies

 12δ3∫Ru(t+x0)e−|t|δdt =12δ3∫Ru(x0)e−|t|δdt−k2u(x0)−f(x0) =1δ2[(1−δ2k2)u(x0)−δ2f(x0)]. (3.16)

From (3.16), we obtain

 u(x0)= 12δ(1−(δk)2)∫Ru(t+x0)e−|t|δdt+δ21−(δk)2f(x0) = 1(1−(δk)2)2uw(x0)+δ21−(δk)2f(x0) = 1(1−(δk)2)2∫RGx0(x)f(x)dx+δ21−(δk)2f(x0) = ∫R(1(1−(δk)2)2Gx0(x)+δ21−(δk)2D(x−x0))f(x)dx, (3.17)

where the second identity uses (3.15), the third identity uses (3.8), and the last identity uses the definition of Dirac delta function. Finally from the identity (3.17), we have that the Green’s function satisfies (3.14). The proof is completed.

###### Remark 2

(a) Since the far field boundary condition for the nonlocal Helmholtz equation is still an open problem, we verify that Green’s function is unique only by the property, and when the nonlocal interaction vanishes, i.e., , the Green’s function converges to the Green’s function of the local Helmholtz equation with Sommerfeld radiation condition

 −∂2xgx0(x)−k2gx0(x) =D(x−x0) x∈R, ∣∣∂rgx0(x)−ikgx0(x)∣∣ →0 as r=|x|→+∞.

It is compatible with the asymptotic convergence property of nonlocal problems and can be proved simply by using the fact that as .

(b) We point out that the result in Theorem 2 is attractive for the study of local Helmholtz equations. For a given positive wavenumber , the solution to the following Helmholtz problem

 −∂2xuloc(x)−~k2uloc(x)= f(x),x∈R, (3.18) ∂ruloc(x)−i~kuloc(x)= o(1),as r=|x|→+∞, (3.19)

can be obtained by solving the nonlocal Helmholtz equation

 Lδu(x)−k2u(x)=f(x), (3.20)

where obtained by (3.3). This is to say, from Eq. (3.17) we have

 (3.21)

In light of Theorem 1, we now present the estimates of nonlocal solutions and .

###### Theorem 3

Assume that is a constant only depending on . Then for the solution of the nonlocal Helmholtz equation (1.1) and its analytic continuation by complex coordinate transform (2.7), we have the following estimates:

• If and , it holds that

 ∥u∥L2(Ω)≤1k(1−(δk)2)(C(1−(δk)2)12+δ2k)∥f∥L2(Ω), (3.22) |u|H1(Ω)≤C(1−(δk)2)2∥f∥L2(Ω)+δ21−(δk)2|f|H1(Ω), (3.23) |~u(x)|≤C1(1−(δk)2)2e−k√1−(δk)2|∫x0σ(t)dt|∥f∥L2(Ω),for |x|>l. (3.24)
• If and , it holds that

 ∥u∥L2(Ω)≤1k((δk)2−1)⎛⎜⎝C((δk)2−1)12+δ2k⎞⎟⎠∥f∥L2(Ω), (3.25) |u|H1(Ω)≤C(1−(δk)2)2∥f∥L2(Ω)+δ2(δk)2−1|f|H1(Ω), (3.26) |u(x)|≤Ck((δk)2−1)32e−k√(δk)2−1(|x|−l)∥f∥L2(Ω),for |x|>l. (3.27)

Proof. (1) From (3.17), can be expressed by

 u(x)=1(1−(δk)2)2uw(x)+δ21−(δk)2f(x). (3.28)

Since , we have if . Combining (