import json import math from pathlib import Path import numpy as np from scipy.linalg import eig, expm from scipy.integrate import solve_ivp # Reproducible toy study of the Fourier Hill construction. np.random.seed(7) T = 2.0 * np.pi omega = 2.0 * np.pi / T def rot(theta): c, s = np.cos(theta), np.sin(theta) return np.array([[c, -s], [s, c]]) def J_of_t(t, gamma=0.0, amp=1.0, rate=1.0): # A genuinely noncommuting periodic linearization. The constant gamma # shifts every Floquet exponent exactly, which gives a sharp prediction. R = rot(rate * t) D = np.diag([gamma + 0.18, gamma - 0.42]) return R @ D @ R.T + amp * np.array([[0.0, np.sin(t)], [0.25*np.cos(t), 0.0]]) def monodromy_exponents(gamma=0.0, amp=1.0, rate=1.0, rtol=2e-10, atol=2e-12): d = 2 def rhs(t, y): Y = y.reshape(d, d) return (J_of_t(t, gamma, amp, rate) @ Y).reshape(-1) sol = solve_ivp(rhs, (0.0, T), np.eye(d).reshape(-1), rtol=rtol, atol=atol, method="DOP853") vals = np.linalg.eigvals(sol.y[:, -1].reshape(d, d)) # Real parts are unambiguous for the dominant multiplier in this toy. mu = np.log(np.abs(vals)) / T return np.sort(mu)[::-1], vals def fourier_coeffs(gamma, amp, rate, M): ts = np.arange(M) * T / M samples = np.stack([J_of_t(t, gamma, amp, rate) for t in ts], axis=0) # np.fft convention: samples[m] = sum_k coeff[k] exp(+i 2pi km/M) coeff = np.fft.fft(samples, axis=0) / M return coeff def hill_exponents(gamma=0.0, amp=1.0, rate=1.0, N=8, M=None, filtered=False): if M is None: M = max(4*N + 1, 64) coeff = fourier_coeffs(gamma, amp, rate, M) modes = np.arange(-N, N+1) d = 2 H0 = np.zeros((len(modes)*d, len(modes)*d), dtype=complex) for ii, n in enumerate(modes): for jj, ell in enumerate(modes): k = (n - ell) % M block = coeff[k].copy() if n == ell: block = block - 1j * ell * omega * np.eye(d) H0[ii*d:(ii+1)*d, jj*d:(jj+1)*d] = block vals, vecs = np.linalg.eig(H0) if filtered: # Finite Hill matrices contain edge-localized spectral pollution. # Physical Fourier modes have small mass in the outermost shell. q = d edge = np.r_[np.arange(q), np.arange(len(modes)*q-q, len(modes)*q)] mass = np.sum(np.abs(vecs[edge, :])**2, axis=0) / np.sum(np.abs(vecs)**2, axis=0) keep = mass < 0.20 if not np.any(keep): keep = np.ones(len(vals), dtype=bool) vals = vals[keep] return np.sort(vals.real)[::-1], vals def penalty(gamma, beta=2.0, N=8): dominant = hill_exponents(gamma=gamma, N=N, filtered=True)[0][0] # Smooth version of max(0, mu+epsilon)^2; here epsilon=0. return beta * max(0.0, dominant)**2, dominant def regularized_parameter_fit(beta=2.0, steps=120, lr=0.08, N=8): # Task loss prefers gamma=+0.20 (unstable), Hill penalty should move it # to the stability boundary or below it. gamma = 0.20 for _ in range(steps): mu = hill_exponents(gamma=gamma, N=N, filtered=True)[0][0] # Numerically, this construction has d(mu)/d(gamma) approximately 1. dpen = 2.0 * beta * max(0.0, mu) grad = 2.0 * (gamma - 0.20) + dpen gamma -= lr * grad return gamma, hill_exponents(gamma=gamma, N=N, filtered=True)[0][0] def main(): report = {"seed": 7, "T": T, "omega": omega} # Prediction 1: adding gamma I shifts every exponent by gamma, slope 1. gs = np.linspace(-0.45, 0.35, 9) exact = np.array([monodromy_exponents(gamma=g)[0][0] for g in gs]) hill = np.array([hill_exponents(gamma=g, N=8, filtered=True)[0][0] for g in gs]) slope_exact, intercept_exact = np.polyfit(gs, exact, 1) slope_hill, intercept_hill = np.polyfit(gs, hill, 1) report["prediction_1_offset_scaling"] = { "predicted": "dominant Floquet exponent shifts as mu(gamma)=mu(0)+gamma; slope=1", "spectrum_filter": "discard eigenvectors with >20% mass in outer Fourier shell", "observed_exact_slope": float(slope_exact), "observed_hill_slope": float(slope_hill), "max_exact_shift_error": float(np.max(np.abs((exact-exact[4])-(gs-gs[4])))), "max_hill_vs_exact": float(np.max(np.abs(hill-exact))), } # Prediction 2: stability boundary occurs where dominant exponent crosses 0. # Estimate by linear interpolation of the sweep and compare to exact root. def crossing(x, y): for i in range(len(x)-1): if y[i] * y[i+1] <= 0: return float(x[i] - y[i]*(x[i+1]-x[i])/(y[i+1]-y[i])) return float("nan") report["prediction_2_boundary"] = { "predicted": "rho_F=1 exactly when max Re(mu)=0", "exact_crossing_gamma": crossing(gs, exact), "hill_N8_crossing_gamma": crossing(gs, hill), "exact_rho_at_crossing": float(np.exp(T * 0.0)), } # Prediction 3: increasing Fourier truncation should converge to monodromy. exact0 = monodromy_exponents()[0][0] trunc = {} for N in [1, 2, 4, 8, 12]: h = hill_exponents(N=N, filtered=True)[0][0] trunc[str(N)] = {"hill_mu": float(h), "abs_error": float(abs(h-exact0))} report["prediction_3_truncation"] = { "predicted": "smooth periodic coefficients give improving Hill approximation; N=8 close to N=4", "spectrum_filter": "discard edge-localized finite-section roots", "monodromy_mu": float(exact0), "values": trunc, "N4_to_N8_relative_change": float(abs(trunc["8"]["hill_mu"]-trunc["4"]["hill_mu"]) / max(abs(trunc["8"]["hill_mu"]), 1e-12)), } # Secondary mini experiment: task-only fit remains at the unstable target; # regularized fit should sacrifice task objective to reduce growth. base_gamma = 0.20 reg_gamma, reg_mu = regularized_parameter_fit() report["mini_optimization"] = { "task_target_gamma": base_gamma, "baseline_gamma": base_gamma, "baseline_mu": float(hill_exponents(gamma=base_gamma, filtered=True)[0][0]), "regularized_gamma": float(reg_gamma), "regularized_mu": float(reg_mu), "baseline_rho": float(np.exp(T * hill_exponents(gamma=base_gamma, filtered=True)[0][0])), "regularized_rho": float(np.exp(T * reg_mu)), } Path("results.json").write_text(json.dumps(report, indent=2)) print(json.dumps(report, indent=2)) if __name__ == "__main__": main()