On the numerical solutions of two-dimensional scattering problems for an open arc

For the scattering problems of acoustic wave for an open arc in two dimensions, we give a uniqueness and existence analysis via the single layer potential approach leading to a system of integral equations that contains a weakly singular operator. For its numerical solutions, we describe an O(h3)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$O(h^{3})$\end{document} order quadrature method based on the specific integral formula including convergence and stability analysis. Moreover, the asymptotic expansion of errors with odd power O(h3)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$O(h^{3})$\end{document} is got and the Richardson extrapolation algorithm (EA) is used to improve the accuracy of numerical solutions. The efficiency of the method is illustrated by a numerical example.


Introduction
The scattering problems of electromagnetic waves or time-harmonic acoustic waves by an infinitely long semi-cylindrical obstacle with a smooth open contour cross-section Γ ⊂ R 2 are reduced to the Helmholtz equation with unbounded boundary value problems where Γ is a smooth open arc, k > 0 is the wave number, and |x| is the Euclidean distance. The solution of (1) is got by the single layer potential theory in the following form [1]: where K 0 |x -y| is the fundamental solution of (1), which is described as and H (1) 0 = iN 0 + J 0 is the Hankel function of order zero and of the first kind.
© The Author(s) 2020. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. N 0 is the Neumann function of order zero, which is given by N 0 (z) = 2 π ln z 2 + γ J 0 (z) + 2 π ∞ n=1 n m=1 1 m (-1) (n+1) (n!) 2 z 2 2n (4) with γ = 0.57721 . . . denoting the Euler constant, and J 0 is the Bessel function of order zero, which is given by Equations (4) and (5) are valid for small |z|. For large |z|, they are inaccurate, and other series expressions are given in [2]. We decompose the integral kernel to analyze its properties as follows: where k 1 (x, y) = -1/2π ln |x -y| is a logarithmic singular function and The singe layer potential (2) solves the unbounded Dirichlet problem provided the density v is the solution of the boundary integral equation which is the first kind boundary integral equation with weakly singular kernel. Assume that C Γ denotes the logarithmic capacity of Γ , when C Γ = 1, the solution of (7) exists and is unique. Once the solution v(x) is solved by (7), the value u(x) (x ∈ R 2 \Γ ) can be got by (2). Boundary element method (BEM) [1,3], finite element method (FEM) [4], and some meshless methods [5,6] have been developed for numerical solutions of the Helmholtz equation. Huang and Wang [7] considered numerical solutions of the boundary integral equation of the Helmholtz equation with mixed boundary value problem on a smooth closed curve by the mechanical quadrature method [8] and extrapolation algorithm. In this case, the solution of the equation is smooth. The boundary element-free method was proposed by Chen, Liu, and Li [2] for solving two-dimensional exterior and interior Helmholtz equations with mixed boundary value problem. The large wave numbers Helmholtz equation was solved by Wang, Zhai, and Zhang [9] using a new weak Galerkin mixed finite element method. Dastour and Liao [10] presented the finite difference method for solving the two-dimensional Helmholtz problems. A radial basis function-generated finite difference scheme was introduced by Londono and Montegranario [11] to get the numerical solutions of the two-dimensional Helmholtz equation. Celiker and Lin [12] solved the Helmholtz equation with the Dirichlet boundary value problem by using a conforming finite element method and obtained the high precision numerical solutions, etc. In the above numerical methods, the boundary element method is a good choice due to the reduction of the dimension of boundary value problems. Boundary element method can be used to transform boundary value problems into boundary integral equations which include singular integral [1]. There are numerical methods to solve the boundary integral equations (7), such as Galerkin method [13], collocation method [14], and meshless methods [15]. However, these methods are either of low accuracy or complicated or unstable.
The focus of the paper is to investigate the numerical solutions of the Helmholtz equation with Dirichlet problems for the smooth open arc by using a high-accuracy quadrature method which has the characteristics of high-accuracy, low computational complexity, and good stability. We firstly convert the problem into a weakly singular boundary integral equation by using the potential theory. Because the boundary integral equation is defined on the smooth open arc, its solution is singular at the end of the open arc, which seriously affects the accuracy of the numerical solutions. In order to obtain high accuracy numerical solutions, we use Sidi transform to eliminate the singularity of solution at the end of open arc. Further, we construct a high-accuracy quadrature method according to the specific quadrature formula and discrete the equation to get a system of linear equations. Then, we prove the convergence of the numerical solutions and the stability of the method, and prove that the error has a single parameter asymptotic expansion with odd power O(h 3 ), and the discrete matrix has a good condition number. The higher accuracy can be obtained by the EA. In the end, numerical results support our theoretical analysis.
The remainder of the paper is organized as follows. We analyze the singularities of the solutions and integral kernels in Sect. 2. A quadrature method with high-accuracy is constructed to discrete the integral equation by using singular integral formula in Sect. 3. Furthermore, we prove the convergence and stability of the method. In Sect. 4, we get the single parameter asymptotic expansion of the error and construct the EA based on the expansion to improve the accuracy of the numerical solutions. A numerical example is given in Sect. 5 to verify our theoretical analysis.

Singularity analysis of integral kernels and solutions
For the convenience of discussion, we assume that the boundary Γ is described by the parameter mapping: 1], m ∈ N and also assume y = x(s) = (x 1 (s), x 2 (s)), s ∈ [0, 1]. Under the parameter mapping, (7) can be converted into the following integral equation: with g(s) = g(x(s)). Since the function v is singular at t = 0 and t = 1, we use Sidi transformation [16] to eliminate singularity, which is defined as follows: where where x(τ ) = x(ψ p (τ )), v(τ ) = v(x(ψ p (τ ))), and ψ p (τ ) has zero points with degree p at τ = 0 and τ = 1.
Some integral operators on [0, 1] are defined as follows: with k 2 (s, τ ) = k 2 (x(s), x(τ )). From (11), (12), and (13), we transform (10) into the following matrix operator equation: (14) is equivalent to where E is the identity operator and G = T -1G . Next, we analyze the singularity of the solution v of (7), Let Q 1 and Q 2 be the end of Γ , and assume the endpoints correspond to the parameter values t m (m = 1, 2). It is well known that the solution v has singularity at points t m (m = 1, 2) [17,18]. The solution v can be described as the following form: where g(0) = 0, g(1) = 0, and β m ≥ -1 2 . Further, we get and where d 1 and d 2 are constants. From (17) and (18), it is easy to know that w(τ ) is bounded in [0, 1] and has zeros at τ = 0 and τ = 1.
Then we get the following results hold:

High-accuracy quadrature method
There exists a quadrature formula as follows: and the error is given by where ζ (t) is the derivative of the Riemann zeta function and h is the mesh widths.
Let us assume that h = 1/n (n ∈ N ) denotes mesh widths, and s j = τ j = (j -1/2)h, j = 1, . . . , n, are middle nodes. For an integral operator M with continuous kernel m(s, τ ) as K 1 and K 2 , by the midpoint or the trapezoidal rule [19], the Nyström approximation M h is defined with the error bounds For the logarithmic singular operator T, by the quadrature formula (19), the Fredholm approximation T h is defined with the error bounds From (21) and (23), we can get the approximate equation of (14) as follows: where w h = (w(τ 1 ), . . . , w(τ n )) andG h = (G(τ 1 ), . . . ,G(τ n )). It is obvious that (25) is a linear equation with n unknowns. The solution w h is solved by (25), then the value u(y) (y ∈ R 2 \Γ ) can be got by the following form: To discuss a unique solution existing for (25), we firstly obtain that the operator T h is invertible. From (23) we know that T h is the circular matrix, which is expressed by Lemma 2 There are two constants c 1 > 0 and c 2 > 0 independent of h, which make the eigenvalues λ α (α = 1, . . . , n) of matrix T h satisfy c 1 ≥ λ α ≥ c 2 h.

Based on the above corollary, (25) is equivalent to
In the following, we will discuss the convergence of the numerical solutions of (32). First, some special operators and subspaces are introduced.
A subspace C 0 [0, 1] ⊂ C[0, 1] is defined as with a norm v * = max 0≤t≤1 |v(t)/ sin 3 (πt)|. A piecewise linear function subspace is defined as where the basis function e j (t) satisfies e j (t i ) = δ ji and {t i } n m -1 i=0 denotes the basis points. A prolongation operator is defined as L h : R n → S h and it satisfies A restricted operator is defined as F h : C 0 [0, 1] → R n and it satisfies Next, we provide some lemmas to prove the convergence.
We get e = (T h ) -1 η with e = O(h 2 ) by Corollary 1, and is the collectively compact convergent toM: Proof According to (35), we get and Since d dsd (s, τ ) is a continuous function on [0, 1] 2 , (36) can be obtained.
Further the following corollary can be got. (21), we have

Corollary 2 Let the Nyström approximation M h of M be defined by
Proof Under the parameter transformations (9), the kernel m(s, τ ) of M is continuous, and its high order derivatives are also continuous. From Also using the results of [21] and Lemma 3, we conclude that {T -1 M h : C[0, 1] → C 3 [0, 1]} must be collectively compact convergent to T -1 M. Hence, the proof is completed.
We construct an operator equation Obviously, ifŵ h is the solution of (38), then F hŵh must be the solution of (32); conversely, if w h is the solution of (32), then L h w h must be the solution of (38). Below we prove that there exists a unique solutionŵ h in (38) such thatŵ h converges to w. Finally, the main theorem is proved to complete the proof of convergence.
where M h stands for K h 1 and K h 2 . Thus we get → shows the point of convergence. The proof is completed.
Meanwhile, the following corollary is presented to obtain the stability of the method.

Corollary 3
Assume that K h 1 , K h 2 , and T h are defined by (21) and (23), respectively. The eigenvalues of Υ h = T h +K h 1 +K h 2 are λ i (i = 1, . . . , n). There is the bound of condition number

Asymptotic expansion of error and EA
In what follows, the single parameter asymptotic expansion of error is got to describe the EA.
Theorem 2 There exists a function Φ independent of h such that the following asymptotic expansion with a single parameter holds at nodes: where w h ∈ S, w ∈ C 3 [0, 1].
Based on the single parameter asymptotic expansion (41), the EA can be constructed as follows [22]: By the extrapolation results, the following a posteriori estimate can be obtained: The self-adaptive algorithm can be got by using the above inequality. Obviously, the EA is not complicated, but the accuracy of the numerical solution is improved.

Numerical example
In what follows, we provide a numerical example to verify the proposed method.

Conclusions
In this paper, we discuss the high-accuracy quadrature method and the EA for scattering problems for a smooth open arc, there exist the following advantages: (a) Computing an entry of the discrete matrices is straightforward and simple, without any singular integrals, hence the method is appropriate to solve weakly singularity problems; (b) The method possesses the feature of high accuracy. The O(h 3 ) order accuracy of the error is proved, and the numerical results verify our theoretical analysis. The EA can be used to raise the accuracy of the solutions, see Tables 1 and 2; (c) The condition number of the discrete matrix increases linearly, which shows that the method is stable. Although the method has the characteristics of high accuracy, low computational complexity, and good stability, it also has a limitation, that is, the method is effective for solving weakly singular boundary integral equations, but not for strongly singular and hypersingular boundary integral equations. Based on the research results in [23,24], we intend to apply the method to the Neumann problems for the Helmholtz equation and the exterior acoustic problems with arbitrary and high wave numbers in the future.