1 Introduction↩︎

Geometric discontinuities such as corners, edges, and junctions are fundamental features of models in continuum mechanics. They occur across a wide range of scales, from crack tips and material interfaces to tectonic faults, boundaries between geological strata, and subsurface scatterers in geophysics [1]–[4]. In elastostatics and related elliptic problems, these nonsmooth geometries generate stress fields with corner singularities that may be oscillatory and strongly localized. The resulting behavior influences scattering patterns, delamination, and failure. Consequently, developing a high order numerical scheme is essential for analyzing these singular structures.

Traditional numerical techniques, particularly the finite element method, are widely used in structural analysis, but they struggle to resolve the branch-type singular behavior near corners. Asymptotics at the material junctions often involve fractional and occasionally complex corner exponents, which demand excessive local refinement and can produce unstable artifacts [5]. Beyond volume-based discretizations, there exists the boundary-oriented class of approaches that represent or approximate the solution from boundary data. Popular schemes include the method of fundamental solutions [6] and the Lightning solver [7]–[9], which place sources or poles that cluster exponentially near corners and capture the expected local behavior. Although these methods are meshless and can achieve high-order accuracy in certain geometries, their underlying theoretical foundations are still being established.

In this paper, we focus on boundary integral equation methods (BIEs) formulated via layer potential theory in elastostatics [10]–[12]. BIEs not only reduce dimensionality but also naturally handle the far-field conditions [13]. However, while BIEs are highly accurate on smooth geometries [14]–[16], their application to nonsmooth domains with cracks, interfaces, and sharp re-entrant corners is complicated by two issues. First, the kernel function for the second-kind formulation is a tensor and involves Cauchy singular integrals, which complicate the numerical analysis. Second, one must account for branch singularities in the unknown boundary density, which can be oscillatory and strongly coupled near corners [17], [18].

To address these issues, various numerical strategies have been proposed. Approaches include kernel splitting techniques combined with generalized Gaussian quadrature [19]–[21], which have proven to be effective for handling logarithmic and Cauchy singularities. Another prominent approach for cornered domains involves local mesh refinement, often combined with compression techniques to control the size of the linear system, including the \(L^2\) weighting method [22] and recursive compressed inverse preconditioning [23]. Alternatively, geometric regularization via corner rounding to modify the domain so that the most severe singular features are suppressed has also been proposed in [4]. Recently, Rokhlin and Serkh constructed an effective quadrature method by utilizing explicit representations of the solution’s asymptotic behavior near corners [24], [25]. They observed that solutions to these integral equations can be represented by rapidly convergent singular power series for Laplace problems, fractional order Bessel functions for Helmholtz cases, and expressions of the form \(\sum_j\left(c_j t^{\mu_j} \sin \left(\beta_j \ln (t)\right)+d_j t^{\mu_j} \cos \left(\beta_j \ln (t)\right)\right)\) for Stokes flow, as summarized in Table 1. Numerically, explicit singular basis methods require the prior determination of the characteristic exponents. To ensure stability and accuracy, the resulting bases are orthogonalized by SVD, yielding a well-conditioned set of functions. These expansions are then discretized using generalized Gaussian quadrature rules tailored to the singular behavior.

Although the singularity analysis proposed in [24] has been successfully extended to vectorial systems such as the 2D Stokes equations [26], its application to elastic equations remains underdeveloped. In particular, the Cauchy-type singularities of the elastic Green’s kernel prevent the parity transform from decoupling the system. A fully coupled \(4\times4\) boundary formulation makes both the analysis of the corner exponents and the design of efficient algorithms substantially more involved. This paper fills this gap by establishing a complete analytical and numerical framework for corner singularities in two dimensional elastostatics. Our contributions are threefold:

1. Corner exponents from a factorizable Mellin symbol: On a wedge, we derive a characteristic equation for the boundary density’s singular exponents that factorizes into two branches, identifying real and complex exponents and explaining which modes are compatible with interior vs.exterior conditions. This connects the BIE symbol calculus directly to classical PDE corner analysis.

2. Fredholm criterion and density-to-Taylor mapping: In the weighted Sobolev space \(H^s_\nu(\Gamma)^2\), we obtain a computable Fredholm criterion \(s-\nu\notin\{\Re z_{n,j}\}\cup\mathbb{Z}\) for the interior double-layer operator, with \(z_{n,j}\) denoting the zeros of the Mellin symbol. We also give an explicit density-to-Taylor map \(\mathbf{B}(\theta)\) and prove its invertibility for all but a discrete set of angles, thereby linking corner singular coefficients to smooth boundary data.

3. A high order singularity guided Nyström (SGN) scheme: Rather than spanning the density with singular powers, whose hybrid variant is severely ill-conditioned, we use the computed exponents to drive panel adaptivity via a multi-exponent Legendre-tail indicator. Moreover, we prove an a priori error bound that couples this indicator with exponentially small far-field quadrature error. In numerical experiments, we observe that the SGN scheme achieves higher accuracy than a uniform Nyström baseline. Moreover, on re-entrant corners, the attainable accuracy is crowding-limited by near singular integration.

The remainder of the paper is organized as follows. Section 2 formulates the elastostatic BIEs and develops the well-posedness analysis using the Fredholm index theory in weighted Sobolev spaces. Section 3 presents the mathematical apparatus, including Mellin-type integrals and a characterization of the root distribution for \(\operatorname{sinc}\) functions. Section 4 performs the local corner analysis, derives the characteristic equation, and explicitly constructs the density-to-Taylor mapping. Section 5 presents the singularity guided Nyström scheme, error analysis, and several numerical experiments. Section 6 summarizes the paper and discusses future directions. Technical details are deferred to Appendices 7 and 8.

2 Problem formulation↩︎

Consider a bounded, simply connected planar elastic domain \(\Omega\subset \mathbb{R}^2\) with a piecewise smooth boundary \(\Gamma\). The boundary is assumed to contain a finite number of corners, each locally diffeomorphic to a wedge with an opening angle (cusps excluded). Let \(\boldsymbol{n}=(n_1, n_2)^\top\) and \(\boldsymbol{\tau}=(\tau_1, \tau_2)^\top\) denote the unit outward normal and the positively oriented unit tangent on \(\Gamma\) respectively, where \(\tau_1 = -n_2\) and \(\tau_2 = n_1\), as shown in Figure 1.

Assume the domain \(\Omega\) is filled with a homogeneous, isotropic elastic material. The displacement field \(\boldsymbol{u}=(u_1,u_2)^\top\) under static loading satisfies the homogeneous Lamé system: \[\label{eq:lame95equation} \mathcal{L} \boldsymbol{u} \triangleq \mu \Delta \boldsymbol{u} + (\lambda + \mu) \nabla (\nabla \cdot \boldsymbol{u}) = \boldsymbol{0}, \quad \boldsymbol{x}\in \Omega,\tag{1}\] where \(\lambda\) and \(\mu\) are the Lamé parameters. The conditions \(\mu > 0\) and \(\lambda + \mu > 0\) ensure the strong ellipticity of \(\mathcal{L}\) [17]. Throughout, we assume that all material parameters are constant and focus on static elasticity, although the methodology extends naturally to quasi-static or time-harmonic elastic wave scattering [21].

Two standard types of boundary conditions are considered. Under Dirichlet boundary conditions, the displacement field satisfies \[\label{eq:HC95bc} \boldsymbol{u}(\boldsymbol{x}) = \boldsymbol{f}(\boldsymbol{x}), \quad \boldsymbol{x}\in \Gamma\tag{2}\] with \(\boldsymbol{f} \in H^{1/2}(\Gamma)^2\). Under Neumann boundary conditions, the traction field satisfies \[\label{eq:SF95bc} \mathcal{T} \boldsymbol{u}(\boldsymbol{x}) = \boldsymbol{g}(\boldsymbol{x}), \quad \boldsymbol{x}\in \Gamma\tag{3}\] with \(\boldsymbol{g} \in H^{-1/2}(\Gamma)^2\). The traction operator \(\mathcal{T}: H^{1}(\Omega)^2 \to H^{-1/2}(\Gamma)^2\) is defined by \(\mathcal{T} \boldsymbol{u} \triangleq \boldsymbol{\sigma}(\boldsymbol{u}) \boldsymbol{n}\), where the stress tensor \(\boldsymbol{\sigma}(\boldsymbol{u})\) satisfies \[\label{eq:stress95tensor} \sigma_{ij}(\boldsymbol{u}) = \lambda(\nabla \cdot \boldsymbol{u})\delta_{ij} + \mu\left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right), \quad i,j=1,2.\tag{4}\] Notably, the analytic framework admits natural extensions to various boundary conditions, which are omitted here for clarity. The central difficulty in such problems lies in the singular behavior of the stress field near corners, which is essential for predicting failure and ensuring structural integrity.

2.1 Boundary integral formulation↩︎

To efficiently resolve the singularity near corners, we employ a boundary integral formulation based on layer potential theory. Let \(\boldsymbol{\phi}:\Gamma \to \mathbb{R}^2\) be a vector-valued density function defined on the boundary. The fundamental solution (Kelvin matrix [27]) of the 2D Lamé system in plane strain is given by \[\label{eq:fundamental95solution} \mathbf{G}(\boldsymbol{x}, \boldsymbol{y}) = \frac{1}{4\pi\mu} \left[ -\frac{\lambda + 3\mu}{\lambda + 2\mu} \ln r \cdot \mathbf{I} + \frac{\lambda + \mu}{\lambda + 2\mu} \frac{\boldsymbol{r}\otimes \boldsymbol{r}}{r^2} \right],\tag{5}\] where \(\boldsymbol{r}= \boldsymbol{x}- \boldsymbol{y}\) and \(r = \|\boldsymbol{r}\|\). This kernel defines the single-layer potential (SLP): \[\label{eq:SLP} S[\boldsymbol{\phi}](\boldsymbol{x}) = \int_\Gamma \mathbf{G}(\boldsymbol{x}, \boldsymbol{y}) \boldsymbol{\phi}(\boldsymbol{y}) \, \mathrm{d}s_{\boldsymbol{y}}.\tag{6}\] The corresponding traction fundamental solution \(\mathbf{D}(\boldsymbol{x},\boldsymbol{y})\), derived via the traction operator \(\mathcal{T}_{\boldsymbol{y}}\) acting on \(\mathbf{G}\), is \[\label{eq:fundamental95solution95traction} \begin{align} \mathbf{D}(\boldsymbol{x},\boldsymbol{y})&=\frac{\mu}{\lambda+2\mu}\left[\frac{\boldsymbol{n}(\boldsymbol{y})\cdot\boldsymbol{r}}{2\pi r^2}\mathbf{I}+\frac{\boldsymbol{\tau}(\boldsymbol{y})\cdot\boldsymbol{r}}{2\pi r^2}\mathbf{L}\right] +\frac{\lambda+\mu}{\lambda+2\mu}\frac{(\boldsymbol{n}(\boldsymbol{y})\cdot\boldsymbol{r})\boldsymbol{r}\otimes\boldsymbol{r}}{\pi r^4}, \end{align}\tag{7}\] where \(\mathbf{L}\triangleq\begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix}\). The double-layer potential (DLP) is subsequently defined by \[\label{eq:DLP} D[\boldsymbol{\phi}](\boldsymbol{x})=\int_\Gamma\mathbf{D}(\boldsymbol{x},\boldsymbol{y})\boldsymbol{\phi}(\boldsymbol{y})\mathrm{d}s_{\boldsymbol{y}}.\tag{8}\]

The single and double-layer potentials satisfy Lamé equations in \(\Omega\) and \(\mathbb{R}^2 \setminus \overline{\Omega}\), and their limiting values on the boundary exhibit the following well-known jump behaviors.

Based on these potentials, we formulate the second-kind boundary integral equations. Throughout, we focus on the interior Dirichlet problem (IDP) for the Lamé system. Using the double-layer potential representation \(\boldsymbol{u}(\boldsymbol{x})=D[\boldsymbol{\phi}](\boldsymbol{x})\) and the jump relation from Lemma 1, we obtain \[~\label{eq:bie95original} \left(-\frac{\theta(\boldsymbol{x})}{2\pi}\mathcal{I}+ \mathcal{K}\right)\boldsymbol{\phi}(\boldsymbol{x}) = \boldsymbol{f}(\boldsymbol{x}), \quad \boldsymbol{x}\in \Gamma.\tag{9}\]

2.2 Well-posedness analysis↩︎

The well-posedness of the boundary integral equation 9 is determined by the spectral properties of the associated double-layer potential operator \(\mathcal{K}\), also known as the Neumann-Poincaré (NP) operator, within an appropriate function space. For domains with corners, the analysis is non-trivial and requires specialized techniques beyond classical theories. Taking the Laplace equation as an example, the operator \(\mathcal{K}\) is compact in \(L^2(\Gamma)\) on a \(C^{1,\alpha}\) boundary, where the Fredholm alternative readily ensures the existence and uniqueness of solutions. However, when the boundary \(\Gamma\) includes corners, the NP operator exhibits a continuous spectrum, and classical Fredholm theory requires local regularization techniques [12], [13]. The situation becomes even more subtle in the Lamé system, since even on smooth boundaries, the NP operator \(\mathcal{K}\) is non-compact and admits the decomposition: \[\label{eq:decompNP95lame} \mathcal{K}[\boldsymbol{\phi}](\boldsymbol{x}) = \int_{\Gamma} \left(\mathcal{T}_{\boldsymbol{y}} \mathbf{G}(\boldsymbol{x},\boldsymbol{y})\right)^{\top} \boldsymbol{\phi}(\boldsymbol{y}) \mathrm{d} s_{\boldsymbol{y}} :=\mathcal{K}_{0}[\boldsymbol{\phi}](\boldsymbol{x})+\mathcal{K}_{1}[\boldsymbol{\phi}](\boldsymbol{x}), \quad \boldsymbol{x}\in \Gamma,\tag{10}\] where the kernels of \(\mathcal{K}_{0}\) and \(\mathcal{K}_{1}\) are given by \[\begin{align} K_{0}(\boldsymbol{x},\boldsymbol{y}) &= \frac{\mu}{\lambda+2 \mu} \frac{\boldsymbol{\tau}(\boldsymbol{y}) \cdot \boldsymbol{r}}{2 \pi r^{2}} \mathbf{L}, \\ K_{1}(\boldsymbol{x},\boldsymbol{y}) &= \frac{\boldsymbol{n}(\boldsymbol{y}) \cdot \boldsymbol{r}}{2 \pi r^{2}} \mathbf{I} -\frac{\lambda+\mu}{\lambda+2 \mu}\left(-\frac{(\boldsymbol{n}(\boldsymbol{y}) \cdot \boldsymbol{r}) \boldsymbol{r}\otimes \boldsymbol{r}}{\pi r^{4}}+\frac{\boldsymbol{n}(\boldsymbol{y}) \cdot \boldsymbol{r}}{2 \pi r^{2}} \mathbf{I}\right). \end{align}\] Here, \(K_1\) is only weakly singular, whereas \(K_0\) contains the Cauchy-type singularity responsible for non-compactness. More precisely, using Stokes’ theorem and the Hilbert transform [28], it is known that \(\mathcal{K}_0^2-k_0^2 \mathcal{I}\) is compact in \(H^{1/2}(\Gamma)^2\) with \(k_0=\frac{\mu}{2(\lambda+2\mu)}\). Although specific regularization techniques enable well-posedness for smooth domains, they become insufficient when the boundary contains corners.

To address this challenge, we adopt a framework developed for elliptic boundary value problems in nonsmooth domains. This theory, initiated by Kondratiev [29] and significantly extended by Costabel, Dauge, and others [30], [31], provides a general operator-theoretic setting for boundary integral formulations on domains with corners. The key component of this framework involves weighted Sobolev spaces, specifically designed to capture the singular behavior of solutions near corners.

Equipped with this functional framework, we can now state the main result of this section, which connects the well-posedness of the operator to the singular exponents derived in Section 4.

We now arrive at the well-posedness of the boundary integral equation for the interior Dirichlet problem.

The sequence of results in this subsection provides a theoretical framework for the analysis of the Lamé system on 2D domains with corners. It confirms that by choosing an appropriate weighted Sobolev space which avoids the corner exponents and the integers, the boundary integral formulation of the interior Dirichlet problem is well-posed, guaranteeing a unique solution.

3 Mathematical Preliminaries↩︎

3.1 Three types of singular integrals↩︎

In this subsection, we compute three Mellin-type singular integrals, referred to as the Cauchy, Laplace, and Stokes cases, respectively. These integrals play a crucial role in the calculations presented in Section 4.

To begin with, we construct two contour paths, as shown in Figure 2. Unless indicated otherwise, we fix the branch cut of the logarithm along the positive real axis, i.e., \(\ln s = \ln|s| + i\,\arg(s)\) with \(\arg(s)\in[0,2\pi)\). Consider a keyhole contour \(\Gamma_{A}=\Gamma_1 \cup \Gamma_2 \cup \Gamma_R \cup \Gamma_r\), as shown in Figure 2 (a), where \(\theta_0=\arcsin(r/R)\), \(r<R\), and \[\begin{align} \Gamma_1 & =\left\{s+i r, 0 \leq s \leq \sqrt{R^2-r^2}\right\}, &\Gamma_2 & = \left\{s-i r, 0 \leq s \leq \sqrt{R^2-r^2}\right\}, \\ \Gamma_R & =\left\{R e^{i \theta}, \theta_0<\theta<2 \pi-\theta_0\right\}, &\Gamma_{r} & =\left\{-r e^{i \theta},-\frac{\pi}{2}<\theta<\frac{\pi}{2}\right\}. \end{align}\] Given \(t>0\), we modify the contours \(\Gamma_1\) and \(\Gamma_2\) by making small semi-circular detours \(\gamma_{\varepsilon}^1\) and \(\gamma_{\varepsilon}^2\) around the point \(s=t\): \[\begin{align} \gamma_{\varepsilon}^1 & =\left\{t + ir + \varepsilon e^{i \theta}, 0<\theta<\pi \right\}, \\ \gamma_{\varepsilon}^2 & =\left\{t - ir + \varepsilon e^{i \theta}, -\pi<\theta<0 \right\}. \end{align}\] Then we construct \(\Gamma_{B}=\Gamma_1^{\varepsilon} \cup \Gamma_2^{\varepsilon} \cup \Gamma_R \cup \Gamma_r \cup \gamma_{\varepsilon}^1 \cup \gamma_{\varepsilon}^2\), as shown in Figure 2 (b).

3.2 Roots of sinc function in the complex domain↩︎

In order to analyze the exponents of the density functions near corners, we need to consider the transcendental equation \[\label{eq:trans2} G(\omega,\alpha)=\operatorname{sinc}(\omega)+\mathrm{K} \operatorname{sinc}(\alpha)=0,\tag{11}\] where \(\omega\in\mathbb{C}\), \(\alpha\in(0,2\pi)\), and \(\mathrm{K}=\pm 1\) or \(\pm \frac{\lambda+\mu}{\lambda+3\mu}\).

Since \(\operatorname{sinc}(\alpha)\) is real-valued, any solution \(\omega\) of 11 must satisfy \(\Im(\operatorname{sinc}(\omega))=0\). Writing \(\omega=x+iy\), we obtain \[\begin{align} &\Re(\operatorname{sinc}(\omega))=\frac{x\sin{x}\cosh{y}+y\cos{x}\sinh{y}}{x^2+y^2},\\ &\Im(\operatorname{sinc}(\omega))=\frac{x\cos{x}\sinh{y}-y\sin{x}\cosh{y}}{x^2+y^2}. \end{align}\] Thus \(\operatorname{sinc}(\omega)\) is purely real if and only if \[\label{eq:im95sinc95zero} x \cos x \sinh y - y \sin x \cosh y = 0.\tag{12}\] For \(y \neq 0\) and \(x \neq 0\), it can be rearranged into a more compact form \[\label{eq:curve95equation} \frac{\tan x}{x} = \frac{\tanh y}{y}.\tag{13}\] Together with the real and imaginary axes, this relation characterizes the entire level set \(\{\omega\in\mathbb{C}:\Im(\operatorname{sinc}(\omega))=0\}\). As illustrated in Figure 3 (a), the black curves represent the contours where \(\Im(\operatorname{sinc}(\omega))=0\), with the top and front views shown in Figures 3 (b) and 3 (c), respectively.

Using Definition 2, we can precisely characterize the properties of \(\operatorname{sinc}(\omega)\) along the curve 13 , as summarized in Lemma 5.

Note that at \(\alpha = \pi\) we have \(\operatorname{sinc}(\alpha) = 0\), so Eq. 11 reduces to \(\operatorname{sinc}(\omega) = 0\), whose real solutions are precisely \(\omega = n\pi\), \(n\in\mathbb{Z}\). We use these points as anchors when continuing the solution branches in \(\alpha\). Now we analyze how a specific root \(\omega\) travels along these curves as \(\alpha\) varies. Along the curves where \(\Im(\operatorname{sinc}(\omega))=0\), the original transcendental equation reduces to finding the roots of a real-valued equation. To characterize the distribution of these roots and their dependence on the parameter \(\alpha\), we explore the properties of the implicit function \(\omega(\alpha)\) defined by \(G(\omega, \alpha) = \operatorname{sinc}(\omega) + \mathrm{K}\operatorname{sinc}(\alpha) = 0\). The implicit function theorem guarantees that for a given point \((\omega_0, \alpha_0)\) satisfying \(G(\omega_0, \alpha_0)=0\), we can find a unique analytic function \(\omega(\alpha)\) in a neighborhood of \(\alpha_0\) as long as \(\partial G/\partial \omega \neq 0\). The “problematic” points, known as branch points, occur precisely where the following condition holds: \[\frac{\partial G}{\partial \omega} = \frac{d}{d\omega}\operatorname{sinc}(\omega) = \frac{\omega\cos(\omega) - \sin(\omega)}{\omega^2} = 0.\] This is equivalent to the equation \(\tan(\omega)=\omega\), which means branch points in the \(\omega\)-plane occur exactly at \(\omega = \lambda_j\) (as defined in Definition 2). This naturally leads to the definition of branch points in the physical parameter \(\alpha\).

To gain insight into the global behavior of these solution branches, we visualize the trajectories for the representative case \(\mathrm{K}=1\) in Figure 4. The behavior of \(\operatorname{sinc}(\alpha)\) over the interval \((0, 2\pi)\) is illustrated in Figure 4 (a). As the parameter \(\alpha\) sweeps from \(0\) to \(2\pi\), the root \(\omega(\alpha)\) anchored at \(\omega(\pi) = n\pi\) traces a specific path in the complex plane. These trajectories exhibit distinct topological patterns based on the parity of \(n\): for even indices \(n=2m\) (\(m\ge 1\)), the roots follow the paths indicated by the red arrows in Figure 4 (b). Conversely, for odd indices \(n=2m+1\) (\(m\ge 0\)), the roots follow the trajectories marked by black arrows.

Formalizing the existence and analyticity of these observed branches, we present the following theorem.

4 Integral equations near corners↩︎

With all the analytic tools prepared above, we are ready to capture the singularities of Eq. 9 on cornered boundaries. By localization and a partition-of-unity technique, the boundary integral operator near each corner reduces to the corresponding operator on a wedge. Suppose the wedge \(\Gamma\subset\mathbb{R}^2\) is composed of two edges \(OA\) and \(OB\), parameterized by \[\label{eq:WedgePara} \gamma(t)=\left\{ \begin{array}{ll} (t, 0), & \;0 \leq t \leq 1, \\ (-t\cos\theta, -t\sin\theta), & \;-1 \leq t \leq 0. \\ \end{array}\right.\tag{15}\] By extending the sides of the wedge to infinity, we divide \(\mathbb{R}^2\) into two open sets \(D_1\) and \(D_2\), as shown in Figure 5.

We will identify a countable set of exponents \(\{z_{n,j}\}_{n\ge1,\,j=1,\dots,4}\subset\mathbb{C}\) together with associated coefficient vectors \(\mathbf{p}^{n,j}=[p_1^{n,j},p_2^{n,j},p_3^{n,j},p_4^{n,j}]^\top\), depending on \(\theta\), such that, for any fixed \((n,j)\), if the density function \(\boldsymbol{\phi}\) is defined by \[\boldsymbol{\phi}(t)= \begin{cases} [p_1^{n, j}, p_2^{n, j}]^\top \cdot t^{z_{n, j}}, &0\leq t\leq 1,\\ [p_3^{n, j}, p_4^{n, j}]^\top \cdot (-t)^{z_{n, j}}, &-1\leq t\leq 0, \end{cases}\] then the left-hand side of Eq. 9 is smooth near \(t=0\).

We will also prove the converse. Specifically, let \(N\) be a positive integer, and suppose further that \(\beta_{n,j} \in \mathbb{C}\), for \(j=1,2,3,4\), and \(n=0,1,2, \ldots, N-1\), with \(\boldsymbol{f}\) satisfying \[\boldsymbol{f}(t)= \begin{cases} \sum\limits_{n=0}^{N-1}[\beta_{n, 1} , \beta_{n,2} ]^\top \cdot t^n, &0\leq t\leq 1,\\ \sum\limits_{n=0}^{N-1}[\beta_{n, 3} , \beta_{n,4} ]^\top \cdot (-t)^n, &-1\leq t\leq 0. \end{cases}\] Then, for all but countably many values of \(0 < \theta < 2\pi\), there exist unique coefficients \(\alpha_{n,j} \in \mathbb{C}\), for \(j=1,2,3,4\) and \(n=1, \ldots N\), such that the density function defined by \[\boldsymbol{\phi}(t)=\begin{cases} \sum\limits_{n=1}^N\sum\limits_{j=1}^4 \alpha_{n,j}[p_1^{n, j}, p_2^{n, j}]^\top \cdot t^{z_{n, j}}, &0\leq t\leq 1,\\ \sum\limits_{n=1}^N\sum\limits_{j=1}^4 \alpha_{n,j}[p_3^{n, j}, p_4^{n, j}]^\top \cdot (-t)^{z_{n, j}}, &-1\leq t\leq 0, \end{cases}\] satisfies Eq. 9 with error \(O\left(|t|^{N}\right)\).

4.1 Boundary integral equations on the wedge↩︎

To derive the boundary integral equations, we place the point \(x\) on the edges \(OA\) and \(OB\) respectively, from which we can derive the corresponding boundary integral operators for \(x\in OA\) and \(x\in OB\). Assume that \(\Re(z)>-1\) and \(z\notin\mathbb{Z}\), and set \[\label{eq:density} \boldsymbol{\phi}(t)= \begin{cases} [p_1,p_2]^\top \cdot t^z, &0\leq t\leq 1,\\ [p_3,p_4]^\top \cdot (-t)^z, &-1\leq t\leq 0. \end{cases}\tag{16}\] To maintain brevity, denote \[c_1 \triangleq \frac{\mu}{2(\lambda+2\mu)}, \quad c_2 \triangleq \frac{\lambda+\mu}{\lambda+2\mu},\] \[I_{OA} \triangleq \int_{OA}\mathbf{D}(\boldsymbol{x},\boldsymbol{y})\boldsymbol{\phi}(\boldsymbol{y})\mathrm{d}s_{\boldsymbol{y}}, \quad I_{OB} \triangleq \int_{OB}\mathbf{D}(\boldsymbol{x},\boldsymbol{y})\boldsymbol{\phi}(\boldsymbol{y})\mathrm{d}s_{\boldsymbol{y}}.\] For \(x\in OA\), we have \(x=(t,0)\) where \(0\leq t\leq 1\), then \[\begin{align} I_{OA} &=\frac{c_1}{\pi}\int_0^1 \frac{s^z}{t-s}\begin{bmatrix}p_2\\-p_1\end{bmatrix}\mathrm{d}s,\\ I_{OB} &= \frac{c_1}{\pi}\left\{\int_{-1}^0\frac{-t\sin\theta\cdot(-s)^z}{ t^2+s^2+2ts\cos\theta}\begin{bmatrix}p_3\\p_4\end{bmatrix}\mathrm{d}s + \int_{-1}^0\frac{(-t\cos\theta-s)\cdot(-s)^z}{t^2+s^2+2ts\cos\theta}\begin{bmatrix}p_4\\-p_3\end{bmatrix}\mathrm{d}s\right\} \\ & + \frac{c_2}{\pi} \int_{-1}^0\frac{-t\sin\theta\cdot(-s)^z}{(t^2+s^2+2ts\cos\theta)^2} \begin{bmatrix} (t+s\cos\theta)^2 & (t+s\cos\theta)s\sin\theta \\ (t+s\cos\theta)s\sin\theta &s^2\sin^2\theta \end{bmatrix} \begin{bmatrix} p_3\\p_4 \end{bmatrix}\mathrm{d}s, \end{align}\] which further gives the boundary double-layer operator on \(OA\): \[\begin{align} \mathcal{D}_{OA}(t) = c_1 \left\{ - I_0(z,t)\cdot \begin{bmatrix}p_2\\-p_1\end{bmatrix} + I_3 \begin{bmatrix}p_3\\p_4\end{bmatrix} + I_4 \begin{bmatrix}p_4\\-p_3\end{bmatrix}\right\} + c_2 \begin{bmatrix} I_{11} & I_{12}\\ I_{21} & I_{22} \end{bmatrix} \begin{bmatrix} p_3\\p_4 \end{bmatrix}, \end{align}\] where \(0\leq t\leq 1\) and \[\begin{align} &I_3 \triangleq - t\sin\theta \cdot I_1(z,\theta,t),\\ &I_4 \triangleq -t\cos\theta\cdot I_1(z,\theta,t)+I_1(z+1,\theta,t),\\ &I_{11} \triangleq -t^3\sin\theta\cdot I_2(z,\theta,t)+2t^2\sin\theta\cos\theta\cdot I_2(z+1,\theta,t)-t\sin\theta\cos^2\theta\cdot I_2(z+2,\theta,t),\\ &I_{12} \triangleq I_{21} = t^2\sin^2\theta\cdot I_2(z+1,\theta,t) - t\cos\theta\sin^2\theta\cdot I_2(z+2,\theta,t),\\ &I_{22} \triangleq -t\sin^3\theta\cdot I_2(z+2,\theta,t). \end{align}\] For \(x\in OB\), we have \(x=-t(\cos\theta,\sin\theta)\) where \(-1\leq t\leq 0\), then \[\begin{align} I_{OB} &= \frac{c_1}{\pi}\int_{-1}^0 \frac{(-s)^z}{t-s} \begin{bmatrix}p_4\\-p_3\end{bmatrix}\mathrm{d}s, \\ I_{OA} &= \frac{c_1}{\pi}\left\{\int_0^1\frac{t\sin\theta\cdot s^z}{ t^2+s^2+2ts\cos\theta}\begin{bmatrix}p_1\\p_2\end{bmatrix}\mathrm{d}s + \int_0^1\frac{(-t\cos\theta-s)\cdot s^z}{t^2+s^2+2ts\cos\theta}\begin{bmatrix}p_2\\-p_1\end{bmatrix}\mathrm{d}s\right\} \\ & + \frac{c_2}{\pi} \int_0^1\frac{t\sin\theta\cdot s^z}{(t^2+s^2+2ts\cos\theta)^2} \begin{bmatrix} (t\cos\theta+s)^2 & (t\cos\theta+s)t\sin\theta \\ (t\cos\theta+s)t\sin\theta &t^2\sin^2\theta \end{bmatrix} \begin{bmatrix} p_1\\p_2 \end{bmatrix}\mathrm{d}s. \end{align}\] Replace \(t\) by \(-t\), \(0\leq t\leq 1\), then we can obtain the boundary double-layer operator on \(OB\): \[\begin{align} \mathcal{D}_{OB}(-t) =& c_1\left\{I_0(z,t)\begin{bmatrix} p_4\\-p_3 \end{bmatrix}+ I_3 \begin{bmatrix}p_1\\p_2\end{bmatrix} - I_4 \begin{bmatrix}p_2\\-p_1\end{bmatrix}\right\} + c_2 \begin{bmatrix} I_{11}^* & I_{12}^*\\ I_{21}^* & I_{22}^* \end{bmatrix} \begin{bmatrix} p_1\\p_2 \end{bmatrix}, \end{align}\] where \(0\leq t\leq 1\) and \[\begin{align} I_{11}^* &\triangleq -t^3\sin\theta\cos^2\theta\cdot I_2(z,\theta,t)+2t^2\sin\theta\cos\theta\cdot I_2(z+1,\theta,t)-t\sin\theta\cdot I_2(z+2,\theta,t),\\ I_{12}^* &\triangleq I_{21}^* = -t^3\cos\theta\sin^2\theta\cdot I_2(z,\theta,t) + t^2\sin^2\theta\cdot I_2(z+1,\theta,t),\\ I_{22}^* &\triangleq -t^3\sin^3\theta\cdot I_2(z,\theta,t). \end{align}\]

For the sake of brevity, denote the left-hand side of Eq. 9 by \(\boldsymbol{h}(\boldsymbol{x})\). Then we derive \[\nonumber \boldsymbol{h}(\boldsymbol{x}(t)) = -\frac{1}{2}\boldsymbol{\phi}(\boldsymbol{x}(t)) + \mathcal{D}[\boldsymbol{\phi}](t),\;-1\leq t\leq 1.\] Replacing \(t\) by \(-t\) on the interval \([-1,0]\), we obtain a linear system of four equations involving the coefficient vector \(\mathbf{p}=[p_1, p_2, p_3, p_4]^\top\): \[\begin{cases} \boldsymbol{h}(\boldsymbol{x}(t)) = -\frac{1}{2}[p_1,p_2]^\top\cdot t^z + \mathcal{D}_{OA}(t), \\ \boldsymbol{h}(\boldsymbol{x}(-t)) = -\frac{1}{2}[p_3,p_4]^\top\cdot t^z + \mathcal{D}_{OB}(-t), \end{cases} 0\leq t\leq 1.\] More specifically, for \(0\leq t\leq 1\), we have \[\nonumber \tilde{\boldsymbol{h}}(t)\triangleq \begin{bmatrix} h_1(t)\\h_2(t)\\h_1(-t)\\h_2(-t) \end{bmatrix}= -\frac{1}{2} \begin{bmatrix} p_1\\p_2\\p_3\\p_4 \end{bmatrix} t^z + \mathbf{G}_{\mathrm{DL}}(z,\theta,t) \begin{bmatrix} p_1\\p_2\\p_3\\p_4 \end{bmatrix},\] where \[\mathbf{G}_{\mathrm{DL}}(z,\theta,t)\triangleq \begin{bmatrix} 0&-c_1 I_0&c_1 I_3+c_2 I_{11}&c_1 I_4+c_2 I_{12}\\ c_1 I_0&0&-c_1 I_4+c_2 I_{12}&c_1 I_3+c_2 I_{22}\\ c_1 I_3+c_2 I_{11}^*&-c_1 I_4+c_2 I_{12}^*&0&c_1 I_0\\ c_1 I_4+c_2 I_{12}^*&c_1 I_3+c_2 I_{22}^*&-c_1 I_0&0 \end{bmatrix}.\] Upon inserting the explicit formulas of \(I_0,I_1,I_2,I_3,I_4,I_{11},I_{12},I_{22},I_{11}^*,I_{12}^*,I_{22}^*\) via Lemmas 2, 3 and 4, we obtain \[\tilde{\boldsymbol{h}}(t) = \mathbf{A}(z,\theta)\mathbf{p} \cdot t^z + \sum_{k=0}^{\infty}\mathbf{F}(k,z,\theta) \mathbf{p} \cdot t^k,\] where the explicit closed forms of the matrix \(\mathbf{A}(z,\theta)\) and \(\mathbf{F}(k,z,\theta)\) are collected in Appendix 7. The above results are summarized in the following theorem.

4.2 Existence of the implicit function \(z(\theta)\)↩︎

Now we analyze the implicit functions \(z(\theta)\) defined by the transcendental equation ?? , which is closely connected to Eq. 11 , as shown in the following lemma.

Consider the case \(H_1(z, \theta)\). Combining Lemma 6 with Theorem 3, we can establish the existence of an implicit function \(z_{n,1}(\theta)\) such that \(H_1\bigl(z_{n,1}(\theta),\theta\bigr)=0\).

4.3 Invertibility of mapping matrix↩︎

Suppose \(\mathbf{B}(\theta)\) is the mapping matrix transforming the coefficients of the truncated basis functions \[(p_1^{n,j} t^{z_{n,j}}, p_2^{n,j} t^{z_{n,j}}, p_3^{n,j} t^{z_{n,j}}, p_4^{n,j} t^{z_{n,j}}),\;n=1,\ldots,N,\;j=1,2,3,4\] to Taylor expansion coefficients of the vector \(\tilde{\boldsymbol{h}}(t)\). The explicit form of \(\mathbf{B}(\theta)\) is established in the following theorem.

The following two lemmas establish the invertibility of \(\mathbf{B}(\theta)\) for almost every angle.

Now we obtain the final result of this section, which shows that for most opening angles \(\theta\), the density function \(\boldsymbol{\phi}\) can be constructed to match arbitrary polynomial boundary data up to order \(N-1\).

4.4 Hybrid basis: limitations on sharp corners↩︎

It seems natural to use the functions obtained from Theorem 6 as a basis for the numerical solutions of BIEs on cornered domains. However, we will show the analytic basis functions are unsuitable for direct use in numerical algorithms. For that purpose, we design a hybrid semi-analytic boundary element method (called HybridBEM) on an equilateral triangle to explore the possibilities of using corner singularities as basis functions. Specifically, on each corner panel \(\Gamma_c\) the density is expanded in the first \(M\) singular powers \(\{t^{z_k}\}_{k=1}^M\) with graded collocation, while the remaining boundary \(\Gamma_s\) uses \(p\)-th order Legendre expansion. Panel interactions are assembled by factoring the kernel-density product into a “singular weight \(\times\) smooth factor” and evaluating the resulting panel moments semi-analytically. Particularly, cross-corner terms follow the explicit mapping series in subsection 4.3.

Figure 6 tests the performance as we vary the number of corner modes \(M\). As \(M\) increases, the number of unknowns grows linearly, yet the conditioning deteriorates: \(\kappa_1(A)\) rises from \(\mathcal{O}(10^2)\) at \(M=1\) to \(\mathcal{O}(10^{19})\) at \(M=20\). The accuracy does not benefit from the extra corner modes, and even in the best case, the relative error levels off near \(O(10^{-5})\), which is already worse than the uniform-panel baseline \(O(10^{-6})\). This breakdown is algebraic rather than quadrature-driven. In particular, sampling the functions \(\{t^{z_k}\}\) on a finite corner panel produces columns with a Vandermonde-like structure. When several real parts \(\Re z_k\) are close, or when a conjugate pair has a very small imaginary part \(\Im z\), these columns become nearly linearly dependent, and the local Gram matrices are severely ill-conditioned, which explains the poor stability and lack of accuracy improvement. To overcome these issues, we introduce a singularity guided Nyström method in the next section.

5 Singularity guided Nyström method↩︎

Solving the elastic integral equation 9 in domains with corners requires specialized techniques to handle both kernel and geometric singularities effectively. Based on the analytical results above, we propose a singularity guided high-order Nyström (SGN) scheme. Specifically, corner singular powers are extracted from the transcendental equation \(H(z,\theta) = 0\) via an analytic root-finding algorithm in the complex plane [34] and then used to guide panelization and accuracy control, while the unknown density is still discretized at Gauss-Legendre (GL) nodes. Near or self interactions are evaluated by kernel splitting and a forward-stable Legendre recurrence [19], [21], [35], effectively addressing Cauchy-type singularities. This yields a square and well-conditioned linear system, so a dense linear solver suffices to obtain stable results. Details of the SGN scheme are given as follows.

5.1 Algorithm design↩︎

We parameterize each boundary arc by arclength \(t\in[0, L]\) originating from the corner at \(t=0\). For the opening angle \(\theta\), let \(z_1, \ldots, z_n\) be the leading singularity exponents with \(\Re z > 0\), computed via the contour integration method [34]. To resolve the corner singularities, we construct an adaptive panel mesh by recursive bisection. On each panel \(I\), we define the multi-exponent Legendre-tail indicator at degree \(p\) by \[\label{eq:legendre95indicator} \eta(I;p):=\max_{1\leq k\leq n}\Bigg(\sum_{j=p}^{2p-1}\big|\alpha_j(z_k;I)\big|^2\Bigg)^{1/2},\tag{17}\] where \(\alpha_j(z_k;I)\) is the \(j\)-th Legendre coefficient of the singular function \(t^{z_k}\) restricted to \(I\). A panel is accepted if \[\label{acceptrule} \eta(I;p)\le\varepsilon_{\mathrm{pan}},\tag{18}\] where \(\varepsilon_{\mathrm{pan}}\) is a given error tolerance. Otherwise, it is bisected. Panels away from corners are kept uniform with the same degree \(p\). In this way, the mesh clusters automatically wherever the computed singular family \(\{t^{z_k}\}\) demands resolution.

Now we consider the quadrature rules for the boundary integral: \[\label{eq:double95layer} \mathcal{K}[\boldsymbol{\phi}](\boldsymbol{x}) = \int \mathbf{D}\left(\boldsymbol{x}, \boldsymbol{y}\right) \boldsymbol{\phi}(\boldsymbol{y}) \mathrm{d}s_{\boldsymbol{y}},\tag{19}\] where \(\mathbf{D}\left(\boldsymbol{x}, \boldsymbol{y}\right)\) is given in 7 . Far-field interactions are computed via the standard \(p\)-point Gauss-Legendre rule. For near-field and self interactions, we utilize the kernel splitting technique: \[\mathbf{D}(t,s) = \frac{1}{t-s} \mathbf{D}_1(t,s) + \mathbf{D}_2(t,s),\] where \(\mathbf{D}_1, \mathbf{D}_2\) are analytic on \([-1,1]^2\): \[\begin{align} \mathbf{D}_1(t, s) &= \frac{\mu}{\lambda+2\mu}\frac{(t-s)\,(\boldsymbol{\tau}(\boldsymbol{x})\cdot\boldsymbol{r})}{2\pi r^2}\mathbf{L},\\ \mathbf{D}_2(t, s) &= \frac{\mu}{\lambda+2\mu}\frac{(\boldsymbol{\tau}(\boldsymbol{y})-\boldsymbol{\tau}(\boldsymbol{x}))\cdot\boldsymbol{r}}{2\pi r^2}\mathbf{L}+\frac{\mu}{\lambda+2\mu}\frac{\boldsymbol{n}(\boldsymbol{y})\cdot\boldsymbol{r}}{2\pi r^2}\mathbf{I} +\frac{\lambda+\mu}{\lambda+2\mu}\frac{(\boldsymbol{n}(\boldsymbol{y})\cdot\boldsymbol{r})\boldsymbol{r}\otimes\boldsymbol{r}}{\pi r^4}. \end{align}\] Moreover, let \(J(s) = |\gamma'(s)|\) be the Jacobian of the panel parameterization at the variable \(s\), i.e., \(d s_{\boldsymbol{y}}=J(s)\,\mathrm{d}s\). Interpolating \(\mathbf{D}_j\,\boldsymbol{\phi}\,J\), \(j=1,2\) by Legendre polynomials \[(\mathbf{D}_1\,\boldsymbol{\phi}\,J)(t,s)\approx\sum_{j=0}^{p-1}\boldsymbol{d}_j(t)\,P_j(s),\qquad (\mathbf{D}_2\,\boldsymbol{\phi}\,J)(t,s)\approx\sum_{j=0}^{p-1}\boldsymbol{f}_j(t)\,P_j(s),\] leads to panel moments \[C_j(t):=\text{p.v.}\int_{-1}^{1}\frac{P_j(s)}{t-s}\,\mathrm{d}s,\qquad I_j:=\int_{-1}^{1}P_j(s)\,\mathrm{d}s .\] Since the Gauss-Legendre nodes never include the endpoints, the Cauchy moments \(C_j\) can be evaluated via the three-term recurrence \[\bigg\{ \begin{align} &C_0(t)=-\ln\Big|\frac{t-1}{t+1}\Big|,\quad C_1(t)=-2+t\,C_0(t), \\ &j\,C_j=(2j-1)\,t\,C_{j-1}-(j-1)\,C_{j-2}, \qquad j\ge2, \end{align} \bigg.\] for any real \(t\neq\pm1\). Thus, near-field or self interactions can be approximated by \(\sum_{j=0}^{p-1}\big[C_j(t)\) \(\boldsymbol{d}_j(t)+I_j\,\boldsymbol{f}_j(t)\big]\). Enforcing the equation at all Gauss-Legendre nodes produces a square Nyström system, which is then solved by a direct solver.

5.2 Error analysis↩︎

Having detailed the algorithmic framework, we now analyze the error in the weighted Sobolev space \(X := H_\nu^s(\Gamma)^2\) with \(s\in(1/2,3/2)\) and \(s<\nu<s+1/2\). Assuming \(s-\nu\notin\{\Re z_{n,j}\} \cup \mathbb{Z}\), Theorem 2 guarantees that the interior double-layer operator \[\mathcal{A}:=\Big(-\frac{\theta(\boldsymbol{x})}{2\pi}\mathcal{I}+\mathcal{K}\Big):X\to X\] is an isomorphism with \(C_{\mathcal{A}}:=\|\mathcal{A}^{-1}\|_{\mathcal{L}(X)}<\infty\), where \(\|\cdot\|_{\mathcal{L}(X)}\) denotes the operator norm from \(X\) to \(X\). Throughout this subsection, \(C>0\) denotes a generic constant independent of \(h\) and \(\varepsilon_{\rm pan}\).

Combining the stability of invertible operators, we derive the main error bound.

5.3 Numerical experiments↩︎

We now test the performance of the SGN scheme through several numerical experiments. Unless stated otherwise, we set the Lamé parameters \((\lambda,\mu)=(1,2)\), polynomial degree \(p=16\), and length factor \(\sigma=0.1\). Boundary data are synthesized using a few exterior sources placed away from the domain. The relative error is measured at interior target points \(\boldsymbol{x}_m\) by \[E=\Bigg(\frac{\sum_m\|\boldsymbol{u}_{\text{num}}(\boldsymbol{x}_m)-\boldsymbol{u}_{\text{ex}}(\boldsymbol{x}_m)\|^2}{\sum_m\|\boldsymbol{u}_{\text{ex}}(\boldsymbol{x}_m)\|^2}\Bigg)^{1/2},\] where \(\boldsymbol{u}_{\text{num}}\) and \(\boldsymbol{u}_{\text{ex}}\) denote the numerical and exact displacement fields, and \(\|\cdot\|\) is the Euclidean norm in \(\mathbb{R}^2\). All the experiments are implemented in Matlab on a personal computer equipped with an ARM processor.

5.3.1 Singularity guided panelization↩︎

In this example, we examine how SGN distributes panels near corners and how the error depends on the panel tolerance. We consider three representative geometries: a droplet, a triangle, and an L-shaped domain. For each geometry, we apply SGN varying the tolerance \(\varepsilon_{\rm pan}\). Figure 7 illustrates the resulting graded boundary discretizations and their parameter-space distributions.

For the same three geometries, we vary the threshold exponent \(N\) so that \(\varepsilon_{\rm pan} = 10^{-N}\), and solve the scattering problem by SGN. The results are plotted in Figure 8 (a), which shows the relative error versus \(N\) for each shape in logarithmic scale. As \(\varepsilon_{\text{pan}}\) decreases, the relative error decays almost linearly with the threshold in the regime where panel refinement dominates, which is in agreement with the error estimate given by Corollary 1.

5.3.2 Comparison with a uniform mesh↩︎

Guided by the first example, we choose \(\varepsilon_{\rm pan}\) such that the accuracy stops increasing for each geometry and compare SGN with a uniform-panel Nyström method (UM). In UM, each edge is discretized by uniform panels with the same Gauss–Legendre degree \(p\). Figure 9 shows the corresponding scattering fields for the droplet, triangle, and L-shaped domains, which provides a qualitative figure of the accuracy gains obtained by singularity guided panelization. Table 3 summarizes the number of degrees of freedom (DoFs) and the relative error \(E\) for a fixed target point. Across all shapes, SGN achieves substantially higher accuracy than UM. In particular, for the droplet and triangle, SGN reaches near machine precision, while UM stagnates at \(10^{-8}\)–\(10^{-6}\). For the re-entrant L-shaped domain, SGN still improves the accuracy by several orders, but the error is noticeably higher. We call this the re-entrant angle effect. The detailed reason for this accuracy loss will be studied in the next example.

5.3.3 Re-entrant angle effect↩︎

To rigorously isolate the role of the corner geometry, we consider a family of single-corner droplet domains. Figure 8 (b) shows \(\log_{10}E\) versus the threshold exponent \(N\) for several opening angles \(\theta\). For any fixed \(N\) in the refinement-dominated regime, the error increases monotonically as \(\theta\) moves from convex to re-entrant. Specifically, while the error remains near machine precision for acute and right angles, it rises sharply as \(\theta\) approaches \(2\pi\). This degradation occurs despite the SGN algorithm automatically clustering panels near the singularity, suggesting an intrinsic difficulty in the re-entrant regime that simple mesh refinement cannot easily overcome.

To investigate the underlying cause, we select two representative cases for a detailed comparison: a convex corner with \(\theta=\pi/2\) and a re-entrant corner with \(\theta=3\pi/2\). Both runs use SGN with the four-branch singular sets and an adaptive \(\varepsilon_{\rm pan}\). We also list the dominant corner exponents to explore the differences, as shown in Table 4. The results show that the re-entrant case requires markedly more panels near the corner and still exhibits a larger error. A possible explanation is that the performance gap is driven by complex roots: in the re-entrant case the characteristic set contains \(z=a\pm ib\), which may induce logarithmic oscillations and make geometric clustering less effective. However, we observe that the same pattern also appears for Laplace problems, where all singular exponents are real, suggesting that complex exponents are not the primary cause.

To visualize what happens near the corner, we plot for \(\theta=\pi/2\) the two density components in Figure 10. Even in the convex case, the innermost panels show mild numerical artifacts—slight flattening and tiny ripples—which are consistent with parameter-space crowding. Locally, a wedge of opening \(\theta\) is conformally flattened by the power map \(w=z^{\pi/\theta}\), whose Jacobian scales like \(|dw/dz|\sim r^{\,\pi/\theta-1}\). When \(\theta>\pi\), the Jacobian blows up as \(r\to0\), forcing many parameter nodes to squeeze into a very short arc in physical space. This amplifies near singular quadrature error and the condition number of the discretized system, forcing the error to plateau earlier than expected given the refined clustering.

6 Conclusions and future work↩︎

This paper investigates the two-dimensional Lamé system on domains with corners through boundary integral equations. On the analytical side, a local Mellin analysis on a wedge yields a factorizable characteristic equation for the boundary density’s singular exponents, cleanly separating physical branches and explaining which exponents are compatible with the corresponding problems. In weighted Sobolev spaces \(H_\nu^s(\Gamma)^2\), we derive a computable Fredholm criterion for the interior double-layer formulation by excluding the real parts of these exponents and the integer poles. Moreover, we prove the invertibility of an explicit density-to-Taylor mapping that links singular density to smooth boundary data. On the numerical side, we design a singularity guided high-order Nyström (SGN) scheme: corner exponents drive panel adaptivity through a multi-exponent Legendre-tail indicator, and near/self interactions are evaluated with a Cauchy-Legendre three-term recurrence. We prove an a priori error bound that couples the panel indicator with exponentially small far-field quadrature error. Numerical experiments on cornered geometries show that SGN achieves higher accuracy than a uniform Nyström baseline, identify re-entrant corners as crowding-limited, and demonstrate that hybrid singular-basis expansions are severely ill-conditioned.

Based on the present results, we propose several directions for future research. We intend to apply the SGN framework to time-harmonic and quasi-static elasticity to quantify corner behaviors, as well as to transmission problems involving interfacial cracks. Geometrically, the method shows promise for generalization to three-dimensional axisymmetric models and polyhedral structures, where edge and vertex singularities interact. Finally, we will investigate shape derivatives in cornered domains for optimal design, and incorporate the close evaluation [38] to ensure accuracy near boundaries.

References↩︎