import json, math
import numpy as np

M=12
D=0.73
L=np.diag([1]+[2]*(M-2)+[1])-np.diag(np.ones(M-1),1)-np.diag(np.ones(M-1),-1)
evals,evecs=np.linalg.eigh(L)
lmax=float(evals[-1])
alpha_c=2/(D*lmax)
# First nonconstant mode: exact predicted geometric decay.
v=evecs[:,1]
alpha=0.5/(D*lmax)
q=v.copy(); amps=[]
for _ in range(100):
    amps.append(abs(q@v)); q=q-alpha*D*(L@q)
slope=float(np.polyfit(np.arange(1,50),np.log(amps[1:50]),1)[0])
pred=math.log(abs(1-alpha*D*evals[1]))
# Slightly above the claimed boundary: exact highest mode grows.
au=1.01*alpha_c
q=evecs[:,-1].copy(); norms=[]
for _ in range(20):
    norms.append(float(np.linalg.norm(q)))
    q=q-au*D*(L@q)
result={
  'lambda_max':lmax,
  'alpha_critical_predicted':alpha_c,
  'decay_slope_observed':slope,
  'decay_slope_predicted':pred,
  'decay_ratio':slope/pred,
  'above_boundary_alpha':au,
  'highest_mode_growth_observed':norms[-1]/norms[0],
  'highest_mode_growth_predicted':abs(1-au*D*lmax)**19
}
with open('math_verification.json','w') as f: json.dump(result,f,indent=2)
print(json.dumps(result,indent=2))