#!/usr/bin/env python3 """ test_347 -- TOTAL ABSOLUTE CURVATURE: SIGNED AGAINST UNSIGNED, THE REVOLUTION REDUCTION, AND J ON THE BERGER FAMILY For MB ยง5.1 (the companion of MA). Conventions as in MA: K_ab = -N . d_a d_b x (the unit sphere with outward normal has K = g); M = -G = adj K; J(g) = int sqrt(det M) dV_g, defined where det M >= 0. PRE-REGISTERED: P1 (the revolution reduction, both ways) For g = dsigma^2 + rho(sigma)^2 (dtheta^2 + sin^2 theta dphi^2), rho an arbitrary function, the intrinsic chain (Christoffels, Ricci, M) gives M = diag((1 - rho'^2)/rho^2, -rho''/rho, -rho''/rho), symbolically. With rho' = cos psi, det M rho^4 = sin^4 psi psi'^2, so sqrt(det M) dV integrates to 4 pi int sin^2 psi |psi'| dsigma. Extrinsically, the hypersurface (rho(sigma) omega, z(sigma)) with rho' = cos psi, z' = sin psi has principal curvatures psi' (meridian) and sin psi / rho (twice), with respect to the normal (sin psi omega, -cos psi), so det K rho^2 = psi' sin^2 psi: the signed integral is 4 pi int sin^2 psi psi'. (Symbolic, with psi an arbitrary function.) P2 (corrugation) For psi = sigma + A sin(2k sigma) on [0, pi]: (a) the signed integral is 2 pi^2 for every A and every integer k >= 1 (exact: the antiderivative of 4 pi sin^2 psi dpsi between psi = 0 and psi = pi), and numerically to 1e-25 at the samples below; (b) the unsigned integral equals 2 pi^2 (to 1e-25) when 2kA <= 1, and exceeds it by more than 1e-3 when 2kA > 1, at (k, A) = (2, 1/5), (2, 1/4), (2, 3/10), (2, 2/5), (3, 3/20), (3, 1/5), (3, 3/10); (c) each profile closes, rho(pi) = 0 (to 1e-25), with rho(sigma)/sin(sigma) > 0.8 on a grid of (0, pi), so each is a closed immersed hypersurface of revolution diffeomorphic to S^3; (d) for 2kA > 1, det M vanishes at the zeros of psi' (to 1e-25 at the first one), so these realizations are not in U; and J computed from the intrinsic M equals the unsigned integral to 1e-20 at (2, 3/10). P3 (J on the Berger family; scale invariance) With Milnor constants (1/r, 1/r, r) (the metric diag(1, 1, r^2) in the frame of the round (1,1,1), whose volume is 16 pi^2): J = 2 pi^2 r^3 sqrt(4 - 3r^2) on 0 < r <= 2/sqrt3; its only critical point there is r = 1, a maximum, so J < 2 pi^2 for every non-round Berger metric in U. Under g -> c^2 g, J is unchanged (det M scales by c^-6, dV by c^3; symbolic, general Milnor frame). """ import sys import sympy as sp import mpmath as mp 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) R3 = range(3) def coord_pipeline(g, X): gi = g.inv() Gam = [[[sp.simplify(sum(gi[c, d] * (sp.diff(g[d, b], X[a]) + sp.diff(g[d, a], X[b]) - sp.diff(g[a, b], X[d])) for d in R3) / 2) for b in R3] for a in R3] for c in R3] def Riem(c, d, a, b): return (sp.diff(Gam[c][b][d], X[a]) - sp.diff(Gam[c][a][d], X[b]) + sum(Gam[c][a][e] * Gam[e][b][d] - Gam[c][b][e] * Gam[e][a][d] for e in R3)) Ric = sp.Matrix(3, 3, lambda d, b: sp.simplify(sum(Riem(c, d, c, b) for c in R3))) Rs = sp.simplify(sum(gi[i, j] * Ric[i, j] for i in R3 for j in R3)) M = (Rs / 2 * sp.eye(3) - gi * Ric).applyfunc(sp.simplify) return M def nodes_for(k, A): """[0, pi] split at the zeros of psi' = 1 + 2kA cos(2k sigma), where |psi'| has kinks.""" pts = [mp.mpf(0), mp.pi] + [m * mp.pi / (2 * k) for m in range(1, 2 * k)] if 2 * k * A > 1: a0 = mp.acos(-1 / (2 * k * A)) for m in range(k + 1): for sgn in (1, -1): v = (sgn * a0 + 2 * mp.pi * m) / (2 * k) if 0 < v < mp.pi: pts.append(v) return sorted(pts) # ------------------------------------------------------------------------------------------------------------- P1 print("P1 the revolution reduction, intrinsic and extrinsic") s, th, ph = sp.symbols("sigma theta phi", real=True) rho = sp.Function("rho")(s) g = sp.diag(1, rho**2, rho**2 * sp.sin(th)**2) M = coord_pipeline(g, [s, th, ph]) r1, r2 = sp.diff(rho, s), sp.diff(rho, s, 2) ok1 = (M - sp.diag((1 - r1**2) / rho**2, -r2 / rho, -r2 / rho)).applyfunc(sp.simplify) == sp.zeros(3, 3) psi = sp.Function("psi")(s) sub = {r2: -sp.sin(psi) * sp.diff(psi, s), r1: sp.cos(psi)} detM = sp.simplify((M.det() * rho**4).subs(sub)) ok1 &= sp.simplify(detM - sp.sin(psi)**4 * sp.diff(psi, s)**2) == 0 check("M = diag((1-rho'^2)/rho^2, -rho''/rho, -rho''/rho); det M rho^4 = sin^4 psi psi'^2 when rho' = cos psi", ok1) # extrinsic: X = (rho omega, z), omega = (cos th, sin th cos ph, sin th sin ph), rho' = cos psi, z' = sin psi rh = sp.Function("rho")(s) zz = sp.Function("z")(s) om = sp.Matrix([sp.cos(th), sp.sin(th) * sp.cos(ph), sp.sin(th) * sp.sin(ph)]) Xv = sp.Matrix([*(rh * om), zz]) Nv = sp.Matrix([*(sp.sin(psi) * om), -sp.cos(psi)]) Y = [s, th, ph] ext = {sp.diff(rh, s): sp.cos(psi), sp.diff(zz, s): sp.sin(psi), sp.diff(rh, s, 2): -sp.sin(psi) * sp.diff(psi, s), sp.diff(zz, s, 2): sp.cos(psi) * sp.diff(psi, s)} tang = [sp.diff(Xv, v).subs(ext) for v in Y] okn = all(sp.simplify(Nv.dot(T)) == 0 for T in tang) and sp.simplify(Nv.dot(Nv) - 1) == 0 gX = sp.Matrix(3, 3, lambda a, b: sp.simplify(tang[a].dot(tang[b]))) KX = sp.Matrix(3, 3, lambda a, b: sp.simplify(-Nv.dot(sp.diff(Xv, Y[a], Y[b]).subs(ext)))) shape = (gX.inv() * KX).applyfunc(sp.simplify) ok1e = okn and gX == sp.diag(1, rh**2, rh**2 * sp.sin(th)**2) ok1e &= (shape - sp.diag(sp.diff(psi, s), sp.sin(psi) / rh, sp.sin(psi) / rh)).applyfunc(sp.simplify) == sp.zeros(3, 3) ok1e &= sp.simplify(shape.det() * rh**2 - sp.diff(psi, s) * sp.sin(psi)**2) == 0 check("extrinsic: principal curvatures psi', sin psi/rho (twice); det K rho^2 = psi' sin^2 psi", ok1e) # ------------------------------------------------------------------------------------------------------------- P2 print("P2 corrugation") k_, A_, sg = sp.symbols("k A s", real=True) anti = 4 * sp.pi * (sg / 2 - sp.sin(2 * sg) / 4) # antiderivative of 4 pi sin^2 in psi ok2a = sp.simplify(anti.subs(sg, sp.pi) - anti.subs(sg, 0) - 2 * sp.pi**2) == 0 psi_end = (sg + A_ * sp.sin(2 * k_ * sg)).subs(sg, sp.pi) ok2a &= sp.simplify(psi_end.subs(k_, 2) - sp.pi) == 0 and sp.simplify(psi_end.subs(k_, 3) - sp.pi) == 0 mp.mp.dps = 40 samples = [(2, mp.mpf(1) / 5), (2, mp.mpf(1) / 4), (2, mp.mpf(3) / 10), (2, mp.mpf(2) / 5), (3, mp.mpf(3) / 20), (3, mp.mpf(1) / 5), (3, mp.mpf(3) / 10)] ok2b, ok2c = True, True rows = [] for k, A in samples: psi_f = lambda u, k=k, A=A: u + A * mp.sin(2 * k * u) dpsi = lambda u, k=k, A=A: 1 + 2 * k * A * mp.cos(2 * k * u) nodes = nodes_for(k, A) sgn = 4 * mp.pi * mp.quad(lambda u: mp.sin(psi_f(u))**2 * dpsi(u), nodes) uns = 4 * mp.pi * mp.quad(lambda u: mp.sin(psi_f(u))**2 * abs(dpsi(u)), nodes) ok2a &= abs(sgn - 2 * mp.pi**2) < mp.mpf("1e-25") if 2 * k * A <= 1: ok2b &= abs(uns - 2 * mp.pi**2) < mp.mpf("1e-25") else: ok2b &= uns - 2 * mp.pi**2 > mp.mpf("1e-3") rho_f = lambda u, psi_f=psi_f, k=k: mp.quad(lambda v: mp.cos(psi_f(v)), mp.linspace(0, u, 2 * k + 1)) ok2c &= abs(rho_f(mp.pi)) < mp.mpf("1e-25") ok2c &= min(rho_f(u) / mp.sin(u) for u in mp.linspace(mp.mpf("0.01"), mp.pi - mp.mpf("0.01"), 60)) > mp.mpf("0.8") rows.append(f"(k={k}, A={mp.nstr(A, 3)}, 2kA={mp.nstr(2 * k * A, 3)}): unsigned/2pi^2 = {mp.nstr(uns / (2 * mp.pi**2), 10)}") check("(a) signed integral = 2 pi^2 for every amplitude (exact antiderivative; 1e-25 at the samples)", ok2a) check("(b) unsigned = 2 pi^2 iff 2kA <= 1, exceeding it beyond", ok2b, "; ".join(rows)) check("(c) each profile closes, rho(pi) = 0, and rho/sin > 0.8 inside", ok2c) k, A = 2, mp.mpf(3) / 10 psi_f = lambda u: u + A * mp.sin(2 * k * u) dpsi = lambda u: 1 + 2 * k * A * mp.cos(2 * k * u) d2psi = lambda u: -4 * k**2 * A * mp.sin(2 * k * u) s0 = mp.acos(-1 / (2 * k * A)) / (2 * k) # first zero of psi' rho_f = lambda u: mp.quad(lambda v: mp.cos(psi_f(v)), [0, u]) def detM_num(u): # from P1's intrinsic M, rho' = cos psi r = rho_f(u) return (1 - mp.cos(psi_f(u))**2) / r**2 * (mp.sin(psi_f(u)) * dpsi(u) / r)**2 Jint = mp.quad(lambda u: 4 * mp.pi * mp.sqrt(detM_num(u)) * rho_f(u)**2, nodes_for(k, A)) uns = 4 * mp.pi * mp.quad(lambda u: mp.sin(psi_f(u))**2 * abs(dpsi(u)), nodes_for(k, A)) check("(d) det M = 0 at the first zero of psi' (outside U); J from the intrinsic M = the unsigned integral", abs(detM_num(s0)) < mp.mpf("1e-25") and abs(Jint - uns) < mp.mpf("1e-20"), f"det M(s0) = {mp.nstr(detM_num(s0), 3)}, J/2pi^2 = {mp.nstr(Jint / (2 * mp.pi**2), 12)}") # ------------------------------------------------------------------------------------------------------------- P3 print("P3 J on the Berger family; scale invariance") r, c = sp.symbols("r c", positive=True) detMb = r**4 * (4 - 3 * r**2) / 64 # test_346 P4 Jb = sp.sqrt(detMb) * 16 * sp.pi**2 * r ok3 = sp.simplify(Jb - 2 * sp.pi**2 * r**3 * sp.sqrt(4 - 3 * r**2)) == 0 crit = sp.solve(sp.diff(r**6 * (4 - 3 * r**2), r), r) ok3 &= [v for v in crit if v.is_positive and v < 2 / sp.sqrt(3)] == [1] ok3 &= sp.diff(r**6 * (4 - 3 * r**2), r, 2).subs(r, 1) < 0 and Jb.subs(r, 1) == 2 * sp.pi**2 ok3 &= all(Jb.subs(r, v) < 2 * sp.pi**2 for v in (sp.Rational(1, 2), sp.Rational(9, 10), sp.Rational(11, 10))) # Berger volume factor: diag(1,1,r^2) in the (1,1,1) frame; the (1,1,1) metric is the sphere of radius 2 ok3 &= sp.sqrt(sp.diag(1, 1, r**2).det()) == r l1, l2, l3 = sp.symbols("l1 l2 l3", real=True) lam = [l1, l2, l3] def milnor_M(lams): mu = [sum(lams) / 2 - lv for lv in lams] return [sp.expand(mu[i] * lams[i] - mu[(i + 1) % 3] * mu[(i + 2) % 3]) for i in R3] dm = sp.prod(milnor_M(lam)) dmc = sp.prod(milnor_M([lv / c for lv in lam])) ok3s = sp.simplify(dmc - dm / c**6) == 0 check("J(Berger r) = 2 pi^2 r^3 sqrt(4 - 3r^2), maximal (= 2 pi^2) only at r = 1", ok3, f"critical points {crit}") check("J is scale invariant: det M(c^2 g) = c^-6 det M(g), dV(c^2 g) = c^3 dV", ok3s) print() if FAILURES: print(f"test_347: {len(FAILURES)} FAILED") sys.exit(1) print("test_347: ALL PASS")