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.
Figures & tables
Figure 1: Singular behaviour at a D → N junction with ω=2π/3 , λ=3/4 . Left: the leading mode rλsin(λθ) is bounded at the vertex. Right: the gradient magnitude ∣∇u∣∼rλ−1=r−1/4 diverges there. The enrichment term wλ represents this mode directly.
Junction
Local solution
∣∇u∣
Term Eλ
Dirichlet to Neumann
rλsin(λθ) , λ=(2k−1)π/(2ω)
rλ−1
wλ
Both Dirichlet, value jump Δg
(Δg/ω)θ
r−1
argw
Both Neumann, flux jump, ω∈πN
ℜ[wLogw]
∣lnr∣
wLogw
Both Neumann, equal flux, ω>π
rλcos(λθ) , λ=kπ/ω
rλ−1
wλ
Both Dirichlet, equal values, ω>π
rλsin(λθ) , λ=kπ/ω
rλ−1
wλ
Otherwise
C∞ away from vertex
bounded
none
Table 1: Enrichment dictionary. Each junction matches one row and contributes the term Eλ with a trainable coefficient cq . Local coordinates are centred at the junction zc , with w=z−zc=reiθ ; ω is the interior angle and k=1,2,… indexes the ladder. Every term is holomorphic, so the representation remains exactly harmonic. Figure 5 is a visualization of this table.
Figure 2: Forward pass of the 3D model. The point is projected onto K complex coordinates ζk=ak⋅x^+bk 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.
Problem
MLP PINN
Plain
Enriched
Closed-form solution
Sector 3π/2
9.5×10−2
7.7×10−4
3.3×10−12
Sector 7π/4
1.8×10−1
1.9×10−2
1.4×10−16
L-shape
7.4×10−2
7.5×10−4
2.9×10−16
Slit disk
1.8×10−1
3.2×10−1
2.8×10−12
Series reference
Table 2: Errors in 2D. All entries report ϵu , median over three seeds at the best budget of Appendix F , chosen by the same rule for every method. Plain and Enriched are the same network without and with the terms of Table 1 .
Figure 3: Enriched model on the 2D CZT detector. Left: normalized potential. Right: error against the reference
Problem
PINN
Deep Ritz
WAN
Ballini et al.
Enriched
Closed-form solution
L-prism
8.7×10−4
2.7×10−3
5.6×10−3
1.3×10−3
6.9×10−9
uA
8.9×10−4
4.5×10−3
1.2×10−2
1.9×10−3
1.8×10−4
uB
1.1×10−3
5.9×10−3
1.2×10−2
1.6×10−3
1.6×10−4
Smooth cube
1.5×10−5
5.9×10−4
2.2×10−3
2.0×10−4
4.6×10−6
Series reference
Table 3: Relative error ϵ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.
Figure 4: Enriched model on the 3D CZT detector, plane z=21 mm. Left: normalized potential. Right: error against the reference, scaled by the root mean square of the reference solution.
Appendix figures & tables12 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 5: Enrichment terms and branch cuts for the rows of Table 1 . (a) A Dirichlet to Neumann junction gives wλ . (b) A Dirichlet value jump gives argw , which carries the jump. (c) A Neumann flux jump gives wLogw . (d) A re-entrant corner with equal Dirichlet values gives wλ ; the equal-flux Neumann row has the same term. Each cut leaves the site along the exterior direction.
Problem
Domain
λ
Reference
Sector 3π/2
r<1 , 0<θ<3π/2
2/3
exact
Sector 7π/4
r<1 , 0<θ<7π/4
4/7
exact
L-shape
(−1,1)2∖[0,1]×[−1,0]
2/3
exact
Slit disk
r<1 cut along 0≤x<1 , y=0
1/2
exact
Motz
(−1,1)×(0,1)
1/2
D0=401.162453745234416
Cracked beam
(−1,1)×(0,1)
1/2
D0=540.56512271362749
Appendix
Table 4: The named benchmarks. λ is the leading exponent at the singular site. The reference coefficients are from Lu et al. (2004) .
Motz
Cracked beam
n
ϵu
D0 error
ϵu
D0 error
1
1.9×10−1
7.9×10−2
3.0×10−1
1.5×10−1
4
5.5×10−3
4.4×10−3
2.2×10−2
1.1×10−2
8
1.9×10−4
3.5×10−5
4.4×10−4
8.6×10−6
12
1.6×10−5
5.1×10−7
2.6×10−5
2.8×10−7
Appendix
Table 5: Least-squares fit of the first n singular terms without the network: ϵu on 20,000 random interior points and the relative error of D0 . The model uses n=4 .
Figure 6: The six benchmarks of Table 4 , with the exact solution drawn where one exists.
Model
MSE
RMSE
MLP PINN
2.40×10−1
4.90×10−1
Plain
1.18×10−4
1.09×10−2
Enriched
4.19×10−5
6.47×10−3
Appendix
Table 6: Potential errors on the CZT cross-section, in units of 350 V, median over three seeds.
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 Table 1 .
Configuration
Boundary loss
Ratio
Converged
Time (s)
starting configuration
6.4×10−5
1.00
3/18
402
depth 2
3.3×10−6
0.03
7/18
280
learning rate 0.01
8.1×10−6
0.14
7/18
381
Adam steps 40 000
1.4×10−5
0.41
5/18
644
width 16
1.4×10−5
0.42
7/18
275
points per segment 1 024
8.5×10−6
0.43
6/18
645
Appendix
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 below 10−6 . The configuration of Section 4 combines the five best changes.
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.
Terms
∣λ1−λ⋆∣
Smallest gap
Spurious ∣cq∣
ϵu
1
4.6×10−11
–
0
5.7×10−11
2
2.3×10−10
6.3×10−1
1.9×10−9
3.4×10−11
3
7.4×10−10
9.5×10−8
3.7×10−2
2.2×10−11
4
8.8×10−10
1.0×10−7
1.8×10−1
3.0×10−11
6
2.0×10−8
2.5×10−6
1.0×10−1
1.3×10−10
8
1.9×10−9
4.7×10−7
2.2×10−1
1.1×10−10
Appendix
Table 8: Ladder depth on the 7π/4 sector, λ⋆=4/7 . λ1 is the trained exponent closest to λ⋆ , the gap is the smallest separation between trained exponents, and the spurious amplitude is the largest ∣cq∣ on a non-leading term.
Setting
Value
Branches / hidden layers / width
16/2/8
Adam steps
120,000, learning rate 0.01 , cosine decay
Maximum L-BFGS iterations
4,000
Uniform boundary samples
1,536 per face (1,024 on the detector)
Near-site samples
128 per site, 10−4 to 10−1 of the local length
Restarts
3 per seed, selected on a fresh boundary sample
Appendix
Table 9: Settings of the enriched model in 3D.
Problem
Seed 0
Seed 1
Seed 2
L-prism
1.3×10−8
1.7×10−9
6.9×10−9
uA
1.8×10−4
1.9×10−4
1.6×10−4
uB
1.8×10−4
1.6×10−4
1.6×10−4
Smooth cube
4.4×10−6
4.9×10−6
4.6×10−6
Conical notch
4.5×10−9
6.0×10−9
6.1×10−9
Extruded Motz
3.4×10−6
2.1×10−6
4.6×10−6
Appendix
Table 10: Per-seed relative error ϵu of the enriched model in Table 3 .
Neural PDE solvers use trainable models to approximate solutions or solution operators for partial differential equations, often with physical constraints.
We present a neural-network-based framework for the solution of three-dimensional boundary value problems where the solution is expressible in terms of harmonic potentials. The approach leverages the Whittaker integral formula, which allows representing the solution through functions that are holomorphic with respect to a suitable complex variable. These functions are subsequently approximated using holomorphic neural networks, which guaranty fulfillment of the holomorphicity requirement. A key feature of the proposed formulation is that the governing partial differential equations (PDEs) are satisfied exactly by construction. Therefore, in contrast to standard physics-informed neural networks, no residual minimization of PDEs is required in the interior of the domain, and training is based exclusively on boundary collocation points. The method is validated against three-dimensional Laplace and linear elasticity problems, where, in the latter case, displacement and stress fields are expressed via the Papkovich-Neuber potentials. The numerical results show an accurate approximation of both scalar and vector fields, with errors remaining controlled throughout the domain. Overall, the work demonstrates that the incorporation of analytical structures into neural network architectures provides a natural and effective framework for the meshless approximation of three-dimensional boundary value problems while preserving the underlying properties of the governing equations.
Enrico Ballini, Allan Peter Engsig-Karup, Tito Andriollo
Department of Mechanical and Production Engineering, Aarhus University, Aarhus, Denmark · Department of Applied Mathematics and Computer Science, Technical University of Denmark, Kongens Lyngby, Denmark
We investigate the use of neural network solvers for infinity and p-Laplace problems, which are fundamental in nonlinear analysis and have practical applications. Our approach employs Physics-Informed Neural Networks (PINNs) and Deep Operator Networks (DeepONets) to address computational challenges associated with large p values, ranging from 2 to 1000, on various 2D and 3D domains. Our method offers advantages over traditional physics-based solvers, especially in three dimensions where mesh-based solvers become very costly for these problems. We also establish conditional convergence results for PINN approximations of both problems and a universal approximation result for DeepONet on the parametric p-Poisson problem. We demonstrate the effectiveness of these neural network solvers through numerical experiments and compare their performance with conventional methods.
Tak Shing Au Yeung, Ka Chun Cheung, Hannah Potgieter +2
NVIDIA AI Technology Center, NVIDIA, Hong Kong, China · Department of Mathematics, Simon Fraser University, Burnaby, British Columbia, Canada · NVIDIA AI Technology Center, NVIDIA, Singapore, Singapore
We introduce Neural Harmonic Measure Operator (NHMO), a neural solver for elliptic PDE problems on variable-shape domains. The harmonic measure of a domain is the boundary probability distribution that, integrated against any boundary data, returns the Dirichlet Laplace solution. It depends only on the geometry, not on the boundary data. NHMO parameterizes the density of this measure as a transformer-based boundary kernel supervised by Walk-on-Spheres exit samples, so one trained kernel handles different boundary values on a shape with no retraining. We extend it to Poisson via a classical decomposition, with an auxiliary network amortizing the source-induced correction and avoiding the singular volume quadrature that breaks direct evaluation. At inference, new boundary values and new sources both yield PDE solutions by re-integration against the fitted kernel and lift, with no retraining. NHMO improves over four prior baselines on the MCB-B 3D variable-shape Poisson benchmark across all five categories, and is competitive with major neural-operator baselines on a controlled 2D testbed.