BEACONS: Bounded-Error, Algebraically-Composable Neural Solvers for Partial Differential Equations

Jonathan Gorard; Ammar Hakim; James Juno

cs.LGFeb 16, 2026

BEACONS: Bounded-Error, Algebraically-Composable Neural Solvers for Partial Differential Equations

Authors: Jonathan Gorard, Ammar Hakim, James Juno

Organizations: Princeton University Princeton, NJ 08544, USA · Princeton Plasma Physics Laboratory Princeton, NJ 08540, USA

Abstract

The traditional limitations of neural networks in reliably generalizing beyond the convex hulls of their training data present a significant problem for computational physics, in which one often wishes to solve PDEs in regimes far beyond anything which can be experimentally or analytically validated. In this paper, we show how it is possible to circumvent these limitations by constructing formally-verified neural network solvers for PDEs, with rigorous convergence, stability, and conservation properties, whose correctness can therefore be guaranteed even in extrapolatory regimes. By using the method of characteristics to predict the analytical properties of PDE solutions a priori (even in regions arbitrarily far from the training domain), we show how it is possible to construct rigorous extrapolatory bounds on the worst-case L^inf errors of shallow neural network approximations. Then, by decomposing PDE solutions into compositions of simpler functions, we show how it is possible to compose these shallow neural networks together to form deep architectures, based on ideas from compositional deep learning, in which the large L^inf errors in the approximations have been suppressed. The resulting framework, called BEACONS (Bounded-Error, Algebraically-COmposable Neural Solvers), comprises both an automatic code-generator for the neural solvers themselves, as well as a bespoke automated theorem-proving system for producing machine-checkable certificates of correctness. We apply the framework to a variety of linear and non-linear PDEs, including the linear advection and inviscid Burgers' equations, as well as the full compressible Euler equations, in both 1D and 2D, and illustrate how BEACONS architectures are able to extrapolate solutions far beyond the training data in a reliable and bounded way. Various advantages of the approach over the classical PINN approach are discussed.

Figures & tables

Panel of Figure 1 — Figure 1: On the left, an example of the extended Racket data structure for representing a system of hyperbolic PDEs (in this case, the 1D compressible Euler equations), a collection of simulation parameters (in this case, representing a 1D Sod-type shock tube problem), and a collection of hyperparameters for a shallow BEACONS architecture (6 layers, 64 neurons per layer, 10,000 maximum training steps). On the right, the resulting optimized C code output by the automatic code-generator for running the specified simulation (using a Lax-Friedrichs solver with second-order flux extrapolation), generating the necessary training data, and training the BEACONS architecture accordingly.

Panel of Figure 2 — Figure 2: On the left, the final few steps of an executable Racket proof of flux continuity for a Roe solver for the 1D inviscid Burgers’ equation. On the right, the final few steps of an executable Racket proof of a bound on the worst-case ${L^{\infty}}$ error (for the non-smooth parts of the solution) for a shallow BEACONS architecture for solving the 1D inviscid Burgers’ equation.

Panel of Figure 3 — Figure 3: Results for the 1D linear advection Riemann problem at times ${t=0.5}$ (left) and ${t=0.8}$ (right), obtained using a formally-verified numerical solver (blue), a 6-layer neural network (red), an 8-layer neural network (yellow), a 6-layer BEACONS architecture (purple), and an 8-layer BEACONS architecture (green). The 6-layer neural network slightly overestimates the advection speed, and the 8-layer neural network significantly underestimates it; both neural network solutions exhibit substantial conservation errors. Both BEACONS solutions track the numerical solution more-or-less perfectly.

Architecture	${L^{\infty}}$ Error (Final)	${L^{2}}$ Error (Final)	${L^{\infty}}$ Error (All)	${L^{2}}$ Error (All)
6-layer NN	1.033641	2.536744	1.076130	3.002461
8-layer NN	0.976984	9.030859	1.002587	4.807408
6-layer BEACONS	0.612160	1.132093	0.782192	1.061877
8-layer BEACONS	0.605036	1.211529	0.633319	1.189765

Table 1:

{L^{2}}

and

{L^{\infty}}

error analysis for the 1D linear advection Riemann problem, comparing the 6-layer neural network, 8-layer neural network, 6-layer BEACONS, and 8-layer BEACONS solutions against the formally-verified numerical solution, both for the final predicted frame, and across all predicted frames.

Architecture	Conservation Error (Final)	Conservation Error (Total)
6-layer NN	-21.169941	498.706193
8-layer NN	-130.422471	-4391.640256
6-layer BEACONS	-3.770073	-155.321387
8-layer BEACONS	1.881948	49.570335

Table 2: Conservation analysis for the 1D linear advection Riemann problem, comparing the 6-layer neural network, 8-layer neural network, 6-layer BEACONS, and 8-layer BEACONS solutions against the formally-verified numerical solution, both for the final predicted frame, and across all frames.

Panel of Figure 4 — Figure 4: Results for the 2D linear advection disk problem at times ${t=0.33}$ (left), ${t=0.66}$ (middle), and ${t=0.99}$ (right), obtained using a formally-verified numerical solver (top), an 8-layer neural network (middle), and an 8-layer BEACONS architecture (bottom). The 8-layer neural network progressively distorts the shape of the disk as it approaches the top right of the domain, causing it to become increasingly “egg-shaped”, while also underestimating the advection speed overall. The BEACONS architecture successfully preserves the shape of the disk, and shows near-perfect agreement with the numerical solution.

Architecture	${L^{\infty}}$ Error (Final)	${L^{2}}$ Error (Final)	${L^{\infty}}$ Error (All)	${L^{2}}$ Error (All)
8-layer NN	1.031098	40.481536	1.031098	22.772318
8-layer BEACONS	0.864784	9.363477	0.938123	6.747000

Table 3:

{L^{2}}

and

{L^{\infty}}

error analysis for the 2D linear advection disk problem, comparing the 8-layer neural network and 8-layer BEACONS solutions against the formally-verified numerical solution, both for the final predicted frame, and across all predicted frames.

Panel of Figure 5 — Figure 5: Results for the 1D inviscid Burgers’ “top-hat” initial value problem at times ${t=0.5}$ (left) and ${t=0.8}$ (right), obtained using a formally-verified numerical solver (blue), a 6-layer neural network (red), an 8-layer neural network (yellow), a 6-layer BEACONS architecture (purple), and an 8-layer BEACONS architecture (green). Both the 6-layer and 8-layer neural networks significantly overestimate the speed of the right-moving shock and underestimate the speed of the left-moving rarefaction; both neural network solutions fail to conserve the advected quantity, and exhibit significant diffusion of the overall qualitative structure of the solution. Both BEACONS solutions track the numerical solution more-or-less perfectly, without significant instability, wave-speed mis-prediction, or conservation errors.

Architecture	${L^{\infty}}$ Error (Final)	${L^{2}}$ Error (Final)	${L^{\infty}}$ Error (All)	${L^{2}}$ Error (All)
6-layer NN	1.310069	15.058654	1.311499	8.144540
8-layer NN	1.208962	12.127816	1.307407	6.691768
6-layer BEACONS	0.986979	2.165441	0.973603	1.554421
8-layer BEACONS	1.028023	2.404486	1.014091	1.636904

Table 4:

{L^{2}}

and

{L^{\infty}}

error analysis for the 1D inviscid Burgers’ “top-hat” initial value problem, comparing the 6-layer neural network, 8-layer neural network, 6-layer BEACONS, and 8-layer BEACONS solutions against the formally-verified numerical solution, both for the final predicted frame, and across all frames.

Architecture	Conservation Error (Final)	Conservation Error (Total)
6-layer NN	69.348586	2254.019031
8-layer NN	30.707755	1442.490291
6-layer BEACONS	4.459920	-0.569301
8-layer BEACONS	12.788327	-35.540890

Table 5: Conservation analysis for the 1D inviscid Burgers’ “top-hat” initial value problem, comparing the 6-layer neural network, 8-layer neural network, 6-layer BEACONS, and 8-layer BEACONS solutions against the formally-verified numerical solution, both for the final predicted frame, and across all frames.

Panel of Figure 6 — Figure 6: Results for the 2D inviscid Burgers’ disk problem at times ${t=0.33}$ (left), ${t=0.66}$ (middle), and ${t=0.99}$ (right), obtained using a formally-verified numerical solver (top), an 8-layer neural network (middle), and an 8-layer BEACONS architecture (bottom). The 8-layer neural network mis-predicts the shape of the solution, causing it to appear much more “pointy” than it is, while also overestimating the propagation speed significantly, causing the advected quantity to hit the boundary of the domain before the end of the simulation. The BEACONS architecture successfully predicts both the correct shape and the correct propagation speed of the solution, matching the numerical solution more-or-less perfectly.

Architecture	${L^{\infty}}$ Error (Final)	${L^{2}}$ Error (Final)	${L^{\infty}}$ Error (All)	${L^{2}}$ Error (All)
8-layer NN	0.976685	37.415877	0.653738	17.396335
8-layer BEACONS	0.319521	2.688830	0.595719	1.888049

Table 6:

{L^{2}}

and

{L^{\infty}}

error analysis for the 2D inviscid Burgers’ disk problem, comparing the 8-layer neural network and 8-layer BEACONS solutions against the formally-verified numerical solution, both for the final predicted frame, and across all predicted frames.

Panel of Figure 7 — Figure 7: Results for the 1D compressible Euler Sod-type shock tube problem at times ${t=0.1}$ (left) and ${t=0.2}$ (right), obtained using a high-resolution numerical solver (blue), a 6-layer neural network (red), an 8-layer neural network (yellow), a 6-layer BEACONS architecture (purple), and an 8-layer BEACONS architecture (green). Both the 6-layer and 8-layer neural networks completely diffuse the qualitative structure of the contact discontinuity and the right-moving shock, and significantly underestimate the speed of the left-moving rarefaction. Both BEACONS solutions successfully predict the speed of the left-moving rarefaction; in the 6-layer BEACONS solution again the qualitative structure of the contact discontinuity and the right-moving shock is diffused, but the 8-layer BEACONS solution tracks to numerical solution more-or-less perfectly.

Architecture	${L^{\infty}}$ Error (Final)	${L^{2}}$ Error (Final)	${L^{\infty}}$ Error (All)	${L^{2}}$ Error (All)
6-layer NN	0.230453	4.522916	0.444256	2.075029
8-layer NN	0.159231	3.002094	0.452337	1.546402
6-layer BEACONS	0.076894	1.004591	0.422681	0.758668
8-layer BEACONS	0.080181	0.411615	0.523623	0.328464

Table 7:

{L^{2}}

and

{L^{\infty}}

error analysis for the 1D compressible Euler Sod-type shock tube problem, comparing the 6-layer neural network, 8-layer neural network, 6-layer BEACONS, and 8-layer BEACONS solutions against the high-resolution numerical solution, both for the final predicted frame, and across all frames.

Architecture	Conservation Error (Final)	Conservation Error (Total)
6-layer NN	133.878340	6219.332405
8-layer NN	83.996340	5127.036345
6-layer BEACONS	3.343495	281.148077
8-layer BEACONS	-0.991722	-211.287987

Table 8: Conservation analysis for the 1D compressible Euler Sod-type shock tube problem, comparing the 6-layer neural network, 8-layer neural network, 6-layer BEACONS, and 8-layer BEACONS solutions against the high-resolution numerical solution, both for the final predicted frame, and across all frames.

Panel of Figure 8 — Figure 8: Results for the 2D compressible Euler quadrants problem at times ${t=0.27}$ (left), ${t=0.53}$ (middle), and ${t=0.8}$ (right), obtained using a high-resolution numerical solver (top), an 8-layer neural network (middle), and an 8-layer BEACONS architecture (bottom). The 8-layer neural network completely diffuses the qualitative structure of the solution and mis-predicts many of the relevant wave-speeds; the protruding “ray” and surrounding structure on the top right of the domain is completely absent, the propagating wave front on the bottom left of the domain is incorrectly predicted as being concave rather than convex, etc. The BEACONS architecture predicts both the correct overall shape and the correct propagation speeds of the constituent waves of the solution, matching the numerical solution more-or-less perfectly (albeit with an absence of certain fine structures and instabilities due to increased numerical diffusion).

Architecture	${L^{\infty}}$ Error (Final)	${L^{2}}$ Error (Final)	${L^{\infty}}$ Error (All)	${L^{2}}$ Error (All)
8-layer NN	0.546258	25.843028	0.536434	12.315916
8-layer BEACONS	0.198093	2.958570	0.310002	1.539846

Table 9:

{L^{2}}

and

{L^{\infty}}

error analysis for the 2D compressible Euler quadrants problem, comparing the 8-layer neural network and 8-layer BEACONS solutions against the high-resolution numerical solution, both for the final predicted frame, and across all predicted frames.

Neural PDE Solvers

Neural PDE solvers use trainable models to approximate solutions or solution operators for partial differential equations, often with physical constraints.

See all 261 papers

Rigorous Error Certification for Neural PDE Solvers: From Empirical Residuals to Solution Guarantees
Mar 19, 2026Amartya Mukherjee, Maxwell Fitzsimmons, David C. Del Rey Fernández +1Neural Network Generalization Neural PDE Solvers
Department of Applied Mathematics, University of Waterloo, Waterloo, ON, Canada
Error-Conditioned Neural Solvers
Jun 25, 2026Haina Jiang, Liam Wang, Peng-Chen Chen +4Neural PDE Solvers PDE Surrogate Modeling
University of Michigan · KAIST AI · Los Alamos National Laboratory
Breakeven complexity: A new perspective on neural partial differential equation solvers
May 14, 2026Yijing Zhang, Nicholas Roberts, Tanya Marwah +1Neural PDE Solvers Benchmark Design
University of Wisconsin–Madison · Google DeepMind

Mar 19, 2026J/K move · Enter open · S save

Rigorous Error Certification for Neural PDE Solvers: From Empirical Residuals to Solution Guarantees

Amartya Mukherjee, Maxwell Fitzsimmons, David C. Del Rey Fernández, Jun Liu

Department of Applied Mathematics, University of Waterloo, Waterloo, ON, Canada

Neural Network Generalization Neural PDE Solvers PDE Solving Physics-Informed ML Generalization Bounds

Open paper AskPDF arXiv

Uncertainty quantification for partial differential equations is traditionally grounded in discretization theory, where solution error is controlled via mesh/grid refinement. Physics-informed neural networks fundamentally depart from this paradigm: they approximate solutions by minimizing residual losses at collocation points, introducing new sources of error arising from optimization, sampling, representation, and overfitting. As a result, the generalization error in the solution space remains an open problem. Our main theoretical contribution establishes generalization bounds that connect residual control to solution-space error. We prove that when neural approximations lie in a compact subset of the solution space, vanishing residual error guarantees convergence to the true solution. We derive deterministic and probabilistic convergence results and provide certified generalization bounds translating residual, boundary, and initial errors into explicit solution error guarantees.

Abstract

Neural PDE Solvers

Explore similar work

Rigorous Error Certification for Neural PDE Solvers: From Empirical Residuals to Solution Guarantees

Error-Conditioned Neural Solvers

Breakeven complexity: A new perspective on neural partial differential equation solvers