Discovery of Unstable Singularities

Abstract

Whether singularities can form in fluids remains a foundational unanswered question in mathematics. This phenomenon occurs when solutions to governing equations, such as the 3D Euler equations, develop infinite gradients from smooth initial conditions. Historically, numerical approaches have primarily identified stable singularities. However, these are not expected to exist for key open problems, such as the boundary-free Euler and Navier-Stokes cases, where unstable singularities are hypothesized to play a crucial role. Here, we present the first systematic discovery of new families of unstable singularities. A stable singularity is a robust outcome, forming even if the initial state is slightly perturbed. In contrast, unstable singularities are exceptionally elusive; they require initial conditions tuned with infinite precision, being in a state of instability whereby infinitesimal perturbations immediately divert the solution from its blow-up trajectory. In particular, we present multiple new, unstable self-similar solutions for the incompressible porous media equation and the 3D Euler equation with boundary, revealing a simple empirical asymptotic formula relating the blow-up rate to the order of instability. Our approach combines curated machine learning architectures and training schemes with a high-precision Gauss-Newton optimizer, achieving accuracies that significantly surpass previous work across all discovered solutions. For specific solutions, we reach near double-float machine precision, attaining a level of accuracy constrained only by the round-off errors of the GPU hardware. This level of precision meets the requirements for rigorous mathematical validation via computer-assisted proofs. This work provides a new playbook for exploring the complex landscape of nonlinear partial differential equations (PDEs) and tackling long-standing challenges in mathematical physics.

Figures (6)

• Figure 1: Research flowchart. Our research methodology consists of two main stages. a) Solution discovery. We start with a candidate solution that allows us to search for self-similar spatial profiles of blow-up solutions characterized by a self-similar scaling rate $\lambda$, illustrated in (i) with Burgers' Equation. We then refine our machine learning pipeline (ii) and solution accuracy using an iterative approach. Empirical results from candidate solutions (iii) and their accuracy guide the mathematical modeling and neural network architecture. Mathematical modeling in turn guides the inductive biases built into the network architecture, such as input coordinate transforms and the shaping of the output fields. We use a Physics-informed Neural Network (PINN) with a Gauss-Newton optimiser and a multi-stage refinement training scheme to generate a highly accurate candidate solution while finding the right scaling rate $\lambda$. b) Solution analysis. We analyze the stability of a high-precision candidate solution by linearising the PDE around it. We discover unstable modes, along which any slight perturbation causes the system to deviate from the blow-up trajectory. Thus, we characterize the degree of stability and confirm the discovery of highly accurate stable and unstable singularities.
• Figure 2: Self-similar singularities to IPM and Boussinesq.a, b, Spatial profile of the vorticity $\Omega$ for the third unstable solutions of the IPM (in a) and Boussinesq equations (in b) near the origin. c, d, The cross section of the vorticity profiles along the $y_1$ axis, namely $\Omega(y_1, 0)$ for the stable to third unstable solutions of the IPM (in c) and Boussinesq equations (in d). e, Inverse scaling rate versus the instability order, the plot reveals a linear trend for both systems (shaded area denotes one standard deviation). f, The scaling parameters of the unstable solutions discovered in the paper (green and yellow) and the previously discovered solutions (blue). The significant digits of the scaling parameters are marked in bold-face. The 4th unstable solution of the Boussinesq equation is un-validated.
• Figure 3: Solution accuracy.a, Spatial profiles of the equation residual for the 1st unstable solution to the IPM equation after multistage training. b, List of the log-10 maximum equation residuals for the best solutions discovered for the CCF, IPM and Boussinesq equations.
• Figure 4: High-precision training.a, The log-10 maxium residual value of the solution evaluated on the evaluation grid for the CCF first unstable solution when training with different optimizers. The Gaussian-Newton optimizer shows better performance and significantly faster convergence than the classical Adam + L-BFGS optimizers. Multi-stage training is able to further improve the results by five orders of magnitude. b, The two residuals of the CCF first unstable solution obtained by multi-stage training plotted against the self-similar spatial coordinate $y$.
• Figure 5: Smoothness signal in $\lambda$. Maximum residuals over the input space as a function of $\lambda - \lambda_i$ for admissible $\lambda_i$ of the Boussinesq equation with error bars. The plotted values and error bars are produced based on five experiments with varying seeds. Each subplot is centered around a different admissible lambda value, where $\lambda_0$ is the self-similar scaling parameter for the stable solution, and $\lambda_n$ with $n \geq 1$ is that for the $n$-th unstable solution. The x-axis shows the difference between the lambda value and the closest admissible lambda value ($\lambda - \lambda_i$), while the y-axis shows the mean and standard deviation of the log 10 maximum residual at 400k training steps. The plots show how the residuals increase as the lambda value deviates from $\lambda_i$ for each $i$, indicating that $\lambda_i$ is close to the local minimum of smoothness errors.