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

TL;DR

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.

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 (16)

• 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 C e^{-\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 schwab2011sparse 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.