# Power-balanced drift regulation for scalar auxiliary variable methods: application to real-time simulation of nonlinear string vibrations

Thomas Risse, Thomas Hélie, Stefan Bilbao

► To cite this version:

Thomas Risse, Thomas Hélie, Stefan Bilbao. Power-balanced drift regulation for scalar auxiliary variable methods: application to real-time simulation of nonlinear string vibrations. DAFx25, Sep 2025, Ancona, Italy. hal-05222856

HAL Id: hal-05222856

https://hal.science/hal-05222856v1

Submitted on 28 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# POWER-BALANCED DRIFT REGULATION FOR SCALAR AUXILIARY VARIABLE METHODS: APPLICATION TO REAL-TIME SIMULATION OF NONLINEAR STRING VIBRATIONS

Thomas Risse

S3AM Team
IRCAM, Sorbonne université
Paris, FR

thomas.risse@ircam.fr

Thomas Hélie

S3AM Team
IRCAM, CNRS
Paris, FR

thomas.helie@ircam.fr

Stefan Bilbao

Acoustics and Audio Group
University of Edinburgh
Edinburgh, UK

s.bilbao@ed.ac.uk

## ABSTRACT

Efficient stable integration methods for nonlinear systems are of great importance for physical modeling sound synthesis. Specifically, a number of musical systems of interest, including vibrating strings, bars or plates may be written as port-Hamiltonian systems with quadratic kinetic energy and non-quadratic potential energy. Efficient schemes have been developed for such systems through the introduction of a scalar auxiliary variable. As a result, the stable real-time simulations of nonlinear musical systems of up to a few thousands of degrees of freedom is possible, even for nearly lossless systems. However, convergence rates can be slow and seem to be system-dependent. Specifically, at audio rates, they may suffer from numerical drift of the auxiliary variable, resulting in dramatic unwanted effects on audio output, such as pitch drifts after several impacts on the same resonator.

In this paper, a novel method for mitigating this unwanted drift while preserving power balance is presented, based on a control approach. A set of modified equations is proposed to control the drift artefact by rerouting energy through the scalar auxiliary variable and potential energy state. Numerical experiments are run in order to check convergence on simulations in the case of a cubic nonlinear string. A real-time implementation is provided as a Max/MSP external. 60-note polyphony is achieved on a laptop, and some simple high level control parameters are provided, making the proposed implementation suitable for use in artistic contexts. All code is available in a public repository, along with compiled Max/MSP externals<fnref n="1" />.

## 1. INTRODUCTION

In the context of physical modelling sound synthesis, efficient numerical integration of differential equations is a great challenge. In order to achieve real-time performance, algorithms must be computationally cheap and numerically stable. A large variety of methods is available in the literature, including digital waveguides [1], modal synthesis [2] or finite difference methods [3]. One of the most difficult problems is in ensuring numerical stability for strongly nonlinear systems with an efficient numerical design.

Several numerical methods preserving an energy-like invariant [4, 5] are interesting to address this issue. They apply to the large

class of passive systems and allow to derive clear stability conditions. Many existing schemes, like discrete gradients methods [6] or projection methods [7, Section 5.3.3] rely on solving implicit nonlinear equations, for which iterative solvers such as Newton-Raphson are required. These are computationally prohibitive for real-time musical use, fundamentally serial, and lead to additional concerns in terms of convergence, stopping criteria, etc. For systems of cubic nonlinear type, including some models of string vibration, as well as plate vibration at high amplitudes, it is possible to arrive at schemes satisfying an energy balance, and with provable numerical stability conditions, See [3, chapter 8]. Such schemes are in general *linearly implicit*—and require the construction and solution of a linear system at each time step. This avoids the use of nonlinear iterative solvers but is still computationally costly.

In recent years, some effort has been directed toward establishing energy preserving schemes using energy quadratisation methods with auxiliary variables (see below) or without [8] (see also [9] for an alternative to quadratisation). The Invariant Energy Quadratisation (IEQ) method [10] is based on the simulation of the state space augmented by a distributed variable storing the square root of the energy density. Slightly later, the Scalar Auxiliary Variable (SAV) method [11] was introduced, in which the state space is augmented by a lumped variable storing the square root of the total energy of the system. The strength of these methods lies in the fact that the energy of the system may be expressed as a quadratic functional of the state variables, making any quadratic invariant preserving scheme a good candidate for stable simulations [12]. The musical instrument modelling community has shown great interest in these methods, for simulation of the geometrically nonlinear string [13]. Such methods remain, however, linearly implicit, and still reliant on linear system solutions in the time loop. More recently, using time-interleaved variables, it has been possible to arrive at *fully explicit* methods satisfying an energy balance for a general class of Hamiltonian systems [14]. In this case, matrix structure is exploited in order to avoid the linear system solution entirely (through the use of the Sherman-Morrison identity [15]), and such schemes can be viewed as stabilised forms of the classic explicit Störmer-Verlet method, with additional state required for storing the auxiliary variable. See [16, 17, 18] for musical applications, including in real time [19]. Other energy-conserving explicit formulations have been proposed, but without a numerical stability guarantee [9].

Even though SAV methods can be computationally efficient, they may suffer from unwanted drifts of the auxiliary variable, resulting in bad long term behaviour of the simulations: reference [20] theoretically discusses this issue as well as convergence prop-

<footnote n="1">¹https://github.com/thomas-risse/SAV-string-simulations</footnote>

<footnote>Copyright: © 2025 Thomas Risse et al. This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, adaptation, and reproduction in any medium, provided the original author and source are credited.</footnote>
erties of different schemes for the musical string, clearly indicating that this problem is not trivial. In [21], alternative schemes are proposed to address this issue in the context of collision problems by constraining the sign of the auxiliary variable. Outside of the musical acoustics community, this problem is also well known: reference [22] introduces a relaxation step correcting the value of the auxiliary variable by solving an optimization problem, for which an efficient solution is derived in [23]. While these modifications ensure stability, they introduce numerical dissipation, which is undesirable for simulating nearly lossless musical instruments.

In this paper, we present a novel method to address the auxiliary variable drift problem. Our approach is based on a modification of the continuous time equations and is suitable for port-Hamiltonian systems (see [24] for an introduction) with quadratic kinetic energy and general non-quadratic potential energy. This method is applied to a nonlinear string model, a relevant case for musical acoustics. The scheme is also implemented in a Max/MSP external which is publicly available¹, along with example sounds and an interactive browser-based simulation on the companion webpage.

The paper is organized as follows: section 2 presents the string model and its semi-discrete approximation using finite differences. Section 3 reviews the Scalar Auxiliary Variable transformation. As a main contribution, Section 4 proposes a modification of the dynamical system in the continuous time domain, based on a control approach, to reduce any SAV drift (due to e.g. inconsistent initial conditions or, numerical errors). An interleaved time stepping scheme is provided in section 5 and evaluated through numerical experiments in section 6. Section 7 discusses the real-time implementation and is followed by a general conclusion in section 8.

## 2. ENERGY CONSISTENT MODELLING

### 2.1. String equations (PDE)

Consider the dynamics of a string under geometrically nonlinear conditions. The transverse displacement $u(x, t)$ of a stiff string is governed by the following PDE (see [25], model $S_{T,4}$)

$$ \partial_t^2 u = \frac{1}{\mu} [(T \partial_x^2 - EI \partial_x^4)u - 2\mu(\eta_0 - \eta_1 \partial_x^2)\partial_t u + \frac{EA - T}{2} \partial_x (\partial_x u)^3 + \delta(x - x_{in}) f_{in}], \quad (1) $$

Here, $x \in [0, l_0]$, for some string length $l_0$ in m, and $t \in \mathbb{R}$. The various parameters, all assumed constant, are the radius $R$ in m, cross section $A = \pi R^2$ in m<fnref n="2" />, moment of inertia $I = \pi R^4 / 4$ in m⁴, density $\rho$ in kg·m⁻³, linear mass density $\mu = \rho A$ in kg·m⁻¹, string tension $T$ in N, Young's modulus $E$ in Pa, and dissipation coefficients $\eta_0$ and $\eta_1$. $f_{in}(t)$ in N, is a driving term at position $x = x_{in}$, where $\delta$ is a Dirac delta function. Simply-supported boundary conditions $u = \partial_x^2 u = 0$ are employed at the domain endpoints $x = 0$ and $x = l_0$. Zero initial conditions are assumed (so that $u(x, 0) = \partial_t u(x, 0) = 0$ for $x \in [0, l_0]$). This simple model results from a 4th order Taylor expansion and reduction to transverse displacement only of the geometrically exact string model [26]. It has been chosen for the following reasons:

* It is a problem that the baseline SAV method struggles to handle properly (see e.g. Figure 5 of [17] and associated discussion),

* It contains some key elements of musical string modelling (stiffness, frequency dependent losses),

* It is relatively simple (no longitudinal displacement, no coupling to resonators).

For an overview of more involved musical string models, see [27, 28]. Note that the time integration scheme presented in the next sections may be adapted to most of these other models after spatial discretisation.

### 2.2. Finite difference discrete model as a port-Hamiltonian system

Consider the class of system written in the port-Hamiltonian framework (see [24] for an introduction) as

$$ \underbrace{\left[ \frac{d}{dt} \mathbf{p} \right]}_{\mathbf{f}} = \underbrace{\left( \begin{bmatrix} 0 & -\mathbf{J}_0 \\ \mathbf{J}_0^\top & 0 \end{bmatrix} \right)}_{\mathbf{J}} - \underbrace{\left( \begin{bmatrix} \mathbf{R}_0 & 0 \\ 0 & 0 \end{bmatrix} \right)}_{\mathbf{R}} \underbrace{\left[ \begin{bmatrix} \mathbf{M}^{-1} \mathbf{p} \\ \mathbf{K} \mathbf{q} + \mathbf{f}_{nl}(\mathbf{q}) \end{bmatrix} \right]}_{\mathbf{e}} + \begin{bmatrix} \mathbf{G}_p \\ 0 \end{bmatrix} \mathbf{u}, \quad (2) $$

with state $\alpha = [\mathbf{p}, \mathbf{q}]^\top$ and associated flow $f$, effort $e$, input $u$, skew-symmetric interconnection matrix $\mathbf{J}$ and symmetric positive semi-definite dissipation matrix $\mathbf{R}$.

The effort vector $e$, including the nonlinear driving force $\mathbf{f}_{nl}(\mathbf{q}) = \frac{\partial E_{nl}}{\partial \mathbf{q}}(\mathbf{q})$ is defined as the gradient of the Hamiltonian (total stored energy of the system)

$$ H(\mathbf{p}, \mathbf{q}) = \frac{1}{2} (\mathbf{p}^\top \mathbf{M}^{-1} \mathbf{p} + \mathbf{q}^\top \mathbf{K} \mathbf{q}) + E_{nl}(\mathbf{q}). \quad (3) $$

Moreover, in this work, $\mathbf{M}$ and $\mathbf{K}$ are assumed to be symmetric positive definite matrices and the nonlinear potential energy $E_{nl}(\mathbf{q})$ is required to be bounded from below with $E_{nl}(0) = 0$. This ensures that $H(\mathbf{p}, \mathbf{q})$ is a Lyapunov function of the autonomous system. (The separation of the potential energy into a quadratic and non-quadratic part allows additional flexibility in the design of numerical methods satisfying an energy balance [14].) From standard port-Hamiltonian theory, system (2) satisfies the power balance

$$ \frac{d}{dt} H(\mathbf{p}, \mathbf{q}) = -\mathbf{p}^\top \mathbf{M}^{-1} \mathbf{R} \mathbf{M}^{-1} \mathbf{p} + \mathbf{p}^\top \mathbf{M}^{-1} \mathbf{G}_p \mathbf{u}. \quad (4) $$

A semi-discrete formulation of the string equations (1) is obtained using structure preserving finite differences. The general method is presented in [29] and not detailed here for brevity. It results in a semi-discrete port-Hamiltonian system of the form (2) with $\mathbf{p} \equiv [p_1, \dots, p_{N-1}]^\top$ momenta and $\mathbf{q} \equiv [q_1, \dots, q_{N-1}]^\top$ transverse displacements at the nodes of a spatial grid, with grid spacing $h$ such that $l_0 = Nh$. System matrices and operators are 2

$$ \mathbf{J}_0 = \frac{\mathbb{I}}{h}, \quad \mathbf{M}^{-1} = h \frac{\mathbb{I}}{\mu}, \quad \mathbf{K} = h(-T\mathbb{I} + EI\mathbf{D}^2)\mathbf{D}^2, \quad (5a) $$

$$ \mathbf{R}_0 = -\frac{2\mu}{h}\eta_1 \mathbf{D}^2 + \frac{2\mu}{h}\eta_0 \mathbb{I}, \quad (5b) $$

$$ E_{nl}(\mathbf{q}) = h \frac{EA - T_0}{8} ((\mathbf{D}^{-}\mathbf{q})^{\circ 2})^{\top} (\mathbf{D}^{-}\mathbf{q})^{\circ 2}, \quad (5c) $$

$$ \mathbf{f}_{nl}(\mathbf{q}) = -h \frac{EA - T_0}{2} \mathbf{D}^{+} [\mathbf{D}^{-}\mathbf{q}]^{\circ 3}, \quad (5d) $$

$$ \mathbf{u} = f_{in} \quad (5e) $$

<footnote n="2">²$\mathbf{x}^{\circ N}$ denotes element-wise exponentiation (Hadamard operation).</footnote>
with finite difference matrices

$$ \mathbf{D}^{-} = \frac{1}{h} \begin{bmatrix} 1 & & & 0 \\ -1 & \ddots & & \\ & \ddots & \ddots & 1 \\ 0 & & -1 & \end{bmatrix} \in \mathbb{R}^{N, N-1}, \quad (6a) $$

$$ \mathbf{D}^{+} = -(\mathbf{D}^{-})^{\mathrm{T}}, \quad \mathbf{D}^{2} = \mathbf{D}^{+}\mathbf{D}^{-}. \quad (6b) $$

### 3. SCALAR AUXILIARY VARIABLE TRANSFORMATION

Scalar auxiliary variable methods are based on a change of state to describe the system dynamics, used to design efficient energy stable numerical schemes. The large literature available proposes a variety of possible state transformations, presented in [30], that may be applied to port-Hamiltonian systems (also in a larger class than (2)). This work is focused on the original SAV method from [11] which has been used for the simulation of systems like (2) in the context of musical acoustics in recent years.

The scalar auxiliary variable is introduced as

$$ r(t) = \sqrt{2E_{nl}(\mathbf{q}) + C_0} \quad (7) $$

and is appended to the system state to yield $\alpha_{sav} = [\mathbf{p}, \mathbf{q}, r]^{\mathrm{T}}$. Here, $C_0 > 0$ is a gauge constant. Extended matrices

$$ \hat{\mathbf{J}} = \begin{bmatrix} \mathbf{J} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix}, \quad \hat{\mathbf{R}} = \begin{bmatrix} \mathbf{R} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} \quad (8) $$

are defined. The dynamics of system (2) may then be rewritten equivalently, using the chain rule

$$ \underbrace{\left[ \begin{array}{c} \frac{d}{dt} \mathbf{p} \\ \frac{d}{dt} \mathbf{q} \\ \frac{d}{dt} r \end{array} \right]}_{\mathbf{f}_{\mathrm{sav}}} = (\hat{\mathbf{J}} + \mathbf{J}_{\mathrm{sav}}(\alpha_{\mathrm{sav}}) - \hat{\mathbf{R}}) \underbrace{\begin{bmatrix} \mathbf{M}^{-1} \mathbf{p} \\ \mathbf{K} \mathbf{q} \\ r \end{bmatrix}}_{\mathbf{e}_{\mathrm{sav}}} + \begin{bmatrix} \mathbf{G}_p \\ \mathbf{0} \\ \mathbf{0} \end{bmatrix} u, \quad (9a) $$

where the skew-symmetric matrix $\mathbf{J}_{sav}$ is of the form

$$ \mathbf{J}_{\mathrm{sav}}(\alpha_{\mathrm{sav}}) = \begin{bmatrix} \mathbf{0} & \mathbf{0} & -\mathbf{g}(\alpha_{\mathrm{sav}}) \\ \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{g}^{\mathrm{T}}(\alpha_{\mathrm{sav}}) & \mathbf{0} & \mathbf{0} \end{bmatrix}, \quad (9b) $$

with $\mathbf{g}(\alpha_{sav}) := \mathbf{g}_{\mathrm{std}}(\mathbf{q}) = \frac{\mathbf{J}_0 f_{nl}(\mathbf{q})}{\sqrt{2E_{nl}(\mathbf{q})+C_0}}$. The effort vector $\mathbf{e}_{sav}$ is now the gradient of the Hamiltonian of the system in the new state space

$$ H_{\mathrm{sav}}(\alpha_{\mathrm{sav}}) = \frac{1}{2} (\mathbf{p}^{\mathrm{T}} \mathbf{M}^{-1} \mathbf{p} + \mathbf{q}^{\mathrm{T}} \mathbf{K} \mathbf{q} + r^2), \quad (10) $$

and the power balance remains unchanged. Systems (2) and (9) are equivalent in the sense that for consistent initial conditions satisfying $r(t_0) = \sqrt{2E_{nl}(\mathbf{q}(t_0)) + C_0}$, they yield the same trajectories in the original phase space $[\mathbf{p}, \mathbf{q}]$. Moreover, any trajectory of (9) generated from consistent initial conditions satisfies relation (7) for all $t \ge t_0$.

### 4. DRIFT CONTROL

System (9) presents a key advantage for efficient energy stable simulation, as the corresponding Hamiltonian is quadratic, as discussed in previous work [14]. However, designing efficient schemes typically requires explicit evaluation of the matrix $\mathbf{J}_{sav}(\alpha_{sav})$. As a result, drift between numerical values of $r$ and $\sqrt{2E_{nl}(\mathbf{q}) + C_0}$ may appear, even for consistent initial conditions. This can have a dramatic effect on simulations, as shown in the numerical experiments section, as the simulated SAV system is then no longer equivalent with the original system. This has been documented in the literature in e.g. [17, 23, 21]. On the studied string model, the auxiliary variable drift causes a progressive degradation of the numerical results for successive hits of the string, as can be seen on blue traces of Figure 4.

The main contribution of this paper consists of the introduction of a new method reducing the impact of the numerical drift problem. Our proposed approach to tackle this problem lies in the addition of a control term to system (9) before time discretization, guided by a set of requirements:

* (R0): The (scalar) drift measure

$$ \epsilon(\mathbf{q}, r) = r - \sqrt{2E_{nl}(\mathbf{q}) + C_0}, \quad (11) $$

is governed by the control law

$$ \dot{\epsilon} = -\lambda_0 \epsilon, \quad (12) $$

* (R1): The drift controller must not change the global power-balance. Consequently, it must be expressible as a skew-symmetric modification of the interconnection matrix $\mathbf{J}_{sav}$ in (9a),

* (R2): The dynamics of the modified system must be that of system (9) if $\epsilon = 0$.

Using (R0), the modified dynamics of $r$ may be written as

$$ \dot{\epsilon} = -\lambda_0 \epsilon, \quad (13a) $$

$$ \Leftrightarrow \dot{r} = -\lambda_0 \epsilon(\mathbf{q}, r) + \frac{\partial_q E_{nl}(\mathbf{q})^{\mathrm{T}} \dot{\mathbf{q}}}{\sqrt{2E_{nl}(\mathbf{q}) + C_0}}, \quad (13b) $$

$$ \Leftrightarrow \dot{r} = -\lambda_0 \epsilon(\mathbf{q}, r) + \frac{\partial_q E_{nl}(\mathbf{q})^{\mathrm{T}} \mathbf{J}_0^{\mathrm{T}} \mathbf{M}^{-1} \mathbf{p}}{\sqrt{2E_{nl}(\mathbf{q}) + C_0}}, \quad (13c) $$

where the second right-hand term is already present in system (9). The addition of the control term $-\lambda_0 \epsilon(q, r)$ under requirement (R1) is however not trivial as $\epsilon$ does not appear in the effort vector $\mathbf{e}_{sav}$. Here, we build an approximation to unity from the velocity $\mathbf{M}^{-1} \mathbf{p}$ as $1 = \frac{\|\mathbf{M}^{-1} \mathbf{p}\|_{L^1}}{\|\mathbf{M}^{-1} \mathbf{p}\|_{L^1}}$, with $L^1$ norm $\|\mathbf{y}\|_{L^1} = \mathrm{sign}(\mathbf{y})^{\mathrm{T}} \mathbf{y}$. Using this choice, the dynamics of the modified systems is given by system (9) with modified auxiliary variable coupling term

$$ \mathbf{g}(\alpha_{\mathrm{sav}}) := \mathbf{g}_{\mathrm{std}}(\mathbf{q}) + \mathbf{g}_{\mathrm{mod}}(\alpha_{\mathrm{sav}}), \quad (14a) $$

$$ \mathbf{g}_{\mathrm{mod}}(\alpha_{\mathrm{sav}}) = -\lambda_0 \epsilon(\mathbf{q}, r) \frac{\mathrm{sign}(\mathbf{M}^{-1} \mathbf{p})}{\|\mathbf{M}^{-1} \mathbf{p}\|_{L^1}}. \quad (14b) $$

The modified system naturally fulfills requirement (R2) as $\mathbf{g}(\alpha_{sav})$ is a first order perturbation of $\mathbf{g}_{\mathrm{std}}(\mathbf{q})$ with parameter $\epsilon(\mathbf{q}, r)$ and as $\dot{\epsilon} = 0$ if $\epsilon = 0$. In order to fulfill requirement (R1), note that the dynamics of $\mathbf{p}$ for the modified system has however been altered in cases where $\epsilon \ne 0$. To summarize, the method consists of adding $\mathbf{g}_{\mathrm{mod}}$ to the standard $\mathbf{g}_{\mathrm{std}}$ term in (9b) and modifying the time-stepping scheme accordingly (see section 5).
# 5. ENERGY-PRESERVING NUMERICAL SCHEME

An efficient energy preserving time-interleaved scheme for (9) is built using a Stormer-Verlet scheme for the linear part of the wave equation coupled to a midpoint scheme for the SAV part [14].

## 5.1. Presentation

For improved efficiency, it is customary to treat some dissipative terms implicitly and others explicitly. In the following, the dissipation matrix is decomposed as $R = R_{sv} + R_{mid}$ with $R_{mid}$ diagonal and positive semidefinite and $R_{sv}$ positive semidefinite. State variables are approximated by time series defined on interleaved grids as $p^n$, $q^{n+\frac{1}{2}}$ and $r^n$. These time series approximate the values of the continuous time functions at times $t^n = n dt$ and $t^{n+\frac{1}{2}} = (n + \frac{1}{2})dt$, where $dt = \frac{1}{sr}$ is the sample period. The forward difference operator $\delta_{t+}w^n = (w^{n+1} - w^n)/dt$ and forward averaging operator $\mu_{t+}w^n = (w^{n+1} + w^n)/2$ are used to define the scheme as

$$ \delta_{t+} \mathbf{p}^n = - \mathbf{J}_0 \mathbf{K} \mathbf{q}^{n+\frac{1}{2}} - \mathbf{R}_{sv} \mathbf{M}^{-1} \mathbf{p}^n - \mathbf{R}_{mid} \mathbf{M}^{-1} \mu_{t+} \mathbf{p}^n - \bar{\mathbf{g}}^n \mu_{t+} r^n + \mathbf{G}_p \mathbf{u}^{n+\frac{1}{2}}, \quad (15a) $$

$$ \delta_{t+} \mathbf{q}^{n-\frac{1}{2}} = \mathbf{J}_0^\top \mathbf{M}^{-1} \mathbf{p}^n, \quad (15b) $$

$$ \delta_{t+} r^n = (\bar{\mathbf{g}}^n)^\top \mathbf{M}^{-1} \mu_{t+} \mathbf{p}^n, \quad (15c) $$

where $\bar{\mathbf{g}}^n$ is a consistent explicit evaluation of $\mathbf{g}$ at timestep $n$. In this work, we use expressions given in table 1.

<table><caption>Table 1: Evaluation of $\bar{\mathbf{g}}^n$ without and with drift control.</caption><thead><tr><th>Standard SAV</th><th>SAV with drift control</th></tr></thead><tbody><tr><td>$\bar{\mathbf{g}}^n = \mathbf{g}_{\text{std}}(\mathbf{q}^{n+1/2})$</td><td>$\bar{\mathbf{g}}^n = \mathbf{g}_{\text{std}}(\mathbf{q}^{n+1/2}) + \mathbf{g}_{\text{mod}}(\mathbf{p}^n, \mu_{t+} \mathbf{q}^{n-1/2}, r^n)$</td></tr></tbody></table>

Independently of this choice, the scheme satisfies the discrete power balance

$$ \frac{E^{n+1} - E^n}{dt} = - \underbrace{(\mathbf{M}^{-1} \mu_{t+} \mathbf{p}^n)^\top \mathbf{R} \mathbf{M}^{-1} \mu_{t+} \mathbf{p}^n}_{P_{\text{diss}}^{n+1/2} \ge 0} + \underbrace{(\mathbf{M}^{-1} \mu_{t+} \mathbf{p}^n)^\top \mathbf{G}_p \mathbf{u}^{n+\frac{1}{2}}}_{\mathbf{g}^n = \mathbf{g}_{\text{std}}(\mathbf{q}^{n+1/2})} \quad (16) $$

with pseudo-energy

$$ E^n = \frac{1}{2} \left( (\mathbf{p}^n)^\top \tilde{\mathbf{M}}^{-1} \mathbf{p}^n + (\mu_{t+} \mathbf{q}^{n-\frac{1}{2}})^\top \mathbf{K} \mu_{t+} \mathbf{q}^{n-\frac{1}{2}} + (r^n)^2 \right), \quad (17) $$

and modified mass matrix

$$ \tilde{\mathbf{M}}^{-1} = [\mathbb{I} - \frac{dt}{4}(dt\mathbf{M}^{-1}\mathbf{J}_0\mathbf{K}^\top\mathbf{J}_0^\top + 2\mathbf{M}^{-1}\mathbf{R}_{sv})](\mathbf{M}^{-1}). \quad (18) $$

The scheme is stable under the condition that the pseudo-energy is non-negative, which holds if $\tilde{\mathbf{M}}^{-1}$ is positive definite. The stability condition thus depends only on a spectral analysis of the linear part of the system. Note that dissipative contributions evaluated with the midpoint effort estimate (corresponding to matrix $R_{mid}$) are purely dissipative as $P_{\text{diss}}^{n+1/2} \ge 0$, whereas dissipative contributions evaluated explicitly (corresponding to matrix $R_{sv}$) change the expression of the conserved pseudo-energy as well as the stability condition through modification of $\tilde{\mathbf{M}}^{-1}$ in (18). This is the same as the case of backward difference approximations used in the case of linear string vibration [3].

## 5.2. Update form

Even though (15) is semi-explicit, efficient update is in fact possible for systems having diagonal mass matrix $\mathbf{M}$. Indeed, given $\mathbf{p}^n$, $\mathbf{q}^{n-1/2}$ and $r^n$, time stepping may be performed following the following operations

$$ \mathbf{q}^{n+1/2} = \mathbf{q}^{n-1/2} + dt\mathbf{J}_0^\top \mathbf{M}^{-1} \mathbf{p}^n \quad (19a) $$

$$ \bar{\mathbf{g}}^n = \mathbf{g}_{\text{std}}(\mathbf{q}^{n+1/2}) + \mathbf{g}_{\text{mod}}(\mathbf{p}^n, \mu_{t+} \mathbf{q}^{n-1/2}, r^n) \quad (19b) $$

$$ \mathbf{A}^n \mathbf{p}^{n+1} = -\mathbf{J}_0 \mathbf{K} \mathbf{q}^{n+1/2} + \mathbf{B}^n \mathbf{p}^n - \bar{\mathbf{g}}^n r^n + \mathbf{G}_p \mathbf{u}^{n+1/2} \quad (19c) $$

$$ r^{n+1} = r^n + dt(\bar{\mathbf{g}}^n)^\top \mathbf{M}^{-1} \mu_{t+} \mathbf{p}^n \quad (19d) $$

with

$$ \mathbf{A}^n = \frac{\mathbb{I}}{dt} + \frac{dt}{4} \bar{\mathbf{g}}^n (\bar{\mathbf{g}}^n)^\top \mathbf{M}^{-1} + \frac{\mathbf{R}_{\text{mid}}}{2} \mathbf{M}^{-1}, \quad (19e) $$

$$ \mathbf{B}^n = \frac{\mathbb{I}}{dt} - \frac{dt}{4} \bar{\mathbf{g}}^n (\bar{\mathbf{g}}^n)^\top \mathbf{M}^{-1} - \frac{\mathbf{R}_{\text{mid}}}{2} \mathbf{M}^{-1} - \mathbf{R}_{\text{sv}} \mathbf{M}^{-1}, \quad (19f) $$

and where the linear system (19c) can be solved efficiently by noting that $\mathbf{A}^n$ is a rank 1 perturbation of a diagonal matrix. Leveraging the Sherman-Morrison inversion formula [15], one obtains

$$ (\mathbf{A}^n)^{-1} = \mathbf{M} \left[ \mathbf{A}_0^{-1} - \frac{dt}{4} \frac{\mathbf{A}_0^{-1} \bar{\mathbf{g}}^n (\bar{\mathbf{g}}^n)^\top \mathbf{A}_0^{-1}}{1 + \frac{dt}{4} (\bar{\mathbf{g}}^n)^\top \mathbf{A}_0^{-1} \bar{\mathbf{g}}^n} \right], \quad (19g) $$

$$ \mathbf{A}_0^{-1} = 2dt (2\mathbf{M} + \mathbf{R}_{\text{mid}} dt)^{-1}. \quad (19h) $$

Reduction to an update on $\mathbf{q}^{n+1/2}$ and $r^n$ uniquely is possible if matrix $\mathbf{J}_0$ is diagonal. This corresponds to cases in which (2) may be reduced to an equivalent second order ODE in $\mathbf{q}$ and is the case for the string model.

## 5.3. Stability condition for the string

The modified mass matrix for the discrete string model (6) reads

$$ \tilde{\mathbf{M}}^{-1} = \frac{h}{\mu} \left[ \mathbb{I} + \left( \frac{dt^2}{4\mu} (T\mathbb{I} - EID^2) + dt\eta_1 \mathbb{I} \right) \mathbf{D}^2 \right]. \quad (20) $$

The stability condition is directly obtained from noticing that the lowest eigenvalue of $-\mathbf{D}^2$ is equal to $-\max(\lambda_{\mathbf{D}^2}) = -4/n^2$ and writes

$$ h \ge h_{\min} = \sqrt{\frac{\gamma + \sqrt{\gamma^2 + 16\mu EIDt^2}}{2\mu}}, \quad \gamma = dt^2T + 4dt\mu\eta_1. \quad (21) $$

## 5.4. Reference scheme

In the following, a reference energy preserving semi-explicit scheme for cubic nonlinearities is used for comparison [25]. Using the same notations and operators as before, as well as $\delta_{tt}w^n = 1/dt^2(w^{n+1} - 2w^n + w^{n-1})$ and $\delta_{t\bullet}w^n = 1/2dt(w^{n+1} - w^{n-1})$, it writes

$$ \delta_{tt} \mathbf{q}^n = \frac{1}{\mu} \left[ (T\mathbf{D}^2 - EID^4) \mathbf{q}^n - 2\mu(\eta_0 \mathbb{I} - \eta_1 \mathbf{D}^2) \delta_{t\bullet} \mathbf{q}^n + \frac{EA - T_0}{2} \mathbf{D}^+ (\mathbf{D}^- \mathbf{q}^n)^2 (\mathbf{D}^- \mu_{t+} \mathbf{q}^n) \right] \quad (22) $$

for the isolated cubic nonlinear string. The stability condition is also given by (21). Updating for this scheme requires the solution of a tridiagonal linear system at each timestep.
# 6. NUMERICAL EXPERIMENTS

This section presents results of numerical experiments performed on the string model using the proposed drift-controlled SAV method. First, the numerical power balance and convergence of the scheme are illustrated. Second, the effect of the new control term is studied and the choice of parameter $\lambda_0$ discussed.

## 6.1. Power balance and convergence of the scheme

**Configuration:** A string with fixed physical parameters given in table 2 is used.

<table><caption>Table 2: Fixed physical parameters.</caption><tr><td>l<sub>0</sub> = 1.1 m, ρ = 8000 kg.m<sup>-3</sup>, T = 60 N, E = 2 × 10<sup>11</sup> Pa</td></tr><tr><td>s<sub>0</sub> = 0.9 s<sup>-1</sup>, s<sub>1</sub> = 4 × 10<sup>-4</sup> m<sup>2</sup>s<sup>-1</sup>, A = π(0.4 × 10<sup>-3</sup>)<sup>2</sup> m<sup>2</sup></td></tr></table>

The external forcing term and initial velocity are set to zero. The initial displacement is set to correspond to the first modal shape of the linear string, or

$$ \mathbf{u}^{-1/2} = \mathbf{u}^{1/2} = U_0 \left[ \sin\left(\frac{\pi}{N}\right), \dots, \sin\left(\frac{(N-1)\pi}{N}\right) \right], \quad (23) $$

with amplitude $U_0$. The auxiliary variable is initialized accordingly from the value of $\mathbf{u}^0 = \frac{1}{2}(\mathbf{u}^{-1/2} + \mathbf{u}^{1/2})$. For a given sampling frequency, $N$ is determined using the formula

$$ N = 2 \left\lfloor \frac{\kappa l_0}{2h_{min}} \right\rfloor, \quad (24) $$

with $0 < \kappa \le 1$ a tuning parameter. Note that $N$ is forced to be even, as this ensures that the displacement at the middle of the string can be observed as $\hat{u} = \mathbf{u}_{N/2}$ (scalar).

**Power balance:** A simulation with $U_0 = 0.01$ m is run at $sr = 44100$ Hz and $h = h_{min}$. Figure 1 shows the relative error on the energy balance, obtained as the error on (16) multiplied by $dt/E^0$. The error is on the order of machine precision in double-precision floating point arithmetic (error values closer to $10^{-15}$ could probably be obtained with proper scaling of the equations, which was not performed in this work), which shows that the scheme is indeed energy stable for quantities defined in the previous section. From (16), this result does not depend on the new control parameter $\lambda_0$ or on the initial amplitude $V_0$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Numerical error on relative energy balance with $\lambda_0 = 1000$ s$^{-1}$.</figcaption></figure>

**Convergence:** A set of $B$ sampling frequencies $\{sr_0, 2sr_0, \dots, 2^{B-1}sr_0\}$ Hz is used for simulations, with $sr_0 = 20000$ Hz. A reference solution $\mathbf{u}_{ref}$ is built using the reference

scheme (22) with $sr = 4 \times 2^{B-1}sr_0$ Hz. For each sampling frequency, $N$ is computed using (24), with $\kappa = 0.9$ or $\kappa = 1$. Solutions are computed for a 0.1s duration using the reference algorithm and the proposed algorithm with $\lambda_0 = \{0, 1000\} \text{ s}^{-1}$. Note that the case $\lambda_0 = 0$ corresponds to the explicit time-interleaved SAV scheme recently reported in the literature [14]. The relative $L^2$ error (intended on time series of scalar quantities) on the displacement of the middle of the string is used and computed as

$$ e(\hat{u}, \hat{u}_{ref}) = \frac{\|\hat{u} - \hat{u}_{ref}\|^2}{\|\hat{u}_{ref}\|^2}, \quad (25) $$

for each solutions. Figure 2 displays convergence curves for $U_0 = 8 \times 10^{-3}$m and $\kappa = 0.9$. In this case, all schemes seem to converge with a second order slope. However, this result should be proved (or disproved) since the numerical scheme is second order accurate when there is no drift but the drift rejection is operated with a formally first order accurate scheme.

<figure><img src="image_3.png" /><figcaption>Figure 2: Convergence study for $U_0 = 8 \times 10^{-3}$m and $\kappa = 0.9$, with $\lambda_0 = 0s^{-1}$ and $\lambda_0 = 1000s^{-1}$. The grey dashed line indicates the second order convergence slope.</figcaption></figure>

Figure 3 displays the same convergence curves for $\kappa = 1$. In this case, convergence is not observed for the studied range of samplersate values.

<figure><img src="image_4.png" /><figcaption>Figure 3: Convergence study for $U_0 = 8 \times 10^{-3}$m and $\kappa = 1$, with $\lambda_0 = 0s^{-1}$ and $\lambda_0 = 1000s^{-1}$. The grey dashed line indicates the second order convergence slope.</figcaption></figure>

For linear systems, it is known that $\kappa = 1$ maximizes output bandwidth. However, for nonlinear systems, this is not the case as the time stepping scheme may introduce some spectral aliasing, which is reduced if the eigenvalues of the discrete system are farther from Nyquist's frequency [3]. This indicates that it might be preferable to slightly reduce $\kappa$ and over-sample the simulation such that eigenfrequencies of the linear part are represented up to the hearing limit $f = 20$ kHz, while reducing unwanted aliasing.
## 6.2. Long term behaviour and effect of the new control term

The effect of the control term $\lambda \ge 0$ is not clearly visible on convergence curves. However, it greatly influences the long term behaviour of the simulations. Figure 4 presents an analysis of the output of the real-time implementation in Max/MSP. The string is set to have a fundamental frequency of $f_0 = 80$ Hz and is excited every second by a half cosine force of period 1ms and maximum amplitude 1N. The top row displays a time-frequency analysis of the output velocities for $\lambda_0 = 0s^{-1}$ and $\lambda_0 = 1000s^{-1}$. The bottom row shows the evolution of the relative SAV drift value.

The original explicit SAV algorithm from [14] leads to a severe drift of the scalar auxiliary variable which is clearly observed on the bottom row and reaches up to 10000 relative error. This results in clear audible artifacts that progressively degrade after each consecutive excitation. These are audible as spurious pitch bends and parasitic high frequency noise. The proposed method successfully keeps the relative drift value low. As a result, synthesis quality does not degrade after several hits: pitch and harmonic content remain stable. Note that the slight up bend observed at each impact is an expected effect due to the nonlinearity of the string model.

The tuning the value of $\lambda_0$ remains an open question. Indeed, it must ideally be as high as possible to be more reactive to a sudden increase in the value of $\epsilon$, possible during high amplitude transients. However, it is important to note that the numerical realization of the control law (12) is not guaranteed to be converging, even though the simulation itself is stable in practice, divergence of $\epsilon$ was only observed for $\lambda_0 \ge sr$.

## 7. REAL-TIME APPLICATION

Due to the explicit interleaved design, the algorithm presented here is efficient enough to be run in real-time audio environments. This section presents the implementation of the model as a Max/MSP external. It allows the instantiation of an object representing a single string with given physical parameters. Using multi-threading, 60 note polyphony was obtained for sample-rate $sr = 96000$ Hz. First, the framework used for this implementation is quickly presented and then, the computational efficiency is studied.

### 7.1. Framework and general presentation

The algorithm is implemented in C++ using the Eigen library for vector arithmetic. All matrix products are reduced to vector operations, noting that that difference matrices are at most pentagonal. The MinDevKit C++ API is used to define and build a Max/MSP object, represented in Figure 5. An example user interface built around this object is shown in Figure 6. It presents the user with 5 audio-rate input signals and 3 audio-rate output signals.

The string model is excited locally by the scalar force signal $F_{in}$ at position $x = x_{ex}/l_0$. $v_L := \dot{u}(x_L, t)$ and $v_R := \dot{u}(x_R, t)$ are velocity signals drawn from the string at two independent listening points ("pickups"), providing an artificial stereo audio output. First-order interpolation between two adjacent grid points is performed, both for excitation and listening positions. A pitch bend signal is also provided to the object. The scalar auxiliary variable drift $\epsilon$ is also generated as an output, mostly for monitoring and demonstration purposes and qualitative assessment of the effect of $\lambda_0$.

<figure><img src="image_5.png" /><figcaption>Figure 5: Max/MSP object diagram.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 6: Screenshot of the proposed Max/MSP interface for interacting with a single string.</figcaption></figure>

### 7.2. Physical parameters

The physical constants of the string are not presented directly to the user. Instead, they are derived from a set of higher level perceptual parameters. The fundamental frequency $f_0$ and inharmonicity factor $\beta$ are related to physical parameters using the relation

$$f_n \approx (n + 1) \frac{1}{2l_0} \sqrt{\frac{T}{\mu}} \sqrt{1 + \beta(n + 1)^2}, \quad \beta = \frac{\pi^2 EI}{Tl_0^2}, \quad (26)$$

which comes from an analysis of the eigen-frequencies of the linear stiff string model. $\rho$, $R$ and $E$ are fixed to a set of values corresponding to the string material. From (26), $T$ and $l_0$ are obtained as

$$T = 2\pi f_0 \sqrt{\frac{EI\mu}{\beta(1 + \beta)}}, \quad l_0 = \frac{\pi}{2} \sqrt{\frac{EI}{\beta T}}. \quad (27)$$
<figure><img src="image_7.png" /><figcaption>Figure 4: Top: evolution of partial frequencies for repeated strikes on the string. Bottom: evolution of the relative SAV drift value. Orange: with anti-drift ($\lambda_0 = 1000s^{-1}$, see (12)), blue: without anti-drift ($\lambda_0 = 0s^{-1}$).</figcaption></figure>

The two dissipation coefficients $\eta_0$ and $\eta_1$ are computed from decay times given at two frequencies $t_{60}(f_{60}^0)$ and $t_{60}(f_{60}^1)$ using formulas given in [3, Chapter 7].

### 7.3. Modulations

The first dissipation coefficient $\eta_0$ may be modulated in real-time without interfering with the stability condition as it does not intervene in (21). Modulation of this parameter allows the modification of bulk general damping by shifting the decay time curve vertically without changing its shape.

Real-time modulation of the other physical parameters is not straightforward as the stability condition depends on their values. In this work, small modulations of $l_0$ modulations (and consequently $h$) are used to implement the pitch bend. An objective fundamental frequency is derived from the pitch bend signal (in cents) as $f_{bend} = f_0 (\frac{\text{bend}}{1200})^2$. From this modified frequency and other physical parameters, a modified length $l_{bend}$ is determined using (26). Finally, $h$ is computed as

$$h = \max \left( \frac{l_{bend}}{N}, h_{\min} \right), \quad (28)$$

ensuring the stability of the simulation. Note that the maximum up-bend depends on the value of $\kappa$ chosen at initialization.

### 7.4. Computational efficiency

Table 3 presents the real-time ratios of the C\* implementation for different system sizes and floating point number representations. Values are obtained as the ratio of the computing time and of the effective computed output duration. Simulations were run on a MacBook Pro M3 with $sr = 44100$ Hz for 20 s of effective output signal. Results for other sampling rates are not provided here

for brevity but may be obtained simply by executing the timing script of the provided code. Compared to the original explicit SAV scheme from [14], the proposed method is near 40% more expensive, which is expected as the nonlinear function must be evaluated twice to obtain $\bar{g}^n$ in equation (19b). However, the computational cost is still well below the real time limit.

<table><caption>Table 3: Real time ratios of the C\* Eigen implementation run on a MacBook Pro M3 for $sr = 44100$ Hz. Different fundamental frequencies and sampling rates are used. Ratios are given for both 32-bit floating point and 64-bit floating point numbers. The last column presents the results for the original explicit SAV method from [14], without the cost of the additional control term computation.</caption><thead><tr><th>$f_0$ (Hz)</th><th>N</th><th>32 bit</th><th>64 bits</th><th>64 bits original SAV</th></tr></thead><tbody><tr><td>20</td><td>360</td><td>3.6 %</td><td>7.3 %</td><td>5.3 %</td></tr><tr><td>50</td><td>215</td><td>2.2 %</td><td>4.1 %</td><td>3 %</td></tr><tr><td>100</td><td>140</td><td>1.5 %</td><td>2.7 %</td><td>2 %</td></tr><tr><td>200</td><td>84</td><td>0.99 %</td><td>1.6 %</td><td>1.2 %</td></tr></tbody></table>

## 8. CONCLUSION

The proposed method extends the explicit time-interleaved SAV scheme presented in [14] to overcome the auxiliary variable drift problem. It is based on the new idea of the introduction of a additional control term in the continuous time equations.

We have presented the results of an application of this method to the simulation of a nonlinear musical string. As with the original SAV scheme, general convergence results remain to be proved. From a practical point of view, the proposed modification enhances the long-term behaviour of the simulation. Specifically, it prevents spurious pitch drifts from building up as the string is struck multiple times. The computational efficiency remains relatively high, making the algorithm useful for real-time synthesis. The proposed
Max/MSP implementation demonstrates a simple digital instru-
ment based on the simulation of 60 uncoupled nonlinear strings
excited by synthetic force signals. A more realistic setting may be
achieved by modelling the exciter (i.e. hammer or bow) and cou-
pling the strings together through a bridge. Note that these com-
ponents may be formulated considering a passive coupling under
a port-Hamiltonian representation, such that the general method-
ology proposed in this paper would still be valid.

Future work includes more general numerical experiments on
a broader category of dynamical systems. In the context of musical
instrument modelling, collision problems are usually the most dif-
ficult to deal with and would present the greatest challenge for the
proposed method. Tests on two-dimensional systems, like nonlin-
ear plate vibrations are also of interest for e.g. drum or gong sound
synthesis (see [19], where the explicit SAV method has been used
for this purpose, and where drift does pose a major problem). From
a theoretical point of view, this paper presents one of many possi-
ble solutions for the discrete evaluation of the proposed additional
control term for which careful analysis may lead to refined choices.
Other control laws may also be considered (by e.g. adding an in-
tegral control term).

9. ACKNOWLEDGMENTS

This research was funded, in whole or in part, by l'Agence Na-
tionale de la Recherche (ANR), project AVATARS (ANR-22-
CE48-0014). We wish to thank Benjamin Matuszewski for his
help on setting up the github project, Victor Paredes for his partic-
ipation in the design of the example Max/MSP patches and Charles
Picasso for helping to produce the top row of Figure 4.

10. REFERENCES

[1] J. O. Smith III, "Principles of digital waveguide models of musical instruments," in *Applications of digital signal processing to audio and acoustics*, pp. 417–466. Springer, 1998.

[2] R. Caussé, J. Bensoam, and N. Ellis, "Modalys, a physical modeling synthesizer: More than twenty years of researches, developments, and musical uses," *J. Acoust. Soc. Am.*, vol. 130, no. 4_Supplement, pp. 2365–2365, 2011.

[3] S. Bilbao, *Numerical sound synthesis: finite difference schemes and simulation in musical acoustics*, John Wiley & Sons, 2009.

[4] E. Hairer, C. Lubich, and G. Wanner, "Structure-preserving algorithms for ordinary differential equations," *Geometric numerical integration*, vol. 31, 2006.

[5] G. Quispel and G. Turner, "Discrete gradient methods for solving odes numerically while preserving a first integral," *J. Physics A: Math. Gen.*, vol. 29, no. 13, pp. L341, 1996.

[6] R. McLachlan, G. Quispel, Reinout W, and N. Robidoux, "Geometric integration using discrete gradients," *Phil. Trans. Royal Soc. London. Series A: Math., Phys. Eng. Sci.*, vol. 357, no. 1754, pp. 1021–1045, 1999.

[7] E. Eich-Soellner and C. Führer, *Numerical methods in multibody dynamics*, vol. 45, chapter Implicit ordinary differential equations, pp. 139–192, Springer, 1998.

[8] N. Lopes, T. Hélie, and A. Falaize, “Explicit second-order accurate method for the passive guaranteed simulation of port-Hamiltonian systems,” in *5th IFAC Workshop on Lagrangian and Hamiltonian Methods for Nonlinear Control LHMNC 2015*, Lyon, France, July 2015, vol. 48, pp. 1–.

[9] F. Marazzato, A. Ern, C. Mariotti, and L. Monasse, "An explicit pseudo-energy conserving time-integration scheme for hamiltonian

dynamics,” *Comp. Meth. Appl. Mech. Eng.*, vol. 347, pp. 906–927, 2019.

[10] J. Zhao, Q. Wang, and X. Yang, “Numerical approximations for a phase field dendritic crystal growth model based on the invariant energy quadratization approach,” *Int. J. Num. Meth. Eng.*, vol. 110, no. 3, pp. 279–300, 2017.

[11] J. Shen, J. Xu, and J. Yang, “The scalar auxiliary variable (sav) approach for gradient flows,” *J. Comp. Phys.*, vol. 353, pp. 407–416, 2018.

[12] S. Sato, "High-order linearly implicit exponential integrators conserving quadratic invariants with application to scalar auxiliary variable approach," *Num. Alg.*, vol. 96, no. 3, pp. 1295-1329, 2024.

[13] M. Ducceschi and S. Bilbao, "Simulation of the geometrically exact nonlinear string via energy quadratisation," *J. Sound Vib.*, vol. 534, pp. 117021, 2022.

[14] S. Bilbao, M. Ducceschi, and F. Zama, "Explicit exactly energy-conserving methods for hamiltonian systems," *J. of Comp. Phys.*, vol. 472, pp. 111697, 2023.

[15] J. Sherman and W. J. Morrison, "Adjustment of an inverse matrix corresponding to a change in one element of a given matrix," *Ann. Math. Stat.*, vol. 21, pp. 124–127, 1950.

[16] M. Ducceschi, M. Hamilton, and R. Russo, "Simulation of the snare-membrane collision in modal form using the scalar auxiliary variable (sav) method," in *Forum Acusticum 2023*, 2023, pp. 1–7.

[17] R. Russo, S. Bilbao, and M. Ducceschi, "Scalar auxiliary variable techniques for nonlinear transverse string vibration," *IFAC-PapersOnLine*, vol. 58, no. 6, pp. 160–165, 2024.

[18] M. van Walstijn, V. Chatzioannou, and N. Athanasopoulos, "An explicit scheme for energy-stable simulation of mass-barrier collisions with contact damping and dry friction," *IFAC-PapersOnLine*, vol. 58, no. 6, pp. 214–219, 2024.

[19] S. Bilbao, C. Webb, Z. Wang, and M. Ducceschi, "Real-time gong synthesis," in *Proc. 26th Int. Conf. Digital Audio Effects*, 2023, pp. 1–8.

[20] G. Castera and J. Chabassier, "Numerical analysis of quadratized schemes. application to the simulation of the nonlinear piano string," Tech. Rep., 2023.

[21] Maarten Van Walstijn, Vasileios Chatzioannou, and Abhiram Bhanuprakash, "Implicit and explicit schemes for energy-stable simulation of string vibrations with collisions: Refinement, analysis, and comparison," *Journal of Sound and Vibration*, vol. 569, pp. 117968, 2024.

[22] M. Jiang, Z. Zhang, and J. Zhao, "Improving the accuracy and consistency of the scalar auxiliary variable (sav) method with relaxation," *J. Comp. Phys.*, vol. 456, pp. 110954, 2022.

[23] Q.-A. Huang, C. Yuan, G. Zhang, and L. Zhang, "A computationally optimal relaxed scalar auxiliary variable approach for solving gradient flow systems," *Comp. & Math. with Appl.*, vol. 156, pp. 64–73, 2024.

[24] A. van Der Schaft, "Port-hamiltonian systems: an introductory survey," in *Proc. Int. Congr. Math.*, 2006, vol. 3, pp. 1339–1365.

[25] S. Bilbao, "Conservative numerical methods for nonlinear strings," *J. Acoust. Soc. Am.*, vol. 118, pp. 3316–3327, 11 2005.

[26] P. Morse and U. Ingard, *Theoretical Acoustics*, Princeton University Press, Princeton, New Jersey, 1968.

[27] S. Bilbao and M. Ducceschi, "Models of musical string vibration," *Acoust. Sci. Tech.*, vol. 44, no. 3, pp. 194–209, 2023.

[28] J. Chabassier, *Modélisation et simulation numérique d'un piano par modèles physiques*, Ph.D. thesis, Ecole Polytechnique X, 2012.

[29] V. Trenchant, H. Ramirez, Y. Le Gorrec, and P. Kotyczka, “Finite differences on staggered grids preserving the port-hamiltonian structure with application to an acoustic duct,” *J. Comp. Phys.*, vol. 373, pp. 673–697, 2018.

[30] Q. Cheng, C. Liu, and J. Shen, "Generalized sav approaches for gradient systems," *J. Comp. Appl. Math.*, vol. 394, pp. 113532, 2021.