PaperlayerView exact revision on arXiv

Source-linked AI summary

Discovery of Singularities in Partial Differential Equations Using Adaptive Numerical Discretizations

Allen Yang, Lucas Amoudruz, Petros Koumoutsakos

arXiv:2610.30492v1math.NAcs.NAmath.AP

TL;DR

Self-similar profiles are used to study finite-time singularities, but unbounded domains and localized sharp features challenge a single spectral interpolant. This paper uses adaptive spectral ODIL with overlapping local and mapped asymptotic patches to compute the profiles. It reports high accuracy for Burgers and generalized De Gregorio profiles and residuals below 1.2 × 10^-8 in the near field for two-dimensional Boussinesq profiles.

  • Problem

    A single Chebyshev interpolant cannot efficiently resolve profiles on an unbounded domain that may have sharp localized features and decay slowly or grow algebraically.

    Source: §2.1 · p. 3

  • Method

    The method optimizes residuals over nodal values on overlapping Chebyshev patches, including a mapped asymptotic patch, and uses coefficient decay to guide adaptive refinement.

    Source: §2.1 · p. 4, §2.3 · p. 5, §4 · p. 15

  • Results

    Burgers errors approach floating-point precision, generalized De Gregorio residuals reach machine-epsilon levels for CLM and stable and first unstable CCF branches, and Boussinesq near-field residuals are below 1.2 × 10^-8.

    Source: §3.1 · p. 8, §3.2 · p. 10, §3.3 · p. 14

  • Takeaways & Limitations

    The computations show that adaptive spectral ODIL can produce accurate self-similar profiles on unbounded domains without domain truncation or prescribed far-field conditions.

    Source: §4 · p. 15

    For profiles with sharp interior features, high-degree differentiation across narrow overlaps may limit attainable accuracy through ill-conditioning.

    Source: §3.2 · p. 10

Abstract

from arXiv · show

We compute self-similar singularity profiles of partial differential equations (PDEs) using adaptive numerical discretizations based on optimizing a discrete loss (ODIL). The fields are represented in overlapping local Chebyshev patches joined with a mapped asymptotic patch. After factoring out the leading power laws of solution ansatzes, their limiting amplitudes at infinity are determined by the optimization, eliminating domain truncation and prescribed far-field conditions. Adaptive refinement is guided by Chebyshev-coefficient decay, and the nonlinear least-squares problem is solved using automatically differentiated Gauss--Newton iterations with Levenberg--Marquardt damping. The method achieves O(10^{-14})-level accuracy for Burgers profiles, residuals of O(10^{-14}) to O(10^{-10}) for the generalized De Gregorio family, and of O(10^{-8}) in a local region for the two-dimensional Boussinesq equations. These results demonstrate that ODIL can identify singularity profiles for PDEs with high accuracy using adaptive numerical methods, complementing ongoing efforts that rely on artificial intelligence.

1 Introduction

Self-similar profiles provide a framework for studying finite-time PDE singularities, and residual-based optimization has been used to discover them. This paper develops a classical adaptive spectral formulation of ODIL for computing such profiles and assessing their numerical resolution.

1 Introduction

  • Prior inverse methods including PINNs have recovered smooth self-similar profiles and unstable branches that are difficult to obtain through forward time evolution.Earlier numerical profiles have also contributed to computer-assisted blowup proofs for the two-dimensional Boussinesq and three-dimensional axisymmetric Euler equations with boundaries.Source: §1 · p. 1
  • ODIL is formulated with nodal field values as optimization variables, so the discrete representation adds no nonlinearity beyond the profile equations while retaining residual-based optimization.The paper presents this as a classical numerical formulation with explicitly controlled resolution and differentiation properties.Source: §1 · pp. 1–2
  • The method supports profile initialization, transfer between resolutions, direct constraints, and overdetermined least-squares fitting without increasing the representation size.Chebyshev-coefficient decay, independent residual evaluations, and comparisons across refinements help assess resolution and asymptotic behavior.Source: §1 · p. 2
  • The framework is assessed on Burgers, generalized De Gregorio, and two-dimensional Boussinesq profiles, spanning an exact-solution benchmark, nonlocal coupling, and a multidimensional system.The introduction also notes reduced training time relative to PINNs and the potential for efficient continuation and sensitivity studies.Source: §1 · p. 2

2 Methodology

The methodology uses adaptive spectral ODIL to solve stationary self-similar profile equations on unbounded domains. Overlapping Chebyshev patches represent local and far-field behavior, while coefficient decay guides one-dimensional refinement and the same optimization framework extends to two dimensions.

2 Methodology

  • ODIL minimizes profile-equation residuals over discretized field values on overlapping Chebyshev patches, using collocation residuals and damped Gauss–Newton iterations.The method also adapts patch decomposition and polynomial degrees using Chebyshev-coefficient decay.Source: §2 · p. 2

2.1 Spectral formulation of ODIL using Chebyshev polynomial interpolation

  • Overlapping patches with a smooth partition of unity address slowly decaying or growing far fields and localized sharp features while keeping the composite profile smooth.A mapped asymptotic patch factors out the leading power law and approximates the normalized field; residuals include oversampled local collocation, local residual-derivative penalties, and overlap matching.Source: §2.1 · pp. 3–4

2.2 Extension to two dimensions

  • In two dimensions, tensor-product Chebyshev patches cover the near field, while overlapping angular patches represent the polar-coordinate far field and are blended smoothly with the near-field patches.The resulting composite representation is smooth, and the one-dimensional residual and Gauss–Newton/Levenberg–Marquardt components carry over.Source: §2.2 · pp. 4–5

2.3 Adaptivity

  • Chebyshev-coefficient decay indicates patch resolution and guides alternating optimization with split, merge, and degree-reduction passes.In this paper, refinement is applied only to local patches in one-dimensional problems; two-dimensional experiments use no adaptive refinement.Source: §2.3 · pp. 5–6

3 Numerical Results

The numerical results assess ODIL on self-similar profiles for Burgers, the generalized De Gregorio family, and the two-dimensional Boussinesq equations, reporting sampled residuals and profile errors. For Boussinesq, they also examine near- and far-field residuals without an imposed far-field condition and compare local results with a PINN-based approach.

3 Numerical Results

  • Validation uses 10^6 random points per region and reports sampled ℓ∞ error and RMSE; Burgers profile errors are relative in the infinite-domain region because its profiles grow without bound.Burgers residuals are absolute in both regions, while profile errors are absolute locally and relative on the infinite domain.Source: §3 · p. 7, Table 1 · p. 8

3.1 The Burgers equation

  • Burgers branch inference recovers λ = 1/2, 1/4, and 1/6 with absolute errors 1.8×10^-9, 1.1×10^-8, and 4.7×10^-7, respectively.Maximum sampled residuals of the extended profiles are of order 10^-13 or below; the authors describe systematic study of λ inference as future work.Source: §3.1 · p. 8
  • For all three Burgers profiles, ODIL achieves smaller errors than PINNs while using substantially fewer parameters and requiring several orders of magnitude less training-iteration time.The computed residuals and errors approach floating-point precision.Source: §3.1 · p. 8, Table 2 · p. 9

3.2 The generalized De Gregorio family

  • Table 3 reports residuals and profile errors for CCF and CLM over 106 points per validation region; infinite-domain samples use y=sinh(z), with z uniform on [-40,40], excluding points near Hilbert-transform quadrature nodes.Sharp features of the second unstable CCF profile require high polynomial degrees and narrow patch overlaps; these may amplify differentiation errors and limit attainable residual accuracy.Source: Table 3 · p. 10, §3.2 · p. 10
  • Resolving the sharp second unstable CCF profile requires high polynomial degrees and narrow patch overlaps, which may amplify differentiation errors and limit attainable residual accuracy.The passage proposes representing derivatives as primary unknowns and recovering fields by integration as a possible conditioning remedy.Source: §3.2 · p. 10

3.3 The two-dimensional Boussinesq equations

  • For Boussinesq, each field uses a local patch on [0, 15]^2 and an asymptotic patch starting at r_start = 10, without adaptivity; residuals are evaluated at prescribed λ = 1.9205599746927.The residuals are visualized on the parity-extended fields.Source: §3.3 · p. 13, Figure 4 · p. 14
  • Within the local region, residual levels are comparable to the PINN-based approach in, with limited discretization tuning and no adaptive refinement.The authors identify extending adaptive refinement to two dimensions as a goal for computing high-precision self-similar profiles.Source: §3.3 · p. 14
  • Preliminary results for singularity formation in the unforced three-dimensional Euler equations [14] use Chebyshev ODIL on a truncated domain.This is presented in Supplementary Section S7.Source: §3.3 · p. 14
  • Near-field residuals have maxima below 1.2 × 10^-8 and RMSE below 4 × 10^-10; infinite-domain maxima are at most 1.0 × 10^-7, attained by f5.The infinite-domain assessment extends to |y| ≈ 3.4 × 10^8, with no far-field condition imposed.Source: §3.3 · p. 14, Table 4 · p. 15

4 Conclusions

The paper computes self-similar singularity profiles on unbounded domains with adaptive spectral ODIL, determining far-field amplitudes through optimization rather than truncating the domain or prescribing far-field conditions. Results include near-machine-precision Burgers profiles, generalized De Gregorio residuals of O(10^-14) to O(10^-10), and Boussinesq residuals of O(10^-8).

4 Conclusions

  • Burgers profiles reach near machine precision, generalized De Gregorio residuals range from O(10^-14) to O(10^-10), and the two-dimensional Boussinesq residual is O(10^-8).Where timing comparisons are available, the computational cost is substantially below PINN-based approaches.Source: §4 · p. 15
  • Sharp interior features can limit accuracy through conditioning of high-degree differentiation across narrow overlaps; an integral reformulation may alleviate this.Adaptive refinement has so far been applied only in one dimension, with extension to two dimensions identified as a next step.Source: §4 · p. 15
  • The approach targets stationary similarity profiles whose normalized fields are smooth in suitable mapped coordinates, while its numerical techniques are not specific to self-similar blowup.Potential extensions include heat and thin-film profiles, as well as dispersive profiles that may require factoring out phase as well as amplitude.Source: §4 · p. 15

S1 Methodology

The supplementary methodology details Chebyshev interpolation, overlapping local and asymptotic patches, and the residual-based optimization used to represent profiles on unbounded domains. It also specifies the solvers and coefficient-based adaptive refinement.

S1.1 Chebyshev interpolation

  • Chebyshev interpolation supplies polynomial approximations and derivatives; coefficient decay reflects smoothness and supports resolution assessment.The interpolation error decreases faster than any fixed power of the node count for sufficiently smooth functions.Source: S1.1 Chebyshev interpolation · pp. 19–20

S1.2 One-dimensional multi-patch representation

  • The one-dimensional representation blends overlapping local patches with a mapped asymptotic patch, factoring out an envelope when the field grows at infinity.The mapping exponent is chosen to avoid fractional-power endpoint behavior that limits regularity.Source: S1.2 One-dimensional multi-patch… · pp. 20–21

S1.3 Two-dimensional representation

  • The two-dimensional representation blends tensor-product Cartesian patches near the origin with angular polar patches in the far field using a radial transition.The polar far-field construction avoids regularity loss from mapping both coordinates separately when both tend to infinity.Source: S1.3 Two-dimensional representation · p. 21

S1.4 Loss and linear solvers

  • The residual loss combines equation residuals, local residual-gradient penalties, and patch-matching conditions, with block scaling to make contributions from different resolutions comparable.The asymptotic patch uses collocation residuals but no gradient penalty.Source: S1.4 Loss and linear solvers · p. 22
  • Damped Gauss–Newton subproblems use dense QR for most problems and preconditioned LSMR for the finest Boussinesq and Euler computations.The latter avoids storing the Jacobian; unpreconditioned LSMR did not reach the required accuracy within a practical iteration count.Source: S1.4 Loss and linear solvers · p. 22

S1.5 Adaptive refinement

  • Adaptive refinement uses normalized Chebyshev coefficient tails to split, merge, or reduce patch degree, but computed tails provide a resolution diagnostic rather than a rigorous error bound.For Burgers and De Gregorio, the procedure holds the local-region endpoint fixed and does not adapt the asymptotic patch.Source: S1.5 Adaptive refinement · pp. 22–23

S2 Parity extension

The computed half-domain fields are extended across symmetry axes by blending them with their parity reflections near each axis. This construction enforces parity exactly while matching the computed field away from the blend band.

S2 Parity extension

  • A narrow blend between the computed field and its parity reflection avoids derivative discontinuities that direct reflection could introduce.The blend uses the computed field near the axis and its reflected counterpart outside the band.Source: S2 Parity extension · p. 23
  • The extension satisfies F(−x) = sF(x) exactly and its first and second derivatives satisfy parity relations to within rounding.The reported deviations in derivative parity are at most 4.5 × 10^-16.Source: S2 Parity extension · p. 23
  • For λ-inference, the origin patch is re-interpolated across the axis as a single polynomial before the residual is evaluated.The re-interpolation uses parity-reflected values and twice as many Chebyshev nodes.Source: S2 Parity extension · p. 23

S3 Far-field and mapping exponents

Far-field envelopes and the mapping exponent are determined from the asymptotic structure of the profile equations. The leading power is prescribed, while limiting amplitudes remain optimization unknowns.

S3 Far-field and mapping exponents

  • The far-field exponent σq follows from dominant balance, and only this leading power is prescribed; the limiting amplitude and two-dimensional angular dependence remain optimization unknowns.A power law q ∼ y^σ annihilates the dominant linear combination when σ = −c/(1 + λ).Source: S3 Far-field and mapping exponents · p. 24
  • The mapping exponent α is chosen to make correction exponents integer-valued in the reference variable, avoiding fractional-power branch points that limit Chebyshev coefficient decay.The selection aims to make as many correction exponents as possible integer multiples of α.Source: S3 Far-field and mapping exponents · p. 24
  • For generalized De Gregorio and Boussinesq, the normalized equations are analytic in s = y^-1/(1+λ), or its radial analogue, so α = 1/(1 + λ).The normalized initial value at s = 0 is free.Source: S3 Far-field and mapping exponents · p. 24
  • For Burgers, the normalized profile is analytic in y^-1/(1+λ), and the smooth branches satisfy the integer-power criterion with α = λ/(1 + λ).The stated branch condition is 1/λ = 2i + 2.Source: S3 Far-field and mapping exponents · p. 24

S4 The Burgers equation

The Burgers profiles are computed on adaptive local patches plus an asymptotic patch, with the three reported configurations requiring at most two local patches. The section also describes a PINN implementation used for comparison.

S4 The Burgers equation

  • The Burgers profiles for λ = 1/2, 1/4, and 1/6 use adaptive local patches and one asymptotic patch, with at most two local patches in each final discretization.The configurations and final discretizations are listed in Table S3.Source: S4 The Burgers equation · p. 25, Table S3 · p. 25

S4.1 Comparison with PINNs

  • The PINN comparison uses a two-stage tanh-network implementation, adding fixed random Fourier features in the second stage and optimizing with stochastic Gauss–Newton.The second-stage correction is scaled using the first-stage maximum residual; the Fourier-feature scale is fixed at σ = 2π.Source: S4.1 Comparison with PINNs · pp. 25–26
  • The PINN validation evaluates maximum residuals and relative maximum errors over the infinite-domain validation region across different network sizes.The supplied table caption specifies these comparison metrics but not their values.Source: Table S4 · p. 26

S5 The generalized De Gregorio family

The generalized De Gregorio profiles are computed with a two-stage optimization that fixes their scaling, alongside a quadrature designed to evaluate the Hilbert transform accurately near the kernel singularity and in the far field.

S5.1 Normalization and two-stage solve

  • A two-stage solve fixes the profile scaling first with a penalty and then by constraining the far-field amplitude.The ccf2 solve instead fixes scaling only in stage 2, without a stage-1 penalty.Source: S5.1 Normalization and two-stage solve · p. 27

S5.2 The Hilbert transform

  • The Hilbert-transform integral is rewritten to remove its kernel singularity, with a finite diagonal limit and a separate treatment of the origin self-term.The self-term at y = 0 uses K(0, 0) = −2 ∂yΩ(0).Source: S5.2 The Hilbert transform · p. 28
  • Clenshaw–Curtis quadrature uses a separate grid with target-specific panels, local segments, and a mapped segment for the tail.Targets are grouped by magnitude so each panel places y0 within its resolved region; near-node cancellation is avoided by interpolating local contributions to targets.Source: S5.2 The Hilbert transform · p. 28
  • The mapped-tail contribution is integrated analytically and differentiated with respect to the far-field amplitude, while quadrature accuracy is checked by doubling panel node counts.The tail construction requires α ≤ 1, which holds for all cases in Table S5.Source: S5.2 The Hilbert transform · p. 28

S6 The two-dimensional Boussinesq equations

The Boussinesq profile is formulated in vorticity form and computed with local and asymptotic patches, transferring a coarse solution to finer grids for the reported profile.

S6 The two-dimensional Boussinesq equations

  • The Boussinesq computation uses one local and one asymptotic patch per field and does not use adaptive refinement.The reported profile comes from continuation from a coarse to a fine discretization.Source: S6 The two-dimensional Boussinesq… · p. 28

S6.1 Formulation

  • The five-field vorticity formulation treats temperature gradients and velocity components as independent unknowns, enforcing gradient, divergence-free, and vorticity relations through residuals.It avoids a stream-function, Poisson, or Biot–Savart evaluation.Source: S6.1 Formulation · pp. 28–29

S6.2 Discretization

  • Each field uses a Cartesian local patch and a polar asymptotic patch whose unknown boundary values represent the angular profile at infinity.The patches blend radially over 10 ≤ r ≤ 15, and residuals are collocated on a grid refined by κ = 2.Source: S6.2 Discretization · p. 29

S6.3 Optimization

  • The reported profile is obtained by solving a coarse problem, transferring its solution to fine grids, and using that transfer as the fine solve’s initial guess.The coarse solve starts from zero and uses dense QR; the fine solve uses sketch-preconditioned LSMR.Source: S6.3 Optimization · p. 29

S7 The three-dimensional axisymmetric Euler equations

A preliminary axisymmetric Euler profile is computed with mapped, truncated domains and a far-field penalty; its residuals are reported as comparable to or better than those in, despite resolution and tail limitations.

Motivation and scope

  • The computation is a preliminary step toward extending the method to singularity formation in the unforced three-dimensional Euler equations [15].The Euler computation uses a mapped, truncated domain and imposes far-field behavior as a penalty rather than using an asymptotic patch.Source: S7 The three-dimensional axisymmetric… · pp. 29–30

Formulation

  • For λ = 1/2, the far-field decay rates are U ∼ R^-2, Ω ∼ R^-3, and Ψ ∼ R^-1, with an unknown axial translation constant C.The profile variables are U = uθ/r, Ω = ωθ/r, and Ψ = ψθ/r.Source: S7 The three-dimensional axisymmetric… · p. 30

Discretization

  • The fields use overlapping tensor-product Chebyshev patches on a mapped computational box truncated at sinh(30) ≈ 5.3 × 10^12.Each patch has 64 × 64 nodes; residuals are collocated on a grid refined by κ = 2.Source: S7 The three-dimensional axisymmetric… · p. 30

Loss and optimization

  • The least-squares loss combines profile-equation residuals and their gradients with patch matching, axis parity, far-field, and normalization terms.The unknowns include U, Ω, Ψ, and C; normalization removes scaling freedom and makes the zero field a nonstationary point of the loss.Source: S7 The three-dimensional axisymmetric… · pp. 30–31
  • The solve uses LM iterations with sketched LSMR and converges to translation constant C = −6.609.The normalization permits initialization from q ≡ 0 and C = 0.Source: S7 The three-dimensional axisymmetric… · p. 31

Results

Loading 2610.30492v1…