06-ilang-source / 07-dressed-fall
07dressed fallverified
Determines the weight of a body with its own mark in the field of another source and whether such a body falls in the same way for all its clock levels and branches, and to which order.
check_weight.py7 KB · raw
"""Check of item 1 of 06-ilang-source/07-dressed-fall, derivation v2.
N_f = 1 field object on rings of L = 8 and L = 13 places, t_f = 1, source at x_s = 0.
Held body: point term g_b = sigma*gamma*M_k at x; source: point term g_s at x_s.
Exact: W(x) - W(x_s) = eps0(x) - eps0(x_s), eps0 the lowest eigenvalue of
C_1 + g_b|x><x| + g_s|x_s><x_s| (7.1), (7.2). Compared with the operator form (7.3), the ring
form (7.5) and the first-order term (4.4) alone; the deviation is split into its parts
odd and even in g_s, whose scaling in gamma and g_s is fitted. Also checks the slope
factor [W'/W'_lead - 1]/g_b -> -((L-2 r0)^2 + 2)/(6 L t_f) at one place, and, for (7.7), (7.8),
the packet acceleration -<[H,[H,X]]> of H = (E''/4)(1 - cos 2P) + W (E'' = 1) for triangular
envelopes symmetric about x0: factor c_phi and geometric factor ((L-2 r0)^2 + 2 + 12 s^2)/(6 L t_f).
"""
from pathlib import Path
import numpy as np
OUT = Path(__file__).parent / "output"
OUT.mkdir(exist_ok=True)
T_F, XS = 1.0, 0
LINES = []
def log(s=""):
LINES.append(s)
def hop(L):
c = np.zeros((L, L))
for x in range(L):
c[x, (x + 1) % L] = c[(x + 1) % L, x] = -T_F
return c
def w_exact(L, gb, gs):
c1 = hop(L)
e = []
for x in range(L):
h = c1.copy()
h[x, x] += gb
h[XS, XS] += gs
e.append(np.linalg.eigvalsh(h)[0])
e = np.array(e)
return e - e[XS]
def green(L):
lap = 2 * np.eye(L) + hop(L) / T_F
return np.linalg.pinv(lap)
def w_operator(L, gb, gs, order2=True):
n, G = 1.0 / L, green(L)
G2 = G @ G
x = np.arange(L)
gr = G[x, XS]
w = -2 * n * gs * gb / T_F * gr
if order2:
w = w + n * gs * gb**2 / T_F**2 * (2 * G[x, x] * gr + gr**2 - 2.0 / L * G2[x, XS])
return w - w[XS]
def u_of(L):
r = (np.arange(L) - XS) % L
return r * (L - r)
def w_ring(L, gb, gs):
n, u = 1.0 / L, u_of(L)
return n * gs * gb / (L * T_F) * ((1 - (L**2 - 2) * gb / (6 * L * T_F)) * u + gb / (3 * L * T_F) * u**2)
def dev(L, gb, gs, f):
return np.max(np.abs(w_exact(L, gb, gs) - f(L, gb, gs)))
def dev_parts(L, gb, gs):
dp = w_exact(L, gb, gs) - w_ring(L, gb, gs)
dm = w_exact(L, gb, -gs) - w_ring(L, gb, -gs)
return np.max(np.abs(dp - dm)) / 2, np.max(np.abs(dp + dm)) / 2
def slope(xs, ys):
return np.polyfit(np.log(xs), np.log(ys), 1)[0]
for L in (8, 13):
om = 4 * np.sin(np.pi / L) ** 2
log(f"=== L = {L}, N_f = 1, t_f = 1, x_s = 0, omega_min = {om:.6f} (gap t_f*omega_min) ===")
G = green(L)
u = u_of(L)
G0 = (L**2 - 1) / (12 * L)
e1 = np.max(np.abs(G[np.arange(L), XS] - (G0 - u / (2 * L))))
G2 = (G @ G)[np.arange(L), XS]
p = G2 - (-u / (12 * L) - u**2 / (24 * L))
e2 = np.max(np.abs(p - p[0]))
e3 = max(np.max(np.abs(w_operator(L, g, 1e-3) - w_ring(L, g, 1e-3))) for g in (0.03, -0.03))
log(f"closed forms (7.4): max|G_r - (G_0 - u/2L)| = {e1:.1e}; (G^2)_r - (-u/12L - u^2/24L) const to {e2:.1e};"
f" max|(7.3) - (7.5)| = {e3:.1e}")
log("grid: max_x |W_exact - W_(7.5)| and relative to max_x |W_exact|, two levels, both branches")
log(" M_k sigma gamma g_s g_b/gap dev rel.dev dev of (4.4) alone")
for M in (1.0, 2.5):
for sg in (+1, -1):
for gam in (0.002, 0.01):
for gs in (0.001, 0.004):
gb = sg * gam * M
we = w_exact(L, gb, gs)
d = dev(L, gb, gs, w_ring)
d1 = dev(L, gb, gs, lambda L_, a, b: w_operator(L_, a, b, False))
log(f" {M:3.1f} {sg:+d} {gam:6.3f} {gs:6.3f} {gb / om:+8.4f} {d:10.3e} "
f"{d / np.max(np.abs(we)):10.3e} {d1:10.3e}")
gs = 0.002 * om
gbs = om * np.array([0.02, 0.04, 0.08, 0.16])
odd = [dev_parts(L, g, gs)[0] for g in gbs]
first = []
for g in gbs:
dp = w_exact(L, g, gs) - w_operator(L, g, gs, False)
dm = w_exact(L, g, -gs) - w_operator(L, g, -gs, False)
first.append(np.max(np.abs(dp - dm)) / 2)
log(f"gamma scaling at g_s = 0.002 gap, g_b/gap = {[float(round(v, 3)) for v in gbs / om]}:")
log(f" g_s-odd deviation of (7.5): {['%.3e' % v for v in odd]}, slope in g_b = {slope(gbs, odd):.3f}"
" (expected 3)")
log(f" g_s-odd deviation of (4.4) alone: {['%.3e' % v for v in first]}, slope in g_b = "
f"{slope(gbs, first):.3f} (expected 2)")
gb = 0.05 * om
gss = om * np.array([0.002, 0.004, 0.008, 0.016])
even = [dev_parts(L, gb, g)[1] for g in gss]
log(f"g_s scaling at g_b = 0.05 gap, g_s/gap = {[float(round(v, 3)) for v in gss / om]}:")
log(f" g_s-even deviation of (7.5): {['%.3e' % v for v in even]}, slope in g_s = {slope(gss, even):.3f}"
" (expected 2)")
log("local slope factor (s = 0 part of the geometric factor of (7.8)): [W'(x0)/W'_lead(x0) - 1]/g_b, odd part in g_s and"
" in g_b (both branches), Richardson in g_b")
gs, g = 2e-3 * om, 0.01 * om
for r0 in range(1, L):
if 2 * r0 == L:
continue
x0 = (XS + r0) % L
def ratio(gb):
wo = (w_exact(L, gb, gs) - w_exact(L, gb, -gs)) / 2
wp = 0.5 * (wo[(x0 + 1) % L] - wo[(x0 - 1) % L])
return wp / ((1.0 / L) * gs * gb * (L - 2 * r0) / (L * T_F))
c1 = (ratio(g) - ratio(-g)) / (2 * g)
c2 = (ratio(g / 2) - ratio(-g / 2)) / g
c = (4 * c2 - c1) / 3 # Richardson step in g_b
pred = -((L - 2 * r0) ** 2 + 2) / (6 * L * T_F)
log(f" r0 = {r0:2d}: numeric {c:+.6f}, predicted {pred:+.6f}, diff {c - pred:+.1e}")
log("packet acceleration (7.7), (7.8): -<[H,[H,X]]>, H = (1/4)(1 - cos 2P) + W_exact (E'' = 1), triangular"
" envelope of half-width h symmetric about x0, arc not containing x_s, cut 3 places beyond the arc")
log(" h c_phi 12 s^2 max|<[K,[K,X]]>| max|lead ratio - 1| max|geom. factor num - pred| over r0")
T2 = np.roll(np.eye(L), 2, axis=0) # T^2 |x> = |x+2>
K = 0.25 * np.eye(L) - 0.125 * (T2 + T2.T)
for h in range(1, (L - 5) // 2 + 1):
phi = np.array([h + 1 - abs(m) for m in range(-h, h + 1)], dtype=float)
phi /= np.linalg.norm(phi)
w = np.array([phi[m + 1] * phi[m - 1] if 0 < m < 2 * h else 0.0 for m in range(2 * h + 1)])
c_phi = w.sum()
s2 = (w * (np.arange(-h, h + 1)) ** 2).sum() / c_phi
kk, lead, geo = 0.0, 0.0, 0.0
for x0 in range(h + 1, L - h):
if 2 * x0 == L:
continue
psi = np.zeros(L)
psi[x0 - h:x0 + h + 1] = phi
X = np.diag(((np.arange(L) - (x0 + h + 3)) % L).astype(float))
def acc(H):
c = H @ X - X @ H
return -psi @ (H @ c - c @ H) @ psi
kk = max(kk, abs(acc(K)))
D0 = L - 2 * x0
def ratio(gb):
wo = (w_exact(L, gb, gs) - w_exact(L, gb, -gs)) / 2
return (acc(K + np.diag(wo)) - acc(K)) / (-c_phi * (1.0 / L) * gs * gb * D0 / (L * T_F))
lead = max(lead, abs((ratio(g / 2) + ratio(-g / 2)) / 2 - 1))
c1 = (ratio(g) - ratio(-g)) / (2 * g)
c2 = (ratio(g / 2) - ratio(-g / 2)) / g
pred = -(D0**2 + 2 + 12 * s2) / (6 * L * T_F)
geo = max(geo, abs((4 * c2 - c1) / 3 - pred))
log(f" {h} {c_phi:.4f} {12 * s2:7.3f} {kk:10.1e} {lead:10.1e} {geo:10.1e}")
log()
(OUT / "check_weight.txt").write_text("\n".join(LINES) + "\n", encoding="utf-8")
print("\n".join(LINES))
output/check_weight.txt7 KB · raw
=== L = 8, N_f = 1, t_f = 1, x_s = 0, omega_min = 0.585786 (gap t_f*omega_min) ===
closed forms (7.4): max|G_r - (G_0 - u/2L)| = 5.0e-16; (G^2)_r - (-u/12L - u^2/24L) const to 1.8e-15; max|(7.3) - (7.5)| = 5.1e-21
grid: max_x |W_exact - W_(7.5)| and relative to max_x |W_exact|, two levels, both branches
M_k sigma gamma g_s g_b/gap dev rel.dev dev of (4.4) alone
1.0 +1 0.002 0.001 +0.0034 3.116e-10 6.244e-04 9.366e-10
1.0 +1 0.002 0.004 +0.0034 4.986e-09 2.502e-03 7.486e-09
1.0 +1 0.010 0.001 +0.0171 1.505e-09 6.060e-04 1.713e-08
1.0 +1 0.010 0.004 +0.0171 2.461e-08 2.483e-03 8.711e-08
1.0 -1 0.002 0.001 -0.0034 3.125e-10 6.245e-04 3.126e-10
1.0 -1 0.002 0.004 -0.0034 5.000e-09 2.503e-03 2.500e-09
1.0 -1 0.010 0.001 -0.0171 1.527e-09 6.073e-04 1.411e-08
1.0 -1 0.010 0.004 -0.0171 2.496e-08 2.487e-03 3.758e-08
2.5 +1 0.002 0.001 +0.0085 7.725e-10 6.203e-04 4.679e-09
2.5 +1 0.002 0.004 +0.0085 1.242e-08 2.498e-03 2.804e-08
2.5 +1 0.010 0.001 +0.0427 3.104e-09 5.048e-04 1.008e-07
2.5 +1 0.010 0.004 +0.0427 5.841e-08 2.379e-03 4.490e-07
2.5 -1 0.002 0.001 -0.0085 7.779e-10 6.208e-04 3.129e-09
2.5 -1 0.002 0.004 -0.0085 1.251e-08 2.500e-03 3.128e-09
2.5 -1 0.010 0.001 -0.0427 3.257e-09 5.134e-04 9.450e-08
2.5 -1 0.010 0.004 -0.0427 6.065e-08 2.394e-03 3.304e-07
gamma scaling at g_s = 0.002 gap, g_b/gap = [0.02, 0.04, 0.08, 0.16]:
g_s-odd deviation of (7.5): ['1.221e-10', '9.700e-10', '7.755e-09', '6.208e-08'], slope in g_b = 2.997 (expected 3)
g_s-odd deviation of (4.4) alone: ['2.504e-08', '9.980e-08', '3.963e-07', '1.562e-06'], slope in g_b = 1.988 (expected 2)
g_s scaling at g_b = 0.05 gap, g_s/gap = [0.002, 0.004, 0.008, 0.016]:
g_s-even deviation of (7.5): ['6.164e-09', '2.466e-08', '9.863e-08', '3.945e-07'], slope in g_s = 2.000 (expected 2)
local slope factor (s = 0 part of the geometric factor of (7.8)): [W'(x0)/W'_lead(x0) - 1]/g_b, odd part in g_s and in g_b (both branches), Richardson in g_b
r0 = 1: numeric -0.791667, predicted -0.791667, diff +1.2e-07
r0 = 2: numeric -0.375000, predicted -0.375000, diff -3.9e-07
r0 = 3: numeric -0.124999, predicted -0.125000, diff +9.3e-07
r0 = 5: numeric -0.124999, predicted -0.125000, diff +7.3e-07
r0 = 6: numeric -0.374999, predicted -0.375000, diff +9.1e-07
r0 = 7: numeric -0.791666, predicted -0.791667, diff +1.9e-07
packet acceleration (7.7), (7.8): -<[H,[H,X]]>, H = (1/4)(1 - cos 2P) + W_exact (E'' = 1), triangular envelope of half-width h symmetric about x0, arc not containing x_s, cut 3 places beyond the arc
h c_phi 12 s^2 max|<[K,[K,X]]>| max|lead ratio - 1| max|geom. factor num - pred| over r0
1 0.1667 0.000 0.0e+00 2.8e-06 9.3e-07
=== L = 13, N_f = 1, t_f = 1, x_s = 0, omega_min = 0.229088 (gap t_f*omega_min) ===
closed forms (7.4): max|G_r - (G_0 - u/2L)| = 1.1e-15; (G^2)_r - (-u/12L - u^2/24L) const to 5.3e-15; max|(7.3) - (7.5)| = 3.4e-21
grid: max_x |W_exact - W_(7.5)| and relative to max_x |W_exact|, two levels, both branches
M_k sigma gamma g_s g_b/gap dev rel.dev dev of (4.4) alone
1.0 +1 0.002 0.001 +0.0087 5.260e-10 1.062e-03 1.584e-09
1.0 +1 0.002 0.004 +0.0087 8.417e-09 4.260e-03 1.265e-08
1.0 +1 0.010 0.001 +0.0437 2.453e-09 9.988e-04 2.890e-08
1.0 +1 0.010 0.004 +0.0437 4.107e-08 4.193e-03 1.468e-07
1.0 -1 0.002 0.001 -0.0087 5.286e-10 1.062e-03 5.292e-10
1.0 -1 0.002 0.004 -0.0087 8.458e-09 4.263e-03 4.227e-09
1.0 -1 0.010 0.001 -0.0437 2.519e-09 1.004e-03 2.393e-08
1.0 -1 0.010 0.004 -0.0437 4.211e-08 4.209e-03 6.368e-08
2.5 +1 0.002 0.001 +0.0218 1.294e-09 1.048e-03 7.905e-09
2.5 +1 0.002 0.004 +0.0218 2.090e-08 4.245e-03 4.734e-08
2.5 +1 0.010 0.001 +0.1091 3.941e-09 6.521e-04 1.692e-07
2.5 +1 0.010 0.004 +0.1091 9.246e-08 3.837e-03 7.536e-07
2.5 -1 0.002 0.001 -0.0218 1.310e-09 1.050e-03 5.301e-09
2.5 -1 0.002 0.004 -0.0218 2.116e-08 4.252e-03 5.294e-09
2.5 -1 0.010 0.001 -0.1091 4.379e-09 6.870e-04 1.609e-07
2.5 -1 0.010 0.004 -0.1091 9.906e-08 3.898e-03 5.621e-07
gamma scaling at g_s = 0.002 gap, g_b/gap = [0.02, 0.04, 0.08, 0.16]:
g_s-odd deviation of (7.5): ['8.885e-12', '7.049e-11', '5.618e-10', '4.474e-09'], slope in g_b = 2.992 (expected 3)
g_s-odd deviation of (4.4) alone: ['2.537e-09', '1.012e-08', '4.025e-08', '1.592e-07'], slope in g_b = 1.991 (expected 2)
g_s scaling at g_b = 0.05 gap, g_s/gap = [0.002, 0.004, 0.008, 0.016]:
g_s-even deviation of (7.5): ['6.270e-10', '2.508e-09', '1.003e-08', '4.013e-08'], slope in g_s = 2.000 (expected 2)
local slope factor (s = 0 part of the geometric factor of (7.8)): [W'(x0)/W'_lead(x0) - 1]/g_b, odd part in g_s and in g_b (both branches), Richardson in g_b
r0 = 1: numeric -1.576923, predicted -1.576923, diff -2.1e-07
r0 = 2: numeric -1.064097, predicted -1.064103, diff +6.0e-06
r0 = 3: numeric -0.653849, predicted -0.653846, diff -2.6e-06
r0 = 4: numeric -0.346152, predicted -0.346154, diff +2.2e-06
r0 = 5: numeric -0.141038, predicted -0.141026, diff -1.2e-05
r0 = 6: numeric -0.038337, predicted -0.038462, diff +1.2e-04
r0 = 7: numeric -0.038469, predicted -0.038462, diff -7.8e-06
r0 = 8: numeric -0.140966, predicted -0.141026, diff +6.0e-05
r0 = 9: numeric -0.346157, predicted -0.346154, diff -2.9e-06
r0 = 10: numeric -0.653856, predicted -0.653846, diff -1.0e-05
r0 = 11: numeric -1.064105, predicted -1.064103, diff -2.3e-06
r0 = 12: numeric -1.576922, predicted -1.576923, diff +6.7e-07
packet acceleration (7.7), (7.8): -<[H,[H,X]]>, H = (1/4)(1 - cos 2P) + W_exact (E'' = 1), triangular envelope of half-width h symmetric about x0, arc not containing x_s, cut 3 places beyond the arc
h c_phi 12 s^2 max|<[K,[K,X]]>| max|lead ratio - 1| max|geom. factor num - pred| over r0
1 0.1667 0.000 0.0e+00 1.4e-06 1.2e-04
2 0.5263 7.200 0.0e+00 1.2e-06 4.1e-05
3 0.7045 15.484 0.0e+00 1.0e-06 1.4e-05
4 0.8000 26.118 0.0e+00 7.6e-07 9.2e-06