#!/usr/bin/env python3
"""Recompute A02 examples; analytic bounds are proved in the linked pages.
Requires mpmath; 70 decimal digits. Sampling checks formulas, not a supremum proof.
Run: python foundations-singular-perturbation-check.py [output.json]
"""
import json,sys
from pathlib import Path
import mpmath as mp
mp.mp.dps=70
checks=[]
def check(label,ok):
 if not ok:raise AssertionError(label)
 checks.append(label)
def st(x):return mp.nstr(x,28)
def grid(a,b,n=400):return [a+(b-a)*k/n for k in range(n+1)]
# Regular expansion: exact geometric remainder, nonnegative.
regular=[]
for eps in map(mp.mpf,['0.1','0.25','0.5']):
 for t in grid(mp.mpf(0),mp.mpf(3),80):
  v=mp.exp(-t)*(1+eps*(1-mp.exp(-t)))
  y=mp.exp(-t)/(1-eps*(1-mp.exp(-t)))
  rem=eps**2*mp.exp(-t)*(1-mp.exp(-t))**2/(1-eps*(1-mp.exp(-t)))
  check('regular exact remainder',abs(y-v-rem)<mp.mpf('1e-65'))
  check('regular bound',0<=rem<=2*eps**2)
 regular.append({'epsilon':st(eps),'proven_bound':st(2*eps**2)})
# Boundary layer and composite: exact sup norm = q, attained at right endpoint.
boundary=[]
for eps in map(mp.mpf,['0.2','0.1','0.05','0.025']):
 q=mp.exp(-1/eps)
 def exact(x):return x-(-mp.expm1(-x/eps))/(-mp.expm1(-1/eps))
 def comp(x):return x-1+mp.exp(-x/eps)
 check('boundary endpoints',abs(exact(0))<mp.mpf('1e-65') and abs(exact(1))<mp.mpf('1e-65'))
 for x in grid(mp.mpf(0),mp.mpf(1),200):
  err=exact(x)-comp(x); identity=-q*(-mp.expm1(-x/eps))/(1-q)
  check('boundary error identity',abs(err-identity)<mp.mpf('1e-65'))
  check('boundary global bound sampled',abs(err)<=q+mp.mpf('1e-65'))
 for x in map(mp.mpf,['0','0.13','0.61','1']):
  check('boundary equation residual',abs(eps*mp.diff(exact,x,2)+mp.diff(exact,x)-1)<mp.mpf('1e-60'))
 check('bound attained at right endpoint',abs(abs(exact(1)-comp(1))-q)<mp.mpf('1e-65'))
 boundary.append({'epsilon':st(eps),'exact_sup_error':st(q),'outer_error_at_left':'1'})
# Coupled slow-fast model, quantitative theorem, T fixed at 1.
slowfast=[];T=mp.mpf(1)
for eps in [mp.mpf(1)/n for n in [4,8,16,32]]:
 rp=2/(1+mp.sqrt(1+4*eps));rm=-(1+mp.sqrt(1+4*eps))/(2*eps);D=rp-rm
 def x(t):return (-rm*mp.exp(rp*t)+rp*mp.exp(rm*t))/D
 def z(t):return -rm*rp*(mp.exp(rp*t)-mp.exp(rm*t))/D
 K=mp.exp(T)*(1+T*mp.exp(T));threshold=eps*mp.log(1/eps)
 ex=[];ez=[];ez_after=[]
 for t in grid(mp.mpf(0),T):
  xe=abs(x(t)-mp.exp(t));ze=abs(z(t)-mp.exp(t));track=abs(z(t)-x(t))
  check('slow-fast invariant region',-mp.mpf('1e-65')<=z(t)<=x(t)+mp.mpf('1e-65') and x(t)<=mp.exp(t)+mp.mpf('1e-65'))
  check('slow whole-window certificate',xe<=eps*K)
  check('fast moving-target certificate',track<=mp.exp(-t/eps)+eps*mp.exp(T))
  check('fast reduced-target certificate',ze<=mp.exp(-t/eps)+eps*(mp.exp(T)+K))
  ex.append(xe);ez.append(ze)
 for t in grid(threshold,T):
  ze=abs(z(t)-mp.exp(t));ez_after.append(ze)
  check('fast after-layer certificate',ze<=eps*(1+mp.exp(T)+K))
 for t in map(mp.mpf,['0','0.17','0.59','1']):
  check('slow differential residual',abs(mp.diff(x,t)-z(t))<mp.mpf('1e-60'))
  check('fast differential residual',abs(eps*mp.diff(z,t)-x(t)+z(t))<mp.mpf('1e-60'))
 check('initial fast error is not small',abs(z(0)-1)==1)
 slowfast.append({'epsilon':st(eps),'T':st(T),'layer_cutoff':st(threshold),'sampled_max_slow_error':st(max(ex)),'slow_proven_bound':st(eps*K),'sampled_max_fast_error_all':st(max(ez)),'sampled_max_fast_error_after_layer':st(max(ez_after)),'fast_after_layer_proven_bound':st(eps*(1+mp.exp(T)+K))})
# Weak damping: full long-window phase-error certificate.
oscillator=[];Ts=mp.mpf(2)
for eps in map(mp.mpf,['0.2','0.1','0.05','0.025']):
 omega=mp.sqrt(1-eps**2)
 def exact(t):return mp.exp(-eps*t)*mp.cos(omega*t)
 def approx(t):return mp.exp(-eps*t)*mp.cos(t)
 errors=[];bound=Ts*eps/(1+omega)
 for t in grid(mp.mpf(0),Ts/eps,1600):
  err=abs(exact(t)-approx(t));errors.append(err)
  check('oscillator pointwise phase bound',err<=mp.exp(-eps*t)*eps**2*t/(1+omega)+mp.mpf('1e-65'))
  check('oscillator long-window bound',err<=bound)
 check('oscillator initial value',exact(0)==1)
 check('oscillator initial derivative',abs(mp.diff(exact,0)+eps)<mp.mpf('1e-65'))
 for t in map(mp.mpf,['0','0.3','4']):check('oscillator equation residual',abs(mp.diff(exact,t,2)+2*eps*mp.diff(exact,t)+exact(t))<mp.mpf('1e-60'))
 oscillator.append({'epsilon':st(eps),'physical_horizon':st(Ts/eps),'sampled_max_error':st(max(errors)),'proven_bound':st(bound)})
secular=[]
for n in [5,20,80]:
 eps=1/(2*mp.pi*n);t=2*mp.pi*n;truth=mp.exp(-1)*mp.cos(mp.sqrt(1-eps**2)*t)
 secular.append({'n':n,'epsilon':st(eps),'time':st(t),'ordinary_first_order':'0','truth':st(truth),'multiple_scales':st(mp.exp(-1))})
# Periodic averaging and a period-depending-on-epsilon counterexample.
averaging=[]
for eps in map(mp.mpf,['0.2','0.1','0.05']):
 errs=[]
 for t in grid(mp.mpf(0),Ts/eps,800):
  x=mp.exp(eps*t+eps*mp.sin(t));y=mp.exp(eps*t);err=abs(x-y);errs.append(err)
  check('averaging relative bound',abs(x/y-1)<=mp.expm1(eps)+mp.mpf('1e-65'))
  check('averaging absolute long-window bound',err<=mp.exp(Ts)*mp.expm1(eps)+mp.mpf('1e-65'))
 averaging.append({'epsilon':st(eps),'sampled_max_error':st(max(errs)),'proven_bound':st(mp.exp(Ts)*mp.expm1(eps))})
check('changing-period counterexample',abs(mp.sin(mp.pi/2)-1)<mp.mpf('1e-65'))
out={'status':'passed','precision_decimal_digits':70,'assertions':len(checks),'scope':'Exact identities and sampled verification of analytic bounds; grid maxima are NOT certified suprema. Boundary sup error has independent exact proof.','regular':regular,'boundary':boundary,'slow_fast':slowfast,'oscillator':oscillator,'secular_failure':secular,'averaging':averaging,'negative_cases':{'initial_fast_error':1,'varying_period_error_at_pi_over_2epsilon':1}}
path=Path(sys.argv[1]) if len(sys.argv)>1 else Path(__file__).with_name('foundations-singular-perturbation-results.json');path.write_text(json.dumps(out,ensure_ascii=False,indent=2)+'\n');print(json.dumps(out,ensure_ascii=False,indent=2))