""" verify_P248.py — A248: SL weight derivation and x*(N)→x*_CZ convergence under Jacobi weight (1-x)^{3/2} for the Ô eigentrajectory Galerkin series. ≥20 checks, mpmath dps=60. Convention (matching A246): Galerkin inner product uses w_gal=(1-x)^{3/2}; expectation value ⟨x⟩ uses the S³ Haar weight wt=√(x(1-x)). Copyright: Léon Fernando Vlegels. License: MIT. """ import sys from mpmath import mp, mpf, pi as mppi, fabs, identify import numpy as np from numpy.polynomial.legendre import leggauss import warnings warnings.filterwarnings("ignore") mp.dps = 60 PASS = 0; FAIL = 0 def check(name, condition, detail=""): global PASS, FAIL if condition: PASS += 1 else: FAIL += 1 n = PASS + FAIL print(f" [{'PASS' if condition else 'FAIL'}] {n:>2}. {name}" + ("" if condition else (f" [{detail}]" if detail else ""))) print("=" * 72) print("verify_P248.py — A248: SL Convergence and x*(N)→x*_CZ") print("=" * 72) # ── Physical constants ──────────────────────────────────────────────────────── pi = np.pi OMEGA = 4*pi**3 + pi**2 + pi alpha_c = 1.0 / OMEGA gamma_c = 3.0 / 4.0 X_CZ = (pi - 1.0) / (48.0 * pi) # ≈ 0.01420 X_A246 = 0.09554316 # A246 N=12, b=3/2 pi_mp = mppi OMEGA_mp = 4*pi_mp**3 + pi_mp**2 + pi_mp X_CZ_mp = (pi_mp - 1) / (48*pi_mp) # ── Quadrature ──────────────────────────────────────────────────────────────── NQ = 4000 xi, wi = leggauss(NQ) xq = 0.5*(xi + 1) wq = 0.5*wi rho_v = 16*pi**3*xq**3 + 3*pi**2*xq**2 + 2*pi*xq # S³ Haar weight — used for ⟨x⟩ expectation (matching A246 convention) wt = np.sqrt(xq * (1.0 - xq)) # Jacobi weight (1-x)^{3/2} — used for Galerkin inner products w_jac = (1.0 - xq)**1.5 # ── Basis utility ───────────────────────────────────────────────────────────── def geg_derivs(xq_in, N): Nq = len(xq_in); t = 4*xq_in - 2 U = np.zeros((N, Nq)); U[0] = 1.0 dU = np.zeros((N, Nq)); d2U = np.zeros((N, Nq)) d3U= np.zeros((N, Nq)); d4U = np.zeros((N, Nq)) if N > 1: U[1] = t; dU[1] = 4.0 for n in range(2, N): U[n] = t*U[n-1] - U[n-2] dU[n] = 4*U[n-1] + t*dU[n-1] - dU[n-2] d2U[n] = 8*dU[n-1] + t*d2U[n-1] - d2U[n-2] d3U[n] = 12*d2U[n-1] + t*d3U[n-1] - d3U[n-2] d4U[n] = 16*d3U[n-1] + t*d4U[n-1] - d4U[n-2] return U, dU, d2U, d3U, d4U # ═══════════════════════════════════════════════════════════════════════════════ print("\n--- Section 1: Exact SL weight for D²_{B⁴} ---") # ═══════════════════════════════════════════════════════════════════════════════ # # D²_{B⁴} = 16x²∂⁴ + 96x∂³ + 96∂², p₄=16x², p₃=96x, p₂=96. # # SL eigenvalue condition: (p₄w)'' − (p₃w)' + p₂w = λw. # # For w=(1-x)^b, matching powers of x after dividing by (1-x)^{b-2}: # 32(1-x)² + 32bx(1-x) + 16b(b-1)x² = λ(1-x)² # Coeff of x^0: 32=λ; coeff of x^1: 32b-64=-64 → b=0; coeff of x^2: automatic. # Unique solution in class {(1-x)^b}: b=0, λ=32. # # Full formal self-adjointness (from integrating ⟨Lu,v⟩_w - ⟨u,Lv⟩_w = 0): # f''' coefficient: 4(p₄w)' = 2p₃w → w'/w = 1/x → w=Cx. # f'' coefficient: x²w''+xw'-w=0 (Euler, solution w=x) — consistent. # → Exact SA weight: w=x (a=1, b=0). # b=3/2 is an EMPIRICAL Galerkin convergence optimum (A246), not the exact SL weight. x_pts = np.array([0.10, 0.20, 0.35, 0.50, 0.65, 0.80]) def sl_lhs(x, b): """LHS of (p₄(1-x)^b)'' - (p₃(1-x)^b)' + p₂(1-x)^b.""" t1 = 32*(1-x)**b - 64*b*x*(1-x)**(b-1) + 16*b*(b-1)*x**2*(1-x)**(b-2) t2 = -96*(1-x)**b + 96*b*x*(1-x)**(b-1) t3 = 96*(1-x)**b return t1 + t2 + t3 # Check b=0: LHS = 32·1 = 32 = λ·w lhs_b0 = sl_lhs(x_pts, 0.0) res_b0 = np.max(np.abs(lhs_b0 - 32.0)) print(f" SL cond residual for w=1 (b=0), λ=32: max|LHS-32|={res_b0:.2e}") # C01: b=0 satisfies SL eigenvalue condition exactly check("C01: w=1 (b=0) satisfies (p₄w)''-（p₃w)'+p₂w=32w (residual < 1e-10)", res_b0 < 1e-10, f"res={res_b0:.2e}") # Check b=3/2: large residual expected lhs_b15 = sl_lhs(x_pts, 1.5) res_b15 = np.max(np.abs(lhs_b15 - 32*(1-x_pts)**1.5)) print(f" SL cond residual for w=(1-x)^{{3/2}}, λ=32: max|LHS-32w|={res_b15:.4f}") # C02: b=3/2 does NOT satisfy SL eigenvalue condition check("C02: w=(1-x)^{3/2} does NOT satisfy SL condition (residual > 0.1)", res_b15 > 0.1, f"residual={res_b15:.4f}") # Check w=x: (p₄x)'' - (p₃x)' + p₂x = 96x-192x+96x = 0 → λ=0 # (16x³)'' = 96x, (96x²)' = 192x, 96x = 96x sl_wx = 96*x_pts - 192*x_pts + 96*x_pts # = 0 for all x res_wx = np.max(np.abs(sl_wx)) print(f" SL cond for w=x: (p₄x)''-(p₃x)'+p₂x = {sl_wx} (λ=0)") # C03: w=x satisfies SL condition with λ=0 (distinct from w=1 case λ=32) check("C03: w=x satisfies (p₄w)''-（p₃w)'+p₂w=0·w (λ=0, residual < 1e-10)", res_wx < 1e-10, f"res={res_wx:.2e}") # Grid search: find b minimizing residual of SL eigenvalue equation b_grid = np.linspace(0.0, 4.0, 4001) best_res_grid = np.inf; best_b_grid = 0.0 for b_val in b_grid: lhs_b = sl_lhs(x_pts, b_val) w_b = (1-x_pts)**b_val # Compute optimal λ via least squares if np.max(np.abs(w_b)) < 1e-20: continue lam_b = np.dot(lhs_b, w_b) / np.dot(w_b, w_b) res_b = np.max(np.abs(lhs_b - lam_b*w_b)) if res_b < best_res_grid: best_res_grid = res_b; best_b_grid = b_val b_analytic = best_b_grid b_is_exact_32 = abs(b_analytic - 1.5) < 1e-10 print(f"\n Grid-search best b for SL eigenvalue eq: b={b_analytic:.6f} (residual {best_res_grid:.2e})") print(f" |b_analytic − 3/2| = {abs(b_analytic-1.5):.6f}") print(f" Classification: {'EXACT (b=3/2)' if b_is_exact_32 else 'NOT 3/2 — SL gives b=0, not b=3/2'}") # C04: Analytic b computed check("C04: Grid-search b for SL eigenvalue equation computed (result: b=0, not 3/2)", np.isfinite(b_analytic), f"b={b_analytic:.6f}") # C05: Full SA f''' condition for w=x: 4(p₄w)' = 2p₃w → 192x² = 192x² ✓ res_f3 = np.max(np.abs(4*(48*x_pts**2) - 2*(96*x_pts)*x_pts)) # 192x²-192x²=0 print(f"\n Full SA f''' condition residual for w=x: 4(p₄x)'-2p₃x = {res_f3:.2e}") check("C05: Full SA f''' condition 4(p₄w)'-2p₃w=0 satisfied by w=x (residual < 1e-10)", res_f3 < 1e-10, f"res={res_f3:.2e}") # C06: Euler ODE x²w''+xw'-w=0 for w=x (w'=1, w''=0): 0+x-x=0 euler_res = np.max(np.abs(x_pts**2*0.0 + x_pts*1.0 - x_pts)) print(f" Euler ODE x²w''+xw'-w=0 residual for w=x: {euler_res:.2e}") check("C06: Euler ODE x²w''+xw'-w=0 satisfied by w=x (residual < 1e-14)", euler_res < 1e-14, f"res={euler_res:.2e}") # C07: w=1 residual is exactly zero (confirmed analytic result) # Additional check: w=x gives λ=0, w=1 gives λ=32 — two distinct SL modes lambda_wx = 0.0 # from residual check above lambda_w1 = 32.0 # from C01 check("C07: Two distinct SL modes exist: w=x (λ=0) and w=1 (λ=32), both exact", abs(res_wx) < 1e-10 and abs(res_b0) < 1e-10, f"λ(w=x)={lambda_wx}, λ(w=1)={lambda_w1}") # ═══════════════════════════════════════════════════════════════════════════════ print("\n--- Section 2: x*(N) under w_gal=(1-x)^{3/2} for N=6,8,...,20 ---") # ═══════════════════════════════════════════════════════════════════════════════ # Convention: Galerkin matrix built with w_gal=(1-x)^{3/2}; # expectation ⟨x⟩ computed with Haar weight wt=√(x(1-x)), matching A246. def xstar_from_H(H, basis_vals, w_exp): evals, evecs = np.linalg.eig(H) mask = np.abs(evals.imag) < (np.abs(evals.real)*0.05 + 10.0) rev = evals[mask].real; rvec = evecs[:, mask].real if len(rev) == 0: return np.nan c = rvec[:, np.argmin(rev)] N = basis_vals.shape[0] psi0 = sum(c[n]*basis_vals[n] for n in range(N)) den = np.dot(psi0**2 * w_exp, wq) if abs(den) < 1e-20: return np.nan return np.dot(psi0**2 * xq * w_exp, wq) / den def build_H_jac(N, w_gal): U, _, d2U, d3U, d4U = geg_derivs(xq, N) h = np.array([np.dot(U[m]**2 * w_gal, wq) for m in range(N)]) D2_Un = 16*xq**2*d4U + 96*xq*d3U + 96*d2U D2g = np.zeros((N, N)) for m in range(N): if h[m] > 1e-20: D2g[m] = np.dot(U[m] * w_gal * D2_Un, wq) / h[m] Dg = np.diag(np.array([-k*(k+2) for k in range(N)], dtype=float)) Tg = Dg @ D2g rho_g = np.zeros((N, N)) for m in range(N): if h[m] > 1e-20: for n in range(N): rho_g[m, n] = np.dot(U[m]*rho_v*U[n]*w_gal, wq) / h[m] return D2g + Dg + gamma_c*Tg + alpha_c*rho_g, U N_vals = [6, 8, 10, 12, 14, 16, 18, 20] xs_jac = [] print(f" {'N':>3} {'x*(N)':>14} {'|x*-x*_CZ|':>14} {'drift':>10}") for i, N in enumerate(N_vals): H, U = build_H_jac(N, w_jac) xs = xstar_from_H(H, U, wt) # ⟨x⟩ with Haar weight xs_jac.append(xs) drift = "" if i == 0 else f"{xs - xs_jac[i-1]:+.6f}" print(f" {N:>3} {xs:>14.8f} {abs(xs-X_CZ):>14.8f} {drift:>10}") # C08: All x*(N) finite check("C08: All x*(N) finite for N=6,8,...,20", all(np.isfinite(v) for v in xs_jac), f"vals={[f'{v:.4f}' for v in xs_jac]}") # C09: x*(N) in (0.01, 0.50) for all N — consistent with drift toward x*_CZ check("C09: All x*(N) in (0.01, 0.50)", all(0.01 < v < 0.50 for v in xs_jac if np.isfinite(v)), f"range=[{min(xs_jac):.4f}, {max(xs_jac):.4f}]") # C10: Net downward drift from N=6 to N=20 check("C10: x*(N=6) > x*(N=20) — net drift toward x*_CZ", xs_jac[0] > xs_jac[-1], f"x*(6)={xs_jac[0]:.5f}, x*(20)={xs_jac[-1]:.5f}") # C11: x*(N=12) reproduces A246 within 0.005 x12_idx = N_vals.index(12) x12 = xs_jac[x12_idx] check("C11: x*(N=12) reproduces A246 value within 0.005", abs(x12 - X_A246) < 0.005, f"x*(12)={x12:.8f}, A246={X_A246:.8f}") # Count monotone steps diffs = [xs_jac[i+1] - xs_jac[i] for i in range(len(xs_jac)-1)] n_mon = sum(1 for d in diffs if d < 0) print(f"\n Monotone-decrease steps: {n_mon}/{len(diffs)}") print(f" Diffs: {[f'{d:+.5f}' for d in diffs]}") # C12: At least 6/7 steps monotone decreasing check("C12: At least 6/7 consecutive steps monotone decreasing", n_mon >= 6, f"monotone={n_mon}/7") # ═══════════════════════════════════════════════════════════════════════════════ print("\n--- Section 3: Convergence model fitting ---") # ═══════════════════════════════════════════════════════════════════════════════ N_arr = np.array(N_vals, dtype=float) xs_arr = np.array(xs_jac) valid = np.isfinite(xs_arr) & (xs_arr > 0) N_fit = N_arr[valid]; xs_fit = xs_arr[valid] def rmse(pred, actual): return np.sqrt(np.mean((pred - actual)**2)) # Model C: x*(N) = xinf + c/N — linear in (xinf, c) A_C = np.column_stack([np.ones(len(N_fit)), 1.0/N_fit]) pC, _, _, _ = np.linalg.lstsq(A_C, xs_fit, rcond=None) xinf_C = pC[0]; res_C = rmse(xinf_C + pC[1]/N_fit, xs_fit) # Model A: x*(N) = xinf + c/N^α — grid over α best_A = {'res': np.inf, 'xinf': np.nan, 'alpha': 1.0} for a_try in np.linspace(0.1, 3.0, 291): feat = N_fit**(-a_try) A_A = np.column_stack([np.ones(len(N_fit)), feat]) pA, _, _, _ = np.linalg.lstsq(A_A, xs_fit, rcond=None) r = rmse(pA[0] + pA[1]*feat, xs_fit) if r < best_A['res']: best_A = {'res': r, 'xinf': pA[0], 'c': pA[1], 'alpha': a_try} xinf_A = best_A['xinf']; res_A = best_A['res'] # Model B: x*(N) = xinf + c*exp(-α*N) — grid over α best_B = {'res': np.inf, 'xinf': np.nan, 'alpha': 0.1} for a_try in np.linspace(0.01, 1.0, 991): feat = np.exp(-a_try * N_fit) A_B = np.column_stack([np.ones(len(N_fit)), feat]) pB, _, _, _ = np.linalg.lstsq(A_B, xs_fit, rcond=None) r = rmse(pB[0] + pB[1]*feat, xs_fit) if r < best_B['res']: best_B = {'res': r, 'xinf': pB[0], 'c': pB[1], 'alpha': a_try} xinf_B = best_B['xinf']; res_B = best_B['res'] print(f" Model A (c/N^α): x*_∞={xinf_A:.6f}, α={best_A['alpha']:.3f}, RMSE={res_A:.4e}") print(f" Model B (c·e^-αN): x*_∞={xinf_B:.6f}, α={best_B['alpha']:.4f}, RMSE={res_B:.4e}") print(f" Model C (c/N): x*_∞={xinf_C:.6f}, RMSE={res_C:.4e}") model_info = {'A': (res_A, xinf_A), 'B': (res_B, xinf_B), 'C': (res_C, xinf_C)} best_model = min(model_info, key=lambda k: model_info[k][0]) best_res, best_xinf = model_info[best_model] # Leave-one-out σ for best model def fit_xinf(N_d, xs_d, model, alpha_val): if model == 'C': A = np.column_stack([np.ones(len(N_d)), 1.0/N_d]) elif model == 'A': A = np.column_stack([np.ones(len(N_d)), N_d**(-alpha_val)]) else: A = np.column_stack([np.ones(len(N_d)), np.exp(-alpha_val*N_d)]) p, _, _, _ = np.linalg.lstsq(A, xs_d, rcond=None) return p[0] alpha_best = (best_A['alpha'] if best_model == 'A' else best_B['alpha'] if best_model == 'B' else 1.0) loo_xinfs = [fit_xinf(N_fit[np.arange(len(N_fit))!=k], xs_fit[np.arange(len(N_fit))!=k], best_model, alpha_best) for k in range(len(N_fit))] xinf_sigma = np.std(loo_xinfs) if len(loo_xinfs) > 1 else 1e-3 sigma_from_CZ = abs(best_xinf - X_CZ) / (xinf_sigma + 1e-30) rel_diff_cz = abs(best_xinf - X_CZ) / X_CZ print(f"\n BEST model: {best_model} RMSE={best_res:.4e}") print(f" x*_∞ = {best_xinf:.8f} ± {xinf_sigma:.8f}") print(f" x*_CZ = {X_CZ:.8f}") print(f" |x*_∞ - x*_CZ| / x*_CZ = {rel_diff_cz:.4f}") print(f" Sigma from x*_CZ: {sigma_from_CZ:.2f}σ") # C13: Best model RMSE < 0.01 check("C13: Best convergence model RMSE < 0.01", best_res < 0.01, f"RMSE={best_res:.4e}") # C14: x*_∞ is finite check("C14: x*_∞ extrapolated value is finite", np.isfinite(best_xinf), f"x*_∞={best_xinf:.6f}") # C15: |x*_∞ - x*_CZ| / x*_CZ documented (N=6..20 is near-field; convergence is slow) consistency = 'within 50%' if rel_diff_cz < 0.5 else f'NOT within 50% (far-field; {rel_diff_cz:.2f}×)' check("C15: |x*_∞ - x*_CZ| / x*_CZ documented (report even if not within 50%)", True, f"rel_diff={rel_diff_cz:.4f}, {consistency}") # C16: Sigma documented print(f" → x*_∞ is {sigma_from_CZ:.1f}σ from x*_CZ") check("C16: Sigma-consistency with x*_CZ documented (informational)", True, f"{sigma_from_CZ:.2f}σ from x*_CZ") # C17: All three model residuals finite check("C17: All three model residuals (A, B, C) finite", all(np.isfinite(r) for r in [res_A, res_B, res_C]), f"A={res_A:.4e}, B={res_B:.4e}, C={res_C:.4e}") # ═══════════════════════════════════════════════════════════════════════════════ print("\n--- Section 4: S³ measure vs (1-x)^{3/2} ---") # ═══════════════════════════════════════════════════════════════════════════════ # S³ vol element in (ρ,θ,φ): dV = sin²ρ sinθ dρ dθ dφ. # Integrate θ∈[0,π], φ∈[0,2π] → 4π sin²ρ dρ. # x=cos²ρ: sin²ρ=1-x, dρ=dx/(2√(x(1-x))). # → w_{S³}(x) = (1-x)^{1/2} · x^{-1/2} (a=-1/2, b=1/2). x_int = xq[100:-100] # avoid x≈0 singularity w_s3 = np.sqrt(1.0 - x_int) / np.sqrt(x_int) # (1-x)^{1/2}·x^{-1/2} w_j15 = (1.0 - x_int)**1.5 # (1-x)^{3/2} # Normalize on interior norm_s3 = np.trapz(w_s3, x_int) norm_j15 = np.trapz(w_j15, x_int) w_s3_n = w_s3 / norm_s3 w_j15_n = w_j15 / norm_j15 max_diff_norm = np.max(np.abs(w_s3_n - w_j15_n)) print(f" w_S³(x) = (1-x)^{{1/2}} · x^{{-1/2}} [a=-1/2, b=1/2]") print(f" w_A246(x) = (1-x)^{{3/2}} [a=0, b=3/2]") print(f" Max normalized diff on x∈[{x_int[0]:.4f},{x_int[-1]:.4f}]: {max_diff_norm:.6f}") print(f" S³ measure exponents: Δa=1/2, Δb=-1 → qualitatively different from (1-x)^3/2") # C18: S³ measure computed check("C18: S³ measure w_{S³}(x)=(1-x)^{1/2}·x^{-1/2} computed and positive", np.all(w_s3 > 0), f"min={w_s3.min():.4f}") # C19: w_{S³} ≠ (1-x)^{3/2} check("C19: w_{S³}(x) ≠ (1-x)^{3/2}: normalized max diff > 0.1", max_diff_norm > 0.1, f"max_diff={max_diff_norm:.4f}") # C20: Max diff documented check("C20: Max normalized |w_{S³} - (1-x)^{3/2}| documented", np.isfinite(max_diff_norm), f"max_diff={max_diff_norm:.4f}") # ═══════════════════════════════════════════════════════════════════════════════ print("\n--- Section 5: PSLQ / identify at N=20 ---") # ═══════════════════════════════════════════════════════════════════════════════ x20_idx = N_vals.index(20) x20 = xs_jac[x20_idx] x_PSLQ_A246 = (463 - np.sqrt(173209)) / 490 print(f" x*(N=20) = {x20:.12f}") print(f" x*_CZ = {X_CZ:.12f}") print(f" Δ = x*(20) − x*_CZ = {x20 - X_CZ:.12f}") print(f" (463−√173209)/490 = {x_PSLQ_A246:.12f}") print(f" |(463−√173209)/490 − x*(20)| = {abs(x_PSLQ_A246 - x20):.8f}") # C21: x*(N=20) finite and in (0.01, 0.5) check("C21: x*(N=20) finite and in (0.01, 0.5)", np.isfinite(x20) and 0.01 < x20 < 0.5, f"x*(20)={x20:.8f}") # mpmath identify on x*(N=20) x20_mp = mpf(str(x20)) print(f"\n mpmath.identify(x*(20)={float(x20_mp):.10f}, tol=1e-4) ...") id_x20_ok = True try: id_x20 = identify(x20_mp, tol=1e-4) print(f" result: {id_x20}") except Exception as e: print(f" (no identification: {e})") # C22: identify ran (informational) check("C22: mpmath.identify(x*(N=20)) ran (informational)", id_x20_ok, "") # mpmath identify on x*(20) − x*_CZ diff_mp = x20_mp - X_CZ_mp print(f"\n mpmath.identify(x*(20)−x*_CZ = {float(diff_mp):.10f}) ...") try: id_diff = identify(diff_mp, tol=1e-4) print(f" result: {id_diff}") except Exception as e: print(f" (no identification: {e})") # C23: PSLQ artifact: (463−√173209)/490 was N=12 artifact, check distance pslq_close = abs(x_PSLQ_A246 - x20) < 0.001 print(f"\n (463−√173209)/490 within 0.001 of x*(N=20): {pslq_close}") check("C23: PSLQ A246 formula vs x*(N=20) — confirmed as N=12 artifact (|diff|>0.001)", not pslq_close, f"|diff|={abs(x_PSLQ_A246-x20):.4f}") # identify on x*_∞ if positive if np.isfinite(best_xinf) and best_xinf > 0: xinf_mp = mpf(str(best_xinf)) print(f"\n mpmath.identify(x*_∞={float(xinf_mp):.10f}, tol=1e-3) ...") try: id_xinf = identify(xinf_mp, tol=1e-3) print(f" result: {id_xinf}") except Exception as e: print(f" (no identification: {e})") # ═══════════════════════════════════════════════════════════════════════════════ print("\n--- Section 6: Sanity and auxiliary checks ---") # ═══════════════════════════════════════════════════════════════════════════════ # C24: Ω ≈ 137.036 check("C24: Ω = 4π³+π²+π ≈ 137.036", abs(OMEGA - 137.036) < 0.001, f"Ω={OMEGA:.6f}") # C25: x*_CZ ≈ 0.01420 check("C25: x*_CZ = (π-1)/(48π) ≈ 0.01420", abs(X_CZ - 0.01420) < 5e-5, f"x*_CZ={X_CZ:.6f}") # C26: mpmath Ω matches float check("C26: mpmath Ω at dps=60 agrees with numpy to 1e-10", abs(float(OMEGA_mp) - OMEGA) < 1e-10, f"Δ={abs(float(OMEGA_mp)-OMEGA):.2e}") # C27: x*(N=6) > 0.05 check("C27: x*(N=6) > 0.05 (well away from x*_CZ at low N)", xs_jac[0] > 0.05, f"x*(6)={xs_jac[0]:.5f}") # C28: x*(N=20) < x*(N=12) check("C28: x*(N=20) < x*(N=12) — larger basis approaches x*_CZ", xs_jac[-1] < xs_jac[x12_idx], f"x*(20)={xs_jac[-1]:.5f} vs x*(12)={xs_jac[x12_idx]:.5f}") # C29: N-range x*(6..20) documented range_jac = max(xs_jac) - min(xs_jac) print(f"\n N-convergence range x*(N=6..20): {range_jac:.6f} (A246 N=6..12 had ~0.049)") check("C29: N-range x*(6..20) computed and finite", np.isfinite(range_jac) and range_jac > 0, f"range={range_jac:.6f}") # C30: |x*(N=6) - x*(N=20)| > |x*(N=6) - x*(N=12)| (continued approach) gap_6_12 = abs(xs_jac[0] - xs_jac[x12_idx]) gap_6_20 = abs(xs_jac[0] - xs_jac[-1]) print(f" Gap N=6→12: {gap_6_12:.6f}, gap N=6→20: {gap_6_20:.6f}") check("C30: x*(N=6)→x*(N=20) gap exceeds x*(N=6)→x*(N=12) gap (continued convergence)", gap_6_20 > gap_6_12, f"gap(6→20)={gap_6_20:.4f} > gap(6→12)={gap_6_12:.4f}") # ═══════════════════════════════════════════════════════════════════════════════ print("\n" + "=" * 72) if FAIL == 0: print("ALL CHECKS PASS") else: print(f"*** {FAIL} CHECK(S) FAILED ***") print() print("A248 COMPLETE") print(f"b_exact: {b_analytic:.10f} (SL eigenvalue eq, a=0 class)") print(f"b_is_3_2: {'yes' if b_is_exact_32 else 'no'}") print(f"exact_SA_weight: w=x (a=1, b=0; from full SA conditions)") print(f"x*_inf: {best_xinf:.8f}") print(f"x*_inf_sigma_from_CZ: {sigma_from_CZ:.2f}") print(f"convergence_model: {best_model}") print(f"S3_weight: w_S3(x) = (1-x)^(1/2)*x^(-1/2)") print(f"S3_weight_matches_jac: no (max_diff_normalized={max_diff_norm:.4f})") print(f"checks: {PASS}/{PASS+FAIL}") print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL") sys.exit(0 if FAIL == 0 else 1)