validate_branch.py — src/21-gravity

1"""The nonlinear completion (Theorem: no self-sourcing, exponential reading).23Validates the two exact arguments determining the nonlinear completion:45(A) Ward identity / exact Gauss law: in the full SINE-coupled (nonlinear) static6    model on Z^3, the synchronisation flux through every closed surface equals the7    enclosed source demand exactly -- for one source, for two sources, and at8    strong gradients -- and the deviation of the nonlinear potential from the9    linear one scales as the CUBE of the source (field-strength corrections,10    (r_*/r)^4 type), never as the square (no u^2 self-sourcing).1112(B) The reading characterisation (symbolic): isotropic Schwarzschild is the13    exponential reading of the self-sourced potential psi = 2 artanh(U/2)14    (series identity verified), psi is not harmonic (Laplacian ~ (Gm)^3/r^5),15    while the exponential branch reads the harmonic u; and the exponential is16    the unique solution of the composition law R(a+b) = R(a)R(b), which the17    Schwarzschild reading violates.18"""19import numpy as np20import sympy as sp2122ok = True2324# ---------- (A) exact Gauss law in the nonlinear sine model ----------25N = 31                                     # lattice [0,N)^3, sources interior26def solve(sources, sine=True, iters=30000, tol=1e-13):27    """Relax kappa * sum_y sin(u_y - u_x) + m_x = 0 with u = 0 on the boundary."""28    u = np.zeros((N, N, N))29    m = np.zeros((N, N, N))30    for (pos, mass) in sources:31        m[pos] = mass32    for it in range(iters):33        nb = (np.roll(u, 1, 0), np.roll(u, -1, 0), np.roll(u, 1, 1),34              np.roll(u, -1, 1), np.roll(u, 1, 2), np.roll(u, -1, 2))35        if sine:36            drive = sum(np.sin(v - u) for v in nb) + m37        else:38            drive = sum(v - u for v in nb) + m39        drive[0, :, :] = drive[-1, :, :] = 0   # Dirichlet boundary40        drive[:, 0, :] = drive[:, -1, :] = 041        drive[:, :, 0] = drive[:, :, -1] = 042        u += 0.15*drive43        if it % 500 == 0 and np.max(np.abs(drive)) < tol:44            break45    return u, np.max(np.abs(drive))4647def flux_box(u, c, R, sine=True):48    """Outward sine-flux through the surface of the cube of half-size R about c."""49    F = 0.050    lo = [ci - R for ci in c]; hi = [ci + R for ci in c]51    couple = np.sin if sine else (lambda x: x)52    for ax in range(3):53        for face, outw in ((lo[ax], -1), (hi[ax], +1)):54            sl_in = [slice(lo[a], hi[a] + 1) for a in range(3)]55            sl_out = list(sl_in)56            sl_in[ax] = face57            sl_out[ax] = face + outw58            F += np.sum(couple(u[tuple(sl_in)] - u[tuple(sl_out)]))59    return F6061c = (N//2, N//2, N//2)62m0 = 4.0                                    # strong: sin argument ~ 0.55 at first shell63u1, res = solve([(c, m0)])64print(f"  (A) single source m={m0}, solver residual {res:.1e}")65for R in (2, 4, 6, 9):66    F = flux_box(u1, c, R)67    print(f"      flux through box R={R}: {F:.10f}  (demand {m0})")68    if abs(F - m0) > 1e-6: ok = False6970c2 = (N//2 + 5, N//2, N//2)71u2, _ = solve([(c, m0), (c2, 2.0)])72F2 = flux_box(u2, (N//2 + 2, N//2, N//2), 12)73print(f"      two sources, enclosing box: flux = {F2:.10f}  (demand {m0 + 2.0})")74if abs(F2 - (m0 + 2.0)) > 1e-6: ok = False7576# no u^2 self-sourcing: ||u_sine - u_linear|| scales as m^3 (arcsin x = x + x^3/6 + ...)77devs = []78for mm in (1.0, 2.0):79    us, _ = solve([(c, mm)])80    ul, _ = solve([(c, mm)], sine=False)81    devs.append(np.max(np.abs(us - ul)))82p = np.log(devs[1]/devs[0])/np.log(2.0)83print(f"      nonlinear deviation scaling: ||u_sin - u_lin|| ~ m^{p:.2f}  (target 3: cubic, no m^2 term)")84if abs(p - 3) > 0.25: ok = False8586# ---------- (B) symbolic characterisation of the two readings ----------87U = sp.symbols('U', positive=True)88psi = 2*sp.atanh(U/2)89schw = ((1 - U/2)/(1 + U/2))**290ident = sp.simplify(sp.powsimp(sp.expand(sp.exp(-2*psi.rewrite(sp.log)))) - schw)91print(f"  (B) g00(Schwarzschild isotropic) == exp(-2 psi), psi = 2 artanh(U/2): residual {ident}")92if ident != 0: ok = False9394r, Gm = sp.symbols('r Gm', positive=True)95psi_r = (Gm/r) + (Gm/r)**3/12               # psi to the order that decides96lap = sp.simplify(sp.diff(r**2*sp.diff(psi_r, r), r)/r**2)97print(f"      Laplacian(psi) = {lap}  (nonzero: psi is self-sourced; u = Gm/r is harmonic)")98if sp.simplify(lap) == 0: ok = False99100a, b = sp.symbols('a b', positive=True)101expo_ok = sp.simplify(sp.exp(-(a + b)) - sp.exp(-a)*sp.exp(-b)) == 0102f = lambda x: (1 - x/2)/(1 + x/2)103schw_viol = sp.simplify(f(a + b) - f(a)*f(b))104print(f"      composition law: exponential satisfies R(a+b)=R(a)R(b): {expo_ok}; "105      f"Schwarzschild reading violates it (mismatch {sp.expand(schw_viol)} != 0)")106if not expo_ok or sp.simplify(schw_viol) == 0: ok = False107108print("PASS" if ok else "FAIL", "- branch decision (Ward identity + composition law)")109