Geometric-integration tools for the simulation of musical sounds

During the last decade, much attention has been given to sound rendering and the simulation of acoustic phenomena by solving appropriate models described by Hamiltonian partial differential equations. In this contribution, we introduce a procedure to develop appropriate tools inspired from geometric integration in order to simulate musical sounds. Geometric integrators are numerical integrators of excellent quality that are designed exclusively for Hamiltonian ordinary differential equations. The introduced procedure is a combination of two techniques in geometric integration: the semi-discretization method by Celledoni et al. (J Comput Phys 231:6770–6789, 2012) and symplectic partitioned Runge–Kutta methods. This combination turns out to be a right procedure that derives numerical schemes that are effective and suitable for computation of musical sounds. By using this procedure we derive a series of explicit integration algorithms for a simple model describing piano sounds as a representative example for virtual instruments. We demonstrate the advantage of the numerical methods by evaluating a variety of numerical test cases.


Introduction
We introduce a systematic procedure inspired from geometric integration in order to simulate musical sounds. In this context, we consider a piano model as a representative example for virtual instruments. The piano sounds are computed by integrating a Hamiltonian partial differential equation (PDE) model describing the oscillations of the string and an ordinary differential equation (ODE) model describing the dynamics of the hammer. The procedure has its basis on the semi-discretization method by Celledoni et al. [13], which semi-discretizes the PDE to a system of ODEs while preserving the Hamiltonian structure. This method allows the application of geometric integrators, which are numerical integrators with excellent quality for Hamiltonian ODEs. Symplectic partitioned Runge-Kutta methods are particularly focused on because they yield explicit numerical schemes for a certain class of Hamiltonian systems. The Hamiltonian systems in the class are called separable Hamiltonian systems, whose Hamiltonian is nicely divided into a sum of two functions. Fortunately, most DE models in sound synthesis are of this type. One example of such models is a model of bar vibrations where u is the transverse displacement of the bar, ψ is the rotation of the bar crosssection relative to the normal, v and φ are the corresponding velocity to u and ψ and γ l , γ s , ε ∈ R. Another example is the Webster equation [30] which is a model of sound waves in vocal tracts or bodies of wind instruments where the pressure in the tube are denoted with p, the volume velocity in it with u, the function of x describing the cross-section area of the tube with S and γ ∈ R.
Other examples are introduced, for example, in [9]. In this contribution, we illustrate that the combination of the above two techniques in geometric integration is a right procedure for designing numerical schemes for computation of sound waves, in that the procedure indeed facilitates the design of stable numerical schemes for the simulation of musical instruments. First, we briefly summarize the recent developments in the field of musical sound synthesis as well as the difficulties, and illustrate their connections to geometric integration.
In the past decade, large efforts have been devoted to the simulation of acoustic effects and sounds. In the context of special effects or more general in computergenerated movies, this is simply motivated by the fact that traditional computational physics simulations usually lead to silent movies, because no practical algorithms existed for synthesizing synchronized sounds automatically. Instead, sound recordings were edited manually during the animation process or triggered automatically in interactive applications. Since the former is inflexible and labor intensive and the latter one produces dreary and repetitive results, researchers have investigated on this; see e.g. [32,34]. Furthermore, the simulation of sounds is well motivated due to the interest in the development of virtual instruments. Such digital devices would be superior to the conventional real musical instruments. For example according to [52] they would be less expensive because different instruments would be able to share a common input device; e.g. a virtual flute would be able to produce sounds of any kind of wind instruments. This makes it affordable for a variety of people enriching their creative work. Also, tuning and any other kind of labor intensive maintaining would not be necessary and the transportation of large and sensitive instruments can be avoided and location-based constraints therefore easily resolved-people from different places can join a common virtual orchestra.
The conventional approach to sound synthesis of musical instruments is based on signal processing-related techniques (e.g. [1,2,52,55]). This is currently an established way of musical sound synthesis because the produced sounds are fairly well perceptually and the algorithms are computationally efficient, so that digital interactive sound systems working in real-time can be developed. Although this approach has achieved a great success, it comes with significant shortcomings: the models have no definite physical interpretation and the quality of sounds is often less than satisfactory. In particular, the unpredictable sounds produced by the non-linear interaction between the input devices (e.g. the hammers in the case of a piano) and the instruments (e.g. the strings and the bodies of the piano) are not successfully reproduced. These difficulties can be resolved using sound synthesis based on appropriate physical models of virtual musical instruments. One of the most significant approaches is the one where the motion of the fundamental components of the instruments is described by differential equations (e.g. strings, hammers, and bows). Compared to the conventional approaches inspired by signal processing, the parameters in these models directly represent physical features of the instruments (e.g. the material of the body). Appropriate fitting parameters can be integrated, which enables the design of more realistic models. Previous research in this direction includes the modeling and simulation of the hammer [11,12,19,53], the key action [31,44,45,48,49], string vibrations [3][4][5]17,18,54], and the soundboard [20,21,36]. The interactions between the components are also considered in the literature; in particular, Chabassier et al. established a model and a numerical method for simulation of the whole piano [14,16]. Illustration of an example of computed sound waves for 0.01 s generated from a string of a piano (left) and the motion of the whole string for the same period (right). There are approximately three periods observed. Because the motion of the string for each period is almost identical, typically, just one motion period is computed The research on sound synthesis using PDEs is still at the beginning and there are many open problems. In this contribution, we address the challenging development of efficient and high-quality numerical integration algorithms for such models. This is a difficult task, because the auditory area of human beings is approximately between 20 and 20,000 Hz, for which reason the computation of thousands of vibrations is required for a sound of just a few seconds. Assume, that this sound wave is generated from a string with both ends fixed on the body of an instrument so that the vibration of the string is modeled by the wave equation under Dirichlet boundary conditions where u is the amplitude of the vibration and c is the speed of the wave. Then each peak of the sound wave corresponds to one periodic motion of the wave packet; see Fig. 1. In typical numerical computations of the wave equation, just one period of the periodic motion is of interest because the behavior of the waves is almost the same in each of the repetitions. However, in the simulation of sound we need to compute thousands of periods of the motion. In other words, the simulation of sound waves requires a long-time calculation compared to the time scale of the phenomena. In those cases, numerical methods must be carefully designed because conventional ones usually result in unstable or meaningless solutions.
In this context, one of the most successful approaches is the energy-based one by Bilbao et al. (e.g. [7][8][9][10][11]15]), which made a breakthrough in this area. Their numerical schemes are designed in a way, that a discrete approximation of the total energy of the system is exactly preserved. Because the energy often dominates the norm of the solution, the preservation of the energy results in a bound of the numerical solution. Hence this way of construction makes the resulting schemes stable and long-time calculations possible.
The aim of this contribution is to introduce a procedure that automatically derives numerical schemes with such a property. The key tools are from geometric integration, which is briefly explained below.
Long-time computations are also required in other research areas such as electromagnetics, quantum theory, fluid-, electro-and molecular dynamics, plasma transport, and celestial mechanics. In such areas so-called geometric integrators are employed to solve the occurring ordinary, partial, or stochastic differential equations derived from Hamiltonian mechanics. These methods typically discretize the underlying equations while preserving the mechanical and/or the geometric structure of the differential equations. As an example, the discrete gradient method is a method to derive energyconservative and energy-dissipative numerical schemes for the Hamilton equation and the gradient flows, respectively (e.g. [27,28,43,46]). A similar method for PDEs also exists, which is called the discrete variational derivative method (e.g. [22][23][24]26]). Other examples are symplectic integrators, which are numerical methods that preserve the symplecticity of the Hamiltonian flow in the discrete setup. The application of a backward error analysis shows that numerical solutions of these methods are the same as solutions of the Hamilton equation which is an approximation of the original equation [47]. As a consequence, energy conservation laws and other similar conservation laws (e.g. the conservation of the linear and the angular momentum) are approximately preserved by these methods, which leads to a globally accurate behavior. Because of these conservation laws, such algorithms often outclass conventional numerical methods in stability and reproducibility of significant phenomena.
In this regard, the goal of our work is the development of efficient geometric integrators for the models for musical instruments. The key observation is that most PDE models for musical instruments are separable Hamiltonian systems. Therefore, as explained in the first paragraph of this section, symplectic partitioned Runge-Kutta methods give explicit and hence efficient integrators for these systems. In order to apply symplectic partitioned Runge-Kutta methods, the models must be semi-discretized to ODEs while preserving the separable Hamiltonian structure. To achieve this, we focus our attention on the semi-discretization method, which we call the variational semi-discretization, by Celledoni et al. [13]. The variational semi-discretization is originally proposed as a method for deriving a suitable semi-discrete scheme for designing numerical schemes that preserve a certain energy behavior. However, as suggested in [13], this method could be used also for deriving semi-discrete schemes for Hamiltonian systems while preserving the Hamiltonian structure. The procedure introduced in this contribution is a combination of this semi-discretization method and symplectic partitioned Runge-Kutta methods. This procedure automatically derives explicit and symplectic integrators for most models for musical instruments. In this contribution, we illustrate this procedure by applying it to a simple model of the piano to develop symplectic numerical methods.
Remark 1 A similar, but slightly different semi-discretization is obtained by the discrete variational derivative method (DVDM) [22][23][24][25][26][37][38][39][40][41][42]. The DVDM derives energy-preserving or -dissipative numerical schemes for a certain class of PDEs. Taking the limit of the scheme for the Hamilton PDEs by the DVDM as the time step size goes to 0 yields the semi-discretized Hamilton ODEs in principle. The difference between these two approaches is the treatment of the boundary conditions. In the variational semi-discretization, the boundary conditions are included in the definition of the discrete phase space, and hence semi-discretized schemes by this approach are always Hamiltonian. On the other hand, in the DVDM appropriate discrete boundary conditions are assumed; that is, discrete boundary conditions that are compatible with the method must be imposed. This paper is organized as follows. In Sect. 2 the model of the piano is described. In Sect. 3 we explain the variational semi-discretization, which is the technique to derive a semi-discrete scheme while preserving the Hamiltonian structure of the equation. We apply this approach to the piano model for illustration reasons. After that, we develop several symplectic numerical methods by applying symplectic Runge-Kutta methods in Sect. 4.

Mathematical model for virtual pianos
Pianos are composed of many distinct parts, such as strings, hammers, black and white keys, and a sounding board. Although an excellent model that consists of most of these parts was recently proposed in [14,16], we use a rather simplified model, which only consists of a string part and a hammer part. This is because the aim of this paper is not the development of a realistic piano model but the introduction of a way to automatically get a simulation method that comprises an arbitrary geometric integrator.
The model we use in this study consists of a wave-type equation for the motion of a string, and the mass-spring model for the hammer's motion. These models are typical for piano simulations. The string displacement and the velocity are denoted with u(t, x) and v(t, x) respectively, the wave speed with c, the stiffness with κ, the frequency independent damping coefficient with d 1 and the frequency dependent damping coefficient with d 3 . Similarly, the displacement, the velocity and the mass of the hammer are denoted with u h , v h and M h . We assume that all the coefficients are positive. Since the ends of the string are fixed to the piano body, we assume the boundary conditions More practical boundary conditions are provided e.g. in [56]. The non-linear interaction between the hammer and the string is specified using the function where we denote the felt stiffness coefficient with K h , the hammer stiffness exponent with α and the felt loss coefficient with μ. These coefficients also have a positive value.
Here, ε(x) determines the hammer-string collision profile and satisfies ε, 1 = 1 so that ε(x) becomes an approximation to the delta function. The expression [ξ ] + is defined by and ε, u denotes the L 2 inner product The estimated values of the parameters for the C4 tone are listed in Tables 1 and 2.
As it is shown below, when the damping terms are ignored, that is, μ = d 1 = d 3 = 0, the above model is a separable Hamiltonian system, which is a system with a remarkable Hamiltonian structure from a viewpoint of numerical analysis. As explained in Sect. 4, this special Hamiltonian structure allows us to design explicit symplectic numerical methods.
To illustrate the Hamiltonian structure of this system, we introduce p h := M h v h , and an energy function This is a separable Hamiltonian in the sense that H is written as a sum of two functions: Then it is straightforward to check the following theorem.
Moreover, if the damping coefficients μ, d 1 and d 3 vanish, the system is a separable Hamiltonian system.
We show the energy behavior of this model in the following theorem.

Theorem 2 Under the boundary conditions (3), the energy function H is not increasing:
Proof From the integration by parts and the boundary conditions Substitution leads to

Variational semi-discretization and the application to the piano model
In order to make use of geometric integrators, a spatial discretization that preserves the Hamiltonian structure is required. The semi-discretization is typically done by replacing the spatial differential operators in the target equation by spatial difference operators, or using the finite element method. However, the resulting scheme does not always admit the Hamiltonian structure by using such approaches. To avoid the absence of the Hamiltonian structure, we use the variational semi-discretization [13] where the Hamiltonian is first discretized and then a semi-discrete scheme is obtained by variational calculus (see Fig. 2). This process automatically leads to a semi-discrete Hamiltonian scheme.
In what follows, for illustration purpose this procedure is applied to the Hamilton equation describing the piano model without the damping terms. We use a uniform mesh with a step size Δx = L/N so that N + 1 equals to a total number of points in the interval [0, L], and denote the approximated value of p(t, lΔx) with p l (t), or p l by omitting the argument t. We also denote a forward, a backward and a second difference operator with Fig. 2 This diagram shows the core idea of the variational semi-discretization. In typical approaches (shown by the dashed lines) the equation is directly discretized and hence the Hamiltonian structure may not be preserved. Instead, in this approach (shown by the solid lines) the Hamiltonian is discretized first and then the semi-discrete scheme is obtained as the Hamilton equation. This procedure automatically leads to a Hamiltonian semi-discrete scheme respectively. We first consider the phase space: Forũ,ṽ ∈X , we define the discrete HamiltonianH d inX using the trapezoidal rule to approximate the integral in H : where ε ∈X is an approximation of the function ε(x) and ·, · X means ũ,ṽ X =ũ 0ṽ0 Δx forũ,ṽ ∈X .

Remark 2
Regarding the approximation to the term u x in H , we chose which is a typical choice in [26]; however, we can also use, for example, the central difference and define the Hamiltonian as As is pointed out in pp. 90-91 of [26], it is difficult to judge whether a choice of the approximations of H defines a useful scheme or not, because it depends on the equation and possibly on other factors. Hereinafter we mainly useH d as the discrete Hamiltonian because it is found from the numerical tests, which are shown in Fig. 4 in Sect. 4, that the numerical solutions derived by usingH d converge to the exact ones slower than that derived by usingH d .
The boundary conditions corresponding to (3) are imposed by and by using them, we can rewriteH d without the boundary terms u −1 , u 0 , u N and u N +1 . This is equivalent to the restriction ofH d to the subspace ofX : We denote this restricted discrete Hamiltonian with H d . For u, v ∈ X , H d is defined as where ·, · X means We note that this is equivalent to (7) under the boundary condition (8).
We calculate the partial derivatives of H d to obtain the gradient in the Hamilton equation. The partial derivatives of H d with respect to u l and v l (l = 1, . . . , N − 1) are and the partial derivatives with respect to u h and p h are We denote the difference matrices of order N − 1 with The semi-discrete scheme is defined by the following separable Hamiltonian system: where ∇ q H d means the gradient of H d in the q direction associated with the inner product (10). In the following theorem we address the energy behavior of this model.
Remark 3 This theorem is generalized to include the damping terms in Theorem 5, and a proof is given there.

Application of symplectic integrators
We apply symplectic partitioned Runge-Kutta (PRK) methods to (11) for illustration purpose. PRK methods are a series of numerical methods for equations of the form Definition 1 Let a i j , b i , andâ i j ,b i be the coefficients of two Runge-Kutta methods. Then, an s-stage PRK method for (12) with a step size Δt is given by As mentioned above, long-time computations are required for the simulation of musical sounds. For this reason, in addition to accuracy, we need to take long-term stability and computational efficiency into consideration. All these three requirements are fulfilled by the application of a special class of PRK methods. As explained before, if a method is symplectic, the method has superior long-term stability in most cases. The following theorem identifies the condition for PRK methods to be symplectic; see [33,50,51].

The first condition is not necessary if the system is a separable Hamiltonian system.
Although the PRK method is available for any semi-discretized equation, the SPRK method requires the Hamiltonian structure for the equations. The semi-discretization in the way of Sect. 3 allows us to apply the SPRK method to any Hamiltonian PDEs.
In addition, if carefully designed, symplectic and explicit methods can be designed for separable Hamiltonians; see [29]. These methods are not only being explicit. Moreover, they have a favourable implementation where no additional storage is necessary. We use the coefficients of Table 3 to achieve the above properties. A PRK method with Table 3 The Butcher tableau of an s-stage SPRK method ···b s these coefficients becomes symplectic and explicit for separable Hamiltonian systems. Following [50], we write the coefficients of Table 3 as Furthermore, these coefficients enable us to use Algorithm 1 to reduce the amount of storage. In Algorithm 1, for example, Q 0 can be overwritten on q n and this applies also for other Q i 's. In summary, when using model (11) and a PRK method with the coefficients shown in Table 3, the method becomes explicit and symplectic because of the separability of the Hamiltonian system. Furthermore using Algorithm 1 the method is implemented with small amount of storage.
We tested three numerical schemes for the computation of piano sounds. These schemes are obtained by applying the SPRK methods to d dt where δ x . These are the semi-discretized Eq. (11) along with the damping terms added for realistic behavior. The energy behavior of this system is described by the following theorem.

Theorem 5 The discrete energy function H d (9) is not increasing:
Δx Moreover, if μ = d 1 = d 3 = 0, the system is conservative according to Theorem 3.
Proof From the differentiation of H d , we obtain Instead of the integration by parts, we introduce the summation by parts [22]: Using the summation by parts, we obtain We used the boundary conditions (8) at the last equality. Substituting the semidiscretized scheme (15), we get In the following numerical tests, we confirm that sounds generated by the numerical schemes derived by the above procedure certainly have basic characteristics of piano tones. We tested the 3-stage 3rd-order, the 4-stage 4th-order, and the 6-stage 4th-order SPRK methods with the following coefficients: where ω = (2 + 2 1/3 + 2 −1/3 )/3 and ν = 1 − 2ω; see [51]. The initial conditions are We employed the approximated delta function ε(l h ) used in [9]: (otherwise) According to [9], this is a 3rd-order approximation of δ(x − (β h + l h )Δx). Because in (11) the central difference operators are used in the spatial direction, the entire scheme becomes 2nd-order in space. We used the values of u(t, 0.7L) as the computed sound waves and also set l h = 0.2N and β h = 0.3. Before testing the piano sound, we investigate the validity of the discretization method for H . We excluded the damping terms and the hammer, and only consider the Hamilton PDE that describes the string in this validation. Figure 3 We set two constants above to A m = B m = 1/15. We also used this exact solution at t = 0 as the initial condition. The graph in Fig. 3 shows that the numerical solution indeed converges to the exact solution as Δt → 0 and Δx → 0. This is in fact quantitatively confirmed in Table 4, where the L 2 -norm and the L ∞ -norm are defined by  We also show the numerical solutions of the scheme derived by usingH d in Fig. 4. Although we used the same or less Δt and Δx, the waveform of that byH d is closer to the exact solution. Figure 5 shows the gap between the computed value and the exact value of the discrete energy with A m = B m = 10 −5 and N = 100. The exact value of the discrete energy is approximately equal to 0.110. We used Δt = 1/ (44,100 · 200) and N = 1000 hereinafter if it is not specifically noted. Figure 6 shows the numerical solutions obtained by the 3-stage 3rd-order, the 4stage 4th-order, and the 6-stage 4th-order SPRK schemes for (14). Any significant difference is not observed between these figures. We also compare the notes calculated by each method by carefully listening to them; however we did not notice a difference again. Hence, concerning the computation time, we conclude that the 3-stage or the 4-stage method is practical enough. Figure 7 is the enlarged figure of the waveform of the 4-stage 4th-order SPRK method. We find that this waveform is formed by repeating several kinds of waves with different amplitudes one after the other. This result gives an expectation that this waveform is a superposition of the wave of 261.63 Hz, which is the frequency of C4, and integer multiples of it. To confirm this, we show the spectrum of the waveform in Fig. 8.
There are large peaks expectedly near the positive integer multiples of 261.63 Hz. The notes of a real piano are indeed a superposition of such frequency components. Actually the spectrum shown in Fig. 8 is similar to those reported in the literature; see [18]. Figure 9 shows the gap of the energy H d (9) between the value of the numerical solution by the 4-stage 4th-order SPRK method and the exact value, which is approximately equal to 4.5496 × 10 −2 . We excluded the damping terms in this numerical test so that the energy is preserved.
The graph on the top shows that the displacements are within a certain range despite a large amount of the calculation steps ( 44,100,000 steps for 10 s). This energy behav- ior is due to the symplectic property of the method and shows that the proposed scheme certainly has a superior property regarding stability. Moreover, a similar result to Fig. 3 in which the recovery of the energy conservation law as Δt → 0 is shown is again observed in this test. Figures 10, 11, 12 show the result when the number of points N is changed from 1000 to 50. We used l h = 0.2N and β h = 0.00945 so that the hammer strikes the same position (x ≈ 0.126 L) of the string as in the previous experiments. In the first two experiments, the damping terms are included. Compared to Fig. 6, the waveform in Fig. 10 is smoother, which implies suppression of high-frequency tones. By carefully listening to the calculated notes, we in fact noticed that the sound was slightly blurred; Fig. 9 The evolutions of the error of the energy of the numerical solutions for the first 10 s computed with the 4-stage 4th-order SPRK method with N = 1000 and Δt = 1/(44,100 · 10), 1/(44,100 · 100). A similar result to Fig. 3 which shows the convergence of the energy as Δt → 0 is observed on the other hand, as shown in Fig. 11, the power and the peak of the spectrum in the low-frequency zone are almost unchanged. The gap between the computed and the exact energy, which is approximately equal to 4.5496 × 10 −2 , is shown in Fig. 12. The values of H d are still within a certain fixed range and converge to the exact value by Δt → 0 as well as in the case illustrated in Fig. 9.

Conclusion
Recently, much attention has been paid to novel approaches to the development of virtual musical instruments, where the PDE models of the components of the instruments Fig. 10 The waveforms of the string computed with the 4-stage 4th-order SPRK method. We changed the number of points N from 1000 to 50. The shape of the outline of the waveform is smoother than that in Fig. 6 are solved numerically. Since extensively long-time calculations are required to reproduce notes even for a few seconds, the computation time is significantly large and the accumulation of errors is not negligible. Hence numerical schemes for the musical simulations must be carefully designed-not only accurate and stable, but also efficient.
In this contribution we have introduced a procedure for deriving numerical schemes for models of musical instruments. The procedure is a combination of the variational semi-discretization by Celledoni et al. and the symplectic Runge-Kutta methods. The outline of the variational semi-discretization is illustrated in Fig. 2. This technique automatically derives a semi-discrete scheme while preserving the Hamiltonian structure. Thereby, geometric integrators can be immediately applied without any additional steps. Geometric integrators are numerical integrators of ODEs that preserve a significant property of the equations, typically energy conservation or symplecticity. By preserving one of these properties, the exact or approximated energy is accurately conserved. Since with this discrete conservation law numerical schemes often have excellent stability properties, the above procedure facilitates the design of several stable numerical schemes for musical simulations. We focus our attention on the observation that most PDE models of musical instruments are separable Hamiltonian systems and also on the fact that a class of SPRK methods yields explicit schemes for this type of Hamiltonian systems. Based on these facts, we have shown that the combination of the variational semi-discretization and SPRK methods is a right procedure for deriving numerical schemes that are suitable for simulations of musical instruments; indeed this procedure automatically yields explicit and symplectic schemes of a high order of accuracy for most of the models for musical instruments.
For illustration purposes, we have applied this procedure to a simple piano model and have derived a series of symplectic integrators by the application of SPRK methods. In absence of the damping terms, the model is shown to be a separable Hamiltonian system, so that the schemes are explicit and computationally efficient for computing Fig. 11 The spectrum in the law-frequency zone of the numerical solution computed with the 4-stage 4th-order SPRK method with N = 50. The power and the peak of the spectrum in the low-frequency zone are almost unchanged piano sounds. We tested the 3-stage 3rd-order, the 4-stage 4th-order, and the 6-stage 4th-order PRK methods numerically and all of them are shown to be sufficiently stable. Although we used higher order schemes (in time), the 3-stage 3rd-order or the 4-stage 4th-order method may be practical enough; almost no difference is observed between the waveforms computed by these methods. In particular, the 6-stage method needs more computational time but the result is almost the same compared to the other methods used in our numerical experiments.
Since we only took the consideration of the accuracy in the time direction into account, and only used the 2nd order difference operators in the spatial direction, in our future work we plan to improve accuracy in the space direction. In particular, the use of higher order compact schemes, which are known to be effective in the calculations Fig. 12 The gap between the energy of the numerical solution computed with the 4-stage 4th-order SPRK method with N = 50 and that of the exact solution. The values are within a fixed range and converge to 0 as Δt tends to 0 of sound waves [35], is of importance. Also, this procedure must be tested for more realistic models of musical instruments. In this context, the model of a whole piano by Chabassier et al. (see [14,15]) is important, for which reason we plan to consider it in our future work.
From a theoretical perspective, the effectiveness of the application of symplectic integrators to dissipative systems should be investigated because the model for the piano has the damping terms. Although this is a challenging problem, there exist a few results on analyses on this topic (e.g. [6]). The results of these analyses could give an insight on the qualitative acoustical analyses of computations of musical sounds.