Enriched Holomorphic Neural Networks for 2D and 3D Laplace Problems

# Enriched Holomorphic Neural Networksfor 2D and 3D Laplace Problems
Source: [https://arxiv.org/html/2610.12541](https://arxiv.org/html/2610.12541)
Andrei Vlad NiculaAffiliation:University of AmsterdamEmail:[jan\.jelinek@student\.uva\.nl](mailto:)Mihnea Andrei BârsanAffiliation:University of AmsterdamEmail:[avama@dtu\.dk](mailto:)Robert\-S,tefan SofroniAffiliation:University of AmsterdamAgneya Bhardwaj Jan Jelínek Alejandro Valverde MahouAffiliation:University of AmsterdamAffiliation:Technical University of Denmark\{kalp\.shah, andrei\.nicula, mihnea\.barsan, robert\.sofroni

###### Abstract

Boundary\-value problems for the Laplace equation develop singular behaviour near corners, edges and junctions that smooth approximators resolve poorly\. We combine holomorphic neural networks with singular functions derived from the local geometry and boundary conditions\. The representation is harmonic by construction, so no interior residual is needed: the network fits the smooth remainder from boundary data alone, and the enrichment terms capture the singular behaviour\. Singular exponents are trainable, so the construction also applies when the local expansion is not known in advance\. We extend the approach to three dimensions using independent networks evaluated along complex null directions\. The enriched model is the most accurate tested method on every problem\. Its margin is several orders of magnitude where the enrichment captures the singularity exactly, and a factor of about 1\.7–2\.1 over the next\-best method on a CdZnTe radiation detector with many singular sites, in both two and three dimensions\.

## 1Introduction

SolvingPartial Differential Equations \(PDEs\)is difficult on domains with non\-trivial geometries or boundary conditions\.Physics\-Informed neural networks \(PINNs\)\([Raissi et al\., 2019](https://arxiv.org/html/2610.12541#bib.bib15)\)offer an alternative by incorporating the governing equations directly into the training objective\. The solution is represented by a neural network, so it is continuous and its derivatives come from automatic differentiation\. The method is therefore mesh\-free and differentiable end\-to\-end, which motivates its use for irregular geometries and inverse problems\([Luo et al\., 2025](https://arxiv.org/html/2610.12541#bib.bib16)\)\.

Unfortunately, the resulting objective suffers from competing loss terms, ill\-conditioning and imbalanced gradients, which slow convergence or degrade accuracy\([Wang et al\., 2021](https://arxiv.org/html/2610.12541#bib.bib17);[Wang et al\., 2022](https://arxiv.org/html/2610.12541#bib.bib18);[Rahaman et al\., 2019](https://arxiv.org/html/2610.12541#bib.bib19)\)\. Much work has therefore improved the optimization and representation of PINNs\. Adaptive loss weighting balances the terms of the objective\([Kendall et al\., 2018](https://arxiv.org/html/2610.12541#bib.bib22);[Heydari et al\., 2019](https://arxiv.org/html/2610.12541#bib.bib23);[Wang et al\., 2022](https://arxiv.org/html/2610.12541#bib.bib18);[McClenny and Braga\-Neto, 2023](https://arxiv.org/html/2610.12541#bib.bib20)\), and adaptive collocation changes where the residual is evaluated\([Wu et al\., 2023](https://arxiv.org/html/2610.12541#bib.bib21)\)\. Other approaches change the representation, through periodic activations or Fourier features\([Sitzmann et al\., 2020](https://arxiv.org/html/2610.12541#bib.bib24);[Liu et al\., 2025](https://arxiv.org/html/2610.12541#bib.bib25)\), or by replacing the multilayer perceptron with a Kolmogorov–Arnold network\([Liu et al\., 2024](https://arxiv.org/html/2610.12541#bib.bib26)\)\. Domain\-decomposition approaches such as conservative PINNs \(cPINNs\)\([Jagtap et al\., 2020](https://arxiv.org/html/2610.12541#bib.bib28)\)and extended PINNs \(XPINNs\)\([Jagtap and Karniadakis, 2020](https://arxiv.org/html/2610.12541#bib.bib29)\)partition the domain into subdomains coupled through interface conditions\.

We address the representation side of these difficulties\. We focus on the Laplace equation,Δ​u=0\\Delta u=0, and build on thePhysics\-Informed Holomorphic Neural Networks \(PIHNNs\)formulations introduced by[Calafà et al\. \(2024\)](https://arxiv.org/html/2610.12541#bib.bib11)for two\-dimensional problems and further developed by[Ballini et al\. \(2026\)](https://arxiv.org/html/2610.12541#bib.bib13)in three dimensions\. These representations have also been shown to benefit from enrichment with known singularities for specific problems, such as crack tips in plane elasticity\([Calafà et al\., 2025](https://arxiv.org/html/2610.12541#bib.bib12)\)\. Our approach addresses the enrichment for a set of local boundary configurations rather than for one particular problem: we construct a dictionary mapping each local boundary configuration to its enrichment family, together with a rule for orienting the branch cuts so that the representation remains single\-valued inside the domain\. We develop the construction in two and three dimensions\.

## 2Related work

Representations that satisfy the governing equation in the interior reduce a boundary value problem to determining model parameters from boundary information\. Classical examples include the boundary element method, the method of fundamental solutions\([Fairweather and Karageorghis, 1998](https://arxiv.org/html/2610.12541#bib.bib8)\), harmonic kernel collocation\([Schaback, 2008](https://arxiv.org/html/2610.12541#bib.bib9)\), and rational approximations such as the lightning Laplace solver\([Gopal and Trefethen, 2019](https://arxiv.org/html/2610.12541#bib.bib32)\)\. These methods construct the solution from fixed analytical families, such as harmonic polynomials or fundamental solutions centred outside the domain\([Axler et al\., 2001](https://arxiv.org/html/2610.12541#bib.bib1)\)\. Neural representations instead learn part of the approximating family, at the cost of a nonconvex fit\.[Ghosh et al\. \(2023\)](https://arxiv.org/html/2610.12541#bib.bib10)represent harmonic functions as the real part of a complex\-valued neural network, while[Calafà et al\. \(2024\)](https://arxiv.org/html/2610.12541#bib.bib11)developPIHNNsfor two\-dimensional problems\. In 3D, where the correspondence with complex analysis is lost,[Ballini et al\. \(2026\)](https://arxiv.org/html/2610.12541#bib.bib13)write the solution as an integral over projection directions, approximate it by a weighted sum over a finite set of directions, and share one network across them\.

Exact satisfaction of the governing equation, however, does not guarantee an efficient representation of solutions with reduced regularity\. At corners and boundary\-condition junctions, local singularities arise whose exponents depend on the geometry and boundary conditions\([Grisvard, 1985](https://arxiv.org/html/2610.12541#bib.bib3)\)\. Enriching the approximation space with known singular functions is a classical remedy in numericalPDEmethods\([Fix et al\., 1973](https://arxiv.org/html/2610.12541#bib.bib4);[Li and Lu, 2000](https://arxiv.org/html/2610.12541#bib.bib5)\)\. Among neural solvers, the singularity\-enrichedPINN\(SEPINN\) of[Hu et al\. \(2024\)](https://arxiv.org/html/2610.12541#bib.bib14)augments a real\-valuedPINNwith singular functions for polygonal corners, mixed boundary conditions and edges in 3D, while[Calafà et al\. \(2025\)](https://arxiv.org/html/2610.12541#bib.bib12)enrichPIHNNswith the square\-root term from Williams’ expansion at crack tips\. In these approaches, the singularity and its exponent are prescribed for the problem at hand\.

## 3Method

### 3\.1Problem formulation

We consider the Laplace boundary\-value problem on a bounded domainΩ⊂ℝd\\Omega\\subset\\mathbb\{R\}^\{d\}, withd∈\{2,3\}d\\in\\\{2,3\\\}\. The boundary is partitioned into disjoint Dirichlet and Neumann portions,∂Ω=ΓD∪ΓN\\partial\\Omega=\\Gamma\_\{D\}\\cup\\Gamma\_\{N\}, and the solutionu:Ω→ℝu:\\Omega\\rightarrow\\mathbb\{R\}satisfies

Δ​u​\(𝒙\)\\displaystyle\\Delta u\(\\bm\{x\}\)=0,\\displaystyle=0,𝒙∈Ω,\\displaystyle\\bm\{x\}\\in\\Omega,\(1\)u⁡\(𝒙\)\\displaystyle u\(\\bm\{x\}\)=gD​\(𝒙\),\\displaystyle=g\_\{D\}\(\\bm\{x\}\),𝒙∈ΓD,\\displaystyle\\bm\{x\}\\in\\Gamma\_\{D\},∇u​\(𝒙\)⋅𝒏​\(𝒙\)\\displaystyle\\nabla u\(\\bm\{x\}\)\\cdot\\bm\{n\}\(\\bm\{x\}\)=gN​\(𝒙\),\\displaystyle=g\_\{N\}\(\\bm\{x\}\),𝒙∈ΓN,\\displaystyle\\bm\{x\}\\in\\Gamma\_\{N\},wheregDg\_\{D\}andgNg\_\{N\}denote the prescribed Dirichlet and Neumann data, respectively, and𝒏\\bm\{n\}is the outward unit normal\. We assume thatΓD\\Gamma\_\{D\}has positive measure, so that the solution is unique\.

We represent the boundary as a collection of segments in two dimensions or faces in three dimensions\. We refer to points at which adjacent boundary pieces meet, or at which the boundary condition or its prescribed value changes, as*junctions*\. These local boundary configurations determine the singular terms used to enrich the neural representation\.

### 3\.2Holomorphic functions

A functionf:U→ℂf:U\\to\\mathbb\{C\}on an open setU⊂ℂU\\subset\\mathbb\{C\}is holomorphic if it is complex differentiable at every point ofUU\. Writingz=x\+i​yz=x\+iyandf⁡\(z\)=p⁡\(x,y\)\+i​q​\(x,y\)f\(z\)=p\(x,y\)\+iq\(x,y\), and using the subscript notation for partial derivatives throughout so thatpx=∂p∂xp\_\{x\}=\\frac\{\\partial p\}\{\\partial x\}, for continuously differentiableppandqq, complex differentiability is equivalent to the Cauchy–Riemann equationspx=qyp\_\{x\}=q\_\{y\}andpy=−qxp\_\{y\}=\-q\_\{x\}\([Axler et al\., 2001](https://arxiv.org/html/2610.12541#bib.bib1);[Evans, 2010](https://arxiv.org/html/2610.12541#bib.bib2)\)\. Differentiating these once more givesΔ​p=px​x\+py​y=qy​x−qx​y=0\\Delta p=p\_\{xx\}\+p\_\{yy\}=q\_\{yx\}\-q\_\{xy\}=0, and the same forqq, so both parts of a holomorphic function are harmonic\. Conversely, on a simply connected domain every harmonic function is the real part of a holomorphic one\. A function that is holomorphic on all ofℂ\\mathbb\{C\}is called*entire*\.

Two properties of holomorphic functions are central to this work\. First, they are closed under composition and linear combination: a network assembled from holomorphic layers is therefore holomorphic, and adding further holomorphic terms to its output keeps it so\. Second, the real part of any holomorphic function is harmonic so a holomorphic network represents a harmonic function for every value of its weights\. Together these mean the network is harmonic by construction, with no interior residual required\. The Cauchy–Riemann equations also give the complex derivative asf′​\(z\)=px−i​pyf^\{\\prime\}\(z\)=p\_\{x\}\-ip\_\{y\}, from which the normal derivative follows directly\.

### 3\.3Enrichment functions

Near corners and junctions, solutions to the Laplace equation develop singular behaviour that smooth approximators resolve poorly: the solution may be bounded, but its gradient can blow up\. Rather than relying on the network to learn this from boundary samples, we augment the representation with enrichment functions drawn from local singular expansions for the given geometry and boundary conditions\([Grisvard, 1985](https://arxiv.org/html/2610.12541#bib.bib3);[Hu et al\., 2024](https://arxiv.org/html/2610.12541#bib.bib14);[Calafà et al\., 2025](https://arxiv.org/html/2610.12541#bib.bib12)\)\.

Consider a planar wedge with interior angleω\\omegaand polar coordinates\(r,θ\)\(r,\\theta\)centred at its vertex\. For a homogeneous Dirichlet condition atθ=0\\theta=0and a homogeneous Neumann condition atθ=ω\\theta=\\omega, the leading singular mode isrλ​sin⁡\(λ​θ\)r^\{\\lambda\}\\sin\(\\lambda\\theta\)withλ=π/\(2​ω\)\\lambda=\\pi/\(2\\omega\)\. The solution stays bounded asr→0r\\to 0, but its gradient grows asrλ−1r^\{\\lambda\-1\}and diverges wheneverλ<1\\lambda<1, i\.e\. whenω\>π/2\\omega\>\\pi/2\. This is illustrated in[Figure1](https://arxiv.org/html/2610.12541#S3.F1)forω=2​π/3\\omega=2\\pi/3,λ=3/4\\lambda=3/4, where the field looks smooth but the gradient concentrates at the vertex\.

In the local complex coordinatew=r​ei​θw=re^\{i\\theta\}, the singular modes are the real and imaginary parts ofwλ=exp⁡\(λ​Log⁡w\)w^\{\\lambda\}=\\exp\(\\lambda\\operatorname\{Log\}w\)\. We place the branch cut ofLog\\operatorname\{Log\}along a fixed exterior direction of the junction \([AppendixA](https://arxiv.org/html/2610.12541#A1)\)\. The complex coefficient of each term then absorbs the phase of the wall from whichθ\\thetais measured and the order of the two conditions\. At each junction there is not one singular exponent but a countably infinite sequenceλ1<λ2<⋯\\lambda\_\{1\}<\\lambda\_\{2\}<\\cdots, each indexing a valid local mode\. For a D→\\toN \(Dirichlet to Neumann\) junction at angleω\\omega, for instance,λk=\(2​k−1\)​π/\(2​ω\)\\lambda\_\{k\}=\(2k\-1\)\\pi/\(2\\omega\)\. The leading term \(k=1k=1\) is the most singular and dominates near the vertex, and higher terms are progressively smoother\. The number of terms per junction is a hyperparameter we call the*ladder depth*\. We use ladder depth one except on Motz’s problem and the cracked beam, where the first four exponents are included \([AppendixF](https://arxiv.org/html/2610.12541#A6)\)\.[Table1](https://arxiv.org/html/2610.12541#S3.T1)lists the enrichment terms\.

![Refer to caption](https://arxiv.org/html/2610.12541v1/enrichment_plot.png)Figure 1:Singular behaviour at a D→\\toN junction withω=2​π/3\\omega=2\\pi/3,λ=3/4\\lambda=3/4\. Left: the leading moderλ​sin⁡\(λ​θ\)r^\{\\lambda\}\\sin\(\\lambda\\theta\)is bounded at the vertex\. Right: the gradient magnitude\|∇u\|∼rλ−1=r−1/4\|\\nabla u\|\\sim r^\{\\lambda\-1\}=r^\{\-1/4\}diverges there\. The enrichment termwλw^\{\\lambda\}represents this mode directly\.Table 1:Enrichment dictionary\. Each junction matches one row and contributes the termEλE\_\{\\lambda\}with a trainable coefficientcqc\_\{q\}\. Local coordinates are centred at the junctionzcz\_\{c\}, withw=z−zc=r​ei​θw=z\-z\_\{c\}=re^\{i\\theta\};ω\\omegais the interior angle andk=1,2,…k=1,2,\\dotsindexes the ladder\. Every term is holomorphic, so the representation remains exactly harmonic\.[Figure5](https://arxiv.org/html/2610.12541#A1.F5)is a visualization of this table\.JunctionLocal solution\|∇u\|\|\\nabla u\|TermEλE\_\{\\lambda\}Dirichlet to Neumannrλ​sin⁡\(λ​θ\)r^\{\\lambda\}\\sin\(\\lambda\\theta\),λ=\(2​k−1\)​π/\(2​ω\)\\lambda=\(2k\-1\)\\pi/\(2\\omega\)rλ−1r^\{\\lambda\-1\}wλw^\{\\lambda\}Both Dirichlet, value jumpΔ​g\\Delta g\(Δ​g/ω\)​θ\(\\Delta g/\\omega\)\\,\\thetar−1r^\{\-1\}arg⁡w\\arg wBoth Neumann, flux jump,ω∈π​ℕ\\omega\\in\\pi\\mathbb\{N\}ℜ⁡\[w​Log⁡w\]\\Re\[w\\operatorname\{Log\}w\]\|ln⁡r\|\|\\ln r\|w​Log⁡ww\\operatorname\{Log\}wBoth Neumann, equal flux,ω\>π\\omega\>\\pirλ​cos⁡\(λ​θ\)r^\{\\lambda\}\\cos\(\\lambda\\theta\),λ=k​π/ω\\lambda=k\\pi/\\omegarλ−1r^\{\\lambda\-1\}wλw^\{\\lambda\}Both Dirichlet, equal values,ω\>π\\omega\>\\pirλ​sin⁡\(λ​θ\)r^\{\\lambda\}\\sin\(\\lambda\\theta\),λ=k​π/ω\\lambda=k\\pi/\\omegarλ−1r^\{\\lambda\-1\}wλw^\{\\lambda\}OtherwiseC∞C^\{\\infty\}away from vertexboundednone#### Trainable exponents\.

The exponents in[Table1](https://arxiv.org/html/2610.12541#S3.T1)and[Figure5](https://arxiv.org/html/2610.12541#A1.F5)can be fixed at their theoretical values or trained\. Sincewλw^\{\\lambda\}is holomorphic for everyλ\>0\\lambda\>0, varying the exponent does not affect harmonicity, but only the exponents of the table satisfy the local boundary conditions, so a trainedλ\\lambdais a fitted parameter rather than the singular exponent itself\. Each term carries one real parameterτ\\tau, withλ=ε\+softplus⁡\(τ\)\\lambda=\\varepsilon\+\\operatorname\{softplus\}\(\\tau\)initialised at the predicted value, except the conical degree, which starts at1/21/2\. This allows the construction to be applied when the exponent is not known in advance, whereas prior enrichment schemes fixλ\\lambda\([Hu et al\., 2024](https://arxiv.org/html/2610.12541#bib.bib14);[Calafà et al\., 2025](https://arxiv.org/html/2610.12541#bib.bib12)\)\.

#### Extension to 3D\.

For a straight edge with unit tangentttand transverse orthonormal framee1,e2e\_\{1\},e\_\{2\}, we set

w=\(𝒙−𝒙e\)⋅\(e1\+i​e2\),Eλ​\(𝒙\)=wλ\.w=\(\\bm\{x\}\-\\bm\{x\}\_\{e\}\)\\cdot\(e\_\{1\}\+ie\_\{2\}\),\\qquad E\_\{\\lambda\}\(\\bm\{x\}\)=w^\{\\lambda\}\.\(2\)The directione1\+i​e2e\_\{1\}\+ie\_\{2\}satisfiesa⋅a=0a\\cdot a=0, soEλE\_\{\\lambda\}is harmonic wherever its branch is holomorphic\.

### 3\.4Enriched Holomorphic neural networks

PIHNNsexploit the compositional property of holomorphic functions to construct a holomorphic neural representation\. Letfη:ℂ→ℂf\_\{\\eta\}:\\mathbb\{C\}\\to\\mathbb\{C\}denote such a network, with complex\-valued parametersη\\eta\. Since holomorphic functions are closed under composition,fηf\_\{\\eta\}is holomorphic by construction, and its real part is therefore harmonic\([Ghosh et al\., 2023](https://arxiv.org/html/2610.12541#bib.bib10);[Calafà et al\., 2024](https://arxiv.org/html/2610.12541#bib.bib11)\)\.

ThePIHNNrepresents the remaining part of the solution, andQQenrichment termsEλqE\_\{\\lambda\_\{q\}\}, one per singular site and ladder index, represent the behaviour near the sites\. The output of the model is

uη,𝒄​\(x,y\)=ℜ⁡\[fη​\(z\)\+∑q=1Qcq​Eλq​\(z\)\],u\_\{\\eta,\\bm\{c\}\}\(x,y\)=\\Re\\left\[f\_\{\\eta\}\(z\)\+\\sum\_\{q=1\}^\{Q\}c\_\{q\}E\_\{\\lambda\_\{q\}\}\(z\)\\right\],\(3\)wherecqc\_\{q\}determine the amplitude and phase of each termEλqE\_\{\\lambda\_\{q\}\}\.

Since the representation is harmonic by construction, no interior residual is needed, and the parameters are fitted to the boundary conditions by minimizing

ℒ⁡\(η,𝒄\)=1\|XD\|​∑𝒙∈XD‖uη,𝒄​\(𝒙\)−gD​\(𝒙\)‖2\+1\|XN\|​∑𝒙∈XN‖∂nuη,𝒄​\(𝒙\)−gN​\(𝒙\)‖2,\\mathcal\{L\}\(\\eta,\\bm\{c\}\)=\\frac\{1\}\{\|X\_\{D\}\|\}\\sum\_\{\\bm\{x\}\\in X\_\{D\}\}\\left\\\|u\_\{\\eta,\\bm\{c\}\}\(\\bm\{x\}\)\-g\_\{D\}\(\\bm\{x\}\)\\right\\\|^\{2\}\+\\frac\{1\}\{\|X\_\{N\}\|\}\\sum\_\{\\bm\{x\}\\in X\_\{N\}\}\\left\\\|\\partial\_\{n\}u\_\{\\eta,\\bm\{c\}\}\(\\bm\{x\}\)\-g\_\{N\}\(\\bm\{x\}\)\\right\\\|^\{2\},\(4\)whereXDX\_\{D\}andXNX\_\{N\}denote the corresponding boundary sample sets, and a term is omitted when the corresponding boundary condition is absent\.

#### Extension to 3D

In 3D there is no single complex coordinate, so we use complex projections of the forma⋅𝒙a\\cdot\\bm\{x\}, wherea∈ℂ3a\\in\\mathbb\{C\}^\{3\}\. For a holomorphic functionffand a fixed directionaa, the chain rule gives

Δ𝒙​f​\(a⋅𝒙\)=\(a⋅a\)​f′′​\(a⋅𝒙\),a⋅a=∑k=13ak2\.\\Delta\_\{\\bm\{x\}\}f\(a\\cdot\\bm\{x\}\)=\(a\\cdot a\)f^\{\\prime\\prime\}\(a\\cdot\\bm\{x\}\),~~~~~a\\cdot a=\\sum\_\{k=1\}^\{3\}a\_\{k\}^\{2\}\.\(5\)Thus, choosing a null direction satisfyinga⋅a=0a\\cdot a=0makes the projected function harmonic\.

We useKKsuch projections and sum their outputs:

u𝜼,𝒄​\(𝒙\)=ℜ⁡\[∑k=1Kfηk​\(ζk\)\+∑q=1Qcq​Eλq​\(𝒙\)\],u\_\{\\bm\{\\eta\},\\bm\{c\}\}\(\\bm\{x\}\)=\\Re\\left\[\\sum\_\{k=1\}^\{K\}f\_\{\\eta\_\{k\}\}\(\\zeta\_\{k\}\)\+\\sum\_\{q=1\}^\{Q\}c\_\{q\}E\_\{\\lambda\_\{q\}\}\(\\bm\{x\}\)\\right\],\(6\)
whereζk=ak⋅𝒙^\+bk\\zeta\_\{k\}=a\_\{k\}\\cdot\\hat\{\\bm\{x\}\}\+b\_\{k\}, with a complex offsetbkb\_\{k\}and𝒙^=\(𝒙−𝒙0\)/L0\\hat\{\\bm\{x\}\}=\(\\bm\{x\}\-\\bm\{x\}\_\{0\}\)/L\_\{0\}\([Section4](https://arxiv.org/html/2610.12541#S4.SS0.SSS0.Px2)\)\. Each branch is harmonic, and so is their sum\. The normal derivative follows directly from the chain rule\. For an outward unit normal𝒏∈ℝ3\\bm\{n\}\\in\\mathbb\{R\}^\{3\},

∂nu𝜼,𝒄​\(𝒙\)=ℜ⁡\[1L0​∑k=1K\(𝒏⋅ak\)​fηk′​\(ζk\)\+∑q=1Qcq​𝒏⋅∇𝒙Eλq​\(𝒙\)\]\.\\partial\_\{n\}u\_\{\\bm\{\\eta\},\\bm\{c\}\}\(\\bm\{x\}\)=\\Re\\left\[\\frac\{1\}\{L\_\{0\}\}\\sum\_\{k=1\}^\{K\}\(\\bm\{n\}\\cdot a\_\{k\}\)f^\{\\prime\}\_\{\\eta\_\{k\}\}\(\\zeta\_\{k\}\)\+\\sum\_\{q=1\}^\{Q\}c\_\{q\}\\bm\{n\}\\cdot\\nabla\_\{\\bm\{x\}\}E\_\{\\lambda\_\{q\}\}\(\\bm\{x\}\)\\right\]\.\(7\)The model is fitted with the same boundary loss as in 2D\.

Since each network branch is entire and the enrichment terms are harmonic away from their branch cuts, the full representation is harmonic inΩ\\Omega\([AppendixA](https://arxiv.org/html/2610.12541#A1)\)\. For pure Dirichlet problems,ΓD=∂Ω\\Gamma\_\{D\}=\\partial\\Omega, this bounds the interior error by the boundary misfit\. The errore=u𝜼,𝒄−u⋆e=u\_\{\\bm\{\\eta\},\\bm\{c\}\}\-u^\{\\star\}is harmonic, and the junctionsJJ, points in 2D and edges in 3D, form a polar set\. Ifeeis bounded and continuous up to∂Ω∖J\\partial\\Omega\\setminus J, the extended maximum principle \([Ransford, 1995](https://arxiv.org/html/2610.12541#bib.bib33), Thm\. 3\.6\.9;[Armitage and Gardiner, 2001](https://arxiv.org/html/2610.12541#bib.bib34)in 3D\) gives

supΩ\|u𝜼,𝒄−u⋆\|≤sup∂Ω∖J\|u𝜼,𝒄−gD\|\.\\sup\_\{\\Omega\}\\left\|u\_\{\\bm\{\\eta\},\\bm\{c\}\}\-u^\{\\star\}\\right\|\\;\\leq\\;\\sup\_\{\\partial\\Omega\\setminus J\}\\left\|u\_\{\\bm\{\\eta\},\\bm\{c\}\}\-g\_\{D\}\\right\|\.\(8\)The error is bounded when every enrichment term is, which rules out aln⁡\|w\|\\ln\|w\|component in the value\-jump terms \([AppendixA](https://arxiv.org/html/2610.12541#A1)\)\. We evaluate the right\-hand side on a finite boundary sample, which can only underestimate the supremum, so the reported value is an error estimate rather than a certified bound\.

## 4Model Setup

In both dimensions the model has the same three parts: a holomorphic core, the enrichment terms that[Table1](https://arxiv.org/html/2610.12541#S3.T1)assigns to the junctions, and a real part that turns their sum into a potential\. The core is a small complex\-valued network in which each layer applies a complex affine map followed by the complex exponential,

h0=ζ,hℓ\+1=exp⁡\(Wℓ​hℓ\+bℓ\),fη​\(ζ\)=Wout​hm\+bout,h\_\{0\}=\\zeta,\\qquad h\_\{\\ell\+1\}=\\exp\(W\_\{\\ell\}h\_\{\\ell\}\+b\_\{\\ell\}\),\\qquad f\_\{\\eta\}\(\\zeta\)=W\_\{\\mathrm\{out\}\}h\_\{m\}\+b\_\{\\mathrm\{out\}\},\(9\)whereζ\\zetais a complex scalar,mmis the number of hidden layers, the exponential acts on each component, andη\\etacollects the complex weights and biases\. The exponential is entire, sofηf\_\{\\eta\}is entire for any value ofη\\eta\.

#### 2D model

In two dimensions the input is the normalized coordinatez^=\(z−z0\)/L0\\hat\{z\}=\(z\-z\_\{0\}\)/L\_\{0\}, withz0z\_\{0\}the centre of the domain andL0L\_\{0\}its characteristic radius\. The core hasm=2m=2hidden layers of width 16 with complex exponential activations\. The enrichment terms are evaluated in the physical coordinate, since their sites, cut directions and interior angles come from the geometry of the problem, and one term is placed per junction and ladder index\. All 2D problems have unit length scale, sowwneeds no separate reference length\.

#### 3D model

The 3D model first projects the point onto several complex coordinates\.[Figure2](https://arxiv.org/html/2610.12541#S4.F2)shows how a point𝒙\\bm\{x\}becomes a potentialu⁡\(𝒙\)u\(\\bm\{x\}\)\.

𝒙\\bm\{x\}centre, scaleand projectfη1​\(ζ1\)f\_\{\\eta\_\{1\}\}\(\\zeta\_\{1\}\)⋮\\vdotsfηK​\(ζK\)f\_\{\\eta\_\{K\}\}\(\\zeta\_\{K\}\)Σ\\Sigma\+\+enrichment∑qcq​Eλq​\(𝒙\)\\sum\_\{q\}c\_\{q\}E\_\{\\lambda\_\{q\}\}\(\\bm\{x\}\)ℜ⁡\[⋅\]\\Re\[\\,\\cdot\\,\]u⁡\(𝒙\)u\(\\bm\{x\}\)ζ1\\zeta\_\{1\}ζK\\zeta\_\{K\}

Figure 2:Forward pass of the 3D model\. The point is projected ontoKKcomplex coordinatesζk=ak⋅𝒙^\+bk\\zeta\_\{k\}=a\_\{k\}\\cdot\\hat\{\\bm\{x\}\}\+b\_\{k\}along null directions, each coordinate is passed through its own holomorphic network, and the outputs are summed\. The enrichment terms are evaluated at the physical point and added before the real part is taken\.The point is centred and scaled,𝒙^=\(𝒙−𝒙0\)/L0\\hat\{\\bm\{x\}\}=\(\\bm\{x\}\-\\bm\{x\}\_\{0\}\)/L\_\{0\}, and projected ontoKKcomplex coordinatesζk=ak⋅𝒙^\+bk\\zeta\_\{k\}=a\_\{k\}\\cdot\\hat\{\\bm\{x\}\}\+b\_\{k\}, with a complex offsetbkb\_\{k\}\. By[Equation5](https://arxiv.org/html/2610.12541#S3.E5), a holomorphic function ofζk\\zeta\_\{k\}is harmonic in𝒙\\bm\{x\}whenak⋅ak=0a\_\{k\}\\cdot a\_\{k\}=0, so we write each direction as

a⁡\(ξ,ρ\)=ρ⁡\(1−ξ22,i⁡\(1\+ξ2\)2,ξ\),ξ,ρ∈ℂ,a\(\\xi,\\rho\)=\\rho\\left\(\\frac\{1\-\\xi^\{2\}\}\{2\},\\;\\frac\{i\(1\+\\xi^\{2\}\)\}\{2\},\\;\\xi\\right\),\\qquad\\xi,\\rho\\in\\mathbb\{C\},\(10\)which is null for everyξ\\xiandρ\\rho\. The directions can therefore be trained freely\. The same length scale is used for all coordinates, since a separate scale per coordinate would break the null condition\. Eachζk\\zeta\_\{k\}is passed through its own network of the form[Equation9](https://arxiv.org/html/2610.12541#S4.E9), with two hidden layers of width 8, and theK=16K=16outputs are added\. This sum is a discrete counterpart of the integral over directions used by[Ballini et al\. \(2026\)](https://arxiv.org/html/2610.12541#bib.bib13), with separate parameters for each direction instead of one network shared across them\. The number and size of the branches are the same for all problems \([AppendixB](https://arxiv.org/html/2610.12541#A2)\)\.

The enrichment terms are evaluated at the physical point, since their location is part of the geometry\. A straight edge uses the termwλw^\{\\lambda\}of[Equation2](https://arxiv.org/html/2610.12541#S3.E2), with the exponent that[Table1](https://arxiv.org/html/2610.12541#S3.T1)gives for its dihedral angle and the conditions on its two faces\. A circular conical vertex uses the separated mode

Sν​\(R,ϑ\)=Rν​Pν​\(cos⁡ϑ\),S\_\{\\nu\}\(R,\\vartheta\)=R^\{\\nu\}P\_\{\\nu\}\(\\cos\\vartheta\),\(11\)whereRRis the distance from the vertex,ϑ\\varthetathe angle from the cone axis, andν\\nuthe first positive root ofPν​\(cos⁡ϑwall\)=0P\_\{\\nu\}\(\\cos\\vartheta\_\{\\mathrm\{wall\}\}\)=0\. A trainedν\\nuenters both factors, andPνP\_\{\\nu\}is evaluated by a 400\-term hypergeometric series whose tail is below10−2710^\{\-27\}on the domain\. Site positions, edge frames and branch cuts are fixed by the geometry and are not trained, since moving them could place a singular line inside the domain\.[AppendixH](https://arxiv.org/html/2610.12541#A8)writes out the model for the L\-shaped prism\.

### 4\.1Parameter fitting

The model is fitted to boundary samples only\. Both loss terms of[Equation4](https://arxiv.org/html/2610.12541#S3.E4)are made dimensionless by a characteristic potentialU0U\_\{0\}and lengthL0L\_\{0\}: the Dirichlet term is divided byU02U\_\{0\}^\{2\}and the Neumann term multiplied byL02/U02L\_\{0\}^\{2\}/U\_\{0\}^\{2\}; both scales are one on the sampled suite\. Enrichment coefficients start at zero and exponents at the values of[Table1](https://arxiv.org/html/2610.12541#S3.T1)\(the conical degree at1/21/2\); exponents are fixed on Motz’s problem, the cracked beam and the sampled suite\. All parameters are optimized jointly withAdamfollowed byL\-BFGS\.

Every sample carries gradient onλ\\lambdathrough∂λwλ=wλ​Log⁡w\\partial\_\{\\lambda\}w^\{\\lambda\}=w^\{\\lambda\}\\operatorname\{Log\}w, but far from a site the change is smooth and the rest of the model can absorb it, so we expect uniform boundary samples to carry little information about the exponent\. Since the coefficients start at zero, an exponent also receives no gradient until its coefficient is nonzero\. We therefore add, except on the sampled suite, samples whose distance to each junction or edge is drawn log\-uniformly between10−410^\{\-4\}and10−110^\{\-1\}of the local length\. Step counts and sample sizes are given in[AppendicesF](https://arxiv.org/html/2610.12541#A6)and[J](https://arxiv.org/html/2610.12541#A10)\.

### 4\.2Evaluation

Where a reference solutionu⋆u^\{\\star\}is available, we report the relativeL2L^\{2\}error

ϵu​\(X\)=\(∑𝒙∈X\|u⁡\(𝒙\)−u⋆​\(𝒙\)\|2∑𝒙∈X\|u⋆​\(𝒙\)\|2\)1/2,\\epsilon\_\{u\}\(X\)=\\left\(\\frac\{\\sum\_\{\\bm\{x\}\\in X\}\|u\(\\bm\{x\}\)\-u^\{\\star\}\(\\bm\{x\}\)\|^\{2\}\}\{\\sum\_\{\\bm\{x\}\\in X\}\|u^\{\\star\}\(\\bm\{x\}\)\|^\{2\}\}\\right\)^\{1/2\},\(12\)estimated on a setXXof interior points\. For the CZT detector, where the reference is numerical, we report the root\-mean\-square error instead\. Where no closed\-form solution exists,u⋆u^\{\\star\}is a singular\-series or numerical reference, described with each problem\. Unless stated otherwise, results are medians over three random seeds\.

## 5Experiments in two dimensions

### 5\.1Problems

We use six named benchmarks, each with a single singular site \([AppendixC](https://arxiv.org/html/2610.12541#A3)\)\. Four have a closed\-form solution: two circular sectors, an L\-shaped domain and a slit disk, whose solution is the leading wedge mode of the opening angle\([Grisvard, 1985](https://arxiv.org/html/2610.12541#bib.bib3)\)\. The other two, Motz’s problem\([Motz, 1947](https://arxiv.org/html/2610.12541#bib.bib6)\)and the cracked beam\([Lu et al\., 2004](https://arxiv.org/html/2610.12541#bib.bib7)\), are scored against their truncated singular expansion, whose leading coefficientD0D\_\{0\}is known to 17 significant digits\([Lu et al\., 2004](https://arxiv.org/html/2610.12541#bib.bib7)\)\.

As an application we use a CdZnTe \(CZT\) drift strip radiation detector\([Kuvvetli et al\., 2010](https://arxiv.org/html/2610.12541#bib.bib35)\)geometry\. CZT detectors are semiconductor gamma\-ray detectors that operate at room temperature, and the electric potential inside the crystal determines how the charge released by a photon drifts to the electrodes\. The electrodes on the crystal faces are held at fixed potentials and separated by insulating gaps, so every electrode edge is a junction between a Dirichlet and a Neumann boundary\. A single detector therefore contains many singular sites of the same type\. In 2D we use a cross\-section of the detector, scored by the root\-mean\-square error against a numerical reference \([AppendixD](https://arxiv.org/html/2610.12541#A4)\)\. To test many sites on simpler geometries, we also use a suite of 24 randomly sampled problems with 4 to 12 sites each \([AppendixE](https://arxiv.org/html/2610.12541#A5)\)\.

### 5\.2Results

Baselines We compare three representations\. TheMultiLayer Perceptron \(MLP\)PINNis a real\-valued network trained on the Laplace residual at interior points and on the boundary misfit, with the weighting tuned per problem\. The plain holomorphic model is[Equation3](https://arxiv.org/html/2610.12541#S3.E3)withQ=0Q=0, and the enriched model adds the terms of[Table1](https://arxiv.org/html/2610.12541#S3.T1)\. The last two are the same network trained with the same budget and restarts, so they differ only in the enrichment\.

Table 2:Errors in 2D\. All entries reportϵu\\epsilon\_\{u\}, median over three seeds at the best budget of[AppendixF](https://arxiv.org/html/2610.12541#A6), chosen by the same rule for every method\.*Plain*and*Enriched*are the same network without and with the terms of[Table1](https://arxiv.org/html/2610.12541#S3.T1)\.The plain holomorphic model improves on theMLPPINNby a factor of 9 to 315 on five of the six benchmarks \([Table2](https://arxiv.org/html/2610.12541#S5.T2)\)\. The two configurations also differ in precision, architecture and boundary sampling \([AppendixF](https://arxiv.org/html/2610.12541#A6)\), so this compares solver configurations rather than isolating exact harmonicity\. The exception is the slit disk, where the plain model is worse than theMLPPINN: its solution vanishes on both slit faces but its transverse Cartesian derivative has opposite one\-sided limits\. An entire core extends smoothly across the slit and cannot reproduce this discontinuity, but the branch enrichment represents it directly\.

![Refer to caption](https://arxiv.org/html/2610.12541v1/enriched_new.png)

![Refer to caption](https://arxiv.org/html/2610.12541v1/error_enriched_new.png)

Figure 3:Enriched model on the 2D CZT detector\. Left: normalized potential\. Right: error against the referenceThe enrichment gives a further 4 to 14 orders of magnitude, and the enriched model reaches roundoff on the7​π/47\\pi/4sector and the L\-shape\. On these closed\-form benchmarks the solution is itself a dictionary term, so they test recovery of the supplied mode\. On Motz’s problem and the cracked beam it recoversD0D\_\{0\}to relative errors of8\.2×10−78\.2\\times 10^\{\-7\}and2\.0×10−62\.0\\times 10^\{\-6\}, with the network supplying the correction that reduces the error of the same four singular terms fitted alone,5\.5×10−35\.5\\times 10^\{\-3\}and2\.2×10−22\.2\\times 10^\{\-2\}\([AppendixC](https://arxiv.org/html/2610.12541#A3)\), to the reported values\. With the same network and budget, tripling theAdamsteps and boundary samples improves the plain model by less than an order of magnitude, while the enriched model improves by several \([AppendixG](https://arxiv.org/html/2610.12541#A7)\)\. On the sampled suite, the enrichment improves all 72 paired runs and reduces the medianϵu\\epsilon\_\{u\}by a factor of3\.2×1063\.2\\times 10^\{6\}on the four closed\-form instances\.

On the detector \([Figure3](https://arxiv.org/html/2610.12541#S5.F3)\), the plain model reduces the error of theMLPPINNby a factor of 46, and the enrichment reduces it by a further factor of 1\.7\. The enrichment gain is much smaller than on the benchmarks: the terms of many junctions must be fitted together with the smooth field between the electrodes\. The enriched model is nevertheless 78 times more accurate than theMLPPINNon this detector cross\-section\.

On the pure Dirichlet benchmarks, the sampled right\-hand side of[Equation8](https://arxiv.org/html/2610.12541#S3.E8)is close to the true error: wherever the error lies above roundoff, the largest interior error of the enriched model reaches 94% to 100% of it, making it a useful estimate when no reference solution is available\.

## 6Experiments in three dimensions

### 6\.1Problems

We group the problems as in 2D \([AppendixI](https://arxiv.org/html/2610.12541#A9)\)\. Four have a closed\-form solution\. The L\-shaped prism has the cross\-section of a re\-entrant wedge and the solutiong⁡\(x,y\)=r2/3​sin⁡\(2​θ/3\)g\(x,y\)=r^\{2/3\}\\sin\(2\\theta/3\), which does not depend on the axial coordinate\. Two further fields on the same prism vary along the edge,

uA​\(x,y,z\)=\(1\+0\.5​z\)​g​\(x,y\)\+0\.1​x​y​z,uB​\(𝒙\)=uA​\(𝒙\)\+0\.1‖𝒙−\(−2,−2,2\)‖,u\_\{A\}\(x,y,z\)=\(1\+0\.5z\)g\(x,y\)\+0\.1xyz,\\qquad u\_\{B\}\(\\bm\{x\}\)=u\_\{A\}\(\\bm\{x\}\)\+\\frac\{0\.1\}\{\\\|\\bm\{x\}\-\(\-2,\-2,2\)\\\|\},\(13\)and a smooth harmonic field on the unit cube has no singular site at all, so it tests the core alone\. Two problems are scored against a series reference: the conical notch, whose solution is the vertex mode of[Equation11](https://arxiv.org/html/2610.12541#S4.E11), and the extruded Motz problem, which adds zero\-flux end faces to the planar problem and uses a 60\-term singular series\. As an application we use the full 3D CZT detector, posed in millimetres and volts, whose 238 electrode edges are all Dirichlet to Neumann junctions\. It is scored against a graded finite\-element reference whose own error is about10−410^\{\-4\}\.

Table 3:Relative errorϵu\\epsilon\_\{u\}in 3D, median over three converged seeds at each method’s best schedule\. The detector is scored against the finite\-element reference\.‡Extension of WAN to mixed boundary data\.
### 6\.2Results

#### Baselines

We compare against four published mesh\-free methods that impose the equation in different ways\. ThePINN\([Raissi et al\., 2019](https://arxiv.org/html/2610.12541#bib.bib15)\)minimizes the strong\-form residual, the Deep Ritz method\([E and Yu, 2018](https://arxiv.org/html/2610.12541#bib.bib30)\)minimizes the Dirichlet energy, and the weak adversarial network \(WAN\) of[Zang et al\. \(2020\)](https://arxiv.org/html/2610.12541#bib.bib31)minimizes the weak form against an adversarially trained test function\. The holomorphic network of[Ballini et al\. \(2026\)](https://arxiv.org/html/2610.12541#bib.bib13)is exactly harmonic but has no singular terms\. Each baseline follows its published protocol, with the authors’ code where it runs on our problems, and is also run with its steps and boundary samples multiplied by 1\.5, 2 and 3\. We report the best implementation and budget of each method \([AppendixJ](https://arxiv.org/html/2610.12541#A10)\)\.

None of the four baselines reaches an error below8×10−48\\times 10^\{\-4\}on any problem with a singularity \([Table3](https://arxiv.org/html/2610.12541#S6.T3)\), and the exactly harmonic network of[Ballini et al\. \(2026\)](https://arxiv.org/html/2610.12541#bib.bib13)is not consistently more accurate than thePINN\. As in 2D, exact harmonicity is not sufficient without the singular terms\. With them, the enriched model outperforms the best baseline by 5 or more orders of magnitude on the L\-prism and the conical notch, whose solution is the supplied mode, and by three orders of magnitude on the extruded Motz problem\.

![Refer to caption](https://arxiv.org/html/2610.12541v1/depth_slice_u_new.png)

![Refer to caption](https://arxiv.org/html/2610.12541v1/depth_slice_error_new.png)

Figure 4:Enriched model on the 3D CZT detector, planez=21z=21mm\. Left: normalized potential\. Right: error against the reference, scaled by the root mean square of the reference solution\.The gain is smaller when the solution varies along the edge\. The edge terms carry a constant amplitude, so the network must carry the axial variation ofuAu\_\{A\}anduBu\_\{B\}, and the margin falls to a factor of 5 to 7\. On the smooth cube there is nothing to enrich, and the enriched model is three times more accurate than thePINN\.

On the detector, the enriched model reaches1\.3×10−21\.3\\times 10^\{\-2\}, against2\.7×10−22\.7\\times 10^\{\-2\}for Deep Ritz and7\.5×10−27\.5\\times 10^\{\-2\}for thePINN\. The network of[Ballini et al\. \(2026\)](https://arxiv.org/html/2610.12541#bib.bib13)reaches2\.4×10−22\.4\\times 10^\{\-2\}\. The smaller margin may reflect the 238 edge terms fitted together with the field of the whole crystal, the unenriched strip ends, or conditioning\. It is still the lowest potential error of the tested methods on this geometry\.

## 7Conclusion

We have shown that building the analytical structure of the Laplace equation directly into the hypothesis space is more effective than enforcing it through a residual penalty\. By leveraging holomorphic network whose real part is harmonic by construction and using a set of singular enrichment terms we remove the need to use the interior points residual during training and handle difficult local behaviours that smooth neural methods cannot handle correctly\. Furthermore, we extend the formulation to three dimensions through null complex projections\. Together they reduce errors by orders of magnitude over both residual\-based and exactly harmonic baselines across a range of 2D and 3D problems, with the largest gains where the dictionary captures the solution’s singular behaviour exactly\.

The advantage shrinks when many junctions share the domain, when edge amplitudes vary axially, or when the domain geometry falls outside the tested configurations\. Several extensions follow naturally from the construction\. A dictionary for general polyhedral vertices would extend the 3D results beyond the edge and conical cases\. Edge terms with an affine axial amplitude would close the remaining gap on problems where the singular strength varies along the edge\. More broadly, since the Laplace equation is linear, a single model withNNoutput heads sharing the same holomorphic core and enrichment structure could representNNsolutions on the same geometry simultaneously, as the junction structure and branch cuts are identical for all of them and only the smooth remainder and enrichment coefficients differ\. For the CZT detector, where 100 weighting potentials share the same electrode geometry, this would reduce training cost by an order of magnitude while exploiting the shared singular structure across all outputs\.

### Acknowledgements

This work originated as a course project at the University of Amsterdam\. We thank Chris Hoes, Remco Hogerwerf, Finn Govers, Storm Koper, Gerrit Janssen and Stef de Wildt, who were part of that project and contributed to the initial ideation and code development\. This work was partially supported by thei\-RASEproject, funded by the European Union under the Horizon Europe EIC Pathfinder Open programme \(Grant Agreement No\. 101130550\)\.

## References

- Armitage and Gardiner \(2001\)D\. H\. Armitage and S\. J\. GardinerClassical potential theory\.Springer Monographs in Mathematics,Springer,London\.External Links:[Document](https://dx.doi.org/10.1007/978-1-4471-0233-5)Cited by:[§3\.4](https://arxiv.org/html/2610.12541#S3.SS4.SSS0.Px1.p4.1)\.
- Axleret al\.\(2001\)S\. Axler, P\. Bourdon, and W\. RameyHarmonic function theory\.2 edition,Graduate Texts in Mathematics, Vol\.137,Springer\.External Links:ISBN 9780387952185Cited by:[§2](https://arxiv.org/html/2610.12541#S2.p1.1),[§3\.2](https://arxiv.org/html/2610.12541#S3.SS2.p1.1)\.
- Balliniet al\.\(2026\)E\. Ballini, A\. P\. Engsig\-Karup, and T\. AndriolloA holomorphic neural network framework for 3D boundary value problems governed by harmonic potentials\.arXiv preprint arXiv:2605\.31231\.Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p3.1),[§2](https://arxiv.org/html/2610.12541#S2.p1.1),[§4](https://arxiv.org/html/2610.12541#S4.SS0.SSS0.Px2.p2.2),[§6\.2](https://arxiv.org/html/2610.12541#S6.SS2.SSS0.Px1.p1.1),[§6\.2](https://arxiv.org/html/2610.12541#S6.SS2.SSS0.Px1.p2.1),[§6\.2](https://arxiv.org/html/2610.12541#S6.SS2.SSS0.Px1.p4.1)\.
- Calafàet al\.\(2024\)M\. Calafà, E\. Hovad, A\. P\. Engsig\-Karup, and T\. AndriolloPhysics\-informed holomorphic neural networks \(PIHNNs\): solving 2D linear elasticity problems\.Computer Methods in Applied Mechanics and Engineering432,pp\. 117406\.External Links:[Document](https://dx.doi.org/10.1016/j.cma.2024.117406)Cited by:[Appendix F](https://arxiv.org/html/2610.12541#A6.p1.1),[§1](https://arxiv.org/html/2610.12541#S1.p3.1),[§2](https://arxiv.org/html/2610.12541#S2.p1.1),[§3\.4](https://arxiv.org/html/2610.12541#S3.SS4.p1.1)\.
- Calafàet al\.\(2025\)M\. Calafà, H\. M\. Jensen, and T\. AndriolloSolving plane crack problems via enriched holomorphic neural networks\.Engineering Fracture Mechanics322,pp\. 111133\.External Links:[Document](https://dx.doi.org/10.1016/j.engfracmech.2025.111133)Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p3.1),[§2](https://arxiv.org/html/2610.12541#S2.p2.1),[§3\.3](https://arxiv.org/html/2610.12541#S3.SS3.SSS0.Px1.p1.1),[§3\.3](https://arxiv.org/html/2610.12541#S3.SS3.p1.1)\.
- E and Yu \(2018\)W\. E and B\. YuThe Deep Ritz method: a deep learning\-based numerical algorithm for solving variational problems\.Communications in Mathematics and Statistics6\(1\),pp\. 1–12\.External Links:[Document](https://dx.doi.org/10.1007/s40304-018-0127-z)Cited by:[§6\.2](https://arxiv.org/html/2610.12541#S6.SS2.SSS0.Px1.p1.1)\.
- Evans \(2010\)L\. C\. EvansPartial differential equations\.2 edition,American Mathematical Society\.Cited by:[§3\.2](https://arxiv.org/html/2610.12541#S3.SS2.p1.1)\.
- Fairweather and Karageorghis \(1998\)G\. Fairweather and A\. KarageorghisThe method of fundamental solutions for elliptic boundary value problems\.Advances in Computational Mathematics9,pp\. 69–95\.External Links:[Document](https://dx.doi.org/10.1023/A%3A1018981221740)Cited by:[§2](https://arxiv.org/html/2610.12541#S2.p1.1)\.
- Fixet al\.\(1973\)G\. J\. Fix, S\. Gulati, and G\. I\. WakoffOn the use of singular functions with finite element approximations\.Journal of Computational Physics13\(2\),pp\. 209–228\.External Links:[Document](https://dx.doi.org/10.1016/0021-9991%2873%2990023-5)Cited by:[§2](https://arxiv.org/html/2610.12541#S2.p2.1)\.
- Ghoshet al\.\(2023\)A\. Ghosh, A\. A\. Gentile, M\. Dagrada, C\. Lee, S\. S\. Kim, H\. Cha, Y\. Choi, D\. Kim, J\. Kye, and V\. E\. ElfvingHarmonic neural networks\.InProceedings of the 40th International Conference on Machine Learning,Proceedings of Machine Learning Research, Vol\.202,pp\. 11340–11359\.External Links:[Link](https://proceedings.mlr.press/v202/ghosh23b.html)Cited by:[§2](https://arxiv.org/html/2610.12541#S2.p1.1),[§3\.4](https://arxiv.org/html/2610.12541#S3.SS4.p1.1)\.
- Gopal and Trefethen \(2019\)A\. Gopal and L\. N\. TrefethenSolving Laplace problems with corner singularities via rational functions\.SIAM Journal on Numerical Analysis57\(5\),pp\. 2074–2094\.External Links:[Document](https://dx.doi.org/10.1137/19M125947X)Cited by:[§2](https://arxiv.org/html/2610.12541#S2.p1.1)\.
- Grisvard \(1985\)P\. GrisvardElliptic problems in nonsmooth domains\.Pitman\.Note:Wedge exponents in chapters 4–5Cited by:[§2](https://arxiv.org/html/2610.12541#S2.p2.1),[§3\.3](https://arxiv.org/html/2610.12541#S3.SS3.p1.1),[§5\.1](https://arxiv.org/html/2610.12541#S5.SS1.p1.1)\.
- Heydariet al\.\(2019\)A\. A\. Heydari, C\. A\. Thompson, and A\. MehmoodSoftAdapt: techniques for adaptive loss weighting of neural networks with multi\-part loss functions\.External Links:1912\.12355,[Link](https://arxiv.org/abs/1912.12355)Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p2.1)\.
- Huet al\.\(2024\)T\. Hu, B\. Jin, and Z\. ZhouSolving Poisson problems in polygonal domains with singularity enriched physics informed neural networks\.SIAM Journal on Scientific Computing46\(4\),pp\. C369–C398\.External Links:[Document](https://dx.doi.org/10.1137/23M1601195)Cited by:[§2](https://arxiv.org/html/2610.12541#S2.p2.1),[§3\.3](https://arxiv.org/html/2610.12541#S3.SS3.SSS0.Px1.p1.1),[§3\.3](https://arxiv.org/html/2610.12541#S3.SS3.p1.1)\.
- Jagtap and Karniadakis \(2020\)A\. D\. Jagtap and G\. E\. KarniadakisExtended physics\-informed neural networks \(xpinns\): a generalized space\-time domain decomposition based deep learning framework for nonlinear partial differential equations\.Communications in Computational Physics28\(5\),pp\. 2002–2041\.Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p2.1)\.
- Jagtapet al\.\(2020\)A\. D\. Jagtap, E\. Kharazmi, and G\. E\. KarniadakisConservative physics\-informed neural networks on discrete domains for conservation laws: applications to forward and inverse problems\.Computer Methods in Applied Mechanics and Engineering365,pp\. 113028\.Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p2.1)\.
- Kendallet al\.\(2018\)A\. Kendall, Y\. Gal, and R\. CipollaMulti\-task learning using uncertainty to weigh losses for scene geometry and semantics\.External Links:1705\.07115,[Link](https://arxiv.org/abs/1705.07115)Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p2.1)\.
- Kuvvetliet al\.\(2010\)I\. Kuvvetli, C\. Budtz\-Jørgensen, E\. Caroli, and N\. AuricchioCZT drift strip detectors for high energy astrophysics\.Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment624\(2\),pp\. 486–491\.Note:New Developments in Radiation DetectorsExternal Links:ISSN 0168\-9002,[Document](https://dx.doi.org/https%3A//doi.org/10.1016/j.nima.2010.03.172),[Link](https://www.sciencedirect.com/science/article/pii/S0168900210008247)Cited by:[§5\.1](https://arxiv.org/html/2610.12541#S5.SS1.p2.1)\.
- Li and Lu \(2000\)Z\. C\. Li and T\. T\. LuSingularities and treatments of elliptic boundary value problems\.Mathematical and Computer Modelling31\(8\),pp\. 97–145\.External Links:[Document](https://dx.doi.org/10.1016/S0895-7177%2800%2900062-5)Cited by:[§2](https://arxiv.org/html/2610.12541#S2.p2.1)\.
- Liu and Nocedal \(1989\)D\. C\. Liu and J\. NocedalOn the limited memory BFGS method for large scale optimization\.Mathematical Programming45\(1–3\),pp\. 503–528\.Cited by:[Appendix F](https://arxiv.org/html/2610.12541#A6.p1.1)\.
- Liuet al\.\(2025\)Y\. Liu, H\. Gu, X\. Yu, and P\. QinDiminishing spectral bias in physics\-informed neural networks using spatially\-adaptive Fourier feature encoding\.Neural Networks182,pp\. 106886\.External Links:[Document](https://dx.doi.org/10.1016/j.neunet.2024.106886)Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p2.1)\.
- Liuet al\.\(2024\)Z\. Liu, Y\. Wang, S\. Vaidya, F\. Ruehle, J\. Halverson, M\. Soljacić, T\. Y\. Hou, and M\. TegmarkKAN: Kolmogorov–Arnold networks\.External Links:2404\.19756,[Link](https://arxiv.org/abs/2404.19756)Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p2.1)\.
- Luet al\.\(2004\)T\. T\. Lu, H\. Y\. Hu, and Z\. C\. LiHighly accurate solutions of Motz’s and the cracked beam problems\.Engineering Analysis with Boundary Elements28\(11\),pp\. 1387–1403\.External Links:[Document](https://dx.doi.org/10.1016/j.enganabound.2004.03.005)Cited by:[Table 4](https://arxiv.org/html/2610.12541#A3.T4),[§5\.1](https://arxiv.org/html/2610.12541#S5.SS1.p1.1)\.
- Luoet al\.\(2025\)K\. Luo, J\. Zhao, Y\. Wang, J\. Li, J\. Wen, J\. Liang, H\. Soekmadji, and S\. LiaoPhysics\-informed neural networks for PDE problems: a comprehensive review\.Artificial Intelligence Review58\.External Links:[Document](https://dx.doi.org/10.1007/s10462-025-11322-7)Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p1.1)\.
- McClenny and Braga\-Neto \(2023\)L\. D\. McClenny and U\. M\. Braga\-NetoSelf\-adaptive physics\-informed neural networks\.Journal of Computational Physics474,pp\. 111722\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2022.111722)Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p2.1)\.
- Motz \(1947\)H\. MotzThe treatment of singularities of partial differential equations by relaxation methods\.Quarterly of Applied Mathematics4\(4\),pp\. 371–377\.External Links:[Link](http://www.jstor.org/stable/43633577)Cited by:[§5\.1](https://arxiv.org/html/2610.12541#S5.SS1.p1.1)\.
- Rahamanet al\.\(2019\)N\. Rahaman, A\. Baratin, D\. Arpit, F\. Draxler, M\. Lin, F\. Hamprecht, Y\. Bengio, and A\. CourvilleOn the spectral bias of neural networks\.InProceedings of the 36th International Conference on Machine Learning,Proceedings of Machine Learning Research, Vol\.97,pp\. 5301–5310\.Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p2.1)\.
- Raissiet al\.\(2019\)M\. Raissi, P\. Perdikaris, and G\. E\. KarniadakisPhysics\-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations\.Journal of Computational Physics378,pp\. 686–707\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2018.10.045)Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p1.1),[§6\.2](https://arxiv.org/html/2610.12541#S6.SS2.SSS0.Px1.p1.1)\.
- Ransford \(1995\)T\. RansfordPotential theory in the complex plane\.London Mathematical Society Student Texts, Vol\.28,Cambridge University Press,Cambridge\.External Links:[Document](https://dx.doi.org/10.1017/CBO9780511623776)Cited by:[§3\.4](https://arxiv.org/html/2610.12541#S3.SS4.SSS0.Px1.p4.1)\.
- Schaback \(2008\)R\. SchabackSolving the Laplace equation by meshless collocation using harmonic kernels\.Advances in Computational Mathematics31,pp\. 457–470\.External Links:[Document](https://dx.doi.org/10.1007/s10444-008-9078-3)Cited by:[§2](https://arxiv.org/html/2610.12541#S2.p1.1)\.
- Sitzmannet al\.\(2020\)V\. Sitzmann, J\. N\. P\. Martel, A\. W\. Bergman, D\. B\. Lindell, and G\. WetzsteinImplicit neural representations with periodic activation functions\.Advances in Neural Information Processing Systems33\.Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p2.1)\.
- Wanget al\.\(2021\)S\. Wang, Y\. Teng, and P\. PerdikarisUnderstanding and mitigating gradient flow pathologies in physics\-informed neural networks\.SIAM Journal on Scientific Computing43\(5\),pp\. A3055–A3081\.External Links:[Document](https://dx.doi.org/10.1137/20M1318043)Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p2.1)\.
- Wanget al\.\(2022\)S\. Wang, X\. Yu, and P\. PerdikarisWhen and why PINNs fail to train: a neural tangent kernel perspective\.Journal of Computational Physics449,pp\. 110768\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2021.110768)Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p2.1)\.
- Wuet al\.\(2023\)C\. Wu, M\. Zhu, Q\. Tan, Y\. Kartha, and L\. LuA comprehensive study of non\-adaptive and residual\-based adaptive sampling for physics\-informed neural networks\.Computer Methods in Applied Mechanics and Engineering403,pp\. 115671\.External Links:[Document](https://dx.doi.org/10.1016/j.cma.2022.115671)Cited by:[§1](https://arxiv.org/html/2610.12541#S1.p2.1)\.
- Zanget al\.\(2020\)Y\. Zang, G\. Bao, X\. Ye, and H\. ZhouWeak adversarial networks for high\-dimensional partial differential equations\.Journal of Computational Physics411,pp\. 109409\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2020.109409)Cited by:[§6\.2](https://arxiv.org/html/2610.12541#S6.SS2.SSS0.Px1.p1.1)\.

## Appendix AHarmonicity and branch cuts

This section states when the model of[Section3](https://arxiv.org/html/2610.12541#S3)is harmonic and how the enrichment terms are placed so that the conditions hold\.

###### Lemma 1\(Harmonicity\)\.

Assume that each network branch is entire and that every enrichment term is twice continuously differentiable and harmonic inΩ\\Omega\. In 3D, also assumeak⋅ak=0a\_\{k\}\\cdot a\_\{k\}=0for every projection direction\. Then the representations in[Equations3](https://arxiv.org/html/2610.12541#S3.E3)and[6](https://arxiv.org/html/2610.12541#S3.E6)are harmonic inΩ\\Omega\.

###### Proof\.

In 2D, the real part of each holomorphic term is harmonic\. In 3D, each branch has zero Laplacian by[Equation5](https://arxiv.org/html/2610.12541#S3.E5), and the centring and scaling of the input only multiply it byL0−2L\_\{0\}^\{\-2\}\. The enrichment terms are harmonic by assumption, and taking real parts and summing preserves harmonicity by linearity of the Laplacian\. ∎

The network branches satisfy the lemma automatically; the enrichment terms satisfy it only if their branch cuts stay outside the domain\. We therefore orient each cut along a direction pointing out of the domain: the outward normal on a segment or face, and the exterior bisector at a corner \([Figure5](https://arxiv.org/html/2610.12541#A1.F5)\)\. On the convex domains \(the rectangles, the square and disk of the suite, the detector\) such a cut meets the closure of the domain only at its site\. On the sectors and the L\-shape it runs through the removed wedge, on the slit disk it lies on the slit, and on the L\-prism it lies in the removed quadrant\. The cone mode has no cut and is singular only on the rayϑ=π\\vartheta=\\pi, outside the domain\. Site positions and orientations are fixed during training\.

The cone mode is evaluated by the seriesPν​\(cos⁡ϑ\)=F12​\(−ν,ν\+1,1,\(1−cos⁡ϑ\)/2\)P\_\{\\nu\}\(\\cos\\vartheta\)=\{\}\_\{2\}F\_\{1\}\(\-\\nu,\\nu\+1;1;\(1\-\\cos\\vartheta\)/2\)truncated at 400 terms\. Its argument is at most0\.8540\.854on the domain, so the truncation error is below10−2710^\{\-27\}\.

The value\-jump term of[Table1](https://arxiv.org/html/2610.12541#S3.T1)isarg⁡w=ℜ⁡\(−i​Log⁡w\)\\arg w=\\Re\(\-i\\operatorname\{Log\}w\), which carries the jump in the data and stays bounded at the junction\. A general complex coefficient onLog⁡w\\operatorname\{Log\}wwould add a radial part, sinceℜ⁡\(c​Log⁡w\)=ℜ⁡\(c\)​ln\|w\|−ℑ⁡\(c\)​arg​w\\Re\(c\\operatorname\{Log\}w\)=\\Re\(c\)\\ln\|w\|\-\\Im\(c\)\\arg w, andln⁡\|w\|\\ln\|w\|is unbounded at the junction\. This is the component excluded in[Section3](https://arxiv.org/html/2610.12541#S3): with it, the error would be unbounded and[Equation8](https://arxiv.org/html/2610.12541#S3.E8)would not apply\.

Figure 5:Enrichment terms and branch cuts for the rows of[Table1](https://arxiv.org/html/2610.12541#S3.T1)\.\(a\)A Dirichlet to Neumann junction giveswλw^\{\\lambda\}\.\(b\)A Dirichlet value jump givesarg⁡w\\arg w, which carries the jump\.\(c\)A Neumann flux jump givesw​Log⁡ww\\operatorname\{Log\}w\.\(d\)A re\-entrant corner with equal Dirichlet values giveswλw^\{\\lambda\}; the equal\-flux Neumann row has the same term\. Each cut leaves the site along the exterior direction\.
## Appendix BChoice of architecture

The 2D network size and optimizer settings come from the ablation of[Table7](https://arxiv.org/html/2610.12541#A6.T7)\. The 3D size, 16 branches of two hidden layers of width eight, was fixed after a preliminary single\-seed screen of 135 configurations on four development problems, varying depth, width, activation, initialization, the number of branches and the treatment of the directions\. We keep one configuration per dimension for all problems\.

## Appendix CNamed benchmarks in two dimensions

[Table4](https://arxiv.org/html/2610.12541#A3.T4)defines the six benchmarks\. Each has a single singular site; every other junction has an integer exponent, so[Table1](https://arxiv.org/html/2610.12541#S3.T1)emits no term there\.

Table 4:The named benchmarks\.λ\\lambdais the leading exponent at the singular site\. The reference coefficients are from[Lu et al\. \(2004\)](https://arxiv.org/html/2610.12541#bib.bib7)\.The first four carry Dirichlet data from the exact solution on the whole boundary\. On the sectors, the data are zero on both radii andrλ​sin⁡\(λ​θ\)r^\{\\lambda\}\\sin\(\\lambda\\theta\),λ=π/ω\\lambda=\\pi/\\omega, on the arc\. On the L\-shape, they are zero on the two edges meeting at the re\-entrant corner andr2/3​sin⁡\(2​θ/3\)r^\{2/3\}\\sin\(2\\theta/3\)on the outer edges\. On the slit disk, they are zero on both faces of the slit andsin⁡\(θ/2\)\\sin\(\\theta/2\)on the rim, so thatu=r1/2​sin⁡\(θ/2\)u=r^\{1/2\}\\sin\(\\theta/2\)\.

Motz’s problem and the cracked beam have mixed data\. Both haveu=0u=0on the bottom edge forx<0x<0and zero flux there forx\>0x\>0\. Motz’s problem has zero flux on the top and left edges andu=500u=500on the right edge; the cracked beam has zero flux on the left and right edges andu=500u=500on the top edge\. Both are scored against the series∑i=059Di​ri\+1/2​cos⁡\(\(i\+12\)​θ\)\\sum\_\{i=0\}^\{59\}D\_\{i\}r^\{i\+1/2\}\\cos\(\(i\+\\frac\{1\}\{2\}\)\\theta\), which satisfies the bottom\-edge conditions exactly, with coefficients fitted by least squares on the other three edges\. It reproduces the publishedD0D\_\{0\}to relative errors of7×10−167\\times 10^\{\-16\}and1\.1×10−151\.1\\times 10^\{\-15\}\. Since the model’s exponents are fixed on these problems, its leading coefficient compares directly withD0D\_\{0\}\.[Table5](https://arxiv.org/html/2610.12541#A3.T5)shows how far the singular terms go without the network\.

Table 5:Least\-squares fit of the firstnnsingular terms without the network:ϵu\\epsilon\_\{u\}on 20,000 random interior points and the relative error ofD0D\_\{0\}\. The model usesn=4n=4\.![Refer to caption](https://arxiv.org/html/2610.12541v1/named_benchmarks_overview.png)Figure 6:The six benchmarks of[Table4](https://arxiv.org/html/2610.12541#A3.T4), with the exact solution drawn where one exists\.
## Appendix DCZT detector cross\-section

The 2D detector is a cross\-section, 40 mm wide and 5 mm thick, of the crystal of[AppendixI](https://arxiv.org/html/2610.12541#A9)\. One face carries 24 anode strips at 0 V, 0\.2 mm wide at a pitch of 1\.6 mm, with 25 groups of three drift strips at−50\-50,−100\-100and−50\-50V between them\. The opposite face is a single cathode at−350\-350V, and potentials are divided by−350\-350V\. The gaps and side walls carry zero flux, so each of the 198 strip edges is a Dirichlet to Neumann junction and receives one termwλw^\{\\lambda\}, withλ\\lambdatrained from1/21/2\. The reference is a finite\-difference solution, and[Table6](https://arxiv.org/html/2610.12541#A4.T6)gives the errors on its grid\.

Table 6:Potential errors on the CZT cross\-section, in units of 350 V, median over three seeds\.
## Appendix ESampled suite

The suite tests many sites on simple geometries\. It has 24 instances, twelve on the unit square and twelve on the unit disk, in three groups of eight: Dirichlet data with value jumps, mixed Dirichlet and Neumann data, and configurations with many segments \([Figure7](https://arxiv.org/html/2610.12541#A5.F7)\)\. The instances have 4 to 12 sites, with a median of 6\.5, and are drawn with seed 20260713; the tuning suite of[Table7](https://arxiv.org/html/2610.12541#A6.T7)is drawn separately with seed 20260804\. The four Dirichlet\-jump instances on the disk have a closed\-form solution and are scored byϵu\\epsilon\_\{u\}, the eight Dirichlet\-jump instances also by the sampled right\-hand side of[Equation8](https://arxiv.org/html/2610.12541#S3.E8), and all instances by their boundary loss\. Exponents are fixed at their dictionary values and no near\-site samples are used\.

Figure 7:The 24 instances of the sampled suite\. Dirichlet segments are coloured by their prescribed value, and Neumann segments are dashed\. Markers show the sites the dictionary emits, with the shape giving the row of[Table1](https://arxiv.org/html/2610.12541#S3.T1)\.
## Appendix FTraining in two dimensions

The model is implemented in JAX in double precision, since the reported errors go well below the10−710^\{\-7\}resolution of single precision\. Weights are drawn from the scaled complex He initialization of[Calafà et al\. \(2024\)](https://arxiv.org/html/2610.12541#bib.bib11)\. Training runs4×1044\\times 10^\{4\}full\-batch Adam steps at learning rate10−210^\{\-2\}with cosine decay, followed by at most6×1036\\times 10^\{3\}L\-BFGS iterations with a history of 100\([Liu and Nocedal, 1989](https://arxiv.org/html/2610.12541#bib.bib27)\), and keeps the best of three restarts\. Each segment receives 1,024 uniform samples, and each endpoint 128 further samples at log\-uniform distances between10−410^\{\-4\}and10−110^\{\-1\}of the segment length\. The sectors, the L\-shape, the slit disk and the detector train their exponents; Motz’s problem and the cracked beam fix them at1/21/2,3/23/2,5/25/2and7/27/2\.

The budget ladder multiplies the Adam steps and the samples per segment by 1, 1\.5, 2 and 3, and each rung is run with seeds 0, 1 and 2\.[Table2](https://arxiv.org/html/2610.12541#S5.T2)reports the best rung of each method\. All other settings were chosen on the tuning suite \([Table7](https://arxiv.org/html/2610.12541#A6.T7)\), which is disjoint from the reported problems\.

#### Baseline\.

The MLP PINN is a real\-valued tanh network with one hidden layer of width 128, trained in single precision on 5,000 interior points and400​b400bboundary points at budgetbb\. Its loss is the Laplace residual plus the boundary misfit, with boundary weight 10 on the sectors and 100 elsewhere, minimized by Adam at learning rate10−310^\{\-3\}followed by 6,000 L\-BFGS iterations, on the same budget ladder and seeds\.

Table 7:Ablation on the tuning suite \(six instances, three seeds, 18 runs per row\), changing one factor at a time from a starting configuration of width 32, depth 3, learning rate 0\.003, 20 000 Adam steps, 6 000 L\-BFGS steps, 512 points per segment, ladder depth 1 and fixed exponents\.*Ratio*is the median boundary loss relative to the starting configuration; below one is an improvement\.*Converged*counts runs with boundary loss below10−610^\{\-6\}\. The configuration of[Section4](https://arxiv.org/html/2610.12541#S4.SS0.SSS0.Px1)combines the five best changes\.

## Appendix GPlain against enriched

The plain and enriched models share the network, budget, restarts and seeds, and differ only in the singular terms\. On the named benchmarks \([Figure8](https://arxiv.org/html/2610.12541#A7.F8)\), tripling the budget improves the plain model by less than an order of magnitude and the enriched model by several\. On the sampled suite \([Figure9](https://arxiv.org/html/2610.12541#A7.F9)\), all 72 paired runs improve with enrichment\. The median paired ratio is1\.5×1031\.5\\times 10^\{3\}for the boundary loss,1\.1×1071\.1\\times 10^\{7\}for the sampled right\-hand side of[Equation8](https://arxiv.org/html/2610.12541#S3.E8)on the eight Dirichlet\-jump instances, and3\.2×1063\.2\\times 10^\{6\}forϵu\\epsilon\_\{u\}on the four closed\-form instances\. The gain is largest on value jumps and smallest on the many\-segment instances\.

To compare the error with the sampled bound of[Equation8](https://arxiv.org/html/2610.12541#S3.E8), each Dirichlet segment is sampled at2×1042\\times 10^\{4\}points, excluding points within10−310^\{\-3\}of a site\. On the two sectors, the largest error of the enriched model on a400×400400\\times 400interior grid reaches 94% to 100% of this bound; on the L\-shape and the7​π/47\\pi/4sector both are at roundoff\. For the plain model the ratio is 83% to 90% on all four\.

Figure 8:The plain and enriched models across the budget ladder on the named benchmarks, median over three seeds with the spread shaded\.Figure 9:Plain against enriched over the sampled suite, one point per instance and seed\. Points below the diagonal improved with enrichment; annotations give the median paired ratio\.#### Ladder depth\.

[Table8](https://arxiv.org/html/2610.12541#A7.T8)adds terms on the7​π/47\\pi/4sector, whose solution is the leading mode alone\. The leading exponent stays within2×10−82\\times 10^\{\-8\}of its true value at every depth, while the surplus exponents cluster and carry amplitudes that cancel in pairs\. The error stays near10−1010^\{\-10\}, so extra terms neither help nor harm here\. These runs use budget one, which is why the depth\-one error is above that of[Table2](https://arxiv.org/html/2610.12541#S5.T2)\.

Table 8:Ladder depth on the7​π/47\\pi/4sector,λ⋆=4/7\\lambda^\{\\star\}=4/7\.λ1\\lambda\_\{1\}is the trained exponent closest toλ⋆\\lambda^\{\\star\}, the gap is the smallest separation between trained exponents, and the spurious amplitude is the largest\|cq\|\|c\_\{q\}\|on a non\-leading term\.

## Appendix HThe 3D model on the L\-shaped prism

On the L\-shaped prism the only singular site is the re\-entrant edge, withω=3​π/2\\omega=3\\pi/2and zero Dirichlet data on both walls\.[Table1](https://arxiv.org/html/2610.12541#S3.T1)assigns it the termw2/3w^\{2/3\}, and the model is

u⁡\(𝒙\)=ℜ⁡\[∑k=116fηk​\(ζk\)\+c​w2/3\]\.u\(\\bm\{x\}\)=\\Re\\left\[\\sum\_\{k=1\}^\{16\}f\_\{\\eta\_\{k\}\}\(\\zeta\_\{k\}\)\+c\\,w^\{2/3\}\\right\]\.The exact solutionr2/3​sin⁡\(2​θ/3\)r^\{2/3\}\\sin\(2\\theta/3\)equalsℜ⁡\[−i​w2/3\]\\Re\[\-i\\,w^\{2/3\}\], soc=−ic=\-irepresents it exactly and the branches need only represent zero\. OnuAu\_\{A\}anduBu\_\{B\}the singular amplitude varies along the edge, and the branches must carry that variation\.

## Appendix IThree\-dimensional problems

#### Closed\-form solution\.

The L\-prism is the L\-shape of[Table4](https://arxiv.org/html/2610.12541#A3.T4)extruded over0<z<10<z<1, with Dirichlet data from the exact solution on every face;uAu\_\{A\}anduBu\_\{B\}use the same geometry\. The smooth cube is\[0,1\]3\[0,1\]^\{3\}with Dirichlet data from1/‖𝒙−\(1\.6,1\.6,1\.6\)‖1/\\\|\\bm\{x\}\-\(1\.6,1\.6,1\.6\)\\\|, whose pole lies outside the cube\.

#### Series reference\.

The conical notch is the unit ball with a cone of half\-angleπ/4\\pi/4removed around the positivezzaxis, so that, measured from the negative axis, the wall lies atϑ=3​π/4\\vartheta=3\\pi/4\. The data are zero on the wall and exact on the sphere, and the solution is[Equation11](https://arxiv.org/html/2610.12541#S4.E11)withν=0\.4630985617801066\\nu=0\.4630985617801066; the trained degree, started at1/21/2, reaches0\.4630990\.463099\. The extruded Motz problem has the planar Motz data on its cross\-section and zero flux on its end faces\. Its 60\-term series reference changes by5\.5×10−135\.5\\times 10^\{\-13\}when reduced to 40 terms and matches the publishedD0D\_\{0\}to7×10−167\\times 10^\{\-16\}\.

#### CZT detector\.

The detector is a40×5×4040\\times 5\\times 40mm crystal, posed in millimetres and volts\. One face carries the anode and drift strips of[AppendixD](https://arxiv.org/html/2610.12541#A4), running the full length of the crystal, and the opposite face 20 cathode strips at−350\-350V, 1\.8 mm wide at a pitch of 2 mm\. The gaps and the four sides carry zero flux, so each of the 119 electrodes has two Dirichlet to Neumann edges on a flat face, 238 in total, each with one edge term andλ\\lambdatrained from1/21/2; the strip ends are not enriched\. Zero flux idealizes an insulating surface, and space charge is neglected, so this is a numerical study of the operating potential rather than a validation of the device\. The reference is a graded finite\-element solution whose own error is about10−410^\{\-4\}\.

## Appendix JTraining and evaluation in three dimensions

[Table9](https://arxiv.org/html/2610.12541#A10.T9)lists the settings of the enriched model\. The core has 3,104 real parameters; with the directions and enrichment terms the model has between 3,200 and 4,390, depending on the number of sites\.[Table10](https://arxiv.org/html/2610.12541#A10.T10)gives the per\-seed errors behind[Table3](https://arxiv.org/html/2610.12541#S6.T3)\.

Table 9:Settings of the enriched model in 3D\.Table 10:Per\-seed relative errorϵu\\epsilon\_\{u\}of the enriched model in[Table3](https://arxiv.org/html/2610.12541#S6.T3)\.#### Baselines\.

Each baseline runs under its own published protocol: the architecture, optimizer, sampling and step counts of its paper, with the authors’ code where it runs on our problems\. It is also run with its steps and boundary samples multiplied by 1\.5, 2 and 3, and[Table3](https://arxiv.org/html/2610.12541#S6.T3)reports its lowest three\-seed median\. Deep Ritz uses three residual blocks of width ten with cubic ReLU activations and a boundary penalty of10310^\{3\}\. WAN is formulated for Dirichlet data; on the extruded Motz problem and the detector, the zero\-flux condition enters through the weak form, and we mark these entries as an extension\. Baselines that require unit scales are fitted on the detector with lengths divided by 5 mm and potentials by 350 V, and mapped back before scoring\.