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
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.
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∞ error (for the non-smooth parts of the solution) for a shallow BEACONS architecture for solving the 1D inviscid Burgers’ equation.
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∞ Error (Final)
L2 Error (Final)
L∞ Error (All)
L2 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: L2 and L∞ 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.
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∞ Error (Final)
L2 Error (Final)
L∞ Error (All)
L2 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: L2 and L∞ 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.
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∞ Error (Final)
L2 Error (Final)
L∞ Error (All)
L2 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: L2 and L∞ 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.
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∞ Error (Final)
L2 Error (Final)
L∞ Error (All)
L2 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: L2 and L∞ 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.
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∞ Error (Final)
L2 Error (Final)
L∞ Error (All)
L2 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: L2 and L∞ 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.
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∞ Error (Final)
L2 Error (Final)
L∞ Error (All)
L2 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: L2 and L∞ 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 use trainable models to approximate solutions or solution operators for partial differential equations, often with physical constraints.
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.
Amartya Mukherjee, Maxwell Fitzsimmons, David C. Del Rey Fernández +1
Department of Applied Mathematics, University of Waterloo, Waterloo, ON, Canada
Neural surrogate models offer fast approximate mappings from PDE parameters to solutions, but they typically treat solving as a purely statistical task: once trained, they struggle to correct their own constraint violations and extrapolate beyond the training distribution. Recent hybrid methods promote physical correctness by targeting the PDE residual via gradient descent or Gauss--Newton steps, but inherit the compute cost and instability of the underlying classical optimizers. We show, theoretically and empirically, that numerically minimizing the PDE residual can be an unreliable proxy for reconstruction accuracy in ill-conditioned systems, explaining why these methods often do not make accurate predictions despite achieving low residuals. We propose error-conditioned Neural Solvers (ENS), built on a different principle: rather than an optimization target, the PDE residual field is passed as a direct input to the network at each iteration, enabling it to read the spatial structure of its own errors and learn an update policy to iteratively correct its predictions. Across four PDE families, ENS attains the highest prediction accuracy in the large majority of settings, with gains reaching 10× on turbulent Kolmogorov flow, while avoiding the expensive compute cost of hybrid methods. ENS's learned correction policy generalizes under distribution shift, including zero-shot parameter changes and cross-equation transfer, where its relative advantage is largest in the ill-conditioned regimes where residual minimization is least reliable. Project website: https://neuralsolver.github.io/.
Haina Jiang, Liam Wang, Peng-Chen Chen +4
University of Michigan · KAIST AI · Los Alamos National Laboratory
Neural surrogate solvers of partial differential equations (PDEs) promise dramatic speedups over numerical methods, especially in scenarios requiring many solves. However, current accuracy-based evaluations do not fully consider two central issues: (1) neural solvers incur substantial up-front costs for data generation, training, and tuning; and (2) classical solvers can also generate low-fidelity solutions at a sufficiently low simulation cost. To explicitly account for these realities and fully incorporate end-to-end costs, we propose an evaluation framework centered on breakeven complexity, a metric that counts the forward solves before a learned solver is cost-effective relative to an error-equivalent traditional solver. To evaluate this measure, we apply scaling laws to determine how much training budget to allocate to data generation and discuss how to achieve smooth error-matching in diverse settings. We evaluate the breakeven complexity of multiple neural PDE solvers on three PDEs on 2D periodic domains from APEBench and a novel benchmark of flows past multiple obstacles generated by the GPU-native PyFR code. Among other findings, our results suggest that neural PDE solvers become more effective as problems get harder in terms of cost, dimension, rollout, physics regime (e.g. higher Reynolds number), etc.