baba784Gong wire: FPUT chain with bending stiffness and a dome, experiments, web appJeremy Magland 1"""A struck beaded wire: the FPUT chain generalized toward a cymbal or gong.
3N unit masses on a wire between two pinned ends move sideways, with displacements y_j (y_0 = y_{N+1} = 0).
4The wire may be built as a shallow arch (a "dome") of height Hd: rest shape z_j = 4 Hd x (1 - x), x = j / (N + 1).
5Energy:
6 links: sum_j [ T D_j^2 / 2 + beta (z'_j D_j + D_j^2 / 2)^2 ], D_j = y_{j+1} - y_j, z'_j = z_{j+1} - z_j
7 bending: sum_j Kb (y_{j+1} - 2 y_j + y_{j-1})^2 / 2 (pinned ends: y_{-1} = -y_1, y_{N+2} = -y_N)
8 on-site: mu^2 sum_j y_j^2 / 2
9With Hd = 0 and Kb = 0 this is the FPUT-beta chain (sideways motion of a pre-tensioned beaded string; the link
10stretch is about D^2 / 2, so the stretching energy is quartic in D). The bending term makes the spectrum
11dispersive, omega ~ k^2 at large Kb, like a plate. The dome makes the link stretch linear in D as well
12(z' D), which adds stiffness (as curvature does in a shell) and a quadratic nonlinearity (z' D^3), the
13three-wave coupling of gongs and cymbals.
15The linear part is K = T L + 2 beta diag-ish(z'^2) + Kb B + mu^2 I; with a dome its normal modes are not sines, so
16they are computed numerically (eigh). Damping: friction c0 and neighbour dashpots c1 (Laplacian of the
17velocity), applied per step with Crank-Nicolson (tridiagonal). Time stepping: Stormer-Verlet (leapfrog).
18"""
20from __future__ import annotations
22from dataclasses import dataclass
24import numpy as np
25from scipy.linalg import eigh, solve_banded
29class Wire:
30 N: int = 64
31 T: float = 1.0
32 beta: float = 1.0
33 Kb: float = 0.0
34 Hd: float = 0.0
35 mu: float = 0.0
36 c0: float = 0.0
37 c1: float = 0.0
40 def zp(self):
41 x = np.arange(self.N + 2) / (self.N + 1)
42 z = 4 * self.Hd * x * (1 - x)
43 return np.diff(z) # (N + 1,) slopes of the links
45 def force(self, y):
46 """y: (M, N) -> force (M, N)."""
47 M, N = y.shape
48 yp = np.pad(y, ((0, 0), (1, 1)))
49 D = np.diff(yp, axis=1) # (M, N + 1)
50 zp = self.zp
51 e = zp * D + 0.5 * D * D # link stretch
52 fl = self.T * D + 2 * self.beta * e * (zp + D) # dV/dD
53 f = fl[:, 1:] - fl[:, :-1]
54 if self.Kb:
55 # bending: curvature c_j = y_{j+1} - 2 y_j + y_{j-1}, j = 0..N+1 with pinned ghosts
56 yg = np.concatenate([-y[:, :1], np.zeros((M, 1)), y, np.zeros((M, 1)), -y[:, -1:]], axis=1) # j = -1..N+2
57 c = yg[:, 2:] - 2 * yg[:, 1:-1] + yg[:, :-2] # j = 0..N+1
58 # dU/dy_j = Kb (c_{j-1} - 2 c_j + c_{j+1}), j = 1..N
59 f -= self.Kb * (c[:, :-2] - 2 * c[:, 1:-1] + c[:, 2:])
60 f -= self.mu**2 * y
61 return f
63 def stiffness(self):
64 """Linear stiffness matrix K (N x N) at y = 0."""
65 N = self.N
66 kl = self.T + 2 * self.beta * self.zp**2 # link stiffness at rest
67 K = np.zeros((N, N))
68 for j in range(N + 1): # link j joins y_j and y_{j+1} (1-based; 0 and N+1 fixed)
69 a, b = j - 1, j
70 if a >= 0:
71 K[a, a] += kl[j]
72 if b < N:
73 K[b, b] += kl[j]
74 if a >= 0 and b < N:
75 K[a, b] -= kl[j]
76 K[b, a] -= kl[j]
77 if self.Kb:
78 L = -2 * np.eye(N) + np.eye(N, k=1) + np.eye(N, k=-1) # Dirichlet Laplacian; pinned bending = Kb L^2
79 K += self.Kb * L @ L
80 K += self.mu**2 * np.eye(N)
81 return K
83 def modes(self):
84 w2, V = eigh(self.stiffness())
85 return np.sqrt(np.maximum(w2, 0)), V
87 def energy(self, y, v):
88 M, N = y.shape
89 yp = np.pad(y, ((0, 0), (1, 1)))
90 D = np.diff(yp, axis=1)
91 e = self.zp * D + 0.5 * D * D
92 U = (0.5 * self.T * D * D + self.beta * e * e).sum(1)
93 if self.Kb:
94 yg = np.concatenate([-y[:, :1], np.zeros((M, 1)), y, np.zeros((M, 1)), -y[:, -1:]], axis=1)
95 c = yg[:, 2:] - 2 * yg[:, 1:-1] + yg[:, :-2]
96 U += 0.5 * self.Kb * (c[:, 1:-1] ** 2).sum(1)
97 U += 0.5 * self.mu**2 * (y * y).sum(1)
98 return 0.5 * (v * v).sum(1) + U
101def strike(w: Wire, E, center, width, M=1, scale=None):
102 j = np.arange(1, w.N + 1)
103 prof = np.exp(-0.5 * ((j - center) / width) ** 2)
104 prof *= np.sqrt(2 * E / (prof**2).sum())
105 v = np.tile(prof, (M, 1))
106 if scale is not None:
107 v *= np.asarray(scale)[:, None]
108 return np.zeros_like(v), v
111def run(w: Wire, y0, v0, t_end, dt, record_every=100, lyap=False):
112 """Stormer-Verlet (kick-drift-kick), damping by Crank-Nicolson on v after each step."""
113 y, v = np.array(y0, float), np.array(v0, float)
114 M, N = y.shape
115 om, V = w.modes()
116 damp = bool(w.c0 or w.c1)
117 if damp:
118 # (I + dt/2 Dm) v+ = (I - dt/2 Dm) v, Dm = c0 I - c1 L
119 a = 0.5 * dt
120 ab = np.zeros((3, N))
121 ab[0, 1:] = -a * w.c1
122 ab[1, :] = 1 + a * (w.c0 + 2 * w.c1)
123 ab[2, :-1] = -a * w.c1
125 def apply_damp(x):
126 rhs = (1 - a * (w.c0 + 2 * w.c1)) * x
127 rhs[:, 1:] += a * w.c1 * x[:, :-1]
128 rhs[:, :-1] += a * w.c1 * x[:, 1:]
129 return solve_banded((1, 1), ab, rhs.T).T
130 if lyap:
131 rng = np.random.default_rng(0)
132 dy, dv = rng.standard_normal(y.shape), rng.standard_normal(v.shape)
133 nrm = np.sqrt((dy**2).sum(1) + (dv**2).sum(1))[:, None]
134 dy /= nrm
135 dv /= nrm
136 eps = 1e-7
137 logsum = np.zeros(M)
138 n = int(round(t_end / dt))
139 rec = dict(t=[], Ek=[], E=[], logsum=[])
140 f = w.force(y)
141 for i in range(n):
142 v = v + 0.5 * dt * f
143 y = y + dt * v
144 f = w.force(y)
145 v = v + 0.5 * dt * f
146 if lyap:
147 # tangent by finite differences of the force (keeps the model in one place)
148 g = (w.force(y + eps * dy) - f) / eps
149 dv = dv + dt * g
150 dy = dy + dt * dv
151 if damp:
152 v = apply_damp(v)
153 if lyap:
154 dv = apply_damp(dv)
155 if lyap and (i + 1) % 10 == 0:
156 nrm = np.sqrt((dy**2).sum(1) + (dv**2).sum(1))
157 logsum += np.log(nrm)
158 dy /= nrm[:, None]
159 dv /= nrm[:, None]
160 if (i + 1) % record_every == 0:
161 q, p = y @ V, v @ V
162 rec["t"].append((i + 1) * dt)
163 rec["Ek"].append(0.5 * (p * p + (om * q) ** 2))
164 rec["E"].append(w.energy(y, v))
165 rec["logsum"].append(logsum.copy())
166 out = {k: np.array(x) for k, x in rec.items()}
167 out.update(y=y, v=v, omega=om)
168 return out