Neural PDEs: Methods and Operator Learning

• Neural PDEs are a spectrum of techniques that merge neural networks with traditional PDE structures to model dynamics, solve equations, and learn operator mappings.
• They utilize methods such as residual-based solvers, continuous-time formulations, and weak-form discretizations to achieve improved accuracy and efficiency.
• Neural PDE frameworks leverage convolutional, graph, and transformer architectures to handle complex geometries, irregular domains, and stochastic processes.

Searching arXiv for recent and foundational papers on Neural PDEs. Partial differential equations can be modeled, solved, or approximated with neural networks in several distinct but overlapping senses. In current research usage, “Neural PDEs” most commonly denotes methods that either learn PDE-governed dynamics from data, solve PDEs by embedding the equation into neural training, or learn solution operators that map functions such as coefficients, boundary data, or initial conditions to PDE solutions (Dulny et al., 2021). The field therefore spans continuous-depth dynamical formulations, residual-based and weak-form solvers, operator-learning architectures, graph- and transformer-based models, and stochastic or geometry-aware extensions. A unifying theme is the replacement or augmentation of hand-crafted discretizations, constitutive operators, or solvers with trainable neural components while preserving some part of PDE structure such as semigroup evolution, boundary integral representations, variational formulations, or the Method of Lines (Dulny et al., 2021).

1. Historical scope and terminological range

The term “Neural PDEs” is not monosemous. One research line uses neural networks as continuous function approximators trained to satisfy PDE residuals and boundary or initial conditions, as in residual-based formulations for Poisson, heat, Burgers, and Navier–Stokes equations (Pedro et al., 2019). Another line treats PDE solving as sequence modeling on spatiotemporal grids, replacing time-stepping rules with recurrent architectures such as bidirectional LSTMs that predict future states from past fields (Hu et al., 2020). A third line learns PDE dynamics directly from data by parameterizing spatial differential operators with convolutional networks and integrating the resulting semi-discrete evolution with differentiable ODE solvers (Dulny et al., 2021).

A separate but closely related trajectory is operator learning. In that setting, the goal is not to solve one PDE instance by residual minimization, but to learn a map between function spaces, for example from coefficient fields, source terms, or boundary traces to the solution field. This perspective underlies neural operators such as graph transformer neural operators on arbitrary geometries (Bryutkin et al., 2024), multiplicative linear–nonlinear fusion operators on regular and irregular domains (Wu et al., 25 Mar 2026), and foundation-style models that ingest the symbolic structure of a PDE through a computational graph and decode mesh-free solutions with an implicit neural representation (Ye et al., 2024).

The literature also includes problem-class-specific extensions. Stochastic forcing motivates Neural SPDEs based on mild solutions and spectral semigroup evolution (Salvi et al., 2021). Elliptic-type PDEs admit reinforcement-learning-inspired formulations such as Q-PDE, which is mesh-free and uses a biased Q-learning-style update with Neural Tangent Kernel analysis (Cohen et al., 2022). Weak-form methods such as ParticleWNN localize test functions to particle-centered compact supports and use an R-adaptive strategy (Zang et al., 2023). Geometry-aware solvers extend Neural PDEs to irregular immersed domains (Khara et al., 2022) and to surface PDEs posed directly on neural shape representations (Welschinger et al., 24 Dec 2025).

This diversity implies that “Neural PDEs” can refer either to a modeling paradigm for dynamical systems, to a class of neural numerical solvers, or to neural operator surrogates. A plausible implication is that the term now names a research area rather than a single method family.

2. Continuous-time dynamics and the Method of Lines perspective

A central formulation treats evolution PDEs as semi-discrete neural dynamical systems. In the NeuralPDE model, the target class is written as

$$\frac{\partial u(x,t)}{\partial t} \;=\; F\big(u, \nabla u, \nabla^2 u, x, t; \theta\big),$$

with $u(x,t)$ possibly multi-channel and $F$ allowed to be nonlinear and depend on spatial derivatives up to any order (Dulny et al., 2021). Applying the Method of Lines to a regular grid produces a continuous-time ODE in the tensor of grid values,

$$\frac{d}{dt}\,\mathcal{U}(t) \;=\; f^\ast\big(t, \mathcal{U}(t)\big),$$

and under time-homogeneity, $f^\ast(t,\mathcal{U}) = f^\ast(\mathcal{U})$ (Dulny et al., 2021).

The key observation is that finite-difference stencils are discrete convolutions. Centered first and second derivatives in one dimension and the standard five-point Laplacian in two dimensions can all be written as convolution kernels scaled by $\Delta x$ or $\Delta y$ (Dulny et al., 2021). NeuralPDE therefore replaces fixed stencils with a learnable convolutional operator $\mathcal{N}_\theta$,

$$\frac{d}{dt}\,\mathcal{U}(t) \;=\; \mathcal{N}_\theta\big(\mathcal{U}(t)\big),$$

so that arbitrary spatially local, translation-equivariant, possibly nonlinear operators can be represented with standard CNN layers (Dulny et al., 2021). This ties the neural architecture directly to the numerical structure of PDE discretization rather than treating the problem as generic sequence prediction.

This continuous-time viewpoint extends earlier residual-network interpretations of deep models. Residual blocks can be read as forward Euler discretizations of ODE or PDE dynamics in a pseudo-time variable, and convolutional operators can be interpreted as reaction, convection, and diffusion terms in parabolic or hyperbolic CNNs (Ruthotto et al., 2018). The connection is architectural rather than application-specific, but it established a PDE lens for stability, reversibility, and discretization design in deep networks (Ruthotto et al., 2018).

Time integration in NeuralPDE is performed by differentiable ODE solvers. The reported experiments use fixed-step explicit Euler,

$$\mathcal{U}^{n+1} \;=\; \mathcal{U}^{n} + \Delta t\;\mathcal{N}_\theta\big(\mathcal{U}^{n}\big),$$

while noting that Runge–Kutta or Dormand–Prince methods are also compatible (Dulny et al., 2021). Training is end-to-end through the solver, with gradients computed via the adjoint sensitivity method from Neural ODEs (Dulny et al., 2021). Stability depends on step size and the Jacobian spectrum of the learned operator, and the reported experiments evaluate RMSE over horizons up to 16 steps, with expected degradation at longer rollouts (Dulny et al., 2021).

A distinct sequence-modeling alternative appears in Neural-PDE, where time-dependent PDEs are cast as sequence-to-sequence learning with bidirectional LSTMs. There, past sampled fields on a grid are encoded and the next $u(x,t)$0 time steps are predicted, with parameter sharing across time steps and an analogy to explicit time marching (Hu et al., 2020). This suggests a conceptual bifurcation within Neural PDEs: one branch models a continuous right-hand side and integrates it, while another directly learns discrete evolution rules.

3. Solver paradigms: strong form, weak form, and integral formulations

Residual-based neural solvers represent the PDE solution $u(x,t)$1 directly with a neural network $u(x,t)$2 and train by minimizing interior residuals together with boundary and initial-condition mismatches. This formulation appears in early general treatments of Poisson, heat, Burgers, and Navier–Stokes problems (Pedro et al., 2019), in small-network collocation studies for Poisson and steady Navier–Stokes (Dockhorn, 2019), and in acceleration-oriented variants inspired by NeRF and integral forms of Maxwell equations (Rzepecki et al., 2022). In these methods, automatic differentiation supplies $u(x,t)$3, $u(x,t)$4, and $u(x,t)$5, while collocation points are sampled in the domain and on the boundary (Pedro et al., 2019).

Weak-form approaches replace pointwise residuals with variational or Petrov–Galerkin conditions. ParticleWNN defines the trial space as neural networks and the test space as compactly supported functions localized to extremely small regions centered at particles (Zang et al., 2023). For a Poisson-type problem with Dirichlet boundary data, the weak residual under a localized test function $u(x,t)$6 is

$u(x,t)$7

and training minimizes the mean squared local residuals plus a boundary penalty (Zang et al., 2023). Its R-descending strategy adaptively reduces support radii during training, and the design lowers derivative order requirements via integration by parts (Zang et al., 2023).

ULGNet takes a spectral Galerkin route. It expands the solution in compact Legendre bases chosen to satisfy homogeneous Dirichlet or Neumann conditions exactly and lets a CNN predict the spectral coefficients from forcing terms or other PDE inputs (Choi et al., 2022). The loss is the squared norm of Galerkin residuals against the basis functions, so no labeled solution data are required (Choi et al., 2022). For singularly perturbed convection–diffusion equations, the basis is enriched by a boundary-layer corrector

$u(x,t)$8

which addresses sharp boundary layers that challenge strong-form PINNs (Choi et al., 2022).

Integral formulations provide another structurally informed family. Q-PDE for elliptic-type PDEs uses a boundary-enforcing ansatz

$u(x,t)$9

where $F$0 vanishes on the boundary, and trains with a truncated residual loss together with a biased Q-learning-style gradient estimator that differentiates only through $F$1, not through $F$2 (Cohen et al., 2022). The theoretical analysis derives an infinite-width ODE

$F$3

and proves $F$4-convergence to the PDE solution for strongly monotone operators (Cohen et al., 2022).

NEKM similarly embeds classical PDE structure into the architecture, but for evolutionary PDEs. It decomposes the elliptic subproblem produced by implicit time stepping into a source-driven volume potential and a boundary-driven double-layer potential, learns the volume operator in a supervised operator-learning fashion, and learns the boundary density through a self-supervised second-kind boundary integral equation (Ling et al., 13 Feb 2026). A plausible implication is that the strongest neural PDE solvers often derive their inductive bias from the algebraic or integral structure of classical numerical methods rather than from generic network flexibility alone.

4. Operator learning, foundation models, and mesh-free representations

Operator-learning approaches seek maps between functions rather than per-instance solutions. LNF-NO represents the operator map as a multiplicative fusion of a linear branch and a nonlinear branch,

$F$5

with separate encoders for multiple functional inputs and an optional decoder for spatial consistency (Wu et al., 25 Mar 2026). It is designed to support regular grids and irregular geometries, and the reported benchmarks include Laplace, Burgers, Darcy, Poisson–Boltzmann, Poisson–Nernst–Planck, Navier–Stokes, and 3D Poisson–Boltzmann problems (Wu et al., 25 Mar 2026). On the tested 3D Poisson–Boltzmann case, it attains the best accuracy among the compared models and trains approximately $F$6 faster than a 3D FNO baseline (Wu et al., 25 Mar 2026).

HAMLET adopts a graph transformer architecture with modular input encoders and sparse neighborhood attention on arbitrary discretizations (Bryutkin et al., 2024). Nodes concatenate parameter values and positions, a radius-based graph defines neighborhoods, and CrossFormer-style cross-attention fuses encoded PDE data with arbitrary query locations (Bryutkin et al., 2024). The model is positioned as discretization-invariant and suitable for uniform grids, irregular meshes, and point clouds, with empirical results on Darcy flow, Shallow Water, Diffusion Reaction, and Airfoil benchmarks (Bryutkin et al., 2024).

PDEformer and PDEformer-1 take a different route by representing the PDE itself as a computational graph whose nodes encode unknown fields, scalar coefficients, coefficient fields, initial conditions, derivatives, algebraic operators, and equality constraints (Ye et al., 2024). A graph Transformer produces layer-wise modulation embeddings that condition an implicit neural representation, yielding mesh-free predictions $F$7 at arbitrary coordinates (Ye et al., 2024). Pretraining on up to three million one-dimensional PDE samples enables zero-shot or few-shot transfer across diffusion–convection–reaction families and beyond, including trigonometric nonlinearities and wave equations after finetuning (Ye et al., 2024).

These models differ from per-instance neural solvers in that they amortize computation across a family of PDEs. They also differ from NeuralPDE-style dynamical models, which learn a single right-hand side from observed trajectories on a fixed grid (Dulny et al., 2021). This suggests a useful distinction between “neural PDE solvers” and “neural PDE operator learners,” even though both fall within the broader Neural PDEs label.

A further extension appears in neural representations of geometry. “Learning to Solve PDEs on Neural Shape Representations” develops a mesh-free surface PDE solver that operates directly on implicit neural fields, spherical neural surfaces, point clouds, Gaussian splats, and meshes (Welschinger et al., 24 Dec 2025). It learns a local update operator enforcing closest-point normal-constancy in a narrow band and then applies standard Cartesian finite differences for heat and Poisson equations on surfaces (Welschinger et al., 24 Dec 2025). The method is trained once on a single representative shape and generalizes across shapes and topologies (Welschinger et al., 24 Dec 2025).

5. Geometry, boundaries, stochasticity, and irregular domains

Boundary conditions and geometry handling are major fault lines in Neural PDE research. NeuralPDE, in the convolutional Method-of-Lines sense, assumes periodic boundary conditions and uniform Cartesian grids, implemented by circular padding (Dulny et al., 2021). The same paper explicitly notes that irregular meshes, non-equidistant sampling, and complex geometries are not supported by its CNN parameterization (Dulny et al., 2021). This is an important limitation because many physically relevant PDEs are posed on irregular domains.

IBN addresses this by combining a regular background grid with a differentiable interior–exterior classifier. Two alternatives are given: a generalized winding number

$F$8

and a signed-distance formulation derived from a stabilized Eikonal equation (Khara et al., 2022). The resulting mask is used in FEM-Galerkin losses over the background grid, while Dirichlet or Neumann conditions are imposed on the irregular boundary through sampled boundary points and normals (Khara et al., 2022). The reported comparisons show markedly lower errors and shorter training times than PINNs on several 2D Poisson benchmarks (Khara et al., 2022).

Boundary handling is also explicit in Q-PDE through a smooth auxiliary function $F$9 that vanishes on $$\frac{d}{dt}\,\mathcal{U}(t) \;=\; f^\ast\big(t, \mathcal{U}(t)\big),$$0, ensuring exact Dirichlet satisfaction by construction (Cohen et al., 2022). ParticleWNN handles Dirichlet conditions by a penalty on sampled boundary points and accommodates natural boundary conditions through weak-form boundary integrals on supports that intersect the boundary (Zang et al., 2023). NEKM uses a second-kind boundary integral equation for the boundary-driven component, thereby reducing the learning target from an interior field to a boundary density (Ling et al., 13 Feb 2026).

Stochasticity adds another dimension. Neural SPDEs formulate semilinear SPDEs through their mild solution,

$$\frac{d}{dt}\,\mathcal{U}(t) \;=\; f^\ast\big(t, \mathcal{U}(t)\big),$$1

with $$\frac{d}{dt}\,\mathcal{U}(t) \;=\; f^\ast\big(t, \mathcal{U}(t)\big),$$2 the semigroup of the linear generator (Salvi et al., 2021). In periodic domains, this becomes diagonal in Fourier space, allowing drift and diffusion terms to be parameterized on spectral modes and evaluated either by a spectral Galerkin solver or by a fixed-point iteration (Salvi et al., 2021). The architecture conditions jointly on the initial condition and the driving-noise realization, which extends deterministic neural operators and continuous-time sequence models (Salvi et al., 2021).

Irregular domains also motivate graph-based operator learning. HAMLET uses sparse graph attention over radius-defined neighborhoods and cross-attention to arbitrary query points (Bryutkin et al., 2024), while LNF-NO can operate decoder-free on node-based outputs for irregular meshes (Wu et al., 25 Mar 2026). A plausible implication is that domain geometry has become one of the principal organizing variables of Neural PDE method design: Cartesian convolution for regular grids, graph or transformer mechanisms for unstructured discretizations, and integral or immersed formulations for complex boundaries.

6. Performance, limitations, and research directions

Reported empirical performance is heterogeneous across paradigms and problem classes. NeuralPDE evaluates on four toy PDEs and three real-world climate or ocean datasets using RMSE averaged over 16-step prediction horizons (Dulny et al., 2021). It attains state-of-the-art accuracy on Burgers, Gas Dynamics, and Wave among the tested toy systems, while PDE-Net 2.0 performs best on the simplest Advection–Diffusion case (Dulny et al., 2021). On real-world data, NeuralPDE is second-best on Ocean Wave with RMSE $$\frac{d}{dt}\,\mathcal{U}(t) \;=\; f^\ast\big(t, \mathcal{U}(t)\big),$$3 versus $$\frac{d}{dt}\,\mathcal{U}(t) \;=\; f^\ast\big(t, \mathcal{U}(t)\big),$$4 for the best model and on WeatherBench with $$\frac{d}{dt}\,\mathcal{U}(t) \;=\; f^\ast\big(t, \mathcal{U}(t)\big),$$5 versus $$\frac{d}{dt}\,\mathcal{U}(t) \;=\; f^\ast\big(t, \mathcal{U}(t)\big),$$6, while ResNet performs best on PlaSim (Dulny et al., 2021).

The RNN-based Neural-PDE reports low MSE values across wave, heat, Burgers, and NSCH examples, including $$\frac{d}{dt}\,\mathcal{U}(t) \;=\; f^\ast\big(t, \mathcal{U}(t)\big),$$7 for the 1D wave equation, $$\frac{d}{dt}\,\mathcal{U}(t) \;=\; f^\ast\big(t, \mathcal{U}(t)\big),$$8 for the 2D heat equation, $$\frac{d}{dt}\,\mathcal{U}(t) \;=\; f^\ast\big(t, \mathcal{U}(t)\big),$$9 for 2D inviscid Burgers, and $f^\ast(t,\mathcal{U}) = f^\ast(\mathcal{U})$0 for NSCH (Hu et al., 2020). It uses only two biLSTM layers of 20 units each plus a small dense head and emphasizes parameter sharing across time (Hu et al., 2020).

PDEformer-1 reports zero-shot relative $f^\ast(t,\mathcal{U}) = f^\ast(\mathcal{U})$1 errors competitive with dedicated expert models on PDEBench, for example $f^\ast(t,\mathcal{U}) = f^\ast(\mathcal{U})$2 on Burgers with $f^\ast(t,\mathcal{U}) = f^\ast(\mathcal{U})$3, $f^\ast(t,\mathcal{U}) = f^\ast(\mathcal{U})$4 on Burgers with $f^\ast(t,\mathcal{U}) = f^\ast(\mathcal{U})$5, and $f^\ast(t,\mathcal{U}) = f^\ast(\mathcal{U})$6 on Advection with $f^\ast(t,\mathcal{U}) = f^\ast(\mathcal{U})$7, with further gains after finetuning (Ye et al., 2024). ParticleWNN reports strong results on high-frequency Poisson, Allen–Cahn, Kovasznay flow, inverse Poisson, and 5D Poisson, often with lower error and lower wall-clock time than vanilla PINN baselines (Zang et al., 2023). Q-PDE demonstrates relative error below $f^\ast(t,\mathcal{U}) = f^\ast(\mathcal{U})$8 in a 1D elliptic example and operator loss below $f^\ast(t,\mathcal{U}) = f^\ast(\mathcal{U})$9 in a 20D case, with approximately $\Delta x$0 faster training than DGM in that setup (Cohen et al., 2022). NEKM reports rel $\Delta x$1 errors around $\Delta x$2–$\Delta x$3 on heat and wave equations and fast batched inference for $\Delta x$4 heat equations in parallel (Ling et al., 13 Feb 2026).

These successes coexist with clear limitations. NeuralPDE is confined to uniform Cartesian grids and periodic boundaries (Dulny et al., 2021). Residual-based neural solvers can be slower than classical methods, sensitive to random initialization, and exhibit weak correlation between residual loss and true solution error, as documented for Poisson in a collocation study (Dockhorn, 2019). Strong-form methods may struggle with singular perturbations and boundary layers, which motivates enriched spectral bases in ULGNet (Choi et al., 2022). Operator learners can lose accuracy on transport-dominated or sharply discontinuous regimes, as seen for LNF-NO on Navier–Stokes and piecewise-coefficient Darcy (Wu et al., 25 Mar 2026). Irregular-domain and neural-geometry methods inherit sensitivity to normals, curvatures, medial-axis ambiguities, or occupancy estimation (Khara et al., 2022).

Several recurring future directions appear across the literature. NeuralPDE explicitly proposes learning boundary conditions from data, exploring spatial modeling beyond discrete convolutions, and using adaptive solvers such as Dormand–Prince (Dulny et al., 2021). ParticleWNN points toward more complex domains and higher-dimensional weak-form solvers (Zang et al., 2023). PDEformer highlights the extension from 1D to multi-dimensional PDEs, more operator vocabularies, and coupled multi-physics systems (Ye et al., 2024). NEKM suggests surrogate boundary kernels and fast integral approximations for variable coefficients and higher dimensions (Ling et al., 13 Feb 2026). This suggests that the field is moving toward hybridization: numerical structure where it is reliable, neural flexibility where classical models are unavailable or too rigid.

7. Conceptual synthesis

Neural PDEs are best understood not as a single algorithmic family but as a spectrum of ways to combine PDE structure with neural function approximation. At one end are residual-based solvers that approximate one solution field directly and rely on automatic differentiation plus collocation or weak forms (Pedro et al., 2019). In the middle are structured dynamical models that learn the right-hand side of a PDE-like evolution and integrate it in continuous time, with NeuralPDE providing a canonical Method-of-Lines realization on regular grids (Dulny et al., 2021). At the other end are operator-learning systems that amortize inference across entire PDE families, from graph transformers on arbitrary discretizations (Bryutkin et al., 2024) to computational-graph foundation models (Ye et al., 2024) and multiplicative fusion operators (Wu et al., 25 Mar 2026).

Across these variants, the most durable technical ideas are structural. Convolutions inherit finite-difference semantics on regular grids (Dulny et al., 2021); graph or transformer encoders represent discretization-agnostic geometry and PDE context (Bryutkin et al., 2024); weak forms reduce derivative order and improve numerical stability (Zang et al., 2023); spectral semigroups support resolution-invariant stochastic evolution (Salvi et al., 2021); and boundary integral or immersed formulations offer principled boundary handling on irregular domains (Ling et al., 13 Feb 2026). This suggests that progress in Neural PDEs is driven less by generic scale than by increasingly precise alignments between neural parameterizations and classical analytical or numerical structure.

The resulting field is therefore both broader and more technically differentiated than simplified accounts imply. “Neural PDEs” encompasses learned continuous-time vector fields, neural weak-form solvers, boundary-integral operator learners, graph and transformer architectures for parametric PDE families, stochastic evolution models, and mesh-free solvers on neural geometry representations. Its central question is no longer whether neural networks can solve or emulate PDEs, but which parts of PDE structure should remain explicit and which should be learned for a given problem class (Dulny et al., 2021).