1import json
 2import numpy as np
 3
 4# Harmonic oscillator H=(q^2+p^2)/2, dq=p, dp=-q.
 5def leapfrog(q, p, h):
 6    p = p - 0.5 * h * q
 7    q = q + h * p
 8    p = p - 0.5 * h * q
 9    return q, p
10
11def map_jacobian(h):
12    # Exact linear map from differentiating kick-drift-kick.
13    return np.array([
14        [1.0 - h*h/2.0, h],
15        [-h + h**3/4.0, 1.0 - h*h/2.0],
16    ])
17
18def run(h, total_time=200.0):
19    q, p = 1.0, 0.0
20    values = []
21    for _ in range(round(total_time / h)):
22        values.append(0.5 * (q*q + p*p))
23        q, p = leapfrog(q, p, h)
24    values = np.asarray(values)
25    return float(values.max() - values.min()), float(abs(values[-1] - values[0]))
26
27def main():
28    hs = np.array([0.4, 0.2, 0.1, 0.05, 0.025])
29    ranges = np.array([run(float(h))[0] for h in hs])
30    slope = float(np.polyfit(np.log(hs), np.log(ranges), 1)[0])
31    determinants = [float(np.linalg.det(map_jacobian(float(h)))) for h in hs]
32    output = {
33        'step_sizes': hs.tolist(),
34        'energy_ranges': ranges.tolist(),
35        'loglog_energy_range_slope': slope,
36        'jacobian_determinants': determinants,
37        'symplectic_volume_preservation': bool(max(abs(d - 1.0) for d in determinants) < 1e-14),
38        'second_order_energy_scaling_observed': bool(1.7 < slope < 2.3),
39    }
40    print(json.dumps(output, indent=2))
41    with open('verification.json', 'w') as f:
42        json.dump(output, f, indent=2)
43
44if __name__ == '__main__':
45    main()