1"""The struck FPUT-beta chain, vectorized over an ensemble of strikes.
3N unit masses between two fixed walls (u_0 = u_{N+1} = 0), joined by springs with force law
4 f(Delta) = Delta + beta Delta^3 (Fermi-Pasta-Ulam-Tsingou, beta model)
5and optionally held by on-site springs of stiffness mu^2 (which makes the spectrum inharmonic). The
6energy is
7 H = sum_j v_j^2 / 2 + sum_{j=0}^{N} (Delta_j^2 / 2 + beta Delta_j^4 / 4) + mu^2 sum_j u_j^2 / 2,
8with Delta_j = u_{j+1} - u_j. The linear normal modes are sin(j k pi / (N + 1)), k = 1..N, with
9 omega_k^2 = mu^2 + 4 sin^2(k pi / (2 (N + 1))).
11Damping: friction c0 on every mass and dashpots c1 between neighbours,
12 dv_j/dt += -c0 v_j + c1 (v_{j+1} - 2 v_j + v_{j-1}),
13so that mode k decays (in amplitude) at (c0 + c1 (omega_k^2 - mu^2)) / 2. The dashpot term is the chain's version
14of the plate's sigma1 gradient damping.
16Integration: 4th-order symplectic composition (Yoshida) of leapfrog for the conservative part; after each
17step the velocities are multiplied by the exact propagator exp(-dt D) of the damping (an N x N matrix,
18precomputed). The tangent dynamics (largest Lyapunov exponent) uses the same splitting, linearized.
19"""
21from __future__ import annotations
23from dataclasses import dataclass
25import numpy as np
26from scipy.linalg import expm
28_C = 2 ** (1 / 3)
29YOSHIDA = (1 / (2 - _C), -_C / (2 - _C), 1 / (2 - _C))
33class Chain:
34 N: int = 32
35 beta: float = 1.0
36 mu: float = 0.0
37 c0: float = 0.0
38 c1: float = 0.0
41 def k(self):
42 return np.arange(1, self.N + 1)
45 def omega(self):
46 return np.sqrt(self.mu**2 + 4 * np.sin(self.k * np.pi / (2 * (self.N + 1))) ** 2)
49 def modes(self):
50 """Orthonormal mode matrix S[j, k] = sqrt(2 / (N + 1)) sin(j k pi / (N + 1)), j, k = 1..N."""
51 j = np.arange(1, self.N + 1)
52 return np.sqrt(2 / (self.N + 1)) * np.sin(np.outer(j, self.k) * np.pi / (self.N + 1))
54 def damping_matrix(self):
55 N = self.N
56 L = -2 * np.eye(N) + np.eye(N, k=1) + np.eye(N, k=-1)
57 return self.c0 * np.eye(N) - self.c1 * L
59 def force(self, u):
60 """u: (M, N) -> force (M, N)."""
61 up = np.pad(u, ((0, 0), (1, 1)))
62 d = np.diff(up, axis=1) # (M, N + 1) stretches
63 f = d + self.beta * d**3
64 return f[:, 1:] - f[:, :-1] - self.mu**2 * u
66 def tangent_force(self, u, du):
67 up = np.pad(u, ((0, 0), (1, 1)))
68 dup = np.pad(du, ((0, 0), (1, 1)))
69 d = np.diff(up, axis=1)
70 dd = np.diff(dup, axis=1)
71 g = (1 + 3 * self.beta * d**2) * dd
72 return g[:, 1:] - g[:, :-1] - self.mu**2 * du
74 def energy(self, u, v):
75 up = np.pad(u, ((0, 0), (1, 1)))
76 d = np.diff(up, axis=1)
77 return 0.5 * (v**2).sum(1) + (0.5 * d**2 + 0.25 * self.beta * d**4).sum(1) + 0.5 * self.mu**2 * (u**2).sum(1)
79 def mode_energies(self, u, v):
80 S = self.modes
81 q, p = u @ S, v @ S
82 return 0.5 * (p**2 + (self.omega * q) ** 2)
85def strike(chain: Chain, E, center, width, M=1, scale=None):
86 """Initial state of a strike: at rest at u = 0, velocity profile exp(-((j - center) / width)^2 / 2)
87 normalized to kinetic energy E. scale (M,) multiplies the velocity of each member (for ensembles)."""
88 j = np.arange(1, chain.N + 1)
89 prof = np.exp(-0.5 * ((j - center) / width) ** 2)
90 prof *= np.sqrt(2 * E / (prof**2).sum())
91 v = np.tile(prof, (M, 1))
92 if scale is not None:
93 v *= np.asarray(scale)[:, None]
94 return np.zeros_like(v), v
97def run(chain: Chain, u0, v0, t_end, dt, lyap=False, record_every=None, renorm_every=10):
98 """Integrate an ensemble (M, N). Returns final state, the finite-time Lyapunov exponent per member
99 (if lyap), and records (t, mode energies, total energy, Lyapunov log-sum) every record_every steps."""
100 u, v = np.array(u0, float), np.array(v0, float)
101 M = u.shape[0]
102 D = chain.damping_matrix()
103 damp = expm(-dt * D) if (chain.c0 or chain.c1) else None
104 n = int(round(t_end / dt))
105 if lyap:
106 rng = np.random.default_rng(0)
107 du = rng.standard_normal(u.shape)
108 dv = rng.standard_normal(v.shape)
109 nrm = np.sqrt((du**2).sum(1) + (dv**2).sum(1))[:, None]
110 du /= nrm
111 dv /= nrm
112 logsum = np.zeros(M)
113 rec = dict(t=[], Ek=[], E=[], logsum=[])
114 for i in range(n):
115 for c in YOSHIDA:
116 h = c * dt
117 v = v + 0.5 * h * chain.force(u)
118 if lyap:
119 dv = dv + 0.5 * h * chain.tangent_force(u, du)
120 du = du + h * dv
121 u = u + h * v
122 v = v + 0.5 * h * chain.force(u)
123 if lyap:
124 dv = dv + 0.5 * h * chain.tangent_force(u, du)
125 if damp is not None:
126 v = v @ damp.T
127 if lyap:
128 dv = dv @ damp.T
129 if lyap and (i + 1) % renorm_every == 0:
130 nrm = np.sqrt((du**2).sum(1) + (dv**2).sum(1))
131 logsum += np.log(nrm)
132 du /= nrm[:, None]
133 dv /= nrm[:, None]
134 if record_every and (i + 1) % record_every == 0:
135 rec["t"].append((i + 1) * dt)
136 rec["Ek"].append(chain.mode_energies(u, v))
137 rec["E"].append(chain.energy(u, v))
138 rec["logsum"].append(logsum.copy())
139 out = dict(u=u, v=v)
140 if lyap:
141 out["lyap"] = logsum / (n * dt)
142 if record_every:
143 out.update({k: np.array(x) for k, x in rec.items()})
144 return out
147def spectral_entropy(Ek):
148 """Normalized effective number of modes exp(S) / N, S = -sum e ln e, e = E_k / sum E_k (last axis).
149 1 at equipartition, 1/N when all energy is in one mode."""
150 e = Ek / Ek.sum(-1, keepdims=True)
151 S = -(e * np.log(np.where(e > 0, e, 1))).sum(-1)
152 return np.exp(S) / Ek.shape[-1]