Nonsmooth contact dynamics for the numerical simulation of collisions in musical string instruments

Collisions in musical string instruments play a fundamental role in explaining the sound production in various instruments such as sitars, tanpuras and electric basses. Contacts occuring during the vibration provide a nonlinear eﬀect which shapes a speciﬁc tone due to energy transfers and enriches the hearing experience. As such, they must be carefully simulated for the purpose of physically-based sound synthesis. Most of the numerical methods presented in the literature rely on a compliant modeling of the contact force between the string and the obstacle. In this contribution, numerical methods from nonsmooth contact dynamics are used to integrate the problem in time. A Moreau-Jean time-stepping scheme is combined with an exact scheme for phases with no contact, thus controlling the numerical dispersion. Results for a two-point bridge mimicking a tanpura and an electric bass are presented, showing the ability of the method to deal eﬃciently with such problems while invoking, as compared to a compliant approach, less modelling parameters and a reduced computational burden.


I. INTRODUCTION
Collisions are of prime importance in musical acoustics for explaining the particular timbre of a number of instruments ranging from strings (a typical example being that of indian instruments such as sitar, tanpura and veena) to drums (e.g.snare drum) (Bilbao, 2012;Fletcher and Rossing, 1998;Raman, 1921).In these examples, the role of contacts is to alter the frequency content due to a nonlinear, nonsmooth interaction that generates high frequencies and contributes to enrich the hearing experience.This effect is particularly prominent in the case of the sitar, where a curved bridge contributes to significantly modify the frequency content of the string vibration, see e.g.(Bilbao et al., 2015;Mandal and Wahi, 2015;Siddiq, 2012;Vyasarayani et al., 2009), and in the case of the tanpura and its particular bridge (Chatziioannou and van Walstijn, 2015;Issanchou et al., 2017;Valette et al., 1991).Other examples that attracted interest in the recent years concern the case of string/frets interactions (for e.g.guitar or electric bass) (Bilbao a Also at: IMSIA, ENSTA ParisTech-CNRS-EDF-CEA, Université Paris Saclay, 828 Boulevard des Maréchaux, 91762 Palaiseau Cedex, France b) cyril.touze@ensta-paristech.fr; Corresponding author.and Torin, 2015; Issanchou et al., 2018;Trautmann and Rabenstein, 2004) and snare drum where metal wires are in contact with a vibrating membrane (Bilbao, 2012;Bilbao et al., 2015).
Earlier studies on vibrating strings with contacts derive a number of analytical results, mostly in the 1980s, see e.g.(Amerio, 1978;Cabannes, 1987;Citrini, 1991;Schatzman, 1980).To overcome the limitations of these approaches, where the string needs to be perfect without stiffness, thus discarding the dispersive effect which has been shown to be of prime importance in the case of the sitar and tanpura (Chatziioannou and van Walstijn, 2015;Issanchou et al., 2017;Siddiq, 2012), recent research efforts concentrate toward the development of efficient, robust and accurate numerical methods in order to simulate musical strings encountering an obstacle during their vibration.Most of the methods presented in the last years use a regularisation in order to treat numerically the contact force, see e.g. the energy-conserving schemes proposed by Bilbao et al. (Bilbao et al., 2015;Desvages and Bilbao, 2015;Ducceschi et al., 2016) and by van Walstijn et al. (Chatziioannou and van Walstijn, 2015;van Walstijn and Bridges, 2016;van Walstijn et al., 2016), the modal approach proposed in (Issanchou et al., 2017) or the approach followed in (Inácio et al., 2006) to model the interaction between a puja (exciting stick) and a Tibetan bowl.In all these studies the contact force is modeled using a power-law method with two parameters defining the stiffness of the repelling force, thus allowing to cover a wide range of contact laws, from soft collisions, see e.g.contact between felt and hammer in piano (Boutillon, 1988), between mallet and membrane in a kettledrum (Rhaouti et al., 1999) or between finger and fretboard (Bilbao and Torin, 2015), to hard contacts.In this article, we will name as compliant approach the methods using a regularisation to express the repelling force.In this framework, power-law are often used but other functions may also be selected.Also, the parameters defining the contact force may have a physical basis, see e.g.(Goldsmith, 2001) for some examples.On the other hand, ad hoc values may be used as numerical free parameters, so that a penalty approach is at hand.
On the other hand, a large body of research has been dedicated to the development of nonsmooth numerical methods in order to deal efficiently with numerical challenges posed by contact and friction forces.These methods rely on specific assumptions (e.g.no interpenetration is allowed between the contacting bodies) and use mathematical tools from the measure theory, differential inclusions and complementarity systems.The first developments have been pioneered by Jean and Moreau, see e.g.(Jean, 1999;Jean and Moreau, 1987), continued by numerous investigations (Doyen et al., 2011;Janin and Lamarque, 2001;Paoli and Schatzman, 2002), and are now summarized in reference books (Acary and Brogliato, 2008;Studer, 2009).Nonsmooth methods have been succesfully applied in a variety of contexts ranging from granular media (Renouf et al., 2004), geomaterials (Jean, 1995), multibody dynamics (Chen et al., 2013) to realistic simulations of hair motions and living systems (Acary et al., 2014;Bertails-Descoubes et al., 2011).In the field of vibration, the method has been applied to rotor/casing contacts in (Meingast et al., 2014) as well as to string vibrations in (Ahn, 2007), where the study was however limited to the case of a perfect string without stiffness and a frictionless contact.In the area of musical acoustics, the action of a grand piano has been recently simulated efficiently by using a nonsmooth approach (Thorin et al., 2017).
As remarked by a number of investigators, numerical integration for contact dynamics is generally timeconsuming due to the high-frequency content generated (Doyen et al., 2011;Issanchou et al., 2017), leading to consider very small time steps in order to achieve convergence.In this context, nonsmooth numerical methods, as avoiding the costly step of finding the zeros of a nonlinear function with a Newton-Raphson approach and using efficient numerical methods to solve a linear complementarity problem, should decrease the computational burden (Acary and Brogliato, 2008).On the modeling point of view, nonsmooth methods are particularly appealing when one only needs an efficient numerical method to repel a vibrating structure from an obstacle with a simple representation of the dissipation phenomena at contact.As such, nonsmooth methods describe the contact law with a single parameter: the coef-ficient of restitution, instead of the two parameters used in the power-law approach as done in e.g.(Bilbao et al., 2015;Chatziioannou and van Walstijn, 2015;Issanchou et al., 2017) in the conservative case, and additional parameters to model the dissipation following a Hunt and Crossley approach.Let us note that it is also difficult with the power-law approach to obtain large dissipation rates at contact.For this purpose, more complicated compliant models must be introduced including plasticity effects (Nguyen and Brogliato, 2014).Less modelling parameters could be seen as an advantage, but also as a drawback if one prefers to have more degrees of freedom in order to represent different contact laws with a large variety of stiffness and/or damping to account for a physical reality having perceptual effects.A typical example is that of the guitar where very different parameters are used for the string/fret and finger/fingerboard collisions (Bilbao and Torin, 2015).Note however that nonsmooth methods can also be conjugated with a material description of a soft structure including nonlinear stiffness and damping to account for these effects and tune them at ease, see e.g. the recent modeling of felt in piano action proposed in (Thorin et al., 2017).
In this contribution, a nonsmooth numerical approach based on a Moreau-Jean time-stepping scheme is adapted to the case of a string vibrating against a stiff obstacle with hard contacts, for the specific purpose of musical acoustics.The examples have been purposely chosen to restrict our discussion to hard contacts where the nonsmooth method should be an interesting alternative to compliant approach.Standard nonsmooth methods are known to produce numerical dispersion (Yoong et al., 2018), which is a specific issue in the field of musical acoustics.During non-contacting phases, an exact scheme, introduced in (Bilbao, 2009) and used for a compliant contact approach in (Issanchou et al., 2017;van Walstijn and Bridges, 2016), is considered in order to ensure a numerical integration method with controlled dispersion.The efficiency of the scheme is demonstrated on two examples which are compared either to a compliant approach or to experimental results.First, the case of a two-point bridge mimicking a tanpura is shown.Results are then extended to a fretted electric bass having at most 20 contact points.The accuracy of the two different modeling options (nonsmooth vs compliant approach) are discussed and the computational burdens of these methods are compared.

II. MODEL A. Vibrating string against a unilateral obstacle
We consider a stiff string of length L (m), mass per unit length µ (kg• m −1 ) and tension T (N).Its Young's modulus E (Pa) and moment of inertia I define its stiffness.The string vibrates against a unilateral obstacle, the profile of which is described by g(x) (see Fig. 1).Eq. ( 1) describes the displacement of the string along (Oz): where the subscript t (respectively x) refers to a partial derivative with respect to time (respectively space).The right-hand side f (x, t) refers to the contact force per unit length.
The string is simply supported at its endpoints, so that ∀t ∈ R + : u(0, t) = u(L, t) = u xx (0, t) = u xx (L, t) = 0.In order to achieve a fine representation of eigenfrequencies and damping parameters, we employ the following modal description of the string: where φ j (x) = 2 L sin jπx L and q j is the j th modal amplitude.
Inserting the modal expansion (2) into the equation of motion (1), multiplying by another mode shape and integrating over the string length, and then adding losses to the final modal equations, one obtains: where q = [q 1 , q 2 , ...q Nm ] T is the vector of unknown modal amplitudes.Ω and Υ are diagonal matrices such that Ω contains radian frequencies Ω jj = ω j = 2πν j , ν j being the j th eigenfrequency, and Υ contains damping coefficients Υ jj = σ j .These quantities follow a model presented in (Issanchou et al., 2017 , 2018; Paté  et al., 2014; Valette and Cuesta, 1993).We briefly recall the expression of eigenfrequencies and damping parameters in Appendix.The vector F contains modal forces and represents the projection of the contact force onto mode shapes, its entries read

B. Contact force
We choose to model the contact condition by a unilateral constraint in order to avoid the interpenetration of the solids.Defining η(x, t) = u(x, t) − g(x) as the distance between the string and the obstacle, this constraint is given by: η In order to satisfy Eq. ( 4), the reaction force introduced in Eq. ( 1) has to be positive and vanishes only if the constraint is not active (open contact).The complete model is the Signorini law which can be simply written with the following notation (Acary and Brogliato, 2008;Signorini, 1933): meaning that either η(x, t) = 0 and f (x, t) ≥ 0, or f (x, t) = 0 and η(x, t) ≥ 0. Second order dynamics with unilateral constraints involves velocity jumps.The usual setting is to consider that the velocity map t → u t (x, t) is a right continuous function, i.e. u + t (x, t) = u t (x, t), and of bounded variation.The notation u + t (x, t) (respectively u − t (x, t)) stands for the right limit (respectively the left limit) of Since the function is of bounded variation, these limits exist.
In order to define the solution, especially in a discrete (finite-dimensional) system, a condition on the velocity after an impact must be specified.To this purpose, we choose the Newton impact law (Acary, 2016): where the coefficient of restitution ∈ [0, 1] defines the string behaviour at impact instants.
In the sequel, the Signorini law is formulated at the velocity level, which allows one to explicitly control contact losses in taking into account the Signorini condition together with the impact law.Therefore, the contact law writes: otherwise.
(7) The viability Lemma of Moreau (Moreau, 1999) ensures that the condition at the velocity level (7) implies the condition at the position level (5) if the constraint ( 4) is satisfied at the initial time.
Since the system has some discontinuities in the velocity w = q, the acceleration ẇ is not everywhere defined in the classical sense.A differential measure dw is associated with the velocity w and plays the role of the acceleration.If we assume furthermore that w is of special bounded variation (i.e.w may be decomposed into a sum of an absolutely continuous function and a jump function), the differential measure can be decomposed with respect to the Lebesgue measure dt as: where dν is a discrete measure of the form i a i δ si with given sequences {a i } and {s i } of real numbers.The notation δ s refers to the Dirac measure supported at time s (see e.g.(Moreau, 1988) for details on differential measures).In other words, ẇ(t) = q(t) almost everywhere and at impact instants t * we have dw = (w Similarly, the reaction force in the modal space is defined by a vector measure: where the vector F corresponds to a modal continuous contact force and the vector P is a modal contact impulse corresponding to velocity jumps.In terms of differential measure, the modal equations read as: By substituting ( 8) and ( 9) in ( 10), we remark that Eq. ( 3) is satisfied dt-almost everywhere.At the instants of discontinuities, we obtain the impact equation: The contact condition at the velocity level is also reformulated in terms of measures as follows: where di is the reaction force in the physical space.The relation between η, η and q, q is given by the relations between the quantities in the physical space and the modal space.Let us introduce the column vector φ(x) containing the N m first string modes, defined by φ j = φ j (x), ∀j ∈ {1, ..., N m }.From ( 2), we get that u(x, t) = φ T (x)q(t) where q = [q 1 (t), ..., q Nm (t)] T and therefore: η(x, t) = φ T (x)q(t) − g(x) and η(x, t) = φ T (x) q(t).
(13) By duality, we also have: Altogether, the dynamics is given by the following measure differential complementarity problem: In this section, we present a discretisation of the problem presented in Section II which is able to han-dle nonsmooth contacts.The distinctive feature of the proposed scheme is to combine an exact method for the linear (non-contacting) part of the equations of motion with a Moreau-Jean time-stepping approach to handle impulses and velocity jumps.The resulting scheme thus prevents dispersion when the string freely vibrates, which represents an improvement of already existing nonsmooth time-stepping methods.
In the rest of the paper, the variable t n denotes the discrete time t n = n∆t, where ∆t is the time step.The spatial grid is defined by x i = i∆x, ∀i ∈ {0, ..., N }, where ∆x = L N is the spatial step.Due to selected boundary conditions, u(x 0 , t) = 0 and u(x N , t) = 0 ∀t ∈ R + , so that only the interior points of the grid are considered.
As prescribed in a Moreau-Jean scheme, Eq. ( 15b) is integrated over (t n , t n+1 ].One thus obtains: The collision term on the right-hand side of Eq. ( 16) is treated with the Moreau-Jean scheme, leading to the following discrete approximation: where T is a contact impulse over the time interval, which will be defined later through a complementarity condition.The Moreau-Jean scheme is generally used with a θ-method to approximate the stiffness and damping terms in Eq. ( 16), as prescribed for example in (Acary and Brogliato, 2008;Jean, 1999;Jean and Moreau, 1987;Moreau, 1999).However, for mechanical vibratory systems expressed in the modal basis and with application to musical acoustics where the problem of numerical dispersion is particularly stringent, the exact scheme used in (Bilbao, 2009;Issanchou et al., 2017) shall be considered in order to improve the discretisation of the stiffness and damping terms.Using the following notations for discrete approximations, w n i ≈ w + i (t n ) , and q n i ≈ q i (t n ), it reads: where the introduced terms γ i and σ * i are such that: with: Finally, one obtains the following time-stepping scheme from ( 16) for the update of the modal velocity vector using ( 18) and ( 19): In matrix form, it writes: where C and Č are diagonal matrices with entries: .
The update of the modal amplitudes is given by the discretisation of (15)(a) as: We introduce a matrix S containing the N m first string modes, defined by S ij = φ j (x i ), ∀(i, j) ∈ {1, ..., N − 1} × {1, ..., N m }.Defining the physical discrete displacement vector u = [u(x 1 , t), ..., u(x N −1 , t)] T , one easily obtains the following relationship: u = Sq (Issanchou et al., 2017).Denoting the displacement of the string at contact points by u c = S c q, its velocity by v c = S c w and introducing the velocity without contact w n+1 free as: where the Delassus' matrix W cc is computed as: The complementarity condition (15g) is discretised in a fully implicit way as prescribed for the Moreau-Jean scheme.For a contact point indexed by c i , one has to solve: Let us denote the index set of active contact points as c = {c i | η n ci ≤ 0}.The impulse p n+1 c is the solution of: which is a Linear Complementarity Problem (LCP).A LCP consists in finding the vectors z, λ such that: for some given matrix W and given vector a.Let us recall that the LCP has a unique solution for all a if the matrix W is a positive definite matrix (Cottle et al., 1992).It is straightforward to identify in (32) the vector a as a = S cw n+1 free + v n c and the matrix W = W cc .In our application, the Delassus matrix W is obviously symmetric positive definite and a unique solution is ensured.If the string vibrates against a single point obstacle, the Delassus matrix reduces to a positive scalar and the LCP can be solved directly as: In the case of a distributed obstacle, a LCP solver has to be employed.In this article, we employ the numerical solvers provided by the Siconos software (Acary et al., 2016) which is dedicated to the modelling and simulation of nonsmooth dynamical systems.In Siconos, several numerical algorithms are implemented to solve LCPs.Depending on the number of constraints and the required accuracy, pivoting or projected successive overrelaxation (PSOR) techniques may be used (Cottle et al., 1992).Simulations in the present paper are based on the use of a standard pivoting technique for solving LCP, known as Lemke's method (Cottle et al., 1992) to get a high precision solution at the machine accuracy.

IV. NUMERICAL RESULTS
A. One contact point: the case of tanpura In this section, the case of a two point bridge is considered, which consists in placing a point obstacle near a boundary.The string parameters are selected as in (Issanchou et al., 2017) for the purpose of comparison.Geometric and material relevant properties are recalled in Table I.The point obstacle is located at 6 mm from the boundary x = 0, in agreement with the range of positions mentioned in (Valette and Cuesta, 1993) to mimick the tanpura bridge with a two point approximation.
L (m) d (mm) T (N) µ (kg.m −1 ) 1.002 0.43 180.5 1.17 × 10 −3 The initial position of the string is a smoothed centered triangle with a maximal amplitude u 0,max = 1.8 mm and no velocity.The nonsmooth numerical procedure is compared to the power-law approach described in (Issanchou et al., 2017), where the contact force is modeled as with K = 10 13 and α = 1.5 for simulations.In order to be in line with the assumptions retained in (Issanchou et al., 2017), we select an equal number of modes N m and grid points N , specifically N m = N − 1, and set N = 1002 for the simulations.For completeness, the outcomes of the nonsmooth method are compared for both = 0 and = 1.
A convergence study has been led using the stringent criterion introduced in (Issanchou et al., 2017) to control the long-term behaviour.The criterion inspects relative errors of the numerical results on a 3 seconds simulation and imposes the error to be less than 10 −1 .Following that, the sampling frequency is selected as F s = 2 MHz for both compliant and nonsmooth cases, since it has been found that the convergence is obtained at that sampling rate for the two methods.
The linear modal characteristics of the string are also selected as in (Issanchou et al., 2017), eigenfrequencies and damping ratios resulting from an experimental identification on a real guitar string.The identification has been realised up to 7200 Hz, then models are employed to input linear characteristics into the simulation up to the desired maximal frequency.Models used for linear characteristics are recalled in Appendix, and selected values resulting from model fitting shown in (Issanchou et al., 2017) are given in Table II.Fig. 2 shows the string displacement at a point located at 1 cm from x = L. Associated sound files corresponding to the displacement resampled at 44.1 kHz are available as supplementary materials 1 .Its appears that selecting = 0 or = 1 has no significant macroscopic effect on the simulation result.Both values give a temporal displacement extremely similar to that obtained with the penalty approach.In all cases, a crenel shape appears which evolves in time due to effects of dispersion and damping.As analysed for the power-law approach in (Issanchou et al., 2017) where simulation outcomes are compared to experimental signals, we can thus conclude that the nonsmooth method also retrieves the essential features of the string vibration against a point obstacle, with the same global accuracy as that obtained with the penalty approach.
Spectrograms of numerical signals with the penalty approach and the nonsmooth method with = 1 are shown in Fig. 3.No significant difference is observed when = 0, therefore we do not show the spectrogram in this case.In both cases, there is no missing mode despite the centered initial condition, due to energy transfers induced by collisions.Moreover, both methods accurately recover the descending formant observed experimentally.Such formants constitute a distinctive feature of musical instruments such as tanpuras.
We now explore the local behaviour of the simulated string displacement at the obstacle position u c in Fig. 4, with the nonsmooth method, = 0 and = 1, as well as with the penalty approach.Zooms focus on the first two contact periods.Despite similar global waveshapes, the string has local specific and distinguishable behaviours, around the obstacle, when colliding with it.When = 0, as expected since no back velocity is enforced, the string sticks to the obstacle for a short period of time when colliding, and then leaves it when the wave propagates back (i.e.around 3.8 × 10 −3 s for the first contact time, then around 8.8×10 −3 s for the second contact).When = 1, the string bounces off the obstacle with a magnitude of about 10 −7 m and the signal takes the form of one time step oscillations until the string leaves the obstacle for a longer time around 3.8 × 10 −3 s.With the penalty approach, smoother oscillations appear with similar amplitudes.Interestingly, the envelope of oscillations of the nonsmooth signal closely fits oscillations of the compliant signal.The detailed behaviour of the string, computed with the nonsmooth method and different values of F s , is presented in Fig. 5. Increasing the sampling frequency does not significantly change the global string behaviour, however the local behaviour at the obstacle position is affected.Note also that a very small penetration occurs with the nonsmooth method, a numerical artifact linked to the time discretisation.This numerical penetration is reduced and small bounces can be observed when = 1, after the string reaches the obstacle.When = 0, the string sticks to the obstacle at the height given by the first discrete point where contact condition is activated, defining the depth of this numerical penetration.
Finally, nonsmooth and power-law approaches give very similar results in the case of a point obstacle, and no significant difference is observed between simulations with = 0 and = 1 in the nonsmooth case, except when focusing on the obstacle position.This can be related to the analytic solutions for a string with a single contact point where a sticking time is obtained as long as the wave has not come back, see e.g.(Cabannes, 1984 , 1987)  and (Doyen et al., 2011;Yoong et al., 2018) for a similar case with bars.
We now discuss time costs of the nonsmooth ( = 1) and power-law (α = 1.5, K = 10 13 ) approaches.Com- putation costs, excluding initialisation, saving data and energy computation, are presented in Table III.For the compliant method, they are presented either with a spatial grid which is fully included in the matrix S, meaning that its dimensions are (N − 1) 2 , or with a matrix S which only includes the obstacle position, meaning that its dimensions are 1 × N − 1.The former case strictly corresponds to the scheme presented in (Issanchou et al., 2017), while the latter case (modified compliant method) contains an improvement where the number of grid points N is not anymore related to the number of modes N m .This further step is described in (van Walstijn and Bridges, 2016), and is also implemented for the nonsmooth method.All computations presented in Table III were led with MATLAB on a single CPU with a clock at 2.4 GHz.
One can observe that taking into account only the contact point in the modified compliant method allows one to obtain a significant gain in computational times since a factor two is present as compared to the power-law method with a full S matrix.A more significant gain can be obtained due to the nonsmooth method: computation costs are 12 to 37 times smaller.This is due in particular to the reduced number of matrix vector products as well as a reduced time for solving the LCP (straightforward in the point obstacle case) compared to Newton-Raphson iterations.Note also that computation times has been found to be once again divided by a factor two by running the same code with the SICONOS software in a optimized C/C++ implementation (Acary et al., 2016).
In this section, we discussed the behaviour of a string vibrating against a point obstacle, mimicking the case of 6.338 6.339 6.34 6.341 6.342 6.343 6.344 a tanpura.We observed that the nonsmooth simulation in the case of a two point bridge was very close to the simulation led with the power-law method, both being consistent with experimental data presented in (Issanchou et al., 2017).Moreover, it seems that the global behaviour of the string does not rely on the coefficient of restitution, even though the local string behaviour at the contact point does.In the following, we investigate a case implying multiple contact points.

B. Multiple contact points: application to the electric bass
In this section, the complexity is increased by considering a string vibrating against the fretboard of an electric bass with up to 20 contact points.Electric basses are known for their modern playing techniques such as pop and slap where numerous contacts are intentionally provoked in the transient attack, giving a peculiar bright and percussive sound (Bacon , 2013;Issanchou et al., 2018)  Fig. 6 shows the guitar neck profile with its 20 frets together with the initial condition imposed to the string.A G-string of an electric bass is considered for simulations, its properties are given in Table IV.The string is plucked at 64 cm from the nut with a maximal initial amplitude u 0,max = 3.6 mm.As in the previous section and according to the study led in (Issanchou et al., 2018), measured values of the linear characteristics are selected to fit with the model up to about 3400 Hz, then models are employed with parameters reported in Table V.In order to perform comparisons with the power-law approach and with experimental results obtained in (Issanchou et al., 2018), 863 modes are retained for simulations.The convergence of data is determined according the criterion presented in (Issanchou et al., 2017 , 2018)    Temporal results are shown in Fig. 7, with a comparison of the nonsmooth approach with the experimental result already shown in (Issanchou et al., 2018), and then with the compliant method.A spectral comparison is provided due to the spectrograms of the displacements shown in Fig. 8. Related sound files corresponding to the displacement resampled at 44.1 kHz are available as supplementary materials 2 .The presence of multiple contacts makes the dynamics more complex but leads once again to similar results for the two limit values of the restitution coefficient, namely = 0 and = 1.Simulation results show that the temporal signals are almost coincident.Consequently only the case = 1 is presented.It thus appears that the dynamics from a macroscopic point of view, is still dominated by waves going back and forth in inter-fret intervals.
Comparing outcomes of the displacement obtained with the nonsmooth method with experimental results shows an excellent agreement, both in the global shape of the time series and in the details shown in zooms in Fig. 7.The occurrence of a substantial high frequency content from the very first collisions is clearly ascertained and retrieved, and the long term behaviour remains similar.A few discrepancies appear, however as explained in (Issanchou et al., 2018) they are probably mainly due to (i) uncertainties related to the experimental measurement of the neck profile and, to a lesser extent, to the measured properties of the string; (ii) the assumption that the fretboard is a rigid obstacle.Comparing the two numerical methods shows that they behave very similarly, with only slight discrepancies.
Focusing on spectrograms, a similar structure between experimental and numerical data can be observed.In particular, energy transfers appear in the transient attack between 0 and 0.09 s, this time window corresponding to the occurrence of contacts between the string and the fretboard.Once again, this collision period of time, together with specific reinforced spectral zones (e.g.around 4500 Hz), is particularly well retrieved by numerical simulations.From a perceptive point of view, associ-  Finally, the global dynamics of the string colliding with frets is well described by both compliant and nonsmooth approaches.In the studied configuration, both numerical methods produce very similar results at the measurement point considered, selected from experiments.Furthermore, the value of the restitution coefficient does not substantially affect simulations, meaning that, in the present case, the nonsmooth approach would be free of parameters to be adjusted.
The simulation of 0.1 s of signal (time interval recovering all contacts) at F s = 1 MHz lasts 5 minutes with Matlab while it lasts only 30 s with Siconos.Numerical costs are thus found to be divided by a factor 10: once again, an important gain in computation can be achieved with the nonsmooth method on Siconos.

V. CONCLUSION
A nonsmooth approach to the numerical simulation of musical strings colliding with an obstacle has been presented.The main feature consists in combining an exact scheme for the vibratory (linear) part and a Moreau-Jean scheme for the contact force, embedded in a modal description of the dynamics.This results in a computationally efficient scheme preventing numerical dispersion during free flight phases.The method has been tested on two cases involving either one contact point (a two point bridge mimicking a tanpura) or multiple contacts (the case of an electric bass).In both cases, it has been found that the nonsmooth approach is very accurate and compares well with experiments and simulations led with a compliant method.A significant gain in computational times is obtained with the proposed nonsmooth method.Interestingly, for the two test cases, simulation results are relatively insensitive to the choice of the value of the restitution coefficient, so that the model could finally be considered as free of parameters to describe collisions.This particular behaviour should however be limited to a few number of pointwise contacts where, from a macroscopic point of view, the behaviour during collisions is governed by waves going back and forth in the string.In contrast, the coefficient of restitution should have an effect when considering a large number of contact points or continuous contact regions (e.g. a fretless bass), where sticking may occur on a whole interval.Note also that the nonsmooth method has been tested in specific examples where a hard collision is at hand and no specific physical detail on the contact is needed to have efficient simulations accounting for the rich contact dynamics of the string.For other cases encountered in musical acoustics where the contact is softer, e.g.mallet/membrane or finger/string interactions, the nonsmooth approach may also be used, either with a refined continuous description of the deformable bodies, involving an increase of the computational burden to take into account the whole dynamics; or with a localised behaviour law to account for the felt as done for example in (Thorin et al., 2017).These points are thus left for further studies where a more complete comparison of nonsmooth and compliant methods could be given.

FIG. 1 .
FIG.1.Scheme of a string vibrating against a unilateral obstacle.

FIG. 2 .
FIG.2.Displacement of the string vibrating against a point obstacle near a boundary, taken at 1 cm from the extremity x = L. Comparison between the penalty approach, Fs = 2 MHz (green line) and the nonsmooth method, Fs = 2 MHz, = 0 (black dashed line) and = 1 (red line).
FIG. 4. Displacement of the string vibrating at the obstacle position.Comparison between the numerical simulation with the penalty approach, Fs = 2 MHz (green line) and the nonsmooth method with = 0 (black dashed line, black line in zooms) and = 1 (red line), Fs = 2 MHz.
FIG. 7. Displacement of the string vibrating against the neck of a bass guitar, taken at 9 mm from the extremity x = L. Comparison between the experimental signal (blue line) and the nonsmooth method, = 1, Fs = 4 MHz (red line), and between the nonsmooth method, = 1, Fs = 4 MHz (red line) and the compliant method, Fs = 8 MHz (green line).

TABLE I .
Electric guitar string properties.

TABLE II .
Model parameters for the electric guitar string. .

TABLE III .
Computation times with the compliant and nonsmooth methods (implemented in Matlab), in seconds, for the simulation of a one second signal in the case of a point obstacle.
applied over a period of time recovering all string/neck collisions.This implies large values of F s , nevertheless smaller values may be selected for other applications.For instance, a convergence study on the spectral content, associated to a perception evaluation of signals with different F s , may lead to smaller values of F s with a satisfactory sound rendition.This is however not the purpose of the present paper.

TABLE IV .
Electric bass string properties.

TABLE V .
Models parameters for the electric bass string installed on an electric bass.