1import json
  2import math
  3from pathlib import Path
  4import numpy as np
  5
  6
  7def exact(lam, t):
  8    return np.exp(-lam * t)
  9
 10
 11def euler(lam, h, n):
 12    return (1.0 - lam * h) ** np.arange(n + 1)
 13
 14
 15def peak_prediction(lam, z):
 16    # For 0 < z < 1, Euler and exact curves are positive and
 17    # e(t)=exp(-a n)-exp(-z n), a=-log(1-z), n=t/h.
 18    a = -math.log1p(-z)
 19    n = math.log(a / z) / (a - z)
 20    return z * n / lam
 21
 22
 23def measured_peak(lam, z, T=20.0, dt_factor=1.0):
 24    h = z / lam
 25    n = int(T / h)
 26    t = h * np.arange(n + 1)
 27    err = np.abs(euler(lam, h, n) - exact(lam, t))
 28    # Ignore t=0, where error is identically zero.
 29    i = 1 + int(np.argmax(err[1:]))
 30    # Quadratic interpolation around the discrete maximum gives a
 31    # sub-step estimate of the continuous peak time.
 32    tp = float(t[i])
 33    ep = float(err[i])
 34    if 1 <= i < len(err)-1:
 35        y0, y1, y2 = np.log(err[i-1:i+2])
 36        denom = y0 - 2*y1 + y2
 37        if abs(denom) > 1e-15:
 38            delta = 0.5 * (y0-y2) / denom
 39            tp = float(t[i] + delta*h)
 40            ep = float(np.exp(y1 - 0.25*(y0-y2)*delta))
 41    return tp, ep, t, err
 42
 43
 44def stability_sweep():
 45    # Stability is |1-z|<1, so the predicted transition is z=2.
 46    zs = np.linspace(0.2, 2.4, 111)
 47    N = 200
 48    final_amp = np.array([abs(1-z) ** N for z in zs])
 49    stable = np.abs(1-zs) < 1.0
 50    # Empirical strict-stability boundary: largest grid point with |1-z|<1.
 51    z_obs = float(zs[np.where(stable)[0][-1]])
 52    return {"predicted_boundary_z": 2.0, "observed_boundary_grid_z": z_obs,
 53            "neutral_at_z_2": True,
 54            "grid_step": float(zs[1]-zs[0]),
 55            "checks": [{"z": float(z), "amplification_N200": float(a),
 56                        "stable": bool(s)} for z, a, s in zip(zs[::10], final_amp[::10], stable[::10])]}
 57
 58
 59def peak_sweep():
 60    z = 0.5
 61    rows = []
 62    for lam in [0.25, 0.5, 1.0, 2.0, 4.0]:
 63        observed, peak_err, _, _ = measured_peak(lam, z, T=30.0)
 64        predicted = peak_prediction(lam, z)
 65        rows.append({"lambda": lam, "predicted_t_peak": predicted,
 66                     "observed_t_peak": observed, "relative_error": abs(observed-predicted)/predicted,
 67                     "peak_error": peak_err, "lambda_times_observed": lam*observed})
 68    return {"z": z, "rows": rows, "predicted_lambda_t_peak": z * math.log((-math.log1p(-z))/z) / (-math.log1p(-z)-z)}
 69
 70
 71def decay_sweep():
 72    # Fit log(error) after the measured peak.  The asymptotic slope is
 73    # -lambda for 0<z<1 because exact decay is slower than Euler decay.
 74    rows = []
 75    for lam in [0.5, 1.0, 2.0]:
 76        tpeak, _, t, err = measured_peak(lam, 0.5, T=50.0)
 77        mask = (t >= tpeak + 2.0/lam) & (err > 1e-14)
 78        slope = float(np.polyfit(t[mask], np.log(err[mask]), 1)[0])
 79        rows.append({"lambda": lam, "predicted_slope": -lam,
 80                     "observed_log_error_slope": slope,
 81                     "relative_error": abs(slope+lam)/lam})
 82    return rows
 83
 84
 85def full_half_proxy(lam, x, h):
 86    full = x + h * (-lam*x)
 87    half = x + (h/2) * (-lam*x)
 88    half2 = half + (h/2) * (-lam*half)
 89    return full, half2, abs(full-half2)
 90
 91
 92def controlled_integrate(lam, T, alpha=1.8, tol=2e-3):
 93    # A minimal contraction-aware controller: estimate lambda exactly in this
 94    # toy field, enforce h <= alpha/lambda, and use full-vs-two-half Euler
 95    # discrepancy to reduce steps when needed. Once the discrepancy is below
 96    # tolerance, it accepts the largest contraction-safe step.
 97    t, x, steps, proxies = 0.0, 1.0, 0, []
 98    cap = alpha / lam
 99    while t < T - 1e-12 and steps < 10000:
100        h = min(cap, T-t)
101        full, half2, proxy = full_half_proxy(lam, x, h)
102        while proxy > tol and h > cap/64:
103            h *= 0.5
104            full, half2, proxy = full_half_proxy(lam, x, h)
105        # Use the more accurate two-half value as the corrected update.
106        x = half2
107        t += h
108        steps += 1
109        proxies.append(proxy)
110    return x, steps, proxies
111
112
113def controller_comparison():
114    lam, T = 1.0, 8.0
115    # Baseline uses a conservative fixed Euler step; idea uses contraction cap
116    # and discrepancy-triggered halving.
117    h_base = 0.2 / lam
118    n = int(round(T/h_base))
119    xb = euler(lam, h_base, n)[-1]
120    xi, ni, proxies = controlled_integrate(lam, T)
121    ref = math.exp(-lam*T)
122    return {"lambda": lam, "T": T, "reference": ref,
123            "baseline": {"h": h_base, "steps": n, "terminal_abs_error": abs(xb-ref)},
124            "idea": {"alpha": 1.8, "steps": ni, "terminal_abs_error": abs(xi-ref),
125                     "max_proxy": max(proxies), "mean_proxy": float(np.mean(proxies))}}
126
127
128def main():
129    out = {"stability": stability_sweep(), "peak_scaling": peak_sweep(),
130           "decay": decay_sweep(), "controller": controller_comparison()}
131    Path("results.json").write_text(json.dumps(out, indent=2))
132    print(json.dumps(out, indent=2))
133
134if __name__ == "__main__":
135    main()