Solving Poisson Problems in Polygonal Domains with Singularity Enriched
  Physics Informed Neural Networks
Physics-Informed Neural Networks (PINNs) are a powerful class of numerical
solvers for partial differential equations, employing deep neural networks with
successful applications across a diverse set of problems. However, their
effectiveness is somewhat diminished when addressing issues involving
singularities, such as point sources or geometric irregularities, where the
approximations they provide often suffer from reduced accuracy due to the
limited regularity of the exact solution. In this work, we investigate PINNs
for solving Poisson equations in polygonal domains with geometric singularities
and mixed boundary conditions. We propose a novel singularity enriched PINN
(SEPINN), by explicitly incorporating the singularity behavior of the analytic
solution, e.g., corner singularity, mixed boundary condition and edge
singularities, into the ansatz space, and present a convergence analysis of the
scheme. We present extensive numerical simulations in two and three-dimensions
to illustrate the efficiency of the method, and also a comparative study with
several existing neural network based approaches.
Tianhao Hu, Bangti Jin, Zhi Zhou
Solving Poisson Problems in Polygonal Domains with Singularity Enriched Physics Informed Neural Networks+

###### Abstract

Physics informed neural networks (PINNs) represent a very powerful class of numerical solvers for partial differential equations using deep neural networks, and have been successfully applied to many diverse problems. However, when applying the method to problems involving singularity, e.g., point sources or geometric singularities, the obtained approximations often have low accuracy, due to limited regularity of the exact solution. In this work, we investigate PINNs for solving Poisson equations in polygonal domains with geometric singularities and mixed boundary conditions. We propose a novel singularity enriched PINN (SEPINN), by explicitly incorporating the singularity behavior of the analytic solution, e.g., corner singularity, mixed boundary condition and edge singularities, into the ansatz space, and present a convergence analysis of the scheme. We present extensive numerical simulations in two and three-dimensions to illustrate the efficiency of the method, and also a comparative study with existing neural network based approaches.

**Key words:** Poisson equation, corner singularity, edge singularity, physics informed neural network, singularity enrichment

## 1 Introduction

Partial differential equations (PDEs) represent a class of mathematical models that occupy a vital role in physics, science and engineering. Many traditional PDE solvers have been developed, e.g., finite difference method, finite element method and finite volume method. These methods have been maturely developed over the past decades and efficient implementations and mathematical guarantees are also available. In the last few years, motivated by the great success in computer vision, speech recognition and natural language processing etc, neural solvers for PDEs using deep neural networks (DNNs) have received much attention [42]. The list of neural solvers includes physics informed neural networks (PINNs) [53], deep Ritz method (DRM) [60], deep Galerkin method [56], weak adversarial network [62] and deep least-squares method [12] etc. Compared with traditional methods, neural PDE solvers have shown very promising results in several direct and inverse problems [36, 22, 17].

In these neural solvers, one employs DNNs as ansatz functions to approximate the solution to the PDE either in strong, weak or Ritz formulations. The existing approximation theory of DNNs [30] indicates that the accuracy of the DNN approximations depends crucially on the Sobolev regularity of the solution (also suitable stability of the mathematical formulation). Thus, these methods might be ineffective or even fail completely when applied to problems with irregular solutions [58, 41], e.g., convection-dominated problems, transport problem, high-frequency wave propagation, problems with geometric singularities (cracks / corner singularity) and singular sources. All these settings lead to either strong directional behavior or solution singularities or highly oscillatory behavior, which are challenging for standard DNNs to approximate effectively.

Thus, there is an imperative need to develop neural solvers for PDEs with nonsmooth solutions. Indeed, several recent efforts have been devoted to addressing the issue, including self-adaptive PINN (SAPINN) [32], failure-informed PINN (FIPINN) [26, 27] and singularity splitting DRM (SSDRM) [31]. SAPINN extends the standard PINN by splitting out the regions with singularities and then setting different weights to compensate the effect of singular regions. FIPINN [27] is inspired by the classical
adaptive FEM, using the PDE residual as the indicator to aid judicious selection of sampling points for training. SSDRM [31] exploits analytic insights into the exact solution, by approximating only the regular part using DNNs whereas extracting the singular part explicitly. These methods have shown remarkable performance for problems with significant singularities. One prime example is point sources, whose solutions involve localized singularities that can be extracted using fundamental solutions.

In this work, we continue this line of research for Poisson problems on polygonal domains, which involve geometric singularities, including corner singularities and mixed boundary condition in the two-dimensional (2D) case, and edge singularity in the three-dimensional (3D) case. This represents an important setting in practical applications that has received enormous attention; see [29, 39, 40, 48] for the solution theory. We shall develop a class of effective neural solvers for Poisson problems with geometric singularities based on the idea of singularity enrichment, building on known analytic insights of the problems, and term the proposed method singularity enriched PINN (SEPINN).

### Problem setting

First, we state the mathematical formulation of the problem. Let \(\Omega\in\mathbb{R}^{d}\) (\(d=2,3\)) be an open, bounded polygonal domain with a boundary \(\partial\Omega\), and \(\Gamma_{D}\) and \(\Gamma_{N}\) be a partition of the boundary \(\partial\Omega\) such that \(\Gamma_{D}\cup\Gamma_{N}=\partial\Omega\) and \(\Gamma_{D}\cap\Gamma_{N}=\emptyset\), with a nonempty \(\Gamma_{D}\) (i.e., Lebesgue measure \(|\Gamma_{D}|\neq 0\)). Let \(n\) denote the unit outward normal vector to the boundary \(\partial\Omega\), and \(\partial_{n}u\) denote taking the outward normal derivative. Given a source \(f\in L^{2}(\Omega)\), consider the following Poisson problem

\[\left\{\begin{aligned} -\Delta u=& f,& \text{in }\Omega,\\ u=& 0,&\text{on }\Gamma_{D},\\ \partial_{n}u=& 0,&\text{on }\Gamma_{N}.\end{aligned}\right. \tag{1.1}\]

We focus on the zero boundary conditions, and nonzero ones can be easily transformed to (1.1) using the trace theorem. Due to the existence of corners, cracks or edges in \(\Omega\), the solution \(u\) of problem (1.1) typically exhibits singularities, even if the source \(f\) is smooth. The presence of singularities in the solution \(u\) severely deteriorates the accuracy of standard numerical methods for constructing approximations, including neural solvers, and more specialized techniques are needed in order to achieve high numerical efficiency. Next we briefly review existing techniques for resolving the singularities.

In the 2D case, there are several classes of traditional numerical solvers based on FEM, including singularity representation based approaches [23, 57, 13], mesh grading [54, 3, 2], generalized FEM [25] and adaptive FEM [28] etc. These methods require different amount of knowledge about the analytic solution. The methods in the first class exploits a singular representation of the solution \(u\) as a linear combination of singular and regular parts [29, 18], and can be further divided into four groups. (i) The singular function method augments singular functions to both trial and test spaces [23, 57]. However, the convergence of the coefficients (a.k.a. stress intensity factors) sometimes is poor [21], which may lead to low accuracy near singularities. (ii) The dual singular function method [10] employs the dual singular function to extract the coefficients as a postprocessing strategy of FEM, which can also achieve the theoretical rate in practical computation. (iii) The singularity splitting method [13] splits the singular part from the solution \(u\) and approximates the smooth part with the Galerkin FEM, and enjoys \(H^{1}(\Omega)\) and \(L^{2}(\Omega)\) error estimates. It can improve the accuracy of the approximation, and the stress intensity factor can also be obtained from the extraction formula, cf. (3.6). (iv) The singular complement method [4] is based on an orthogonal decomposition of the solution \(u\) into a singular part and a regular part, by augmenting the FEM trial space with specially designed singular functions.

For 3D problems with edges, the singular functions belong to an infinite-dimensional space and their coefficients are functions defined along the edges [29]. Thus, their computation involves approximating functions defined along edges, and there are relatively few numerical methods, and numerical investigations are strikingly lacking. The methods developed for the 2D case do not extend directly to the 3D case. In fact, in several existing studies, numerical algorithms and error analysis have been provided, but the methods are nontrivial to implement [52, 57]. For example, the approach in [52] requires evaluating a few dozens of highly singular integrals at each step, which may lead to serious numerical issues.

### Our contributions

In this work, following the paradigm of building analytic knowledge of the problem into numerical schemes, we construct a novel numerical method using PINN to solve Poisson problems with geometric
singularities. The key analytic insight is that the solution \(u\) has a singular function representation as a linear combination of a singular function \(S\) and a regular part \(w\)[18, 19, 29, 38]: \(u=S+w\), with \(w\in H^{2}(\Omega)\). The singular function \(S\) is determined by the domain \(\Omega\), truncation functions and their coefficients. Using this fact, we develop, analyze and test SEPINN, and make the following contributions:

1. develop a novel class of SEPINNs for corner singularities, mixed boundary conditions and edge singularity, for the Poisson problem and modified Helmholtz problem.
2. provide error bounds for the SEPINN approximation.
3. present numerical experiments for several challenging scenarios, including 3D problems with edge singularities and the eigenvalue problem on an L-shaped domain, to illustrate the flexibility and accuracy of SEPINN. We also include a comparative study with existing approaches.

The rest of the paper is organized as follows. In Section 2 we recall preliminaries on DNNs and its use in PINNs. Then we develop the singularity enriched PINN in Section 3 for 2D case (corner singularity and mixed boundary conditions) and 3D case (edge singularity) separately, and also extend it to the Helmholtz case. In Section 4, we discuss the convergence analysis of SEPINN. In Section 5, we present extensive numerical experiments to illustrate the performance of SEPINN, including a comparative study with PINN and its variants, and give further discussions in Section 6.

## 2 Preliminaries

### Deep neural networks

We employ standard fully connected feedforward DNNs, i.e., functions \(f_{\theta}:\mathbb{R}^{d}\rightarrow\mathbb{R}\), with the DNN parameters \(\theta\in\mathbb{R}^{N_{\theta}}\) (\(N_{\theta}\) is the dimensionality of the DNN parameters). Given a sequence of integers \(\{n_{\ell}\}_{\ell=0}^{L}\), with \(n_{0}=d\) and \(n_{L}=1\), \(f_{\theta}(x)\) is defined recursively by

\[f^{(0)} =x,\] \[f^{(\ell)} =\rho(A_{\ell}f^{(\ell-1)}+b_{\ell}),\quad\ell=1,2,\cdots,L-1,\] \[f_{\theta}(x) :=f^{(L)}(x)=A_{L}f^{(L-1)}+b_{L},\]

where \(A_{\ell}\in\mathbb{R}^{n_{\ell}\times n_{\ell-1}}\) and \(b_{\ell}\in\mathbb{R}^{n_{\ell}},\,\ell=1,2,\cdots,L\), are the weight matrix and bias vector at the \(\ell\)th layer. The nonlinear activation function \(\rho:\mathbb{R}\rightarrow\mathbb{R}\) is applied componentwise to a vector. The integer \(L\) is called the depth, and \(W:=\max\{n_{\ell},\ell=0,1,\cdots,L\}\) the width of the DNN. The set of parameters \(\{A_{\ell},\,b_{\ell}\}_{\ell=1}^{L}\) of the DNN is trainable and stacked into a big vector \(\theta\). \(f^{(0)}\) is called the input layer, \(f^{(\ell)}\), \(\ell=1,2,\cdots,L-1\), are called the hidden layer and \(f_{\theta}(x)\) is the output layer.

There are many possible choices of \(\rho\). The most frequently used one in computer vision is rectified linear unit (ReLU), \(\rho(x)=\max(x,0)\). However, it is not smooth enough for PINN, since PINN requires thrice differentiability of \(\rho\): two spatial derivatives in the loss, and another derivative in the DNN parameters \(\theta\) (for the optimizer). For neural PDE solvers, hyperbolic tangent \(\rho(x)=\frac{\mathrm{e}^{x}-\mathrm{e}^{-x}}{\mathrm{e}^{x}+\mathrm{e}^{-x}}\) and logistic / sigmoid \(\rho(x)=\frac{1}{1+\mathrm{e}^{-x}}\) are often used [53, 17]. We employ the hyperbolic tangent. We denote the collection of DNN functions of depth \(L\), with \(N_{\theta}\) nonzero parameters, and each the parameter bounded by \(R\), with the activation function \(\rho\) by \(\mathcal{N}_{\rho}(L,N_{\theta},R)\), i.e., \(\mathcal{N}_{\rho}(L,N_{\theta},R)=\{w_{\theta}:w_{\theta}\) has a depth \(L\), \(|\theta|_{0}\leq N_{\theta},\,|\theta|_{\ell^{\infty}}\leq R\}\), where \(|\cdot|_{\ell^{0}}\) and \(|\cdot|_{\ell^{\infty}}\) denote the number of nonzero entries in and the maximum norm of a vector, respectively. We also use the notation \(\mathcal{A}\) to denote this collection of functions.

### Physics informed neural networks

Physics informed neural networks (PINNs) [53] represent one popular neural solver based on the principle of PDE residual minimization. For problem (1.1), the continuous loss \(\mathcal{L}(u)\) is given by

\[\mathcal{L}_{\mathbf{\sigma}}(u)=\|\Delta u+f\|_{L^{2}(\Omega)}^{2}+\sigma_{d}\| u\|_{L^{2}(\Gamma_{D})}^{2}+\sigma_{n}\|\partial_{n}u\|_{L^{2}(\Gamma_{N})}^{2},\]

where the tunable penalty weights \(\sigma_{d},\sigma_{n}>0\) are to approximately enforce the boundary conditions, and \(\mathbf{\sigma}=(\sigma_{d},\sigma_{n})\). We approximate the solution \(u\) by an element \(u_{\theta}\in\mathcal{A}\), and then discretize relevant integrals using quadrature, e.g., Monte Carlo method. Let \(U(D)\) be the uniform distribution over a set \(D\), and \(|D|\) denotes taking its Lebesgue measure. Then we can rewrite the loss \(\mathcal{L}(u_{\theta})\) as

\[\mathcal{L}_{\mathbf{\sigma}}(u_{\theta})=|\Omega|\mathbb{E}_{U(\Omega)}[(\Delta u _{\theta}(X)+f(X))^{2}]+\sigma_{d}|\Gamma_{D}|\mathbb{E}_{U(\Gamma_{D})}[(u_{ \theta}(Y))^{2}]+\sigma_{n}|\Gamma_{N}|\mathbb{E}_{U(\Gamma_{N})}[(\partial_{ n}u_{\theta}(Z))^{2}],\]
where \(\mathbb{E}_{\mathbf{\nu}}\) denotes taking expectation with respect to a distribution \(\nu\). Let the sampling points \(\{X_{i}\}_{i=1}^{N_{\pi}}\), \(\{Y_{j}\}_{j=1}^{N_{\pi}}\) and \(\{Z_{k}\}_{k=1}^{N_{\pi}}\) be identically and independently distributed (i.i.d.), uniformly on the domain \(\Omega\) and the boundaries \(\Gamma_{D}\) and \(\Gamma_{N}\), respectively, i.e., \(\{X_{i}\}_{i=1}^{N_{\pi}}\sim U(\Omega)\), \(\{Y_{j}\}_{j=1}^{N_{\pi}}\sim U(\Gamma_{D})\) and \(\{Z_{k}\}_{k=1}^{N_{\pi}}\sim U(\Gamma_{N})\). Then the empirical loss function \(\widehat{\mathcal{L}}_{\mathbf{\sigma}}(u_{\theta})\) is given by

\[\widehat{\mathcal{L}}_{\mathbf{\sigma}}(u_{\theta})=\frac{|\Omega|}{N_{r}}\sum_{ i=1}^{N_{\pi}}(\Delta u_{\theta}(X_{i})+f(X_{i}))^{2}+\frac{\sigma_{d}|\Gamma_{D}|}{N _{d}}\sum_{j=1}^{N_{d}}(u_{\theta}(Y_{j}))^{2}+\frac{\sigma_{n}|\Gamma_{N}|}{ N_{n}}\sum_{k=1}^{N_{n}}(\partial_{n}u_{\theta}(Z_{k}))^{2}.\]

Note that the resulting optimization problem \(\widehat{\mathcal{L}}_{\mathbf{\sigma}}(u_{\theta})\) over \(\mathcal{A}\) is well posed due to the box constraint on the DNN parameters \(\theta\), i.e., \(|\theta|_{\ell\sim\infty}\leq R\) for suitable \(R\), which induces a compact set in \(\mathbb{R}^{N_{\theta}}\). Meanwhile the empirical loss \(\widehat{\mathcal{L}}_{\mathbf{\sigma}}(u_{\theta})\) is continuous in \(\theta\), when \(\rho\) is smooth. In the abscence of the box constraint, the optimization problem might not have a finite minimizer.

The loss \(\widehat{\mathcal{L}}_{\mathbf{\sigma}}(u_{\theta})\) is minimized with respect to the DNN parameters \(\theta\). This is often achieved by gradient type algorithms, e.g., Adam [37] or limited memory BFGS [11], which returns an approximate minimizer \(\theta^{*}\). The DNN approximation to the PDE solution \(u\) is given by \(u_{\theta^{*}}\). Note that the major computational effort, e.g., gradient of \(\widehat{\mathcal{L}}_{\mathbf{\sigma}}(u_{\theta})\) with respect to the DNN parameters \(\theta\) and the DNN \(u_{\theta}(x)\) with respect to the input \(x\) can both be computed efficiently via automatic differentiation [7], which is available in many software platforms, e.g., PyTorch or Tensorflow. Thus, the method is very flexible and easy to implement, and applicable to a wide range of direct and inverse problems for PDEs [36].

The population loss \(\mathcal{L}_{\mathbf{\sigma}}(u_{\theta})\) and empirical loss \(\widehat{\mathcal{L}}_{\mathbf{\sigma}}(u_{\theta})\) have different minimizers, due to the presence of quadrature errors. The analysis of these errors is known as generalization error analysis in statistical learning theory [1]. The theoretical analysis of PINNs has been investigated in several works under different settings [55, 33, 49, 46]. These mathematical theories require that the solutions to the problems be smooth, e.g., \(C^{2}(\overline{\Omega})\) in order to achieve consistency [55] and even stronger regularity for convergence rates [33]. Such conditions unfortunately cannot be met for problem (1.1), due to the inherently limited solution regularity. Thus, it is not _a priori_ clear that one can successfully apply PINNs to problem (1.1). This is also confirmed by the numerical experiments in Section 5. In this work, we shall develop an effective strategy to overcome the challenge.

## 3 Singularity enriched PINN

Now we develop a class of singularity enriched PINN (SEPINN) for solving Poisson problems with geometric singularities, including 2D problems with mixed boundary conditions or on polygonal domains, and 3D problems with edge singularity. The key is the singular function representation. We discuss the 2D case in Section 3.1, and the more involved 3D case in Section 3.2. Finally we extend the approach to the modified Helmholtz equation in Section 3.3.

### Two-dimensional problem

#### 3.1.1 Singular function representation

First we develop a singular function representation in the 2D case. We determine the analytic structure of the solution \(u\) of problem (1.1) by means of the Fourier method. Consider a vertex \(\mathbf{v}_{j}\) of the polygonal domain \(\Omega\) with an interior angle \(\omega_{j}\). We denote by \((r_{j},\theta_{j})\) the local polar coordinate of the vertex \(\mathbf{v}_{j}\) so that the interior angle \(\omega_{j}\) is spanned counterclockwise by two rays \(\theta_{j}=0\) and \(\theta_{j}=\omega_{j}\). Then consider the local behavior of \(u\) near the vertex \(\mathbf{v}_{j}\), or in the sector \(G_{j}=\{(r_{j},\theta_{j}):0<r_{j}<R_{j},0<\theta_{j}<\omega_{j}\}.\) Let the edge overlapping with \(\theta_{j}=0\) be \(\Gamma_{j_{1}}\) and with \(\theta_{j}=\omega_{j}\) be \(\Gamma_{j_{2}}\). We employ the system of orthogonal and complete set of basis functions in \(L^{2}(0,\omega_{j})\) in Table 3.1, according to the boundary conditions (b.c.) on the edges \(\Gamma_{j_{1}}\) and \(\Gamma_{j_{2}}\)[51]. We denote the representation of \(u\) in the local polar coordinate by \(\widetilde{u}(r_{j},\theta_{j})\), i.e., \(\tilde{u}(r_{j},\theta_{j})=u(r_{j}\cos\theta_{j},r_{j}\sin\theta_{j})\).

To study the behavior of \(u\) in the sector \(G_{j}\), we assume \(\widetilde{u}(R_{j},\theta_{j})=0\) and \(|\widetilde{u}(0,\theta_{j})|<\infty\). Since \(u\) and \(f\) belong to \(L^{2}(G_{j})\), they can be represented in Fourier series with respect to \(\{\phi_{j,k}\}_{k=1}^{\infty}\):

\[u(x,y)=\widetilde{u}(r_{j},\theta_{j})=\sum_{k=1}^{\infty}u_{k}(r_ {j})\phi_{j,k}(\theta_{j}),\quad\text{in }G_{j}, \tag{3.1}\] \[f(x,y)=\widetilde{f}(r_{j},\theta_{j})=\sum_{k=1}^{\infty}f_{k}( r_{j})\phi_{j,k}(\theta_{j}),\quad\text{in }G_{j}, \tag{3.2}\]
with the Fourier coefficients \(u_{k}(r_{j})\) and \(f_{k}(r_{j})\) given respectively by \(u_{k}(r_{j})=\frac{2}{\omega_{j}}\int_{0}^{\omega_{j}}\widetilde{u}(r_{j},\theta_ {j})\phi_{j,k}(\theta_{j})\mathrm{d}\theta_{j}\) and \(f_{k}(r_{j})=\frac{2}{\omega_{j}}\int_{0}^{\omega_{j}}\widetilde{f}(r_{j}, \theta_{j})\phi_{j,k}(\theta_{j})\mathrm{d}\theta_{j}\). Substituting (3.1) and (3.2) into (1.1) gives the following two-point boundary value problem:

\[\begin{cases}\quad-u_{j,k}^{\prime\prime}-r^{-1}u_{j,k}^{\prime}+\lambda_{j,k} ^{2}r^{-2}u_{j,k}=f_{k},\quad 0<r<R_{j},\\ |u_{j,k}(r)|_{r=0}<\infty,\quad u_{j,k}(r)|_{r=R_{j}}=0,\end{cases}\]

Solving the ODE yields [51, (2.10)]

\[u_{k}(r_{j})= r_{j}^{\lambda_{j,k}}\left(\frac{1}{2\lambda_{j,k}}\int_{r_{j}}^ {R_{j}}f_{k}(\tau)\tau^{1-\lambda_{j,k}}\mathrm{d}\tau-\frac{1}{2\lambda_{j,k} R_{j}^{2\lambda_{j,k}}}\int_{0}^{R_{j}}f_{k}(\tau)\tau^{1+\lambda_{j,k}} \mathrm{d}\tau\right)+\frac{r_{j}^{-\lambda_{j,k}}}{2\lambda_{j,k}}\int_{0}^{ r_{j}}f_{k}(\tau)\tau^{1+\lambda_{j,k}}\mathrm{d}\tau.\]

Since the factor in front of \(r_{j}^{\lambda_{j,k}}\) is generally nonzero, and \(r_{j}^{\lambda_{j,k}}\notin H^{2}(G)\) for \(\lambda_{j,k}<1\), there exist singular terms in the representation (3.1). This shows the limited regularity of the solution \(u\), which is the culprit of low efficiency of standard numerical schemes.

Next we list all singularity points and the corresponding singular functions. Let \(\mathbf{v}_{j},j=1,2,\cdots,M\), be the vertices of \(\Omega\) whose interior angles \(\omega_{j},j=1,2,\cdots,M\), satisfy

\[\begin{cases}\quad\pi<\omega_{j}<2\pi,&\text{b.c. doesn't change its type},\\ \pi/2<\omega_{j}<2\pi,&\text{b.c. changes its type}.\end{cases} \tag{3.3}\]

Table 3.2 gives the index set \(\mathbb{I}_{j}\) and the associated singularity functions [14, p. 2637]. Note that \(s_{j,1}\in H^{1+\frac{\pi}{2\omega}-\epsilon}(\Omega)\), \(s_{j,\frac{1}{2}}\in H^{1+\frac{\pi}{2\omega}-\epsilon}(\Omega)\), and \(s_{j,\frac{3}{2}}\in H^{1+\frac{\pi}{2\omega}-\epsilon}(\Omega)\) for small \(\epsilon>0\). Upon letting

\[\omega^{*}=\max_{1\leq j\leq M}\hat{\omega_{j}},\quad\text{with }\hat{\omega_{j}}= \begin{cases}\omega_{j},&\text{b.c. does't change its type at }\mathbf{v}_{j},\\ 2\omega_{i},&\text{b.c. changes its type at }\mathbf{v}_{j},\end{cases}\]

the solution \(u\) belongs to \(H^{1+\frac{\pi}{2\omega}-\epsilon}(\Omega)\), just falling short of \(H^{1+\frac{\pi}{2\omega}}(\Omega)\). If \(\omega^{*}>\pi\), \(u\) fails to belong to \(H^{2}(\Omega)\). Hence, it is imperative to develop techniques to resolve such singularities.

We employ a smooth cut-off function. For \(\rho\in(0,2]\), we define a cut-off function \(\eta_{\rho}\) by

\[\eta_{\rho}(r_{j})=\begin{cases}\quad 1,&0<r_{j}<\frac{\rho R}{2},\\ \frac{15}{16}\left(\frac{8}{15}-\left(\frac{4r_{j}}{\rho R}-3\right)+\frac{2} {3}\left(\frac{4r_{j}}{\rho R}-3\right)^{3}-\frac{1}{5}\left(\frac{4r_{j}}{ \rho R}-3\right)^{5}\right),&\frac{\rho R}{2}\leq r_{j}<\rho R,\\ 0,&r_{j}\geq\rho R,\end{cases} \tag{3.4}\]

where \(R\in\mathbb{R}_{+}\) is a fixed number so that \(\eta_{\rho}\) vanishes identically on \(\partial\Omega\). In practice, we take \(R\) to be small enough so that when \(i\neq j\), the support of \(\eta_{\rho}(r_{i})\) does no intersect with that of \(\eta_{\rho}(r_{j})\). By construction, we have \(\eta_{\rho}\in C^{2}([0,\infty))\). Then the solution \(u\) of problem (1.1) has the following singular function representation [5, 18, 43]

\[u=w+\sum_{j=1}^{M}\sum_{i\in\mathbb{I}_{j}}\gamma_{j,i}\eta_{\rho_{j}}(r_{j})s_ {j,i}(r_{j},\theta_{j}),\quad\text{with }w\in H^{2}(\Omega)\cap H^{1}_{0}(\Omega), \tag{3.5}\]

where the scalars \(\gamma_{j,i}\in\mathbb{R}\) are known as stress intensity factors and given by the following extraction formulas [29, Lemma 8.4.3.1]:

\[\gamma_{j,i}=\frac{1}{i\pi}\left(\int_{\Omega}f\eta_{\rho_{j}}s_{j,-i}\mathrm{ d}x+\int_{\Omega}u\Delta(\eta_{\rho_{j}}s_{j,-i})\mathrm{d}x\right), \tag{3.6}\]
where \(s_{j,-i}\) denotes dual singular functions of \(s_{j,i}\). Specifically, if \(s_{j,i}(r_{j},\theta_{j})=r_{j}^{\frac{i\pi}{\nu_{j}}}\sin\frac{i\pi\theta_{j}}{ \omega_{j}}\), the dual function \(s_{j,-i}\) is given by \(s_{j,-i}(r_{j},\theta_{j})=r_{j}^{\frac{i\pi}{\nu_{j}}}\sin\frac{i\pi\theta_{j}} {\omega_{j}}\)[14, p. 2639]. Moreover, the following regularity estimate on the regular part \(w\) holds

\[\|w\|_{H^{2}(\Omega)}+\sum_{j=1}^{M}\sum_{i\in\mathbb{I}_{j}}|\gamma_{j,i}|\leq c \|f\|_{L^{2}(\Omega)}. \tag{3.7}\]

#### 3.1.2 Singularity enriched physics-informed neural network

Now we propose singularity enriched PINN (SEPINN) for problem (1.1), inspired by the representation (3.5). We discuss only the case with one singular function in \(\Omega\) satisfying the condition (3.3). The case of multiple singularities can be handled similarly. Upon letting \(S=\gamma\eta_{\rho}s\), the regular part \(w\) satisfies

\[\left\{\begin{aligned} -\Delta w&=f+\gamma\Delta(\eta_{ \rho}s),&\text{in }\Omega,\\ w&=0,&\text{on }\Gamma_{D},\\ \partial_{n}w&=0,&\text{on }\Gamma_{N}, \end{aligned}\right. \tag{3.8}\]

where the parameter \(\gamma\) is unknown. Since \(w\in H^{2}(\Omega)\), it can be well approximated using PINN. The parameter \(\gamma\) can be either learned together with \(w\) or extracted from \(w\) via (3.6).

Based on the principle of PDE residual minimization, the solution \(w^{*}\) of problem (3.8) and the exact parameter \(\gamma^{*}\) in the splitting (3.5) is a global minimizer of the following loss

\[\mathcal{L}_{\mathbf{\sigma}}(w;\gamma)=\|\Delta w+f+\gamma\Delta(\eta_{\rho}s) \|_{L^{2}(\Omega)}^{2}+\sigma_{d}\|w\|_{L^{2}(\Gamma_{D})}^{2}+\sigma_{n}\left\| \partial_{n}w\right\|_{L^{2}(\Gamma_{N})}^{2}, \tag{3.9}\]

where the penalty weights \(\mathbf{\sigma}=(\sigma_{d},\sigma_{n})\in\mathbb{R}_{+}^{2}\) are tunable. Following the PINN paradigm in Section 2.2, we employ a DNN \(w_{\theta}\in\mathcal{A}\) to approximate \(w^{*}\in H^{2}(\Omega)\), and treat the parameter \(\gamma\) as a trainable parameter and learn it along with the DNN parameters \(\theta\). This leads to an empirical los

\[\widehat{\mathcal{L}}_{\mathbf{\sigma}}(w_{\theta};\gamma)= \frac{|\Omega|}{N_{r}}\sum_{i=1}^{N_{r}}(\Delta w_{\theta}(X_{i})+f(X _{i})+\gamma\Delta(\eta_{\rho}s)(X_{i}))^{2}+\sigma_{d}\frac{|\Gamma_{D}|}{N_ {d}}\sum_{j=1}^{N_{d}}w_{\theta}^{2}(Y_{j})+\sigma_{n}\frac{|\Gamma_{N}|}{N_{ n}}\sum_{k=1}^{N_{n}}\left(\partial_{n}w_{\theta}(Z_{k})\right)^{2}, \tag{3.10}\]

with i.i.d. sampling points \(\{X_{i}\}_{i=1}^{N_{r}}\sim U(\Omega)\), \(\{Y_{j}\}_{j=1}^{N_{d}}\sim U(\Gamma_{D})\) and \(\{Z_{k}\}_{k=1}^{N_{n}}\sim U(\Gamma_{N})\). Let \((\widehat{\theta}^{*},\widehat{\gamma}^{*})\) be a minimizer of the empirical loss \(\widehat{L}(w_{\theta};\gamma)\). Then \(w_{\widehat{\theta}^{*}}\in\mathcal{A}\) is the DNN approximation of the regular part \(w^{*}\), and the approximation \(\hat{u}\) to \(u\) is given by \(\hat{u}=w_{\theta^{*}}+\hat{\gamma}^{*}\eta_{\rho}s\).

\begin{table}
\begin{tabular}{c|c|c} \hline \hline \(\Gamma_{j2}\backslash\Gamma_{j1}\) & Dirichlet & Neumann \\ \hline Dirichlet & \(s_{j,1}(r_{j},\theta_{j})=r_{j}^{\frac{\pi}{\nu_{j}}}\sin\frac{\pi\theta_{j}}{ \omega_{j}}\), & \(s_{j,\frac{1}{2}}(r_{j},\theta_{j})=r_{j}^{\frac{\pi}{\nu_{j}}}\cos\frac{\pi \theta_{j}}{2\omega_{j}}\), if \(\frac{\pi}{2}<\omega_{j}\leq\frac{3\pi}{2}\), \\  & \(\mathbb{I}_{j}=\left\{\frac{1}{2}\right\}\). & \(\mathbb{I}_{j}=\left\{\frac{1}{2}\right\}\). \\ \cline{2-3}  & \(\mathbb{I}_{j}=\left\{1\right\}\). & \(s_{j,\frac{1}{2}}(r_{j},\theta_{j})=r_{j}^{\frac{\pi}{\nu_{j}}}\cos\frac{\pi \theta_{j}}{2\omega_{j}}\) and \\  & \(s_{j,\frac{3}{2}}(r_{j},\theta_{j})=r_{j}^{\frac{3\pi}{2\omega_{j}}}\cos\frac{3 \pi\theta_{j}}{2\omega_{j}}\), if \(\frac{3\pi}{2}<\omega_{j}\leq 2\pi\), \\  & \(\mathbb{I}_{j}=\left\{\frac{1}{2},\frac{3}{2}\right\}\). & \\ \hline Neumann & \(s_{j,\frac{1}{2}}(r_{j},\theta_{j})=r_{j}^{\frac{\pi}{2\omega_{j}}}\sin\frac{ \pi\theta_{j}}{2\omega_{j}}\), if \(\frac{\pi}{2}<\omega_{j}\leq\frac{3\pi}{2}\), \\  & \(\mathbb{I}_{j}=\left\{\frac{1}{2}\right\}\). & \(s_{j,1}(r_{j},\theta_{j})=r_{j}^{\frac{\pi}{2\omega_{j}}}\cos\frac{\pi\theta_{ j}}{\omega_{j}}\), \\  & \(s_{j,\frac{1}{2}}(r_{j},\theta_{j})=r_{j}^{\frac{3\pi}{2\omega_{j}}}\sin\frac{3 \pi\theta_{j}}{2\omega_{j}}\) and \\  & \(s_{j,\frac{3}{2}}(r_{j},\theta_{j})=r_{j}^{\frac{3\pi}{2\omega_{j}}}\sin\frac{3 \pi\theta_{j}}{2\omega_{j}}\), if \(\frac{3\pi}{2}<\omega_{j}\leq 2\pi\), \\  & \(\mathbb{I}_{j}=\left\{\frac{1}{2},\frac{3}{2}\right\}\). & \\ \hline \hline \end{tabular}
\end{table}
Table 3.2: Singularity functions \(s_{j,i}\), depending on the boundary condition. The tuple \((r_{j},\theta_{j})\) refers to the local polar coordinate of the vertex \(\mathbf{v}_{j}\), and \(\mathbb{I}_{j}\) is an index set for leading singularities associated with \(\mathbf{v}_{j}\).
Now we discuss the training of the loss \(\widehat{\mathcal{L}}_{\mathbf{\sigma}}(w_{\theta};\gamma)\). One can minimize \(\widehat{\mathcal{L}}_{\mathbf{\sigma}}(w_{\theta};\gamma)\) directly with respect to \(\theta\) and \(\gamma\), which works reasonably. However, the DNN approximation \(\widehat{u}\) tends to have larger errors on the boundary \(\partial\Omega\) than in the domain \(\Omega\), but the estimated \(\widehat{\gamma}^{*}\) is often accurate. Thus we adopt a two-stage training procedure.

1. At Stage 1, minimize the loss \(\widehat{\mathcal{L}}_{\mathbf{\sigma}}(w_{\theta};\gamma)\) for a fixed \(\mathbf{\sigma}\), and obtain the minimizer \((\widehat{\theta}^{*},\widehat{\gamma}^{*})\).
2. At Stage 2, fix \(\gamma\) in \(\widehat{\mathcal{L}}_{\mathbf{\sigma}}(w_{\theta};\gamma)\) at \(\widehat{\gamma}^{*}\), and learn \(\theta\) via a path-following strategy [45, 31].

Now we describe the path-following strategy on updating \(\mathbf{\sigma}\). We start with small values \(\mathbf{\sigma}^{(1)}=\mathbf{\sigma}\). After each loop (i.e., finding one minimizer \(\widehat{\theta}^{*}_{k}\)), we update \(\mathbf{\sigma}\) geometrically: \(\mathbf{\sigma}^{(k+1)}=q\mathbf{\sigma}^{(k)}\), with \(q>1\). By updating \(\mathbf{\sigma}^{(k)}\), the minimizer \(\widehat{\theta}^{*}_{k}\) of the loss \(\widehat{\mathcal{L}}_{\mathbf{\sigma}^{(k)}}(w_{\theta};\widehat{\gamma}^{*})\) also approaches that of problem (3.9), and the path-following strategy enforces the boundary conditions progressively, which is beneficial to obtain good approximations, since when \(\mathbf{\sigma}\) are large, the optimization problem is known to be stiff (and thus numerically challenging). Note that the minimizer \(\widehat{\theta}^{*}_{k+1}\) of the \(\mathbf{\sigma}^{(k+1)}\)-problem (i.e., minimizing \(\widehat{\mathcal{L}}_{\mathbf{\sigma}^{(k+1)}}(w_{\theta};\widehat{\gamma}^{*})\)) can be initialized to \(\widehat{\theta}^{*}_{k}\) of the \(\mathbf{\sigma}^{(k)}\)-problem to warm start the training process. Hence, for each fixed \(\mathbf{\sigma}^{(k)}\) (except \(\mathbf{\sigma}^{(1)}\)), the initial parameter configuration is close to the optimal one, and the training only requires few iterations to reach convergence. The overall procedure is shown in Algorithm 1 for 2D problems with corner singularities and / or mixed boundary conditions.

```
1:Set \(\mathbf{\sigma}^{(1)}\), and obtain the minimizer \((\widehat{\theta}^{*},\widehat{\gamma}^{*})\) of the loss \(\widehat{\mathcal{L}}_{\mathbf{\sigma}^{(1)}}(w_{\theta};\gamma)\).
2:Set \(k=1\), \(\widehat{\theta}^{*}_{0}=\widehat{\theta}^{*}\), and increasing factor \(q>1\).
3:while Stopping condition not met do
4: Find a minimizer \(\widehat{\theta}^{*}_{k}\) of the loss \(\widehat{\mathcal{L}}_{\mathbf{\sigma}^{(k)}}(w_{\theta};\widehat{\gamma}^{*})\) (initialized to \(\widehat{\theta}^{*}_{k-1}\)).
5: Update \(\mathbf{\sigma}\) by \(\mathbf{\sigma}^{(k+1)}=q\mathbf{\sigma}^{(k)}\), and \(k\gets k+1\).
6:Output the SEPINN approximation \(\widehat{u}=w_{\widehat{\theta}^{*}_{k-1}}+\widehat{\gamma}^{*}\eta s\).
```

**Algorithm 1** SEPINN for 2D problems

### Three-dimensional problem

Now we develop SEPINN for the 3D Poisson problem with edge singularities.

#### 3.2.1 Singular function representation

There are several different types of geometric singularities in the 3D case, and each case has to be dealt with separately [19]. We only study edge singularities, which cause strong solution singularities. Indeed, the \(H^{2}(\Omega)\)-regularity of the solution \(u\) of problem (1.1) is not affected by the presence of conic points [38, 29, 40]. Now we state the precise setting. Let \(\Omega_{0}\subset\mathbb{R}^{2}\) be a polygonal domain as in Section 3.1, and \(\Omega=\Omega_{0}\times(0,l)\). Like before, let \(\mathbf{v}_{j},j=1,2,\cdots,M\), be the vertices of \(\Omega_{0}\) whose interior angles \(\omega_{j},j=1,2,\cdots,M\), satisfy (3.3), and let \(\Gamma_{z_{1}}=\Omega_{0}\times\{0\}\) and \(\Gamma_{z_{2}}=\Omega_{0}\times\{l\}\). We employ a complete orthogonal system \(\{Z_{j,n}\}_{n=0}^{\infty}\) of \(L^{2}(0,l)\), given in Table 3.3.

The functions \(u\in L^{2}(\Omega)\) and \(f\in L^{2}(\Omega)\) from problem (1.1) can then be represented by the following
convergent Fourier series in the 3D wedge \(G_{j}=\mathbb{K}_{j}\times(0,l)\) (suppressing the subscript \(j\)):

\[u(x,y,z) =\frac{1}{2}u_{0}(x,y)Z_{0}(z)+\sum_{n=1}^{\infty}u_{n}(x,y)Z_{n}(z), \tag{3.11}\] \[f(x,y,z) =\frac{1}{2}f_{0}(x,y)Z_{0}(z)+\sum_{n=1}^{\infty}f_{n}(x,y)Z_{n}( z), \tag{3.12}\]

where the Fourier coefficients \(\{u_{n}\}_{n\in\mathbb{N}}\) and \(\{f_{n}\}_{n\in\mathbb{N}}\) are defined on the 2D domain \(\Omega_{0}\) by

\[u_{n}(x,y)=\frac{2}{l}\int_{0}^{l}u(x,y,z)Z_{n}(z)\mathrm{d}z\quad\text{and} \quad f_{n}(x,y)=\frac{2}{l}\int_{0}^{l}f(x,y,z)Z_{n}(z)\mathrm{d}z.\]

By substituting (3.11) and (3.12) into (1.1), we get countably many 2D elliptic problems:

\[\begin{cases}-\Delta u_{n}+\xi_{j,n}^{2}u_{n}=f_{n},&\text{in }\Omega_{0},\\ \phantom{-\Delta u_{n}+\xi_{j,n}^{2}u_{n}=}{}0,&\text{on }\Gamma_{D},\\ \phantom{-\Delta u_{n}+\xi_{j,n}^{2}u_{n}=}{}0,&\text{on }\Gamma_{N}.\end{cases} \tag{3.13}\]

Then problem (1.1) can be analyzed via the 2D problems. Below we describe the edge behavior of the weak solution \(u\in H^{1}(\Omega)\)[52, Theorem 2.1]. The next theorem gives a crucial decomposition of \(u\in H^{1+}\overrightarrow{\omega}^{*}\leftarrow\) for every \(\epsilon>0\). The functions \(\Psi_{j,i}\) are the so-called edge flux intensity functions.

**Theorem 3.1**.: _For any fixed \(f\in L^{2}(\Omega)\), let \(u\in H^{1}(\Omega)\) be the unique weak solution to problem (1.1). Then there exist unique functions \(\Psi_{j,i}\in H^{1-\lambda_{j,i}}(0,l)\) of the variable \(z\) such that \(u\) can be split into a sum of a regular part \(w\in H^{2}(\Omega)\) and a singular part \(S\) with the following properties:_

\[u =w+S,\quad S=\sum_{j=1}^{M}\sum_{i\in\mathbb{I}_{j}}S_{j,i}(x_{j}, y_{j},z), \tag{3.14a}\] \[S_{j,i}(x_{j},y_{j},z) =(T_{j}(r_{j},z)*\Psi_{j,i}(z))\eta_{\rho_{j}}(r_{j})s_{j,i}(r_{j},\theta_{j}), \tag{3.14b}\]

_where the functions \(T_{j}\) are fixed Poisson's kernels, and the symbol \(*\) denotes the convolution in \(z\in(0,l)\). Moreover, there exists a constant \(C>0\) independent of \(f\in L^{2}(\Omega)\) such that \(\|w\|_{H^{2}(\Omega)}\leq C\|f\|_{L^{2}(\Omega)}\)._

Now we assume that near each edge \(\mathbf{v}_{j}\times(0,l)\), the domain \(\Omega\) coincides with a 3D wedge \(G_{j}\) defined by \(G_{j}=\{(x_{j},y_{j},z)\in\Omega:0<r_{j}<R_{j},0<\theta_{j}<\omega_{j},0<z<l\}\), where \((r_{j},\theta_{j})\) are local polar coordinates linked with the local Cartesian coordinates \((x_{j},y_{j})\). The explicit form of the function \(T_{j}(r_{j},z)*\Psi_{j,i}(z)\) and the formula for the coefficients \(\gamma_{j,i,n}\) are given below [52, Theorem 2.2][51, pp. 179-182].

**Theorem 3.2**.: _The coefficients \(\Phi_{j,i}(x_{j},y_{j},z)=T_{j}(r_{j},z)*\Psi_{j,i}(z)\) of the singularities in (3.14b) can be represented by Fourier series in \(z\) and with respect to the orthogonal system \(\{Z_{n}(z)\}_{n=0}^{\infty}\):_

\[T_{j}(r_{j},z) =\frac{1}{2}Z_{0}(z)+\sum_{n=1}^{\infty}\mathrm{e}^{-\xi_{j,n}r_{ j}}Z_{n}(z),\] \[\Psi_{j,i}(z) =\frac{1}{2}\gamma_{j,i,0}Z_{0}(z)+\sum_{n=1}^{\infty}\gamma_{j, i,n}Z_{n}(z),\] \[\Phi_{j,i}(x_{j},y_{j},z) =\frac{1}{2}\gamma_{j,i,0}Z_{0}(z)+\sum_{n=1}^{\infty}\gamma_{j, i,n}\mathrm{e}^{-\xi_{j,n}r_{j}}Z_{n}(z), \tag{3.15}\]

_where the coefficients \(\gamma_{j,i,n}\) are given explicitly by_

\[\gamma_{j,i,n} =\frac{2}{l\omega_{j}\lambda_{j}}\int_{G_{j}}f_{j}^{*}\mathrm{e}^ {\xi_{j,n}r_{j}}s_{j,-k}(r_{j},\theta_{j})Z_{n}(z)\mathrm{d}x\mathrm{d}y\mathrm{ d}z, \tag{3.16}\] \[f_{j}^{*} =f\eta_{\rho_{j}}-u\left(\frac{\partial^{2}\eta_{\rho_{j}}}{ \partial r_{j}^{2}}+\left(2\xi_{j,n}+\frac{1}{r_{j}}\right)\frac{\partial\eta _{\rho_{j}}}{\partial r_{j}}+\left(2\xi_{j,n}^{2}+\frac{\xi_{j,n}}{r_{j}} \right)\eta_{\rho_{j}}\right)-2\frac{\partial u}{\partial r_{j}}\left(\frac{ \partial\eta_{\rho_{j}}}{\partial r_{j}}+\xi_{j,n}\eta_{\rho_{j}}\right).\]

_Moreover there exists a constant \(C>0\) independent of \(f\) such that_

\[|\gamma_{j,i,0}|^{2}+\sum_{n=1}^{\infty}\xi_{j,n}^{2(1-\lambda_{j,i})}|\gamma_{j, i,n}|^{2}\leq C\|f\|_{L^{2}(G_{j})}.\]
Note that the formula (3.16) is a singular integral and numerically inconvenient to evaluate. We discuss the case of only one edge with one singular function below. Then the solution \(u\) can be split into

\[u=w+\Phi\eta s. \tag{3.17}\]

Since the functions \(\Phi\eta_{\rho}s\) and \(\partial_{n}(\Phi\eta_{\rho}s)\) vanish on \(\Gamma_{D}\) and \(\Gamma_{N}\), respectively, \(w\) solves

\[\begin{cases}-\Delta w=f+\Delta(\Phi\eta_{\rho}s),&\text{in }\Omega,\\ w=0,&\text{on }\Gamma_{D},\\ \partial_{n}w=0,&\text{on }\Gamma_{N}.\end{cases} \tag{3.18}\]

This forms the basis of SEPINN for 3D problems with edge singularities. Below we describe two strategies for constructing SEPINN approximations, i.e., SEPINN-C based on a cutoff approximation of the infinite series and SEPINN-N based on multiple DNN approximations.

#### 3.2.2 Sepinn - Cutoff approximation

The expansion (3.15) of the singular function \(\Phi\) in (3.18) involves infinitely many unknown scalar coefficients \(\{\gamma_{j,i,n}\}_{=0}^{\infty}\). In practice, it is infeasible to learn all of them. However, since the series is convergent, we may truncate it to a finite number of terms: the edge flux intensity function \(\Phi(x_{j},y_{j},z)\) from (3.14b) is approximated by the truncation

\[\Phi_{j,i}^{N}(x_{j},y_{j},z)=\frac{1}{2}\gamma_{j,i,0}+\sum_{n=1}^{N}\gamma_ {j,i,n}\mathrm{e}^{-\xi_{j,n}r_{j}}Z_{n}(z),\quad\text{with }N\in\mathbb{N}.\]

The approximate singular function \(s_{j,i}^{N}\) is given by

\[S_{j,i}^{N}(x_{j},y_{j},z)=\Phi_{j,i}^{N}(x_{j},y_{j},z)\eta_{\rho_{j}}(r_{j}) s_{j,i}(r_{j},\theta_{j}). \tag{3.19}\]

In view of the splitting (3.14a) and truncation (3.19), the approximate regular part \(w\) satisfies

\[\begin{cases}-\Delta w=f+\gamma_{0}\Delta(\eta_{\rho}s)+\sum_{n=1}^{N}\gamma_ {n}\Delta(\mathrm{e}^{-\xi_{n}r}Z_{n}(z)\eta_{\rho}s),&\text{in }\Omega,\\ \quad\quad w=0,&\text{on }\Gamma_{D},\\ \partial_{n}w=0,&\text{on }\Gamma_{N},\end{cases} \tag{3.20}\]

with \(\gamma_{i},i=0,1,\cdots,N\), being \(N+1\) unknown parameters. Let \(w^{N}\) be the solution of (3.20) and \(w^{*}\) the solution of (3.18). The next result shows that when \(N\) is large enough, the error \(w^{N}-w^{*}\) can be made small, and so is the error between \(u^{N}\) and \(u^{*}\), which underpins the truncation method.

**Theorem 3.3**.: _Let \(w^{N}\) be the solution of (3.20) and \(S^{N}\) be the truncated singular function. Let \(u^{N}=w^{N}+S^{N}\) and let \(u^{*}\) be the solution of (1.1). Then there holds_

\[\|u^{N}-u^{*}\|_{H^{1}(G)}\leq CN^{-1}\|f\|_{L^{2}(G)}.\]

Proof.: By [52, Lemma 3.3], there exists a constant \(C>0\) independent of \(f\) such that

\[\|S^{N}-S^{*}\|_{H^{1}(G)}\leq CN^{-1}\|f\|_{L^{2}(G)}. \tag{3.21}\]

Since \(w^{N}\) solves (3.20) and \(w^{*}\) solves (3.18), \(w^{*}-w^{N}\) satisfies problem (1.1) with the source \(f^{N}=\sum_{n=N+1}^{\infty}\gamma_{n}\Delta(\mathrm{e}^{-\xi_{n}r}Z_{n}(z) \eta_{\rho}s)\). By elliptic regularity theory, we have

\[\|w^{*}-w^{N}\|_{H^{1}(\Omega)}\leq C\|f^{N}\|_{H^{-1}_{D}(\Omega)}\leq C\| \Delta(S^{N}-S^{*})\|_{H^{-1}_{\Gamma_{D}}(\Omega)},\]

where the notation \(H^{-1}_{\Gamma_{D}}(\Omega)\) denotes the dual space of \(H^{1}_{\Gamma_{D}}(\Omega)=\{v\in H^{1}(\Omega):v=0\text{ on }\Gamma_{D}\}\). Upon integration by parts, we get

\[\|\Delta(S^{N}-S^{*})\|_{H^{-1}_{\Gamma_{D}}(\Omega)}=\sup_{v\in H ^{1}_{\Gamma_{D}}(\Omega),\|v\|_{H^{1}(\Omega)}\leq 1}\left|\int_{\Omega} \Delta(S^{N}-S^{*})v\mathrm{d}\mathbf{x}\right|\] \[=\sup_{v\in H^{1}_{\Gamma_{D}}(\Omega),\|v\|_{H^{1}(\Omega)}\leq 1} \left|\int_{\Omega}\nabla(S^{N}-S^{*})\cdot\nabla v\mathrm{d}x\right|\leq\|S^{ N}-S^{*}\|_{H^{1}(\Omega)}.\]

Combining this estimate with (3.21) yields the desired assertion.
Following the PINN paradigm in Section 2.2, we approximate \(w^{N}\) by an element \(w_{\theta}\in\mathcal{A}\). Like in the 2D case, we view the parameters \(\boldsymbol{\gamma}_{N}:=(\gamma_{i},\dots,\gamma_{N})\) as trainable parameters and learn them along with DNN parameters \(\theta\). The population loss \(\mathcal{L}_{\boldsymbol{\sigma}}(w_{\theta};\boldsymbol{\gamma}_{N})\) is given by

\[\mathcal{L}_{\boldsymbol{\sigma}}(w_{\theta};\boldsymbol{\gamma}_{N})=\left\| \Delta w_{\theta}+f+\gamma_{0}\Delta(\eta_{\rho}s)+\sum_{n=1}^{N}\gamma_{n} \Delta(\mathrm{e}^{-\xi_{n}r}Z_{n}(z)\eta_{\rho}s)\right\|_{L^{2}(\Omega)}^{2}+ \sigma_{d}\|w_{\theta}\|_{L^{2}(\Gamma_{D})}^{2}+\sigma_{n}\left\|\partial_{n} w_{\theta}\right\|_{L^{2}(\Gamma_{N})}^{2}.\]

In practice, we employ the following empirical loss

\[\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}(w_{\theta};\boldsymbol {\gamma}_{N})= \frac{|\Omega|}{N_{r}}\sum_{i=1}^{N_{r}}\left(\Delta w_{\theta}(X _{i})+f(X_{i})+\gamma_{0}\Delta(\eta_{\rho}s)(X_{i})+\sum_{n=1}^{N}\gamma_{n} \Delta(\mathrm{e}^{-\xi_{n}r}Z_{n}(z)\eta_{\rho}s)(X_{i})\right)^{2}\] \[+\sigma_{d}\frac{|\Gamma_{D}|}{N_{d}}\sum_{j=1}^{N_{d}}w_{\theta} ^{2}(Y_{j})+\sigma_{n}\frac{|\Gamma_{N}|}{N_{n}}\sum_{k=1}^{N_{n}}\left( \partial_{n}w_{\theta}(Z_{k})\right)^{2},\]

with i.i.d. sampling points \(\{X_{i}\}_{i=1}^{N_{r}}\sim U(\Omega)\), \(\{Y_{j}\}_{j=1}^{N_{d}}\sim U(\Gamma_{D})\) and \(\{Z_{k}\}_{k=1}^{N_{n}}\sim U(\Gamma_{N})\). The empirical loss \(\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}(w_{\theta};\boldsymbol{\gamma}_{N})\) can be minimized using the two-stage procedure as in the 2D case.

#### 3.2.3 SEPINN - Neural networks approximation

There is actually a direct way to resolve the singular part \(S\), i.e., using a DNN to approximate \(\Phi\) in (3.14a). We term the resulting method SEPINN-N. This strategy eliminates the necessity of explicitly knowing the expansion basis and relieves us from lengthy derivations. Thus, compared with SEPINN-C, it is more direct and simpler to implement. The downside is an increase in the number of parameters that need to be learned. Specifically, let \(\mathcal{B}\) to be a DNN set with a fixed architecture (possibly different from \(\mathcal{A}\)) and \(\zeta\) its parameterization. Then a DNN \(\Phi_{\zeta}\in\mathcal{B}\) is employed to approximate \(\Phi\) in (3.17), where the DNN parameters \(\zeta\) are also learned. The splitting is then given by

\[u=w+\Phi_{\zeta}\eta s.\]

Since we cannot guarantee \(\Phi_{\zeta}=0\) on \(\Gamma_{D}\) or \(\partial_{n}\Phi_{\zeta}=0\) on \(\Gamma_{N}\), the boundary conditions of \(w\) have to be modified accordingly (noting \(\partial_{n}(\eta_{\rho}s)=0\)):

\[\begin{cases}-\Delta w=f+\Delta(\Phi_{\zeta}\eta_{\rho}s),&\text{in }\Omega,\\ w=-\Phi_{\zeta}\eta_{\rho}s,&\text{on }\Gamma_{D},\\ \partial_{n}w=-\partial_{n}(\Phi_{\zeta})\eta_{\rho}s,&\text{on }\Gamma_{N}. \end{cases}\]

Like before, we can obtain the following empirical loss

\[\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}(w_{\theta};\Phi_{ \zeta})= \frac{|\Omega|}{N_{r}}\sum_{i=1}^{N_{r}}\left(\Delta w_{\theta}(X _{i})+f(X_{i})+\Delta(\Phi_{\zeta}\eta_{\rho}s)(X_{i})\right)^{2}+\sigma_{d} \frac{|\Gamma_{D}|}{N_{d}}\sum_{j=1}^{N_{d}}(w_{\theta}(Y_{j})+\Phi_{\zeta} \eta_{\rho}s(Y_{j}))^{2}\] \[+\sigma_{n}\frac{|\Gamma_{N}|}{N_{n}}\sum_{k=1}^{N_{n}}\left( \partial_{n}w_{\theta}(Z_{k})+\partial_{n}(\Phi_{\zeta})\eta_{\rho}s(Z_{k}) \right)^{2},\]

with i.i.d. sampling points \(\{X_{i}\}_{i=1}^{N_{r}}\sim U(\Omega)\), \(\{Y_{j}\}_{j=1}^{N_{d}}\sim U(\Gamma_{D})\) and \(\{Z_{k}\}_{k=1}^{N_{n}}\sim U(\Gamma_{N})\). The implementation of SEPINN-N is direct, since both DNNs \(w_{\theta}\) and \(\Psi_{\zeta}\) are learned, and the resulting optimization problem can be minimized using the path-following strategy.

### Modified Helmholtz equation

Now we extend SEPINN to the modified Helmholtz equation equipped with mixed boundary conditions:

\[\begin{cases}-\Delta u+A^{2}u=f,&\text{in }\Omega,\\ \qquad\qquad u=0,&\text{on }\Gamma_{D},\\ \partial_{n}u=0,&\text{on }\Gamma_{N},\end{cases} \tag{3.22}\]

where \(A>0\) is the wave number. It is also known as the screened Poisson problem in the literature.
#### 3.3.1 Two-dimensional case

Following the discussion in Section 3.1, we use the basis \(\{\phi_{j,k}\}_{k=1}^{\infty}\) to expand the solution \(u\), find the singular term and split it from the solution \(u\) of (3.22). Since \(u\) cannot be expressed in terms of an elementary function, the derivation of the singularity splitting differs from that for the Poisson equation. Nonetheless, a similar decomposition holds [20, Theorem 3.3].

**Theorem 3.4**.: _Let \(\mathbf{v}_{j},j=1,2,\cdots,M\), be the vertices of \(\Omega\) whose interior angles \(\omega_{j},j=1,2,\cdots,M\), satisfy (3.3). Then the unique weak solution \(u\in H^{1}(\Omega)\) of (3.22) can be decomposed into_

\[u=w+\sum_{j=1}^{M}\sum_{i\in\mathbb{I}_{j}}\gamma_{j,i}\eta_{\rho_{j}}(r_{j}) \mathrm{e}^{-Ar_{j}}s_{j,i}(r_{j},\theta_{j}),\quad\text{with }w\in H^{2}(\Omega).\]

_Moreover, there holds \(|w|_{H^{2}(\Omega)}+A|w|_{H^{1}(\Omega)}+A^{2}\|w\|_{L^{2}(\Omega)}\leq C\|f \|_{L^{2}(\Omega)}\)._

Once having the singularity solutions \(s_{j,i}\), one can learn the parameters \(\gamma_{j,i}\) and the DNN parameters \(\theta\) of the regular part \(w_{\theta}\) as Algorithm 1. Thus, SEPINN applies equally well to the (modified) Helmholtz case. In contrast, the FEM uses an integral formula and numerical integration when calculating the flux intensity factors [57, 51], which can be rather cumbersome. Further, such extraction formulas are still unavailable for the Helmholtz equation. In SEPINN, training the coefficients directly along with the DNN not only reduces manual efforts, but also increases the computational efficiency.

#### 3.3.2 Three-dimensional case

Following the discussion in Section 3.2, we employ the orthogonal basis \(\{Z_{n}\}_{n=0}^{\infty}\) in Table 3.3 to expand \(u\) and \(f\), cf. (3.11) and (3.12), and substitute them into (3.22):

\[\begin{cases}-\Delta u_{n}+(\xi_{j,n}^{2}+A^{2})u_{n}=f_{n},&\text{in }\Omega_{0},\\ \phantom{-\Delta u_{n}+(\xi_{j,n}^{2}+A^{2})u_{n}=0,\quad\text{on }\Gamma_{D},}\\ \phantom{-\Delta u_{n}+(\xi_{j,n}^{2}+A^{2})u_{n}=0,\quad\text{on }\Gamma_{N}.} \end{cases} \tag{3.23}\]

The only difference of (3.23) from (3.13) lies in the wave number. Upon letting \(\widetilde{A}_{j,n}=(\xi_{j,n}^{2}+A^{2})^{1/2}\), we can rewrite Theorems 3.1 and 3.2 with \(\xi_{j,n}=\widetilde{A}_{j,n}\), and obtain the following result.

**Theorem 3.5**.: _For each \(f\in L^{2}(\Omega)\), let \(u\in H^{1}(\Omega)\) be the unique weak solution to problem (3.22). Then \(u\) can be split into a sum of a regular part \(w\) and a singular part \(S\) as_

\[u=w+S,\quad\text{with }w\in H^{2}(\Omega),\quad S=\sum_{j=1}^{M} \sum_{i\in\mathbb{I}_{j}}S_{j,i}(x_{j},y_{j},z),\] \[S_{j,i}(x_{j},y_{j},z)=\widetilde{\Phi}_{j,i}(r_{j},z)\eta_{\rho _{j}}(r_{j})s_{j,i}(r_{j},\theta_{j}),\quad\text{with }\widetilde{\Phi}_{j,i}(x_{j},y_{j},z)=\frac{1}{2}\gamma_{j,i,0}+\sum_{n=1}^{ \infty}\gamma_{j,i,n}\mathrm{e}^{-\widetilde{A}_{j,n}\tau_{j}}Z_{n}(z).\]

Note that the expressions for the coefficients \(\gamma_{j,i,n}\) differ from that in Section 3.2, which however are not needed for implementing SEPINN and thus not further discussed. Nonetheless, we can still use both SEPINN-C and SEPINN-N to solve 3D Helmholtz problems.

## 4 Error analysis

Now we discuss the error analysis of SEPINN developed in Section 3, following the strategies established in the recent works [33, 46, 31], in order to provide theoretical guarantee of SEPINN. We only analyze the 2D problem in Section 3.1. Let \((w^{*},\gamma^{*})\) be the global minimizer of the loss \(\mathcal{L}_{\mathbf{\sigma}}(w;\gamma)\), cf. (3.9), and the exact solution \(u^{*}\) to problem (1.1) be given by \(u^{*}=w^{*}+\gamma^{*}\eta_{\rho}s.\) Moreover, we assume that \(w^{*}\in H^{3}(\Omega)\) and \(|\gamma^{*}|\leq B\), cf. (3.7); otherwise, we can split out additional singular function(s); see Table 3.1 and the argument in Section 3.1.1.

The following approximation property holds [30, Proposition 4.8].
**Lemma 4.1**.: _Let \(s\in\mathbb{N}\cup\{0\}\) and \(p\in[1,\infty]\) be fixed, and \(v\in W^{k,p}(\Omega)\) with \(k\geq s+1\). Then for any tolerance \(\epsilon>0\), there exists at least one \(v_{\theta}\) of depth \(\mathcal{O}\big{(}\log(d+k)\big{)}\), with \(|\theta|_{\epsilon^{0}}\) bounded by \(\mathcal{O}\big{(}\epsilon^{-\frac{d}{k-s-\mu(\epsilon=2)}}\big{)}\) and \(|\theta|_{\epsilon^{\infty}}\) bounded by \(\mathcal{O}(\epsilon^{-2-\frac{2(d\epsilon/\epsilon+\epsilon+\epsilon+\epsilon (\epsilon=2))+d\epsilon/\epsilon+d}{k-s-\mu(\epsilon=2)}})\), where \(\mu>0\) is arbitrarily small, such that_

\[\|v-v_{\theta}\|_{W^{s,p}(\Omega)}\leq\epsilon.\]

For any \(\epsilon>0\), with \(d=2\), \(k=3\) and \(s=2\), Lemma 4.1 implies that there exists a DNN

\[v_{\theta}\in\mathcal{N}(C,C\epsilon^{-\frac{2}{1-\epsilon}},C\epsilon^{-2- \frac{16+2\epsilon}{1-\epsilon}})=:\mathcal{W}_{\epsilon},\]

such that \(\|w^{*}-v_{\theta}\|_{H^{2}(\Omega)}\leq\epsilon\). Also let \(I_{\gamma}=[-B,B]\). Let \((\widehat{w}_{\theta},\widehat{\gamma})\) be a minimizer of \(\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}(w_{\theta};\gamma)\), cf. (3.10) over \(\mathcal{W}_{\epsilon}\times I_{\gamma}\), and set \(\widehat{u}=\widehat{w}_{\theta}+\widehat{\gamma}\eta_{\rho}s.\) The next lemma gives a decomposition of the error \(\|u^{*}-\widehat{u}\|_{L^{2}(\Omega)}\).

**Lemma 4.2**.: _For any \(\epsilon>0\), let \((\widehat{w}_{\theta},\widehat{\gamma})\in\mathcal{W}_{\epsilon}\times I_{\gamma}\) be a minimizer to the loss \(\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}(w_{\theta};\gamma)\). Then there exists a constant \(c(\boldsymbol{\sigma})\) such that_

\[\|u^{*}-\widehat{u}\|_{L^{2}(\Omega)}^{2}\leq c(\boldsymbol{\sigma})\mathcal{ L}_{\boldsymbol{\sigma}}(\widehat{w}_{\theta};\widehat{\gamma})\leq c( \boldsymbol{\sigma})\Big{(}\epsilon^{2}+\sup_{(w_{\theta},\gamma)\in\mathcal{ W}_{\epsilon}\times I_{\gamma}}|\mathcal{L}_{\boldsymbol{\sigma}}(w_{ \theta};\gamma)-\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}(w_{\theta}; \gamma)|\Big{)}.\]

Proof.: For any \((w_{\theta},\gamma)\in\mathcal{W}_{\epsilon}\times I_{\gamma}\), let \(u=w_{\theta}+\gamma\eta_{\rho}s\), and \(e=u^{*}-u\). By the trace theorem, we have

\[\mathcal{L}_{\boldsymbol{\sigma}}(w_{\theta};\gamma)=\|\Delta e\|_{L^{2}( \Omega)}^{2}+\sigma_{d}\|e\|_{L^{2}(\Gamma_{D})}^{2}+\sigma_{n}\|\partial_{n}e \|_{L^{2}(\Gamma_{N})}^{2}\leq c(\boldsymbol{\sigma})\|e\|_{H^{2}(\Omega)}^{2}.\]

To treat the nonzero boundary conditions of \(w_{\theta}\), we define the harmonic extension \(\zeta\) by

\[\begin{cases}-\Delta\zeta=0,&\text{in }\Omega,\\ \quad\zeta=w_{\theta},&\text{on }\Gamma_{D},\\ \quad\partial_{n}\zeta=\partial_{n}w_{\theta},&\text{on }\Gamma_{N}.\end{cases}\]

Then the following elliptic regularity estimate holds [9, Theorem 4.2, p. 870]

\[\|\zeta\|_{L^{2}(\Omega)}\leq c\big{(}\|w_{\theta}\|_{L^{2}(\Gamma_{D})}+\| \partial_{n}w_{\theta}\|_{L^{2}(\Gamma_{N})}\big{)}. \tag{4.1}\]

Let \(\tilde{e}=e+\zeta\). Then it satisfies

\[\begin{cases}-\Delta\tilde{e}=-\Delta e,&\text{in }\Omega,\\ \quad\tilde{e}=0,&\text{on }\Gamma_{D},\\ \quad\partial_{n}\tilde{e}=0,&\text{on }\Gamma_{N}.\end{cases}\]

Since \(\Delta e\in L^{2}(\Omega)\), the standard energy argument and Poincare inequality imply

\[\|\tilde{e}\|_{L^{2}(\Omega)}\leq c\|\nabla\tilde{e}\|_{L^{2}(\Omega)}\leq c \|\Delta e\|_{L^{2}(\Omega)}.\]

This, the stability estimate (4.1) and the triangle inequality lead to

\[\|e\|_{L^{2}(\Omega)}^{2}\leq c\big{(}\|\tilde{e}\|_{L^{2}(\Omega)}+\|\zeta\|_ {L^{2}(\Omega)}\big{)}\leq c\big{(}\|\Delta e\|_{L^{2}(\Omega)}^{2}+\|w_{ \theta}\|_{L^{2}(\Gamma_{D})}^{2}+\|\partial_{n}w_{\theta}\|_{L^{2}(\Gamma_{N })}^{2}\big{)}\leq c(\boldsymbol{\sigma})\mathcal{L}_{\boldsymbol{\sigma}}(w_{ \theta};\gamma).\]

This proves the first inequality of the lemma. Next, by Lemma 4.1 and the assumption \(w^{*}\in H^{3}(\Omega)\), there exists \(w_{\theta}\in\mathcal{W}_{\epsilon}\) such that \(\|w_{\theta}-w^{*}\|_{H^{2}(\Omega)}\leq\epsilon\). Let \(\bar{u}=w_{\theta}+\gamma^{*}\eta_{\rho}s\). Then we derive

\[\mathcal{L}_{\boldsymbol{\sigma}}(w_{\bar{\theta}};\gamma^{*}) =\|\Delta(w_{\bar{\theta}}-w^{*})\|_{L^{2}(\Omega)}^{2}+\sigma_{d}\|w_{ \bar{\theta}}-w^{*}\|_{L^{2}(\Gamma_{D})}^{2}+\sigma_{n}\|\partial_{n}(w_{ \bar{\theta}}-w^{*})\|_{L^{2}(\Gamma_{N})}^{2}\] \[\leq c(\boldsymbol{\sigma})\|w_{\bar{\theta}}-w^{*}\|_{H^{2}( \Omega)}^{2}\leq c(\boldsymbol{\sigma})\epsilon^{2}.\]

Thus by the minimizing property of \((\widehat{w}_{\theta},\widehat{\gamma})\), we arrive at

\[\mathcal{L}_{\boldsymbol{\sigma}}(\widehat{w}_{\theta};\widehat{ \gamma}) =\mathcal{L}_{\boldsymbol{\sigma}}(\widehat{w}_{\theta};\widehat{ \gamma})-\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}(\widehat{w}_{\theta}; \widehat{\gamma})+\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}(\widehat{w}_{\theta}; \widehat{\gamma})-\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}(w_{\bar{\theta}}; \gamma^{*})+\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}(w_{\bar{\theta}}; \gamma^{*})-\mathcal{L}_{\boldsymbol{\sigma}}(w_{\bar{\theta}};\gamma^{*})+ \mathcal{L}_{\boldsymbol{\sigma}}(w_{\bar{\theta}};\gamma^{*})\] \[\leq\big{|}\mathcal{L}_{\boldsymbol{\sigma}}(\widehat{w}_{\theta}; \widehat{\gamma})-\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}(\widehat{w}_{ \theta};\widehat{\gamma})\big{|}+\big{|}\widehat{\mathcal{L}}_{\boldsymbol{ \sigma}}(w_{\bar{\theta}};\gamma^{*})-\mathcal{L}_{\boldsymbol{\sigma}}(w_{\bar{ \theta}};\gamma^{*})\big{|}+\mathcal{L}_{\boldsymbol{\sigma}}(w_{\bar{\theta}}; \gamma^{*})\] \[\leq 2\sup_{(w_{\theta},\gamma)\in\mathcal{W}_{\epsilon}\times I_{ \gamma}}|\mathcal{L}_{\boldsymbol{\sigma}}(w_{\theta};\gamma)-\widehat{ \mathcal{L}}_{\boldsymbol{\sigma}}(w_{\theta};\gamma)|+c(\boldsymbol{\sigma}) \epsilon^{2}.\]

This completes the proof of the lemma.
Next, we bound the statistical error \(\mathcal{E}_{stat}=\sup_{(w_{\theta},\gamma)\in\mathcal{W}_{\kappa}\times I_{\gamma}} \big{|}\mathcal{L}_{\mathbf{\sigma}}(w_{\theta};\gamma)-\widehat{\mathcal{L}}_{\bm {\sigma}}(w_{\theta};\gamma)\big{|}\), which arises from approximating the integrals by Monte Carlo. By the triangle inequality, we have the splitting

\[\mathcal{E}_{stat}\leq \sup_{(w_{\theta},\gamma)\in\mathcal{W}_{\kappa}\times I_{\gamma}} |\Omega|\Big{|}\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}h_{r}(X_{i};w_{\theta},\gamma) -\mathbb{E}_{X}(h_{r}(X;w_{\theta},\gamma))\Big{|}\] \[+\sup_{w_{\theta}\in\mathcal{W}_{\epsilon}}\sigma_{d}|\Gamma_{D}| \Big{|}\frac{1}{N_{d}}\sum_{j=1}^{N_{d}}h_{d}(Y_{j};w_{\theta})-\mathbb{E}_{Y }(h_{d}(Y;w_{\theta}))\Big{|} \tag{4.2}\] \[+\sup_{w_{\theta}\in\mathcal{W}_{\epsilon}}\sigma_{n}|\Gamma_{N} |\Big{|}\frac{1}{N_{n}}\sum_{k=1}^{N_{n}}h_{n}(Z_{k};w_{\theta})-\mathbb{E}_{ Y}(h_{n}(Z;w_{\theta}))\Big{|},\]

with \(h_{r}(x;w_{\theta},\gamma)=(\Delta w_{\theta}+f+\gamma\Delta(\eta_{\theta}s) )^{2}(x)\) for \(x\in\Omega\), \(h_{d}(y;w_{\theta})=|w_{\theta}(y)|^{2}\) for \(y\in\Gamma_{D}\) and \(h_{n}(z;w_{\theta})=|\partial_{n}w_{\theta}(z)|^{2}\) for \(z\in\Gamma_{N}\). Thus, we define the following three function classes

\[\mathcal{H}_{r}=\{h_{r}(w_{\theta},\gamma):w_{\theta}\in\mathcal{W}_{\epsilon },\gamma\in I_{\gamma}\},\ \ \mathcal{H}_{d}=\{h_{d}(w_{\theta}):w_{\theta}\in\mathcal{W}_{\epsilon}\}\ \ \text{and}\ \ \mathcal{H}_{n}=\{h_{n}(w_{\theta}):w_{\theta}\in\mathcal{W}_{\epsilon}\}.\]

To bound the errors in (4.2), we employ Rademacher complexity [1, 6], which measures the complexity of a collection of functions by the correlation between function values with Rademacher random variables.

**Definition 4.1**.: _Let \(\mathcal{F}\) be a real-valued function class defined on the domain \(D\) and \(\xi=\{\xi_{j}\}_{j=1}^{n}\) be i.i.d. samples from the distribution \(\mathcal{U}(D)\). Then the Rademacher complexity \(\mathfrak{R}_{n}(\mathcal{F})\) of \(\mathcal{F}\) is defined by_

\[\mathfrak{R}_{n}(\mathcal{F})=\mathbb{E}_{\xi,\omega}\bigg{[}\sup_{f\in \mathcal{F}}\ n^{-1}\bigg{|}\ \sum_{j=1}^{n}\omega_{j}f(\xi_{j})\ \bigg{|}\bigg{]},\]

_where \(\omega=\{\omega_{j}\}_{j=1}^{n}\) are i.i.d Rademacher random variables with probability \(P(\omega_{j}=1)=P(\omega_{j}=-1)=\frac{1}{2}\)._

Then we have the following PAC-type generalization bound [50, Theorem 3.1].

**Lemma 4.3**.: _Let \(X_{1},\dots,X_{n}\) be a set of i.i.d. random variables, and let \(\mathcal{F}\) be a function class defined on \(D\) such that \(\sup_{f\in\mathcal{F}}\|f\|_{L^{\infty}(D)}\leq M_{\mathcal{F}}<\infty\). Then for any \(\tau\in(0,1)\), with probability at least \(1-\tau\):_

\[\sup_{f\in\mathcal{F}}\bigg{|}n^{-1}\sum_{j=1}^{n}f(X_{j})-\mathbb{E}[f(X)] \bigg{|}\leq 2\mathfrak{R}_{n}(\mathcal{F})+2M_{\mathcal{F}}\sqrt{\frac{\log \frac{1}{\tau}}{2n}}.\]

To apply Lemma 4.3, we bound Rademacher complexities of the function classes \(\mathcal{H}_{r}\), \(\mathcal{H}_{d}\) and \(\mathcal{H}_{n}\). This follows from Dudley's formula in Lemma 4.6. The next lemma gives useful boundedness and Lipschitz continuity of the DNN function class in terms of \(\theta\); see [34, Lemma 3.4 and Remark 3.3] and [35, Lemma 5.3]. The estimates also hold when \(L^{\infty}(\Omega)\) is replaced with \(L^{\infty}(\Gamma_{D})\) or \(L^{\infty}(\Gamma_{N})\).

**Lemma 4.4**.: _Let \(L\), \(W\) and \(R\) be the depth, width and maximum weight bound of a DNN function class \(\mathcal{W}\), with \(N_{\theta}\) nonzero weights. Then for any \(v_{\theta}\in\mathcal{W}\), the following estimates hold_

* \(\|v_{\theta}\|_{L^{\infty}(\Omega)}\leq WR\)_,_ \(\|v_{\theta}-v_{\tilde{\theta}}\|_{L^{\infty}(\Omega)}\leq 2LW^{L}R^{L-1}| \theta-\tilde{\theta}|_{\ell^{\infty}}\)_;_
* \(\|\nabla v_{\theta}\|_{L^{\infty}(\Omega;\mathbb{R}^{d})}\leq\sqrt{2}W^{L-1}R^{L}\)_,_ \(\|\nabla(v_{\theta}-v_{\tilde{\theta}})\|_{L^{\infty}(\Omega;\mathbb{R}^{d})} \leq\sqrt{2}L^{2}W^{2L-2}R^{2L-2}|\theta-\tilde{\theta}|_{\ell^{\infty}}\)_;_
* \(\|\Delta v_{\theta}\|_{L^{\infty}(\Omega)}\leq 2LW^{2L-2}R^{2L}\)_,_ \(\|\Delta(v_{\theta}-v_{\tilde{\theta}})\|_{L^{\infty}(\Omega)}\leq 8N_{\theta}L^{2}W^{3L-3}R^{3L-3}| \theta-\tilde{\theta}|_{\ell^{\infty}}\)_._

Lemma 4.4 implies boundedness and Lipschitz continuity of functions in \(\mathcal{H}_{r}\), \(\mathcal{H}_{d}\) and \(\mathcal{H}_{n}\).

**Lemma 4.5**.: _There exists a constant \(c\) depending on \(\|f\|_{L^{\infty}(\Omega)}\), \(\|\Delta(\eta_{\rho}s)\|_{L^{\infty}(\Omega)}\) and \(B\) such that_

\[\|h(w_{\theta},\gamma)\|_{L^{\infty}(\Omega)}\leq cL^{2}W^{4L-4}R^{L},\quad\forall h\in\mathcal{H}_{r},\] \[\|h(w_{\theta})\|_{L^{\infty}(\Gamma_{D})}\leq cW^{2R}R^{2},\quad\forall h\in\mathcal{H}_{d},\] \[\|h(w_{\theta})\|_{L^{\infty}(\Gamma_{N})}\leq cW^{2L-2}R^{2L},\quad\forall h\in\mathcal{H}_{n}.\]

_Moreover, the following Lipschitz continuity estimates hold_

\[\|h(w_{\theta},\gamma)-\tilde{h}(w_{\tilde{\theta}},\gamma)\|_{L^{ \infty}(\Omega)}\leq cN_{\theta}L^{3}W^{5L-5}R^{5L-3}(|\theta-\tilde{\theta}|_{ \ell^{\infty}}+|\gamma-\tilde{\gamma}|),\quad\forall h,\tilde{h}\in\mathcal{H}_{r},\] \[\|h(w_{\theta})-h(w_{\tilde{\theta}})\|_{L^{\infty}(\Gamma_{D})}\leq cLW^{L+1}R^{L}| \theta-\tilde{\theta}|_{\ell^{\infty}},\quad\forall h,\tilde{h}\in\mathcal{H}_{d},\] \[\|h(w_{\theta})-\tilde{h}(w_{\tilde{\theta}})\|_{L^{\infty}( \Gamma_{N})}\leq cL^{2}W^{3L-3}R^{3L-2}|\theta-\tilde{\theta}|_{\ell^{\infty}}, \quad\forall h,\tilde{h}\in\mathcal{H}_{n}.\]
Proof.: The uniform bounds follow directly from Lemma 4.5. The Lipschitz estimates follow similarly, and we show it only for the set \(H_{r}\). Indeed, for any \(h=h_{r}(w_{\theta},\gamma),\tilde{h}=h_{r}(w_{\tilde{\theta}},\tilde{\gamma}) \in\mathcal{H}_{r}\), we have

\[\|h-\tilde{h}\|_{L^{\infty}(\Omega)}=\|(\Delta w_{\theta}+f+\gamma\Delta(\eta_{ \rho}s))^{2}-(\Delta w_{\tilde{\theta}}+f+\tilde{\gamma}\Delta(\eta_{\rho}s))^ {2}\|_{L^{\infty}(\Omega)}.\]

By Lemma 4.5, we have

\[\|(\Delta w_{\theta}+f+\gamma\Delta(\eta_{\rho}s))+(\Delta w_{ \tilde{\theta}}+f+\tilde{\gamma}\Delta(\eta_{\rho}s))\|_{L^{\infty}(\Omega)}\] \[\leq 2(2LW^{2L-2}R^{2L}+\|f\|_{L^{\infty}(\Omega)}+B\|\Delta(\eta_{ \rho}s)\|_{L^{\infty}(\Omega)}\leq cLW^{2L-2}R^{2L},\] \[\|\Delta(w_{\theta}-w_{\tilde{\theta}})\|_{L^{\infty}(\Omega)}+| \gamma-\tilde{\gamma}|\|\Delta(\eta_{\rho}s))^{2}\|_{L^{\infty}(\Omega)}\] \[\leq 8N_{\theta}L^{2}W^{3L-3}R^{3L-3}|\theta-\tilde{\theta}|_{\ell \infty}+|\gamma-\tilde{\gamma}||\Delta(\eta_{\rho}s)\|_{L^{\infty}(\Omega)}\] \[\leq cN_{\theta}L^{2}W^{3L-3}R^{3L-3}(|\theta-\tilde{\theta}|_{\ell ^{\infty}}+|\gamma-\tilde{\gamma}|).\]

Combining the last three estimates yield the bound on \(\|h-\tilde{h}\|_{L^{\infty}(\Omega)}\). 

Next, we state Dudley's lemma ([47, Theorem 9], [59, Theorem 1.19]), which bounds Rademacher complexities using covering number.

**Definition 4.2**.: _Let \((\mathcal{M},\rho)\) be a metric space of real valued functions, and \(\mathcal{G}\subset\mathcal{M}\). A set \(\{x_{i}\}_{i=1}^{n}\subset\mathcal{G}\) is called an \(\epsilon\)-cover of \(\mathcal{G}\) if for any \(x\in\mathcal{G}\), there exists a \(x_{i}\) such that \(\rho(x,x_{i})\leq\epsilon\). The \(\epsilon\)-covering number \(\mathcal{C}(\mathcal{G},\rho,\epsilon)\) is the minimum cardinality among all \(\epsilon\)-covers of \(\mathcal{G}\) with respect to \(\rho\)._

**Lemma 4.6** (Dudley's lemma).: _Let \(M_{\mathcal{F}}:=\sup_{f\in\mathcal{F}}\|f\|_{L^{\infty}(\Omega)}\), and \(\mathcal{C}(\mathcal{F},\|\cdot\|_{L^{\infty}(\Omega)},\epsilon)\) be the covering number of \(\mathcal{F}\). Then the Rademacher complexity \(\mathfrak{R}_{n}(\mathcal{F})\) is bounded by_

\[\mathfrak{R}_{n}(\mathcal{F})\leq\inf_{0<s<M_{\mathcal{F}}}\left(4s\ +\ 12n^{-\frac{1}{2}}\int_{s}^{M_{\mathcal{F}}}\left(\log\mathcal{C}(\mathcal{F}, \|\cdot\|_{L^{\infty}(\Omega)},\epsilon)\right)^{\frac{1}{2}}\,\mathrm{d} \epsilon\right).\]

**Lemma 4.7**.: _For any small \(\tau\), with probability at least \(1-3\tau\), there holds_

\[\sup_{(w_{\theta},\gamma)\in\mathcal{W}_{\epsilon}\times I_{\gamma}}|\mathcal{ L}_{\mathbf{\sigma}}(w_{\theta};\gamma)-\widehat{\mathcal{L}}_{\mathbf{\sigma}}(w_{ \theta};\gamma)|\leq c(e_{r}+\sigma_{d}e_{d}+\sigma_{n}e_{n}),\]

_where \(c\) depends on \(\|f\|_{L^{\infty}(\Omega)}\) and \(\|\Delta(\eta_{\rho}s)\|_{L^{\infty}(\Omega)}\), and \(e_{r}\), \(e_{d}\) and \(e_{n}\) are respectively defined by_

\[e_{r} \leq c\frac{L^{2}R^{4L}N_{\theta}^{4L-4}\big{(}N_{\theta}^{\frac{1}{2 }}\left(\log^{\frac{1}{2}}R+\log^{\frac{1}{2}}N_{\theta}+\log^{\frac{1}{2}}N_{r }\right)+\log^{\frac{1}{2}}\frac{1}{\tau}\big{)}}{\sqrt{N_{r}}},\] \[e_{d} \leq c\frac{R^{2}N_{\theta}^{2}\big{(}N_{\theta}^{\frac{1}{2}}\left( \log^{\frac{1}{2}}R+\log^{\frac{1}{2}}N_{\theta}+\log^{\frac{1}{2}}N_{d} \right)+\log^{\frac{1}{2}}\frac{1}{\tau}\big{)}}{\sqrt{N_{d}}},\] \[e_{n} \leq c\frac{R^{2L}N_{\theta}^{2L-2}\big{(}N_{\theta}^{\frac{1}{2}} \left(\log^{\frac{1}{2}}R+\log^{\frac{1}{2}}N_{\theta}+\log^{\frac{1}{2}}N_{n }\right)+\log^{\frac{1}{2}}\frac{1}{\tau}\big{)}}{\sqrt{N_{n}}}.\]

Proof.: Fix \(m\in\mathbb{N}\), \(R\in[1,\infty)\), \(\epsilon\in(0,1)\), and \(B_{R}:=\{x\in\mathbb{R}^{m}:\ |x|_{\ell^{\infty}}\leq R\}\). Then by [16, Prop. 5],

\[\log\mathcal{C}(B_{R},|\cdot|_{\ell^{\infty}},\epsilon)\leq m\log(4R\epsilon^{-1}).\]

The Lipschitz continuity estimates in Lemmas 4.4 and 4.5 imply

\[\log\mathcal{C}(\mathcal{H}_{r},\|\cdot\|_{L^{\infty}(\Omega)},\epsilon)\leq \log\mathcal{C}(\Theta\times I_{\gamma},\max(|\cdot|_{\ell^{\infty}},|\cdot|), \Lambda_{r}^{-1}\epsilon)\leq cN_{\theta}\log(4R\Lambda_{r}\epsilon^{-1}),\]

with \(\Lambda_{r}=cN_{\theta}L^{3}W^{5L-5}R^{5L-3}\). Moreover, by Lemma 4.5, we have \(M_{r}\equiv M_{\mathcal{H}_{r}}=cL^{2}W^{4L-4}R^{4L}\). Then setting \(s=n^{-\frac{1}{2}}\) in Lemma 4.6 and using the estimates \(1\leq R\), \(1\leq L\) and \(1\leq W\leq N_{\theta}\), \(1\leq L\leq c\log 5\) (with \(d=2\) and \(k=3\)), cf. Lemma 4.1, lead to

\[\mathfrak{R}_{n}(\mathcal{H}_{r})\leq 4n^{-\frac{1}{2}}+12n^{-\frac{1}{2}}\int_{n^{-\frac{1}{2}}}^{M_{r}} \big{(}cN_{\theta}\log(4R\Lambda_{r}\epsilon^{-1})\big{)}^{\frac{1}{2}}\,\mathrm{d}\epsilon\] \[\leq 4n^{-\frac{1}{2}}+12n^{-\frac{1}{2}}M_{r}\big{(}cN_{\theta}\log( 4R\Lambda_{r}n^{\frac{1}{2}})\big{)}^{\frac{1}{2}}\] \[\leq 4n^{-\frac{1}{2}}+cn^{-\frac{1}{2}}W^{4L-4}R^{4L}N_{\theta}^{ \frac{1}{2}}\left(\log^{\frac{1}{2}}R+\log^{\frac{1}{2}}\Lambda_{r}+\log^{\frac{1} {2}}n\right)\]
\[\leq cn^{-\frac{1}{2}}R^{4L}N_{\theta}^{4L-\frac{1}{2}}\big{(}\log^{\frac{1}{2}}R +\log^{\frac{1}{2}}N_{\theta}+\log^{\frac{1}{2}}n\big{)}.\]

Similarly, repeating the preceding argument leads to

\[\mathfrak{R}_{n}(\mathcal{H}_{d}) \leq cn^{-\frac{1}{2}}R^{2}N_{\theta}^{\frac{5}{2}}(\log^{\frac{1}{2 }}R+\log^{\frac{1}{2}}N_{\theta}+\log^{\frac{1}{2}}n),\] \[\mathfrak{R}_{n}(\mathcal{H}_{n}) \leq cn^{-\frac{1}{2}}R^{2L}N_{\theta}^{2L-\frac{3}{2}}(\log^{ \frac{1}{2}}R+\log^{\frac{1}{2}}N_{\theta}+\log^{\frac{1}{2}}n).\]

Finally, the desired result follows from the PAC-type generalization bound in Lemma 4.3. 

Then combining Lemma 4.2 with Lemma 4.7 yields the following error estimate. Thus, by choosing the numbers \(N_{r}\), \(N_{d}\) and \(N_{n}\) of sampling points sufficiently large, the \(L^{2}(\Omega)\) error of the SEPINN approximation can be made about \(O(\epsilon^{2})\).

**Theorem 4.1**.: _Fix a tolerance \(\epsilon>0\), and let \((\widehat{w}_{\theta},\widehat{\gamma})\in\mathcal{W}_{\epsilon}\times I_{\gamma}\) be a minimizer to the empirical loss \(\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}(w_{\theta};\gamma)\) in (3.10). Then for any small \(\tau\), with the statistical errors \(e_{r}\), \(e_{d}\) and \(e_{n}\) from Lemma 4.7, we have with probability at least \(1-3\tau\), the following error estimate holds_

\[\|u^{*}-\widehat{u}\|_{L^{2}(\Omega)}^{2}\leq c(\boldsymbol{\sigma})\big{(} \epsilon^{2}+e_{r}+e_{d}+e_{n}\big{)}.\]

## 5 Numerical experiments

Now we present numerical examples (with zero Dirichlet and / or Neumann boundary condition) to illustrate SEPINN, and compare it with existing PINN type solvers. In the training, \(N_{r}=10,000\) points in the domain \(\Omega\) and \(N_{b}=800\) points on the boundary \(\partial\Omega\) are selected uniformly at random to form the empirical loss \(\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}\), unless otherwise specified. In the path-following (PF) strategy, we take an increasing factor \(q=1.5\). All numerical experiments were conducted on a personal laptop (Windows 10, with RAM 8.0GB, Intel(R) Core(TM) i7-10510U CPU, 2.3 GHz), with Python 3.9.7, with the software framework PyTorch. The gradient of the DNN output \(w_{\theta}(x)\) with respect to the input \(x\) (i.e., spatial derivative) and that of the loss \(\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}\) to \(\theta\) are computed via automatic differentiation [7], as implemented by torch.autograd.

For SEPINN and SEPINN-C (based on cutoff), we minimize the loss \(\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}\) in two stages: first determine the coefficients \(\boldsymbol{\gamma}_{N}\), and then reduce the boundary error and refine the DNN approximation \(w_{\theta}\) (of the regular part \(w\)) by the PF strategy. We use different optimizers at these two stages. First, we minimize the loss \(\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}(w_{\theta};\boldsymbol{\gamma})\) in both \(\theta\) and \(\boldsymbol{\gamma}\) using Adam [37] (from the SciPy library), with the default setting (tolerance: 1.0e-8, no box constraint, maximum iteration number: 1000); then, we minimize \(\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}(w_{\theta};\widehat{\gamma}^{*})\) (fixing \(\boldsymbol{\gamma}\) at \(\widehat{\gamma}^{*}\) fixed) using limited memory BFGS (L-BFGS) [11], with the default setting (tolerance: 1.0e-9, no box constraint, maximum iteration number: 2500). For SEPINN-N, we employ only L-BFGS [11]. To measure the accuracy of an approximation \(\hat{w}\) of \(w^{*}\), we use the relative \(L^{2}(\Omega)\)-error \(e=\|w^{*}-\hat{w}\|_{L^{2}(\Omega)}/\|w^{*}\|_{L^{2}(\Omega)}\). The stopping condition of the PF strategy is set to \(e<1.00\)e-3 and \(\sigma_{d}^{(k)},\sigma_{n}^{(k)}\leq\sigma^{*}\), for some fixed \(\sigma^{*}>0\). The first condition ensures that \(\hat{w}\) can achieve the desired accuracy and the second one terminates the iteration after a finite number of loops. The detailed hyper-parameter setting of the PF strategy is listed in Table 5.1. The Python code for reproducing the numerical experiments will be made available at [https://github.com/hhjc-web/SEPINN.git](https://github.com/hhjc-web/SEPINN.git).

First, we showcase the approach on an L-shape domain [15, Example 5.2].

**Example 5.1**.: _The domain \(\Omega=(-1,1)^{2}\backslash\left([0,1)\times(-1,0]\right)\). Set \(\rho=1\) and \(R=\frac{1}{2}\) in (3.4), the source_

\[f=\begin{cases}\sin(2\pi x)\left[2\pi^{2}\left(y^{2}+2y\right)\left(y^{2}-1 \right)-\left(6y^{2}+6y-1\right)\right]-\Delta\left(\eta_{\rho}s\right),&-1 \leq y\leq 0,\\ \sin(2\pi x)\left[2\pi^{2}\left(-y^{2}+2y\right)\left(y^{2}-1\right)-\left(-6 y^{2}+6y+1\right)\right]-\Delta\left(\eta_{\rho}s\right),&0\leq y\leq 1,\end{cases}\]

_with the singular function \(s=r^{\frac{2}{3}}\sin\left(\frac{2\theta}{3}\right)\), and \(\Gamma_{D}=\partial\Omega\). The exact solution \(u\) of the problem is given by \(u=w+\eta_{\rho}s\), with the regular part \(w\) given by_

\[w=\begin{cases}\sin(2\pi x)\left(\frac{1}{2}y^{2}+y\right)\left(y^{2}-1\right),& -1\leq y\leq 0,\\ \sin(2\pi x)\left(-\frac{1}{2}y^{2}+y\right)\left(y^{2}-1\right),&0\leq y\leq 1.\end{cases}\]
Here the regular part \(w\) belongs to \(H^{2}(\Omega)\) but not to \(H^{3}(\Omega)\), and \(u\) lies in \(H^{1}(\Omega)\) but not in \(H^{2}(\Omega)\). In SEPINN, we employ a 2-20-20-20-1 DNN (3 hidden layers, each having 20 neurons). The first stage of the PF strategy gives an estimate \(\widehat{\gamma}^{*}=0.9958\), and the final prediction error \(e\) after the second stage is 3.43e-3. Fig. 1 shows that the pointwise error of the SEPINN approximation is small and the accuracy around the singularity at the reentrant corner is excellent. In contrast, applying PINN and DRM directly fails to yield satisfactory results near the reentrant corner because of the presence of the singular term \(r^{\frac{7}{3}}\sin\frac{2}{3}\theta\), consistent with the approximation theory of DNNs to singular functions [30]. DRM shows larger errors over the whole domain, not just in the vicinity of the corner. By adaptively adjusting the empirical loss, SAPINN and FIPINN can yield more accurate approximations than PINN, but the error around the singularity is still large. Numerically, FIPINN can adaptively add sampling points near the corner but does not give high concentration, which limits the accuracy of the final DNN approximation. This shows the advantages of explicitly incorporating analytic insights into SEPINN.

\begin{table}
\begin{tabular}{|c|c|c c c c|} \hline Example & method & \(\sigma^{*}\) & \(\boldsymbol{\sigma}^{(1)}\) & \(\boldsymbol{\sigma}^{(K)}\) & \(\widehat{\gamma}^{*}\) & \(e\) \\ \hline
5.1 & SEPINN & 800 & 100 & 759.4 & 0.9958 & 3.43e-3 \\ \hline
5.2 & SEPINN & 800 & (100,100) & (759.4,759.4) & 1.0023 & 3.48e-3 \\ \hline
5.3 & SEPINN-C & 7000 & 2000 & 6750 & Table 5.2 & 3.93e-2 \\  & SEPINN-N & 7000 & 400 & 6834.4 & – – & 2.05e-2 \\ \hline
5.4 & SEPINN-C & 1000 & (100,100) & (759.4,759.4) & Table 5.4 & 2.08e-2 \\  & SEPINN-N & 4000 & (100,400) & (759.4,3037.5) & – – & 3.83e-2 \\ \hline
5.5 & SEPINN & 800 & 100 & 759.4 & 0.9939 & 1.02e-2 \\  & SEPINN-C & 2000 & 400 & 759.4 & Table 5.5 & 3.30e-2 \\  & SEPINN-N & 4000 & 400 & 3037.5 & – & 2.55e-2 \\ \hline
5.6 & SEPINN & 600 & 50 & 569.5 & Table 5.6 & – – \\ \hline \end{tabular}
\end{table}
Table 5.1: The hyper-parameters for the PF strategy for SEPINN (2D) and SEPINN-C and SEPINN-N (3D) for Examples 5.1-5.5. The notation \(\widehat{\gamma}^{*}\) denotes the estimated stress intensity factor, and \(e\) the prediction error.

Figure 1: The numerical approximations for Example 5.1 by SEPINN, SAPINN, FIPINN, PINN and DRM. From top to bottom: exact solution, DNN approximation and pointwise error.
To shed further insights into the methods, we show in Fig. 2 the training dynamics of the empirical loss \(\widehat{\mathcal{L}}\) and relative error \(e\), where \(i\) denotes the total iteration index along with the PF loops. For all methods, the loss \(\widehat{\mathcal{L}}\) and error \(e\) both decay steadily as the iteration proceeds, indicating stable convergence of the optimizer, but SEPINN enjoys the fastest decay and smallest error \(e\), due to the improved regularity of \(w^{*}\). The final error \(e\) saturates at around \(10^{-3}\) for SEPINN and \(10^{-2}\) for PINN, SAPINN and FIPINN, but only \(10^{-1}\) for DRM. Indeed the accuracy of neural PDE solvers tends to stagnate at a level of \(10^{-2}\sim 10^{-3}\)[53, 60, 62, 17], which differs markedly from more traditional solvers.

We next investigate a mixed boundary value problem, adapted from [14, Example 1].

**Example 5.2**.: _The domain \(\Omega\) is the unit square \(\Omega=(0,1)^{2}\), \(\Gamma_{N}=\{(x,0):x\in(0,\frac{1}{2})\}\) and \(\Gamma_{D}=\partial\Omega\backslash\Gamma_{N}\). The singular function \(s=r^{\frac{1}{2}}\sin\frac{\theta}{2}\) in the local polar coordinate \((r,\theta)\) at \((\frac{1}{2},0)\). Set \(\rho=1\) and \(R=\frac{1}{4}\) in (3.4), the source \(f=-\sin(\pi x)(-\pi^{2}y^{2}(y-1)+6y-2)-\Delta(\eta_{\rho}s).\) The exact solution \(u=\sin(\pi x)y^{2}(y-1)+\eta_{\rho}s\), with the regular part \(w=\sin(\pi x)y^{2}(y-1)\) and stress intensity factor \(\gamma=1\)._

This problem has a geometric singularity at the point \((\frac{1}{2},0)\), where the boundary condition changes from Dirichlet to Neumann with an interior angle \(\omega=\pi\). We employ a 3-10-10-10-1 DNN (with 3 hidden layers, each having 10 neurons). The first stage of the PF strategy gives an estimate \(\widehat{\gamma}^{*}=1.0023\), and the prediction error \(e\) after the second stage is 3.48e-3. The singularity at the crack point is accurately resolved by SEPINN, cf. Fig. 3, which shows also the approximations by DRM and other PINN techniques. The overall solution accuracy is very similar to Example 5.1, and SEPINN achieves the

Figure 3: The numerical approximations of Example 5.2 by the proposed SEPINN (error: 3.48e-3), SAPINN (error: 4.40e-2), FIPINN (error: 3.65e-2), PINN (error: 7.33e-2) and DRM (error: 1.86e-1). From top to bottom: analytic solution, DNN approximation and pointwise error.

Figure 2: Training dynamics for SEPINN and benchmark methods: (a) the decay of the empirical loss \(\widehat{\mathcal{L}}\) versus the iteration index \(i\) (counted along the path-following trajectory) and (b) the error \(e\) versus \(i\).
smallest error. SAPINN and FIPINN improve the standard PINN, but still suffer from large errors near the singularity.

The next example is a 3D Poisson equation adapted from [52, Example 1].

**Example 5.3**.: _Let \(\Omega_{0}=(-1,1)^{2}\backslash\ ([0,1)\times(-1,0])\), and the domain \(\Omega=\Omega_{0}\times(-1,1)\). Define_

\[\Phi(r,z)=-2\arctan\frac{\mathrm{e}^{-\pi r}\sin\pi z}{1+\mathrm{e}^{-\pi r} \cos\pi z}=2\sum_{n=1}^{\infty}\frac{(-1)^{n}\mathrm{e}^{-\pi nr}\sin n\pi z} {n}, \tag{5.1}\]

_and set \(\rho=1\) and \(R=\frac{1}{2}\) in (3.4), the source \(f=6x(y-y^{3})(1-z^{2})+6y(x-x^{3})(1-z^{2})+2(y-y^{3})(x-x^{3})-\Delta(\Phi\eta _{\rho}s)\), with the singular function \(s=r^{\frac{3}{3}}\sin(\frac{2\Phi}{3})\), and a zero Dirichlet boundary condition. The exact solution \(u\) is given by \(u=(x-x^{3})(y-y^{3})(1-z^{2})+\Phi\eta_{\rho}s\)._

It follows from (5.1) that the coefficients \(\gamma_{n}^{*}\) are given by \(\gamma_{0}^{*}=0\) and \(\gamma_{n}^{*}=(-1)^{n}\frac{2}{n}\) for \(n\in\mathbb{N}\). In SEPINN-C (cutoff), we take a truncation level \(N=20\) to approximate the first \(N+1\) coefficients in the series (5.1), and a 3-10-10-10-1 DNN to approximate \(w\). In the first stage, we employ a learning rate 1.0e-3 for the DNN parameters \(\theta\) and 8.0e-3 for coefficients \(\boldsymbol{\gamma}_{N}\), and the prediction error \(e\) after the second stage is 3.93e-2. We present the slices at \(z=\frac{1}{2}\) and \(z=\frac{1}{4}\) in Fig. 4. The true and estimated values of \(\boldsymbol{\gamma}_{n}\) are given in Table 5.2: The first few terms of the expansion are well approximated, but as the index \(n\) gets larger, the approximations of \(\widehat{\gamma}_{n}^{*}\) becomes less accurate. The error initially indeed decreases as the truncation level \(N\) increases, but it does not consistently decrease with \(N\), cf. Table 5.3. The difference between the empirical observation and the theoretical prediction in Theorem 3.3 is possibly due to the optimization error during the training. Next we present the SEPINN-N approximation. We use two 4-layer DNNs, both of 3-10-10-10-1, for \(w\) and \(\Phi\). The PF strategy is run with a maximum 2500 iterations for each fixed \(\boldsymbol{\sigma}\), and the final prediction error \(e\) is 2.05e-2. The maximum error of the SEPINN-N approximation is slightly smaller, and both can give excellent approximations.
Fig. 5 compares the training dynamics for SEPINN-C and SEPINN-N. Fig. 5 (a) shows the convergence for the fist few flux intensity factors \(\boldsymbol{\gamma}\), all initialized to \(1\). Within hundreds of iterations (\(\leq 1000\)), the iterates converge to the exact one steadily (and hence we set the maximum iteration number for \(\boldsymbol{\gamma}\) to \(1000\)). To accurately approximate \(w\), more iterations are needed. Fig. 5 (b) shows the training dynamics for the DNN \(\Phi_{\zeta}\), where \(e_{\boldsymbol{\Phi}}=\|\boldsymbol{\Phi}^{*}-\Phi_{\zeta}\|_{L^{2}(G)}/\| \Phi^{*}\|_{L^{2}(G)}\) (\(G\) is the support of the cut-off function \(\eta\) and \(\Phi^{*}\) is the exact one). The error \(e_{\boldsymbol{\Phi}}\) eventually decreases to \(10^{-2}\). The entire training process of the two methods is similar, cf. Figs. 5 (c) and (d). SEPINN-N takes more iterations than SEPINN-C, but SEPINN-C actually takes longer training time: SEPINN-C requires evaluating the coefficient \(\gamma_{n}\), which incurs taking Laplacian of singular terms, whereas in SEPINN-N, all parameters are trained together.

The next example involves a combination of four singularities.

**Example 5.4**.: _Let the domain \(\Omega=(-\pi,\pi)^{3}\), \(\Gamma_{D}=\{(x,-\pi,z):x\in(-\pi,0),z\in(-\pi,\pi)\}\cup\{(-\pi,y,z):y\in(- \pi,0),z\in(-\pi,\pi)\}\cup\{(x,\pi,z):x\in(0,\pi),z\in(-\pi,\pi)\}\cup\{(\pi,y,z):y\in(0,\pi),z\in(-\pi,\pi)\}\) and \(\Gamma_{D}=\partial\Omega\backslash\Gamma_{N}\). Let the vertices \(\boldsymbol{v}_{1}:(0,-\pi)\), \(\boldsymbol{v}_{2}:(\pi,0)\), \(\boldsymbol{v}_{3}:(0,\pi)\) and \(\boldsymbol{v}_{4}:(-\pi,0)\). This problem has four geometric singularities at boundary edges \(\boldsymbol{v}_{j}\times(-\pi,\pi)\), where the type of the boundary condition changes from Dirichlet to Neumann with interior angles \(\omega_{j}=\pi\), \(j=1,2,3,4\). Define_

\[\Phi_{j}(r_{j},z)=r-\ln(2\cosh r_{j}-2\cos z)=\sum_{n=1}^{\infty}\frac{2}{n} \mathrm{e}^{-nr_{j}}\cos{nz},\quad j=1,2,3,4, \tag{5.2}\]

_and set \(\rho_{j}=1\) and \(R=\frac{1}{2}\) in (3.4), the source \(f=((\frac{4}{\pi^{2}}-\frac{12\pi^{2}}{\pi^{2}})(1-\frac{y^{2}}{\pi^{2}})^{2}+ (\frac{4}{\pi^{2}}-\frac{12\eta^{2}}{\pi^{4}})(1-\frac{x^{2}}{\pi^{2}})^{2}+( 1-\frac{x^{2}}{\pi^{2}})^{2}(1-\frac{y^{2}}{\pi^{2}})^{2}\cos{z}-\sum_{j=1}^{4 }\Delta(\Phi_{j}\eta_{\rho_{j}}s_{j})\) with the singular functions \(s_{j}=r_{j}^{\frac{1}{2}}\cos(\frac{\theta_{j}}{2})\) for \(j=1,3\) and \(s_{j}=r_{j}^{\frac{1}{2}}\sin(\frac{\theta_{j}}{2})\) for \(j=2,4\). The exact solution \(u\) is given by \(u=(1-\frac{x^{2}}{\pi^{2}})^{2}(1-\frac{y^{2}}{\pi^{2}})^{2}\cos{z}+\sum_{j=1} ^{4}\Phi_{j}\eta_{\rho_{j}}s_{j}\)._

This example requires learning a large number of parameters regardless of the method: for SEPINN-C, we have to expand the four singular functions and learn their coefficients, whereas for SEPINN-N, we employ five networks to approximate \(w\) and \(\Phi_{\zeta,j}\) (\(j=1,2,3,4\)). We take the number of sampling points \(N_{d}=800\) and \(N_{n}=1200\) on the boundaries \(\Gamma_{D}\) and \(\Gamma_{N}\), respectively, and \(N_{r}=10000\) in the domain \(\Omega\).

First we present the SEPINN-C approximation. From (5.2), we have the explicit coefficients \(\gamma_{0}^{*}=0\) and \(\gamma_{n}^{*}=\frac{2}{n}\) for \(n\in\mathbb{N}\). In SEPINN-C, we take a truncation level \(N=15\) for all four singularities. In the first stage, we take a learning rate 2.0e-3 for \(\theta\) and for coefficients \(r_{1}=r_{3}=1.1\)e-2, \(r_{2}=r_{4}=7.0\)e-3. The final prediction error \(e\) is 2.08e-2, and the estimated coefficients \(\boldsymbol{\widehat{\gamma}}^{*}\) are shown in Table 5.4. The first few coefficients are well approximated, but the high-order ones are less accurately approximated.

\begin{table}
\begin{tabular}{c|c c c c} \hline \hline \(N\) & 5 & 10 & 15 & 20 \\ \hline \(e\) & 8.15e-2 & 5.96e-2 & 3.48e-2 & 3.93e-2 \\ \hline rate & & 0.45 & 1.32 & -0.42 \\ \hline \hline \end{tabular}
\end{table}
Table 5.3: Convergence of the SEPINN-C approximation with respect to the truncation level \(N\) for Example 5.3.

Figure 5: The training dynamics of SEPINN-C and SEPINN-N: (a) the variation of first few coefficients versus iteration index \(i\), (b) the error \(e_{\boldsymbol{\Phi}}\) versus iteration index \(i\), (c) the decay of the loss \(\widehat{\mathcal{L}}\) versus iteration index \(i\), (d) the error \(e\) versus iteration index \(i\).
Next we show the SEPINN-N approximation, obtained with five 4-layer 3-10-10-10-1 DNNs to approximate \(w\) and \(\Phi_{\zeta,j}\) (\(j=1,2,3,4\)) separately. The training process suffers from the following problem: using only L-BFGS tends to be trapped into a local minimum of the loss \(\widehat{\mathcal{L}}_{\boldsymbol{\sigma}}\), which persists even after adjusting extensively the hyper-parameters. Therefore, we take the following approach: we first train the DNNs with Adam (learning rate \(r=4.0\)e-3) for 1000 iterations and then switch to L-BFGS (learning rate \(r=0.2\)) for a maximum 4000 iterations. This training strategy greatly improves the accuracy. The prediction error \(e\) after the second stage is 3.83e-2. The approximation is fairly accurate but slightly less accurate than that by SEPINN-C in Fig. 6 in both \(L^{\infty}(\Omega)\) and \(L^{2}(\Omega)\) norms.

The next example illustrates SEPINN on the Helmholtz equation.

**Example 5.5**.: _Consider the following problem in both 2D and 3D domains:_

\[-\Delta u+\pi^{2}u=f,\quad\text{in}\quad\Omega\]

_with the source \(f\) taken to be \(f=(4x^{3}+6x)y\mathrm{e}^{x^{2}-1}(\mathrm{e}^{y^{2}-1}-1)-(4y^{3}+6y)x \mathrm{e}^{y^{2}-1}(\mathrm{e}^{x^{2}-1}-1)+\pi^{2}xy(\mathrm{e}^{x^{2}-1}-1) (\mathrm{e}^{y^{2}-1}-1)-\Delta(\mathrm{e}^{-\pi r}\eta_{\rho}s)+\pi^{2} \mathrm{e}^{-\pi r}\eta_{\rho}s\) for the 2D domain \(\Omega=\Omega_{1}=(-1,1)^{2}\backslash([0,1)\times(-1,0])\), and \(f=-((4x^{3}+6x)y\mathrm{e}^{x^{2}-1}(\mathrm{e}^{y^{2}-1}-1)-(4y^{3}+6y)x \mathrm{e}^{y^{2}-1}(\mathrm{e}^{x^{2}-1}-1)+2\pi^{2}xy(\mathrm{e}^{x^{2}-1}-1 )(\mathrm{e}^{y^{2}-1}-1))\sin(\pi z)-\Delta(\mathrm{e}^{-\sqrt{2}\pi r}\eta_{ \rho}s\sin(\pi z))+\pi^{2}\mathrm{e}^{-\sqrt{2}\pi r}\eta_{\rho}s\sin(\pi z)\) for the 3D domain \(\Omega=\Omega_{2}=\Omega_{1}\times(-1,1)\), with \(\rho=1\) and

\begin{table}
\begin{tabular}{c|c|c c c c} \hline \(n\backslash\gamma_{j,n}\) & exact & \(\gamma_{1,n}\) & \(\gamma_{2,n}\) & \(\gamma_{3,n}\) & \(\gamma_{4,n}\) \\ \hline
0 & 0.0000e0 & 2.4763e-4 & 2.9810e-4 & 1.1018e-3 & -1.0403e-4 \\ \hline
1 & 2.0000e0 & 1.9985e0 & 1.9932e0 & 1.9902e0 & 1.9968e0 \\ \hline
2 & 1.0000e0 & 1.0087e0 & 1.0020e0 & 9.9630e-1 & 1.0005e0 \\ \hline
3 & 6.6667e-1 & 6.7281e-1 & 6.6629e-1 & 6.5206e-1 & 6.5664e-1 \\ \hline
4 & 5.0000e-1 & 5.1538e-1 & 4.8951e-1 & 4.8902e-1 & 4.8838e-1 \\ \hline
5 & 4.0000e-1 & 4.0184e-1 & 3.9272e-1 & 3.9834e-1 & 3.7305e-1 \\ \hline
6 & 3.3333e-1 & 3.1959e-1 & 3.1089e-1 & 3.2622e-1 & 2.9855e-1 \\ \hline
7 & 2.8571e-1 & 2.5984e-1 & 2.7772e-1 & 2.3514e-1 & 2.4714e-1 \\ \hline
8 & 2.5000e-1 & 1.8928e-1 & 2.3612e-1 & 2.1324e-1 & 1.9534e-1 \\ \hline
9 & 2.2222e-1 & 1.3174e-1 & 1.6849e-1 & 3.0520e-1 & 1.5977e-1 \\ \hline
10 & 2.0000e-1 & 4.9766e-2 & 1.3174e-1 & 2.4587e-1 & 1.2393e-1 \\ \hline
11 & 1.8182e-1 & 1.7544e-3 & 1.4534e-1 & 2.0706e-1 & 1.1372e-1 \\ \hline
12 & 1.6667e-1 & -6.0320e-2 & 1.3719e-1 & 2.9289e-1 & 9.5344e-2 \\ \hline
13 & 1.5385e-1 & -8.6117e-2 & 1.0655e-1 & 2.1058e-1 & 5.5091e-2 \\ \hline
14 & 1.4286e-1 & -5.5869e-2 & 1.1616e-1 & 1.2455e-1 & 5.2139e-2 \\ \hline
15 & 1.3333e-1 & -6.2956e-2 & 6.1959e-2 & 2.2816e-1 & 1.4306e-2 \\ \hline \end{tabular}
\end{table}
Table 5.4: The comparison between true and estimated values of the parameters \(\gamma_{j,n}\) for Example 5.4, with 5 significant digits.

Figure 6: The SEPINN-C and SEPINN-N approximations for Example 5.4, slices at \(z=\frac{\pi}{2}\) (top) and \(z=\pi\) (bottom).
\(R=\frac{1}{2}\) in (3.4), and \(s=r^{\frac{2}{3}}\sin(\frac{2\theta}{3})\), \(\Gamma_{D}=\partial\Omega\). The exact solution \(u\) is given by_

\[u=\begin{cases}xy(\mathrm{e}^{x^{2}-1}-1)(\mathrm{e}^{\mathsf{v}^{2}-1}-1)+ \mathrm{e}^{-\pi r}\eta_{\rho}s,\quad\Omega=\Omega_{1},\\ xy(\mathrm{e}^{x^{2}-1}-1)(\mathrm{e}^{\mathsf{v}^{2}-1}-1)\sin(\pi z)+\mathrm{e }^{-\sqrt{2}\pi r}\eta_{\rho}s\sin(\pi z),\quad\Omega=\Omega_{2}.\end{cases}\]

For the 2D problem, we employ a 2-20-20-20-1 DNN. The first stage of the PF strategy yields an estimate \(\widehat{\gamma}^{*}=0.9939\), and the final prediction error \(e\) is 1.02e-2. This accuracy is comparable with that for the Poisson problem in Example 5.1. In the 3D case, we have the coefficients \(\gamma_{1}^{*}=1\) and \(\gamma_{n}^{*}=0\) for \(n\in\mathbb{N}\backslash\{1\}\). In SEPINN-C, we take \(N=10\) terms. In the first stage, we set the learning rate 2.0e-3 for the DNN parameters \(\theta\) and 4.0e-3 for coefficients \(\boldsymbol{\gamma}\). The estimated \(\widehat{\boldsymbol{\gamma}}^{*}\) are shown in Table 5.5, which approximate accurately the exact coefficients, and the prediction error \(e\) after the second stage is 3.30e-2. In SEPINN-N, we use five 4-layer (3-20-20-20-1) DNNs to approximate \(w\) and 3-10-10-10 approximate \(\Phi\), and optimize the loss with L-BFGS with a maximum of 2500 iterations. The final prediction error \(e\) is 2.55e-2. The results are shown in Fig. 7. SEPINN-C and SEPINN-N give very similar pointwise errors, and the singularity at the reentrant corner is accurately resolved, due to singularity enrichment. This shows clearly the flexibility of SEPINN for the screened Poisson equation with geometric singularities.

Last, we illustrate SEPINN on the Laplacian eigenvalue problem in an L-shaped domain.

**Example 5.6**.: _On the domain \(\Omega=(-1,1)^{2}\backslash([0,1)\times(-1,0])\), consider the eigenvalue problem:_

\[\begin{cases}-\Delta u=\mu u,&\text{in }\Omega,\\ u=0,&\text{on }\partial\Omega,\end{cases}\]

_where \(\mu>0\) is the eigenvalue and \(u\not\equiv 0\) is the corresponding eigenfunction._

The eigenvalue problem on an L-shaped domain has been studied extensively [24, 44, 61]. The eigenfunctions may have singularity around the reentrant corner, but the analytic forms appear unavailable: the first eigenfunction \(u_{1}\) has a leading singular term \(r^{\frac{2}{3}}\sin(\frac{2}{3}\theta)\), the second one \(u_{2}\) has \(r^{\frac{1}{3}}\sin(\frac{4}{3}\theta)\)[24], and the third \(u_{3}\) is analytic, given by \(u_{3}(x_{1},x_{2})=\sin(\pi x_{1})\sin(\pi x_{2})\). To illustrate the flexibility of SEPINN, we compute the first two leading eigenpairs \((\mu_{1},u_{1})\) and \((\mu_{2},u_{2})\). The preceding discussions

\begin{table}
\begin{tabular}{c|c|c|c|c|c|c} \hline  & \(\gamma_{0}\) & \(\gamma_{1}\) & \(\gamma_{2}\) & \(\gamma_{3}\) & \(\gamma_{4}\) & \(\gamma_{5}\) \\ \hline estimate & 5.6180e-4 & 9.8857e-1 & 5.8005e-3 & -2.3331e-3 & 3.3624e-3 & -2.3693e-3 \\ \hline  & \(\gamma_{6}\) & \(\gamma_{7}\) & \(\gamma_{8}\) & \(\gamma_{9}\) & \(\gamma_{10}\) & \\ \hline estimate & 2.7073e-3 & -2.9761e-3 & -1.8276e-3 & -9.7370e-4 & 1.1867e-3 & \\ \hline \end{tabular}
\end{table}
Table 5.5: The estimated coefficients \(\gamma_{n}\) in the 3D case for Example 5.5, with five significant digits.

Figure 7: 2D problem with SEPINN (top) and 3D problem with SEPINN-C and SEPINN-N (bottom) for Example 5.5, with slice at \(z=\frac{1}{2}\).
indicate \(u_{1}\in H^{1}(\Omega)\) and \(u_{2}\in H^{2}(\Omega)\), and we split the leading singularity from both functions in order to benefit from SEPINN, i.e., \(u_{i}=w_{i}+\gamma_{i}\eta_{\rho}s\), with \(s=r^{\frac{2}{3}}\sin(\frac{2\theta}{3})\) and \(w_{i}\) is approximated by a DNN. Following the ideas in [8] and SEPINN, we employ the following loss

\[\mathcal{L}_{\mathbf{\sigma}}(w_{1},w_{2};\gamma_{1},\gamma_{2})= \sum_{i=1}^{2}\left(\|\Delta(w_{i}+\gamma_{i}\eta_{\rho}s)+R(u_{i} )(w_{i}+\gamma_{i}\eta_{\rho}s)\|_{L^{2}(\Omega)}^{2}+\sigma_{1}\|w_{i}\|_{L^{ 2}(\partial\Omega)}^{2}\right. \tag{5.3}\] \[+\alpha\big{|}\|w_{i}+\gamma_{i}\eta_{\rho}s\|_{L^{2}(\Omega)}^{2 }-1\big{|}+\nu_{i}R(u_{i})\Big{)}+\beta\left|(w_{1}+\gamma_{1}\eta_{\rho}s,w_{ 2}+\gamma_{2}\eta_{\rho}s)_{L^{2}(\Omega)}\right|.\]

where \(\alpha\), \(\nu_{i}(i=1,2)\) and \(\beta\) are hyper-parameters and \(R(u_{i})\) is the Rayleigh quotient:

\[R(u_{i})=\|\nabla u_{i}\|_{L^{2}(\Omega)}^{2}/\|u_{i}\|_{L^{2}( \Omega)}^{2},\quad i=1,2, \tag{5.4}\]

which estimates the eigenvalue \(\mu_{i}\) using the eigenfunction \(u_{i}\), by Rayleigh's principle. We employ an alternating iteration method: we first approximate the eigenfunction \(u_{i}\) by minimizing the loss (5.3) and then update the eigenvalue \(\mu_{i}\) by (5.4), which is then substituted back into (5.3). These two steps are repeated until convergence.

In SEPINN, we employ two 2-10-10-10-10-10-10-1 DNNs to approximate the regular parts \(w_{1}\) and \(w_{2}\), and take \(\alpha=100\), \(\beta=135\), \(\nu_{1}=0.02\) and \(\nu_{2}=0.01\), and determine the parameter \(\sigma_{d}\) by the PF strategy. We use Adam with a learning rate 2e-3 for all DNN parameters and \(\mathbf{\gamma}\). Table 5.6 shows that singularity enrichment helps solve the eigenvalue problem. Indeed, we can approximate \(u_{1}\) better and get more accurate eigenvalue estimates. Note that during the training process of the standard PINN, the DNN approximation actually directly approaches \(u_{2}\) and cannot capture \(u_{1}\), since it cannot be resolved accurately. Even with a larger penalty \(\nu_{1}\) for the first eigenvalue, this does not improve the accuracy, while \(u_{2}\) can be approximated well since \(u_{2}\in H^{2}(\Omega)\). This is also seen from estimated stress intensity factor in Table 5.6, which is much larger for \(u_{1}\) than for \(u_{2}\), indicating stronger singularity of \(u_{1}\).

## 6 Conclusions

In this work, we have developed a family of neural solvers for solving boundary value problems with geometric singularities, e.g., corner singularity and mixed boundary conditions in the 2D case, and edge singularities in the 3D case. The basic idea is to enrich the ansatz space, which deep neural networks span, by incorporating specially designed singular functions. These singular functions are designed to capture the leading singularities of the exact solution. In so doing, we can achieve a significantly improved convergence rate. Additionally, we provide preliminary theoretical guarantees of the approach. We discuss several variants of the method, depending on the specific scenarios, and discuss the extensions to the Helmholtz equation and eigenvalue problems. To the best of our knowledge, it is the first work systematically exploring the use of singularity enrichment in a neural PDE solver. The extensive numerical experiments indicate that the approach is indeed highly effective and flexible, and works for a broad range of problem settings.
