Linearly-implicit schemes for collisions in musical acoustics based on energy quadratisation^a)

Collisions play a key role in the operation of many musical instruments. The most obvious examples are the hammer-string¹ and mallet-membane interactions,² but there are many others: fret/string interactions in instruments such as the guitar;^3,4 reed-beating effects in wind instruments;^5,6 the sitar⁷ and tanpura;⁸ and wire/membrane collisions in the snare drum.⁹ Some collisions may be modelled as lumped, and considered to act only over a very small portion of a system (e.g., a piano hammer). Others are distributed in spatial extent such as the wire-membrane interaction. Furthermore, some of these collisions involve obstacles that are conveniently modelled as rigid (e.g., a fretboard), while in others, the effects of deformation are critical. The collision force is strongly nonlinear and cannot be approximated through linearisation.

At the numerical level, various approaches are available. Unilaterally-constrained dynamics may be used to model the collision of a vibrating object, such as a string, against a rigid, immovable obstacle;^7,10,11 non-smooth dynamical representations have also been employed for the same purpose.^12,13 In contrast, when the colliding objects are deformable, a common approach is to model the interaction via a suitable potential function: under perfectly elastic conditions, the collision energy is exchanged while remaining conserved overall.¹⁴ The potential function depends on the amount of deformation of the colliding objects, often in power-law fashion.¹⁵ Such models may then be extended to include losses that may take place during the collision.¹⁶ The possibility of modelling a collision via energy methods is particularly attractive from a numerical design perspective since this passivity property can be used as a condition on stability.¹⁷ Thus, such potential-based methods have been extended to cases involving rigid obstacles, though the interpenetration is now interpreted as a penalty.^8,17 At the numerical level, energy conservation may be achieved via schemes involving the solution of a system of nonlinear equations. Though the existence and uniqueness of the underlying solutions have been proven,^8,17 the resulting numerical schemes can only be approached via iterative root-finding routines such as e.g., Newton-Raphson. Furthermore, for collisions taking place in systems of finite spatial extent, approaches based on modal decompositions are impaired due to the implicit character of the update equations,¹⁸ and efficient solutions are only available in the case of linear barrier force.¹⁹ 

In this work, a method is presented such that the resulting numerical schemes maintain a notion of passivity, via energy conservation, while avoiding iterative methods; there is at most the solution of one linear system per update. On top of the reduction in computational cost, in this case, existence and uniqueness of solutions follow in an obvious manner. Such schemes are based on the quadratisation of the collision potential energy, through the introduction of an auxiliary function treated as a new additional state variable. Quadratisation strategies allowing for explicit numerical updates appeared first in the context of Port-Hamiltonian systems,^20,21 for invertible potentials. The introduction of an additional state variable was proposed within the context of the Invariant Energy Quadratisation^22,23 method. In this work, the case of a non-invertible potential is considered. For distributed collisions, the proposed schemes generalise naturally to the case of modal-based discretisations, yielding an efficient update.

The article is organised as follows: Sec. II introduces the proposed schemes for the simple case of the mass-barrier collision. In this section, the proposed scheme is compared against benchmark schemes borrowed from the literature.^8,17 Section III presents the case of the hammer-string interaction using a finite difference discretisation on the spatial operators. Section IV extends the previous example to a fully distributed barrier and by making use of both finite difference and modal schemes. Finally, Secs. V and VI present applications of the proposed schemes for the cases of the snare drum and the tromba marina.

In this section, the case of a point mass colliding against a rigid barrier is presented. This section also introduces the temporal difference operators that will be used throughout the text.

The mass-barrier collision, though not a musical system per se, serves as an introductory test case, from which the properties of the numerical schemes and their main operational principles can be understood. The motion is expressed through an ordinary differential equation (ODE) of the following form:

Here, M represents the mass of the particle, K is the stiffness coefficient of the linear restoring force, and $u=u(t)$ is the displacement measured from the rest position, and dependent on time $t≥0$⁠. η represents the distance between the mass and a barrier, located at z. Time differentiation is here indicated with dots, and the gradient is taken with respect to η, as suggested by the gradient subscript. Here, linear stiffness has been separated out from the general potential, as it may often be advantageous to approach discretisation through such a splitting. Equation (1) must be complemented by two initial conditions $u(0)=u0, u̇(0)=v0$⁠. From these, one has $η0=u0−z$⁠. Finally, $ϕ=ϕ(η)$ is a nonlinear potential, that for collisions takes the form

where $Kη$ is a stiffness parameter, and where $α≥1$⁠.

Using the chain rule, when the barrier height z is constant, and assuming $η≠0,u̇≠0$⁠, one may write Eq. (1) as

Upon multiplication of Eq. (3) by $u̇$⁠, one arrives at

and thus energy remains constant,

Furthermore, boundedness of the solution follows, since $ϕ(u)≥0$⁠. Under this condition, from Eq. (4),

and thus u and $u̇$ are bounded in terms of the initial energy $H0$⁠. As will be seen shortly, an iterative, conservative scheme may be derived as a discretisation of Eq. (3).

Following Ref. 8, the same motion may be described using Hamilton's equations,

The system in Eq. (6) conserves the energy $H=T+V$⁠. Here, $T=p2/2M$ is the kinetic energy, $V=Ku2/2+ϕ(u)$ is the potential energy, and p is the momentum of the particle. This system may also be discretised directly, as will be seen shortly.

Non-iterative time discretisations follow from a change of variables applied to the potential function in Eq. (1). Consider a quadratisation of the potential function $ϕ$⁠, as

One may substitute such form in the expression for the energy [Eq. (4)], obtaining

Notice that, under the condition of non-negativity of $ϕ$⁠, one may always perform such a substitution. This form of the energy includes quadratic terms only. Performing time differentiation of Eq. (8), one obtains the following equation of motion:

Formally, Eqs. (1) and (9) are entirely equivalent. They yield the same solution u(t), as well as the same bounds on the growth of such solution and its time derivative. Notice that bounds of Eq. (5) hold in this case too. Looking towards discrete time implementation (see Sec. II C 3), it is useful to rewrite Eq. (9) as

Quadratisation strategies appeared in various other contexts, such as e.g., Port-Hamiltonian systems,^20,21 and fluid dynamics,²³ and they form the core of the Invariant Energy Quadratisation method,^22,23 which is similar to the method proposed here, in that an extra auxiliary state variable is defined.

The finite difference method is employed for the simulation of the nonlinear equations. Thus, the continuous function u(t) is approximated at the time nk by the time series uⁿ, where $n∈ℕ0$⁠, and where k is the time step (and $1/k$ is the sample rate.) The basic operators in discrete time are the identity and shift operators, defined as

From these, one may define the time difference operators, all approximating the first time derivative, as

An approximation to the second time derivative is constructed from the above as

Averaging operators are also used throughout the text and are

Analogous definitions of the difference and averaging operators hold for time series defined at interleaved time instants $n−1/2$⁠. Finally, an identity used throughout the text is given here as

The finite difference formalism introduced in the previous section is now used to construct three conservative schemes. The first scheme is taken from Chatziioannou and van Walstijn.⁸ This is a method making use of a discretisation of Hamilton's equations given in first order form, and will be labelled IT-1. The second scheme is taken from Bilbao et al.,¹⁷ discretising directly Eq. (1), and labelled IT-2. Finally, the proposed scheme follows from a discretisation of the quadratised Eq. (8), and will be labelled N-IT.

Following Chatziioannou and van Walstijn,⁸ one may discretise Hamilton's Eq. (6) as

$a=un−1/2, b=kpn−1/2/M, c=k2/2M$⁠. Energy conservation arises naturally from Eq. (16), as

where

Notice that this discrete energy is non-negative by definition, reflecting the implicit nature of the discretisation of the linear part of Eq. (16). Existence and uniqueness of the solution may be shown for Eq. (16), and the resulting update Eq. (17) may be approached via a suitable root-finding algorithm, such as Newton-Raphson.⁸ 

A suitable discretisation of Eq. (3) is given in Bilbao et al.¹⁷ as

where

At each time step, the update may be written as a nonlinear function in the unknown $s=un+1−un−1$⁠,

where here $a=un−1, b=2M(δt−un)/k−Kun, c=k2/M$⁠. One then solves G(s) = 0 using a nonlinear root finder such as e.g., Newton-Raphson.

This scheme conserves a discrete energy. To see this, it is enough to multiply Eq. (19) by $δt·u$⁠, to get

where the discrete energy has the form

The discrete energy is not necessarily non-negative. One may easily show that a condition for non-negativity of the total energy is obtained for²⁴ 

This serves as a necessary and sufficient condition for stability for the scheme in Eq. (19) and is independent of the state u and of the particular form of the nonlinear function $ϕ$⁠, provided it is non-negative. Under such a condition, a discrete counterpart to the bounds [Eq. (5)] may be derived as

Existence and uniqueness may be proven for this scheme as well.¹⁷ 

Turning now to Eq. (9), and the form given in Eq. (10), a particular discretisation is given by the following system:

With this in mind, the scheme in Eq. (24) has a completely explicit form. Furthermore, inserting Eq. (24b) into Eq. (24a) and multiplying by $δt·u$ leads to a discrete energy balance,

The discrete Hamiltonian has the form

It is immediate to verify that the nonlinear potential energy is non-negative, and thus the stability condition, Eq. (22), and bounds, Eq. (23), hold in this case too.

One important aspect pertains to the choice of the explicit gradient gⁿ. Previous preliminary works employed the following form:

However, it was observed that some spurious oscillations are obtained under such choice.^27–29 A better approximation, employed in the remainder of this work, is given by

where κ = 1 if $ψn−1/2≥0$⁠, and $κ=−1$ otherwise. Once gⁿ is computed, one must also check that $gnηn−1<4ψn−1/2$⁠. If this condition is violated, then gⁿ is set to zero. This procedure ensures that the collision force is directed outwardly. Furthermore, $η⋆=u⋆−z$ is the update of the system in the absence of the collision potential, i.e.,

The particular form for Eq. (27b) can be derived by considering an implicit realisation for g, as

In the event of no collision at $n+1/2$⁠, then $gimp$ reduces to Eq. (27b), since $ψn+1/2=0$⁠, and since $un+1$ would be obtained as the solution of the system for zero collision force, as per Eq. (3). Notice as well that g in Eq. (27) is given entirely from previous values of the time series u, ψ, making the scheme fully explicit.

The existence and uniqueness of the numerical solution, regardless of the particular form of gⁿ, follow immediately, as the system is solved by simple division.

As a first experiment, consider Fig. 1. In the figure, the three schemes are compared against each other for a barrier of increasing stiffness. It can be appreciated that the schemes return consistent solutions: as the barrier stiffness is increased, the interpenetration becomes smaller. The numerical energy is conserved to the order of machine accuracy for all the three schemes. Notice as well that, as the stiffness of the barrier is increased, more iterations of the Newton-Rapshon algorithm are needed for the iterative schemes, while the computational cost of N-IT remains fixed. This is an important aspect in view of any real-time implementation requiring a precise allocation of computing resources. In fact, while it is possible to estimate the upper bound on the number of iterations required for Newton-Raphson,²⁵ the iterative routine may be affected by poor convergence, or instability in certain cases,^25,26 if, for instance, the initial guess is not carefully estimated, or the if value of the barrier stiffness is too large.

The experiment in Fig. 2 reports the convergence of the numerical schemes, computed against a reference analytic solution for a barrier with linear restoring force. Second-order accuracy is maintaned for lower stiffness values. However, for values of $Kη$ such that $Kη/M k>1$⁠, the schemes become first-order accurate, as proven by Taylor-expanding the schemes about t_n = kn.

As a first example of a collision typical of musical instruments, the hammer-string collision is investigated here. The interpenetration in this case may be interpreted as the compression experienced by the hammer felt during contact with the string. Many works have employed the power law [Eq. (2)] as a model for this case, though not all include a conservative discretisation of the resulting dynamics.^1,30,31

In a basic configuration, the system may be described by the following coupled differential equations:

The string is assumed initially at rest, and is fixed at the two ends, i.e., $u(0,t)=u(L,t)=0∀t$⁠, where it was assumed that $x∈[0,L]$⁠, with L being the string's length. The hammer has initial displacement U₀ and initial velocity V₀.

The system in Eq. (28) is probably insufficient as a musical model as such. It is lacking several important features such as stiffness³³ and losses,³⁴ and perhaps a nonlinearity inherent to the string's geometrical stretching (see, e.g., Bilbao,³⁵ as well as Morse and Ingard,³⁶ Chap. 14), though all such features may be added into the model without substantial changes to the template schemes presented here. For the purpose of illustration, they are therefore neglected at this stage, and one may refer to the case study on the tromba marina given in Sec. VI for a working example of a complete system.

The system in Eq. (28) is conservative, with energy given by

where the L₂ norm notation is used, see, e.g., Ref. 17.

Quadratisation of the energy may be performed in the same fashion as Eq. (8), yielding

where again $ψ2=2ϕ$⁠. The associated equations of motion read

Discrete realisations of both Eqs. (28) and (31) are given here in terms of appropriate finite difference schemes. The temporal finite difference operators and notation are as given in Sec. II B. Here, because of the distributed character of the string, it is convenient to introduce a matrix-vector formalism for the spatial difference operators. Thus, the string is divided into N subintervals by means of N + 1 grid points including the end points. Each subinterval is of length h, the grid spacing. The displacement u(x, t) is approximated by the grid function $umn$⁠, where n is the time step, and m is the grid index.

In a vector notation, one may then denote the grid function as $un∈ℝN−1$⁠, where the dimensionality reflects the fact that, under fixed end conditions, the end points need not be stored or updated. Spatial difference operators may then be realised as matrices. The first difference operator is given as

Thus, $D$ is a $N×(N−1)$ rectangular matrix. From this, the second difference operator is constructed simply as

yielding a square $(N−1)×(N−1)$ matrix.

The impact spatial distribution may as well be given as a vector. Hence $δ(x−xc)→r∈ℝN−1$⁠. Denoting $mh=floor(xc/h), ν=xc/h−mh$⁠, one may set $rmh=(1−ν)/h, rmh+1=ν/h$⁠, thus effectively employing linear interpolation.²⁴ 

Two finite difference schemes are given here for the solutions of Eqs. (28) and (31).

A discretisation of Eq. (28) follows immediately from this formalism as¹⁷ 

which is clearly a discrete counterpart of Eq. (29). While the nonlinear potential energy is non-negative, the linear part of the discrete energy is non-negative only under the following Courant–Friedrichs–Lewy (CFL) condition:^24,37

Non-negativity of the energy overall allows to derive a bound on the growth of the norms of the grid functions, thus effectively ensuring stability.

In order to solve Eq. (34), one first computes the update of the interpenetration η, by multiplying Eq. (34a) by $hr⊺$ and subtracting Eq. (34b), thus effectively projecting the dynamics onto the collision point. This results in an implicit scalar equation of the same form as Eq. (20), i.e., G(s) = 0, where here $s=ηn+1−ηn−1, a=ηn−1, b=2ηn−2ηn−1−T0k2h/ρr⊺D(2)un, c=k2h/ρr⊺r+k2/M$⁠. Existence and uniqueness of this nonlinear algebraic equation are proven using the same arguments as before.^8,17 Once the interpenetration is known, one may update Eqs. (34a) and (34b) explicitly.

The system in Eq. (31) may be approximated by the following scheme:

thus effectively discretising Eq. (30). The same arguments on stability as IT-2 apply here, i.e., the nonlinear energy is clearly non-negative, and for stability, the CFL condition in Eq. (35) is a necessary and sufficient condition.

A solution to Eq. (36) may be found by using the identity from Eq. (15) in both Eqs. (36a) and (36b), and then expressing the time difference of ψ using Eq. (36c), so as to effectively reduce the two equations to a linear system. The system is written as

where $b1=2un−un−1+(T0k2/ρ)D(2)un, b2=2Un−Un−1$⁠, and $a=−(g2/4)ηn−1+gψn−1/2$⁠. This linear system may be solved block-wise. Notice in particular that the matrix yields itself to a fast inversion since it is a rank-1 perturbation of the identity matrix.³⁸ However, efficient inversion techniques of the system in Eq. (37) will not be explored in this article.

Once Eq. (37) is solved, one may compute $ηn+1$ via Eq. (36d) and then update ψ via Eq. (36c).

As a first illustrative example, consider Fig. 3: here, the snapshots of IT-2 and N-IT are plotted. The hammer here has a high enough stiffness to serve as a test case.

The same dynamics may be represented in terms of time series at one output point, rather than as snaphsots. This is done in Fig. 4. It is seen that, for lower values of the hammer stiffness, the solutions of the two schemes are perfectly superimposed. As the stiffness is increased, some small differences are noticed, though the two schemes converge to the same solution in the limit of high sample rate. Energy is again conserved to the order of machine accuracy during collision.

When collisions are distributed (i.e., taking place across a spatially extended portion of a system), it may be convenient to think of a density collision potential. Here, collisions of the string against an immovable, distributed obstacle are considered. One then has $z=z(x)$⁠, and $η(x,t)=z(x)−u(x,t)$⁠. The equation of motion for the string may then be written as

where now $ϕ$ has dimension of J/m (i.e., it is a potential density.) The associated energy now reads

which is a distributed generalisation of Eq. (29). This form of the energy lends itself naturally to quadratisation. Using again $ψ2=2ϕ$⁠, one gets

with an associated equation of motion,

In this section, two applications of the non-iterative scheme are presented: one making use of a time-space finite difference scheme, and one making use of a modal projection for the spatial part. Iterative conservative finite difference schemes for the string in contact with a distributed barrier have been employed,¹⁷ based on the model Eq. (38), and extended to the case of frets in a full model of guitar strings.⁴ Modal schemes in the context of collision dynamics have been successfully presented in other works: an implicit modal update was used in Ref. 13, where the collision force is resolved at each time step by employing a spatial grid, thus effectively employing a finite difference formulation; the special case α = 1 was given in Ref. 19. Here, it is shown that N-IT yields an efficient modal resolution that may be applied directly for all values of the barrier exponent α.

This scheme is a generalisation of Eq. (36). In order to account for a potential density, one may think of the barrier as being composed of N_b discrete points. For a continuous barrier (such as the backboard of a fretless instrument), one may assume that the barrier points are located at the string's grid points, in which case $Nb=N−1$⁠. For other kinds of barrier, one may need to specify points in between the string's grid points (e.g., for fretted instruments). In either case, one may approximate $ψ(x,t)$ by a vector $ψn−1/2∈ℝNb$⁠. The density distribution can be thought of as an $(N−1)×Nb$ sparse matrix $R$⁠. In practice,

where $ri∈ℝN−1$ is the interpolated sparse density vector of the ith barrier point. It is convenient, formally, to introduce a diagonal matrix containing all the explicit gradients $gin,i∈[1,Nb]$⁠. Hence, $G=diag(gin)$⁠. Then, define $T=RG$⁠.

With this notation, a finite difference scheme discretising Eq. (41) is

The discrete energy in this case is

which clearly discretises Eq. (40).

Proceeding in a similar manner as before, one may express this system as

where $a=(h/4)T⊺un−1+ψn−1/2, b=2un−un−1+T0k2/ρD(2)un$⁠. The update matrix in Eq. (44) is a rank-N_b perturbation of the identity matrix, though it is typically a perturbation of much smaller rank, equal to the number of points colliding at the time step n for which $gi≠0$ (in general, just a fraction of the total points.) Efficient inversion strategies such as the Woodbury identity³⁹ may be employed here. Note that, when the barrier points are collocated at the grid locations, the update matrix is in fact diagonal, and the scheme is fully explicit. Once $un+1$ is known, one may update $ψ$ using Eq. (43b). Note that, when the barrier points are collocated at the finite difference grid locations, the update matrix is in fact diagonal, and the scheme is fully explicit.

A suitable modal expansion for the string's displacement under fixed conditions is given by

where $m∈[1,Nm]$ and N_m is the total number of modes. In order to account for an appropriate density of barrier points, the vector d contains the spatial distributions, such as, e.g., delta functions. Thus,

where N_b as before is the total number of barrier points. Thus, $d$ is of length N_b. Then, modal projection is performed by means of the L₂ inner product denoted here for two square integrable functions f, g as

Thus, the projected modal equations for Eq. (41) become

Owing to modal orthogonality, the matrix $⟨X,X⊺⟩$ is the identity matrix times the norm of the modes, which in this case is the same for all the modes, i.e., $||Xm||2=L/2∀m$⁠. Similarly, $⟨X,(X″)⊺⟩=−L/2Λ2$⁠, where $Λ$ is a diagonal matrix with diagonal elements $Λm,m=mπ/L$⁠. The matrix $⟨X,d⊺⟩=R$ is a $Nm×Nb$ a dense matrix containing the projections of all the modes at each barrier point, columnwise. As before, define $T=RG$⁠. Using a finite difference approximation on the time operators, one gets the following modal system:

with modal discrete energy given by

Inspection of the energy yields a stability condition, which in this case reads

The modal update is

where here $a=(1/4)T⊺qn−1+ψn−1/2, b=2qn−qn−1−(k2T0/ρ)Λ2qn$⁠. Hence, the modal update has the same form as Eq. (44), except now the matrix $T$ is dense.

Snapshots of the numerical outcome of the schemes are presented in Figs. 5 and 6. In both cases, the modal scheme and the finite difference scheme return consistent solutions, even after a large number of collisions. Notice that, although the barrier parameters are selected so to simulate a hard collision, the schemes are perfectly stable and compute the solution solving one single linear system per update.

In Fig. 5, a bent obstacle similar to the bridge of instruments such as the tanpura is obtained as a quadratic function of x. In Fig. 6, frets are placed at intermediate string grid points using linear interpolation.

Power-law contact forces may be applied to the case of collisions between moving distributed objects, see, e.g., Ref. 40. An interesting case is represented by the wire-membrane interaction in the snare drum. In this instrument, a set of wires collide against the snare membrane. The system is activated after the batter membrane, at the opposite end of the drum, is set into motion by the player. This produces a vibration of the air cavity which in turn sets the snare membrane into vibration. Computationally, this is a complex system since various subsystems of different wave speeds are coupled in a nonlinear manner. A full model of this system was first offered by Bilbao,⁹ where the wire-membrane collisions are modelled in a semi-conservative manner. Bilbao et al.¹⁷ subsequently used the iterative model to obtain conservation of the discrete energy to machine accuracy. Here, as a test case, only the wire-membrane interaction is shown, using a non-iterative finite difference scheme.

In this basic configuration, the wire-membrane system may be described by two coupled partial differential equations,

Implementation details for the system in Eq. (51) are not given here. In terms of the membrane, one may use a 2D cartesian grid over which the differential operators are discretised. This leads to a staircase representation of the circular boundary, but it has overall numerous beneficial effects in terms of numerical dispersion and ease of implementation compared to, e.g., a polar difference grid, see, e.g., Bilbao²⁴ (Chap. 11). The collision force may be resolved once a suitable interpolation is implemented, to switch between the membrane and string grids. This can be done with a suitable 2D linear interpolator, as the one given by Bilbao²⁴ (Chap. 10). A non-iterative scheme for the collision force can then be implemented easily from the templates given above.

As an illustration, consider Fig. 7. Here, the wire is initialised in its first mode of vibration and released. The collision with the membrane induces a set of waves propagating in the membrane. The system can be updated at each time step by solving one single sparse linear system, as a perturbation of rank lower than the wire's grid points.

The Tromba Marina is a medieval bowed monochord instrument that produces a trumpet-like sound when played (hence the name tromba, meaning trumpet in Italian) (see Fig. 8). This characteristic sound arises from the fact that the string rests on a loose shoe-shaped bridge that collides with the instrument body as the string vibrates (see Fig. 9). This makes it a suitable test case for the method proposed here.

A working simulation of this instrument where the interactions between different components were based on Eq. (3) was previously published by the same authors.²⁹ Below, the details of the simulation using Eq. (27) are shown.

The complete instrument is subdivided into three components: the bowed string, the bridge, and the body. These are modelled as a stiff string, mass, and plate, respectively, all with loss terms.

Consider a damped stiff string of length L, described by displacement $u=u(χ,t)$⁠, with $χ∈[0,L]$⁠. Assume a linear differential operator of the form

with linear density ρ_s, cross-sectional area $A=πr2$⁠, radius r, tension T_s, Young's modulus E, area moment of inertia $I=πr4/4$⁠, and loss coefficients $σs1$ and $σs1$⁠. The equation of motion for the bowed string can then be given as

Here, Dirac delta function δ locates the bowing force at externally supplied bowing position $χb=χb(t)$⁠, and $Fb=Fb(t)$ is the externally supplied bowing force. Finally, $Φ$ is the dimensionless friction characteristic described in Ref. 29, with relative velocity (between the bow and the string at the bowing location) $vrel=∂tu(χb,t)−vb$⁠, and externally supplied bowing force $vb=vb(t)$⁠.

Similar to Eq. (1), the bridge is modelled as a simple point-like mass. Its displacement is $w=w(t)$ and its differential operator is

with mass M, stiffness K, and loss coefficient σ_m.

The body is here modelled as a 2D plate whose flexural displacement is $z=z(x,y,t)$⁠, where $(x,y)∈[0,Lx]×[0,Ly]$ and where L_x and L_y are the side lengths. Thus,

with surface density ρ_p, stiffness coefficient D, and loss coefficients $σp0$ and $σp1$⁠.

Interactions between the components are modelled using N-IT in two different configurations.

The interaction between the bridge and the body is modelled using Eq. (2), such that

is the difference between the state of the body at input location (x_i, y_i) and the bridge. The input spatial distribution is here assumed to be a 2D Dirac's delta $δ(x−xi,y−yi)$⁠.

The interaction between the string and the bridge is using a two-sided version of Eq. (2), allowing the collision potential to act as a connection, and is modelled as²⁸ 

depending on the difference between the state of the bridge and the string at contact location $ζ(t)=w(t)−u(χm,t)$⁠. Here, $Kζ$ is a constant, and $β≥1$⁠.

The effects of the interactions can be added to the respective components to yield the complete system

For most implementation details such as a discrete form of the equations in Eq. (58) and parameter values, one may refer to previous work.²⁹ The improvements will be discussed here.

The main improvement is the use of Eq. (27) as a discrete form of Eq. (2). The change to this new form of gⁿ, as noted above, solves the problem of spurious oscillations experienced at values of α and β larger than 1. These are now changed to be 1.3 and the interactions between components are now nonlinear.

The full system has been implemented in real time using C++ and the JUCE framework.⁴¹ A video showcasing the implementation can be found via.⁴² Snapshots of the bowed system can be seen in Fig. 10. As the string is bowed, it causes the bridge to collide with the body.