1import json
  2from pathlib import Path
  3import numpy as np
  4from scipy.optimize import minimize
  5
  6
  7def mu_factors(N, m, q, weights=None):
  8    weights = np.ones(m) if weights is None else np.asarray(weights)
  9    ell = np.arange(1, m + 1)
 10    cq = np.cos(2 * np.pi * q * ell / N)
 11    return np.array([np.sum(weights * cq * (1 - np.cos(2 * np.pi * k * ell / N))) for k in range(N)])
 12
 13
 14def eigenvalues(N, m, q, kappa=1.0, weights=None):
 15    weights = np.ones(m) if weights is None else np.asarray(weights)
 16    ell = np.arange(1, m + 1)
 17    cq = np.cos(2 * np.pi * q * ell / N)
 18    return np.array([kappa * np.sum(weights * cq * (np.exp(2j * np.pi * k * ell / N) - 1)) for k in range(N)])
 19
 20
 21def rhs(theta, m, kappa=1.0, weights=None):
 22    weights = np.ones(m) if weights is None else np.asarray(weights)
 23    y = np.zeros_like(theta)
 24    for l, a in enumerate(weights, 1):
 25        y += a * np.sin(np.roll(theta, -l) - theta)
 26    return kappa * y
 27
 28
 29def mode_amplitude(x, profile, k):
 30    N = len(x)
 31    u = np.exp(2j * np.pi * k * np.arange(N) / N)
 32    return abs(np.vdot(u, x - profile) / N)
 33
 34
 35def simulate(profile, m, k, kappa, h, steps, weights=None, eps=1e-7):
 36    N = len(profile)
 37    u = np.exp(2j * np.pi * k * np.arange(N) / N)
 38    x = profile + eps * np.real(u)
 39    amps = [mode_amplitude(x, profile, k)]
 40    for _ in range(steps):
 41        x = x + h * rhs(x, m, kappa, weights)
 42        amps.append(mode_amplitude(x, profile, k))
 43    return np.asarray(amps)
 44
 45
 46def fit_slope(a, h):
 47    t = np.arange(len(a)) * h
 48    return float(np.polyfit(t, np.log(np.maximum(a, 1e-30)), 1)[0])
 49
 50
 51def main():
 52    N, m, q, k = 32, 4, 2, 1
 53    profile = 2 * np.pi * q * np.arange(N) / N
 54    mu = mu_factors(N, m, q)
 55    lam = eigenvalues(N, m, q)[k]
 56    # Prediction 1: infinitesimal slope equals -mu.
 57    a = simulate(profile, m, k, 1.0, .002, 1000)
 58    slope = fit_slope(a, .002)
 59    # Prediction 2: slope scales linearly with coupling.
 60    scaling = []
 61    for kap in [.25, .5, 1., 1.5]:
 62        aa = simulate(profile, m, k, kap, .001, 500)
 63        scaling.append({'kappa': kap, 'predicted': float(-kap * mu[k]), 'observed': fit_slope(aa, .001)})
 64    # Prediction 3: explicit-Euler boundary |1+h lambda|=1.
 65    lb = eigenvalues(N, m, 0)[5]
 66    hcrit = float(-2 * lb.real / abs(lb) ** 2)
 67    hs = np.linspace(.2 * hcrit, 1.8 * hcrit, 65)
 68    stable = np.abs(1 + hs * lb) <= 1
 69    crossings = np.where(stable[:-1] & ~stable[1:])[0]
 70    observed_hcrit = float((hs[crossings[0]] + hs[crossings[0] + 1]) / 2)
 71    # Negative-factor prediction.
 72    candidates = [(qq, kk, float(mu_factors(N, m, qq)[kk])) for qq in range(N) for kk in range(1, N)
 73                  if mu_factors(N, m, qq)[kk] < -.05]
 74    qq, kk, negmu = candidates[0]
 75    prof_neg = 2 * np.pi * qq * np.arange(N) / N
 76    an = simulate(prof_neg, m, kk, 1., .001, 1000)
 77    neg_slope = fit_slope(an, .001)
 78    # Baseline versus Fourier-shaped coupling on the same unstable profile/mode.
 79    baseline_w = np.ones(m)
 80    target = .25
 81    objective = lambda w: (mu_factors(N, m, qq, w)[kk] - target) ** 2 + .02 * np.sum((w - 1.) ** 2)
 82    opt = minimize(objective, baseline_w, method='L-BFGS-B', bounds=[(-3., 3.)] * m)
 83    shaped_w = opt.x
 84    base_mu = float(mu_factors(N, m, qq, baseline_w)[kk])
 85    shaped_mu = float(mu_factors(N, m, qq, shaped_w)[kk])
 86    base_amp = simulate(prof_neg, m, kk, 1., .001, 1000, baseline_w)
 87    shaped_amp = simulate(prof_neg, m, kk, 1., .001, 1000, shaped_w)
 88    comparison = {'profile_q': qq, 'mode': kk, 'target_mu': target,
 89      'baseline': {'weights': baseline_w.tolist(), 'mu': base_mu, 'slope': fit_slope(base_amp, .001),
 90                   'amplitude_ratio': float(base_amp[-1] / base_amp[0])},
 91      'fourier_shaped': {'weights': shaped_w.tolist(), 'mu': shaped_mu, 'slope': fit_slope(shaped_amp, .001),
 92                         'amplitude_ratio': float(shaped_amp[-1] / shaped_amp[0])}}
 93    report = {'config': {'N': N, 'm': m, 'q': q, 'mode': k},
 94      'growth_rate_test': {'predicted': float(-mu[k]), 'observed': slope, 'relative_error': abs(slope + mu[k]) / abs(mu[k])},
 95      'coupling_scaling_test': scaling,
 96      'euler_boundary_test': {'lambda': [float(lb.real), float(lb.imag)], 'predicted_hcrit': hcrit,
 97          'observed_hcrit': observed_hcrit, 'relative_error': abs(observed_hcrit - hcrit) / hcrit},
 98      'negative_factor_test': {'q': qq, 'mode': kk, 'mu': negmu, 'predicted_growth': -negmu,
 99          'observed_growth': neg_slope, 'relative_error': abs(neg_slope + negmu) / abs(negmu)},
100      'baseline_vs_shaped': comparison}
101    Path('results.json').write_text(json.dumps(report, indent=2))
102    print(json.dumps(report, indent=2))
103
104if __name__ == '__main__':
105    main()