Partial Differential Equations in the Age of Machine Learning: A Critical Synthesis of Classical, Machine Learning, and Hybrid Methods

Mohammad Nooraiepour; Jakub Wiktor Both; Teeratorn Kadeethum; Saeid Sadeghnejad

cs.LGMar 8, 2026

Partial Differential Equations in the Age of Machine Learning: A Critical Synthesis of Classical, Machine Learning, and Hybrid Methods

Authors: Mohammad Nooraiepour, Jakub Wiktor Both, Teeratorn Kadeethum, Saeid Sadeghnejad

Organizations: Faculty of Mathematics and Natural Sciences, University of Oslo, 1047 Blindern, 0316 Oslo, Norway · Department of Mathematics, University of Bergen, Bergen, Norway · AI Lab, Siemens Energy, Orlando, FL 32826, USA · Institute for Geosciences, Applied Geology, Friedrich-Schiller-University Jena, 07749 Jena, Germany

Abstract

Partial differential equations (PDEs) govern physical phenomena across the full range of scientific scales, yet their computational solution remains one of the defining challenges of modern science. This critical review examines two mature but epistemologically distinct paradigms for PDE solution, classical numerical methods and machine learning approaches, through a unified evaluative framework organized around six fundamental computational challenges. Classical methods are assessed for their structure-preserving properties, rigorous convergence theory, and scalable solver design; their persistent limitations in high-dimensional and geometrically complex settings are characterized precisely. Machine learning approaches are introduced under a taxonomy organized by the degree to which physical knowledge is incorporated and subjected to the same critical evaluation applied to classical methods. Classical methods are deductive -- errors are bounded by quantities derivable from PDE structure and discretization parameters -- while machine learning methods are inductive -- accuracy depends on statistical proximity to the training distribution. This epistemological distinction is the primary criterion governing responsible method selection. We identify three genuine complementarities between the paradigms and develop principles for hybrid design, including a framework for the structure inheritance problem that addresses when classical guarantees propagate through hybrid couplings, and an error budget decomposition that separates discretization, neural approximation, and coupling contributions. We further assess emerging frontiers, including foundation models, differentiable programming, quantum algorithms, and exascale co-design, evaluating each against the structural constraints that determine whether current barriers are fundamental or contingent on engineering progress.

Figures & tables

Figure 1: Hierarchical taxonomy of partial differential equations (PDEs), organized by type (elliptic, parabolic, hyperbolic), linearity, and canonical examples with their governing mathematical forms.

Figure 2: Characterisation of the six fundamental computational challenges in PDE solution. A. Radar chart showing the challenge severity profile (scored 0–10) for six representative application domains across the axes: (I) high dimensionality, (II) nonlinearity and multiplicity, (III) geometric complexity, (IV) discontinuities and non-smoothness, (V) multiscale phenomena, and (VI) multiphysics coupling. B. Symmetric interaction-severity matrix quantifying how each pair of challenges amplifies computational difficulty. The diagonal entries denote inherent challenge severity, and the markers indicate the dominant co-occurring challenge pair for each application.

Figure 3: Illustration of the multiscale challenge for PDEs in spatially heterogeneous media with rapidly oscillating coefficient

a(x,x/\epsilon)

. Shown across four scale-separation ratios

\epsilon\in\{0.20,0.10,0.04,0.015\}

: (A) Coefficient

a(x,x/\epsilon)=1.5+0.8\sin(2\pi x/\epsilon)

for decreasing

\epsilon

. The dashed line indicates the homogenized value

a^{*}=1.5

. Inset: magnification near

x=0

showing increasingly rapid oscillations that direct discretization must resolve. (B) Solutions

u^{\epsilon}(x)

(solid) compared to the smooth homogenized limit

u^{*}(x)

(dashed). Inset: zoom at the peak (

x\approx 0.5

) revealing that oscillatory corrections decay with

\epsilon

but remain visible on the microscale, even where the macroscopic profiles are nearly indistinguishable. (C) Relative

L^{2}

error versus mesh size

h

(decreasing left to right) for direct finite element discretization and the homogenized solver. Direct simulation shows an

\mathcal{O}(\epsilon)

error plateau for

h\gg\epsilon

, with second-order

\mathcal{O}(h^{2})

convergence only when

h\ll\epsilon

— a resolution demand that becomes computationally intractable for small

\epsilon

. The homogenized approach yields uniform

\mathcal{O}(h^{2})

convergence independent of

\epsilon

at greatly reduced cost.

Figure 4: Illustration of four principal discretization paradigms on a common L-shaped domain, highlighting differences in geometric flexibility and approximation properties. (A) Finite difference method (FDM): structured Cartesian grid with the canonical five-point stencil emphasized. It shows the central-difference approximation

\partial^{2}u/\partial x^{2}\approx(u_{i+1}-2u_{i}+u_{i-1})/h^{2}

at the stencil center. (B) Finite element method (FEM): unstructured triangulation with one representative element

K_{e}

highlighted. The global system is governed by the weak form

a(u_{h},v_{h})=\ell(v_{h})

for all test functions

v_{h}\in V_{h}

. (C) Finite volume method (FVM): structured control-volume partition enforcing local conservation. The integral form

\frac{d}{dt}\int_{V_{i}}u\,dV+\oint_{\partial V_{i}}\mathbf{F}\cdot\mathbf{n}\,dS=0

is applied cell-wise; flux arrows indicate the anti-symmetric condition

F_{ij}=-F_{ji}

. (D) Spectral element method (SEM): three macro-elements

\Omega_{1}

–

\Omega_{3}

with Gauss–Lobatto–Legendre (GLL) quadrature points displayed. For smooth solutions, exponential convergence

\|u-u_{N}\|\leq Ce^{-\sigma N}

holds within each element.

Figure 5: Quantitative illustration of the curse of dimensionality for classical grid-based and sampling methods applied to

d

-dimensional PDEs and parametric problems. A. Degrees of freedom (or function evaluations) required for a fixed resolution

N

as a function of dimension

d

, plotted on a logarithmic scale. Full tensor-grid discretisation with

N

points per dimension scales as

N^{d}

(crimson curves;

N=20

dashed,

N=50

solid), crossing the practical ceiling of

10^{12}

degrees of freedom (shaded region, corresponding to terabyte-scale memory) at

d=6

for

N=50

and

d=8

for

N=20

. The Smolyak sparse grid [ 329 ] scales as

\mathcal{O}(N(\log_{2}N)^{d-1})

(blue), providing order-of-magnitude reductions over the full grid but requiring bounded mixed derivatives

\|D^{\alpha}f\|_{L^{2}}\leq C_{|\alpha|_{1}}

; the dashed continuation past

d=14

marks the regime where this smoothness condition generically fails. Monte Carlo (MC, emerald) and quasi-Monte Carlo (QMC, teal) methods are dimension-independent by construction, with sample counts determined solely by the target accuracy

\varepsilon

and independent of

d

; the reduced basis method (gold, dashed) achieves low online cost

\mathcal{O}(N_{\mathrm{rb}}^{3})

when the solution manifold

\mathcal{M}_{\mu}

has low Kolmogorov

n

-width, but its offline snapshot cost grows with the dimension of the parameter space. Background shading identifies four tractability regimes: all methods viable (

d\leq 3

); sparse grid viable (

d=3

–

8

); only MC/QMC viable (

d=8

–

14

); and MC viable at reduced accuracy only (

d>14

). B. Achievable error

\varepsilon

for a fixed total degree-of-freedom budget of

N_{\mathrm{total}}=10^{6}

, showing how accuracy degrades with dimension for each method. Full-grid FD/FEM with

k=2

and sparse-grid

k=2

deliver engineering accuracy (

\varepsilon\leq 10^{-2}

) only up to

d\approx 4

and

d\approx 6

, respectively, within this budget. Monte Carlo achieves a dimension-independent accuracy of

\varepsilon\sim N^{-1/2}=10^{-3}

regardless of

d

. Horizontal reference lines mark engineering (

\varepsilon=10^{-2}

) and scientific (

\varepsilon=10^{-4}

) accuracy thresholds. C. The formal DOF scaling expression for the four method families, the key mathematical assumption required for that scaling to hold, and the condition under which the method breaks down. The reduced basis method occupies a qualitatively different position: its tractability depends not on the spatial dimension

d

but on the intrinsic dimension of the solution manifold

\mathcal{M}_{\mu}

, quantified by the Kolmogorov

n

-width

d_{n}(\mathcal{M}_{\mu})\to 0

; when this condition fails (e.g. for transport-dominated or bifurcating problems), the method loses its efficiency advantage. Taken together, the three panels establish that breaking the curse of dimensionality within the classical paradigm requires either strong a priori regularity assumptions (sparse grid) or a fundamental change in approximation paradigm (MC/QMC, reduced basis), each of which incurs its own structural cost.

Figure 6: Capability assessment matrix for five classical discretization families across the six core computational challenges outlined in Section 2 . Each cell is color-coded based on the method’s inherent mathematical suitability, accompanied by a concise annotation of the key enabling mechanism or limiting factor. Rows: finite difference methods (FDM), finite element methods (FEM), finite volume methods (FVM), spectral/spectral element methods, and meshless/boundary element methods. Columns: (I) high dimensionality, (II) nonlinearity and solution multiplicity, (III) geometric complexity, (IV) discontinuities and non-smoothness, (V) multiscale phenomena, (VI) multiphysics coupling. No classical method excels uniformly across all challenges. The most intractable weaknesses—high dimensionality (I) and multiscale phenomena (V)—persist across the board and provide the primary impetus for exploring machine learning approaches in Section 5 .

Figure 7: Classification of ML models for PDE modeling. Methods are organized by the degree to which physical knowledge is embedded into the learning process, ranging from purely data-driven surrogates to hybrid approaches that tightly couple learned components with classical solvers.

Method Family	Core Principle	Complexity	Data Needs
Pure Data-Driven Surrogates
Generic Deep Surrogates (CNNs, MLPs, U-Nets, RNNs)	Black-box regression from inputs to solution fields	$O(N)$ – $O(N^{2})$	High
Graph Neural Networks	Message passing on meshes and point clouds	$O(\|V\|+\|E\|)$	Moderate
Transformers	Global attention over spatial tokens	$O(N^{2})$ ; $O(N\log N)$ with efficient attention	High
Generative Models	Probabilistic distribution over solutions	Problem-dependent	Moderate
Physics-Embedded Methods

Table 1: Machine learning method families for PDE solving: core principles, computational complexity, and data scaling. Complexity reflects the typical computational cost per forward evaluation; actual performance depends on problem-specific factors including resolution, dimensionality, and architecture choices.

Variant	Key Modification	Use Case	Reference
Variational (VPINN)	Weak formulation via energy functional	Problems with variational structure	[ 197 , 313 ]
Extended (XPINN)	Domain decomposition, interface conditions	Complex geometries, parallelization	[ 178 , 168 ]
Conservative	Architectural conservation enforcement	Long-time integration, transport	[ 179 ]
Deep Ritz	Energy minimization for elliptic PDEs	Variational problems	[ 184 , 410 ]
Weak Adversarial	Min-max formulation avoiding high derivatives	High-order PDEs	[ 419 , 283 ]
Multi-fidelity	Hierarchical data fusion	Mixed-quality data sources	[ 262 , 363 ]

Table 2: PINN variants: modifications and applications.

Figure 8: Spectrum of hybrid ML–PDE methods, ordered by retained physical fidelity (left: fully classical solvers; right: fully data-driven surrogates). Five representative strategy families are positioned along the axis: - ML-accelerated solvers : preserve the full PDE operator, augmenting only numerical kernels with learned components for speed; - residual / error-correcting models : maintain physics baseline while learning discrepancy terms for unresolved scales or missing closures; - hybrid physics-ML solvers : replace costly sub-models (e.g., constitutive laws, fine-scale homogenization, geochemical modules) inside a structure-preserving PDE framework; - physics-enhanced deep surrogates & neural differential equations (PEDS, Neural ODEs): leverage low-fidelity physics as inductive bias, expressing dynamics via adjoint-friendly continuous architectures; - Koopman and foundation models : linearize nonlinear dynamics in learned latent spaces or enable transfer across PDE families with minimal retraining. The upper boxes list provably preserved physical structures (conservation laws, symmetries, boundary conditions); the lower boxes highlight key computational gains over classical or purely data-driven baselines.

Application Scenario	Recommended Approach	Rationale
Unknown or complex material constitutive behavior	Neural constitutive models	Preserves equilibrium; learns material response from data while satisfying thermodynamic constraints
Repeated microscale homogenization at many macroscale points	Homogenization surrogates (FE 2 replacement)	Amortizes expensive offline microscale training across online macroscale queries
Scarce training data with known governing equations	PEDS (physics-enhanced deep surrogate)	Low-fidelity physics solver provides strong inductive bias, reducing labeled data requirements by 10–100 $\times$
High accuracy required for smooth solutions	DeePoly (neural + polynomial refinement)	Achieves near-spectral convergence through polynomial correction; applicable to smooth problems only
Parameter estimation or design optimization from observations	Differentiable physics	Enables gradient-based optimization through the solver via adjoint or implicit differentiation
Highly heterogeneous or high-contrast coefficients	Neural preconditioners	Learns approximate inverse operators that reduce Krylov iteration count in regimes where algebraic multigrid may degrade

Table 3: Hybrid and physics-coupled method selection guidelines. The entries in this table are not mutually exclusive, as many applications benefit from combining multiple strategies.

Figure 9: Computational economics of machine learning methods for PDE solution. All cost and speedup estimates are order-of-magnitude syntheses from representative published benchmarks [ 224 , 244 , 192 , 294 , 350 , 76 , 160 ] and illustrate the relative computational character of each method family. Actual values vary significantly with problem dimensionality, resolution, rollout horizon, implementation quality, and classical baseline. (A) Log-log scatter of training cost (GPU-hours) versus inference speedup (relative to classical solver) for eight representative ML-PDE families. Diagonal iso-lines show break-even query count

N_{\rm be}=C_{\rm train}/(C_{\rm classical}-C_{\rm inference})

(Eq. 50 ): methods above the

N_{\rm be}=100

line recoup training investment after fewer than 100 queries. Shaded quadrants separate the favorable regime (high speedup, low training cost) from the unfavorable one (low speedup, high training cost). (B) Cumulative cost versus query count for a classical solver and three representative neural methods, based on the amortization model of Eq. 50 . Vertical dashed lines indicate break-even crossovers: the minimum query count at which each neural method becomes cheaper overall than repeated classical solves. Shaded regions delineate unfavorable regimes (few-query problems, left) from favorable ones (many-query scenarios such as uncertainty quantification, parametric studies, and real-time control, right).

Method	GPU-Hours	Training Data	Typical GPU Memory Scaling
Standard PINN [ 192 ]	$10^{1}$ – $10^{2}$	Minimal: PDE + BC/IC; optional $O(10^{2})$ sparse observations	$O(B\cdot P)$ for activations; $O(N_{\text{pts}})$ for collocation storage
FNO [ 224 ]	$10^{2}$ – $10^{3}$	$10^{3}$ – $10^{4}$ solution fields or rollouts (supervised)	$O(B\cdot C\cdot N^{d})$ for grid-based operators
DeepONet [ 244 ]	$10^{2}$ – $10^{3}$	$10^{3}$ – $10^{5}$ input–output function pairs	$O(B\cdot(N_{s}+N_{q})\cdot C)$
GNN simulators [ 294 ]	$10^{2}$ – $10^{4}$	$10^{2}$ – $10^{4}$ trajectories (mesh or graph)	$O(B\cdot(\|V\|+\|E\|)\cdot C)$
Diffusion / score-based [ 350 ]	$10^{3}$ – $10^{4}$	$10^{3}$ – $10^{5}$ samples (conditional on parameters or BCs)	$O(B\cdot N\cdot C)$ ; inference cost scales with $T_{\text{steps}}$
Neural ODEs [ 76 ]	$10^{1}$ – $10^{3}$	$10^{2}$ – $10^{4}$ trajectories	$O(B\cdot P)$ via adjoint; independent of integration depth

Table 4: Order-of-magnitude training costs by method family. Compute and memory are indicative; actual costs vary substantially with problem size, rollout horizon, and hardware. GPU memory reflects typical on-device usage per batch during training, not offline dataset storage.

Figure 10: Hierarchical summary of major challenges in applying machine learning to scientific computing, organized into four major categories: solution reliability, data-related challenges, computational challenges, and model and theoretical limitations.

Figure 11: Out-of-distribution (OoD) generalisation and degradation for representative ML-PDE method families. All contours and curves reflect qualitative trends synthesised from published benchmarks and illustrate characteristic degradation behaviour rather than precise, problem-specific performance. Actual rates depend on the PDE type, the parameter-space geometry, the training-set coverage, and the architectural choices. (A) Two-dimensional slice of parameter space (e.g., Reynolds number vs. diffusion coefficient). The white ellipse bounds the training distribution

\mathcal{D}_{\mathrm{train}}

; the background colour field shows normalised relative

L^{2}

error (green: near-zero inside training support; crimson: high far outside). Iso-accuracy contours indicate extrapolation range before unacceptable degradation: - PINNs (steel, solid): minimal extension beyond training boundary (per-instance training). - Neural operators FNO/DeepONet (emerald, dashed): moderate extension via resolution-invariant operator learning. - Hybrid physics-ML methods (gold, dash-dot): substantially larger reach, driven by PDE operator as inductive bias. - Foundation models (crimson, dotted): irregular/asymmetric boundary due to uneven cross-PDE-family transfer. Classical solvers (omitted) show no distribution-shift dependence; error is controlled solely by discretisation parameters. (B) Normalised

L^{2}

error versus OoD distance

d(\xi,\mathcal{D}_{\mathrm{train}})

for the four ML families plus classical reference (dark navy, flat). Vertical dashed line: training boundary. Shaded regions separate in-distribution from OoD regimes. Horizontal band: representative acceptable-accuracy threshold. PINNs degrade most steeply; hybrids degrade most slowly (physics structure persists); foundation models show intermediate rates with slower initial fall-off (benefit of broad pretraining) but higher in-distribution error. No ML family achieves the distribution-independent behaviour of classical solvers.

Panel of Figure 12 — Figure 12: Bayesian optimization for hyperparameter selection. Left: After evaluating validation performance (red curve, unknown to the optimizer) at a sparse set of hyperparameter configurations, a Gaussian process surrogate (black line, with $\pm 1\sigma$ uncertainty in gray) models the objective landscape. The expected improvement acquisition function (green) identifies the next configuration to evaluate by balancing exploitation of the current best region (left peak) against exploration of high-uncertainty regions (right peak). In this iteration, exploration dominates and the next candidate is selected near $\lambda=8.5$ . Right: After observing the objective at the new configuration, the surrogate model is updated with substantially reduced uncertainty in the explored region. The acquisition function now concentrates around the promising region near $\lambda=2.7$ , demonstrating convergence toward the true optimum. For neural PDE solvers, $\lambda$ represents a high-dimensional hyperparameter vector and each evaluation requires training a complete model to convergence, making systematic optimization computationally expensive. Adapted from [ 125 ] .

Figure 13: Tractability assessment of persistent challenges facing ML-based PDE solvers, organised by the four categories of Sections 6.1 – 6.4 . Each row details one failure mode across three columns: Current Status (how the challenge manifests in practice today), Improvement Trajectory (most credible near- to medium-term mitigation pathways from the literature), and Tractability ( Addressable , Partially addressable , or Fundamental ), reflecting whether clear technical solutions exist or the limitation is intrinsic to neural approximation and stochastic optimisation. It also indicates tractability level: green (addressable, solutions already demonstrated), gold (partially addressable, requiring further conceptual progress), red (fundamental, unlikely to yield without new mathematical foundations). The rightmost column uses filled coloured markers (following the taxonomy of Figure 7 ) to show which method families are primarily affected; unfilled circles indicate minor or no impact. No Fundamental challenge is confined to a single architecture: the accuracy ceiling and lack of a priori error bounds affect every ML-PDE family reviewed, confirming these as intrinsic properties of the neural paradigm rather than specific design choices.

Challenge	Current Status	Improvement Trajectory	Tractability
Physical inconsistency	Frequent in pure data-driven models	Hard architectural constraints; energy-based models	Addressable
Training instability	Significant for PINNs	NTK-based weighting; curriculum learning	Addressable
Hyperparameter burden	Pervasive	AutoML; principled initialization	Partially addressable
Multiscale handling	Architecture-dependent	Wavelet operators; hierarchical GNNs	Partially addressable
Data generation cost	Dominates for supervised methods	Multi-fidelity; active learning	Partially addressable
OoD generalization	Rapid degradation	Foundation models; physics-informed training	Partially addressable

Table 5: Assessment of ML-PDE challenge tractability. Challenges are classified by whether improvement trajectories are clear (Addressable), partially available (Partially addressable), or reflect fundamental properties of neural approximation and stochastic optimization (Fundamental).

Method	Architecture	Complexity	UQ	Accuracy †	Best Application
Physics-Informed
Standard PINN	Residual MLP	$O(N_{\text{pts}}\cdot P)$	Limited	$10^{-3}$ – $10^{-1}$	Inverse problems, sparse data
VPINN	Energy functional	$O(N_{\text{quad}}\cdot P)$	Limited	$10^{-4}$ – $10^{-1}$	Variational problems
XPINN	Domain decomp.	$O(K\cdot N\cdot P)$	No	$10^{-3}$ – $10^{-1}$	Complex geometries
Neural Operators
FNO	Fourier layers	$O(N\log N)$	Approx. (ensemble)	$10^{-2}$ – $10^{-1}$	Periodic, parametric

Table 6: ML-PDE method comparison. Accuracy † reflects indicative relative

L^{2}

errors; values vary substantially across problem settings. UQ column indicates native support (Native), approximate support via ensemble or dropout (Approx.), or absence (No/Limited).

Axis	Classical Numerical Methods	Machine Learning Methods
Knowledge encoding	Structural: conservation laws, variational principles, and inf-sup conditions encoded by construction into discretization and solver design	Statistical: prior knowledge enters via training distribution, architecture inductive bias, and loss function regularization
Error characterization	Deductive: a priori and a posteriori bounds derived from PDE regularity and discretization parameters; convergence hierarchy is systematic and refinable	Empirical: accuracy assessed on held-out data; no refinement hierarchy; training loss does not certify deployment accuracy
Generalization guarantee	Convergence for any well-posed input satisfying regularity assumptions, independent of training history	Statistical: degrades for inputs outside training distribution; rate of degradation generally unpredictable without additional assumptions
Failure mode character	Identifiable and quantifiable: spurious oscillations (Gibbs), pressure locking, non-physical modes, CFL violation—each with known cause and remedy	Often opaque: physically inadmissible outputs, silent distributional failure, accumulated rollout errors not detectable from internal model quantities
Human-in-the-loop form	Geometric and mathematical: mesh generation, solver parameter selection, stability analysis, preconditioner design	Statistical and engineering: architecture search, hyperparameter tuning, training pipeline design, validation dataset construction

Table 7: Comparative anatomy of classical numerical methods and machine learning approaches for PDE solution across five discriminating axes. Entries characterize paradigm-level behavior, not every individual method within each paradigm.

Figure 14: Convergent evolution of classical numerical methods and machine learning approaches toward shared mathematical structures for PDE approximation. The figure depicts classical solvers (left) and ML methods (right), both rooted in three common mathematical foundations: operator/functional analysis, variational principles and Sobolev theory, and spectral/harmonic analysis. Four horizontal levels connect independently developed method pairs. Each level is colour-coded by its shared mathematical domain. Central convergence boxes comprise three tiers: the structural parallel, the precise governing equations/mathematical equivalence, and the key distinction that separates the approaches despite their shared architecture. Row 1: Geometric multigrid

\leftrightarrow

graph neural networks (GNNs). Both propagate corrections hierarchically across scales. A multigrid V-cycle is mathematically equivalent to a GNN with fixed aggregation weights

W_{s}

W_{e}

determined by the operator

A

rather than learned. Classical multigrid guarantees convergence via two-grid analysis (

h

-ellipticity, approximation property); GNNs learn the hierarchy from data without certified rates. Row 2: Spectral Galerkin

\leftrightarrow

Fourier Neural Operator (FNO). Both rely on truncated spectral expansions computed in

\mathcal{O}(N\log N)

via FFT, yielding identical Gibbs ringing at discontinuities. Classical methods certify convergence via Jackson/Bernstein theory; FNO’s learned filter

R_{\theta}(k)

lacks prescribed approximation rates. Row 3: Reduced Basis / POD-Galerkin

\leftrightarrow

neural operators (DeepONet, FNO). Both exploit low Kolmogorov

n

-width of parametric solution manifolds. DeepONet’s branch-trunk decomposition nonlinearly generalises the POD ansatz (trunk

\approx

modes, branch

\approx

coefficients), but learned rather than SVD-derived. This enables approximation in regimes requiring many POD modes, at the expense of the greedy certificate and a priori dimension bound from

n

-width theory. Row 4: Sparse grids / Smolyak construction

\leftrightarrow

transformer self-attention. Both capitalise on energy concentration in low-order ANOVA terms of high-dimensional functions. Sparse grids analytically derive interactions from mixed Sobolev smoothness

\|D^{\alpha}f\|_{L^{2}}\leq C_{|\alpha|_{1}}

, reducing cost from

\mathcal{O}(N^{d})

\mathcal{O}(N(\log N)^{d-1})

with certification. Transformer attention learns informative interactions from data, forgoing a priori regularity assumptions at the cost of

\mathcal{O}(n^{2}d)

scaling and no certified dimension reduction. Convergent evolution across these parallels is driven by a single mathematical imperative: efficient PDE approximation—classical or learned—must exploit the low effective dimensionality, hierarchical organisation, and spectral content of solution manifolds.

Figure 15: Complementarity phase space for classical numerical methods and machine learning approaches to PDE solution. A. Two-dimensional deployment space spanned by the degree of physics specification (horizontal axis, from empirically characterised to be fully specified from first principles) and the certification requirement (vertical axis, from statistical accuracy sufficient to deductive certified accuracy required). Three shaded regions demarcate the natural domain of each paradigm: the classical-dominant zone where the governing operator

\mathcal{L}

is fully known and certification is required; the ML-dominant zone where physics is unknown or intractable and statistical accuracy suffices; and the hybrid transition where the problem decomposes naturally into a classical skeleton and a learned component. Two dashed boundary curves mark the certifiability threshold

\mathcal{C}

(above which deductive error bounds are required) and the physics specification threshold

\mathcal{P}

(to the right of which the governing equations are sufficiently characterised for classical discretisation); the hybrid zone is the region between these boundaries. Eighteen representative application domains are plotted as labelled points, colour-coded by discipline cluster: safety and certification, rigorous physical modelling, data-rich and empirical physics, inverse and data assimilation problems, surrogate and real-time deployment, and the research frontier methods. The deployment decisions depend simultaneously on both axes. B. The three genuine complementarities identified in Section 7.3 , each corresponding to a zone of panel A. Each box presents the intrinsic limitation of each paradigm and the open target — the capability that a well-designed hybrid method must achieve to address both limitations simultaneously. Complementarity 1: classical refinement guarantees convergence but costs

\mathcal{O}(h^{-d})

in dimension

d

; ML avoids the curse of dimensionality but provides no certified bound outside the training distribution

\mathcal{D}_{\mathrm{train}}

. Complementarity 2: classical methods enforce conservation and inf-sup conditions by construction, but require geometry-conforming meshes; ML methods are mesh-free but trade geometric error guarantees for distributional assumptions on the training geometry. The open target for Complementarity 2 — simultaneously mesh-free and conservation-certified, enforcing

\nabla\cdot\mathbf{F}=0

\Delta S\geq 0

by architecture on unstructured inputs — is a core unsolved problem of hybrid design. Complementarity 3: classical methods certify solutions to completely and correctly specified problems; ML methods accommodate unknown constitutive relations and empirical closures, but at the cost of conflating physical uncertainty and statistical approximation error with no mechanism for separating the two. These limitations are intrinsic in the sense of Section 7.3 : they follow from the foundational assumptions of each paradigm and cannot be resolved by engineering refinement within that paradigm alone.

Figure 16: Error budget decomposition across hybrid design archetypes (Eq. ( 51 ): terms (1) discretisation, (2) neural approximation, (3) coupling). A. Illustrative relative magnitudes of all three terms across five archetypes, ordered conservative to aggressive. A1 Neural constitutive law : FEM skeleton intact; classical

h

p

certificates preserved; coupling small. A2 Neural preconditioner : coupling structurally near-zero, as the preconditioner affects convergence rate only, not the converged solution. A3 Physics-ML hybrid : all three terms contribute; closest to a balanced budget. A4 Learned dynamics (Neural ODE / PEDS): coupling accumulates per autoregressive step and dominates alongside neural approximation error. A5 End-to-end neural solver (e.g. FNO, DeepONet, PINN): no classical skeleton; coupling is identically zero; effectively all error is neural approximation with no computable bound. Dashed line marks the A3 total budget

\varepsilon

used as normalisation in panel B. B. Balanced error budget for A3. Bars show the fraction of

\varepsilon

per term; dashed line at

\varepsilon/3

is the equal-investment target; arrows indicate over- or under-targeting. Control mechanisms (right): mesh refinement at certified rate

{\sim}\,h^{k}

for (1); architecture, data

N

, regularisation

\lambda

with no certified rate for (2); Schwarz / mortar / DD coupling and flux continuity for (3). Term (3) is chronically under-targeted because no computable a posteriori estimator for

\|\mathcal{N}_{\theta}(u_{h})-u_{\mathrm{hybrid}}\|_{V}

yet exists (cf. open problem P4, § 7.5 ).

Method	Computational	Mesh	UQ	Nonlinearity	Theory	Key Distinguishing Characteristics
	Cost Class	Adaptivity	Support	Handling	Guarantees
FINITE DIFFERENCE METHODS
Standard FDM	Low	None	No	Explicit time	Partial	Optimal memory efficiency; structured grids only; von Neumann stability analysis; CFL restrictions $\Delta t\sim h$ (hyperbolic), $h^{2}$ (parabolic); 2nd-order typical
Compact FD	Low-Medium	None	No	Implicit stencils	Strong	Padé-type schemes; tridiagonal systems; spectral-like resolution; 4th-6th order standard; excellent for DNS/LES
High-Order Compact	Medium	None	No	Implicit stencils	Strong	8th-10th order achievable; requires high solution regularity; recent advances in dispersion-relation-preserving schemes
WENO	Medium	r-adapt	No	Conservative	Strong	Nonlinear adaptive stencils; TVD/ENO properties; robust shock capturing; 5th order standard; improved variants address critical points

Table 8: Comprehensive comparison of classical PDE solution methods: computational characteristics, adaptivity capabilities, nonlinearity handling, theoretical guarantees, implementation requirements, and key distinguishing features.

Appendix figures & tables2 assets

Supplementary material from the paper’s appendix.

Appendix

Method	Adaptation Cost	Expressiveness	Uncertainty	Best For
MAML	Medium	High	No	Maximum quality
Reptile	Medium	High	No	Simple implementation
Neural Processes	Low	Medium	Yes	Sparse observations
Hypernetworks	Low	Medium	No	Real-time adaptation
Transfer Learning	Medium	High	No	Hierarchical tasks

Appendix

Table S.1 : Meta-learning method comparison for PDE applications.

Method	Architecture	Complexity	UQ	Multiscale	Nonlin.	Inference	Accuracy	Notes
Physics-Informed Neural Networks
Standard PINN	Residual MLP	$O(N_{\text{pts}}\cdot P)$	Limited	Moderate	Sensitive	Fast	$10^{-3}$ – $10^{-1}$	AD-based; spectral bias
VPINN	Energy functional	$O(N_{\text{quad}}\cdot P)$	Limited	Good	Better	Fast	$10^{-4}$ – $10^{-1}$	Superior conditioning
XPINN	Domain decomp.	$O(K\cdot N\cdot P)$	No	Good	Moderate	Medium	$10^{-3}$ – $10^{-1}$	Parallel subdomains
Conservative	Conservation arch.	$O(N_{\text{pts}}\cdot P)$	No	Moderate	Physical	Fast	$10^{-3}$ – $10^{-1}$	Long-time stable
Deep Ritz	Variational	$O(N_{\text{quad}}\cdot P)$	No	Limited	Smooth	Fast	$10^{-4}$ – $10^{-1}$	High-dimensional

Appendix

Table S.2: Comprehensive ML-PDE method comparison. Accuracy values are indicative relative

L^{2}

errors from PDEBench [ 366 ] and benchmark studies (2020–2025); results vary substantially with PDE family, coefficient regime (in-distribution vs. out-of-distribution), rollout horizon, resolution, and training protocol.

PDE Solving

PDE solving develops numerical and learned methods to approximate solutions of partial differential equations across scientific applications.

See all 249 papers

Sep 29, 2025

A Greedy PDE Router for Blending Neural Operators and Classical Methods

When solving PDEs, classical numerical solvers are often computationally expensive, while machine learning methods can suffer from spectral bias, failing to capture high-frequency components. Designing an optimal hybrid iterative solver--where, at each iteration, a solver is selected from an ensemble of solvers to leverage their complementary strengths--poses a challenging combinatorial problem. While greedy selection is desirable for its constant-factor approximation guarantee to the optimal solution under Lipschitz assumptions, it requires knowledge of the true error at each step, which is unavailable in practice. We address this by proposing an approximate greedy router that efficiently mimics a greedy approach to solver selection. Empirical results on the Poisson, convection-diffusion, and anisotropic diffusion equations show that our method consistently reduces time-to-solution relative to single-solver baselines and existing hybrid approaches such as HINTS.

Sahana Rayan, Yash Patel, Ambuj Tewari

Department of Statistics University of Michigan Ann Arbor, MI 48104, USA

PDE Solving

Ask PDF

Jun 18, 2026

Learning universal approximations for partial differential equations with Physics-Informed Broad Learning System

Partial differential equations (PDEs) play a central role in modeling complex physical, biological, and engineering systems. While traditional numerical solvers are robust, they often incur prohibitive computational costs due to mesh dependencies, whereas recent Physics-Informed Neural Networks (PINNs) offer a mesh-free alternative but frequently suffer from slow convergence and optimization instability. To bridge this gap, this article proposes the Physics-Informed Broad Learning System (PIBLS), a novel backpropagation-free framework that reformulates PDE solving as a direct least-squares optimization. We improved an algorithm within this framework to handle nonlinear PDEs efficiently and provide a rigorous mathematical proof establishing the universal approximation property of PIBLS for these equations. Experiments on linear and nonlinear PDEs demonstrate that PIBLS is one to three orders of magnitude faster than conventional PINNs while achieving significantly higher solution accuracy. This framework provides a computationally efficient paradigm for scientific machine learning, offering a practical, high-speed alternative for real-time simulation and design optimization tasks.

Zhiwen Yu, Derong Yang, Liujian Zhang +5

School of Computer Science and Engineering, South China University of Technology, Guangzhou 510006, China · Peng Cheng Laboratory, Shenzhen 518000, China · School of Future Technology, South China University of Technology, Guangzhou 510006, China +2

PDE Solving Universal Approximation Neural Network Approximation Theory Physics-Informed ML

Ask PDF

Apr 2, 2026

Flow Learners for PDEs: Toward a Physics-to-Physics Paradigm for Scientific Computing

Partial differential equations (PDEs) govern nearly every physical process in science and engineering, but solving them at scale remains prohibitively expensive. Generative AI has transformed language, vision, and protein science, but learned PDE solvers have not undergone a comparable shift. Existing paradigms each capture part of the problem. Physics-informed neural networks embed residual structure, although they are often difficult to optimize in stiff, multiscale, or large-domain regimes. Neural operators amortize across instances, although they commonly inherit a snapshot-prediction view of solving and can degrade over long rollouts. Diffusion-based solvers model uncertainty, although they are often built on a solver template that still centers on state regression. We argue that the core issue is the abstraction used to train learned solvers. Many models are asked to predict states, while many scientific settings require modeling how uncertainty moves through constrained dynamics. The relevant object is transport over physically admissible futures. This motivates flow learners: models that parameterize transport vector fields and generate trajectories through integration, echoing the continuous dynamics that define PDE evolution. This physics-to-physics alignment supports continuous-time prediction, native uncertainty quantification, and new opportunities for physics-aware solver design. We explain why transport-based learning offers a stronger organizing principle for learned PDE solving and outline the research agenda that follows from this shift.

Yilong Dai, Shengyu Chen, Xiaowei Jia +1

The University of Alabama, Tuscaloosa, Alabama, USA · University of Pittsburgh, Pittsburgh, Pennsylvania, USA

Neural PDE Solvers PDE Solving Flow-Based Generative Modeling

Ask PDF

Abstract

PDE Solving

Explore similar work

A Greedy PDE Router for Blending Neural Operators and Classical Methods

Learning universal approximations for partial differential equations with Physics-Informed Broad Learning System

Flow Learners for PDEs: Toward a Physics-to-Physics Paradigm for Scientific Computing