06-ilang-source / 04-profile
04profileverified
Determines for a class of hopping fields how the mutual cost of two held sources depends on their distance and on the dimension of the field graph, and what ties field places to body positions.
check_control_decay.py3 KB · raw
"""Check (iii) of 06-ilang-source/04-profile.
Control field (4.14): one field object with places {empty} + an open chain of L places,
own term C1 = Delta * sum_x |x><x| - J * sum_x (|x><x+1| + h.c.), coupled operator Phi_x = |x><empty| + h.c.
Second-order mutual cost from the single-object contract (generic form (4.3)):
M2 = -2 g1 g2 Re <phi|Phi_x1 R Phi_x2|phi>,
compared with -g1 g2 exp(-d/xi) / (J sinh(1/xi)), cosh(1/xi) = Delta/(2J); decay length fitted from ln|M2|.
"""
from pathlib import Path
import numpy as np
from scipy.linalg import eigh, solve
OUT = Path(__file__).parent / "output"
OUT.mkdir(exist_ok=True)
LOG = []
def out(s=""):
print(s)
LOG.append(s)
L, J, x1 = 400, 1.0, 100 # chain of L places (indices 1..L), place 0 = empty
out(f"Control field on an open chain of L = {L} places, J = {J}, x1 = {x1}, g1 = g2 = 1, N_f = 1")
for ratio in [2.02, 2.1, 2.5, 3.0]:
Delta = ratio * J
n = L + 1
C = np.zeros((n, n))
for x in range(1, n):
C[x, x] = Delta
if x + 1 < n:
C[x, x + 1] = C[x + 1, x] = -J
e, v = eigh(C)
phi, eps = v[:, 0], e[0]
P0 = np.outer(phi, phi)
A = C - eps * np.eye(n) + (e[1] - e[0]) * P0 # invertible; equals C - eps on phi-perp
def R_apply(u):
"""Reduced resolvent of C at eps applied to u."""
return solve(A, u - phi * (phi @ u), assume_a="sym")
def Phi_apply(x, u):
"""Coupled operator Phi_x = |x><empty| + |empty><x| applied to u."""
w = np.zeros(n)
w[x] += u[0]
w[0] += u[x]
return w
theta = np.arccosh(Delta / (2 * J))
xi = 1 / theta
ds = np.arange(1, 161)
M2 = np.array([-2 * (phi @ Phi_apply(x1, R_apply(Phi_apply(x1 + d, phi)))) for d in ds])
form = -np.exp(-ds * theta) / (J * np.sinh(theta))
ok = np.abs(form) > 1e-12 * abs(form[0])
rel = np.max(np.abs(M2[ok] / form[ok] - 1))
slope = np.polyfit(ds[ok], np.log(np.abs(M2[ok])), 1)[0]
out(f"Delta/J = {ratio}: phi = |empty> (overlap {abs(phi[0]):.15f}), gap = {e[1]-e[0]:.4f}; "
f"xi = 1/arccosh(Delta/2J) = {xi:.5f}; fitted decay length = {-1/slope:.5f}; "
f"max |M2/(4.14) - 1| over d = 1..{ds[ok][-1]} = {rel:.1e}")
out(" d : " + " ".join(f"{d:>9d}" for d in [1, 5, 10, 20, 40, 80]))
out(" M2 : " + " ".join(f"{M2[d-1]: .2e}" for d in [1, 5, 10, 20, 40, 80]))
out(" d^6 |M2| : " + " ".join(f"{d**6*abs(M2[d-1]):.2e}" for d in [1, 5, 10, 20, 40, 80]))
(OUT / "control_decay.txt").write_text("\n".join(LOG) + "\n", encoding="utf-8")
check_exact_ring.py6 KB · raw
"""Check (ii) of 06-ilang-source/04-profile.
One field object on a ring of L places, two sources g1 = g2 = g at distance d.
Exact mutual cost M = E(g,g) - E(g,0) - E(0,g) + E(0,0) from
(a) exact diagonalization of h = C1 + g1 P_x1 + g2 P_x2,
(b) the closed eigenvalue equations (4.11), (4.12),
against the second-order form (4.5) and the third-order term (4.8).
Plus one case g2 = 2 g1 for the determinant equation (4.11).
Plus the mean over placements (4.17): Mbar = (1/L) sum_y M(0, y), y = 0..L-1 (y = 0: both sources on
place 0), so that E(0, d) - mean_y E(0, y) = M(d) - Mbar; (4.17) states that Mbar = O(g1^2 g2^2).
"""
from pathlib import Path
import numpy as np
from scipy.linalg import eigh
from scipy.optimize import brentq
OUT = Path(__file__).parent / "output"
OUT.mkdir(exist_ok=True)
LOG = []
def out(s=""):
print(s)
LOG.append(s)
t = 1.0
def ring(L):
return -t * (np.eye(L, k=1) + np.eye(L, k=-1) + np.eye(L, k=L - 1) + np.eye(L, k=-(L - 1)))
def E_exact(C1, sources):
h = C1.copy()
for x, g in sources:
h[x, x] += g
return eigh(h, eigvals_only=True, subset_by_index=[0, 0])[0]
def G(q, d, L):
"""(4.11): Green function of the ring at eps = -2 t cos q (q complex allowed via q = i theta)."""
return np.cos(q * (L / 2 - d)) / (2 * t * np.sin(q) * np.sin(q * L / 2))
def root(f, L):
"""Lowest eigenvalue: the largest root theta > 0 of f(i theta) (eps0 = -2t cosh theta < -2t) if any,
else the smallest root q in (0, 2 pi/L) (eps0 = -2t cos q)."""
fh = lambda th: np.real(f(1j * th))
ths = np.linspace(1.0, 1e-7, 20001)
v = [fh(th) for th in ths]
for a, b, fa, fb in zip(ths[:-1], ths[1:], v[:-1], v[1:]):
if np.sign(fa) != np.sign(fb):
th = brentq(fh, b, a, xtol=1e-16, rtol=1e-15)
return -2 * t * np.cosh(th)
qs = np.linspace(1e-9, 2 * np.pi / L * (1 - 1e-9), 4001)
v = [np.real(f(q)) for q in qs]
for a, b, fa, fb in zip(qs[:-1], qs[1:], v[:-1], v[1:]):
if np.sign(fa) != np.sign(fb):
q = brentq(lambda x: np.real(f(x)), a, b, xtol=1e-16, rtol=1e-15)
return -2 * t * np.cos(q)
raise RuntimeError("no root")
def E_closed(sources, L):
if len(sources) == 1:
g = sources[0][1]
f = lambda q: np.real(1 - g * G(q, 0, L)) * np.real(np.sin(q) * np.sin(q * L / 2))
return root(f, L)
(x1, g1), (x2, g2) = sources
d = abs(x2 - x1)
f = lambda q: np.real(((1 - g1 * G(q, 0, L)) * (1 - g2 * G(q, 0, L)) - g1 * g2 * G(q, d, L) ** 2)
* (np.sin(q) * np.sin(q * L / 2)) ** 2)
return root(f, L)
out("Exact mutual cost on a ring (one field object, t_f = 1), gamma = g L / t_f, w = d(L-d)/L^2")
out("rel2 = (M_exact - M2)/M2, pred3 = M3/M2 from (4.8), rel3 = (M_exact - M2 - M3)/M2,")
out("pred3_L = -(gamma/5)(1-10w+20w^2)/(1-6w) (leading in 1/L, (4.9)); |closed - eigh| = max deviation of (4.11) roots")
for L in [64, 256]:
C1 = ring(L)
P0 = np.full((L, L), 1.0 / L) # projector on the lowest eigenvector of C1
Lp = np.linalg.inv(C1 + 2 * t * np.eye(L) + t * P0) - P0 / t # reduced resolvent R = L^+/t_f
Lp2 = Lp @ Lp
out(f"-- L = {L}")
out(" gamma d w M_exact M2 rel2 pred3 pred3_L rel3 |closed-eigh|")
for gam in [0.03, 0.1, 0.3, 1.0, 3.0, -0.3, -1.0]:
g = gam * t / L
E00 = -2 * t
E10 = E_exact(C1, [(0, g)])
c10 = E_closed([(0, g)], L)
for d in [1, L // 8, L // 4, L // 2]:
E11 = E_exact(C1, [(0, g), (d, g)])
c11 = E_closed([(0, g), (d, g)], L)
Mx = E11 - 2 * E10 + E00
w = d * (L - d) / L**2
M2 = g * g / (t * L * L) * (d * (L - d) - (L * L - 1) / 6)
l0, ld, q0, qd = Lp[0, 0], Lp[0, d], Lp2[0, 0], Lp2[0, d] # R and R^2 entries (t_f = 1)
M3 = g * g * (2 * g) / L * (2 * l0 * ld + ld**2 - (q0 + 2 * qd) / L)
pL = -(gam / 5) * (1 - 10 * w + 20 * w * w) / (1 - 6 * w)
dev = max(abs(c11 - E11), abs(c10 - E10))
out(f" {gam:5.2f} {d:5d} {w:6.4f} {Mx: .5e} {M2: .5e} {(Mx-M2)/M2: .3e} {M3/M2: .3e} {pL: .3e} "
f"{(Mx-M2-M3)/M2: .2e} {dev:.1e}")
out("")
out("Unequal sources g2 = 2 g1 (determinant equation (4.11)), L = 64, d = 8:")
L, d = 64, 8
C1 = ring(L)
for gam in [0.3, 1.0, -1.0]:
g1 = gam / L
g2 = 2 * g1
ex = E_exact(C1, [(0, g1), (d, g2)])
cl = E_closed([(0, g1), (d, g2)], L)
Mx = ex - E_exact(C1, [(0, g1)]) - E_exact(C1, [(d, g2)]) + (-2 * t)
M2 = g1 * g2 / (t * L * L) * (d * (L - d) - (L * L - 1) / 6)
out(f" gamma1 = {gam:5.2f}: E_eigh = {ex:.15f}, E_(4.11) = {cl:.15f}, (M_exact - M2)/M2 = {(Mx-M2)/M2: .3e}")
out("")
out("Mean over placements (4.17), g1 = g2 = g, all y = 0..L-1 (y = 0: both sources on one place):")
out("Mbar = (1/L) sum_y M(0,y); inv = E(0,d) - mean_y E(0,y) = M(d) - Mbar (computed from total contracts only);")
out("res = max_d |inv - M2 - M3| / max_d |M2| (M2, M3 from (4.5), (4.8)); (4.17): Mbar has no term of order < g^4")
out(" L gamma Mbar*L^2/t_f Mbar*L^2/(t_f gamma^4) |Mbar|/max_d|M2| res")
for L in [64, 256]:
C1 = ring(L)
P0 = np.full((L, L), 1.0 / L)
Lp = np.linalg.inv(C1 + 2 * t * np.eye(L) + t * P0) - P0 / t
Lp2 = Lp @ Lp
ds = np.arange(L)
for gam in [0.1, 0.3, 1.0, -0.3, -1.0]:
g = gam * t / L
E00, E10 = -2 * t, E_exact(C1, [(0, g)])
Ms = np.array([E_exact(C1, [(0, g), (y, g)]) - 2 * E10 + E00 for y in ds])
Mbar = Ms.mean()
M2 = g * g / (t * L * L) * (ds * (L - ds) - (L * L - 1) / 6)
M3 = g * g * (2 * g) / L * (2 * Lp[0, 0] * Lp[0, ds] + Lp[0, ds] ** 2 - (Lp2[0, 0] + 2 * Lp2[0, ds]) / L)
inv = Ms - Mbar
res = np.max(np.abs(inv - M2 - M3)) / np.max(np.abs(M2))
out(f" {L:3d} {gam:6.2f} {Mbar*L*L/t: .4e} {Mbar*L*L/(t*gam**4): .4e} "
f"{abs(Mbar)/np.max(np.abs(M2)):.3e} {res:.2e}")
(OUT / "exact_ring.txt").write_text("\n".join(LOG) + "\n", encoding="utf-8")
check_second_order.py6 KB · raw
"""Check (i) of 06-ilang-source/04-profile.
Second-order mutual cost M2 computed numerically from the single-object contract
C1 = -t_f * (adjacency of the graph) via (4.3):
M2 = -2 N_f g1 g2 phi(x1) phi(x2) R(x1, x2), R = reduced resolvent of C1 at the eigenvalue of phi,
compared with the ring closed form (4.5) and the leading torus forms (4.6), (4.7).
Also: lowest + highest state on even rings (sum over the two ends, (4.13)).
"""
from pathlib import Path
from math import lgamma, log, pi
import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import cg, LinearOperator
OUT = Path(__file__).parent / "output"
OUT.mkdir(exist_ok=True)
LOG = []
def out(s=""):
print(s)
LOG.append(s)
def torus_contract(L, D, t=1.0):
"""C1 = -t * sum over nearest-neighbour edges of the L^D torus (L >= 3)."""
n = L**D
idx = np.arange(n).reshape((L,) * D)
r = np.concatenate([idx.ravel()] * D)
c = np.concatenate([np.roll(idx, -1, axis=ax).ravel() for ax in range(D)])
A = sp.coo_matrix((np.ones(len(r)), (r, c)), shape=(n, n)).tocsr()
return -t * (A + A.T), idx
def phi_and_R(C1, which):
"""Eigenvector phi (which = 0: lowest, -1: highest) and reduced resolvent R of C1 at its eigenvalue."""
e, v = np.linalg.eigh(C1)
j = which % len(e)
w = np.array([0.0 if k == j else 1.0 / (e[k] - e[j]) for k in range(len(e))])
return v[:, j], (v * w) @ v.T
def m2_dense(C1, x1, x2, which, Nf=1, g1=1.0, g2=1.0, pr=None):
"""(4.3): M2 = -2 N_f g1 g2 phi(x1) phi(x2) R(x1, x2)."""
phi, R = pr if pr is not None else phi_and_R(C1, which)
return -2 * Nf * g1 * g2 * phi[x1] * phi[x2] * R[x1, x2]
def resolvent_column_cg(C1, phi, eps, x2, t=1.0):
"""R(., x2) from (C1 - eps + t |phi><phi|) y = e_x2 - phi phi(x2) by conjugate gradients."""
n = C1.shape[0]
op = LinearOperator((n, n), matvec=lambda u: C1 @ u - eps * u + t * phi * (phi @ u))
b = -phi * phi[x2]
b[x2] += 1.0
y, info = cg(op, b, rtol=1e-13, atol=0.0, maxiter=20000)
assert info == 0
return y
# ---------------- ring: closed form (4.5) ----------------
out("== Ring: numerical (4.3) vs closed form (4.5); M2 * t_f * L^2 / (N_f g1 g2) vs d(L-d) - (L^2-1)/6 ==")
for L in [3, 4, 16, 101, 400]:
for t in [1.0, 2.0]:
C1 = -t * (np.eye(L, k=1) + np.eye(L, k=-1) + np.eye(L, k=L - 1) + np.eye(L, k=-(L - 1)))
if L == 3:
C1 = -t * (np.ones((3, 3)) - np.eye(3))
dev = 0.0
pr = phi_and_R(C1, 0)
for d in range(1, L):
M2 = m2_dense(C1, 0, d, 0, pr=pr)
closed = (d * (L - d) - (L * L - 1) / 6) / (t * L * L)
dev = max(dev, abs(M2 - closed) / ((L * L - 1) / (6 * t * L * L)))
out(f"L={L:4d} t_f={t}: max_d |M2 - (4.5)| / max|(4.5)| = {dev:.2e}")
out("Ring L=3, d=1, t_f=1: M2 = %.12f, 2/27 = %.12f" % (m2_dense(-(np.ones((3, 3)) - np.eye(3)), 0, 1, 0), 2 / 27))
out("")
out("== Even rings: lowest + highest state (sum over the two ends, (4.13)); odd ring: top degeneracy ==")
for L in [4, 16, 64, 15]:
C1 = -(np.eye(L, k=1) + np.eye(L, k=-1) + np.eye(L, k=L - 1) + np.eye(L, k=-(L - 1)))
if L % 2 == 0:
p0, p1 = phi_and_R(C1, 0), phi_and_R(C1, -1)
s = max(abs(m2_dense(C1, 0, d, 0, pr=p0) + m2_dense(C1, 0, d, -1, pr=p1)) for d in range(1, L))
out(f"L={L}: max_d |M2_lowest + M2_highest| = {s:.2e}")
else:
e = np.linalg.eigvalsh(C1)
out(f"L={L} (odd): two highest eigenvalues of C1 = {e[-1]:.12f}, {e[-2]:.12f} (degenerate)")
# ---------------- square and cubic tori ----------------
kappa2 = 2 * lgamma(0.25) - log(2) - 0.5 * log(pi) # ln(Gamma(1/4)^2 / (2 sqrt(pi)))
xi3 = -2.837297 / (4 * pi)
out("")
out(f"kappa_2 = {kappa2:.6f}, xi_3 = {xi3:.6f}")
def torus_table(D, Ls, rs, dirs, lead):
for L in Ls:
C1, idx = torus_contract(L, D)
n = L**D
phi = np.full(n, n**-0.5)
eps = phi @ (C1 @ phi)
res = np.linalg.norm(C1 @ phi - eps * phi)
x2 = 0
y = resolvent_column_cg(C1, phi, eps, x2)
out(f"-- D={D}, L={L}: |X|={n}, eps0={eps:.6f}, ||C1 phi - eps0 phi|| = {res:.1e}, phi>0 (Perron)")
out(" d vector |d| m_num m_lead m_num-m_lead")
maxdev = 0.0
for dv in dirs:
for r in rs:
d = tuple(r * a for a in dv)
if max(d) > L // 4:
continue
x1 = idx[d]
M2 = -2 * phi[x1] * phi[x2] * y[x1] # N_f = g1 = g2 = t_f = 1
m_num = M2 * n # M2 * t_f / (n g1 g2)
dist = float(np.linalg.norm(d))
m_lead = lead(dist, L)
out(f" {str(d):15s} {dist:6.2f} {m_num: .6f} {m_lead: .6f} {m_num - m_lead: .2e}")
if 4 <= dist <= L / 8:
maxdev = max(maxdev, abs(m_num - m_lead))
out(f" max |m_num - m_lead| over 4 <= |d| <= L/8: {maxdev:.2e}")
out("")
out("== Square torus: m = M2 t_f/(n g1 g2) vs leading form (4.6): m_lead = -(1/pi)[ln(L/|d|) - kappa_2] ==")
torus_table(2, [64, 128, 256], [2, 4, 8, 16, 32], [(1, 0), (1, 1), (2, 1)],
lambda r, L: -(log(L / r) - kappa2) / pi)
out("")
out("== Cubic torus: m vs leading form (4.7): m_lead = -2[1/(4 pi |d|) + xi_3/L] ==")
torus_table(3, [32, 48, 64], [2, 3, 4, 6, 8], [(1, 0, 0), (1, 1, 0), (1, 1, 1)],
lambda r, L: -2 * (1 / (4 * pi * r) + xi3 / L))
(OUT / "second_order.txt").write_text("\n".join(LOG) + "\n", encoding="utf-8")
output/control_decay.txt2 KB · raw
Control field on an open chain of L = 400 places, J = 1.0, x1 = 100, g1 = g2 = 1, N_f = 1
Delta/J = 2.02: phi = |empty> (overlap 1.000000000000000), gap = 0.0201; xi = 1/arccosh(Delta/2J) = 7.07695; fitted decay length = 7.07695; max |M2/(4.14) - 1| over d = 1..160 = 5.7e-13
d : 1 5 10 20 40 80
M2 : -6.12e+00 -3.48e+00 -1.72e+00 -4.18e-01 -2.48e-02 -8.69e-05
d^6 |M2| : 6.12e+00 5.44e+04 1.72e+06 2.67e+07 1.01e+08 2.28e+07
Delta/J = 2.1: phi = |empty> (overlap 1.000000000000000), gap = 0.1001; xi = 1/arccosh(Delta/2J) = 3.17536; fitted decay length = 3.17536; max |M2/(4.14) - 1| over d = 1..88 = 2.2e-15
d : 1 5 10 20 40 80
M2 : -2.28e+00 -6.47e-01 -1.34e-01 -5.74e-03 -1.06e-05 -3.57e-11
d^6 |M2| : 2.28e+00 1.01e+04 1.34e+05 3.68e+05 4.33e+04 9.37e+00
Delta/J = 2.5: phi = |empty> (overlap 1.000000000000000), gap = 0.5001; xi = 1/arccosh(Delta/2J) = 1.44270; fitted decay length = 1.44270; max |M2/(4.14) - 1| over d = 1..40 = 2.4e-15
d : 1 5 10 20 40 80
M2 : -6.67e-01 -4.17e-02 -1.30e-03 -1.27e-06 -1.21e-12 -1.10e-24
d^6 |M2| : 6.67e-01 6.51e+02 1.30e+03 8.14e+01 4.97e-03 2.89e-13
Delta/J = 3.0: phi = |empty> (overlap 1.000000000000000), gap = 1.0001; xi = 1/arccosh(Delta/2J) = 1.03904; fitted decay length = 1.03904; max |M2/(4.14) - 1| over d = 1..29 = 2.9e-15
d : 1 5 10 20 40 80
M2 : -3.42e-01 -7.27e-03 -5.91e-05 -3.91e-09 -1.71e-17 -3.26e-34
d^6 |M2| : 3.42e-01 1.14e+02 5.91e+01 2.50e-01 7.00e-08 8.55e-23
output/exact_ring.txt8 KB · raw
Exact mutual cost on a ring (one field object, t_f = 1), gamma = g L / t_f, w = d(L-d)/L^2
rel2 = (M_exact - M2)/M2, pred3 = M3/M2 from (4.8), rel3 = (M_exact - M2 - M3)/M2,
pred3_L = -(gamma/5)(1-10w+20w^2)/(1-6w) (leading in 1/L, (4.9)); |closed - eigh| = max deviation of (4.11) roots
-- L = 64
gamma d w M_exact M2 rel2 pred3 pred3_L rel3 |closed-eigh|
0.03 1 0.0154 -3.30464e-08 -3.32326e-08 -5.602e-03 -5.619e-03 -5.625e-03 1.63e-05 2.2e-16
0.03 8 0.1094 -1.25477e-08 -1.25796e-08 -2.535e-03 -2.530e-03 -2.540e-03 -5.34e-06 2.2e-16
0.03 16 0.1875 4.54887e-09 4.58658e-09 -8.221e-03 -8.249e-03 -8.250e-03 2.78e-05 2.2e-16
0.03 32 0.2500 1.82646e-08 1.83195e-08 -2.994e-03 -2.999e-03 -3.000e-03 5.02e-06 4.4e-16
0.10 1 0.0154 -3.62402e-07 -3.69251e-07 -1.855e-02 -1.873e-02 -1.875e-02 1.81e-04 4.4e-16
0.10 8 0.1094 -1.38586e-07 -1.39773e-07 -8.489e-03 -8.432e-03 -8.466e-03 -5.74e-05 4.4e-16
0.10 16 0.1875 4.95764e-08 5.09620e-08 -2.719e-02 -2.750e-02 -2.750e-02 3.06e-04 4.4e-16
0.10 32 0.2500 2.01527e-07 2.03550e-07 -9.940e-03 -9.995e-03 -1.000e-02 5.53e-05 4.4e-16
0.30 1 0.0154 -3.14194e-06 -3.32326e-06 -5.456e-02 -5.619e-02 -5.625e-02 1.63e-03 2.2e-16
0.30 8 0.1094 -1.22553e-06 -1.25796e-06 -2.578e-02 -2.530e-02 -2.540e-02 -4.80e-04 2.2e-16
0.30 16 0.1875 4.22078e-07 4.58658e-07 -7.975e-02 -8.249e-02 -8.250e-02 2.73e-03 2.2e-16
0.30 32 0.2500 1.77793e-06 1.83195e-06 -2.949e-02 -2.999e-02 -3.000e-02 4.96e-04 4.4e-16
1.00 1 0.0154 -3.06651e-05 -3.69251e-05 -1.695e-01 -1.873e-01 -1.875e-01 1.78e-02 4.4e-16
1.00 8 0.1094 -1.27422e-05 -1.39773e-05 -8.836e-02 -8.432e-02 -8.466e-02 -4.04e-03 4.4e-16
1.00 16 0.1875 3.84439e-06 5.09620e-06 -2.456e-01 -2.750e-01 -2.750e-01 2.93e-02 4.4e-16
1.00 32 0.2500 1.84310e-05 2.03550e-05 -9.452e-02 -9.995e-02 -1.000e-01 5.43e-03 4.4e-16
3.00 1 0.0154 -1.94134e-04 -3.32326e-04 -4.158e-01 -5.619e-01 -5.625e-01 1.46e-01 0.0e+00
3.00 8 0.1094 -9.24888e-05 -1.25796e-04 -2.648e-01 -2.530e-01 -2.540e-01 -1.18e-02 2.2e-16
3.00 16 0.1875 1.88729e-05 4.58658e-05 -5.885e-01 -8.249e-01 -8.250e-01 2.36e-01 2.2e-16
3.00 32 0.2500 1.36755e-04 1.83195e-04 -2.535e-01 -2.999e-01 -3.000e-01 4.64e-02 2.2e-16
-0.30 1 0.0154 -3.51540e-06 -3.32326e-06 5.782e-02 5.619e-02 5.625e-02 1.63e-03 4.4e-16
-0.30 8 0.1094 -1.28903e-06 -1.25796e-06 2.470e-02 2.530e-02 2.540e-02 -5.93e-04 0.0e+00
-0.30 16 0.1875 4.97778e-07 4.58658e-07 8.529e-02 8.249e-02 8.250e-02 2.81e-03 0.0e+00
-0.30 32 0.2500 1.88780e-06 1.83195e-06 3.049e-02 2.999e-02 3.000e-02 5.01e-04 4.4e-16
-1.00 1 0.0154 -4.45014e-05 -3.69251e-05 2.052e-01 1.873e-01 1.875e-01 1.79e-02 4.4e-16
-1.00 8 0.1094 -1.50412e-05 -1.39773e-05 7.612e-02 8.432e-02 8.466e-02 -8.20e-03 4.4e-16
-1.00 16 0.1875 6.66095e-06 5.09620e-06 3.070e-01 2.750e-01 2.750e-01 3.21e-02 4.4e-16
-1.00 32 0.2500 2.25039e-05 2.03550e-05 1.056e-01 9.995e-02 1.000e-01 5.62e-03 4.4e-16
-- L = 256
gamma d w M_exact M2 rel2 pred3 pred3_L rel3 |closed-eigh|
0.03 1 0.0039 -2.22219e-09 -2.23535e-09 -5.886e-03 -5.906e-03 -5.906e-03 1.94e-05 4.4e-16
0.03 32 0.1094 -7.84745e-10 -7.86746e-10 -2.543e-03 -2.539e-03 -2.540e-03 -4.26e-06 4.4e-16
0.03 64 0.1875 2.83784e-10 2.86137e-10 -8.225e-03 -8.250e-03 -8.250e-03 2.46e-05 4.4e-16
0.03 128 0.2500 1.14102e-09 1.14444e-09 -2.996e-03 -3.000e-03 -3.000e-03 4.24e-06 4.4e-16
0.10 1 0.0039 -2.43535e-08 -2.48372e-08 -1.948e-02 -1.969e-02 -1.969e-02 2.11e-04 2.2e-16
0.10 32 0.1094 -8.66715e-09 -8.74163e-09 -8.520e-03 -8.464e-03 -8.466e-03 -5.65e-05 2.2e-16
0.10 64 0.1875 3.09285e-09 3.17930e-09 -2.719e-02 -2.750e-02 -2.750e-02 3.07e-04 2.2e-16
0.10 128 0.2500 1.25896e-08 1.27160e-08 -9.944e-03 -1.000e-02 -1.000e-02 5.54e-05 2.2e-16
0.30 1 0.0039 -2.10755e-07 -2.23535e-07 -5.717e-02 -5.906e-02 -5.906e-02 1.89e-03 2.2e-16
0.30 32 0.1094 -7.66398e-08 -7.86746e-08 -2.586e-02 -2.539e-02 -2.540e-02 -4.73e-04 2.2e-16
0.30 64 0.1875 2.63314e-08 2.86137e-08 -7.976e-02 -8.250e-02 -8.250e-02 2.73e-03 2.2e-16
0.30 128 0.2500 1.11068e-07 1.14444e-07 -2.950e-02 -3.000e-02 -3.000e-02 4.97e-04 2.2e-16
1.00 1 0.0039 -2.04534e-06 -2.48372e-06 -1.765e-01 -1.969e-01 -1.969e-01 2.04e-02 2.2e-16
1.00 32 0.1094 -7.96704e-07 -8.74163e-07 -8.861e-02 -8.464e-02 -8.466e-02 -3.97e-03 4.4e-16
1.00 64 0.1875 2.39830e-07 3.17930e-07 -2.457e-01 -2.750e-01 -2.750e-01 2.93e-02 2.2e-16
1.00 128 0.2500 1.15137e-06 1.27160e-06 -9.456e-02 -1.000e-01 -1.000e-01 5.44e-03 2.2e-16
3.00 1 0.0039 -1.28013e-05 -2.23535e-05 -4.273e-01 -5.906e-01 -5.906e-01 1.63e-01 4.4e-16
3.00 32 0.1094 -5.78102e-06 -7.86746e-06 -2.652e-01 -2.539e-01 -2.540e-01 -1.13e-02 2.2e-16
3.00 64 0.1875 1.17749e-06 2.86137e-06 -5.885e-01 -8.250e-01 -8.250e-01 2.37e-01 0.0e+00
3.00 128 0.2500 8.54277e-06 1.14444e-05 -2.535e-01 -3.000e-01 -3.000e-01 4.64e-02 2.2e-16
-0.30 1 0.0039 -2.37166e-07 -2.23535e-07 6.098e-02 5.906e-02 5.906e-02 1.92e-03 0.0e+00
-0.30 32 0.1094 -8.06263e-08 -7.86746e-08 2.481e-02 2.539e-02 2.540e-02 -5.85e-04 4.4e-16
-0.30 64 0.1875 3.10548e-08 2.86137e-08 8.531e-02 8.250e-02 8.250e-02 2.81e-03 0.0e+00
-0.30 128 0.2500 1.17935e-07 1.14444e-07 3.050e-02 3.000e-02 3.000e-02 5.02e-04 4.4e-16
-1.00 1 0.0039 -3.02602e-06 -2.48372e-06 2.183e-01 1.969e-01 1.969e-01 2.15e-02 4.4e-16
-1.00 32 0.1094 -9.41059e-07 -8.74163e-07 7.653e-02 8.464e-02 8.466e-02 -8.11e-03 4.4e-16
-1.00 64 0.1875 4.15579e-07 3.17930e-07 3.071e-01 2.750e-01 2.750e-01 3.21e-02 4.4e-16
-1.00 128 0.2500 1.40593e-06 1.27160e-06 1.056e-01 1.000e-01 1.000e-01 5.64e-03 4.4e-16
Unequal sources g2 = 2 g1 (determinant equation (4.11)), L = 64, d = 8:
gamma1 = 0.30: E_eigh = -1.999791523936965, E_(4.11) = -1.999791523936965, (M_exact - M2)/M2 = -3.896e-02
gamma1 = 1.00: E_eigh = -1.999382359521463, E_(4.11) = -1.999382359521462, (M_exact - M2)/M2 = -1.335e-01
gamma1 = -1.00: E_eigh = -2.000878297773575, E_(4.11) = -2.000878297773575, (M_exact - M2)/M2 = 1.042e-01
Mean over placements (4.17), g1 = g2 = g, all y = 0..L-1 (y = 0: both sources on one place):
Mbar = (1/L) sum_y M(0,y); inv = E(0,d) - mean_y E(0,y) = M(d) - Mbar (computed from total contracts only);
res = max_d |inv - M2 - M3| / max_d |M2| (M2, M3 from (4.5), (4.8)); (4.17): Mbar has no term of order < g^4
L gamma Mbar*L^2/t_f Mbar*L^2/(t_f gamma^4) |Mbar|/max_d|M2| res
64 0.10 1.2774e-08 1.2774e-04 7.666e-06 2.28e-04
64 0.30 9.6820e-07 1.1953e-04 6.456e-05 2.03e-03
64 1.00 9.4761e-05 9.4761e-05 5.687e-04 2.17e-02
64 -0.30 1.1827e-06 1.4601e-04 7.886e-05 2.09e-03
64 -1.00 1.8410e-04 1.8410e-04 1.105e-03 2.37e-02
256 0.10 1.2767e-08 1.2767e-04 7.660e-06 2.29e-04
256 0.30 9.6940e-07 1.1968e-04 6.463e-05 2.04e-03
256 1.00 9.4897e-05 9.4897e-05 5.694e-04 2.18e-02
256 -0.30 1.1840e-06 1.4617e-04 7.893e-05 2.10e-03
256 -1.00 1.8427e-04 1.8427e-04 1.106e-03 2.38e-02
output/second_order.txt7 KB · raw
== Ring: numerical (4.3) vs closed form (4.5); M2 * t_f * L^2 / (N_f g1 g2) vs d(L-d) - (L^2-1)/6 ==
L= 3 t_f=1.0: max_d |M2 - (4.5)| / max|(4.5)| = 5.62e-16
L= 3 t_f=2.0: max_d |M2 - (4.5)| / max|(4.5)| = 5.62e-16
L= 4 t_f=1.0: max_d |M2 - (4.5)| / max|(4.5)| = 1.78e-16
L= 4 t_f=2.0: max_d |M2 - (4.5)| / max|(4.5)| = 1.78e-16
L= 16 t_f=1.0: max_d |M2 - (4.5)| / max|(4.5)| = 1.50e-15
L= 16 t_f=2.0: max_d |M2 - (4.5)| / max|(4.5)| = 1.50e-15
L= 101 t_f=1.0: max_d |M2 - (4.5)| / max|(4.5)| = 6.56e-14
L= 101 t_f=2.0: max_d |M2 - (4.5)| / max|(4.5)| = 6.56e-14
L= 400 t_f=1.0: max_d |M2 - (4.5)| / max|(4.5)| = 1.47e-12
L= 400 t_f=2.0: max_d |M2 - (4.5)| / max|(4.5)| = 1.47e-12
Ring L=3, d=1, t_f=1: M2 = 0.074074074074, 2/27 = 0.074074074074
== Even rings: lowest + highest state (sum over the two ends, (4.13)); odd ring: top degeneracy ==
L=4: max_d |M2_lowest + M2_highest| = 1.39e-17
L=16: max_d |M2_lowest + M2_highest| = 1.30e-15
L=64: max_d |M2_lowest + M2_highest| = 2.78e-15
L=15 (odd): two highest eigenvalues of C1 = 1.956295201468, 1.956295201468 (degenerate)
kappa_2 = 1.310533, xi_3 = -0.225785
== Square torus: m = M2 t_f/(n g1 g2) vs leading form (4.6): m_lead = -(1/pi)[ln(L/|d|) - kappa_2] ==
-- D=2, L=64: |X|=4096, eps0=-4.000000, ||C1 phi - eps0 phi|| = 0.0e+00, phi>0 (Perron)
d vector |d| m_num m_lead m_num-m_lead
(2, 0) 2.00 -0.695059 -0.686022 -9.04e-03
(4, 0) 4.00 -0.469300 -0.465387 -3.91e-03
(8, 0) 8.00 -0.253040 -0.244751 -8.29e-03
(16, 0) 16.00 -0.056436 -0.024116 -3.23e-02
(2, 2) 2.83 -0.573480 -0.575705 2.22e-03
(4, 4) 5.66 -0.358125 -0.355069 -3.06e-03
(8, 8) 11.31 -0.149590 -0.134433 -1.52e-02
(16, 16) 22.63 0.027659 0.086202 -5.85e-02
(4, 2) 4.47 -0.431846 -0.429872 -1.97e-03
(8, 4) 8.94 -0.218862 -0.209237 -9.63e-03
(16, 8) 17.89 -0.027187 0.011399 -3.86e-02
max |m_num - m_lead| over 4 <= |d| <= L/8: 8.29e-03
-- D=2, L=128: |X|=16384, eps0=-4.000000, ||C1 phi - eps0 phi|| = 0.0e+00, phi>0 (Perron)
d vector |d| m_num m_lead m_num-m_lead
(2, 0) 2.00 -0.915339 -0.906658 -8.68e-03
(4, 0) 4.00 -0.688478 -0.686022 -2.46e-03
(8, 0) 8.00 -0.467770 -0.465387 -2.38e-03
(16, 0) 16.00 -0.252726 -0.244751 -7.97e-03
(32, 0) 32.00 -0.056370 -0.024116 -3.23e-02
(2, 2) 2.83 -0.793395 -0.796340 2.95e-03
(4, 4) 5.66 -0.575856 -0.575705 -1.52e-04
(8, 8) 11.31 -0.358750 -0.355069 -3.68e-03
(16, 16) 22.63 -0.149758 -0.134433 -1.53e-02
(32, 32) 45.25 0.027599 0.086202 -5.86e-02
(4, 2) 4.47 -0.650663 -0.650508 -1.55e-04
(8, 4) 8.94 -0.432210 -0.429872 -2.34e-03
(16, 8) 17.89 -0.218948 -0.209237 -9.71e-03
(32, 16) 35.78 -0.027219 0.011399 -3.86e-02
max |m_num - m_lead| over 4 <= |d| <= L/8: 7.97e-03
-- D=2, L=256: |X|=65536, eps0=-4.000000, ||C1 phi - eps0 phi|| = 0.0e+00, phi>0 (Perron)
d vector |d| m_num m_lead m_num-m_lead
(2, 0) 2.00 -1.135885 -1.127294 -8.59e-03
(4, 0) 4.00 -0.908750 -0.906658 -2.09e-03
(8, 0) 8.00 -0.686940 -0.686022 -9.18e-04
(16, 0) 16.00 -0.467447 -0.465387 -2.06e-03
(32, 0) 32.00 -0.252650 -0.244751 -7.90e-03
(2, 2) 2.83 -1.013850 -1.016976 3.13e-03
(4, 4) 5.66 -0.795763 -0.796340 5.77e-04
(8, 8) 11.31 -0.576473 -0.575705 -7.68e-04
(16, 16) 22.63 -0.358907 -0.355069 -3.84e-03
(32, 32) 45.25 -0.149800 -0.134433 -1.54e-02
(4, 2) 4.47 -0.870843 -0.871144 3.00e-04
(8, 4) 8.94 -0.651019 -0.650508 -5.11e-04
(16, 8) 17.89 -0.432288 -0.429872 -2.42e-03
(32, 16) 35.78 -0.218969 -0.209237 -9.73e-03
(64, 32) 71.55 -0.027227 0.011399 -3.86e-02
max |m_num - m_lead| over 4 <= |d| <= L/8: 7.90e-03
== Cubic torus: m vs leading form (4.7): m_lead = -2[1/(4 pi |d|) + xi_3/L] ==
-- D=3, L=32: |X|=32768, eps0=-6.000000, ||C1 phi - eps0 phi|| = 0.0e+00, phi>0 (Perron)
d vector |d| m_num m_lead m_num-m_lead
(2, 0, 0) 2.00 -0.071711 -0.065466 -6.24e-03
(3, 0, 0) 3.00 -0.041074 -0.038940 -2.13e-03
(4, 0, 0) 4.00 -0.026635 -0.025677 -9.58e-04
(6, 0, 0) 6.00 -0.013006 -0.012414 -5.91e-04
(8, 0, 0) 8.00 -0.006579 -0.005783 -7.96e-04
(2, 2, 0) 2.83 -0.042082 -0.042158 7.58e-05
(3, 3, 0) 4.24 -0.023502 -0.023402 -1.01e-04
(4, 4, 0) 5.66 -0.014303 -0.014023 -2.79e-04
(6, 6, 0) 8.49 -0.005343 -0.004645 -6.98e-04
(8, 8, 0) 11.31 -0.001184 0.000044 -1.23e-03
(2, 2, 2) 3.46 -0.031314 -0.031833 5.18e-04
(3, 3, 3) 5.20 -0.016597 -0.016518 -7.90e-05
(4, 4, 4) 6.93 -0.009248 -0.008860 -3.88e-04
(6, 6, 6) 10.39 -0.002170 -0.001203 -9.67e-04
(8, 8, 8) 13.86 0.001009 0.002626 -1.62e-03
max |m_num - m_lead| over 4 <= |d| <= L/8: 9.58e-04
-- D=3, L=48: |X|=110592, eps0=-6.000000, ||C1 phi - eps0 phi|| = 4.3e-14, phi>0 (Perron)
d vector |d| m_num m_lead m_num-m_lead
(2, 0, 0) 2.00 -0.076384 -0.070170 -6.21e-03
(3, 0, 0) 3.00 -0.045711 -0.043644 -2.07e-03
(4, 0, 0) 4.00 -0.031220 -0.030381 -8.38e-04
(6, 0, 0) 6.00 -0.017433 -0.017118 -3.15e-04
(8, 0, 0) 8.00 -0.010770 -0.010487 -2.83e-04
(2, 2, 0) 2.83 -0.046727 -0.046862 1.35e-04
(3, 3, 0) 4.24 -0.028077 -0.028105 2.88e-05
(4, 4, 0) 5.66 -0.018780 -0.018727 -5.23e-05
(6, 6, 0) 8.49 -0.009549 -0.009349 -2.00e-04
(8, 8, 0) 11.31 -0.005031 -0.004660 -3.71e-04
(2, 2, 2) 3.46 -0.035932 -0.036536 6.05e-04
(3, 3, 3) 5.20 -0.021112 -0.021222 1.10e-04
(4, 4, 4) 6.93 -0.013626 -0.013564 -6.21e-05
(6, 6, 6) 10.39 -0.006195 -0.005907 -2.88e-04
(8, 8, 8) 13.86 -0.002602 -0.002078 -5.23e-04
max |m_num - m_lead| over 4 <= |d| <= L/8: 8.38e-04
-- D=3, L=64: |X|=262144, eps0=-6.000000, ||C1 phi - eps0 phi|| = 0.0e+00, phi>0 (Perron)
d vector |d| m_num m_lead m_num-m_lead
(2, 0, 0) 2.00 -0.078728 -0.072522 -6.21e-03
(3, 0, 0) 3.00 -0.048046 -0.045996 -2.05e-03
(4, 0, 0) 4.00 -0.033543 -0.032733 -8.10e-04
(6, 0, 0) 6.00 -0.019720 -0.019470 -2.50e-04
(8, 0, 0) 8.00 -0.013004 -0.012839 -1.65e-04
(2, 2, 0) 2.83 -0.049065 -0.049214 1.49e-04
(3, 3, 0) 4.24 -0.030397 -0.030457 6.04e-05
(4, 4, 0) 5.66 -0.021076 -0.021079 3.44e-06
(6, 6, 0) 8.49 -0.011777 -0.011701 -7.63e-05
(8, 8, 0) 11.31 -0.007166 -0.007012 -1.55e-04
(2, 2, 2) 3.46 -0.038262 -0.038888 6.26e-04
(3, 3, 3) 5.20 -0.023417 -0.023574 1.56e-04
(4, 4, 4) 6.93 -0.015897 -0.015916 1.97e-05
(6, 6, 6) 10.39 -0.008369 -0.008259 -1.10e-04
(8, 8, 8) 13.86 -0.004654 -0.004430 -2.23e-04
max |m_num - m_lead| over 4 <= |d| <= L/8: 8.10e-04