1c822d5Add modal-series and fd-second-order submissions for string-pluck-v1, with results workflowJeremy Magland 1"""Finite-difference solution of string-pluck-v1, second order in space and time.
3With D the standard three-point second difference (zero at the fixed ends) and u^n the grid
4values at step n, the scheme is
6 (u^{n+1} - 2u^n + u^{n-1})/dt^2 = c^2 D u^n / dx^2
7 - sigma0 (u^{n+1} - u^{n-1}) / dt
8 + kappa D (u^{n+1} - u^{n-1}) / (2 dt dx^2).
10The Kelvin-Voigt term is centered in time, so the scheme is implicit in the damping and a
11tridiagonal system is solved at every step. (An earlier explicit variant with a backward
12difference for that term is only first order in time; its error cancelled against the
13second-order dispersion error near m = 8 and the sequence was not monotone.) Stability
14requires the Courant number c dt/dx <= 1.
16The time step is dt = 1/(fs*m) for an integer oversampling factor m, and the grid spacing is
17dx = c dt / lam with lam = 0.9, so the scheme is deliberately dispersive. The pickup velocity is
18(u^{n+1} - u^{n-1})/(2 dt) linearly interpolated to x_m, taken every m-th step.
20Usage: python fd_string.py <challenge.json> <out-dir> [m ...]
21"""
22from __future__ import annotations
24import json
25import sys
26import time
27from pathlib import Path
29import numpy as np
30from scipy.linalg import solve_banded
32sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
33from common import hardware_string, initial_shape, load_challenge, write_wav # noqa: E402
35LAM = 0.9
38def solve(p: dict, fs: int, num_samples: int, m: int, lam: float = LAM):
39 dt = 1.0 / (fs * m)
40 n_int = int(np.floor(p["L"] * lam / (p["c"] * dt))) # intervals; floor keeps dx >= c dt / lam
41 dx = p["L"] / n_int
42 x = np.linspace(0.0, p["L"], n_int + 1)
43 lam2 = (p["c"] * dt / dx) ** 2
44 s = p["sigma0"] * dt
45 beta = p["kappa"] * dt / (2.0 * dx**2)
47 def D(u):
48 out = np.zeros_like(u)
49 out[1:-1] = u[:-2] - 2.0 * u[1:-1] + u[2:]
50 return out
52 # Left-hand operator on interior points: (1 + s) I - beta D, tridiagonal, in solve_banded form.
53 n_in = n_int - 1
54 ab = np.zeros((3, n_in))
55 ab[0, 1:] = -beta
56 ab[1, :] = (1.0 + s) + 2.0 * beta
57 ab[2, :-1] = -beta
59 u_prev = initial_shape(x, p)
60 u_cur = u_prev + 0.5 * dt**2 * (p["c"] ** 2) * D(u_prev) / dx**2 # zero initial velocity, second-order start
61 xi = p["x_m"] / dx
62 i0 = int(np.floor(xi))
63 w1 = xi - i0
64 out = np.zeros(num_samples)
65 total_steps = num_samples * m
66 u_next = np.zeros_like(u_cur)
67 for n in range(1, total_steps + 1):
68 rhs = 2.0 * u_cur + lam2 * D(u_cur) - (1.0 - s) * u_prev - beta * D(u_prev)
69 u_next[1:-1] = solve_banded((1, 1), ab, rhs[1:-1], overwrite_b=True, check_finite=False)
70 if n % m == 0:
71 k = n // m
72 if k < num_samples:
73 v = (u_next - u_prev) / (2.0 * dt)
74 out[k] = (1.0 - w1) * v[i0] + w1 * v[i0 + 1]
75 u_prev, u_cur, u_next = u_cur, u_next, u_prev
76 return out, {"m": m, "dt": dt, "dx": dx, "n_points": n_int + 1, "courant": p["c"] * dt / dx}
79def main():
80 challenge = load_challenge(Path(sys.argv[1]))
81 out_dir = Path(sys.argv[2])
82 ms = [int(a) for a in sys.argv[3:]] or [1, 2, 4, 8, 16]
83 p = challenge["parameters"]
84 out = challenge["output"]
85 meta = []
86 for m in ms:
87 t0 = time.perf_counter()
88 sig, info = solve(p, out["sample_rate_hz"], out["num_samples"], m)
89 dt = time.perf_counter() - t0
90 sid = f"fd-m{m}"
91 write_wav(out_dir / f"{sid}.wav", sig, challenge)
92 meta.append({
93 "solution_id": sid,
94 "output": f"{out_dir.name}/{sid}.wav",
95 "method": "finite differences, second order, implicit damping",
96 "parameters": {"oversampling": m, "n_points": info["n_points"], "dx_m": info["dx"], "dt_s": info["dt"], "courant": round(info["courant"], 6)},
97 "timing_seconds": round(dt, 3),
98 "hardware": hardware_string(),
99 "peak_si": float(np.max(np.abs(sig))),
100 })
101 print(f"m={m} points={info['n_points']} peak={np.max(np.abs(sig)):.4f} m/s t={dt:.2f}s", file=sys.stderr)
102 print(json.dumps(meta, indent=1))
105if __name__ == "__main__":
106 main()