Listing: Listing

from itertools import product
import sympy as s

R = s.Rational
eta = s.diag(1, -1, -1, -1)
pairs = [(i, j) for i in range(4) for j in range(i, 4)]

B = []
for i, j in pairs:
    X = s.zeros(4)
    X[i, j] = 1
    X[j, i] = 1
    B.append(X)

def tr(X):
    return s.trace(X)

def zero_matrix():
    return s.zeros(4)

def lincomb(coefficients, matrices):
    return sum((c * X for c, X in zip(coefficients, matrices)),
               zero_matrix())

def P(X):
    return X - eta * tr(eta * X) / 2

G = s.Matrix(10, 10,
             lambda i, j: tr(eta * B[i] * eta * B[j]))
Gi = G.inv()
dual = [lincomb([Gi[i, j] for j in range(10)], B)
        for i in range(10)]
Z = [-P(X) for X in dual]

assert G == s.diag(1, -2, -2, -2, 1, 2, 2, 1, 2, 1)

a = s.symbols("a0:10")
H = lincomb(a, B)
u = [eta, -eta * H * eta, eta * H * eta * H * eta]
b1 = tr(eta * H) / 2
b2 = tr(eta * H)**2 / 8 - tr(eta * H * eta * H) / 4

# Curvature-atom computation at the residual jet:
# h(0)=H, first derivatives vanish, h_{,11}=2*eta.
# Since Gamma(0)=0, products of Christoffel symbols vanish.
def second_jet(r, m, n, ell):
    return 2 * eta[n, ell] if r == 1 and m == 1 else 0

def dGamma(n, r, rho, m, v):
    return sum(
        u[n-1][rho, ell] * (
            second_jet(r, m, v, ell)
            + second_jet(r, v, m, ell)
            - second_jet(r, ell, m, v)
        ) / 2
        for ell in range(4)
    )

Ric = {}
for n in (1, 2, 3):
    Ric[n] = s.Matrix(4, 4, lambda m, v:
        sum(dGamma(n, rho, rho, m, v)
            - dGamma(n, v, rho, m, rho)
            for rho in range(4)))

rscalar = {}
for n in (1, 2, 3):
    rscalar[n] = s.expand(sum(
        u[i][m, v] * Ric[n-i][m, v]
        for i in range(n)
        for m, v in product(range(4), repeat=2)
    ))

ell_at_jet = s.expand(-2 * (
    rscalar[3] + b1 * rscalar[2] + b2 * rscalar[1]
))
L_atoms = s.hessian(ell_at_jet, a)

N = eta * H
F = s.expand(-4 * (N*N)[1, 1] - tr(N*N) + tr(N)**2 / 2)
L_trace = s.hessian(F, a)

assert s.expand(ell_at_jet - F) == 0
assert L_atoms == L_trace
L = L_atoms

L_expected = s.Matrix([
 [-1,0,0,0,-1,0,0,-1,0,-1],
 [0,12,0,0,0,0,0,0,0,0],
 [0,0,4,0,0,0,0,0,0,0],
 [0,0,0,4,0,0,0,0,0,0],
 [-1,0,0,0,-9,0,0,1,0,1],
 [0,0,0,0,0,-12,0,0,0,0],
 [0,0,0,0,0,0,-12,0,0,0],
 [-1,0,0,0,1,0,0,-1,0,1],
 [0,0,0,0,0,0,0,0,-4,0],
 [-1,0,0,0,1,0,0,1,0,-1]
])
assert L == L_expected

# A general symmetric covariant source jet.
js = s.symbols("J00 J01 J02 J03 J11 J12 J13 J22 J23 J33")
J = lincomb(js, B)

# Source method 1: explicit four-index C2.
def source_indices(X):
    U1 = -eta * X * eta
    U2 = eta * X * eta * X * eta
    t = tr(eta * X)
    c1 = t / 2
    c2 = t*t / 8 - tr(eta * X * eta * X) / 4

    def cross(U, V, m, n, p, q):
        return U[m,n]*V[p,q] - (
            U[m,p]*V[q,n] + U[m,q]*V[p,n]
        ) / 2

    def coefficient(m, n, p, q):
        return (
            c2 * cross(eta, eta, m,n,p,q)
            + c1 * (cross(U1, eta, m,n,p,q)
                    + cross(eta, U1, m,n,p,q))
            + cross(U2, eta, m,n,p,q)
            + cross(eta, U2, m,n,p,q)
            + cross(U1, U1, m,n,p,q)
        )

    return s.Matrix(4, 4, lambda m, n:
        s.expand(2 * sum(
            coefficient(m,n,p,q) * J[p,q]
            for p, q in product(range(4), repeat=2)
        )))

# Source method 2: independent matrix expression.
def source_matrix(X):
    V = eta * X * eta
    W = eta * X * eta * X * eta
    t = tr(eta * X)
    c1 = t / 2
    c2 = t*t / 8 - tr(eta * X * eta * X) / 4
    c, k, j = tr(eta*J), tr(V*J), tr(W*J)
    result = 2 * (
        c2 * (eta*c - eta*J*eta)
        + c1 * (-V*c - eta*k + V*J*eta + eta*J*V)
        + W*c + eta*j + V*k
        - W*J*eta - eta*J*W - V*J*V
    )
    return result.applyfunc(s.expand)

assert (source_indices(H) - source_matrix(H)).applyfunc(
    s.expand) == s.zeros(4)

# Direct double polarization, with dual edge vectors.
Qsingle = [source_indices(X) for X in Z]
K_pair = s.zeros(4)
for i in range(10):
    for j in range(i, 10):
        if L[i, j] == 0:
            continue
        polar = (
            source_indices(Z[i] + Z[j])
            - Qsingle[i] - Qsingle[j]
        ) / 2
        multiplicity = 1 if i == j else 2
        K_pair += multiplicity * L[i, j] * polar
K_pair = K_pair.applyfunc(s.expand)

# Independent contraction by a quadratic differential operator.
Hedge = lincomb(a, Z)
Qedge = source_matrix(Hedge)
origin = {x: 0 for x in a}
K_diff = s.Matrix(4, 4, lambda m, n:
    s.expand(sum(
        L[i,j] * s.diff(Qedge[m,n], a[i], a[j]).subs(origin) / 2
        for i, j in product(range(10), repeat=2)
    )))
assert K_pair == K_diff

j00,j01,j02,j03,j11,j12,j13,j22,j23,j33 = js
K_expected = s.Matrix([
 [14*j11+10*j22+10*j33, -14*j01, -10*j02, -10*j03],
 [-14*j01, 14*(j00-j22-j33), 14*j12, 14*j13],
 [-10*j02, 14*j12, 10*j00-14*j11-10*j33, 10*j23],
 [-10*j03, 14*j13, 10*j23, 10*j00-14*j11-10*j22]
])
assert K_pair == K_expected

print("L:")
for row in L.tolist():
    print(row)
print("K upper triangular:")
for m, n in pairs:
    print((m,n), s.expand(K_pair[m,n]))
print("All exact polynomial checks passed.")