1import json
 2import numpy as np
 3
 4R = np.array([[0.0, 1.0], [-1.0, 0.0]])
 5GAIN = 1.04
 6DAMP = 0.8
 7
 8def g(x):
 9    x = np.asarray(x, float)
10    q = 1.0 + DAMP * (x @ x)
11    return GAIN * (R @ x) / q
12
13def jac(x):
14    x = np.asarray(x, float)
15    q = 1.0 + DAMP * (x @ x)
16    z = R @ x
17    # derivative of GAIN*(R x)/(1+DAMP*x^T x)
18    return GAIN * (R / q - np.outer(z, 2.0 * DAMP * x) / (q*q))
19
20def nonlinear(x0, H):
21    xs = [np.asarray(x0, float)]
22    for _ in range(H):
23        xs.append(g(xs[-1]))
24    return np.asarray(xs)
25
26def frozen(x0, H, radius=None):
27    A = jac(x0)
28    c = g(x0) - A @ x0
29    rawrho = max(abs(np.linalg.eigvals(A)))
30    if radius is not None and rawrho > radius:
31        A *= radius / rawrho
32        c = g(x0) - A @ x0
33    rho = max(abs(np.linalg.eigvals(A)))
34    xs = [np.asarray(x0, float)]
35    for _ in range(H):
36        xs.append(A @ xs[-1] + c)
37    return np.asarray(xs), rawrho, rho
38
39def main():
40    H = 100
41    x0 = np.array([0.1, 0.0])
42    truth = nonlinear(x0, H)
43    raw, rawrho, _ = frozen(x0, H, None)
44    clipped, _, cliprho = frozen(x0, H, 0.98)
45    d = np.array([0.7, -0.3]); h = 1e-6
46    fd = (g(x0 + h*d) - g(x0))/h
47    out = {
48      'initial_jacobian_spectral_radius': float(rawrho),
49      'clipped_spectral_radius': float(cliprho),
50      'horizon': H,
51      'final_norms': {'nonlinear': float(np.linalg.norm(truth[-1])),
52                      'frozen_unclipped': float(np.linalg.norm(raw[-1])),
53                      'frozen_clipped': float(np.linalg.norm(clipped[-1]))},
54      'max_norms': {'nonlinear': float(np.max(np.linalg.norm(truth, axis=1))),
55                    'frozen_unclipped': float(np.max(np.linalg.norm(raw, axis=1))),
56                    'frozen_clipped': float(np.max(np.linalg.norm(clipped, axis=1)))},
57      'finite_difference_jacobian_error': float(np.linalg.norm(fd - jac(x0) @ d))
58    }
59    open('unstable_results.json','w').write(json.dumps(out, indent=2))
60    print(json.dumps(out, indent=2))
61
62if __name__ == '__main__':
63    main()