#!/usr/bin/env python3 """ test_325 -- THE REALIZABILITY OPERATOR'S PRINCIPAL SYMBOL AT EVERY NONDEGENERATE POINT: ITS KERNEL, WHERE IT IS ELLIPTIC (THE FIRST NEWTON TENSOR), ITS FOUR IDENTITIES, AND CARTAN'S TEST For MA ยง4 (the paper this script is deposited with). Independent of tests 228-232 and 246: the symbol is built fresh, exactly where it can be. Set-up. At a point, in an orthonormal frame, F(g) = Cod(K~(g)) with K~ = sqrt(det M) M^{-1}, M = -G = adj K. The third-order principal symbol is sigma(xi) = sigma(Cod)(xi) o (D adj_K)^{-1} o sigma_2(DM)(xi), a map from symmetric 2-tensors h (6) to Codazzi tensors C_abc = -C_bac (9). Write Q_xi = P K P, P the projection off xi, and T_1 = (tr K) g - K, the first Newton transformation. PRE-REGISTERED: P1 Kernel. For every K with det K != 0 and every xi != 0, sigma(xi) annihilates the gauge directions xi (.) v and Q_xi, exactly (symbolic: K diagonal with symbolic entries, xi symbolic), and has rank exactly 2 (exact rational arithmetic at random rational K, definite and indefinite, and random xi), so ker sigma(xi) = gauge + span{Q_xi}, dimension 4. P2 Ellipticity. sigma(xi) restricted to TT_xi (transverse to xi, trace-free) is injective iff xi^T T_1 xi != 0. Exact counterexample K = diag(2,3,-3), xi = e1: h = diag(0,1,-1) is TT and sigma(e1) h = 0. A characteristic covector exists iff k1 + k2 <= 0 <= k2 + k3 (k1 <= k2 <= k3), i.e. iff T_1 is not definite; definite K (and indefinite K with T_1 definite, e.g. diag(2,2,-1)) have none. At fixed xi the fourth direction is conformal iff Q_xi is scalar: K = diag(3,2,2), xi = e1, non-umbilic, has P in the kernel. P3 The count. The nine components of C carry exactly four identities at the symbol level: the cyclic identity C_abc + C_bca + C_cab = 0, and the contracted Bianchi identity div(adj K~) = div M = 0, which for any symmetric K reads B_b(C) = K^c_b C^a_ca - sigma_1 C^a_ba + K^ac C_abc = 0, linear in C. These four functionals are independent at generic K and span the left null space of the order-3 symbol map tau_3 : S^3 V* (x) Sym^2 -> U, whose rank is 5 (as test_232). So F is five genuine conditions per point, not nine (the 7 October review's correction of "nine against six"). P4 Cartan's test (paper 4 sec. 8 used rank tau_4 = 3 rank tau_3, which is not Cartan's test). Here: rank tau_4 = 15 (as test_232), but with g_q = ker tau_q and g_{q,k} its elements free of the first k variables in generic coordinates, dim g_4 = 75 < dim g_3 + dim g_{3,1} + dim g_{3,2} = 55 + 20 + 4 = 79. The symbol is NOT involutive at order 3, and paper 4's "no first-order involutivity failure" is withdrawn. (Registered from paper 4's own numbers before running.) FIRST RUN, AND RE-POSED: the registered claim failed, the other way. Built independently, tau_4 has rank 12, not 15, and the full restriction of tau_3 to a 2-plane of covectors has rank 5, not 4 (test_231's "2, 4, 5" stacked two single covectors, which is not the restriction Cartan's test uses). So dim g_4 = 90 - 12 = 78 and the Cartan sum is 55 + 19 + 4 = 78: EQUALITY, the symbol IS involutive at order 3. The three extra relations at the first prolongation are xi ^ C = 0, the symbol of d_nabla Cod(K) = d_nabla^2 K = R.K, which vanishes identically once Gauss holds (R^i_j = K^i ^ K_j, and K^i ^ K_j ^ K^j = 0). Paper 4's conclusion stands, by Cartan's test; its numbers (rank tau_4 = 15, characters (2,2,1)) are corrected. Where test_232's 15 came from: its prolongation multiplied each block by the monomial's multiplicity (m4.count(a)); in jet coordinates d_a d^alpha g = d^(alpha + e_a) g with coefficient one. Run without that factor, test_232's own construction gives rank 12 and dim g_4 = 78 (checked 26 October). Re-posed P4: rank tau_4 = 12; the left null space of tau_4 is spanned by the twelve prolonged algebraic identities and the three functionals xi ^ C; dim g_4 = 78 = dim g_3 + dim g_{3,1} + dim g_{3,2} = 55 + 19 + 4. P5 The generality. dim g_q grows as 2q^2 + O(q) (second differences equal 4 for q = 6..9), i.e. the general formal solution depends on four functions of three variables: three for diffeomorphisms and one for the normal displacement, matching the per-covector kernel of dimension 4. The joint kernel of sigma(xi) over all xi is span{K}: the normal-variation direction delta g = -2 phi K. P6 (Added 26 October after P1-P5's run; registered before its own run. The sign in (ii) was fixed in a flat-space scratch first.) The relations of P3-P4 are exact identities, not only symbol identities. For every symmetric K on a Riemannian 3-manifold, with C = Cod K, C_abc = nabla_a K_bc - nabla_b K_ac: (i) C_abc + C_bca + C_cab = 0; (ii) nabla^a (adj K)_ab = B_b(C), with adj K = K^2 - sigma_1 K + sigma_2 g; (iii) the cyclic sum over (d,a,b) of nabla_d C_abc equals -(R^e_cda K_be + R^e_cab K_de + R^e_cbd K_ae) (Ricci identity plus first Bianchi), and this vanishes identically when R_ecda = K_ed K_ca - K_ea K_cd (Gauss). Checked exactly at the origin of a random curved polynomial metric with a random polynomial K, and (iii)'s cancellation symbolically. So every left null vector of tau_4 is an exact identity, the prolongation to order 4 adds no condition, and with P4 the system F = 0 is involutive. BUNDLE EDIT (29 October, before the bundle's final run; claims unchanged): sympy's nsimplify removed from the value- at-the-origin helper of P6. On a number with no free symbols nsimplify evaluates numerically and can return a different algebraic number agreeing to about 15 digits; here every value is rational, so it was a no-op. The check now compares exact expressions directly. The opening line's reference to a working document was replaced by the section of the paper. The archive's [M] labels were removed from the claims. """ import sys import itertools import numpy as np import sympy as sp 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) SYM = [(0, 0), (1, 1), (2, 2), (0, 1), (0, 2), (1, 2)] CODIDX = [(a, b, c) for (a, b) in ((0, 1), (0, 2), (1, 2)) for c in range(3)] # C_abc, a < b def sym_basis(): out = [] for (i, j) in SYM: E = sp.zeros(3, 3) E[i, j] = 1 E[j, i] = 1 out.append(E) return out def sigma_M(xi, h): """principal symbol of M' = -G'(h) at a flat point (overall sign immaterial for kernels).""" x2 = (xi.T * xi)[0] trh = h.trace() hx = h * xi Ric = (xi * hx.T + hx * xi.T - x2 * h - xi * xi.T * trh) / 2 R = (xi.T * h * xi)[0] - x2 * trh return -(Ric - R * sp.eye(3) / 2) def dadj(K, X): """derivative of the adjugate at K in direction X (adj K = det K * K^{-1}).""" Ki = K.inv() return K.det() * ((Ki * X).trace() * Ki - Ki * X * Ki) def dadj_inv(K, Y): """solve dadj(K, X) = Y for symmetric X.""" xs = sp.symbols("x0:6") X = sp.zeros(3, 3) for v, (i, j) in zip(xs, SYM): X[i, j] = v X[j, i] = v eqs = [(dadj(K, X) - Y)[i, j] for (i, j) in SYM] sol = sp.solve(eqs, xs, dict=True)[0] return X.subs(sol) def symbol_matrix(K, xi): """9 x 6 matrix of sigma(xi), columns = SYM basis of h.""" cols = [] for E in sym_basis(): L = dadj_inv(K, sigma_M(xi, E)) cols.append([xi[a] * L[b, c] - xi[b] * L[a, c] for (a, b, c) in CODIDX]) return sp.Matrix(cols).T def apply_symbol(K, xi, h): L = dadj_inv(K, sigma_M(xi, h)) return sp.Matrix([xi[a] * L[b, c] - xi[b] * L[a, c] for (a, b, c) in CODIDX]) # ---------------------------------------------------------------- P1 print("P1 the kernel at every nondegenerate point") k1, k2, k3 = sp.symbols("k1 k2 k3", nonzero=True) x, y, z = sp.symbols("x y z", real=True) Ks = sp.diag(k1, k2, k3) xs_ = sp.Matrix([x, y, z]) Pp = sp.eye(3) - xs_ * xs_.T / (x**2 + y**2 + z**2) Q = Pp * Ks * Pp ok = all(sp.simplify(c) == 0 for c in apply_symbol(Ks, xs_, Q)) v1, v2, v3 = sp.symbols("v1 v2 v3") vv = sp.Matrix([v1, v2, v3]) ok &= all(sp.simplify(c) == 0 for c in apply_symbol(Ks, xs_, xs_ * vv.T + vv * xs_.T)) rng = np.random.default_rng(325) ranks = [] for trial in range(12): A = sp.Matrix(3, 3, lambda i, j: sp.Rational(int(rng.integers(-9, 10)), int(rng.integers(1, 6)))) Kr = (A + A.T) / 2 + (sp.eye(3) * 6 if trial < 6 else sp.zeros(3, 3)) # half near-definite, half arbitrary if Kr.det() == 0: continue xi = sp.Matrix([sp.Rational(int(rng.integers(-7, 8)), int(rng.integers(1, 5))) for _ in range(3)]) if (xi.T * xi)[0] == 0: continue ranks.append(symbol_matrix(Kr, xi).rank()) ok &= all(r == 2 for r in ranks) and len(ranks) >= 10 check("sigma(xi) kills xi(.)v and Q_xi = P K P exactly (symbolic K, xi); rank 2 at random rational K and xi", ok, f"exact ranks {ranks}") # ---------------------------------------------------------------- P2 print("P2 where the symbol is elliptic: the first Newton tensor") K0 = sp.diag(2, 3, -3) h0 = sp.diag(0, 1, -1) e1 = sp.Matrix([1, 0, 0]) ok = all(c == 0 for c in apply_symbol(K0, e1, h0)) # the trace of Q_xi is xi^T T1 xi / |xi|^2, with T1 = tr(K) g - K: symbolic identity T1 = Ks.trace() * sp.eye(3) - Ks ok &= sp.simplify(Q.trace() - (xs_.T * T1 * xs_)[0] / (x**2 + y**2 + z**2)) == 0 # Q_xi is transverse; so Q_xi is TT iff tr Q_xi = 0 iff xi^T T1 xi = 0; and kernel = gauge + span{Q} means # the only TT element of the kernel (gauge has no transverse element) is Q itself ok &= all(sp.simplify(c) == 0 for c in (Q * xs_)) # existence of characteristic covectors iff k1 + k2 <= 0 <= k2 + k3 (sorted) agree = True for trial in range(400): kk = np.sort(rng.normal(size=3) * 2) indefinite_T1 = (kk[0] + kk[1]) * (kk[1] + kk[2]) <= 0 t = kk.sum() - kk # T1's eigenvalues k_j + k_l has_null = (t.min() <= 0 <= t.max()) agree &= (indefinite_T1 == has_null) ok &= agree # numerical confirmation that injectivity fails exactly there (float symbol, many K) def fsym(Kf, xf): return np.array(symbol_matrix(sp.Matrix(Kf), sp.Matrix(xf)).evalf(), dtype=float) def tt_basis(xf): xf = xf / np.linalg.norm(xf) u = np.linalg.svd(np.eye(3) - np.outer(xf, xf))[0][:, :2] a, b = u[:, 0], u[:, 1] return [np.outer(a, a) - np.outer(b, b), np.outer(a, b) + np.outer(b, a)] def vec6(Hm): return np.array([Hm[i, j] if i == j else Hm[i, j] for (i, j) in SYM]) def tt_smin(Kf, xf): S = fsym(Kf, xf) # columns of S act on the SYM basis where off-diagonal basis elements carry both (i,j),(j,i): coefficient = H_ij T = np.array([S @ vec6(B) for B in tt_basis(xf)]).T return np.linalg.svd(T, compute_uv=False).min() null_vals, def_vals = [], [] for trial in range(8): kk = rng.normal(size=3) * 2 t = kk.sum() - kk if np.all(t > 0.2) or np.all(t < -0.2): def_vals.append(min(tt_smin(np.diag(kk), rng.normal(size=3)) for _ in range(6))) elif t.min() < -0.2 and t.max() > 0.2: a_, b_ = int(np.argmax(t)), int(np.argmin(t)) xf = np.zeros(3); xf[a_] = np.sqrt(-t[b_]); xf[b_] = np.sqrt(t[a_]) null_vals.append(tt_smin(np.diag(kk), xf)) ok &= all(v < 1e-9 for v in null_vals) and all(v > 1e-3 for v in def_vals) ok &= tt_smin(np.diag([2.0, 2.0, -1.0]), np.array([0.3, -0.5, 0.8])) > 1e-3 # conformal coincidence at fixed xi iff Q scalar: K = diag(3,2,2), xi = e1, P = diag(0,1,1) in the kernel ok &= all(c == 0 for c in apply_symbol(sp.diag(3, 2, 2), e1, sp.diag(0, 1, 1))) ok &= any(c != 0 for c in apply_symbol(sp.diag(3, 2, 1), e1, sp.diag(0, 1, 1))) check("TT-injective at xi iff xi^T T1 xi != 0; diag(2,3,-3) counterexample exact; characteristic covectors iff " "k1+k2 <= 0 <= k2+k3; conformal coincidence at fixed xi iff Q scalar", ok, f"TT singular value at null covectors {max(null_vals) if null_vals else None:.1e}; " f"min over definite-T1 cases {min(def_vals) if def_vals else None:.2e}") # ---------------------------------------------------------------- jets: tau_q def symbol_poly(Kr): """sigma(xi) as a 9 x 6 matrix of cubic polynomials in symbolic xi.""" xs3 = sp.Matrix(sp.symbols("X0:3")) return symbol_matrix(Kr, xs3).applyfunc(sp.expand), list(xs3) def monomials(deg, nvar=3): return [m for m in itertools.product(range(deg + 1), repeat=nvar) if sum(m) == deg] def tau(Spoly, X, q): """the order-q prolonged symbol map: rows (monomial of degree q-3, U), columns (monomial of degree q, W).""" cols = {(M, w): i for i, (M, w) in enumerate((M, w) for M in monomials(q) for w in range(6))} rows = [] polys = [[sp.Poly(Spoly[u, w], *X) for w in range(6)] for u in range(9)] for m in monomials(q - 3): for u in range(9): row = np.zeros(len(cols)) for w in range(6): for mon, coef in polys[u][w].terms(): M = tuple(a + b for a, b in zip(mon, m)) row[cols[(M, w)]] += float(coef) rows.append(row) return np.array(rows), cols def rank(Mx, tol=1e-8): s = np.linalg.svd(Mx, compute_uv=False) return int((s > tol * s.max()).sum()) if s.size else 0 # a generic nondegenerate K in generic coordinates (a random rotation of a definite and of an indefinite K) def rotated(kdiag, seed): Rr = np.linalg.qr(np.random.default_rng(seed).normal(size=(3, 3)))[0] Kf = Rr @ np.diag(kdiag) @ Rr.T return sp.Matrix(3, 3, lambda i, j: sp.Rational(Kf[i, j]).limit_denominator(1000)) def Cfull(vec): C = np.zeros((3, 3, 3)) for val, (a, b, c) in zip(vec, CODIDX): C[a, b, c] = val; C[b, a, c] = -val return C def identities(Kr): """the cyclic identity and the three Bianchi functionals B_b, as functionals on U (rows).""" Kf = np.array(Kr, dtype=float); s1 = np.trace(Kf) basisU = np.eye(9) cyc = np.array([Cfull(e)[0, 1, 2] + Cfull(e)[1, 2, 0] + Cfull(e)[2, 0, 1] for e in basisU]) bianchi = [] for b in range(3): bianchi.append(np.array([np.einsum('c,aca->', Kf[:, b], Cfull(e)) - s1 * np.einsum('aa->', Cfull(e)[:, b, :]) + np.einsum('ac,ac->', Kf, Cfull(e)[:, b, :]) for e in basisU])) return np.array([cyc] + bianchi) # ---------------------------------------------------------------- P3 print("P3 four identities among the nine components; five genuine conditions") ok = True detail = [] for kd, seed in (([1.3, 2.1, 0.8], 1), ([2.0, 3.0, -3.0], 2), ([1.0, 1.7, -0.6], 3)): Kr = rotated(kd, seed) Sp_, X = symbol_poly(Kr) T3, cols = tau(Sp_, X, 3) r3 = rank(T3) # left null space of tau_3 (functionals on U killing its image) U_, s_, Vt = np.linalg.svd(T3.T) left_null = Vt[r3:] # 9 - r3 functionals, as rows on U F4 = identities(Kr) ok &= r3 == 5 and rank(F4) == 4 and np.allclose(F4 @ T3, 0, atol=1e-8) # the four span the whole left null space ok &= rank(np.vstack([F4, left_null])) == 4 detail.append(f"rank tau_3 {r3}") check("rank tau_3 = 5; the cyclic identity and the three Bianchi functionals B_b(C) are independent and span the " "left null space", ok, ", ".join(detail)) # ---------------------------------------------------------------- P4, P5 print("P4 Cartan's test at order 3") ok4 = True; det4 = []; dims_all = {} for kd, seed in (([1.3, 2.1, 0.8], 1), ([2.0, 3.0, -3.0], 2)): Kr = rotated(kd, seed) Sp_, X = symbol_poly(Kr) dims = {} for q in range(3, 10): Tq, cols = tau(Sp_, X, q) rq = rank(Tq) dims[q] = 6 * len(monomials(q)) - rq if q == 4: r4 = rq if q == 3: # g_{3,k}: columns whose monomial has no power of the first k variables gk = [] for kfree in (1, 2): idx = [i for (M, w), i in cols.items() if all(M[j] == 0 for j in range(kfree))] gk.append(len(idx) - rank(Tq[:, idx])) cartan_sum = dims[3] + gk[0] + gk[1] ok4 &= (r4 == 12 and dims[4] == 78 and cartan_sum == 78 and dims[3] == 55) # re-posed (docstring) # the left null space of tau_4: the prolonged four identities (12) and the three xi ^ C functionals T4, cols4 = tau(Sp_, X, 4) F4 = identities(Kr) rows4 = [(m, u) for m in monomials(1) for u in range(9)] def Cvec_fun(fun): out = np.zeros(len(rows4)) for r, (m, u) in enumerate(rows4): d = m.index(1); a, b, c = CODIDX[u] out[r] = fun(d, a, b, c) return out eps = np.zeros((3, 3, 3)) for (i, j, k), sgn in (((0, 1, 2), 1), ((1, 2, 0), 1), ((2, 0, 1), 1), ((0, 2, 1), -1), ((2, 1, 0), -1), ((1, 0, 2), -1)): eps[i, j, k] = sgn wedge = [Cvec_fun(lambda d, a, b, c, cc=cc: 2 * eps[d, a, b] * (c == cc)) for cc in range(3)] # prolonged algebraic identities: (direction d) x (cyclic, B_0, B_1, B_2) as functionals on U def Ufun(vec9): return vec9 alg = [] for dd in range(3): for F in F4: v = np.zeros(len(rows4)) for r, (m, u) in enumerate(rows4): if m.index(1) == dd: v[r] = F[u] alg.append(v) L = np.array(alg + wedge) ok4 &= np.allclose(L @ T4, 0, atol=1e-8) and rank(L) == 15 det4.append(f"K={kd}: rank tau_4 {r4}, dim g_4 {dims[4]}, Cartan sum {dims[3]} + {gk[0]} + {gk[1]} = {cartan_sum}") dims_all[tuple(kd)] = dims check("(re-posed) rank tau_4 = 12, its left null space = 12 prolonged identities + 3 functionals xi ^ C; dim g_4 = 78 = " "55 + 19 + 4: involutive at order 3", ok4, "; ".join(det4)) print("P5 four functions of three variables; the joint kernel is K") ok5 = True; det5 = [] for kd, dims in dims_all.items(): d2 = [dims[q + 1] - 2 * dims[q] + dims[q - 1] for q in range(7, 9)] ok5 &= all(v == 4 for v in d2) det5.append(f"K={list(kd)}: dim g_q (q=3..9) {[dims[q] for q in range(3, 10)]}, second differences {d2}") # joint kernel over many covectors Kr = rotated([1.3, 2.1, 0.8], 1) stack = np.vstack([np.array(symbol_matrix(Kr, sp.Matrix([sp.Rational(int(v), 7) for v in rng.integers(-9, 10, 3)])).evalf(), dtype=float) for _ in range(6)]) _, s_, Vt = np.linalg.svd(stack) kvec = Vt[-1] Kf = np.array(Kr, dtype=float) Kv = np.array([Kf[i, j] for (i, j) in SYM]) cosang = abs(kvec @ Kv) / np.linalg.norm(Kv) ok5 &= rank(stack) == 5 and cosang > 1 - 1e-9 check("dim g_q has second differences 4 (four functions of three variables); joint kernel = span{K}", ok5, "; ".join(det5) + f"; joint-kernel alignment with K {cosang:.12f}") print("P6 the relations are exact identities (cyclic, contracted Bianchi, d_nabla Cod K = R.K = 0 under Gauss)") X3 = sp.symbols("y0:3") rr = np.random.default_rng(3250) def rq(deg): mons = [sp.Integer(1)] + list(X3) + [X3[i] * X3[j] for i in range(3) for j in range(i, 3)] \ + [X3[i] * X3[j] * X3[k] for i in range(3) for j in range(i, 3) for k in range(j, 3)] nm = {0: 1, 1: 4, 2: 10, 3: 20}[deg] return sum(sp.Rational(int(rr.integers(-6, 7)), int(rr.integers(2, 9))) * m for m in mons[:nm]) def trunc(e, deg): """drop monomials of total degree > deg""" p = sp.Poly(sp.expand(e), *X3) return sum(c * sp.prod([X3[i] ** m[i] for i in range(3)]) for m, c in p.terms() if sum(m) <= deg) A = sp.zeros(3, 3) for i in range(3): for j in range(i, 3): a_ij = rq(3) A[i, j] = A[j, i] = sp.expand(a_ij - a_ij.subs({v: 0 for v in X3})) # A(0) = 0, so g(0) = identity g = sp.eye(3) + A ginv = sp.Matrix(3, 3, lambda i, j: trunc((sp.eye(3) - A + A * A)[i, j], 2)) dg = [[[sp.diff(g[i, j], X3[k]) for k in range(3)] for j in range(3)] for i in range(3)] Gam = [[[trunc(sum(ginv[e, l] * (dg[l][a][b] + dg[l][b][a] - dg[a][b][l]) for l in range(3)) / 2, 2) for b in range(3)] for a in range(3)] for e in range(3)] # Gam[e][a][b] = Gamma^e_ab Kp = sp.zeros(3, 3) for i in range(3): for j in range(i, 3): Kp[i, j] = Kp[j, i] = rq(2) + (3 if i == j else 0) at0 = lambda e: sp.expand(e.subs({v: 0 for v in X3})) def nabla2(T): """covariant derivative of a (0,2) field: N[a][b][c] = nabla_a T_bc (truncated)""" return [[[trunc(sp.diff(T[b, c], X3[a]) - sum(Gam[e][a][b] * T[e, c] + Gam[e][a][c] * T[b, e] for e in range(3)), 2) for c in range(3)] for b in range(3)] for a in range(3)] nK = nabla2(Kp) C = [[[nK[a][b][c] - nK[b][a][c] for c in range(3)] for b in range(3)] for a in range(3)] cyc = max(abs(at0(C[a][b][c] + C[b][c][a] + C[c][a][b])) for a in range(3) for b in range(3) for c in range(3)) # (ii) Ke = ginv * Kp # endomorphism K^a_b s1e = trunc(Ke.trace(), 2) s2e = trunc(((Ke.trace()) ** 2 - (Ke * Ke).trace()) / 2, 2) Nadj = sp.Matrix(3, 3, lambda a, b: trunc((Kp * ginv * Kp)[a, b] - s1e * Kp[a, b] + s2e * g[a, b], 2)) nN = nabla2(Nadj) divN = [at0(sum(nN[a][a][b] for a in range(3))) for b in range(3)] # g = delta at 0 K0 = Kp.subs({v: 0 for v in X3}); s10 = K0.trace() C0 = [[[at0(C[a][b][c]) for c in range(3)] for b in range(3)] for a in range(3)] Bb = [sum(K0[c, b] * C0[a][c][a] for a in range(3) for c in range(3)) - s10 * sum(C0[a][b][a] for a in range(3)) + sum(K0[a, c] * C0[a][b][c] for a in range(3) for c in range(3)) for b in range(3)] bian = max(abs(divN[b] - Bb[b]) for b in range(3)) # (iii) Ricci identity: R^e_{f d a} = d_d Gam^e_af - d_a Gam^e_df + Gam^e_dk Gam^k_af - Gam^e_ak Gam^k_df Rm = [[[[at0(sp.diff(Gam[e][a][f], X3[d]) - sp.diff(Gam[e][d][f], X3[a]) + sum(Gam[e][d][k] * Gam[k][a][f] - Gam[e][a][k] * Gam[k][d][f] for k in range(3))) for a in range(3)] for d in range(3)] for f in range(3)] for e in range(3)] # Rm[e][f][d][a] def nabla3(Cf): """nabla_d C_abc at 0""" return [[[[at0(sp.diff(Cf[a][b][c], X3[d]) - sum(Gam[e][d][a] * Cf[e][b][c] + Gam[e][d][b] * Cf[a][e][c] + Gam[e][d][c] * Cf[a][b][e] for e in range(3))) for c in range(3)] for b in range(3)] for a in range(3)] for d in range(3)] nC = nabla3(C) ricci = 0 for d in range(3): for a in range(3): for b in range(3): for c in range(3): lhs = nC[d][a][b][c] + nC[a][b][d][c] + nC[b][d][a][c] rhs = -sum(Rm[e][c][d][a] * K0[b, e] + Rm[e][c][a][b] * K0[d, e] + Rm[e][c][b][d] * K0[a, e] for e in range(3)) ricci = max(ricci, abs(lhs - rhs)) # Gauss cancellation, symbolic kk = sp.symbols("q0:6") Ks6 = sp.zeros(3, 3) for v, (i, j) in zip(kk, SYM): Ks6[i, j] = Ks6[j, i] = v RG = lambda e, c, d, a: Ks6[e, d] * Ks6[c, a] - Ks6[e, a] * Ks6[c, d] gauss = all(sp.expand(sum(RG(e, c, d, a) * Ks6[b, e] + RG(e, c, a, b) * Ks6[d, e] + RG(e, c, b, d) * Ks6[a, e] for e in range(3))) == 0 for d in range(3) for a in range(3) for b in range(3) for c in range(3)) nontriv = max(abs(v) for v in divN) > 0 and max(abs(nC[0][1][2][c]) for c in range(3)) > 0 check("cyclic sum of Cod K = 0; div adj K = B(Cod K); cyclic nabla Cod K = -R.K (Ricci + Bianchi), exactly at a curved point; " "R.K = 0 under Gauss, symbolically", cyc == 0 and bian == 0 and ricci == 0 and gauss and nontriv, f"residuals {cyc}, {bian}, {ricci}; Gauss cancellation {gauss}; div adj K at 0 = {divN}") print() print(f"test_325: {'ALL PASS' if not FAILURES else 'FAILURES: ' + ', '.join(FAILURES)}") sys.exit(1 if FAILURES else 0)