1import json, math
 2import numpy as np
 3
 4M=12
 5D=0.73
 6L=np.diag([1]+[2]*(M-2)+[1])-np.diag(np.ones(M-1),1)-np.diag(np.ones(M-1),-1)
 7evals,evecs=np.linalg.eigh(L)
 8lmax=float(evals[-1])
 9alpha_c=2/(D*lmax)
10# First nonconstant mode: exact predicted geometric decay.
11v=evecs[:,1]
12alpha=0.5/(D*lmax)
13q=v.copy(); amps=[]
14for _ in range(100):
15    amps.append(abs(q@v)); q=q-alpha*D*(L@q)
16slope=float(np.polyfit(np.arange(1,50),np.log(amps[1:50]),1)[0])
17pred=math.log(abs(1-alpha*D*evals[1]))
18# Slightly above the claimed boundary: exact highest mode grows.
19au=1.01*alpha_c
20q=evecs[:,-1].copy(); norms=[]
21for _ in range(20):
22    norms.append(float(np.linalg.norm(q)))
23    q=q-au*D*(L@q)
24result={
25  'lambda_max':lmax,
26  'alpha_critical_predicted':alpha_c,
27  'decay_slope_observed':slope,
28  'decay_slope_predicted':pred,
29  'decay_ratio':slope/pred,
30  'above_boundary_alpha':au,
31  'highest_mode_growth_observed':norms[-1]/norms[0],
32  'highest_mode_growth_predicted':abs(1-au*D*lmax)**19
33}
34with open('math_verification.json','w') as f: json.dump(result,f,indent=2)
35print(json.dumps(result,indent=2))