Numerical solution of scattering problems using a Riemann--Hilbert formulation

TL;DR

This paper introduces a fast, spectrally accurate numerical method for solving scalar and matrix Wiener--Hopf problems by formulating them as Riemann--Hilbert problems on the real line, applicable to diffraction problems.

Contribution

It develops a novel numerical approach leveraging Riemann--Hilbert formulation for Wiener--Hopf problems, enabling high-accuracy solutions for generalized diffraction scenarios.

Findings

Achieved spectrally accurate numerical solutions.

Successfully solved generalized Wiener--Hopf problems.

Demonstrated applicability to classical diffraction problems.

Abstract

A fast and accurate numerical method for the solution of scalar and matrix Wiener--Hopf problems is presented. The Wiener--Hopf problems are formulated as Riemann--Hilbert problems on the real line, and a numerical approach developed for these problems is used. It is shown that the known far-field behaviour of the solutions can be exploited to construct numerical schemes providing spectrally accurate results. A number of scalar and matrix Wiener--Hopf problems that generalize the classical Sommerfeld problem of diffraction of plane waves by a semi-infinite plane are solved using the approach.

Tables

Table 1. Table 1 : Diffraction by a semi-infinite plane and boundary conditions on top and bottom sides of the semi-infinite boundary x > 0 𝑥 0 x>0 , y = 0 𝑦 0 y=0 . Note that in some cases [ 12 , 3 , 25 ] the semi-infinite plane lies along the negative real axis. The parameter μ 𝜇 \mu is defined to be μ = ρ 2 / ρ 1 𝜇 subscript 𝜌 2 subscript 𝜌 1 \mu=\sqrt{\rho_{2}/\rho_{1}} .

(B1)	(B2)	$c_1/c_2$	${\rho }_1/{\rho }_2$	$a$	References
hard	hard	1	1	1/2	Sommerfeld [3, 24]
soft	soft	1	1	1/2	Sommerfeld [3, 24]
soft	hard	1	1	1/4	Rawlins [25]
impedance1	impedance1	1	1	1/2	Senior [3, 26]
impedance1	impedance2	1	1	1/2	Hurd [9], Barton & Rawlins [27]
soft	hard	$\ne 1$	$\ne 1$	${\pi }^{-1}{\tan }^{-1}\mu$	Wickham [21], Abrahams [12]

Equations

$$A(\alpha ){\Phi }_{+}(\alpha )+B(\alpha ){\Phi }_{-}(\alpha )+C(\alpha )=0,$$

$$A(\alpha ){\Phi }_{+}(\alpha )+B(\alpha ){\Phi }_{-}(\alpha )+C(\alpha )=0,\alpha \in \Gamma ,$$

$$\varphi (r,\theta )\sim D(\theta )\frac{{\mathrm{e}}^{\mathrm{i}kr}}{\sqrt{r}},\text{as}r\to \infty ,$$

$$\Phi (\alpha )=\sum _{k=0}^{n-1}{\check{U}}_k{T}_k(x)=\sum _{k=0}^{n-1}{\check{U}}_k{T}_k(M^{-1}(\alpha ))$$

$$\alpha =M(x),x=M^{-1}(\alpha ).$$

$${\check{U}}_k=\sum _{q=1}^n{F}_{kq}\Phi ({\alpha }_q),$$

$${\mathbf{x}}^{\mathbb{I}}=(x_1,\dots ,x_n)=\left(-1,\cos \frac{(n-2)\pi }{n-1},\dots ,\cos \frac{\pi }{n-1},1\right).$$

$$\mathcal{C}\Phi (\alpha )=\frac{1}{2\pi \mathrm{i}}{\int }_{\Gamma }\frac{\Phi ({\alpha }^{\prime })}{{\alpha }^{\prime }-\alpha }\mathrm{d}{\alpha }^{\prime }.$$

$$G(\alpha )=\sum _{k=0}^{n-1}{\check{U}}_k[\mathcal{C}{T}_k](M^{-1}(\alpha ))=\sum _{k=0}^{n-1}{\check{U}}_k{S}_k(M^{-1}(\alpha )),$$

$$G({\alpha }_p)=\sum _{k=0}^{n-1}{\check{U}}_k{S}_k(M^{-1}({\alpha }_p)).$$

$$G({\alpha }_p)=\sum _{k=0}^{n-1}\sum _{q=1}^n{F}_{kq}{S}_k(M^{-1}({\alpha }_p))\Phi ({\alpha }_q)=\sum _{q=1}^n{C}_{pq}\Phi ({\alpha }_q),$$

$$C_{pq}=\sum _{k=0}^{n-1}{F}_{kq}[\mathcal{C}{T}_k](M^{-1}({\alpha }_p))=\sum _{k=0}^{n-1}{F}_{kq}{S}_k(M^{-1}({\alpha }_p)).$$

$$C_{pq}=\sum _{j=1}^d{C}_{pq}^{(j)},$$

$$\alpha =M(x)=\frac{x}{1-x^2}$$

$$M_{1,2}^{-1}(\alpha )=\frac{-1\pm \sqrt{1+4{\alpha }^2}}{2\alpha }\Rightarrow M_1^{-1}(\alpha )=x,M_2^{-1}(\alpha )=-x^{-1}.$$

$$C^{+}_{1}=\frac{1}{{\rm i}\pi}\left(\begin{array}[]{ccc}\frac{\log 2}{2}+\frac{{\rm i}\pi}{2}&&\\ &-2\text{ diag}(\text{arctanh }T_{\downarrow}^{-1}({\bf x}^{\mathbb{I}}))&\\ &&-\frac{\log 2}{2}+\frac{{\rm i}\pi}{2}\end{array}\right)+\frac{1}{{\rm i}\pi}F^{-1}PF,$$

$$T_{\downarrow }^{-1}(x)=x-\mathrm{i}\sqrt{1-x}\sqrt{1+x}.$$

$$C^{+}_{2}=\left(\begin{array}[]{c}\mu^{R}\\ \Psi^{z}_{n}\\ \mu^{L}\end{array}\right)F,$$

$$z=T_{+}^{-1}(M_2^{-1}(\alpha )),\text{where}{T}_{+}^{-1}(x)=x-\sqrt{x-1}\sqrt{1+x}.$$

$$C^{+}=C_1^{+}+C_2^{+}-1_{n\times n}\text{diag}({\mu }^L+{\mu }^R)F.$$

$$\alpha \sim \frac{1}{2(1-x)}$$

$$\alpha =M(x)=\frac{x+x^3}{(1-x^2{)}^2}$$

$$M_1^{-1}(\alpha )=x,M_2^{-1}(\alpha )=x^{-1},$$

$$M_{3,4}^{-1}(\alpha )=\frac{-c(x)\pm \sqrt{c(x{)}^2-4}}{2},\text{where}c(x)=x+x^{-1}-{\alpha }^{-1}.$$

$$\alpha \sim \frac{1}{2(1-x{)}^2}$$

$$C^{+}=C_1^{+}+C_2^{+}+C_3^{+}+C_4^{+}-1_{n\times n}\text{diag}2({\mu }^L+{\mu }^R)F.$$

$$[A({\alpha }_p)C^{+}+B({\alpha }_p)C^{-}]\Phi ({\alpha }_p)+C({\alpha }_p)=0$$

$${\alpha }_p^{\text{rotated}}={\alpha }_p{\mathrm{e}}^{-\mathrm{i}\chi },$$

$$E_n^r={\left[{\int }_{-1}^1|Q(x)-Q_n(x){|}^r\mathrm{d}x\right]}^{1/r}$$

$$E_n^r\approx {\left[\sum _{q=1}^n{w}_q|Q(x_q)-Q_n(x_q){|}^r\right]}^{1/r},$$

Full text

Numerical solution of scattering problems using a Riemann–Hilbert formulation

Stefan G. Llewellyn Smith1,2 and Elena Luca1

1Department of Mechanical and Aerospace Engineering

Jacobs School of Engineering, UCSD

La Jolla, CA 92093-0411, USA.

2Scripps Institution of Oceanography, UCSD

La Jolla, CA 92039-0213, USA.

sgls@ucsd.edu

elouca@eng.ucsd.edu

Abstract

A fast and accurate numerical method for the solution of scalar and matrix Wiener–Hopf problems is presented. The Wiener–Hopf problems are formulated as Riemann–Hilbert problems on the real line, and a numerical approach developed for these problems is used. It is shown that the known far-field behaviour of the solutions can be exploited to construct numerical schemes providing spectrally accurate results. A number of scalar and matrix Wiener–Hopf problems that generalize the classical Sommerfeld problem of diffraction of plane waves by a semi-infinite plane are solved using the approach.

1 Introduction

The Wiener–Hopf (WH) method was originally developed to solve integral equations and mixed boundary value problems [1, 2]. Standard references include [3, 4, 5]. In its direct formulation introduced by Jones [6], the WH method, in combination with the Fourier transform, reduces the solution of a boundary-value problem (in e.g. electrodynamics, acoustics, elasticity, etc…) to the problem of solving a functional equation by finding functions analytic in the upper and lower complex half-planes by factorization. The typical WH functional equation (with $\alpha$ the transform variable) is

[TABLE]

valid in the strip $a_{-}<\text{Im}(\alpha)<a_{+}$, where the functions $\Phi_{+}(\alpha)$, $\Phi_{-}(\alpha)$ are, respectively, analytic in $a_{-}<\text{Im}(\alpha)$ and $\text{Im}(\alpha)<a_{+}$. Further, $A(\alpha)$, $B(\alpha)$ and $C(\alpha)$ are known functions of $\alpha$. In addition, the behaviour of $\Phi_{+}(\alpha)$, $\Phi_{-}(\alpha)$ for large $|\alpha|$ is given by the behaviour of physical variables near the origin.

While it is straightforward to apply the WH method to a scalar problem for which an explicit solution formula based on Cauchy’s integral formula exists [3], there is no general method for carrying out the decomposition for matrix functions. A survey of constructive methods for factorization problems is given by Rogosin & Mishuris [7]. The decomposition can be carried out for matrices of special form, following the ideas of Khrapkov [8], Hurd [9], Daniele [10] and Rawlins & Williams [11], but for other matrices, the only general semi-analytical approach that has been put forward is the Padé decomposition method of Abrahams [12], in which the matrices are approximated by rational functions, permitting a decomposition. The method seems to work well in applications [13, 14] but a lot of manipulation of the equations and algebra is required. A recent study by Kisil [15] proposes an iterative method for triangular matrix problems with exponential factors.

WH problems are related to Riemann–Hilbert (RH) problems. The former connect boundary values of sectionally analytic functions in a common strip of analyticity, while the latter couple these on a contour. The RH problem asks for the construction of a function $\Phi(\alpha)$ that is analytic everywhere in the complex plane except along a given oriented contour $\Gamma$ where it has a prescribed jump

[TABLE]

where $\Phi_{+}(\alpha)$ and $\Phi_{-}(\alpha)$ are the representations of $\Phi(\alpha)$ on the $+$ and $-$ sides of the contour, respectively and $A(\alpha)$, $B(\alpha)$ and $C(\alpha)$ are known functions on $\Gamma$. The functional form (2) is identical to (1), but it is now valid along a given oriented contour rather than a strip. Many RH problems arise in the context of singular integral equations, but RH problems have also been connected to random matrix theory, nonlinear special functions, nonlinear wave equations, and other problems. Olver [16] and Trogdon & Olver [17] have recently developed accurate and efficient numerical algorithms for the solution of RH problems.

The aim of this study is to present fast and accurate numerical schemes for the solution of scalar and matrix WH problems by exploiting the links between the WH and RH problems. The idea is to solve the associated RH problems, adapting the methods of Olver [16] and Trogdon & Olver [17] to take into account the known far-field behaviour of the solutions. In particular, the approach adopted is to use rational mappings with multiple inverses, extending the Möbius mappings discussed in [16, 17], to account for the $\mathcal{O}(\alpha^{-1/2})$ decay of the solutions for large $|\alpha|$. We focus here on problems of diffraction of plane waves by a semi-infinite plane that produce scalar and matrix WH problems. For the acoustic scattering problems considered, the conditions on $A(\alpha),B(\alpha),C(\alpha)$ guaranteeing existence and uniqueness of solutions are satisfied [18, 17]. These schemes are also applicable to other more complicated diffraction problems, such as the diffraction of plane waves by two identical strips [19].

In § 2 we present the numerical scheme for the solution of Riemann–Hilbert problems using rational mappings. In § 3 we give error measures that we use to verify our numerical solutions. This is followed by implementation of the numerical scheme for a number of scalar and matrix problems for diffraction by a semi-infinite plane (§ 4). In § 5 we compute the far-field directivity, $D(\theta)$, which is defined via

[TABLE]

where $\phi$ is the diffracted field, $k$ is the acoustic wave number and $(r,\theta)$ are polar coordinates. Finally we conclude and discuss further applications in § 6.

2 Numerical solution of Riemann–Hilbert problems

Following Olver [16] and Trogdon & Olver [17], we present the numerical scheme for the solution of RH problems (2). Although the method can be used to solve RH problems over any contour $\Gamma$ whose individual pieces can be mapped from the unit interval $\mathbb{I}=[-1,1]$, here we take $\Gamma=\mathbb{R}$, since this is relevant to WH problems.

The aim is to find a function $\Phi(\alpha)$ analytic everywhere in the complex plane except along $\mathbb{R}$ where it has a prescribed jump of the form (2). We represent the function $\Phi(\alpha)$ by

[TABLE]

for $\alpha\in\mathbb{R}$ and by analytic continuation off $\mathbb{R}$. Here $T_{k}(x)$ is the usual $k$th Chebyshev polynomial and $M:\mathbb{I}=[-1,1]\mapsto\mathbb{R}$ is a mapping

[TABLE]

In this study, we shall consider rational mappings $M$ with multiple inverses, extending the Möbius mappings discussed in [16, 17]. (The notation $M^{-1}$ refers to mapping back to $\mathbb{I}$; subscripts for $M^{-1}$ will indicate other branches.) The function $\Phi(\alpha)$ can be scalar, vector or matrix, and the coefficients $\check{U}_{k}$ are of the same nature. The coefficients $\check{U}_{k}$ can be found from the function values via

[TABLE]

where $F_{kq}$ are the elements of the Chebyshev operator $\mathcal{F}$. We seek a collocation method in which the function $\Phi$ and the relation (2) are evaluated at points $\alpha_{q}=M(x_{q})$ along the contour $\Gamma=\mathbb{R}$, where $\{x_{q}|q=1,\dots,n\}$ are the Chebyshev points of the second kind:

[TABLE]

2.1 Collocation method

The Cauchy operator is defined by

[TABLE]

The Cauchy operators $\mathcal{C}^{\pm}$ are defined as the limiting values of the Cauchy operator (8) as $\alpha$ approaches the oriented contour $\Gamma=\mathbb{R}$ from the $+$ and $-$ sides, i.e. $\mathcal{C}^{+}\Phi=\Phi_{+}$ and $\mathcal{C}^{-}\Phi=\Phi_{-}$. To construct a collocation method for solving (2), we must compute the Cauchy matrices $C^{\pm}$ at points $\alpha_{q}=M(x_{q})$, $x_{q}\in{\bf x}^{\mathbb{I}}$, along $\Gamma$. (The notation $\mathcal{C}^{\pm}$ refer to the operators; $C^{\pm}$ are matrices.) Since $C^{+}-C^{-}=I$, it is sufficient to construct $C^{+}$.

Consider

[TABLE]

where $S_{k}\equiv{\cal C}T_{k}$. This is the Cauchy transform of $\Phi(\alpha)$, ${\cal C}\Phi(\alpha)$, up to a constant which requires analysis for large $|\alpha|$. Next, we evaluate $G(\alpha)$ at the point $\alpha_{p}=M(x_{p})$, where $x_{p}$ is a Chebyshev point in (7), giving

[TABLE]

Substituting in for the coefficients $\check{U}_{k}$ (6), we can write this as a linear operator acting on the function values:

[TABLE]

where we have defined

[TABLE]

In the next section, we consider rational mappings with $d$ inverses and, therefore, we must write

[TABLE]

where $C^{(j)}_{pq}$ is given by (12) with $M^{-1}$ replaced by the $j$-th inverse $M_{j}^{-1}$. When the contour $\Gamma$ is unbounded with the endpoints $\pm 1$ of ${\bf x}^{\mathbb{I}}$ the preimages of $|\alpha|=\infty$, the above form of ${\cal C}\Phi(\alpha)$ needs to take into account the limiting behaviour for large $|\alpha|$ [16]; this will be discussed in the next section.

2.2 Rational mappings and Cauchy matrices

Olver [16] and Trogdon & Olver [17] considered mappings $M$ that are Möbius transformations and also discussed the extension to polynomials of degree $d$ with $d$ inverses. In what follows, we extend their formulation to rational functions $M$. This allows us to represent functions that have $\mathcal{O}(\alpha^{-1/2})$ decay for large $|\alpha|$. The use of the Chebyshev expansion (4) and subsequent matrix operators imply spectral convergence of the numerical scheme [20].

2.2.1 2-to-1 mapping

We consider the rational mapping

[TABLE]

which has two inverses:

[TABLE]

The points $x=\pm 1$ are mapped to $\pm\infty$ and, therefore, the images of $M_{j}^{-1}(\alpha)$, $j=1,2$ are both $\mathbb{R}$. A schematic is given in figure 1. The rational mapping (14), which is an odd function of $x$, has the desirable property that it maps the unit $x$-interval to the entire real axis, as needed for WH problems, with $\infty$ having two preimages ($x=\pm 1$).

For each inverse, we construct the Cauchy transform matrix. Following [16, 17], the contribution of $M^{-1}_{1}(\alpha)=\mathbb{I}$ to $C^{+}$ is given by

[TABLE]

where $F$ is the transform matrix from the function values at ${\bf x}^{\mathbb{I}}$ to the coefficients, $P$ is an almost Toeplitz matrix and

[TABLE]

Explicit expressions for $F$ and $P$ are given in [16, 17].

The contour $M^{-1}_{2}(\alpha)=\mathbb{R}\setminus\mathbb{I}$ is unbounded and connected to $\mathbb{I}$ at the points $\pm 1$. The contribution of this contour to $C^{+}$ is

[TABLE]

where $\mu^{L/R}$ are row vectors (where $L$ is identified with $-1$ and $R$ with $+1$) and $\Psi^{z}_{n}$ is an $n\times n$ matrix with the first and last rows removed, so that the matrix in (18) is of dimension $n\times n$. Both are defined in [16, 17] and are associated with the Cauchy transforms of the Chebyshev basis evaluated at points

[TABLE]

As in [16, 17], the orientation of the unbounded contour $M^{-1}_{2}(\alpha)$ (connected to $\mathbb{I}$ at $\pm 1$) implies that the first and last rows of the matrix in (18) must be the row vectors $\mu^{R}$ and $\mu^{L}$, respectively.

The final form of the Cauchy matrix $C^{+}$ for the rational mapping (14) is given by

[TABLE]

The final term is required to ensure that $C^{+}$ has the correct limiting behaviour [16]. Our numerical scheme takes into account the behaviour of the transform functions $\Phi_{\pm}$ for large $|\alpha|$. Their decay at infinity implies that the first and last rows in the matrices $C^{\pm}$ are irrelevant in our calculations.

However, this mapping does not possess the desired structure at the endpoints, since, e.g. for $x\sim 1$,

[TABLE]

and similarly for $x\sim-1$. For the diffraction problems we consider, the functions $\Phi_{\pm}$ are of the form $a_{1}\alpha^{-1/2}+a_{2}\alpha^{-1}+a_{3}\alpha^{-3/2}+\cdots$ and, therefore, their representations in this case will not be spectrally convergent Chebyshev series in $x$.

2.2.2 4-to-1 mapping

We introduce the rational mapping

[TABLE]

which has four inverses:

[TABLE]

and the solutions of a quadratic given by

[TABLE]

A schematic is given in figure 2. Again, the points $x=\pm 1$ are mapped to $\pm\infty$, but now the mapping will also be able to capture exactly the far-field behaviour of the functions $\Phi_{\pm}$. The critical point is that, for $x\sim 1$,

[TABLE]

and similarly for $x\sim-1$. One can hence represent functions of the form $a_{1}\alpha^{-1/2}+a_{2}\alpha^{-1}+a_{3}\alpha^{-3/2}+\cdots$ as spectrally convergent Chebyshev series in $x$.

Again, for each one of the inverses, we construct the Cauchy matrix. The contour $M^{-1}_{1}(\alpha)$ is identical to that presented above for the 2-to-1 mapping and, therefore, its contribution to $C^{+}$, $C^{+}_{1}$, is again given by (16).

The contour $M^{-1}_{2}(\alpha)=\mathbb{R}\setminus\mathbb{I}$ is unbounded and connected to $\mathbb{I}$ at the points $\pm 1$. The orientation of the contour is shown in figure 2. Its contribution to $C^{+}$ is given by (18) evaluated at $z=T^{-1}_{+}(M^{-1}_{2}(\alpha))$. Note that, consistent with [16, 17], the orientation of the unbounded contour $M^{-1}_{2}(\alpha)$ (connected to $\mathbb{I}$ at $\pm 1$) implies that the first and last rows of the matrix must be the row vectors $\mu^{L}$ and $\mu^{R}$, respectively.

The contours $M^{-1}_{3}(\alpha)=\{\mathrm{e}^{{\rm i}t}|t\in[0,\pi]\}$, $M^{-1}_{4}(\alpha)=\{\mathrm{e}^{{\rm i}t}|t\in[-\pi,0]\}$ are bounded and connected to $\mathbb{I}$ at the points $\pm 1$ (their orientation is shown in figure 2). The contributions of these contours to $C^{+}$ are given again by (18), evaluated at $z=T^{-1}_{+}(M^{-1}_{3}(\alpha))$ and $z=T^{-1}_{+}(M^{-1}_{4}(\alpha))$, respectively. Again, consistent with [16, 17], the orientation of the bounded contours $M^{-1}_{3}(\alpha)$, $M^{-1}_{4}(\alpha)$ (connected to $\mathbb{I}$ at $\pm 1$) implies that the first and last rows of the matrix must be the row vectors $\mu^{R}$ and $\mu^{L}$, respectively.

The final form of the Cauchy matrix $C^{+}$ for the rational mapping (22) is given by

[TABLE]

The fifth term has been subtracted to ensure that $C^{+}$ has the correct limiting behaviour (Olver [16]). Our numerical scheme takes into account the behaviour of the transform functions $\Phi_{\pm}$ for large $|\alpha|$. Their decay at infinity implies that the first and last rows in the matrices $C^{\pm}$ are again irrelevant for our calculations.

2.3 Rotation, scaling and evaluation

Once the matrix $C^{+}$ has been constructed (using either the 2-to-1 or 4-to-1 mapping), $C^{-}$ follows from the relation $C^{+}-C^{-}=I$. Therefore, we have a collocation method for the RH equation (2) which can be written as

[TABLE]

for $\alpha_{p}=M(x_{p})$, $p=1,\dots,n$, and $x_{p}$ given in (7). The WH problem requires the contour to lie in the strip of analyticity which encloses $\Gamma=\mathbb{R}$. For numerical purposes it is helpful to take the contour $\Gamma$ far from singularities, for example by taking a rotation of the real axis by angle $-\chi$, e.g. evaluating (27) along

[TABLE]

provided that no branch points are traversed. A RH problem can be then formulated along $\mathbb{R}\mathrm{e}^{-{\rm i}\chi}$; this poses no problem given the decay properties of the functions $\Phi_{+}$ and $\Phi_{-}$ and the fact that the rotation matrix providing (28) is invertible. Specifically, if we multiply (27) by a rotation matrix, then in order for (27) still to hold, we require functions $A(\alpha)$ and $B(\alpha)$ to be well-behaved in the region between the initial and rotated contour.

We also note that for the diffraction problems to be considered in the present study, there is one associated length scale, the diffraction parameter $k$. This implies that the mapping $M$ should be rescaled incorporating the length scale $k$. This is not pursued here, since we shall be considering the case $k=1$. For more complicated problems, e.g. the problem of Wickham and Abrahams [21, 12] discussed below, there are 2 associated length scales ($1$ and $k$), which should be taken into account in the scaling of the mapping.

To evaluate $\mathcal{C}\Phi$ at a point $\alpha=M(x)$ off the contour $\Gamma$, we use (4) and obtain $\mathcal{C}\Phi(\alpha)=\sum_{k=0}^{n-1}\check{U}_{k}\mathcal{C}T_{k}(x)$, where the coefficients $\check{U}_{k}$ are known e.g. from the numerical solution of the Riemann–Hilbert problem. To compute $\mathcal{C}T_{k}(x)$ for each of the inverses $M^{-1}_{j}(\alpha)$, we compute a row vector $\Psi^{z}$, analogous to the matrix $\Psi_{n}^{z}$, corresponding to the Cauchy transform of the Chebyshev polynomials at a single point $\alpha=M(x)$.

3 Error estimate

In the following sections we apply the numerical scheme presented above to solve various scalar and matrix WH problems. To validate our results, we compare computed values for the sectionally analytic functions $\Phi_{\pm}$ against those given by exact solutions. We use the error estimate

[TABLE]

($r>1$; for $r=\infty$ the estimate becomes the maximum error), where $Q(x)$ is the exact solution and $Q_{n}(x)$ is the numerical value at the collocation points. Note that (29) is a measure of the error in the mapped unit interval $\mathbb{I}$ in the $x$-plane. This integral can be evaluated using Clenshaw–Curtis quadrature [22] using the values of the integrand at the Chebyshev points. We write

[TABLE]

where $\{x_{q}|q=1,\dots,n\}$ are the Chebyshev points (7) and $w_{q}$ are weights that can be found in e.g. [20, 23].

Alternatively, the error estimate in the variable $\alpha$ can be used:

[TABLE]

where $R(\alpha)$ is the exact solution and $R_{n}(\alpha)$ is the numerical value for $n$ collocation points in the $\alpha$-plane. The error function $\mathcal{E}^{r}_{n}$ can be equivalently written as

[TABLE]

This can also be computed using Clenshaw–Curtis quadrature.

4 Diffraction by a semi-infinite plane

We analyse the problem of diffraction of a plane wave by a half-plane, where the plane is taken to lie along the positive real axis $x>0$, $y=0$. A schematic is shown in figure 3. The parameters $\rho_{1},\rho_{2}$ and $c_{1},c_{2}$ are the density and sound speeds in the upper (medium 1) and lower (medium 2) half-planes respectively. The boundary conditions along the top and bottom sides of the semi-infinite plane along the positive real axis are denoted by (B1) and (B2).

The governing equations for the fields $\phi_{1}$ and $\phi_{2}$ in the upper and lower half-planes, respectively, are given by

[TABLE]

where $k_{1}=\omega/c_{1}$, $k_{2}=\omega/c_{2}$ with $\omega$ the frequency. We decompose the total field $\phi_{t}$ into an incident wave $\phi_{\text{inc}}$ and a scattered field, with $\phi_{\text{inc}}$ given by

[TABLE]

where $k$ is the acoustic wave number and $\theta_{0}$ is the angle of the incident wave. Note that there is also an underlying time dependence, $\mathrm{e}^{-{\rm i}\omega t}$, that is only needed to determine branch cuts later.

The general form of the boundary conditions (B1) and (B2) on the semi-infinite plane $x>0$, $y=0^{\pm}$ is

[TABLE]

for $\lambda\in\mathbb{C}$ which can be different on the two sides of the semi-infinite plane. An acoustically hard plane corresponds to $\lambda=0$ while an acoustically soft plane corresponds to $|\lambda|\rightarrow\infty$.

Near the origin, we write

[TABLE]

using polar coordinates $x=r\cos\theta$, $y=r\sin\theta$, where $a$ and $h(\theta)$ are determined from the boundary and interface conditions.

Table 1 shows a list of diffraction problems by a half-plane with boundary conditions (B1) and (B2) along the top and bottom sides of the semi-infinite plane respectively, ratio of sound speeds $c_{1}/c_{2}$, ratio of densities $\rho_{1}/\rho_{2}$, the value of $a$ determining the local behaviour near the origin (36), as well as references to the original papers treating each case. Different notations and conventions are used in these references, and we follow them for ease of comparison, rather than recasting all problems in a single unified notation. Since our goal is to show the effectiveness of the numerical method presented in § 2, we do not pursue an exhaustive parameter study. We note that the numerical scheme is, in general, robust with regards to parameter changes, except for the cases $\theta_{0}\rightarrow\pi/2$ and impedance parameter $|\lambda|\rightarrow 0$ in (35) where convergence becomes slower. In the former case, this reduces the scope for using the angle of rotation used to separate the contour from the branch points. In the latter case, the behaviour near the origin which governs convergence properties could possibly change from the singularity structure considered in our numerical scheme near the origin ($a=1/2$).

We define full and half-range Fourier transforms in $x$ according to

[TABLE]

with

[TABLE]

This definition of the Fourier transform will be used in the next sections for the Sommerfeld problem (a) and the Senior problem (b). Note that Noble’s analysis[3], which will be followed to formulate both problems, uses (37) with a different normalising factor. For the Hurd problem (c), a different definition of the Fourier transform pair is used (following Hurd & Przezdziecki [28]).

4.1 The Sommerfeld problem [3, 24, 29]

First, we revisit the Sommerfeld problem. This problem can be formulated as a scalar Wiener–Hopf problem and admits an exact solution, since the associated kernel function can be factorized explicitly into a product of an upper and a lower analytic function.

For the diffraction problem originally solved by Sommerfeld [24], we have $c_{1}=c_{2}$, $\rho_{1}=\rho_{2}$, with hard-hard boundary conditions along $x<0$, $y=0^{\pm}$. (Following [3], the semi-infinite boundary is taken to lie along the negative real axis unlike in the rest of this paper). Then

[TABLE]

We can write $\phi_{t}=\phi_{\text{inc}}+\phi$, where $\phi_{\text{inc}}$ is an incident wave of the form (34) and $\phi$ is the scattered potential everywhere in the complex plane. Here we take $-\pi/2<\theta_{0}<\pi/2$; for $\theta_{0}$ outside this range, the analysis must be adapted. Noble [3] presents the following analysis stating that $0<\theta_{0}<\pi$, but a careful observation of the analyticity properties used indicates that the problem should really be considered separately for the two cases $0<\theta_{0}<\pi/2$ and $\pi/2<\theta_{0}<\pi$, even though the final result for directivity is the same. On the edge of the plate, $\phi$ is bounded and $\partial\phi/\partial y=O(x^{-1/2})$. Substitution of $\phi_{t}$ into (39), gives the boundary conditions for the scattered field $\phi$:

[TABLE]

4.1.1 Exact solution [3]

It is convenient to assign a small positive imaginary part to $k\mapsto k+{\rm i}\epsilon$, $\epsilon>0$, to improve convergence of subsequent Fourier integrals (and then let $\epsilon\rightarrow 0$).

Since $\mathrm{d}\Phi/\mathrm{d}y=\Phi^{\prime}$ is continuous everywhere, we write

[TABLE]

where $\gamma(\alpha)=(\alpha^{2}-k^{2})^{1/2}$ which has branch cuts from $\pm k$ to $\infty$ in the first and third quadrants, as required by causality. The function $A(\alpha)$ has to be determined from the boundary conditions. It follows that

[TABLE]

Expressions with one argument, e.g. $\Phi_{+}(0)$, will always refer to the value of, in this case, $\Phi_{+}(\alpha,y)$ at $y=0$. Taking the half-range Fourier transform of the boundary condition (40) gives

[TABLE]

Elimination of the function $A(\alpha)$ from (42) and use of (43) gives

[TABLE]

where

[TABLE]

The two equations in (44) hold in the strip $-k_{i}<\text{Im }\alpha<k_{i}\cos{\theta}_{0}$ where $k_{i}=\text{Im }k$. Next, we divide the second equation in (44) by $(\alpha+k)^{1/2}$:

[TABLE]

The first term on the left-hand side is upper analytic, while the right-hand side is lower analytic. Applying an additive splitting to the second term gives

[TABLE]

Therefore, (46) can be written in the form

[TABLE]

The function $J(\alpha)$ is regular in the whole plane by analytic continuation and is bounded from edge conditions. Hence, by Liouville’s theorem, $J(\alpha)=0$, implying

[TABLE]

The function $A(\alpha)$ follows from (42).

4.1.2 Numerical solution using the RH formulation

We now solve the hard-hard Sommerfeld problem numerically by expressing it as a RH problem on the (rotated) real line. Recall that

[TABLE]

We divide by $\gamma$ so that

[TABLE]

where $\Phi_{+}^{\prime}(0)=\mathcal{C}^{+}\Phi$ and $D_{-}=\mathcal{C}^{-}\Phi$. This becomes the matrix problem $M\Phi=N$, where

[TABLE]

The Cauchy operator $C^{+}$ is given by (20) or (26) and $C^{-}=C^{+}-I$. The sectionally analytic functions $\Phi^{\prime}_{+}(0)$ and $D_{-}$ can be computed by applying the operators $C^{+}$ and $C^{-}$ to $\Phi$.

Figure 4 shows the error estimates $E^{2}_{n}$, $E^{\infty}_{n}$, $\mathcal{E}^{2}_{n}$ for $\Phi^{\prime}_{+}(0)$ and $D_{-}$ using the 2-to-1 and 4-to-1 mappings, as functions of $n$ (the number of collocation points). We observe that spectral convergence is obtained using the 4-to-1 mapping, since, as discussed earlier, this mapping is able to capture exactly the far-field behaviour of the sectionally analytic functions. On the other hand, the error estimate using the 2-to-1 mapping is much larger, e.g. $E^{2}_{n}=\mathcal{O}(10^{-3})$ using $n=100$; this is to be expected, since the numerical scheme using this mapping does not incorporate the square root singularity near the origin which governs convergence properties [30].

4.2 Senior’s problem [3, 26]

We take $c_{1}=c_{2}$, $\rho_{1}=\rho_{2}$ with impedance boundary conditions along $x>0$, $y=0^{\pm}$:

[TABLE]

where $S$ is the impedance parameter. Again, we can write $\phi_{t}=\phi_{\text{inc}}+\phi$, where $\phi$ is the scattered potential everywhere in the complex plane and $\phi_{\text{inc}}$ is the incident wave (34), but now with $\pi/2<\theta_{0}<3\pi/2$ (this is omitted in [3]).

4.2.1 Analysis [3]

We write

[TABLE]

where $\gamma=(\alpha^{2}-k^{2})^{1/2}$ and $A(\alpha)$, $B(\alpha)$ have to be determined from the boundary conditions. This gives (using the same notation as previously):

[TABLE]

Elimination of $A(\alpha)$, $B(\alpha)$, substitution of $\phi_{t}=\phi_{\text{inc}}+\phi$ in (53) and the definition of the Fourier transform yield

[TABLE]

These equations can be expressed in matrix form as

[TABLE]

where

[TABLE]

Addition and subtraction of (56) and (57) give two independent scalar WH problems:

[TABLE]

4.2.2 Numerical solution using the RH formulation

To solve this problem numerically, we express it as a RH problem. Use of (58) results in the matrix RH problem $M\Phi=N$, where

[TABLE]

and

[TABLE]

The functions $\Phi_{+}$ and $\Phi_{-}$ can be computed by applying the operators $C^{+}$ and $C^{-}$ to $\Phi$. One can also solve the two scalar WH problems given by (60) and (61) numerically.

Figure 5 shows the error estimates $E^{2}_{n}$, $E^{\infty}_{n}$, $\mathcal{E}^{2}_{n}$ for $\Phi^{\prime}_{+}(0^{+}),\Phi^{\prime}_{+}(0^{-}),\Phi^{\prime}_{-}(0),\Phi_{-}(0)$ using the 4-to-1 mapping and both scalar and matrix formulations, as functions of $n$. Again, spectral convergence is observed, since the 4-to-1 mapping is able to capture exactly the far-field behaviour of the sectionally analytic functions. The exact solutions used to compute the error estimates are provided from the highest resolution ($n=257$) numerical solutions.

4.3 Hurd’s problem [9, 28, 31, 27]

Hurd [9] presented a method to deal with matrix WH problems with symmetric branch cuts at $\pm k$ and isolated singularities in the upper (or lower) half-plane. The analysis is based on the so called Wiener–Hopf–Hilbert method which involves the transformation of WH equations into a pair of coupled Hilbert equations that can be solved using Muskhelishvili’s theory [32]. The problem was revisited by Hurd & Przezdziecki [28] who carried out a more detailed analysis. Hurd’s problem [9] is similar to that analysed by Senior [3, 26], but with different impedance boundary conditions on top and bottom sides of the semi-infinite plane:

[TABLE]

where $S_{j}=\sin\theta_{j}$, $0\leq\theta_{j}\leq\pi/2$, for $j=1,2$. Again, we write $\phi_{t}=\phi_{\text{inc}}+\phi$, where $\phi_{\text{inc}}$ is the incident wave (34) with $-\pi/2<\theta_{0}<\pi/2$ (see figures 1 and 2 in [28]) and $\phi$ is the scattered potential everywhere in the complex plane. This problem was also analysed by Rawlins [31] using a different approach and, later, by Barton & Rawlins [27] who modified the Wiener–Hopf–Hilbert method to consider the case when surface waves can propagate along the two sides of the semi-infinite plane.

We write

[TABLE]

where $\beta(\alpha)=(k^{2}-\alpha^{2})^{1/2}$ such that $\beta(0)=k$ and $A(\alpha)$, $B(\alpha)$ have to be determined from the boundary conditions. Following Hurd & Przezdziecki [28], the above problem leads to a matrix WH problem of the form:

[TABLE]

where

[TABLE]

and $S_{0}=\sin\theta_{0}$. Note that the definition of the Fourier transform pair in [28] differs from that given by (37).

4.3.1 Numerical solution using the RH formulation

To solve this problem numerically, we express it as a RH problem. We write $M\Phi=N$, where

[TABLE]

and

[TABLE]

The functions $\Phi_{+}$, $\Phi_{-}$ can be computed by applying the operators $C^{+}$ and $C^{-}$ to $\Phi$.

Figure 6 shows the error estimates $E^{2}_{n}$, $E^{\infty}_{n}$, $\mathcal{E}^{2}_{n}$ for $U_{1},U_{2},L_{1},L_{2}$ using the 4-to-1 mapping, as functions of $n$. We observe spectral convergence, since, as already mentioned, the 4-to-1 mapping is able to capture exactly the far-field behaviour of the sectionally analytic functions. The exact solutions used to compute the error estimates are provided from the highest resolution ($n=257$) numerical solutions. Although closed-form expressions are given by Hurd [9], we have not used those to verify our results since their computation is awkward.

5 Far-field directivity

The full range Fourier transform was written in the form

[TABLE]

(For the Sommerfeld problem: $B(\alpha)=-A(\alpha)$.) Taking the Fourier inverse of (70), the diffracted field $\phi$ can be written as

[TABLE]

in the upper half-plane. To obtain the far-field asymptotics ($x,y\rightarrow\infty$), we deform the integration contour on the steepest descent path [3]. Note that the stationary phase method can alternatively be used [29]. Application of the steepest descent method to (71) gives

[TABLE]

where $(r,\theta)$ are polar coordinates. A similar representation for $\phi(r,\theta)$ in the lower half-plane ($-\pi\leq\theta\leq 0$) can be found. The far-field directivity, $D(\theta)$, is defined by (3).

Note that the formula (72) breaks down near the shadow boundaries. The reason is that the pole $\alpha=k\cos\theta_{0}$ is near the saddle point $\alpha=-k\cos\theta$ and, therefore, the steepest descent method which was employed to obtain the far-field asymptotics (72) fails. To overcome this, we can proceed as in [3]. Suppose that $A(\alpha)$ can be written as

[TABLE]

where function $H(\alpha)$ varies slowly near the saddle point $\alpha=-k\cos\theta$. The form of $H(\alpha)$ depends on the problem under consideration. Then one can break up the solution $\phi(x,y)$ into two terms, one of which is uniformly valid while the other takes the form of a plane wave in the far field. We do not present results from this procedure here, since our goal is to show that the numerical method reproduces the Fourier transform $A(\alpha)$ accurately and we compare to exact results that are not uniformly valid. If desired, one could then use the modification presented in [3] to obtain uniformly valid solutions numerically.

The Sommerfeld problem Using (72) and $B(\alpha)=-A(\alpha)$, it follows that the directivity is given by

[TABLE]

where

[TABLE]

Expression (75) follows from (44). If the argument of $A(\cdot)$ in (74) is in the upper half-plane, then we use the first equality in (75) to compute $A$. In this case, the function $\Phi_{-}^{\prime}(\alpha,0)$ is defined in (43), while we use the numerical solution to the Sommerfeld problem above to compute $\Phi_{+}^{\prime}(\alpha,0)$. If the argument of $A(\cdot)$ in (74) is in the lower half-plane, then we use the latter equality in (75) to compute $A$; in this case, $D_{-}$ is computed from our numerical solution.

The exact far-field directivity for the hard-hard Sommerfeld problem [3] which is used for comparison to our numerical solutions is given by

[TABLE]

Senior’s problem Using (55), the functions $A(\alpha)$ and $B(\alpha)$ can be expressed in terms of the sectionally analytic functions as

[TABLE]

and the directivity can be computed using (72). Similarly to the Sommerfeld problem, different but equivalent representations for the functions $A$ and $B$ using (56)–(57) can be written which are then used to compute the far-field directivity numerically.

Hurd’s problem Following Hurd & Przezdziecki [28], the functions $A(\alpha)$ and $B(\alpha)$ can be expressed as

[TABLE]

with additional equivalent representations provided by (66)–(67), and directivity can be computed from application of the steepest descent method to the analogous expression (71), since the definition of the Fourier transform pair in [28] is different, with $\gamma$ replaced by $\beta$ (consistent with (65)) and different normalising factor. An analogous representation to (72) can be written with $A(-k\cos\theta)$ replaced by $A(k\cos\theta)$.

Bowman [33] gave the far-field asymptotics for different impedance boundary conditions along the top and bottom sides of the semi-infinite plane:

[TABLE]

where $U(\theta,\theta_{0})$ is given by

[TABLE]

where $\psi(x)$ is a known function [33]. Other expressions are given in [28, 34, 35]. The expression (79) is used for comparison to our numerical solutions for Senior’s and Hurd’s problems.

Figure 7 shows the the far-field directivity $|D(\theta)|$ for the hard-hard Sommerfeld, Senior and Hurd problems, as a function of angle $\theta$, computed using the exact and numerical solutions. The numerical results are indistinguishable from the exact directivities; this reflects the excellent approximation properties of our numerical schemes. We are computing approximations to the Fourier transforms at the stationary phase points. Since the transforms are analytic in their respective half-planes, we obtain spectral convergence. We observe that our numerical method recovers the singularities in the far-field directivities associated with the shadow boundary regions.

6 Discussion

We have presented fast and accurate numerical schemes for the solution of scalar and matrix WH problems by exploiting their links with RH problems. The idea is to solve the corresponding RH problems adapting the methods of Olver and Trogdon [16, 17] to take into account the known far-field behaviour of the solutions to construct tailor-made numerical schemes. In particular, we have used rational mappings with multiple inverses, extending the Möbius mappings discussed in [17] to account for the $\mathcal{O}(\alpha^{-1/2})$ decay of the solutions for large $|\alpha|$.

The 4-to-1 mapping was used to capture the $\alpha^{-1/2}$ behaviour of the Fourier transforms that stems from the $r^{1/2}$ behaviour near the origin found in diffraction problems with e.g. hard-hard and impedance-impedance boundary conditions. For soft-hard boundary conditions, $\phi_{j}(r,\theta)\sim r^{1/4}h_{j}(\theta)$, $j=1,2$ and, therefore, the numerical scheme presented here must be adapted (e.g. using an 8-to-1 rational mapping).

We analysed problems of diffraction of plane waves by a semi-infinite plane with different boundary conditions that produce scalar and matrix WH problems and obtained accurate and spectrally convergent solutions using our numerical schemes. Our results were verified against exact solutions for the cases for which these exist, as well as high-resolution numerical solutions. We computed the far-field directivity obtained from the far-field asymptotics, which is the most useful physical prediction from these problems.

The numerical scheme presented in this study can be adapted to solve other more complicated problems some of which we mention below. Scattering by wedges and other simple geometric shapes can lead to WH problems (for a discussion of many scattering problems, see [36]), as can linear elasticity problems with straight-line geometry.

Wickham [21] and Abrahams [12] analysed the problem of scattering of acoustic waves by a semi-infinite screen at the interface between two compressible media with different physical properties, so that $c_{1}\neq c_{2}$, $\rho_{1}\neq\rho_{2}$. The problem was reduced to a matrix WH problem of the form

[TABLE]

where

[TABLE]

and $\gamma(\alpha)=(\alpha^{2}-k^{2})^{1/2}$, $k>1$, $\delta(\alpha)=(\alpha^{2}-1)^{1/2}$, with the same branch cut structure as $\gamma(\alpha)$, and $\mu=\sqrt{\rho_{2}/\rho_{1}}$. Local analysis near the origin [21, 12] gives $\phi_{j}(r,\theta)\sim r^{a}h_{j}(\theta)$, $j=1,2$, with $a=\pi^{-1}\tan^{-1}{\mu}$ . Therefore, to construct accurate numerical schemes for solving problems with solutions exhibiting irrational $a$, we need to consider a different approach, e.g. expanding in Jacobi polynomials rather than the Chebyshev polynomials considered here. This is beyond the scope of this work, but remains a topic for future study.

The present numerical approach can also be adapted to devise accurate numerical schemes for WH problems with jump matrices containing exponential factors [37, 38, 15], e.g. of the form

[TABLE]

Problems with exponentials $\mathrm{e}^{\pm{\rm i}\alpha L}$ in the elements of the jump matrix arise in acoustics when there are different boundary conditions in the physical space $(-\infty,0)$, $(0,L)$ and $(L,\infty)$, for example in the scattering of acoustic waves by multiple slits or poroelastic plate extensions.

Finally, motivated by various applications such as calculating effective properties of periodic composites [39], we aim to adapt the method for solving RH problems with jump matrices defined piecewise along $\mathbb{R}$. To solve such problems, we have to construct numerical schemes and rational mappings with different behaviours at critical points along the contour.

Data Accessibility. This article has no additional data.

Authors’ Contributions. SGLS had the original idea for this work. The calculations were carried out by both authors equally. Both authors contributed to the writing of the manuscript.

Competing Interests. The authors declare that they have no competing interests.

Funding. This research was partly funded by NSF award DMS-1522675.

Acknowledgments. Conversations with David Abrahams, Sheehan Olver and Tom Trogdon were very helpful.