#!/usr/bin/env python3 """ test_338 -- HYPERSURFACES OF R^4 SWEPT BY ROUND 2-SPHERES: THE METRIC, THE TILT, THE DRIFT, FIRST ORDER A one-parameter family of round 2-spheres in R^4: x(r, omega) = x0(r) + Phi(r) R(r) (omega, 0), R in SO(4), omega in S^2. Write R^T x0' = (a, alpha) and R^T R' = [[rho x, -b], [b^T, 0]]. In coordinates (r, p, q) that follow the orthogonal trajectories of the shells, labelled by a Moebius chart z = p + iq with E = (|z - (u1 + i u2)|^2 + u3^2)/(2 u3) and n the unit normal direction in the shell's 3-plane, the claim is g = [H^2 + (alpha + Phi b.n)^2] dr^2 + Phi^2 E^{-2} (dp^2 + dq^2), H = Phi' + Phi W.n, W = (u1', u2', u3')/u3, with a = Phi W. The tilt b and the out-of-plane speed alpha are free. F = 0 is tested in the polynomial (Thomas) form 2 det(M) Cod(adj M) - d(det M) ^ adj M = 0, in exact rational jets. PRE-REGISTERED: P1 (the construction, numerically, independent of F) For freely chosen Phi, alpha, b(r) and u1, u2, u3 (alpha and b not tied to the chart), integrating R' = R Omega and x0' = R(a, alpha) and differentiating the embedding by finite differences gives g_rr = H^2 + (alpha + Phi b.n)^2, g_pp = g_qq = Phi^2/E^2 and zero cross terms, to 1e-7 relative. P2 (tilt, no drift) With u1 = u2 = 0, u3 = 1 (W = 0) and b != 0 on a non-round base, the metric A^2 = Phi'^2 + (s Phi' + Phi b.n)^2 satisfies F = 0 exactly at two points, with M positive definite. P3 The general swept family (drift, independent tilt, independent alpha) satisfies F = 0 exactly at a point with det M > 0. The subfamily alpha = s Phi', b = s W, s^2 = k/(1 - k), has A^2 = H^2/(1 - k) identically (Szekeres's quasi-spherical form), and satisfies F = 0 exactly on a non-round, non-axisymmetric example. P4 (route 1 is the tilt) On the round base (Phi = sin r, alpha = sin r), the tilt b = (eps/sin^2 r) e_3 gives exactly A = 1 + eps cos(theta) + eps^2 [cot^2 r/6 + (cot^2 r/3) P_2(cos theta)] + O(eps^3). P5 (first order, general alpha) About the shell metric with A0^2 = Phi'^2 + alpha^2, delta A / A0 = Phi (Phi' W + alpha b).n / A0^2; with W = 0 the l = 1 part of delta(A^2) is 2 alpha Phi b.n. P6 The axisymmetric tilt b = (0, 0, r/20), W = 0, on the non-round base of P2 satisfies F = 0 exactly, M > 0. P7 Each shell is intrinsically round: the Gauss curvature of Phi^2 E^{-2}(dp^2 + dq^2) is 1/Phi^2 for every u1, u2, u3 (symbolic). P8 (the drift in shell coordinates) For x = (c(r) e_3 + Phi(r) n(theta, phi), Z(r)) with Phi'^2 + Z'^2 = A^2: g_rr = A^2 + 2 c' Phi' cos(theta) + c'^2, g_{r theta} = -c' Phi sin(theta), and the angular block is Phi^2 times the round metric (symbolic). On the round base with radial modulation 2 eps cos(theta) (c' = eps/cos r): g_{r theta} = eps beta sin(theta) with beta = -tan r. This script replaces test_328 for the paper's verification bundle: the computations of P1-P4 and P6 are test_328's, P5 is generalised from s Phi' to any alpha, and P7-P8 are new. """ import sys import numpy as np import sympy as sp from scipy.integrate import solve_ivp FAILURES = [] def check(name, ok, detail=""): print(f" [{'PASS' if ok else 'FAIL'}] {name}" + (f" -- {detail}" if detail else "")) if not ok: FAILURES.append(name) # ---------------------------------------------------------------- P1: the explicit construction (numerical) print("P1 the swept-sphere construction reproduces A^2 = H^2 + (alpha + Phi b.n)^2, round shells, no cross terms") Phi = lambda r: 1 + r / 2 dPhi = lambda r: 0.5 s_ = lambda r: 0.3 + 0.2 * r alpha_ = lambda r: s_(r) * dPhi(r) + 0.05 * r # NOT the Szekeres value b_ = lambda r: np.array([0.2 * r, -0.1, 0.15 * r * r]) # NOT tied to the drift Pf = lambda r: 0.25 * r Qf = lambda r: r * r / 3 Sf = lambda r: 1 + r / 5 dPf = lambda r: 0.25 dQf = lambda r: 2 * r / 3 dSf = lambda r: 0.2 def sig_inv(w): u, v = w.real, w.imag d = u * u + v * v + 1 return np.array([2 * u, 2 * v, u * u + v * v - 1]) / d def omega(r, z): return sig_inv((z - (Pf(r) + 1j * Qf(r))) / Sf(r)) def domega_dr(r, z, h=1e-6): return (omega(r + h, z) - omega(r - h, z)) / (2 * h) def cross_mat(v): return np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]]) def fit_ckf(r): """d_r omega = beta - (beta.omega) omega + gamma x omega, fitted exactly (a CKF) from sample points.""" rows, rhs = [], [] for z in (0.3 + 0.1j, -0.7 + 0.4j, 0.2 - 0.9j, 1.3 + 0.2j): om = omega(r, z) rows.append(np.hstack([np.eye(3) - np.outer(om, om), -cross_mat(om)])) rhs.append(domega_dr(r, z)) Amat, y = np.vstack(rows), np.hstack(rhs) sol, res, *_ = np.linalg.lstsq(Amat, y, rcond=None) return sol[:3], sol[3:], np.abs(Amat @ sol - y).max() def rhs(r, Y): Rm = Y[:16].reshape(4, 4) beta, gamma, _ = fit_ckf(r) a, rho = -Phi(r) * beta, -gamma b = b_(r) Om = np.zeros((4, 4)) Om[:3, :3] = cross_mat(rho) Om[:3, 3] = -b Om[3, :3] = b dR = Rm @ Om dx0 = Rm @ np.hstack([a, alpha_(r)]) return np.hstack([dR.ravel(), dx0]) r0, r1 = 0.4, 0.5 Y0 = np.hstack([np.eye(4).ravel(), np.zeros(4)]) sol = solve_ivp(rhs, (r0, r1 + 0.01), Y0, rtol=1e-12, atol=1e-14, dense_output=True, method="DOP853") def X(r, p, q): Y = sol.sol(r) Rm, x0 = Y[:16].reshape(4, 4), Y[16:] return x0 + Phi(r) * Rm @ np.hstack([omega(r, p + 1j * q), 0.0]) p0, q0, h = 0.3, -0.2, 1e-5 D = [lambda f: (f(r1 + h, p0, q0) - f(r1 - h, p0, q0)) / (2 * h), lambda f: (f(r1, p0 + h, q0) - f(r1, p0 - h, q0)) / (2 * h), lambda f: (f(r1, p0, q0 + h) - f(r1, p0, q0 - h)) / (2 * h)] J = [d(X) for d in D] g = np.array([[J[i] @ J[j] for j in range(3)] for i in range(3)]) z = p0 + 1j * q0 E = (abs(z - (Pf(r1) + 1j * Qf(r1))) ** 2 + Sf(r1) ** 2) / (2 * Sf(r1)) Eh = lambda r: (abs(z - (Pf(r) + 1j * Qf(r))) ** 2 + Sf(r) ** 2) / (2 * Sf(r)) dE = (Eh(r1 + 1e-6) - Eh(r1 - 1e-6)) / 2e-6 H = dPhi(r1) - Phi(r1) * dE / E nvec = omega(r1, z) grr = H ** 2 + (alpha_(r1) + Phi(r1) * b_(r1) @ nvec) ** 2 want = np.array([[grr, 0, 0], [0, Phi(r1) ** 2 / E ** 2, 0], [0, 0, Phi(r1) ** 2 / E ** 2]]) rel = np.abs(g - want).max() / np.abs(want).max() ckf_res = max(fit_ckf(r)[2] for r in (0.4, 0.45, 0.5)) check("explicit embedding: g = diag(H^2 + (alpha + Phi b.n)^2, Phi^2/E^2, Phi^2/E^2), cross terms 0", rel < 1e-7 and ckf_res < 1e-8, f"max relative deviation {rel:.1e}; chart flow is a CKF to {ckf_res:.1e}; g_rr {g[0,0]:.10f} vs {grr:.10f}") # ---------------------------------------------------------------- exact jets for P2-P3 (the polynomial Thomas form) x, y, zz = sp.symbols("x y z") VV = (x, y, zz) DS = 3 def trunc(e, deg=DS): e = sp.expand(e) if e == 0: return sp.Integer(0) Pl = sp.Poly(e, *VV) return sum((c * x**m[0] * y**m[1] * zz**m[2] for m, c in Pl.terms() if sum(m) <= deg), sp.Integer(0)) def jet_rational(expr, pt, rs, deg=DS + 1): e = sp.together(expr.subs({s: v + d for s, v, d in zip(rs, pt, VV)})) n, d = sp.fraction(e) n = trunc(sp.expand(n), deg) dd = sp.expand(d) D0 = dd.subs({x: 0, y: 0, zz: 0}) Ex = trunc(dd - D0, deg) inv = sp.Integer(1) / D0 term = inv for _ in range(deg + 1): term = trunc(-term * Ex / D0, deg) inv = inv + term return trunc(n * trunc(inv, deg), deg) def matmul(A, B): return sp.Matrix(3, 3, lambda i, j: trunc(sum(A[i, k] * B[k, j] for k in range(3)))) def inv_mat(A): A0 = A.subs({x: 0, y: 0, zz: 0}) A0i = A0.inv() Ex = sp.Matrix(3, 3, lambda i, j: trunc(A[i, j] - A0[i, j])) out, Pm = A0i, A0i for _ in range(DS + 1): Pm = matmul(-A0i * Ex, Pm) out = out + Pm return sp.Matrix(3, 3, lambda i, j: trunc(out[i, j])) dd_ = lambda e, i: sp.diff(e, VV[i]) def det3(A): return trunc(A[0, 0] * trunc(A[1, 1] * A[2, 2] - A[1, 2] * A[2, 1]) - A[0, 1] * trunc(A[1, 0] * A[2, 2] - A[1, 2] * A[2, 0]) + A[0, 2] * trunc(A[1, 0] * A[2, 1] - A[1, 1] * A[2, 0])) def adj3(A): Cf = lambda i, j: trunc(A[(i + 1) % 3, (j + 1) % 3] * A[(i + 2) % 3, (j + 2) % 3] - A[(i + 1) % 3, (j + 2) % 3] * A[(i + 2) % 3, (j + 1) % 3]) return sp.Matrix(3, 3, lambda i, j: Cf(j, i)) def thomas_residual(gm): gi = inv_mat(gm) Gam = [[[trunc(sum(gi[a, e] * (dd_(gm[e, b], c) + dd_(gm[e, c], b) - dd_(gm[b, c], e)) for e in range(3)) / 2) for c in range(3)] for b in range(3)] for a in range(3)] Riem = lambda a, b, c, d: trunc(dd_(Gam[a][d][b], c) - dd_(Gam[a][c][b], d) + sum(Gam[a][c][e] * Gam[e][d][b] - Gam[a][d][e] * Gam[e][c][b] for e in range(3))) Ric = sp.Matrix(3, 3, lambda b, d: trunc(sum(Riem(a, b, a, d) for a in range(3)))) Rs = trunc(sum(gi[b, d] * Ric[b, d] for b in range(3) for d in range(3))) G = sp.Matrix(3, 3, lambda a, b: trunc(Ric[a, b] - Rs * gm[a, b] / 2)) M = sp.Matrix(3, 3, lambda a, b: trunc(-sum(gi[a, c] * G[c, b] for c in range(3)))) detM, Pl = det3(M), matmul(gm, adj3(M)) nab = lambda a, b, c: trunc(dd_(Pl[b, c], a) - sum(Gam[e][a][b] * Pl[e, c] + Gam[e][a][c] * Pl[b, e] for e in range(3))) at0 = lambda e: e.subs({x: 0, y: 0, zz: 0}) res = [at0(2 * detM * (nab(a, b, c) - nab(b, a, c)) - (dd_(detM, a) * Pl[b, c] - dd_(detM, b) * Pl[a, c])) for (a, b) in ((0, 1), (0, 2), (1, 2)) for c in range(3)] Mnum = np.array(M.subs({x: 0, y: 0, zz: 0}).evalf(30).tolist(), dtype=float) return res, at0(detM), np.linalg.eigvals(Mnum).real r, p, q = sp.symbols("r p q") def metric_jet(A2, E, Phis, pt): return sp.diag(*[jet_rational(c, pt, (r, p, q)) for c in (A2, Phis**2 / E**2, Phis**2 / E**2)]) # ---------------------------------------------------------------- P2 print("P2 tilt with no chart drift, convex branch, F = 0") E0 = (p**2 + q**2 + 1) / 2 n0 = (p / E0, q / E0, 1 - 1 / E0) PhiS, sS = r, r / 2 + r**2 / 3 bS = (r / 5, r**2 / 7, sp.Rational(1, 3)) A2t = sp.diff(PhiS, r)**2 + (sS * sp.diff(PhiS, r) + PhiS * sum(bi * ni for bi, ni in zip(bS, n0)))**2 ok2, det2 = True, [] for pt in ((sp.Rational(1, 2), sp.Rational(1, 3), sp.Rational(-1, 5)), (sp.Rational(2, 3), sp.Rational(-1, 2), sp.Rational(1, 4))): res, dM, ev = thomas_residual(metric_jet(A2t, E0, PhiS, pt)) ok2 &= all(v == 0 for v in res) and dM > 0 and min(ev) > 0 det2.append(f"point {pt}: residual 0 {all(v == 0 for v in res)}, M eigenvalues {np.round(ev, 6).tolist()}") # the l = 1 part of the radial coefficient is nonzero at first order in b while W = 0 (P5 gives the formula) check("tilt-only family: exact F = 0, M > 0, chart drift W = 0, radial coefficient with an l = 1 part", ok2, "; ".join(det2)) # ---------------------------------------------------------------- P3 print("P3 the general swept family, and its Szekeres subfamily") PfS, QfS, SfS = r / 4, r**2 / 3, 1 + r / 5 ES = ((p - PfS)**2 + (q - QfS)**2 + SfS**2) / (2 * SfS) nS = ((p - PfS) / ES, (q - QfS) / ES, 1 - SfS / ES) HS = sp.diff(PhiS, r) - PhiS * sp.diff(ES, r) / ES alphaS = r**2 / 3 + r / 7 b2 = (r / 3, -r / 6, r**2 / 5) A2g = HS**2 + (alphaS + PhiS * sum(bi * ni for bi, ni in zip(b2, nS)))**2 pt = (sp.Rational(1, 2), sp.Rational(1, 3), sp.Rational(-1, 5)) resg, dMg, evg = thomas_residual(metric_jet(A2g, ES, PhiS, pt)) # Szekeres subfamily: alpha = s Phi', Phi b.n = s Phi W.n = s (H - Phi'), i.e. alpha + Phi b.n = s H kS = sS**2 / (1 + sS**2) A2z = HS**2 / (1 - kS) ident = sp.simplify(HS**2 + (sS * HS)**2 - A2z) == 0 resz, dMz, evz = thomas_residual(metric_jet(A2z, ES, PhiS, pt)) ok3 = all(v == 0 for v in resg) and dMg > 0 and ident and all(v == 0 for v in resz) and dMz > 0 check("general swept family F = 0 (det M > 0); Szekeres subfamily A^2 = H^2/(1-k) identically, and F = 0", ok3, f"general: M eigenvalues {np.round(evg, 6).tolist()}; Szekeres: M eigenvalues {np.round(evz, 6).tolist()}") # ---------------------------------------------------------------- P4 print("P4 route 1 of the second-order removal is the tilt") rr, th, eps = sp.symbols("r theta epsilon", positive=True) A2r = sp.cos(rr)**2 + sp.sin(rr)**2 * (1 + eps / sp.sin(rr)**2 * sp.cos(th))**2 Aser = sp.series(sp.sqrt(A2r), eps, 0, 3).removeO() P2 = (3 * sp.cos(th)**2 - 1) / 2 want4 = 1 + eps * sp.cos(th) + eps**2 * (sp.cot(rr)**2 / 6 + sp.cot(rr)**2 / 3 * P2) ok4 = sp.simplify(sp.expand(Aser - want4).rewrite(sp.cos)) == 0 # and the base: Phi = sin r, s = tan r, alpha = s Phi' = sin r, so A0^2 = cos^2 r + sin^2 r = 1 (the unit S^3) check("round base: A = 1 + eps cos(theta) + eps^2 [cot^2 r/6 + (cot^2 r/3) P_2] + O(eps^3)", ok4) # ---------------------------------------------------------------- P5 print("P5 first order, general alpha") Ph5, dPh5, al5, ep = sp.symbols("Phi Phi_p alpha varepsilon", positive=True) Wn, bn = sp.symbols("Wn bn", real=True) A2f = (dPh5 + ep * Ph5 * Wn)**2 + (al5 + ep * Ph5 * bn)**2 A0f = sp.sqrt(dPh5**2 + al5**2) first = sp.simplify(sp.diff(sp.sqrt(A2f), ep).subs(ep, 0) / A0f) expected5 = Ph5 * (dPh5 * Wn + al5 * bn) / A0f**2 ok5 = sp.simplify(first - expected5) == 0 ok5 &= sp.simplify(sp.diff(A2f.subs(Wn, 0), ep).subs(ep, 0) - 2 * al5 * Ph5 * bn) == 0 check("delta A / A0 = Phi (Phi' W + alpha b).n / A0^2; with W = 0, delta(A^2) = 2 alpha Phi b.n", ok5, f"{first}") print("P6 the axisymmetric case, off the round base") b6 = (0, 0, r / 20) A26 = sp.diff(PhiS, r)**2 + (sS * sp.diff(PhiS, r) + PhiS * sum(bi * ni for bi, ni in zip(b6, n0)))**2 res6, dM6, ev6 = thomas_residual(metric_jet(A26, E0, PhiS, (sp.Rational(1, 2), sp.Rational(1, 3), sp.Rational(-1, 5)))) check("axisymmetric tilt, W = 0, non-round base: F = 0 exactly, M > 0", all(v == 0 for v in res6) and dM6 > 0 and min(ev6) > 0, f"M eigenvalues {np.round(ev6, 6).tolist()}") # ---------------------------------------------------------------- P7 print("P7 each shell is intrinsically round") u1, u2, u3, Ph = sp.symbols("u1 u2 u3 Phi", positive=True) pp, qq = sp.symbols("p q", real=True) Ef = ((pp - u1)**2 + (qq - u2)**2 + u3**2) / (2 * u3) lam = Ph / Ef # conformal factor: metric lam^2 (dp^2 + dq^2) Kg = sp.simplify(-(sp.diff(sp.log(lam), pp, 2) + sp.diff(sp.log(lam), qq, 2)) / lam**2) check("Gauss curvature of Phi^2 E^{-2}(dp^2 + dq^2) = 1/Phi^2", sp.simplify(Kg - 1 / Ph**2) == 0, f"K = {Kg}") # ---------------------------------------------------------------- P8 print("P8 the drift in shell coordinates") rr8, th8, ph8 = sp.symbols("r theta phi", real=True) cf, Phf, Zf = sp.Function("c")(rr8), sp.Function("Phi")(rr8), sp.Function("Z")(rr8) nh = sp.Matrix([sp.sin(th8) * sp.cos(ph8), sp.sin(th8) * sp.sin(ph8), sp.cos(th8)]) Xe = sp.Matrix([cf * 0 + Phf * nh[0], Phf * nh[1], cf + Phf * nh[2], Zf]) Jm = Xe.jacobian([rr8, th8, ph8]) G8 = sp.simplify(Jm.T * Jm) A2s = sp.diff(Phf, rr8)**2 + sp.diff(Zf, rr8)**2 ok8 = (sp.simplify(G8[0, 0] - (A2s + 2 * sp.diff(cf, rr8) * sp.diff(Phf, rr8) * sp.cos(th8) + sp.diff(cf, rr8)**2)) == 0 and sp.simplify(G8[0, 1] + sp.diff(cf, rr8) * Phf * sp.sin(th8)) == 0 and sp.simplify(G8[0, 2]) == 0 and sp.simplify(G8[1, 1] - Phf**2) == 0 and sp.simplify(G8[2, 2] - Phf**2 * sp.sin(th8)**2) == 0 and sp.simplify(G8[1, 2]) == 0) ep8 = sp.symbols("varepsilon") cp_round = ep8 / sp.cos(rr8) # Phi' c' = eps on Phi = sin r: radial modulation 2 eps cos(theta) beta8 = sp.simplify(-cp_round * sp.sin(rr8) / ep8) ok8 &= sp.simplify(beta8 + sp.tan(rr8)) == 0 check("g_rr = A^2 + 2c'Phi' cos th + c'^2, g_{r th} = -c' Phi sin th, angular block round; round base beta = -tan r", ok8) print() print(f"test_338: {'ALL PASS' if not FAILURES else 'FAILURES: ' + ', '.join(FAILURES)}") sys.exit(1 if FAILURES else 0)