08-ilang-propagation / 09-anchor-formation
09anchor formationverified
Determines how the made field of a held source forms in unpinned oscillator cells with an anchor, uniform part included, and whether mark, formation and weighing hold together.
check_formation.py17 KB · raw
"""Check (ii) of 09-anchor-formation: formation of the place-dependent part (item 2).
Exact solution of the linear equations of the means (M -> infinity),
d^2 q/dl^2 = -K_a q - e_{x_s}, q(0) = dq/dl(0) = 0, K_a = kappa_h*Lap + kappa_a*P_h,
in the window below the wrap time (emulating L -> infinity first):
cubic torus L = 64, 96 by fourth-order Runge-Kutta (dt = 0.05; one run repeated with dt = 0.1);
ring L = 1024, 2048 by the eigen-decomposition of K_a.
Cubic torus: Delta_3^perp(rho, l) = da^perp(x,l) / ref - 1 at x = x_s + rho e_1 in order (O1):
da^perp = da - U with the place average U -> 0 at fixed l as L -> infinity (U = O(l^2/L^3) is
printed; the field da(x,l) itself is L-independent in the window, see the size check),
ref = g_inf(x-h) - g_inf(x-x_s) (L -> infinity limit of the place-dependent part of (5.11)),
compared with (9.11): limit -1 (w<1), 1/(b-1) (w>1), b = |x-h|/rho; finite-rho value of the
plateau -g_inf(x-h)/ref. Anchor far (antipode) and near (h = x_s - 8 e_1).
Wakes: max over 1.3 < w <= w_max of |Delta_3^perp - pred(finite rho)| = |u_0 + I|/|ref| and the w of the maximum;
for the far anchor (free field in the window) the free-field wake max|u_0|/g_inf(rho e_1), i.e. the far-from-anchor
limit of |Delta_3^perp| (Step 7), over w > 1.3 and over the common window 1.3 < w <= 2.5; for the near anchor the
image part max|I|/|ref| with I = da(near) - da(far), and |u_0|/|ref| at the w of the wake maximum.
Also printed: the end of the window w_max = c0*l_max/rho (all formation times are finite-window
values), and the lambda-mean of Delta_3^perp in order (O1), (beta-1) g_inf(x-h)/ref with
beta = g_inf(x_s-h)/c_phi (the static field on Z^3 with the anchor), used by the negative part of (9.11).
Above the threshold kappa_h/gamma_3 = 3.957 (localized mode): kappa_a = 8, near anchor, L = 64; the
localized mode (nu_b from 1 + kappa_a gamma_b(0) = 0 on Z^3 and from a sparse eigen-solver, its share
psi_b(x)psi_b(x_s)/nu_b), and the image term I = da(near) - da(far) in the window (the far-anchor field is
the free field there) for kappa_a = 1 and 8, scaled by d_x d_s^(1/6) as in (9.18).
Ring: Delta_1(rho, l) = grad da / grad da_st - 1 between source and anchor, gradient toward the anchor,
reference 1/kappa_h (L -> infinity limit), compared with (9.13): -1, -1/2 for 1<w<2a-1, 0 after;
also max|Delta_1| for w > 1.4 and the range of Delta_1 across the image front (2.6 < w < 3.4).
Uniform form of Step 8: sup over the window of kappa_h |grad da - F(rho) - F(m)|, m = 2 d_s - rho - 1, with
F(n,l) = (2 kappa_h)^-1 int_0^{2 c0 l} J_{2n+1}(u) du (trapezoid rule; checked against the free-ring
gradient), for kappa_a = 0.25, 1, 4, a = 2,
c0 l <= 8 rho; the localized eigenvalue against 2 kappa_h + sqrt(4 kappa_h^2 + kappa_a^2).
"""
from pathlib import Path
import numpy as np
from scipy.integrate import quad
from scipy.special import ive, jv
from scipy.optimize import brentq
import scipy.sparse as sp
from scipy.sparse.linalg import eigsh
OUT = Path(__file__).parent / "output"
OUT.mkdir(exist_ok=True)
KH, KA = 1.0, 1.0
C0 = np.sqrt(KH)
def g_inf(x):
"""Lattice Green function of kappa_h*Lap on Z^3: int_0^inf prod_i e^{-2t} I_{x_i}(2t) dt."""
n = [abs(int(v)) for v in x]
f = lambda u: np.exp(u) * ive(n[0], 2 * np.exp(u)) * ive(n[1], 2 * np.exp(u)) * ive(n[2], 2 * np.exp(u))
umax = 19.0 # ive is unreliable beyond argument ~1e9
pts = [np.log(max(sum(v * v for v in n), 1) / 6.0)]
val = quad(f, -20, umax, points=pts, limit=800, epsabs=1e-14, epsrel=1e-11)[0]
val += 2.0 / (4 * np.pi) ** 1.5 / np.sqrt(np.exp(umax)) # tail (4 pi t)^{-3/2}
return val / KH
def lap3(q):
return 6 * q - sum(np.roll(q, s, axis=a) for a in range(3) for s in (1, -1))
def run_cubic(L, h, rhos, lmax, dt, ka=KA):
q = np.zeros((L, L, L))
p = np.zeros_like(q)
src = np.zeros_like(q)
src[0, 0, 0] = 1.0
hh = tuple(v % L for v in h)
def acc(qq):
a = -KH * lap3(qq) - src
a[hh] -= ka * qq[hh]
return a
nsteps = int(round(lmax / dt))
lam = np.arange(nsteps + 1) * dt
rec = np.zeros((len(rhos), nsteps + 1))
U = np.zeros(nsteps + 1)
for m in range(1, nsteps + 1):
k1q, k1p = p, acc(q)
k2q, k2p = p + 0.5 * dt * k1p, acc(q + 0.5 * dt * k1q)
k3q, k3p = p + 0.5 * dt * k2p, acc(q + 0.5 * dt * k2q)
k4q, k4p = p + dt * k3p, acc(q + dt * k3q)
q = q + dt / 6 * (k1q + 2 * k2q + 2 * k3q + k4q)
p = p + dt / 6 * (k1p + 2 * k2p + 2 * k3p + k4p)
rec[:, m] = [q[r % L, 0, 0] for r in rhos]
U[m] = q.mean()
return lam, rec, U
def form_lambda(lam, Delta, eps):
bad = np.nonzero(np.abs(Delta) > eps)[0]
if len(bad) == 0:
return lam[0]
if bad[-1] == len(lam) - 1:
return np.inf
return lam[bad[-1] + 1]
def cubic_part(lines, csv):
rhos = [3, 4, 6, 8, 12, 16]
lmax = 40.0
lines.append("CUBIC TORUS (window below the wrap time; reference = L->infinity limit)")
lines.append(f" g_inf(0) = {g_inf((0, 0, 0)):.7f} (gamma_3/kappa_h = 0.2527310)")
gr = {r: g_inf((r, 0, 0)) for r in rhos}
cphi = 1.0 / KA + g_inf((0, 0, 0))
store = {}
for L in (64, 96):
for name, h in (("far", (L // 2, L // 2, L // 2)), ("near", (-8, 0, 0))):
lam, rec, U = run_cubic(L, h, rhos, lmax, 0.05)
store[(L, name)] = (lam, rec, U)
beta = g_inf(h) / cphi
lines.append(f" L={L} anchor {name}: h-x_s={h}, |h-x_s|={np.linalg.norm(h):.2f}; "
f"finite-L place average U(lambda_max)={U[-1]:.3e} "
f"(-lambda_max^2/2L^3 = {-lmax**2/(2*L**3):.3e}; -> 0 in order (O1))")
for i, r in enumerate(rhos):
xh = np.array((r, 0, 0)) - np.array(h)
gxh = g_inf(xh)
ref = gxh - gr[r]
b = np.linalg.norm(xh) / r
D = rec[i] / ref - 1 # order (O1): U -> 0
w = C0 * lam / r
sel = (w > 1.6) & (w <= 3.0)
plat = D[sel].mean()
pred_fin = -gxh / ref
mean_O1 = (beta - 1) * gxh / ref
sw = w > 1.3
dev = np.abs(D[sw] - pred_fin) # = |u_0 + I| / |ref|
jw = np.argmax(dev)
wake, w_wake = dev[jw], w[sw][jw]
if name == "far": # in the window the far-anchor field is the free field: u_0 = da + g_inf(rho e_1)
rel = np.abs(rec[i] + gr[r]) / gr[r]
jf = np.argmax(rel[sw])
extra = (f"free-field wake max|u_0|/g_inf(rho e1): w>1.3 {rel[sw][jf]:.3f} at w={w[sw][jf]:.2f}, "
f"1.3<w<=2.5 {rel[sw & (w <= 2.5)].max():.3f}")
else: # image term I = da(near) - da(far), the far-anchor field being the free field in the window
I = rec[i] - store[(L, "far")][1][i]
w_img = (np.linalg.norm(h) + np.linalg.norm(xh)) / r
extra = (f"image part max_(w>1.3)|I|/|ref| {np.max(np.abs(I[sw])) / abs(ref):.3f} "
f"(image front at w={w_img:.2f}); |u_0|/|ref| at the wake maximum {abs(rec[i][sw][jw] - I[sw][jw] + gr[r]) / abs(ref):.3f}")
lf = [form_lambda(lam, D, e) for e in (0.2, 0.3, 0.5)]
lines.append(f" rho={r:2d} b={b:6.2f}: plateau(1.6<w<=min(3,w_max))={plat:+.3f} "
f"pred(finite rho) {pred_fin:+.3f} pred 1/(b-1) {1/(b-1):+.3f} "
f"wake max_(w>1.3)|D-pred| {wake:.3f} at w={w_wake:.2f} {extra} "
f"min D(w<0.6)={D[w < 0.6].min():+.3f} lambda-mean (O1) {mean_O1:+.3f} "
f"w_max={C0*lmax/r:.1f} c0*l_form/rho within window (eps=.2,.3,.5) = " +
", ".join("not formed" if np.isinf(x) else f"{C0*x/r:.2f}" for x in lf))
for j in range(0, len(lam), 4):
csv.append(f"3,{L},{name},{r},{lam[j]:.3f},{D[j]:.6f}")
# size independence below the wrap time and time-step check
for name in ("far", "near"):
a, b_ = store[(64, name)], store[(96, name)]
d = np.max(np.abs(a[1] - b_[1]))
lines.append(f" size check ({name}): max_l,rho |da_64 - da_96| = {d:.2e} "
f"(max |da| = {np.max(np.abs(b_[1])):.2e})")
lam, rec2, _ = run_cubic(64, (-8, 0, 0), rhos, lmax, 0.1)
d = np.max(np.abs(rec2 - store[(64, "near")][1][:, ::2]))
lines.append(f" time-step check (L=64 near): max |da(dt=0.1) - da(dt=0.05)| = {d:.2e}")
above_threshold(lines, csv, store, rhos, lmax, gr)
def gamma_b0(nu):
"""(nu - kappa_h Lap)^{-1}_{00} on Z^3 for nu > 12 kappa_h: int_0^inf e^{-(nu-12kh)t} ive(0,2kh t)^3 dt."""
f = lambda t: np.exp(-(nu - 12 * KH) * t) * ive(0, 2 * KH * t) ** 3
return quad(f, 0, np.inf, limit=800, epsabs=1e-13, epsrel=1e-11)[0]
def sparse_Ka(L, h, ka):
one = sp.identity(L, format="csr")
shift = sp.csr_matrix((np.ones(L), (np.arange(L), (np.arange(L) + 1) % L)), shape=(L, L))
lap1 = 2 * one - shift - shift.T
lap = (sp.kron(sp.kron(lap1, one), one) + sp.kron(sp.kron(one, lap1), one)
+ sp.kron(sp.kron(one, one), lap1))
K = (KH * lap).tolil()
ih =np.ravel_multi_index(tuple(v % L for v in h), (L, L, L))
K[ih, ih] = K[ih, ih] + ka
return K.tocsc()
def above_threshold(lines, csv, store, rhos, lmax, gr):
L, h, ka8 = 64, (-8, 0, 0), 8.0
lines.append(f" ABOVE THRESHOLD: kappa_a = {ka8} > kappa_h/gamma_3 = {KH/g_inf((0,0,0)):.4f}; L={L}, near anchor h-x_s={h}")
lam, rec8, U8 = run_cubic(L, h, rhos, lmax, 0.05, ka=ka8)
nub_inf = brentq(lambda nu: gamma_b0(nu) - 1 / ka8, 12 * KH + 1e-6, 12 * KH + 10 * ka8)
K = sparse_Ka(L, h, ka8)
nu, psi = eigsh(K, k=1, which="LA")
nub, psib = nu[0], psi[:, 0]
lines.append(f" localized mode: nu_b from 1 + kappa_a gamma_b(0) = 0 on Z^3: {nub_inf:.8f}; sparse eigen-solver (L=64): {nub:.8f}")
cphi8 = 1.0 / ka8 + g_inf((0, 0, 0))
beta8 = g_inf(h) / cphi8
lam1, rec_far, _ = store[(L, "far")]
rec_near1 = store[(L, "near")][1]
i0 = np.ravel_multi_index((0, 0, 0), (L, L, L))
for i, r in enumerate(rhos):
xh = np.array((r, 0, 0)) - np.array(h)
dx, ds = np.linalg.norm(xh), np.linalg.norm(h)
gxh = g_inf(xh)
ref = gxh - gr[r]
b = dx / r
D = rec8[i] / ref - 1
w = C0 * lam / r
plat = D[(w > 1.6) & (w <= 3.0)].mean()
I8 = rec8[i] - rec_far[i]
I1 = rec_near1[i] - rec_far[i]
ix = np.ravel_multi_index((r % L, 0, 0), (L, L, L))
share = psib[ix] * psib[i0] / nub
lf = form_lambda(lam, D, 0.5)
lines.append(f" rho={r:2d} b={b:5.2f}: plateau={plat:+.3f} pred(finite rho) {-gxh/ref:+.3f} 1/(b-1) {1/(b-1):+.3f} "
f"lambda-mean (O1) {(beta8-1)*gxh/ref:+.3f} (beta={beta8:.4f}) min D(w<0.6)={D[w < 0.6].min():+.3f} "
f"c0*l_form/rho (eps=.5, window) = {'not formed' if np.isinf(lf) else f'{C0*lf/r:.2f}'}; "
f"max_window|I|*d_x*d_s^(1/6): kappa_a=1 {np.max(np.abs(I1))*dx*ds**(1/6):.4f}, kappa_a=8 {np.max(np.abs(I8))*dx*ds**(1/6):.4f}; "
f"localized share psi_b(x)psi_b(x_s)/nu_b = {share:+.2e}")
for j in range(0, len(lam), 4):
csv.append(f"3,{L},near_ka8,{r},{lam[j]:.3f},{D[j]:.6f}")
def ring_eig(L, h):
lap = np.zeros((L, L))
i = np.arange(L)
lap[i, i] = 2
lap[i, (i + 1) % L] = -1
lap[i, (i - 1) % L] = -1
K = KH * lap
K[h % L, h % L] += KA
return np.linalg.eigh(K)
def ring_field(w2, V, sites, lam):
c = 1 - np.cos(np.sqrt(w2)[:, None] * lam[None, :])
return np.array([-(V[s % len(w2)] * V[0] / w2) @ c for s in sites])
def ring_part(lines, csv):
lines.append("")
lines.append("RING (exact eigen-decomposition; window below the wrap time; reference grad = 1/kappa_h)")
lines.append(" (a) anchor at d_s = 2 rho (a = 2): prediction -1 (w<1), -1/2 (1<w<3), 0 (w>3)")
for L in (1024, 2048):
for r in (20, 40, 80, 160):
ds = 2 * r
w2, V = ring_eig(L, ds)
lam = np.linspace(0, 4 * r / C0, 1601)
f = ring_field(w2, V, (r, r + 1), lam)
D = (f[1] - f[0]) * KH - 1
w = C0 * lam / r
p1 = D[(w > 1.4) & (w < 2.6)].mean()
p2 = D[w > 3.4].mean()
mx = np.max(np.abs(D[w > 1.4]))
tr = D[(w > 2.6) & (w < 3.4)]
lf = [form_lambda(lam, D, e) for e in (0.25, 0.6)]
lines.append(f" L={L} rho={r:3d}: plateau(1.4<w<2.6)={p1:+.3f} plateau(w>3.4)={p2:+.3f} "
f"max|D|(w>1.4)={mx:.3f} D across image front (2.6<w<3.4) in [{tr.min():+.3f}, {tr.max():+.3f}] "
f"c0*l_form/rho: eps=0.25 -> {C0*lf[0]/r:.2f} (pred 3), eps=0.6 -> {C0*lf[1]/r:.2f} (pred 1);"
f" finite-L ref grad*kappa_h = {(L-ds)/L:.3f}")
for j in range(0, len(lam), 8):
csv.append(f"1,{L},a2,{r},{lam[j]:.3f},{D[j]:.6f}")
lines.append(" (b) anchor far, d_s = 400: prediction -1/2 until c0*l = 2 d_s - rho, then 0")
for L in (1024, 2048):
ds = 400
w2, V = ring_eig(L, ds)
lam = np.linspace(0, 900 / C0, 3601)
for r in (10, 20, 40, 80):
f = ring_field(w2, V, (r, r + 1), lam)
D = (f[1] - f[0]) * KH - 1
sel = (lam > 1.4 * r / C0) & (lam < 0.9 * (2 * ds - r) / C0)
lf = form_lambda(lam, D, 0.25)
lines.append(f" L={L} rho={r:3d}: mean D on 1.4 rho < c0 l < 0.9(2d_s-rho) = {D[sel].mean():+.3f}, "
f"max|D+1/2| there = {np.max(np.abs(D[sel] + 0.5)):.3f}; c0*l_form(eps=0.25) = "
f"{C0*lf:.1f} (pred 2d_s-rho = {2*ds-r}), c0*l_form/rho = {C0*lf/r:.1f}")
for j in range(0, len(lam), 8):
csv.append(f"1,{L},far,{r},{lam[j]:.3f},{D[j]:.6f}")
def F_bessel(n, lam):
"""Free gradient on the line: (2 kappa_h)^-1 int_0^{2 c0 l} J_{2n+1}(u) du (trapezoid rule, step 0.004)."""
X = 2 * C0 * lam
u = np.linspace(0, X.max(), int(X.max() / 0.004) + 2)
J = jv(2 * n + 1, u)
cum = np.concatenate(([0.0], np.cumsum(0.5 * (J[1:] + J[:-1]) * np.diff(u))))
return np.interp(X, u, cum) / (2 * KH)
def F_free_ring(n, lam, L):
"""Free gradient at place n on the ring without anchor (zero mode removed), by the mode sum."""
ks = 2 * np.pi * np.arange(1, L) / L
om = 2 * C0 * np.abs(np.sin(ks / 2))
amp = (np.exp(1j * ks) - 1) * np.exp(1j * ks * n) / om ** 2
out = np.zeros_like(lam)
for c in range(0, len(ks), 256):
out -= (amp[c:c + 256, None] * (1 - np.cos(om[c:c + 256, None] * lam[None, :]))).sum(0).real
return out / L
def ring_uniform(lines):
lines.append(" (c) uniform form of Step 8: sup_window kappa_h|grad da - F(rho) - F(m)|, a = 2, m = 3 rho - 1, c0 l <= 8 rho, L = 2048")
L = 2048
Fcache = {}
bess = []
for ka in (0.25, 1.0, 4.0):
row = []
for r in (20, 40, 80, 160):
ds = 2 * r
i = np.arange(L)
K = np.zeros((L, L))
K[i, i] = 2 * KH
K[i, (i + 1) % L] = -KH
K[i, (i - 1) % L] = -KH
K[ds, ds] += ka
w2, V = np.linalg.eigh(K)
lam = np.linspace(0, 8 * r / C0, 3201)
f = ring_field(w2, V, (r, r + 1), lam)
grad = f[1] - f[0]
m = 2 * ds - r - 1
if r not in Fcache:
Fr, Fm = F_bessel(r, lam), F_bessel(m, lam)
dF = max(np.max(np.abs(Fr - F_free_ring(r, lam, L))), np.max(np.abs(Fm - F_free_ring(m, lam, L)))) * KH
Fcache[r] = (Fr, Fm, dF)
Fr, Fm, dF = Fcache[r]
dev = np.max(np.abs(grad - Fr - Fm)) * KH
row.append(f"rho={r}: {dev:.4f}")
if ka == 1.0:
bess.append(f"rho={r}: {dF:.1e}")
if r == 160:
nub = w2[-1]
lines.append(f" kappa_a={ka}: " + "; ".join(row) + f"; localized eigenvalue {nub:.6f} vs 2kh+sqrt(4kh^2+ka^2) = {2*KH+np.sqrt(4*KH**2+ka**2):.6f}")
lines.append(" Bessel form of F (n = rho and m) against the free-ring gradient, max kappa_h|diff|: " + "; ".join(bess))
def main():
lines = ["Check (ii): formation of the place-dependent part (item 2)",
f"kappa_h = {KH}, kappa_a = {KA}, c_0 = {C0}", ""]
csv = ["D,L,anchor,rho,lambda,Delta"]
cubic_part(lines, csv)
ring_part(lines, csv)
ring_uniform(lines)
(OUT / "formation_check.txt").write_text("\n".join(lines) + "\n", encoding="utf-8")
(OUT / "formation_Delta.csv").write_text("\n".join(csv) + "\n", encoding="utf-8")
print("\n".join(lines))
if __name__ == "__main__":
main()
check_made_field.py6 KB · raw
"""Check (i) of 09-anchor-formation: the made field of item 1.
Exact solution of the linear equations of the means (M -> infinity):
d^2 q/dl^2 = -K_a q - e_{x_s}, q(0) = dq/dl(0) = 0,
computed from the eigen-decomposition of K_a = kappa_h*Lap + kappa_a*P_h.
It is compared with the derived expressions of item 1:
(9.4) image representation: free torus kernels S_0, da_0 plus the image at the
anchor with strength sigma_h from the Volterra equation;
(9.5) uniform part (place average) from the force balance,
U(l) = -L^{-D} int_0^l (l-l') [1 - sigma_h(l')] dl',
and the static split (9.6) of (5.11).
Ring L = 64, 96 and cubic torus L = 8, 12; anchor at two distances; three places.
"""
from pathlib import Path
import numpy as np
OUT = Path(__file__).parent / "output"
OUT.mkdir(exist_ok=True)
KH, KA = 1.0, 1.0
def idx(c, L):
return int(np.ravel_multi_index(tuple(int(v) % L for v in c), (L,) * len(c)))
def laplacian(L, D):
N = L ** D
grid = np.arange(N).reshape((L,) * D)
lap = np.zeros((N, N))
rows = np.arange(N)
for ax in range(D):
for sh in (1, -1):
lap[rows, np.roll(grid, sh, axis=ax).ravel()] -= 1.0
lap[rows, rows] += 2 * D
return lap
def modes(L, D):
ks = 2 * np.pi * np.arange(L) / L
kk = np.meshgrid(*([ks] * D), indexing="ij")
om2 = 2 * KH * sum(1 - np.cos(k) for k in kk)
return [k.ravel() for k in kk], om2.ravel()
def free_kernels(y, t, kk, om2, L, D):
"""S_0(y,t) and da_0(y,t) on the torus (zero mode included)."""
ph = np.cos(sum(k * yi for k, yi in zip(kk, y)))
nz = om2 > 1e-12
ph, w2 = ph[nz], om2[nz]
S = np.zeros_like(t)
A = np.zeros_like(t)
for c in range(0, len(w2), 128):
om = np.sqrt(w2[c:c + 128])[:, None]
arg = om * t[None, :]
S += (ph[c:c + 128, None] * np.sin(arg) / om).sum(0)
A += (ph[c:c + 128, None] * (1 - np.cos(arg)) / om ** 2).sum(0)
S = (S + t) / L ** D
A = -(A + t ** 2 / 2) / L ** D
return S, A
def run(D, L, xs, h, places, lines):
N = L ** D
lap = laplacian(L, D)
K = KH * lap
ih, isrc = idx(h, L), idx(xs, L)
K[ih, ih] += KA
w2, V = np.linalg.eigh(K)
lmax = 3.0 * L
dt = 0.005 if D == 3 else 0.01
t = np.arange(0.0, lmax + dt / 2, dt)
n = len(t)
# --- derived: Volterra equation for sigma_h (trapezoid rule) ---
kk, om2 = modes(L, D)
S00, _ = free_kernels(np.zeros(D), t, kk, om2, L, D)
_, A_hs = free_kernels(np.subtract(h, xs), t, kk, om2, L, D)
f = -KA * A_hs
sig = np.zeros(n)
for m in range(1, n):
conv = 0.5 * S00[m] * sig[0] + np.dot(S00[m - 1:0:-1], sig[1:m])
sig[m] = f[m] - KA * dt * conv # S00[0] = 0
# evaluation times
ev = np.linspace(0, n - 1, 241).astype(int)
lam = t[ev]
errs = []
for p in places:
ip = idx(p, L)
exact = -(V[ip] * V[isrc] / w2) @ (1 - np.cos(np.sqrt(w2)[:, None] * lam[None, :]))
S_ph, _ = free_kernels(np.subtract(p, h), t, kk, om2, L, D)
_, A_ps = free_kernels(np.subtract(p, xs), t, kk, om2, L, D)
img = np.array([dt * (np.dot(S_ph[m:0:-1], sig[:m]) + 0.5 * S_ph[0] * sig[m]
- 0.5 * S_ph[m] * sig[0]) if m > 0 else 0.0 for m in ev])
derived = A_ps[ev] + img
errs.append(np.max(np.abs(derived - exact)) / np.max(np.abs(exact)))
# uniform part: exact place average vs force formula
ones_proj = V.sum(0) / N
U_exact = -(ones_proj * V[isrc] / w2) @ (1 - np.cos(np.sqrt(w2)[:, None] * lam[None, :]))
g1 = 1.0 - sig
U_der = []
for m in ev:
if m == 0:
U_der.append(0.0)
continue
integ = (t[m] - t[:m + 1]) * g1[:m + 1]
U_der.append(-dt * (integ.sum() - 0.5 * integ[0] - 0.5 * integ[-1]) / N)
U_der = np.array(U_der)
errU = np.max(np.abs(U_der - U_exact)) / np.max(np.abs(U_exact))
# static split of (5.11)
Ginv = V @ np.diag(1 / w2) @ V.T
da_st = -Ginv[:, isrc]
# g = (kappa_h Lap)^+ from the plane-wave modes (k = 0 excluded)
nz = om2 > 1e-12
coords = np.array(np.unravel_index(np.arange(N), (L,) * D))
def g_of(y): # y: array (D, M) of displacements
ph = np.cos(sum(kk[a][nz, None] * y[a][None, :] for a in range(D)))
return (ph / om2[nz, None]).sum(0) / N
g_xh = g_of(coords - np.array(h)[:, None])
g_xs = g_of(coords - np.array(xs)[:, None])
g0 = g_of(np.zeros((D, 1)))[0]
g_sh = g_of((np.array(h) - np.array(xs))[:, None])[0]
r_sh = 2 * (g0 - g_sh)
U_st_pred = -(1 / KA + r_sh / 2)
perp_pred = g_xh - g_xs
e_split = max(abs(da_st.mean() - U_st_pred),
np.max(np.abs(da_st - da_st.mean() - perp_pred))) / np.max(np.abs(da_st))
# Bohr mean of sigma_h -> 1 (static image strength on the torus)
sig_static = -KA * da_st[ih]
lines.append(f"D={D} L={L} xs={tuple(xs)} h={tuple(h)} |h-xs|={np.linalg.norm(np.subtract(h, xs)):.3f} "
f"lambda_max={lmax:.0f} dt={dt}")
for p, e in zip(places, errs):
lines.append(f" place {tuple(p)}: max|derived(9.4)-exact|/max|exact| = {e:.2e}")
lines.append(f" uniform part (9.5) force formula vs exact place average: rel. err = {errU:.2e}")
lines.append(f" static split (9.6): rel. err = {e_split:.2e}; static image strength -kappa_a*da_st(h) = {sig_static:.12f}")
return max(errs + [errU, e_split])
def main():
lines = ["Check (i): made field of item 1 against the exact solution of the means",
f"kappa_h = {KH}, kappa_a = {KA}, c_0 = sqrt(kappa_h) = {np.sqrt(KH)}", ""]
worst = 0.0
for L in (64, 96):
for h in ((3,), (17,)):
worst = max(worst, run(1, L, (0,), h, [(-6,), (4,), (11,)], lines))
for L in (8, 12):
for h in ((1, 1, 0), (3, 2, 1)):
worst = max(worst, run(3, L, (0, 0, 0), h, [(2, 0, 0), (1, -2, 1), (-3, 1, 2)], lines))
lines += ["", f"largest relative deviation over all cases: {worst:.2e}"]
(OUT / "made_field_check.txt").write_text("\n".join(lines) + "\n", encoding="utf-8")
print("\n".join(lines))
if __name__ == "__main__":
main()
check_uniform.py12 KB · raw
"""Check (iii) of 09-anchor-formation: the uniform part U(l) = place average of da(x,l) (item 3).
Exact solution of the linear equations of the means (M -> infinity),
d^2 q/dl^2 = -K_a q - e_{x_s}, q(0) = dq/dl(0) = 0, K_a = kappa_h*Lap + kappa_a*P_h.
Cubic torus L = 8, 12, 16, 20 (kappa_a = 1; RK4, dt = 0.05) against (9.14):
U(l) = U_st [1 - cos(l/T_U)] + O(1/L) U_st, T_U = sqrt(L^3 c_phi^(L)), c_phi^(L) = 1/kappa_a + g_L(0),
U_st = -[c_phi^(L) - g_L(x_s-h)];
plus two variations at L = 12 (anchor farther, kappa_a = 4), the lowest normal mode of K_a from
a sparse eigen-solver, and the slow image: place-dependent part of psi_0(x)psi_0(x_s)/nu_0 against
-(1-beta_L) g_L(x-h), beta_L = g_L(x_s-h)/c_phi^(L).
Ring L = 128, 256, 512 (kappa_a = 1, d_s = L/4; exact eigen-decomposition) against (9.15):
U(l) = -(L/2kappa_h) [f(alpha) - f(alpha+tau)/2 - f(alpha-tau)/2], tau = c_0 l/L, alpha = d_s/L,
f the odd 2-periodic extension of x(1-x).
Ring, discrete derivation of (9.15): the even in-band labels from the secular equation (9.19),
tan(kL/2) = kappa_a/(2 kappa_h sin k), against the eigenvalues of K_a; the exact mode sum of Step 9
(weights N_j^2 cos(delta_j) cot(k_j/2)/(4 kappa_h L sin^2(k_j/2)), phases sin(k_j d_s + delta_j)) against
the eigen-decomposition (difference = localized mode, bound 2/(kappa_a L)); the remainder
R = U - U_cont for alpha = 1/8, 1/4, 1/2, kappa_a = 1, 4 over tau <= 8, scaled by kappa_a/(1+tau)
(bound C[(1+tau)/kappa_a + tau^(2/3) L^(-1/3)/kappa_h + log L/(kappa_h L)]), and for a stiff anchor
(kappa_a = 1e4) over tau <= 8 and 32 against tau^(2/3) L^(-1/3).
Ring with a weak anchor (kappa_a L << kappa_h; L = 32, 64): the slow-mode form of (9.14) with
c_phi^(L) = 1/kappa_a + (L^2-1)/(12 kappa_h L), T_U = sqrt(L c_phi^(L)).
"""
from pathlib import Path
import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import eigsh
from scipy.optimize import brentq
OUT = Path(__file__).parent / "output"
OUT.mkdir(exist_ok=True)
KH = 1.0
C0 = np.sqrt(KH)
def g_torus(L, ys):
"""g_L(y) = L^-3 sum_{k != 0} cos(k.y) / (kappa_h Lap)(k), for an array ys of shape (3, M)."""
ks = 2 * np.pi * np.arange(L) / L
k1, k2, k3 = [k.ravel() for k in np.meshgrid(ks, ks, ks, indexing="ij")]
om2 = 2 * KH * (3 - np.cos(k1) - np.cos(k2) - np.cos(k3))
nz = om2 > 1e-12
out = []
for c in range(0, ys.shape[1], 512):
y = ys[:, c:c + 512]
ph = np.cos(k1[nz, None] * y[0] + k2[nz, None] * y[1] + k3[nz, None] * y[2])
out.append((ph / om2[nz, None]).sum(0) / L ** 3)
return np.concatenate(out)
def lap3(q):
return 6 * q - sum(np.roll(q, s, axis=a) for a in range(3) for s in (1, -1))
def run_cubic(L, h, ka, lmax, dt=0.05):
q = np.zeros((L, L, L))
p = np.zeros_like(q)
src = np.zeros_like(q)
src[0, 0, 0] = 1.0
hh = tuple(v % L for v in h)
def acc(qq):
a = -KH * lap3(qq) - src
a[hh] -= ka * qq[hh]
return a
n = int(round(lmax / dt))
U = np.zeros(n + 1)
for m in range(1, n + 1):
k1q, k1p = p, acc(q)
k2q, k2p = p + 0.5 * dt * k1p, acc(q + 0.5 * dt * k1q)
k3q, k3p = p + 0.5 * dt * k2p, acc(q + 0.5 * dt * k2q)
k4q, k4p = p + dt * k3p, acc(q + dt * k3q)
q = q + dt / 6 * (k1q + 2 * k2q + 2 * k3q + k4q)
p = p + dt / 6 * (k1p + 2 * k2p + 2 * k3p + k4p)
U[m] = q.mean()
return np.arange(n + 1) * dt, U
def sparse_Ka(L, h, ka):
one = sp.identity(L, format="csr")
shift = sp.csr_matrix((np.ones(L), (np.arange(L), (np.arange(L) + 1) % L)), shape=(L, L))
lap1 = 2 * one - shift - shift.T
lap = (sp.kron(sp.kron(lap1, one), one) + sp.kron(sp.kron(one, lap1), one)
+ sp.kron(sp.kron(one, one), lap1))
K = (KH * lap).tolil()
ih = np.ravel_multi_index(tuple(v % L for v in h), (L, L, L))
K[ih, ih] += ka
return K.tocsc(), ih
def cubic_case(L, h, ka, lines, csv, tag):
coords = np.array(np.unravel_index(np.arange(L ** 3), (L, L, L)))
g0 = g_torus(L, np.zeros((3, 1)))[0]
gsh = g_torus(L, -np.array(h, dtype=float)[:, None])[0]
cphi = 1 / ka + g0
Ust = -(cphi - gsh)
TU = np.sqrt(L ** 3 * cphi)
lam, U = run_cubic(L, h, ka, 2.2 * np.pi * TU)
pred = Ust * (1 - np.cos(lam / TU))
dev = np.max(np.abs(U - pred)) / abs(Ust)
first = lam < 2 * np.pi * TU
lmin = lam[first][np.argmin(U[first])]
# lowest normal mode (independent, sparse shift-invert)
K, ih = sparse_Ka(L, h, ka)
nu, psi = eigsh(K, k=1, sigma=0, which="LM")
nu0, psi0 = nu[0], psi[:, 0]
A0 = psi0 * psi0[0] / nu0
A0perp = A0 - A0.mean()
beta = gsh / cphi
gxh = g_torus(L, coords - np.array(h)[:, None] % L)
img = -(1 - beta) * gxh
dimg = np.max(np.abs(A0perp - img)) / np.max(np.abs(img))
lines.append(f" L={L:2d} kappa_a={ka} h-x_s={h}: c_phi^(L)={cphi:.5f} U_st={Ust:.5f} T_U={TU:.2f}; "
f"1/sqrt(nu_0)={1/np.sqrt(nu0):.2f} (ratio {1/np.sqrt(nu0)/TU:.5f}); "
f"max|U-pred|/|U_st|={dev:.2e} (x L = {dev*L:.4f}); "
f"first minimum at l/(pi T_U)={lmin/(np.pi*TU):.4f}, U_min/U_st={U[first].min()/Ust:.4f}; "
f"U(2 pi T_U)/U_st={np.interp(2*np.pi*TU, lam, U)/Ust:+.4f}; "
f"slow image: max|A0perp+(1-beta)g(x-h)|/max|.|={dimg:.2e} (x L^2 = {dimg*L*L:.2f}), beta_L={beta:.4f}")
for j in range(0, len(lam), 20):
csv.append(f"3,{tag},{L},{ka},{lam[j]:.3f},{U[j]:.6e},{pred[j]:.6e}")
def fext(x):
y = np.mod(x + 1, 2) - 1 # in [-1, 1)
return np.where(y >= 0, y * (1 - y), y * (1 + y))
def ring_case(L, lines, csv, ka=1.0):
ds = L // 4
i = np.arange(L)
K = np.zeros((L, L))
K[i, i] = 2 * KH
K[i, (i + 1) % L] = -KH
K[i, (i - 1) % L] = -KH
K[ds, ds] += ka
w2, V = np.linalg.eigh(K)
tau = np.linspace(0, 4, 2001)
lam = tau * L / C0
U = -((V.sum(0) / L) * V[0] / w2) @ (1 - np.cos(np.sqrt(w2)[:, None] * lam[None, :]))
alpha = ds / L
Uc = -(L / (2 * KH)) * (fext(alpha) - 0.5 * fext(alpha + tau) - 0.5 * fext(alpha - tau))
Ust = -(1 / ka + ds * (L - ds) / (2 * KH * L))
dev = np.max(np.abs(U - Uc)) / abs(Ust)
dev1 = np.max(np.abs(U - Uc)[tau <= 2]) / abs(Ust)
first = tau < 2
tmin = tau[first][np.argmin(U[first])]
lines.append(f" L={L:3d} d_s={ds}: U_st={Ust:.3f}; max|U-U_cont|/|U_st| over tau<=2: {dev1:.4f} (x L = {dev1*L:.2f}), "
f"over tau<=4: {dev:.4f}; first minimum at c0 l/L = {tmin:.4f} (pred 1), U_min/U_st={U[first].min()/Ust:.4f} (pred 2); "
f"U(c0 l = 2L)/U_st = {np.interp(2.0, tau, U)/Ust:+.4f} (pred 0)")
for j in range(0, len(tau), 4):
csv.append(f"1,ring,{L},{ka},{lam[j]:.3f},{U[j]:.6e},{Uc[j]:.6e}")
def ring_eig(L, ds, ka):
"""Anchor at place ds, source at place 0."""
i = np.arange(L)
K = np.zeros((L, L))
K[i, i] = 2 * KH
K[i, (i + 1) % L] = -KH
K[i, (i - 1) % L] = -KH
K[ds % L, ds % L] += ka
return np.linalg.eigh(K)
def U_eig(w2, V, L, lam):
return -((V.sum(0) / L) * V[0] / w2) @ (1 - np.cos(np.sqrt(w2)[:, None] * lam[None, :]))
def even_labels(L, ka):
"""Roots of -2 kappa_h sin k sin(kL/2) + kappa_a cos(kL/2) = 0 in (2j pi/L, (2j+1) pi/L), L even."""
f = lambda k: -2 * KH * np.sin(k) * np.sin(k * L / 2) + ka * np.cos(k * L / 2)
return np.array([brentq(f, max(2 * j * np.pi / L, 1e-12), (2 * j + 1) * np.pi / L, xtol=1e-15)
for j in range(L // 2)])
def U_modesum(L, ka, ds, lam):
k = even_labels(L, ka)
dl = np.arctan(2 * KH * np.sin(k) / ka)
N2 = 1 / (L / 2 + np.sin(2 * dl) * np.cos(k) / (2 * np.sin(k)))
w = N2 * np.cos(dl) / np.tan(k / 2) / (4 * KH * L * np.sin(k / 2) ** 2)
om = 2 * C0 * np.sin(k / 2)
return -(w * np.sin(k * ds + dl)) @ (1 - np.cos(om[:, None] * lam[None, :])), k, dl
def U_cont(L, alpha, tau):
return -(L / (2 * KH)) * (fext(alpha) - 0.5 * fext(alpha + tau) - 0.5 * fext(alpha - tau))
def ring_discrete(L, ka, lines):
tau = np.linspace(0, 8, 4001)
lam = tau * L / C0
out = []
for frac in (8, 4, 2):
ds = L // frac
w2, V = ring_eig(L, ds, ka)
U = U_eig(w2, V, L, lam)
Um, k, dl = U_modesum(L, ka, ds, lam)
if frac == 8:
ev = np.sort(w2[(np.abs(V[ds]) > 1e-9) & (w2 < 4 * KH)])
dq = np.max(np.abs(ev - 4 * KH * np.sin(k / 2) ** 2)) if len(ev) == len(k) else np.nan
dd = np.max(np.abs(k * L - (2 * np.arange(len(k)) + 1) * np.pi + 2 * dl))
nub = w2[-1]
out.append(f" L={L} kappa_a={ka}: {len(k)} even in-band labels from (9.19) vs {len(ev)} eigenvalues, "
f"max|nu_eig - 4 kappa_h sin^2(k_j/2)| = {dq:.1e}; max|k_j L - (2j+1)pi + 2 delta_j| = {dd:.1e}; "
f"top eigenvalue {nub:.6f} (line value 2kh+sqrt(4kh^2+ka^2) = {2*KH+np.sqrt(4*KH**2+ka**2):.6f})")
loc = np.max(np.abs(U - Um))
R = U - U_cont(L, ds / L, tau)
sc = np.max(np.abs(R) * ka / (1 + tau))
out.append(f" alpha=1/{frac}: max|U_eig - U_modesum| = {loc:.2e} (localized mode; bound 2/(kappa_a L) = {2/(ka*L):.2e}); "
f"max_(tau<=8)|R| = {np.max(np.abs(R)):.3f}, max_(tau<=2)|R| = {np.max(np.abs(R[tau <= 2])):.3f}, "
f"max_tau |R| kappa_a/(1+tau) = {sc:.3f}, |R(tau=8)| = {abs(R[-1]):.3f}, |U_st| = {1/ka + ds*(L-ds)/(2*KH*L):.2f}")
lines += out
def ring_stiff(L, ka, lines):
tau = np.linspace(0, 32, 16001)
lam = tau * L / C0
ds = L // 4
w2, V = ring_eig(L, ds, ka)
R = U_eig(w2, V, L, lam) - U_cont(L, 0.25, tau)
r8, r32 = np.max(np.abs(R[tau <= 8])), np.max(np.abs(R))
lines.append(f" stiff anchor kappa_a={ka:.0e}, L={L}, alpha=1/4: max_(tau<=8)|R| = {r8:.3f} "
f"(8^(2/3) L^(-1/3) = {8**(2/3)*L**(-1/3):.3f}, ratio {r8/(8**(2/3)*L**(-1/3)):.3f}); max_(tau<=32)|R| = {r32:.3f} "
f"(32^(2/3) L^(-1/3) = {32**(2/3)*L**(-1/3):.3f}, ratio {r32/(32**(2/3)*L**(-1/3)):.3f})")
def ring_weak(L, ka, lines, csv):
ds = L // 4
i = np.arange(L)
K = np.zeros((L, L))
K[i, i] = 2 * KH
K[i, (i + 1) % L] = -KH
K[i, (i - 1) % L] = -KH
K[ds, ds] += ka
w2, V = np.linalg.eigh(K)
cphi = 1 / ka + (L * L - 1) / (12 * KH * L)
TU = np.sqrt(L * cphi)
Ust = -(1 / ka + ds * (L - ds) / (2 * KH * L))
lam = np.linspace(0, 2.2 * np.pi * TU, 4001)
U = -((V.sum(0) / L) * V[0] / w2) @ (1 - np.cos(np.sqrt(w2)[:, None] * lam[None, :]))
pred = Ust * (1 - np.cos(lam / TU))
dev = np.max(np.abs(U - pred)) / abs(Ust)
lines.append(f" L={L} kappa_a={ka} (kappa_a L/kappa_h = {ka*L/KH:.3f}): T_U={TU:.1f} (wrap L/c0 = {L/C0:.0f}), "
f"1/sqrt(nu_0)/T_U = {1/np.sqrt(w2[0])/TU:.5f}, max|U-pred|/|U_st| = {dev:.4f}")
for j in range(0, len(lam), 8):
csv.append(f"1,ring_weak,{L},{ka},{lam[j]:.3f},{U[j]:.6e},{pred[j]:.6e}")
def main():
lines = ["Check (iii): the uniform part (item 3)", f"kappa_h = {KH}, c_0 = {C0}", "",
"CUBIC TORUS, kappa_a = 1, h - x_s = (2,1,0): prediction (9.14)"]
csv = ["D,case,L,kappa_a,lambda,U_exact,U_pred"]
for L in (8, 12, 16, 20):
cubic_case(L, (2, 1, 0), 1.0, lines, csv, "sizes")
lines.append("CUBIC TORUS, L = 12: anchor farther and stiffer anchor")
cubic_case(12, (5, 4, 3), 1.0, lines, csv, "far_anchor")
cubic_case(12, (2, 1, 0), 4.0, lines, csv, "kappa_a_4")
lines += ["", "RING, kappa_a = 1, d_s = L/4: prediction (9.15)"]
for L in (128, 256, 512):
ring_case(L, lines, csv)
lines += ["", "RING, discrete secular equation (9.19), exact mode sum and remainder of (9.15)"]
for L in (128, 256, 512):
for ka in (1.0, 4.0):
ring_discrete(L, ka, lines)
for L in (128, 256, 512):
ring_stiff(L, 1e4, lines)
lines += ["", "RING, weak anchor kappa_a L << kappa_h: slow-mode form (9.14) with c_phi^(L) = 1/kappa_a + (L^2-1)/(12 kappa_h L)"]
for L, ka in ((32, 0.004), (64, 0.002), (64, 0.0005)):
ring_weak(L, ka, lines, csv)
(OUT / "uniform_check.txt").write_text("\n".join(lines) + "\n", encoding="utf-8")
(OUT / "uniform_U.csv").write_text("\n".join(csv) + "\n", encoding="utf-8")
print("\n".join(lines))
if __name__ == "__main__":
main()
output/formation_Delta.csv340 KB · raw
output/formation_check.txt16 KB · raw
Check (ii): formation of the place-dependent part (item 2)
kappa_h = 1.0, kappa_a = 1.0, c_0 = 1.0
CUBIC TORUS (window below the wrap time; reference = L->infinity limit)
g_inf(0) = 0.2527310 (gamma_3/kappa_h = 0.2527310)
L=64 anchor far: h-x_s=(32, 32, 32), |h-x_s|=55.43; finite-L place average U(lambda_max)=-3.052e-03 (-lambda_max^2/2L^3 = -3.052e-03; -> 0 in order (O1))
rho= 3 b= 17.92: plateau(1.6<w<=min(3,w_max))=+0.119 pred(finite rho) +0.057 pred 1/(b-1) +0.059 wake max_(w>1.3)|D-pred| 0.318 at w=3.27 free-field wake max|u_0|/g_inf(rho e1): w>1.3 0.301 at w=3.27, 1.3<w<=2.5 0.251 min D(w<0.6)=-1.000 lambda-mean (O1) +0.057 w_max=13.3 c0*l_form/rho within window (eps=.2,.3,.5) = 3.82, 3.73, 1.00
rho= 4 b= 13.30: plateau(1.6<w<=min(3,w_max))=+0.086 pred(finite rho) +0.080 pred 1/(b-1) +0.081 wake max_(w>1.3)|D-pred| 0.291 at w=1.54 free-field wake max|u_0|/g_inf(rho e1): w>1.3 0.270 at w=1.54, 1.3<w<=2.5 0.270 min D(w<0.6)=-1.000 lambda-mean (O1) +0.079 w_max=10.0 c0*l_form/rho within window (eps=.2,.3,.5) = 3.55, 2.34, 1.04
rho= 6 b= 8.70: plateau(1.6<w<=min(3,w_max))=+0.113 pred(finite rho) +0.129 pred 1/(b-1) +0.130 wake max_(w>1.3)|D-pred| 0.249 at w=2.34 free-field wake max|u_0|/g_inf(rho e1): w>1.3 0.220 at w=2.34, 1.3<w<=2.5 0.220 min D(w<0.6)=-1.000 lambda-mean (O1) +0.129 w_max=6.7 c0*l_form/rho within window (eps=.2,.3,.5) = 3.49, 3.09, 1.03
rho= 8 b= 6.40: plateau(1.6<w<=min(3,w_max))=+0.176 pred(finite rho) +0.184 pred 1/(b-1) +0.185 wake max_(w>1.3)|D-pred| 0.279 at w=2.73 free-field wake max|u_0|/g_inf(rho e1): w>1.3 0.235 at w=2.73, 1.3<w<=2.5 0.227 min D(w<0.6)=-1.000 lambda-mean (O1) +0.184 w_max=5.0 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, 4.22, 1.02
rho=12 b= 4.12: plateau(1.6<w<=min(3,w_max))=+0.321 pred(finite rho) +0.319 pred 1/(b-1) +0.320 wake max_(w>1.3)|D-pred| 0.233 at w=2.35 free-field wake max|u_0|/g_inf(rho e1): w>1.3 0.177 at w=2.35, 1.3<w<=2.5 0.177 min D(w<0.6)=-1.000 lambda-mean (O1) +0.319 w_max=3.3 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, not formed, 3.12
rho=16 b= 3.00: plateau(1.6<w<=min(3,w_max))=+0.501 pred(finite rho) +0.499 pred 1/(b-1) +0.500 wake max_(w>1.3)|D-pred| 0.236 at w=1.35 free-field wake max|u_0|/g_inf(rho e1): w>1.3 0.157 at w=1.35, 1.3<w<=2.5 0.157 min D(w<0.6)=-1.000 lambda-mean (O1) +0.499 w_max=2.5 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, not formed, not formed
L=64 anchor near: h-x_s=(-8, 0, 0), |h-x_s|=8.00; finite-L place average U(lambda_max)=-3.036e-03 (-lambda_max^2/2L^3 = -3.052e-03; -> 0 in order (O1))
rho= 3 b= 3.67: plateau(1.6<w<=min(3,w_max))=+0.437 pred(finite rho) +0.357 pred 1/(b-1) +0.375 wake max_(w>1.3)|D-pred| 0.408 at w=3.27 image part max_(w>1.3)|I|/|ref| 0.004 (image front at w=6.33); |u_0|/|ref| at the wake maximum 0.408 min D(w<0.6)=-1.000 lambda-mean (O1) +0.354 w_max=13.3 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, not formed, 3.83
rho= 4 b= 3.00: plateau(1.6<w<=min(3,w_max))=+0.496 pred(finite rho) +0.487 pred 1/(b-1) +0.500 wake max_(w>1.3)|D-pred| 0.401 at w=1.54 image part max_(w>1.3)|I|/|ref| 0.005 (image front at w=5.00); |u_0|/|ref| at the wake maximum 0.401 min D(w<0.6)=-1.000 lambda-mean (O1) +0.483 w_max=10.0 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, not formed, 9.91
rho= 6 b= 2.33: plateau(1.6<w<=min(3,w_max))=+0.717 pred(finite rho) +0.742 pred 1/(b-1) +0.750 wake max_(w>1.3)|D-pred| 0.384 at w=2.34 image part max_(w>1.3)|I|/|ref| 0.007 (image front at w=3.67); |u_0|/|ref| at the wake maximum 0.384 min D(w<0.6)=-1.000 lambda-mean (O1) +0.736 w_max=6.7 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, not formed, not formed
rho= 8 b= 2.00: plateau(1.6<w<=min(3,w_max))=+0.980 pred(finite rho) +0.994 pred 1/(b-1) +1.000 wake max_(w>1.3)|D-pred| 0.469 at w=2.73 image part max_(w>1.3)|I|/|ref| 0.010 (image front at w=3.00); |u_0|/|ref| at the wake maximum 0.469 min D(w<0.6)=-1.000 lambda-mean (O1) +0.986 w_max=5.0 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, not formed, not formed
rho=12 b= 1.67: plateau(1.6<w<=min(3,w_max))=+1.494 pred(finite rho) +1.496 pred 1/(b-1) +1.500 wake max_(w>1.3)|D-pred| 0.447 at w=2.35 image part max_(w>1.3)|I|/|ref| 0.015 (image front at w=2.33); |u_0|/|ref| at the wake maximum 0.442 min D(w<0.6)=-1.000 lambda-mean (O1) +1.484 w_max=3.3 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, not formed, not formed
rho=16 b= 1.50: plateau(1.6<w<=min(3,w_max))=+1.991 pred(finite rho) +1.997 pred 1/(b-1) +2.000 wake max_(w>1.3)|D-pred| 0.472 at w=1.35 image part max_(w>1.3)|I|/|ref| 0.020 (image front at w=2.00); |u_0|/|ref| at the wake maximum 0.472 min D(w<0.6)=-1.000 lambda-mean (O1) +1.981 w_max=2.5 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, not formed, not formed
L=96 anchor far: h-x_s=(48, 48, 48), |h-x_s|=83.14; finite-L place average U(lambda_max)=-9.042e-04 (-lambda_max^2/2L^3 = -9.042e-04; -> 0 in order (O1))
rho= 3 b= 27.15: plateau(1.6<w<=min(3,w_max))=+0.097 pred(finite rho) +0.037 pred 1/(b-1) +0.038 wake max_(w>1.3)|D-pred| 0.312 at w=3.27 free-field wake max|u_0|/g_inf(rho e1): w>1.3 0.301 at w=3.27, 1.3<w<=2.5 0.251 min D(w<0.6)=-1.000 lambda-mean (O1) +0.037 w_max=13.3 c0*l_form/rho within window (eps=.2,.3,.5) = 3.80, 2.92, 1.00
rho= 4 b= 20.22: plateau(1.6<w<=min(3,w_max))=+0.058 pred(finite rho) +0.051 pred 1/(b-1) +0.052 wake max_(w>1.3)|D-pred| 0.284 at w=1.54 free-field wake max|u_0|/g_inf(rho e1): w>1.3 0.270 at w=1.54, 1.3<w<=2.5 0.270 min D(w<0.6)=-1.000 lambda-mean (O1) +0.051 w_max=10.0 c0*l_form/rho within window (eps=.2,.3,.5) = 2.98, 1.62, 1.05
rho= 6 b= 13.30: plateau(1.6<w<=min(3,w_max))=+0.065 pred(finite rho) +0.081 pred 1/(b-1) +0.081 wake max_(w>1.3)|D-pred| 0.238 at w=2.34 free-field wake max|u_0|/g_inf(rho e1): w>1.3 0.220 at w=2.34, 1.3<w<=2.5 0.220 min D(w<0.6)=-1.000 lambda-mean (O1) +0.081 w_max=6.7 c0*l_form/rho within window (eps=.2,.3,.5) = 3.47, 3.06, 1.04
rho= 8 b= 9.85: plateau(1.6<w<=min(3,w_max))=+0.105 pred(finite rho) +0.112 pred 1/(b-1) +0.113 wake max_(w>1.3)|D-pred| 0.262 at w=2.73 free-field wake max|u_0|/g_inf(rho e1): w>1.3 0.235 at w=2.73, 1.3<w<=2.5 0.227 min D(w<0.6)=-1.000 lambda-mean (O1) +0.112 w_max=5.0 c0*l_form/rho within window (eps=.2,.3,.5) = 4.24, 2.77, 1.03
rho=12 b= 6.40: plateau(1.6<w<=min(3,w_max))=+0.187 pred(finite rho) +0.185 pred 1/(b-1) +0.185 wake max_(w>1.3)|D-pred| 0.210 at w=2.35 free-field wake max|u_0|/g_inf(rho e1): w>1.3 0.177 at w=2.35, 1.3<w<=2.5 0.177 min D(w<0.6)=-1.000 lambda-mean (O1) +0.185 w_max=3.3 c0*l_form/rho within window (eps=.2,.3,.5) = 3.33, 3.14, 1.02
rho=16 b= 4.69: plateau(1.6<w<=min(3,w_max))=+0.272 pred(finite rho) +0.271 pred 1/(b-1) +0.271 wake max_(w>1.3)|D-pred| 0.200 at w=1.35 free-field wake max|u_0|/g_inf(rho e1): w>1.3 0.157 at w=1.35, 1.3<w<=2.5 0.157 min D(w<0.6)=-1.000 lambda-mean (O1) +0.270 w_max=2.5 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, not formed, 1.23
L=96 anchor near: h-x_s=(-8, 0, 0), |h-x_s|=8.00; finite-L place average U(lambda_max)=-8.996e-04 (-lambda_max^2/2L^3 = -9.042e-04; -> 0 in order (O1))
rho= 3 b= 3.67: plateau(1.6<w<=min(3,w_max))=+0.437 pred(finite rho) +0.357 pred 1/(b-1) +0.375 wake max_(w>1.3)|D-pred| 0.408 at w=3.27 image part max_(w>1.3)|I|/|ref| 0.004 (image front at w=6.33); |u_0|/|ref| at the wake maximum 0.408 min D(w<0.6)=-1.000 lambda-mean (O1) +0.354 w_max=13.3 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, not formed, 3.83
rho= 4 b= 3.00: plateau(1.6<w<=min(3,w_max))=+0.496 pred(finite rho) +0.487 pred 1/(b-1) +0.500 wake max_(w>1.3)|D-pred| 0.401 at w=1.54 image part max_(w>1.3)|I|/|ref| 0.005 (image front at w=5.00); |u_0|/|ref| at the wake maximum 0.401 min D(w<0.6)=-1.000 lambda-mean (O1) +0.483 w_max=10.0 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, not formed, 9.91
rho= 6 b= 2.33: plateau(1.6<w<=min(3,w_max))=+0.717 pred(finite rho) +0.742 pred 1/(b-1) +0.750 wake max_(w>1.3)|D-pred| 0.384 at w=2.34 image part max_(w>1.3)|I|/|ref| 0.007 (image front at w=3.67); |u_0|/|ref| at the wake maximum 0.384 min D(w<0.6)=-1.000 lambda-mean (O1) +0.736 w_max=6.7 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, not formed, not formed
rho= 8 b= 2.00: plateau(1.6<w<=min(3,w_max))=+0.980 pred(finite rho) +0.994 pred 1/(b-1) +1.000 wake max_(w>1.3)|D-pred| 0.469 at w=2.73 image part max_(w>1.3)|I|/|ref| 0.010 (image front at w=3.00); |u_0|/|ref| at the wake maximum 0.469 min D(w<0.6)=-1.000 lambda-mean (O1) +0.986 w_max=5.0 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, not formed, not formed
rho=12 b= 1.67: plateau(1.6<w<=min(3,w_max))=+1.494 pred(finite rho) +1.496 pred 1/(b-1) +1.500 wake max_(w>1.3)|D-pred| 0.447 at w=2.35 image part max_(w>1.3)|I|/|ref| 0.015 (image front at w=2.33); |u_0|/|ref| at the wake maximum 0.442 min D(w<0.6)=-1.000 lambda-mean (O1) +1.484 w_max=3.3 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, not formed, not formed
rho=16 b= 1.50: plateau(1.6<w<=min(3,w_max))=+1.991 pred(finite rho) +1.997 pred 1/(b-1) +2.000 wake max_(w>1.3)|D-pred| 0.472 at w=1.35 image part max_(w>1.3)|I|/|ref| 0.020 (image front at w=2.00); |u_0|/|ref| at the wake maximum 0.472 min D(w<0.6)=-1.000 lambda-mean (O1) +1.981 w_max=2.5 c0*l_form/rho within window (eps=.2,.3,.5) = not formed, not formed, not formed
size check (far): max_l,rho |da_64 - da_96| = 1.95e-07 (max |da| = 3.55e-02)
size check (near): max_l,rho |da_64 - da_96| = 1.93e-07 (max |da| = 3.55e-02)
time-step check (L=64 near): max |da(dt=0.1) - da(dt=0.05)| = 8.74e-06
ABOVE THRESHOLD: kappa_a = 8.0 > kappa_h/gamma_3 = 3.9568; L=64, near anchor h-x_s=(-8, 0, 0)
localized mode: nu_b from 1 + kappa_a gamma_b(0) = 0 on Z^3: 14.78961903; sparse eigen-solver (L=64): 14.78961903
rho= 3 b= 3.67: plateau=+0.437 pred(finite rho) +0.357 1/(b-1) +0.375 lambda-mean (O1) +0.348 (beta=0.0264) min D(w<0.6)=-1.000 c0*l_form/rho (eps=.5, window) = 3.83; max_window|I|*d_x*d_s^(1/6): kappa_a=1 0.0011, kappa_a=8 0.0037; localized share psi_b(x)psi_b(x_s)/nu_b = -9.45e-17
rho= 4 b= 3.00: plateau=+0.496 pred(finite rho) +0.487 1/(b-1) +0.500 lambda-mean (O1) +0.474 (beta=0.0264) min D(w<0.6)=-1.000 c0*l_form/rho (eps=.5, window) = 9.90; max_window|I|*d_x*d_s^(1/6): kappa_a=1 0.0011, kappa_a=8 0.0037; localized share psi_b(x)psi_b(x_s)/nu_b = +1.88e-17
rho= 6 b= 2.33: plateau=+0.717 pred(finite rho) +0.742 1/(b-1) +0.750 lambda-mean (O1) +0.722 (beta=0.0264) min D(w<0.6)=-1.000 c0*l_form/rho (eps=.5, window) = not formed; max_window|I|*d_x*d_s^(1/6): kappa_a=1 0.0011, kappa_a=8 0.0037; localized share psi_b(x)psi_b(x_s)/nu_b = +7.64e-19
rho= 8 b= 2.00: plateau=+0.980 pred(finite rho) +0.994 1/(b-1) +1.000 lambda-mean (O1) +0.968 (beta=0.0264) min D(w<0.6)=-1.000 c0*l_form/rho (eps=.5, window) = not formed; max_window|I|*d_x*d_s^(1/6): kappa_a=1 0.0011, kappa_a=8 0.0037; localized share psi_b(x)psi_b(x_s)/nu_b = +3.18e-20
rho=12 b= 1.67: plateau=+1.480 pred(finite rho) +1.496 1/(b-1) +1.500 lambda-mean (O1) +1.456 (beta=0.0264) min D(w<0.6)=-1.000 c0*l_form/rho (eps=.5, window) = not formed; max_window|I|*d_x*d_s^(1/6): kappa_a=1 0.0011, kappa_a=8 0.0037; localized share psi_b(x)psi_b(x_s)/nu_b = +5.76e-23
rho=16 b= 1.50: plateau=+1.970 pred(finite rho) +1.997 1/(b-1) +2.000 lambda-mean (O1) +1.944 (beta=0.0264) min D(w<0.6)=-1.000 c0*l_form/rho (eps=.5, window) = not formed; max_window|I|*d_x*d_s^(1/6): kappa_a=1 0.0011, kappa_a=8 0.0037; localized share psi_b(x)psi_b(x_s)/nu_b = +4.01e-25
RING (exact eigen-decomposition; window below the wrap time; reference grad = 1/kappa_h)
(a) anchor at d_s = 2 rho (a = 2): prediction -1 (w<1), -1/2 (1<w<3), 0 (w>3)
L=1024 rho= 20: plateau(1.4<w<2.6)=-0.496 plateau(w>3.4)=-0.007 max|D|(w>1.4)=0.577 D across image front (2.6<w<3.4) in [-0.544, +0.069] c0*l_form/rho: eps=0.25 -> 3.08 (pred 3), eps=0.6 -> 1.34 (pred 1); finite-L ref grad*kappa_h = 0.961
L=1024 rho= 40: plateau(1.4<w<2.6)=-0.501 plateau(w>3.4)=-0.002 max|D|(w>1.4)=0.561 D across image front (2.6<w<3.4) in [-0.531, +0.097] c0*l_form/rho: eps=0.25 -> 3.06 (pred 3), eps=0.6 -> 1.20 (pred 1); finite-L ref grad*kappa_h = 0.922
L=1024 rho= 80: plateau(1.4<w<2.6)=-0.500 plateau(w>3.4)=+0.000 max|D|(w>1.4)=0.543 D across image front (2.6<w<3.4) in [-0.522, +0.098] c0*l_form/rho: eps=0.25 -> 3.03 (pred 3), eps=0.6 -> 1.12 (pred 1); finite-L ref grad*kappa_h = 0.844
L=1024 rho=160: plateau(1.4<w<2.6)=-0.500 plateau(w>3.4)=-0.000 max|D|(w>1.4)=0.530 D across image front (2.6<w<3.4) in [-0.515, +0.117] c0*l_form/rho: eps=0.25 -> 3.02 (pred 3), eps=0.6 -> 1.02 (pred 1); finite-L ref grad*kappa_h = 0.688
L=2048 rho= 20: plateau(1.4<w<2.6)=-0.496 plateau(w>3.4)=-0.007 max|D|(w>1.4)=0.577 D across image front (2.6<w<3.4) in [-0.544, +0.069] c0*l_form/rho: eps=0.25 -> 3.08 (pred 3), eps=0.6 -> 1.34 (pred 1); finite-L ref grad*kappa_h = 0.980
L=2048 rho= 40: plateau(1.4<w<2.6)=-0.501 plateau(w>3.4)=-0.002 max|D|(w>1.4)=0.561 D across image front (2.6<w<3.4) in [-0.531, +0.097] c0*l_form/rho: eps=0.25 -> 3.06 (pred 3), eps=0.6 -> 1.20 (pred 1); finite-L ref grad*kappa_h = 0.961
L=2048 rho= 80: plateau(1.4<w<2.6)=-0.500 plateau(w>3.4)=+0.000 max|D|(w>1.4)=0.543 D across image front (2.6<w<3.4) in [-0.522, +0.098] c0*l_form/rho: eps=0.25 -> 3.03 (pred 3), eps=0.6 -> 1.12 (pred 1); finite-L ref grad*kappa_h = 0.922
L=2048 rho=160: plateau(1.4<w<2.6)=-0.500 plateau(w>3.4)=-0.000 max|D|(w>1.4)=0.530 D across image front (2.6<w<3.4) in [-0.515, +0.117] c0*l_form/rho: eps=0.25 -> 3.02 (pred 3), eps=0.6 -> 1.02 (pred 1); finite-L ref grad*kappa_h = 0.844
(b) anchor far, d_s = 400: prediction -1/2 until c0*l = 2 d_s - rho, then 0
L=1024 rho= 10: mean D on 1.4 rho < c0 l < 0.9(2d_s-rho) = -0.500, max|D+1/2| there = 0.111; c0*l_form(eps=0.25) = 793.2 (pred 2d_s-rho = 790), c0*l_form/rho = 79.3
L=1024 rho= 20: mean D on 1.4 rho < c0 l < 0.9(2d_s-rho) = -0.500, max|D+1/2| there = 0.087; c0*l_form(eps=0.25) = 783.2 (pred 2d_s-rho = 780), c0*l_form/rho = 39.2
L=1024 rho= 40: mean D on 1.4 rho < c0 l < 0.9(2d_s-rho) = -0.500, max|D+1/2| there = 0.062; c0*l_form(eps=0.25) = 763.5 (pred 2d_s-rho = 760), c0*l_form/rho = 19.1
L=1024 rho= 80: mean D on 1.4 rho < c0 l < 0.9(2d_s-rho) = -0.500, max|D+1/2| there = 0.045; c0*l_form(eps=0.25) = 723.0 (pred 2d_s-rho = 720), c0*l_form/rho = 9.0
L=2048 rho= 10: mean D on 1.4 rho < c0 l < 0.9(2d_s-rho) = -0.500, max|D+1/2| there = 0.111; c0*l_form(eps=0.25) = 793.2 (pred 2d_s-rho = 790), c0*l_form/rho = 79.3
L=2048 rho= 20: mean D on 1.4 rho < c0 l < 0.9(2d_s-rho) = -0.500, max|D+1/2| there = 0.087; c0*l_form(eps=0.25) = 783.2 (pred 2d_s-rho = 780), c0*l_form/rho = 39.2
L=2048 rho= 40: mean D on 1.4 rho < c0 l < 0.9(2d_s-rho) = -0.500, max|D+1/2| there = 0.062; c0*l_form(eps=0.25) = 763.5 (pred 2d_s-rho = 760), c0*l_form/rho = 19.1
L=2048 rho= 80: mean D on 1.4 rho < c0 l < 0.9(2d_s-rho) = -0.500, max|D+1/2| there = 0.045; c0*l_form(eps=0.25) = 723.0 (pred 2d_s-rho = 720), c0*l_form/rho = 9.0
(c) uniform form of Step 8: sup_window kappa_h|grad da - F(rho) - F(m)|, a = 2, m = 3 rho - 1, c0 l <= 8 rho, L = 2048
kappa_a=0.25: rho=20: 0.4402; rho=40: 0.4055; rho=80: 0.3688; rho=160: 0.3305; localized eigenvalue 4.015564 vs 2kh+sqrt(4kh^2+ka^2) = 4.015564
kappa_a=1.0: rho=20: 0.2152; rho=40: 0.1817; rho=80: 0.1515; rho=160: 0.1248; localized eigenvalue 4.236068 vs 2kh+sqrt(4kh^2+ka^2) = 4.236068
kappa_a=4.0: rho=20: 0.0648; rho=40: 0.0524; rho=80: 0.0421; rho=160: 0.0337; localized eigenvalue 6.472136 vs 2kh+sqrt(4kh^2+ka^2) = 6.472136
Bessel form of F (n = rho and m) against the free-ring gradient, max kappa_h|diff|: rho=20: 1.2e-07; rho=40: 8.3e-08; rho=80: 5.9e-08; rho=160: 4.2e-08
output/made_field_check.txt4 KB · raw
Check (i): made field of item 1 against the exact solution of the means
kappa_h = 1.0, kappa_a = 1.0, c_0 = sqrt(kappa_h) = 1.0
D=1 L=64 xs=(0,) h=(3,) |h-xs|=3.000 lambda_max=192 dt=0.01
place (-6,): max|derived(9.4)-exact|/max|exact| = 4.17e-06
place (4,): max|derived(9.4)-exact|/max|exact| = 7.55e-06
place (11,): max|derived(9.4)-exact|/max|exact| = 4.98e-06
uniform part (9.5) force formula vs exact place average: rel. err = 3.81e-06
static split (9.6): rel. err = 1.97e-14; static image strength -kappa_a*da_st(h) = 1.000000000000
D=1 L=64 xs=(0,) h=(17,) |h-xs|=17.000 lambda_max=192 dt=0.01
place (-6,): max|derived(9.4)-exact|/max|exact| = 1.72e-06
place (4,): max|derived(9.4)-exact|/max|exact| = 1.58e-06
place (11,): max|derived(9.4)-exact|/max|exact| = 2.79e-06
uniform part (9.5) force formula vs exact place average: rel. err = 1.78e-06
static split (9.6): rel. err = 4.09e-14; static image strength -kappa_a*da_st(h) = 1.000000000000
D=1 L=96 xs=(0,) h=(3,) |h-xs|=3.000 lambda_max=288 dt=0.01
place (-6,): max|derived(9.4)-exact|/max|exact| = 4.28e-06
place (4,): max|derived(9.4)-exact|/max|exact| = 7.70e-06
place (11,): max|derived(9.4)-exact|/max|exact| = 4.96e-06
uniform part (9.5) force formula vs exact place average: rel. err = 3.65e-06
static split (9.6): rel. err = 1.27e-13; static image strength -kappa_a*da_st(h) = 1.000000000000
D=1 L=96 xs=(0,) h=(17,) |h-xs|=17.000 lambda_max=288 dt=0.01
place (-6,): max|derived(9.4)-exact|/max|exact| = 1.25e-06
place (4,): max|derived(9.4)-exact|/max|exact| = 1.57e-06
place (11,): max|derived(9.4)-exact|/max|exact| = 2.76e-06
uniform part (9.5) force formula vs exact place average: rel. err = 1.42e-06
static split (9.6): rel. err = 5.19e-14; static image strength -kappa_a*da_st(h) = 1.000000000000
D=3 L=8 xs=(0, 0, 0) h=(1, 1, 0) |h-xs|=1.414 lambda_max=24 dt=0.005
place (2, 0, 0): max|derived(9.4)-exact|/max|exact| = 2.13e-07
place (1, -2, 1): max|derived(9.4)-exact|/max|exact| = 1.75e-07
place (-3, 1, 2): max|derived(9.4)-exact|/max|exact| = 1.75e-07
uniform part (9.5) force formula vs exact place average: rel. err = 1.69e-07
static split (9.6): rel. err = 3.14e-13; static image strength -kappa_a*da_st(h) = 1.000000000000
D=3 L=8 xs=(0, 0, 0) h=(3, 2, 1) |h-xs|=3.742 lambda_max=24 dt=0.005
place (2, 0, 0): max|derived(9.4)-exact|/max|exact| = 1.38e-07
place (1, -2, 1): max|derived(9.4)-exact|/max|exact| = 1.24e-07
place (-3, 1, 2): max|derived(9.4)-exact|/max|exact| = 1.37e-07
uniform part (9.5) force formula vs exact place average: rel. err = 1.25e-07
static split (9.6): rel. err = 1.68e-12; static image strength -kappa_a*da_st(h) = 0.999999999998
D=3 L=12 xs=(0, 0, 0) h=(1, 1, 0) |h-xs|=1.414 lambda_max=36 dt=0.005
place (2, 0, 0): max|derived(9.4)-exact|/max|exact| = 1.71e-07
place (1, -2, 1): max|derived(9.4)-exact|/max|exact| = 1.47e-07
place (-3, 1, 2): max|derived(9.4)-exact|/max|exact| = 1.41e-07
uniform part (9.5) force formula vs exact place average: rel. err = 1.37e-07
static split (9.6): rel. err = 3.30e-12; static image strength -kappa_a*da_st(h) = 0.999999999996
D=3 L=12 xs=(0, 0, 0) h=(3, 2, 1) |h-xs|=3.742 lambda_max=36 dt=0.005
place (2, 0, 0): max|derived(9.4)-exact|/max|exact| = 1.06e-07
place (1, -2, 1): max|derived(9.4)-exact|/max|exact| = 9.34e-08
place (-3, 1, 2): max|derived(9.4)-exact|/max|exact| = 9.09e-08
uniform part (9.5) force formula vs exact place average: rel. err = 9.06e-08
static split (9.6): rel. err = 2.51e-13; static image strength -kappa_a*da_st(h) = 1.000000000000
largest relative deviation over all cases: 7.70e-06
output/uniform_U.csv275 KB · raw
output/uniform_check.txt9 KB · raw
Check (iii): the uniform part (item 3)
kappa_h = 1.0, c_0 = 1.0
CUBIC TORUS, kappa_a = 1, h - x_s = (2,1,0): prediction (9.14)
L= 8 kappa_a=1.0 h-x_s=(2, 1, 0): c_phi^(L)=1.22461 U_st=-1.21518 T_U=25.04; 1/sqrt(nu_0)=25.04 (ratio 1.00006); max|U-pred|/|U_st|=3.24e-04 (x L = 0.0026); first minimum at l/(pi T_U)=0.9998, U_min/U_st=2.0000; U(2 pi T_U)/U_st=+0.0000; slow image: max|A0perp+(1-beta)g(x-h)|/max|.|=5.26e-04 (x L^2 = 0.03), beta_L=0.0077
L=12 kappa_a=1.0 h-x_s=(2, 1, 0): c_phi^(L)=1.23394 U_st=-1.21632 T_U=46.18; 1/sqrt(nu_0)=46.18 (ratio 1.00003); max|U-pred|/|U_st|=1.53e-04 (x L = 0.0018); first minimum at l/(pi T_U)=0.9999, U_min/U_st=2.0000; U(2 pi T_U)/U_st=+0.0000; slow image: max|A0perp+(1-beta)g(x-h)|/max|.|=2.04e-04 (x L^2 = 0.03), beta_L=0.0143
L=16 kappa_a=1.0 h-x_s=(2, 1, 0): c_phi^(L)=1.23863 U_st=-1.21660 T_U=71.23; 1/sqrt(nu_0)=71.23 (ratio 1.00001); max|U-pred|/|U_st|=8.33e-05 (x L = 0.0013); first minimum at l/(pi T_U)=0.9999, U_min/U_st=2.0000; U(2 pi T_U)/U_st=+0.0000; slow image: max|A0perp+(1-beta)g(x-h)|/max|.|=1.06e-04 (x L^2 = 0.03), beta_L=0.0178
L=20 kappa_a=1.0 h-x_s=(2, 1, 0): c_phi^(L)=1.24145 U_st=-1.21670 T_U=99.66; 1/sqrt(nu_0)=99.66 (ratio 1.00001); max|U-pred|/|U_st|=5.45e-05 (x L = 0.0011); first minimum at l/(pi T_U)=1.0001, U_min/U_st=2.0000; U(2 pi T_U)/U_st=+0.0000; slow image: max|A0perp+(1-beta)g(x-h)|/max|.|=6.75e-05 (x L^2 = 0.03), beta_L=0.0199
CUBIC TORUS, L = 12: anchor farther and stiffer anchor
L=12 kappa_a=1.0 h-x_s=(5, 4, 3): c_phi^(L)=1.23394 U_st=-1.23791 T_U=46.18; 1/sqrt(nu_0)=46.18 (ratio 1.00003); max|U-pred|/|U_st|=1.33e-04 (x L = 0.0016); first minimum at l/(pi T_U)=0.9999, U_min/U_st=2.0000; U(2 pi T_U)/U_st=-0.0000; slow image: max|A0perp+(1-beta)g(x-h)|/max|.|=2.31e-04 (x L^2 = 0.03), beta_L=-0.0032
L=12 kappa_a=4.0 h-x_s=(2, 1, 0): c_phi^(L)=0.48394 U_st=-0.46632 T_U=28.92; 1/sqrt(nu_0)=28.92 (ratio 1.00017); max|U-pred|/|U_st|=9.64e-04 (x L = 0.0116); first minimum at l/(pi T_U)=1.0000, U_min/U_st=1.9999; U(2 pi T_U)/U_st=+0.0001; slow image: max|A0perp+(1-beta)g(x-h)|/max|.|=3.71e-04 (x L^2 = 0.05), beta_L=0.0364
RING, kappa_a = 1, d_s = L/4: prediction (9.15)
L=128 d_s=32: U_st=-13.000; max|U-U_cont|/|U_st| over tau<=2: 0.2308 (x L = 29.55), over tau<=4: 0.3846; first minimum at c0 l/L = 1.0320 (pred 1), U_min/U_st=2.0000 (pred 2); U(c0 l = 2L)/U_st = +0.0193 (pred 0)
L=256 d_s=64: U_st=-25.000; max|U-U_cont|/|U_st| over tau<=2: 0.1201 (x L = 30.75), over tau<=4: 0.2002; first minimum at c0 l/L = 1.0160 (pred 1), U_min/U_st=2.0000 (pred 2); U(c0 l = 2L)/U_st = +0.0050 (pred 0)
L=512 d_s=128: U_st=-49.000; max|U-U_cont|/|U_st| over tau<=2: 0.0613 (x L = 31.38), over tau<=4: 0.1022; first minimum at c0 l/L = 1.0080 (pred 1), U_min/U_st=2.0000 (pred 2); U(c0 l = 2L)/U_st = +0.0013 (pred 0)
RING, discrete secular equation (9.19), exact mode sum and remainder of (9.15)
L=128 kappa_a=1.0: 64 even in-band labels from (9.19) vs 64 eigenvalues, max|nu_eig - 4 kappa_h sin^2(k_j/2)| = 2.2e-15; max|k_j L - (2j+1)pi + 2 delta_j| = 9.6e-14; top eigenvalue 4.236068 (line value 2kh+sqrt(4kh^2+ka^2) = 4.236068)
alpha=1/8: max|U_eig - U_modesum| = 1.76e-07 (localized mode; bound 2/(kappa_a L) = 1.56e-02); max_(tau<=8)|R| = 5.499, max_(tau<=2)|R| = 2.500, max_tau |R| kappa_a/(1+tau) = 1.161, |R(tau=8)| = 3.219, |U_st| = 8.00
alpha=1/4: max|U_eig - U_modesum| = 1.43e-10 (localized mode; bound 2/(kappa_a L) = 1.56e-02); max_(tau<=8)|R| = 9.000, max_(tau<=2)|R| = 3.001, max_tau |R| kappa_a/(1+tau) = 1.318, |R(tau=8)| = 4.004, |U_st| = 13.00
alpha=1/2: max|U_eig - U_modesum| = 1.28e-10 (localized mode; bound 2/(kappa_a L) = 1.56e-02); max_(tau<=8)|R| = 14.469, max_(tau<=2)|R| = 3.968, max_tau |R| kappa_a/(1+tau) = 1.706, |R(tau=8)| = 4.000, |U_st| = 17.00
L=128 kappa_a=4.0: 64 even in-band labels from (9.19) vs 64 eigenvalues, max|nu_eig - 4 kappa_h sin^2(k_j/2)| = 1.8e-15; max|k_j L - (2j+1)pi + 2 delta_j| = 7.9e-14; top eigenvalue 6.472136 (line value 2kh+sqrt(4kh^2+ka^2) = 6.472136)
alpha=1/8: max|U_eig - U_modesum| = 2.17e-11 (localized mode; bound 2/(kappa_a L) = 3.91e-03); max_(tau<=8)|R| = 1.407, max_(tau<=2)|R| = 0.633, max_tau |R| kappa_a/(1+tau) = 1.186, |R(tau=8)| = 0.253, |U_st| = 7.25
alpha=1/4: max|U_eig - U_modesum| = 4.25e-11 (localized mode; bound 2/(kappa_a L) = 3.91e-03); max_(tau<=8)|R| = 2.286, max_(tau<=2)|R| = 0.759, max_tau |R| kappa_a/(1+tau) = 1.343, |R(tau=8)| = 0.251, |U_st| = 12.25
alpha=1/2: max|U_eig - U_modesum| = 1.20e-10 (localized mode; bound 2/(kappa_a L) = 3.91e-03); max_(tau<=8)|R| = 3.923, max_(tau<=2)|R| = 1.014, max_tau |R| kappa_a/(1+tau) = 1.840, |R(tau=8)| = 0.252, |U_st| = 16.25
L=256 kappa_a=1.0: 128 even in-band labels from (9.19) vs 128 eigenvalues, max|nu_eig - 4 kappa_h sin^2(k_j/2)| = 1.8e-15; max|k_j L - (2j+1)pi + 2 delta_j| = 1.9e-13; top eigenvalue 4.236068 (line value 2kh+sqrt(4kh^2+ka^2) = 4.236068)
alpha=1/8: max|U_eig - U_modesum| = 8.36e-11 (localized mode; bound 2/(kappa_a L) = 7.81e-03); max_(tau<=8)|R| = 5.511, max_(tau<=2)|R| = 2.502, max_tau |R| kappa_a/(1+tau) = 1.170, |R(tau=8)| = 2.019, |U_st| = 15.00
alpha=1/4: max|U_eig - U_modesum| = 5.19e-10 (localized mode; bound 2/(kappa_a L) = 7.81e-03); max_(tau<=8)|R| = 9.011, max_(tau<=2)|R| = 3.003, max_tau |R| kappa_a/(1+tau) = 1.326, |R(tau=8)| = 2.000, |U_st| = 25.00
alpha=1/2: max|U_eig - U_modesum| = 8.34e-10 (localized mode; bound 2/(kappa_a L) = 7.81e-03); max_(tau<=8)|R| = 15.234, max_(tau<=2)|R| = 3.987, max_tau |R| kappa_a/(1+tau) = 1.781, |R(tau=8)| = 1.999, |U_st| = 33.00
L=256 kappa_a=4.0: 128 even in-band labels from (9.19) vs 128 eigenvalues, max|nu_eig - 4 kappa_h sin^2(k_j/2)| = 1.3e-15; max|k_j L - (2j+1)pi + 2 delta_j| = 2.9e-13; top eigenvalue 6.472136 (line value 2kh+sqrt(4kh^2+ka^2) = 6.472136)
alpha=1/8: max|U_eig - U_modesum| = 3.48e-10 (localized mode; bound 2/(kappa_a L) = 1.95e-03); max_(tau<=8)|R| = 1.403, max_(tau<=2)|R| = 0.633, max_tau |R| kappa_a/(1+tau) = 1.187, |R(tau=8)| = 0.130, |U_st| = 14.25
alpha=1/4: max|U_eig - U_modesum| = 3.00e-10 (localized mode; bound 2/(kappa_a L) = 1.95e-03); max_(tau<=8)|R| = 2.280, max_(tau<=2)|R| = 0.758, max_tau |R| kappa_a/(1+tau) = 1.343, |R(tau=8)| = 0.127, |U_st| = 24.25
alpha=1/2: max|U_eig - U_modesum| = 8.15e-10 (localized mode; bound 2/(kappa_a L) = 1.95e-03); max_(tau<=8)|R| = 3.974, max_(tau<=2)|R| = 1.014, max_tau |R| kappa_a/(1+tau) = 1.866, |R(tau=8)| = 0.126, |U_st| = 32.25
L=512 kappa_a=1.0: 256 even in-band labels from (9.19) vs 256 eigenvalues, max|nu_eig - 4 kappa_h sin^2(k_j/2)| = 1.8e-15; max|k_j L - (2j+1)pi + 2 delta_j| = 5.1e-13; top eigenvalue 4.236068 (line value 2kh+sqrt(4kh^2+ka^2) = 4.236068)
alpha=1/8: max|U_eig - U_modesum| = 1.86e-09 (localized mode; bound 2/(kappa_a L) = 3.91e-03); max_(tau<=8)|R| = 5.514, max_(tau<=2)|R| = 2.503, max_tau |R| kappa_a/(1+tau) = 1.174, |R(tau=8)| = 1.002, |U_st| = 29.00
alpha=1/4: max|U_eig - U_modesum| = 2.64e-09 (localized mode; bound 2/(kappa_a L) = 3.91e-03); max_(tau<=8)|R| = 9.014, max_(tau<=2)|R| = 3.003, max_tau |R| kappa_a/(1+tau) = 1.330, |R(tau=8)| = 1.000, |U_st| = 49.00
alpha=1/2: max|U_eig - U_modesum| = 4.49e-09 (localized mode; bound 2/(kappa_a L) = 3.91e-03); max_(tau<=8)|R| = 15.618, max_(tau<=2)|R| = 3.996, max_tau |R| kappa_a/(1+tau) = 1.832, |R(tau=8)| = 1.000, |U_st| = 65.00
L=512 kappa_a=4.0: 256 even in-band labels from (9.19) vs 256 eigenvalues, max|nu_eig - 4 kappa_h sin^2(k_j/2)| = 1.8e-15; max|k_j L - (2j+1)pi + 2 delta_j| = 4.7e-13; top eigenvalue 6.472136 (line value 2kh+sqrt(4kh^2+ka^2) = 6.472136)
alpha=1/8: max|U_eig - U_modesum| = 1.97e-09 (localized mode; bound 2/(kappa_a L) = 9.77e-04); max_(tau<=8)|R| = 1.399, max_(tau<=2)|R| = 0.631, max_tau |R| kappa_a/(1+tau) = 1.186, |R(tau=8)| = 0.064, |U_st| = 28.25
alpha=1/4: max|U_eig - U_modesum| = 1.95e-09 (localized mode; bound 2/(kappa_a L) = 9.77e-04); max_(tau<=8)|R| = 2.274, max_(tau<=2)|R| = 0.757, max_tau |R| kappa_a/(1+tau) = 1.343, |R(tau=8)| = 0.062, |U_st| = 48.25
alpha=1/2: max|U_eig - U_modesum| = 4.40e-09 (localized mode; bound 2/(kappa_a L) = 9.77e-04); max_(tau<=8)|R| = 3.998, max_(tau<=2)|R| = 1.012, max_tau |R| kappa_a/(1+tau) = 1.879, |R(tau=8)| = 0.062, |U_st| = 64.25
stiff anchor kappa_a=1e+04, L=128, alpha=1/4: max_(tau<=8)|R| = 0.046 (8^(2/3) L^(-1/3) = 0.794, ratio 0.058); max_(tau<=32)|R| = 0.121 (32^(2/3) L^(-1/3) = 2.000, ratio 0.060)
stiff anchor kappa_a=1e+04, L=256, alpha=1/4: max_(tau<=8)|R| = 0.037 (8^(2/3) L^(-1/3) = 0.630, ratio 0.058); max_(tau<=32)|R| = 0.096 (32^(2/3) L^(-1/3) = 1.587, ratio 0.060)
stiff anchor kappa_a=1e+04, L=512, alpha=1/4: max_(tau<=8)|R| = 0.029 (8^(2/3) L^(-1/3) = 0.500, ratio 0.059); max_(tau<=32)|R| = 0.077 (32^(2/3) L^(-1/3) = 1.260, ratio 0.061)
RING, weak anchor kappa_a L << kappa_h: slow-mode form (9.14) with c_phi^(L) = 1/kappa_a + (L^2-1)/(12 kappa_h L)
L=32 kappa_a=0.004 (kappa_a L/kappa_h = 0.128): T_U=89.9 (wrap L/c0 = 32), 1/sqrt(nu_0)/T_U = 1.00001, max|U-pred|/|U_st| = 0.0001
L=64 kappa_a=0.002 (kappa_a L/kappa_h = 0.128): T_U=179.8 (wrap L/c0 = 64), 1/sqrt(nu_0)/T_U = 1.00001, max|U-pred|/|U_st| = 0.0001
L=64 kappa_a=0.0005 (kappa_a L/kappa_h = 0.032): T_U=358.2 (wrap L/c0 = 64), 1/sqrt(nu_0)/T_U = 1.00000, max|U-pred|/|U_st| = 0.0000
nDot Physics Research