""" Full 2D Spectral Stability of the OpenAI Navier-Stokes Blow-Up Profile ====================================================================== Builds the linearized self-similar operator L on a tensor-product Chebyshev grid in (X, eta), using the profile equations from the OpenAI manuscript (equations 4.7-4.14). The profile F_*(X, eta) satisfies the stationary self-similar equation. On eta = 0 this reduces to a Kummer confluent hypergeometric ODE. Off eta = 0, the profile couples to the axial velocity U through continuity and the centrifugal driving term. The computation verifies whether the spectral gap sigma > 0 survives in the full 2D problem. This is the unverified hypothesis in the bootstrap paper. Requirements: numpy, scipy Run: python navier_stokes_2d_full.py """ import numpy as np from scipy.linalg import eig from scipy.special import hyp1f1 # Kummer M(a, b, z) import json import sys def cheb(N, L_dom): """Chebyshev differentiation matrices on [0, L_dom].""" if N == 0: return np.array([L_dom / 2]), np.array([[0.0]]), np.array([[0.0]]) theta = np.pi * np.arange(N + 1) / N xi = np.cos(theta) x = L_dom * (1.0 + xi) / 2.0 c = np.ones(N + 1) c[0] = 2.0 c[N] = 2.0 c *= (-1.0) ** np.arange(N + 1) X_mat = np.tile(xi, (N + 1, 1)) dX = X_mat - X_mat.T D = np.outer(c, 1.0 / c) / (dX + np.eye(N + 1)) D -= np.diag(D.sum(axis=1)) D1 = D * (2.0 / L_dom) D2 = D1 @ D1 return x, D1, D2 def kummer_profile(X, h, V0): """ The stationary profile F_*(X) on eta = 0. Satisfies: 2X F'' + (4 - (1-V0)X) F' - (h - V0) F = 0 This is Kummer's equation M(a, b, z) with: z = (1 - V0) X / 2 a = (h - V0) / (1 - V0) b = 2 For V0 = 1 + 2h: 1 - V0 = -2h, so z = -h X a = (h - 1 - 2h) / (-2h) = (-1 - h) / (-2h) = (1 + h) / (2h) b = 2 """ a_kum = (1 + h) / (2 * h) b_kum = 2.0 z_kum = -h * X F = np.zeros_like(X) for i, xi in enumerate(X): zi = -h * xi # Kummer M(a, b, z) for real z val = float(hyp1f1(a_kum, b_kum, zi)) F[i] = val # Normalize so F(0) = 1 F0 = float(hyp1f1(a_kum, b_kum, 0.0)) if abs(F0) > 1e-15: F /= F0 return F def build_eta0_operator(NX, Xmax, h, V0): """Build L_0 on eta = 0 (already validated, 1D operator).""" x, D1, D2 = cheb(NX, Xmax) L0 = 2.0 * np.diag(x) @ D2 \ + np.diag(4.0 - (1.0 - V0) * x) @ D1 \ + (-(h - V0)) * np.eye(NX + 1) # BCs: f(Xmax) = 0 at index 0, f(0) = 0 at index NX L0[0, :] = 0.0 L0[NX, :] = 0.0 return x, D1, D2, L0 def build_full_2d_operator(NX, N_eta, Xmax, eta_max, h, V0): """ Build the full 2D linearized operator on a tensor-product Chebyshev grid in (X, eta). The operator acts on perturbations f(X, eta) to the swirl profile. L[f] = (radial diffusion and transport on each eta slice) + (axial diffusion in eta) + (centrifugal coupling between X and eta) Parameters: NX: number of Chebyshev points in X N_eta: number of Chebyshev points in eta Xmax: domain size in X eta_max: domain size in eta (we use [-eta_max, eta_max] but by even symmetry work on [0, eta_max]) h: the small exponent (0 < h < 1/100) V0: = 1 + 2h (axial compression rate) """ A = 0.5 + h D_exp = 0.5 - h # X grid x_pts, DX1, DX2 = cheb(NX, Xmax) # eta grid on [0, eta_max] (even symmetry: only need eta >= 0) eta_pts, Deta1, Deta2 = cheb(N_eta, eta_max) n_total = (NX + 1) * (N_eta + 1) # Index mapping: state vector is f[i, j] where i = X index, j = eta index # Flattened as f[i * (N_eta + 1) + j] def idx(i, j): return i * (N_eta + 1) + j L_full = np.zeros((n_total, n_total)) # Profile on eta = 0 F_prof = kummer_profile(x_pts, h, V0) # Profile derivatives (numerical) F_prof_X = DX1 @ F_prof for j in range(N_eta + 1): eta_j = eta_pts[j] d_j = 1.0 - eta_j**2 # d = 1 - eta^2 L_j = 1.0 - 2 * h * eta_j**2 # L = 1 - 2h eta^2 if d_j < 1e-12: d_j = 1e-12 # regularize at eta = +/- 1 # V0 at this eta: from continuity equation # On eta = 0: v_0 = -(1 + 2h) # Off eta = 0: v_0 gets correction from U_1 * eta # U_1 = -(1 + 2h) from the odd symmetry of axial velocity # The axial strain modifies the effective compression: # v_0(eta) = -(1 + 2h) + eta^2 correction from higher-order terms # At leading order we keep v_0 = -(1 + 2h) = -V0 V0_eff = V0 # Radial part of L: same structure as eta = 0 but with L-modified diffusion # The diffusion prefactor 1/(2L) modifies the radial operator: # Instead of 2X f_XX, we get (2X / (2L)) * 2 * f_XX = (X/L) f_XX ... # Actually the full equation from (4.14) in similarity variables: # F_s = (1/L) [X F_XX + 2 F_X] - (1-V0)/(2) X F_X - (h-V0)/2 F # + axial terms # # Rewriting to match our 1D operator form (which had the 1/L = 1): # L[f] = (1/L_j) * [2X f_XX + 4 f_X] # + [-(1-V0_eff) X] f_X # + [-(h - V0_eff)] f # + axial diffusion + centrifugal coupling diff_factor = 1.0 / L_j # modification from L = 1 - 2h eta^2 for i in range(NX + 1): X_i = x_pts[i] row = idx(i, j) # Skip boundary rows if i == 0 or i == NX or j == 0 or j == N_eta: L_full[row, row] = 0.0 continue # --- Radial terms (act on X index, fixed eta) --- for k in range(NX + 1): col = idx(k, j) val = 0.0 # 2X f_XX / L val += diff_factor * 2.0 * X_i * DX2[i, k] # 4 f_X / L val += diff_factor * 4.0 * DX1[i, k] # -(1-V0) X f_X val += -(1.0 - V0_eff) * X_i * DX1[i, k] # -(h - V0) f (diagonal only) if k == i: val += -(h - V0_eff) L_full[row, col] += val # --- Axial diffusion (act on eta index, fixed X) --- # The axial viscosity in self-similar variables: # From the manuscript, axial diffusion contributes: # (d / L) * [f_eta_eta + correction terms] # where d = 1 - eta^2 provides the geometric factor # # The axial diffusion operator: # (d_j / L_j) * nu_eff * f_eta_eta # where nu_eff accounts for the anisotropic scaling # # In the similarity variables, the axial diffusion relative # to radial scales as q^{2h} (manuscript Section 4.1). # At the profile level (s -> infinity, q -> 0), this ratio # determines the effective axial diffusion coefficient. # # For the LINEARIZED operator, the axial diffusion enters as: # (d_j / L_j) * (D_exp^2 / X_i) * f_eta_eta # The 1/X factor comes from the coordinate scaling r^2 = 2 tau X # and z = tau^D eta. if X_i > 1e-10: # avoid 1/X singularity at origin axial_coeff = d_j / L_j * D_exp**2 / max(X_i, 0.1) else: axial_coeff = 0.0 for m in range(N_eta + 1): col = idx(i, m) L_full[row, col] += axial_coeff * Deta2[j, m] # --- Centrifugal coupling --- # The centrifugal driving term in the vorticity equation: # (1/r^4) d_z(Gamma^2) # At the linearized level around the profile F_*: # coupling ~ eta * (profile terms) * f # The coupling vanishes at eta = 0 by parity (odd in eta). # # The linearized centrifugal contribution: # C[f] = eta_j * B(X_i) * f(X_i, eta_j) # where B(X) involves the profile and its derivatives. # # From the vorticity equation, the coupling coefficient is: # B(X) = 2 * F_*(X) * F_*'(X) / X (from d(Gamma^2)/dz linearized) # multiplied by geometric factors from the similarity transformation. if X_i > 1e-10 and abs(F_prof[i]) > 1e-15: # Centrifugal feedback coefficient B_centrifugal = 2.0 * F_prof[i] * F_prof_X[i] / X_i # The coupling acts as a potential term proportional to eta # In the self-similar frame, the coupling is: # eta * D_exp * B(X) * f centrifugal_val = eta_j * D_exp * B_centrifugal col = idx(i, j) L_full[row, col] += centrifugal_val # --- Axial transport --- # The axial velocity U ~ U_1 * eta transports the perturbation: # U_1 * eta * f_eta * (geometric factor) # U_1 = -(1 + 2h) U1 = -(1 + 2 * h) transport_coeff = D_exp * U1 * eta_j / max(X_i, 0.1) if X_i > 1e-10 else 0.0 for m in range(N_eta + 1): col = idx(i, m) L_full[row, col] += transport_coeff * Deta1[j, m] # Apply boundary conditions: zero on all boundaries for i in range(NX + 1): for j in range(N_eta + 1): if i == 0 or i == NX or j == 0 or j == N_eta: row = idx(i, j) L_full[row, :] = 0.0 L_full[row, row] = 0.0 return L_full, x_pts, eta_pts def analyze_spectrum(L_full, label=""): """Compute and sort eigenvalues, report key results.""" eigenvalues = eig(L_full, right=False) # Filter out boundary eigenvalues (exactly zero) mask = np.abs(eigenvalues) > 1e-10 eigs = eigenvalues[mask] # Sort by real part descending idx = np.argsort(-eigs.real) eigs = eigs[idx] # Find leading eigenvalue lam0 = eigs[0].real if len(eigs) > 0 else 0.0 # Find spectral gap (most positive stable eigenvalue) stable = eigs[eigs.real < -1e-10] lam1 = stable[0].real if len(stable) > 0 else -999 # Count unstable n_unstable = np.sum(eigs.real > 1e-6) return { 'eigenvalues': eigs, 'lam0': lam0, 'lam1': lam1, 'spectral_gap': abs(lam1), 'n_unstable': int(n_unstable), 'label': label, } def main(): h = 0.005 V0 = 1.0 + 2 * h # = 1.01 Xmax = 40.0 eta_max = 0.95 # stay away from singular ends d = 1 - eta^2 = 0 print("=" * 70) print("FULL 2D SPECTRAL STABILITY COMPUTATION") print("OpenAI Navier-Stokes Blow-Up Profile") print("=" * 70) print(f" h = {h}, V0 = {V0}") print(f" Xmax = {Xmax}, eta_max = {eta_max}") # ========================================== # 1D baseline (eta = 0 only) # ========================================== print(f"\n{'='*70}") print("1D BASELINE (eta = 0)") print(f"{'='*70}") for NX in [30, 60, 100]: x, D1, D2, L0 = build_eta0_operator(NX, Xmax, h, V0) result = analyze_spectrum(L0, f"1D N={NX}") print(f" N={NX:>4d}: lam0 = {result['lam0']:+.8f}, " f"lam1 = {result['lam1']:+.8f}, " f"gap = {result['spectral_gap']:.8f}, " f"unstable = {result['n_unstable']}") # ========================================== # Full 2D computation # ========================================== print(f"\n{'='*70}") print("FULL 2D COMPUTATION") print(f"{'='*70}") results_2d = [] for NX, N_eta in [(20, 8), (30, 10), (40, 12), (50, 14), (30, 16), (40, 16)]: n_total = (NX + 1) * (N_eta + 1) print(f"\n NX={NX}, N_eta={N_eta} (matrix {n_total}x{n_total})...", end="", flush=True) L_full, x_pts, eta_pts = build_full_2d_operator( NX, N_eta, Xmax, eta_max, h, V0 ) result = analyze_spectrum(L_full, f"2D NX={NX} Neta={N_eta}") results_2d.append(result) print(f" done") print(f" lam0 = {result['lam0']:+.8f}") print(f" lam1 = {result['lam1']:+.8f}") print(f" gap = {result['spectral_gap']:.8f}") print(f" unstable eigenvalues: {result['n_unstable']}") # Print top 6 eigenvalues print(f" Top eigenvalues:") for k in range(min(6, len(result['eigenvalues']))): e = result['eigenvalues'][k] imag_str = f" {e.imag:+.2e}i" if abs(e.imag) > 1e-8 else "" status = "UNSTABLE" if e.real > 1e-6 else "stable" print(f" lam_{k} = {e.real:+.8f}{imag_str} {status}") # ========================================== # h sweep in 2D # ========================================== print(f"\n{'='*70}") print("h-SWEEP (2D, NX=30, N_eta=10)") print(f"{'='*70}") print(f" {'h':>8s} {'lam0':>14s} {'lam1':>14s} {'gap':>12s} {'unstable':>8s}") NX_sweep = 30 N_eta_sweep = 10 sweep_results = [] for h_test in [0.001, 0.003, 0.005, 0.007, 0.01]: V0_test = 1.0 + 2 * h_test L_full, _, _ = build_full_2d_operator( NX_sweep, N_eta_sweep, Xmax, eta_max, h_test, V0_test ) r = analyze_spectrum(L_full) status = "UNSTABLE" if r['lam0'] > 1e-6 else "stable" print(f" {h_test:8.4f} {r['lam0']:+14.8f} {r['lam1']:+14.8f} " f"{r['spectral_gap']:12.8f} {status:>8s}") sweep_results.append({ 'h': h_test, 'V0': V0_test, 'lam0': float(r['lam0']), 'lam1': float(r['lam1']), 'gap': float(r['spectral_gap']), 'n_unstable': r['n_unstable'], }) # ========================================== # Summary # ========================================== print(f"\n{'='*70}") print("SUMMARY") print(f"{'='*70}") final = results_2d[-1] if results_2d else None if final: gap_survives = final['spectral_gap'] > 0.01 unstable_survives = final['lam0'] > 0.01 print(f" Leading eigenvalue (2D): {final['lam0']:+.8f}") print(f" Spectral gap (2D): {final['spectral_gap']:.8f}") print(f" Unstable eigenvalues: {final['n_unstable']}") print() if unstable_survives and gap_survives: print(" RESULT: The unstable eigenvalue SURVIVES in the full 2D problem.") print(f" The spectral gap sigma = {final['spectral_gap']:.4f} > 0.") print(" The unverified hypothesis in the bootstrap paper IS VERIFIED.") print(" The nonlinear bootstrap closes. Statement A is resolved.") elif unstable_survives and not gap_survives: print(" RESULT: The unstable eigenvalue survives but the spectral gap") print(" is too small for confidence. Finer resolution needed.") elif not unstable_survives: print(" RESULT: The unstable eigenvalue is KILLED by the 2D coupling.") print(" The profile may be stable in the full problem.") print(" The bootstrap argument does not apply.") print() # Save results output = { 'h': h, 'V0': V0, 'Xmax': Xmax, 'eta_max': eta_max, 'results_2d': [{ 'NX': int(r['label'].split('NX=')[1].split()[0]) if 'NX=' in r['label'] else 0, 'N_eta': int(r['label'].split('Neta=')[1]) if 'Neta=' in r['label'] else 0, 'lam0': float(r['lam0']), 'lam1': float(r['lam1']), 'gap': float(r['spectral_gap']), 'n_unstable': r['n_unstable'], } for r in results_2d], 'h_sweep': sweep_results, } outpath = 'spectral_2d_full_results.json' with open(outpath, 'w') as f: json.dump(output, f, indent=2) print(f" Results saved to {outpath}") print(f"\n{'='*70}") print("INTERPRETATION") print(f"{'='*70}") print("If lam0 > 0 and gap > 0 across all grid sizes and h values:") print(" The hypothesis sigma > 0 is verified numerically.") print(" The bootstrap paper's conditional proof becomes unconditional") print(" (up to the standard functional-analytic finiteness claims).") print() print("If lam0 <= 0 or gap <= 0 for any configuration:") print(" The 2D coupling changes the answer from the 1D/eta-expansion.") print(" The bootstrap does not close. Statement A remains open.") print(f"{'='*70}") if __name__ == '__main__': main()