1# Cymbal2
3A time-domain physical simulation of a struck cymbal, after Samejima (2021), "Nonlinear physical modeling sound synthesis of cymbals involving dynamics of washers and sticks/mallets" (Acoust. Sci. & Tech. 42, 6). The body is a shallow spherical shell governed by the von Karman equations, so the crash is not modeled but computed: the geometric nonlinearity that couples bending to mid-surface stretching is integrated directly, and the slow migration of energy into the high frequencies after a hard stroke, the pitch glides, and the noise-like bifurcations of a crash emerge from the PDE. The support is not a rigid clamp but a washer of finite impedance, entering as the shear boundary condition at the mounting circle, which is what lets the whole shell bounce and tilt on its felt the way a mounted cymbal visibly does. The striking mechanism is carried in full: a mallet head with a Hertz-type nonlinear contact, on a stick modeled as an elastic beam, held by a hand that grips either rigidly or through a rotational spring.
5The sibling [cymbal](../cymbal) repository models the same instrument by an entirely different route: a modal eigensolve of a profiled shell of revolution, with the nonlinearity carried statistically as a cascade of modal energies and the sound radiated to a listener point. The two models make opposite trade-offs. The modal cymbal renders in seconds, carries the hammered irregularity and a radiation model, but replaces the true nonlinear dynamics with an energy random walk that discards phase. Cymbal2 integrates the true dynamics, so intermodulation, pitch glide, and the attack's nonlinear character are real, the trade-off being that a render takes minutes, the plate is ideally axisymmetric, and the output is the shell's motion at a point rather than a pressure at a listener.
7## Contents
9| Path | What it is |
10|---|---|
11| `matlab/` | The algorithm in MATLAB syntax, run by [numbl](https://numbl.org) |
12| `webapp/` | Interactive browser app that executes the MATLAB code client-side |
14The MATLAB code runs from the command line with `cd matlab && numbl run cymbal2_demo.m`, or `cymbal2_validate.m` to run the checks. The web app is a Vite project; see `webapp/README.md`.
16## The model
18The signal chain, with each stage in its own `.m` file:
20```
21polar grid, azimuthal transform cymbal2_grid.m
22 -> radial operators per azimuthal wavenumber cymbal2_perm.m
23 -> implicit update + boundary rows, inverted cymbal2_system.m
24stick beam FEM, grip, head mass cymbal2_mallet.m
25 -> hand stroke -> Hertz contact -> von Karman
26 shell integration, step by step cymbal2_sim.m
27 -> pickup acceleration -> decimation -> PCM cymbal2_output.m
28```
30Five points are worth stating plainly, because they are where the model commits itself.
32**The shell is a uniform shallow spherical cap.** A real cymbal has a raised bell and a tapered bow whose curvatures differ; the paper (following Bilbao 2010) replaces both by one spherical cap whose rise sets a single radius of curvature, here 66 cm from a 4 cm rise over a 23 cm radius, with uniform 1 mm thickness. That is the price of the von Karman formulation, whose stress-function form is what makes the nonlinear terms computable: the displacement u and an Airy stress function couple through the bilinear operator L, the stretching enters as one Laplacian term in 1/Rs, and both elliptic solves separate over azimuthal wavenumber. The consequence is audible and known: the linear spectrum is that of an idealized cap, not of a hammered instrument, and there is no lathing or hammering irregularity anywhere in the model. The sibling cymbal repository is the converse commitment.
34**The support is an impedance, not a clamp.** The inner circle at the washer radius carries a moment-free condition and a shear balance against a mass-damper-spring per unit length of the mounting circumference (the paper's Eq. 15), with the measured washer stiffness 1.65e5 N/m per meter. Solved to convergence, this places the bounce of the whole shell on its felt near 13 Hz and the tilt near 4 Hz, both of which this discretization reproduces against the closed-form rigid-body values in the validation. These wobble modes are subsonic, but they matter twice: they set the boundary condition the audible modes actually see (the paper's point, visible in its Fig. 8 as a shifted low-frequency response), and under a hard stroke the shell swings away from the stick by centimeters, which is what makes the contact of a follow-through stroke a shove rather than a strike. `support_clamped = 1` restores the rigid clamp of the earlier literature and is worth hearing once: the low end stiffens and the wobble disappears.
36**The nonlinearity is integrated explicitly, and that sets a hardness ceiling.** Following the paper's scheme (Eq. 20), the stiff linear terms are implicit with the free parameter g = 0.2498, unconditionally stable on their own, while L(phi, u) and the membrane coupling are evaluated at the current step. An explicit quadratic term is only conditionally stable, and the condition tightens with amplitude: strokes whose tip arrives above roughly 10 m/s drive peak displacements past forty plate thicknesses and the integration diverges within a few hundred milliseconds, at any rate we tried. This is also the regime where the shallow-shell expansion itself has stopped being valid, so the model refuses more or less where its physics runs out; the renderer detects the divergence and reports it rather than emitting garbage. The practical range of `strike_velocity` is about 0.25 to 7 m/s of hand speed, which spans a fingertip tap to a violent crash. Two discretization choices govern the headroom, and not the ones intuition suggests. The azimuthal resolution matters more than the time step: the cascade carries energy up the nodal-diameter ladder, and if `n_theta` truncates that ladder short, energy piles up at the grid's cutoff and the integration blows up where the full grid rides it out, so a cheap render should drop `oversample` and never `n_theta`. The simulation runs at `oversample` times the output rate; 2 matches the paper's 2x rate and adds margin, but at the paper's grid even 1 carries the default stroke.
38**The striking mechanism is the paper's, with three repairs it needs in practice.** The stick is a 39 cm Euler-Bernoulli beam in cubic Hermite elements with consistent mass, damped in proportion to its mass, gripped at one end either rigidly or against a 140 N m/rad rotational spring, with the 28 g head at the tip and the Hertz-type contact f = K delta^2.54 against the shell. First, the paper's separate head SDOF coupled to the beam by a rigid constraint (its Eqs. 31-39) mixes an explicit and an implicit update through that constraint and the discrete pair is unstable; since the connection is rigid, the head mass is lumped onto the tip node, which is the same physics under substitution and inherits the beam scheme's stability. Second, the implicit weighting printed for the washer and beam schemes (the (1-phi)/2 split with phi = 0.7) is unstable whenever the mass term does not dominate, massless washer included; the weight belongs on the (k+1, k-1) average, and with it phi = 0.7 is unconditionally stable. Third, the paper's hand accelerates to a terminal speed and never reverses, so it cannot play a second stroke and eventually drags the stick through geometry no linear beam describes; here the hand accelerates at an arm-like 250 m/s^2, coasts, and on contact snaps back at 2500 m/s^2, returning to rest so a figure can repeat. Because the stick is compliant and the grip rotates, the tip arrives faster than the hand by a factor that settles around 1.3 when the clearance is sized to the stroke, which `strike_gap_m <= 0` (the default) does automatically. The multiple contacts of a real mallet bounce, which the paper shows as double force peaks, appear here the same way.
40**The output is the shell's motion at a point, not a microphone.** As in the paper, the rendered signal is the normal acceleration of the plate at a pickup point (default at 0.6 of the radius, opposite the strike), DC-blocked and decimated. This is the model's honest output: a radiation computation like the sibling cymbal's dipole sheet would need the full surface velocity history and is not part of the paper's method. The audible consequence is that low-order modes that radiate poorly from a real plate are present here at full strength, so the result is closer to a contact pickup than to a room microphone.
42Numerically, the solver departs from the paper in one structural way: instead of assembling the Kronecker-product operators over the full grid and solving one large system per step, both elliptic solves are block-diagonalized over azimuthal wavenumber (every operator and boundary row of the model is axisymmetric, so this is exact), and each step performs two batches of small precomputed radial solves plus dense transforms. The update it computes is the paper's, at about a hundredth the cost, which is the difference between a render taking minutes and taking hours. On grids coarse enough that a wavenumber's inner boundary layer cannot be resolved, the discretized boundary conditions can hand that wavenumber a negative stiffness; the system builder detects such columns by their spectral radius and truncates them, and at the default grid none are.
44## Validation
46`matlab/cymbal2_validate.m` (also `npm run selftest` in `webapp/`) checks the flat-plate limit against the closed-form free circular plate for one through four nodal diameters, the washer bounce and tilt against rigid-body closed forms in the stiff-plate limit, the clamped-center eigenfrequencies against the linear-FEM values published in the paper's Fig. 9 for the modes both discretizations resolve the same way, the exact invertibility of the azimuthal transform and the agreement of the physical-space derivative pipeline with the per-wavenumber operators, the non-amplification (spectral radius at most one) of every azimuthal column of the assembled update, the stick's static tip deflection under both grips against beam theory, exact reproducibility of a full render, pinned contact and output diagnostics from a reduced configuration, and that the crash decays. The nonlinear trajectory is sensitive to rounding, so the acoustic pins carry percent-level tolerances where the contact pins are tight.
48## Defaults
50The paper's measured instrument: an 18-inch medium crash, 23 cm radius, 4 cm rise, 1 mm bronze at 8860 kg/m^3, on its measured washer stiffness, struck by the paper's mallet (25 mm, 28 g head, on the measured stick) at 0.6 of the radius. The grid is the paper's, 24 radial intervals by 60 azimuthal points, and the damping constants are the paper's fitted sigma0 = 1.34 s^-1 and sigma1 = 1.2e-3, which put T60 near 5 s at the bottom of the spectrum and near 2 s at 5 kHz; the default tail of 6 s carries the ring to the end. The default stroke (5.75 m/s hand, the tip arriving at 8.2 m/s, 311 N peak contact) was chosen against the authors' published demonstration recording: it matches that recording's band balance to within about 2 dB from 100 Hz to 12 kHz, its high-frequency buildup over the first half second, and its decay envelope, the residual being a few dB of extra energy below 700 Hz, where the pickup hears the washer wobble. The same stroke's peak pickup acceleration, 2.4e4 m/s^2, sits at the scale of the paper's Fig. 7 waveforms. A default render integrates about 580,000 steps at 96 kHz and takes a few minutes in MATLAB or under the numbl CLI; `oversample = 1` halves that at some cost in headroom for hard strokes. In a browser with WebGPU the web app runs the same integration on the GPU in well under a minute (see `webapp/README.md`). The mesh resolves the dynamics reliably to roughly 2 kHz by the paper's own estimate; the crash extends far above that, but its top octaves are dispersion-warped and weak, which is the main respect in which this model is duller than a real crash. Raising `n_radial` and `n_theta` together sharpens it at quadratic cost.
52`strike_velocity` is the one control that matters most: 0.5 m/s is a tap that stays nearly linear and rings like the bell of a gong, 2 m/s begins to wash, and 4 m/s and above crash, with the high-frequency energy arriving over hundreds of milliseconds after the attack. `strike_count` with `strike_rate_hz` plays a repeated figure; each stroke's contact speed varies with `velocity_spread`. The washer parameters, the grip (`grip_clamped`, `grip_spring`), and the stick and head properties are all exposed and all audible in the attack.