---
title: 'Neural PDEs: Methods and Operator Learning'
url: https://www.emergentmind.com/topics/neural-partial-differential-equations-neural-pdes
type: topic
---

# Neural PDEs: Methods and Operator Learning

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 [2111.07671]. 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 [2111.07671].

## 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 [1912.04737]. 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 [2009.03892]. 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 [2111.07671].

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 [2402.03541], multiplicative linear–nonlinear fusion operators on regular and irregular domains [2603.24143], 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 [2407.06664].

The literature also includes problem-class-specific extensions. Stochastic forcing motivates Neural SPDEs based on mild solutions and spectral semigroup evolution [2110.10249]. 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 [2203.17128]. Weak-form methods such as ParticleWNN localize test functions to particle-centered compact supports and use an R-adaptive strategy [2305.12433]. Geometry-aware solvers extend Neural PDEs to irregular immersed domains [2211.03241] and to surface PDEs posed directly on neural shape representations [2512.21311].

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 [2111.07671]. 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})\) [2111.07671].

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\) [2111.07671]. 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 [2111.07671]. 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 [1804.04272]. The connection is architectural rather than application-specific, but it established a PDE lens for stability, reversibility, and discretization design in deep networks [1804.04272].

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 [2111.07671]. Training is end-to-end through the solver, with gradients computed via the adjoint sensitivity method from Neural ODEs [2111.07671]. 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 [2111.07671].

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 \(M\) time steps are predicted, with parameter sharing across time steps and an analogy to explicit time marching [2009.03892]. 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\) directly with a neural network \(u_\theta\) 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 [1912.04737], in small-network collocation studies for Poisson and steady Navier–Stokes [1904.07200], and in acceleration-oriented variants inspired by NeRF and integral forms of Maxwell equations [2205.08978]. In these methods, automatic differentiation supplies \(\partial_t u_\theta\), \(\nabla u_\theta\), and \(\nabla^2 u_\theta\), while collocation points are sampled in the domain and on the boundary [1912.04737].

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 [2305.12433]. For a Poisson-type problem with Dirichlet boundary data, the weak residual under a localized test function \(\varphi_i\) is
\[
\mathcal{R}(u_{NN}; \varphi_i) := \int_{B(\mathbf{x}_i^c, R_i)} \nabla u_{NN} \cdot \nabla \varphi_i \, \mathrm{d}\mathbf{x} \;-\; \int_{B(\mathbf{x}_i^c, R_i)} f\, \varphi_i \, \mathrm{d}\mathbf{x},
\]
and training minimizes the mean squared local residuals plus a boundary penalty [2305.12433]. Its R-descending strategy adaptively reduces support radii during training, and the design lowers derivative order requirements via integration by parts [2305.12433].

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 [2207.10241]. The loss is the squared norm of Galerkin residuals against the basis functions, so no labeled solution data are required [2207.10241]. For singularly perturbed convection–diffusion equations, the basis is enriched by a boundary-layer corrector
\[
\phi_{cor}(x) = \exp(−(1 + x)/\epsilon) - \Big[ 1 - (1 - \exp(−2/\epsilon)) (x + 1)/2 \Big],
\]
which addresses sharp boundary layers that challenge strong-form PINNs [2207.10241].

Integral formulations provide another structurally informed family. Q-PDE for elliptic-type PDEs uses a boundary-enforcing ansatz
\[
Q^N(x;\theta) := S^N(x;\theta)\,\eta(x) + f(x),
\]
where \(\eta\) vanishes on the boundary, and trains with a truncated residual loss together with a biased Q-learning-style gradient estimator that differentiates only through \(Q^N\), not through \(\mathcal{L}Q^N\) [2203.17128]. The theoretical analysis derives an infinite-width ODE
\[
\frac{dQ_t}{dt} = \alpha\,\mathcal{B}\,\mathcal{L}(Q_t),
\]
and proves \(L^2\)-convergence to the PDE solution for strongly monotone operators [2203.17128].

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 [2602.12872]. 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,
\[
u_{raw} = \alpha_\theta \cdot ( \mathcal{B}_{L,\theta}(z) \odot \mathcal{B}_{N,\theta}(z) ),
\]
with separate encoders for multiple functional inputs and an optional decoder for spatial consistency [2603.24143]. 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 [2603.24143]. On the tested 3D Poisson–Boltzmann case, it attains the best accuracy among the compared models and trains approximately \(2.7\times\) faster than a 3D FNO baseline [2603.24143].

HAMLET adopts a graph transformer architecture with modular input encoders and sparse neighborhood attention on arbitrary discretizations [2402.03541]. 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 [2402.03541]. 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 [2402.03541].

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 [2402.12652]. A graph Transformer produces layer-wise modulation embeddings that condition an implicit neural representation, yielding mesh-free predictions \( \hat{u}(t,x) \) at arbitrary coordinates [2407.06664]. 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 [2407.06664].

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 [2111.07671]. 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 [2512.21311]. 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 [2512.21311]. The method is trained once on a single representative shape and generalizes across shapes and topologies [2512.21311].

## 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 [2111.07671]. The same paper explicitly notes that irregular meshes, non-equidistant sampling, and complex geometries are not supported by its CNN parameterization [2111.07671]. 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
\[
\chi_w(x) = \sum_{i=1}^n a_i \frac{(p_i - x)\cdot \hat{n}_i}{4\pi \|p_i - x\|^3},
\]
and a signed-distance formulation derived from a stabilized Eikonal equation [2211.03241]. 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 [2211.03241]. The reported comparisons show markedly lower errors and shorter training times than PINNs on several 2D Poisson benchmarks [2211.03241].

Boundary handling is also explicit in Q-PDE through a smooth auxiliary function \(\eta\) that vanishes on \(\partial\Omega\), ensuring exact Dirichlet satisfaction by construction [2203.17128]. 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 [2305.12433]. 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 [2602.12872].

Stochasticity adds another dimension. Neural SPDEs formulate semilinear SPDEs through their mild solution,
\[
u(t) = S(t) u_0 + \int_0^t S(t-s) F(u(s), s)\, ds + \int_0^t S(t-s) G(u(s), s)\, dW_s,
\]
with \(S(t)=e^{tL}\) the semigroup of the linear generator [2110.10249]. 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 [2110.10249]. The architecture conditions jointly on the initial condition and the driving-noise realization, which extends deterministic neural operators and continuous-time sequence models [2110.10249].

Irregular domains also motivate graph-based operator learning. HAMLET uses sparse graph attention over radius-defined neighborhoods and cross-attention to arbitrary query points [2402.03541], while LNF-NO can operate decoder-free on node-based outputs for irregular meshes [2603.24143]. 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 [2111.07671]. 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 [2111.07671]. On real-world data, NeuralPDE is second-best on Ocean Wave with RMSE \(0.435\) versus \(0.427\) for the best model and on WeatherBench with \(0.097\) versus \(0.096\), while ResNet performs best on PlaSim [2111.07671].

The RNN-based Neural-PDE reports low MSE values across wave, heat, Burgers, and NSCH examples, including \(3.54\times 10^{-5}\) for the 1D wave equation, \(2.1551\times 10^{-6}\) for the 2D heat equation, \(2.3070\times 10^{-6}\) for 2D inviscid Burgers, and \(6.1631\times 10^{-7}\) for NSCH [2009.03892]. It uses only two biLSTM layers of 20 units each plus a small dense head and emphasizes parameter sharing across time [2009.03892].

PDEformer-1 reports zero-shot relative \(L^2\) errors competitive with dedicated expert models on PDEBench, for example \(0.0061\) on Burgers with \(\nu=0.1\), \(0.0180\) on Burgers with \(\nu=0.01\), and \(0.0120\) on Advection with \(\beta=0.1\), with further gains after finetuning [2407.06664]. 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 [2305.12433]. Q-PDE demonstrates relative error below \(0.02\%\) in a 1D elliptic example and operator loss below \(10^{-5}\) in a 20D case, with approximately \(5\times\) faster training than DGM in that setup [2203.17128]. NEKM reports rel \(L^2\) errors around \(0.1\%\)–\(0.5\%\) on heat and wave equations and fast batched inference for \(10{,}000\) heat equations in parallel [2602.12872].

These successes coexist with clear limitations. NeuralPDE is confined to uniform Cartesian grids and periodic boundaries [2111.07671]. 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 [1904.07200]. Strong-form methods may struggle with singular perturbations and boundary layers, which motivates enriched spectral bases in ULGNet [2207.10241]. Operator learners can lose accuracy on transport-dominated or sharply discontinuous regimes, as seen for LNF-NO on Navier–Stokes and piecewise-coefficient Darcy [2603.24143]. Irregular-domain and neural-geometry methods inherit sensitivity to normals, curvatures, medial-axis ambiguities, or occupancy estimation [2211.03241].

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 [2111.07671]. ParticleWNN points toward more complex domains and higher-dimensional weak-form solvers [2305.12433]. PDEformer highlights the extension from 1D to multi-dimensional PDEs, more operator vocabularies, and coupled multi-physics systems [2407.06664]. NEKM suggests surrogate boundary kernels and fast integral approximations for variable coefficients and higher dimensions [2602.12872]. 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 [1912.04737]. 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 [2111.07671]. At the other end are operator-learning systems that amortize inference across entire PDE families, from graph transformers on arbitrary discretizations [2402.03541] to computational-graph foundation models [2407.06664] and multiplicative fusion operators [2603.24143].

Across these variants, the most durable technical ideas are structural. Convolutions inherit finite-difference semantics on regular grids [2111.07671]; graph or transformer encoders represent discretization-agnostic geometry and PDE context [2402.03541]; weak forms reduce derivative order and improve numerical stability [2305.12433]; spectral semigroups support resolution-invariant stochastic evolution [2110.10249]; and boundary integral or immersed formulations offer principled boundary handling on irregular domains [2602.12872]. 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 [2111.07671].

Source: https://www.emergentmind.com/topics/neural-partial-differential-equations-neural-pdes