1 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

$$\begin{aligned}&u_t=v, \quad v_{t}=\gamma _s^2(u_{xx}-\psi _x), \\&\psi _t=\epsilon ^2\phi , \quad \epsilon ^2\phi _{t}=\epsilon ^2\gamma _l^2\psi _{xx}+\gamma _s^2(u_x-\psi ), \quad (x\in [0,1]) \\&H(u,\psi ,v,\phi )=\int _0^1 \left[ \frac{1}{2}v^2+\frac{\varepsilon ^2}{2}\phi ^2+\frac{\varepsilon ^2\gamma _l^2}{2}\psi _x^2+\frac{\gamma _s^2}{2}(u_x-\psi )^2\right] \mathrm {d}x, \end{aligned}$$

where u is the transverse displacement of the bar, \(\psi \) is the rotation of the bar cross-section relative to the normal, v and \(\phi \) are the corresponding velocity to u and \(\psi \) and \(\gamma _l\), \(\gamma _s\), \(\varepsilon \in {\mathbb {R}}\). Another example is the Webster equation [30] which is a model of sound waves in vocal tracts or bodies of wind instruments

$$\begin{aligned}&Sp_t=\gamma u_x, \quad \frac{1}{S}\,u_t=-\gamma p_x, \quad (x\in [0,1])\\&H(p,u)=\frac{\gamma }{2} \int _0^1 \left[ Sp^2+\frac{1}{S}\,u^2 \right] \mathrm {d}x, \end{aligned}$$

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 \(\gamma \in {\mathbb {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 computer-generated 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].

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

$$\begin{aligned}&u_{tt} - c^2 u_{xx} = 0, \quad x \in (0, L)\\&u(t, 0) = u(t, L) = 0, \end{aligned}$$

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.

Fig. 1
figure 1

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

Full size image

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 energy-conservative 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.

2 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

$$\begin{aligned}&{\left\{ \begin{array}{ll} u_t = v, \\ v_{t} = c^2 u_{xx} - \kappa ^2 u_{xxxx} - d_1 u_t + d_3 u_{txx} + \varepsilon f_\mathrm {h}, \end{array}\right. } \end{aligned}$$
(1)

for the motion of a string, and the mass-spring model

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle \frac{\mathrm {d}u_\mathrm {h}}{\mathrm {d}t} = v_\mathrm {h}, \\ \\ \displaystyle \frac{\mathrm {d}v_\mathrm {h}}{\mathrm {d}t} = -\frac{1}{M_\mathrm {h}}f_\mathrm {h}({u,u_\mathrm {h}, v,v_\mathrm {h}}), \end{array}\right. } \end{aligned}$$
(2)

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 \(\kappa \), 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_\mathrm {h}\), \(v_\mathrm {h}\) and \(M_\mathrm {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

$$\begin{aligned} u(t, 0) = u(t, L) = u_{xx}(t, 0) = u_{xx}(t, L) = 0. \end{aligned}$$
(3)

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

$$\begin{aligned} f_\mathrm {h}(u,u_\mathrm {h}, v,v_\mathrm {h}) = -K_\mathrm {h}([\langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^\alpha (1+\mu v_\mathrm {r}), \quad v_\mathrm {r} = \left\langle \varepsilon , v \right\rangle - v_\mathrm {h}, \end{aligned}$$

where we denote the felt stiffness coefficient with \(K_\mathrm {h}\), the hammer stiffness exponent with \(\alpha \) and the felt loss coefficient with \(\mu \). These coefficients also have a positive value. Here, \(\varepsilon (x)\) determines the hammer-string collision profile and satisfies \(\langle \varepsilon , 1 \rangle =1\) so that \(\varepsilon (x)\) becomes an approximation to the delta function. The expression \([ \xi ]^+\) is defined by

$$\begin{aligned}{}[ \xi ]^+ = {\left\{ \begin{array}{ll} \xi , &{} (\xi > 0) \\ 0, &{} (\text{ otherwise }) \end{array}\right. } \end{aligned}$$

and \(\langle \varepsilon , u\rangle \) denotes the \(L^2\) inner product

$$\begin{aligned} \langle \varepsilon , u\rangle := \int _0^L \varepsilon (x) u(t, x) \mathrm {d}x. \end{aligned}$$

The estimated values of the parameters for the C4 tone are listed in Tables 1 and 2. For parameters of other tones we refer e.g. to [3, 4, 18].

Table 1 Estimated parameter values used in [4] to describe the hammer model
Full size table
Table 2 Estimated parameter values used in [4] to specifying the C4 tone
Full size table

As it is shown below, when the damping terms are ignored, that is, \(\mu =\mathrm {d}_1=\mathrm {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_\mathrm {h}:= M_\mathrm {h}v_\mathrm {h}\), and an energy function

$$\begin{aligned}&{{H}}(u, v, {u}_\mathrm {h}, {p}_\mathrm {h})\nonumber \\&= \int _0^L \left[ \frac{1}{2} v^2 + \frac{c^2}{2 } u_x^2 + \frac{\kappa ^2}{2 } u_{xx}^2 \right] \mathrm {d}x + \frac{1}{2M_\mathrm {h}} {p_\mathrm {h}}^2 + \frac{K_\mathrm {h}}{\alpha + 1} ([ \langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^{\alpha +1}. \end{aligned}$$
(4)

This is a separable Hamiltonian in the sense that \({{H}}\) is written as a sum of two functions:

$$\begin{aligned} {{H}}(u, v, {u}_\mathrm {h}, {p}_\mathrm {h})&= T(v, p_\mathrm {h}) + U(u, u_\mathrm {h}), \\ T(v, p_\mathrm {h})&= \int _0^L \left[ \frac{1}{2} v^2 \right] \mathrm {d}x + \frac{1}{2M_\mathrm {h}} {p_\mathrm {h}}^2, \\ U(u, u_\mathrm {h})&= \int _0^L \left[ \frac{c^2}{2 } u_x^2 + \frac{\kappa ^2}{2 } u_{xx}^2 \right] \mathrm {d}x + \frac{K_\mathrm {h}}{\alpha + 1} ([ \langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^{\alpha +1}. \end{aligned}$$

Then it is straightforward to check the following theorem.

Theorem 1

The model (1, 2) is equivalent to

$$\begin{aligned} \begin{pmatrix} u_t \\ \dot{u}_\mathrm {h}\\ v_t \\ \dot{p}_\mathrm {h}\end{pmatrix}&= \begin{pmatrix} \phantom {-}0 &{} \phantom {-}0 &{} 1 &{} 0\\ \phantom {-}0 &{} \phantom {-}0 &{} 0 &{} 1\\ -1 &{} \phantom {-}0 &{} 0 &{} 0 \\ \phantom {-}0 &{} -1 &{} 0 &{} 0 \end{pmatrix} \,\, \begin{pmatrix} \frac{\delta {{H}}}{\delta u}&\frac{\partial {{H}}}{\partial u_\mathrm {h}}&\frac{\delta {{H}}}{\delta v}&\frac{\partial {{H}}}{\partial p_\mathrm {h}} \end{pmatrix}^\top \\&\qquad +\begin{pmatrix} 0 \\ 0 \\ -d_1 v + d_3 v_{xx} -\varepsilon K_\mathrm {h}([\langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^\alpha \mu v_\mathrm {r} \\ {K_\mathrm {h}([\langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^\alpha \mu v_\mathrm {r}} \end{pmatrix} \\&= \begin{pmatrix} v \\ \frac{p_\mathrm {h}}{M_\mathrm {h}} \\ c^2 u_{xx} - \kappa ^2 u_{xxxx} - {\varepsilon K_\mathrm {h}} ([\langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^\alpha (1 + \mu v_\mathrm {r}) - {d_1} v + {d_3} v_{xx} \\ K_\mathrm {h}([\langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^\alpha (1+\mu v_\mathrm {r}) \end{pmatrix}. \end{aligned}$$

Moreover, if the damping coefficients \(\mu , 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:

$$\begin{aligned} \frac{\mathrm {d}{{H}}}{\mathrm {d}t} =&- \mu v_\mathrm {r}^2 K_\mathrm {h}([\langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^{\alpha } - d_1 \int _0^L v^2 \,\mathrm {d}x - d_3 \int _0^L v_x^2 \,\mathrm {d}x \le 0. \end{aligned}$$

Moreover, if \(\mu = d_1 = d_3 = 0\), the system is conservative.

Proof

From the integration by parts and the boundary conditions \(u(t, 0)=u_{xx}(t, 0)=u(t, L)=u_{xx}(t, L)=0\), we obtain

$$\begin{aligned} \frac{\mathrm {d}{{H}}}{\mathrm {d}t} =&\int _0^L \left[ v v_t - c^2 u_{xx} u_t + \kappa ^2 u_{xxxx} u_t + K_\mathrm {h}([\langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^{\alpha } \varepsilon u_t \right] \mathrm {d}x \\&+ \left[ c^2 u_x u_t \right] _0^L + \left[ \kappa ^2 u_{xx} u_{xt} \right] _0^L - \left[ \kappa ^2 u_{xxx} u_{t} \right] _0^L + \frac{p_\mathrm {h}}{M_\mathrm {h}}{\dot{p}_\mathrm {h}} - K_\mathrm {h}([\langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^{\alpha } \dot{u}_\mathrm {h}\\ =&\int _0^L \left[ v v_t - c^2 u_{xx} u_t + \kappa ^2 u_{xxxx} u_t + K_\mathrm {h}([\langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^{\alpha } \varepsilon u_t \right] \mathrm {d}x \\&+ \frac{p_\mathrm {h}}{M_\mathrm {h}}{\dot{p}_\mathrm {h}} - K_\mathrm {h}([\langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^{\alpha } \dot{u}_\mathrm {h}. \end{aligned}$$

Substitution leads to

$$\begin{aligned} \frac{\mathrm {d}{{H}}}{\mathrm {d}t} =&\int _0^L \left[ v_t - c^2u_{xx} + \kappa ^2 u_{xxxx}+\varepsilon K_\mathrm {h}([\langle \varepsilon ,u \rangle -u_\mathrm {h}]^+)^\alpha ) \right] v \mathrm {d}x \\&+ \frac{p_\mathrm {h}}{M_\mathrm {h}} K_\mathrm {h}([\langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^{\alpha } (1 + \mu v_\mathrm {r}) - K_\mathrm {h}([\langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^{\alpha } \frac{p_\mathrm {h}}{M_\mathrm {h}} \\ =\,\,&\mu v_\mathrm {r}\frac{p_\mathrm {h}}{M_\mathrm {h}} K_\mathrm {h}([\langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^{\alpha } - \mu v_\mathrm {r}\langle \varepsilon , v \rangle K_\mathrm {h}([\langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^{\alpha } \\&- d_1 \int _0^L v^2 \,\mathrm {d}x + d_3 \int _0^L v v_{xx} \,\mathrm {d}x. \end{aligned}$$

Because \(v_\mathrm {r}= \langle \varepsilon , v\rangle - p_\mathrm {h}/M_\mathrm {h}\) holds and the boundary condition implies \(v(t, 0) = v(t, L) = 0\), we obtain

$$\begin{aligned} \frac{\mathrm {d}{{H}}}{\mathrm {d}t} =&- \mu v_\mathrm {r}^2 K_\mathrm {h}([\langle \varepsilon , u\rangle - u_\mathrm {h}]^+)^{\alpha } - d_1 \int _0^L v^2 \,\mathrm {d}x - d_3 \int _0^L v_x^2 \,\mathrm {d}x \le 0. \end{aligned}$$

\(\square \)

3 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.

Fig. 2
figure 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

Full size image

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 \(\varDelta 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\varDelta 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

$$\begin{aligned} \delta _x^+ q_l = \frac{q_{l+1}-q_l}{\varDelta x}, \quad \delta _x^- q_l = \frac{q_l-q_{l-1}}{\varDelta x}\quad \text{ and } \quad \delta _x^{\langle 2\rangle } q_l = \frac{q_{l+1}-2q_l+q_{l-1}}{\varDelta x^2}, \end{aligned}$$

respectively. We first consider the phase space:

$$\begin{aligned} \tilde{\mathscr {X}}{:=}\{ \tilde{\varvec{u}}\in {\mathbb {R}}^{N+3}\mid \tilde{\varvec{u}} = ({\tilde{u}}_{-1} \ {\tilde{u}}_0 \ {\tilde{u}}_1 \ \cdots \ {\tilde{u}}_{N-1} \ {\tilde{u}}_{N} \ {\tilde{u}}_{N+1} )^\top \}. \end{aligned}$$
(5)

For \(\tilde{\varvec{u}}\), \(\tilde{\varvec{v}} \in \tilde{\mathscr {X}}\), we define the discrete Hamiltonian \({{\tilde{H}}}_\mathrm {d}\) in \(\tilde{\mathscr {X}}\) using the trapezoidal rule to approximate the integral in H:

$$\begin{aligned}&{{\tilde{H}}}_\mathrm {d}(\tilde{\varvec{u}}, \tilde{\varvec{v}}, u_\mathrm {h}, p_\mathrm {h})\nonumber \\&=\left[ \frac{1}{2}v_0^2 + \frac{c^2}{2} \frac{\left( \delta _x^+u_0 \right) ^2 + \left( \delta _x^-u_0 \right) ^2}{2} + \frac{\kappa ^2}{2}\left( \delta _x^{\langle 2\rangle }u_0 \right) ^2 \right] \frac{\varDelta x}{2}\nonumber \\&\quad +\sum _{l=1}^{N-1} \left[ \frac{1}{2}v_l^2 + \frac{c^2}{2} \frac{\left( \delta _x^+u_l \right) ^2 + \left( \delta _x^-u_l \right) ^2}{2} +\frac{\kappa ^2}{2}\left( \delta _x^{\langle 2\rangle }u_l \right) ^2 \right] \varDelta x \nonumber \\&\quad + \left[ \frac{1}{2}v_N^2 + \frac{c^2}{2} \frac{\left( \delta _x^+u_N \right) ^2 + \left( \delta _x^-u_N \right) ^2}{2} + \frac{\kappa ^2}{2}\left( \delta _x^{\langle 2\rangle }u_N \right) ^2 \right] \frac{\varDelta x}{2}\nonumber \\&\quad +\frac{1}{2M_\mathrm {h}}p_\mathrm {h}^2 + \frac{K_\mathrm {h}}{\alpha +1}\left( [\langle \varvec{\varepsilon }, \tilde{\varvec{u}} \rangle _{\tilde{\mathscr {X}}} -u_\mathrm {h}]^{+} \right) ^{\alpha +1}, \end{aligned}$$
(6)

where \(\varvec{\varepsilon }\in \tilde{\mathscr {X}}\) is an approximation of the function \(\varepsilon (x)\) and \(\langle \cdot , \cdot \rangle _{\tilde{\mathscr {X}}}\) means

$$\begin{aligned} \langle \tilde{\varvec{u}}, \tilde{\varvec{v}} \rangle _{\tilde{\mathscr {X}}}={\tilde{u}}_0 {\tilde{v}}_0\frac{\varDelta x}{2}+\sum _{l=1}^{N-1} {\tilde{u}}_l {\tilde{v}}_l \varDelta x + {\tilde{u}}_N {\tilde{v}}_N\frac{\varDelta x}{2} \end{aligned}$$
(7)

for \(\tilde{\varvec{u}}, \tilde{\varvec{v}} \in \tilde{\mathscr {X}}\).

Remark 2

Regarding the approximation to the term \(u_x\) in H, we chose

$$\begin{aligned} u_x(t,l\varDelta x)^2\approx \frac{(\delta _x^+u_l)^2+(\delta _x^-u_l)^2}{2}, \end{aligned}$$

which is a typical choice in [26]; however, we can also use, for example, the central difference and define the Hamiltonian as

$$\begin{aligned}&\tilde{{\tilde{H}}}_\mathrm {d}(\tilde{\varvec{u}}, \tilde{\varvec{v}}, u_\mathrm {h}, p_\mathrm {h})\nonumber \\&=\left[ \frac{1}{2}v_0^2 + \frac{c^2}{2} \left( \frac{\delta _x^++\delta _x^-}{2}u_0\right) ^2 + \frac{\kappa ^2}{2}\left( \delta _x^{\langle 2\rangle }u_0 \right) ^2 \right] \frac{\varDelta x}{2}\nonumber \\&\quad +\sum _{l=1}^{N-1} \left[ \frac{1}{2}v_l^2 + \frac{c^2}{2} \left( \frac{\delta _x^++\delta _x^-}{2}u_l\right) ^2 +\frac{\kappa ^2}{2}\left( \delta _x^{\langle 2\rangle }u_l \right) ^2 \right] \varDelta x \nonumber \\&\quad + \left[ \frac{1}{2}v_N^2 + \frac{c^2}{2} \left( \frac{\delta _x^++\delta _x^-}{2}u_N\right) ^2 + \frac{\kappa ^2}{2}\left( \delta _x^{\langle 2\rangle }u_N \right) ^2 \right] \frac{\varDelta x}{2}\nonumber \\&\quad +\frac{1}{2M_\mathrm {h}}p_\mathrm {h}^2 + \frac{K_\mathrm {h}}{\alpha +1}\left( [\langle \varvec{\varepsilon }, \tilde{\varvec{u}} \rangle _{\tilde{\mathscr {X}}} -u_\mathrm {h}]^{+} \right) ^{\alpha +1}. \end{aligned}$$

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 use \({{\tilde{H}}}_\mathrm {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 using \(\tilde{{{\tilde{H}}}}_\mathrm {d}\) converge to the exact ones slower than that derived by using \({{{\tilde{H}}}}_\mathrm {d}\).

The boundary conditions corresponding to (3) are imposed by

$$\begin{aligned} u_0=u_N=0, \quad u_{1}=-u_{-1}, \quad \text{ and } \quad u_{N-1}=-u_{N+1}, \end{aligned}$$
(8)

and by using them, we can rewrite \({{\tilde{H}}}_\mathrm {d}\) without the boundary terms \(u_{-1}\), \(u_0\), \(u_N\) and \(u_{N+1}\).

This is equivalent to the restriction of \({{\tilde{H}}}_\mathrm {d}\) to the subspace of \(\tilde{\mathscr {X}}\):

$$\begin{aligned} {\mathscr {X}} = \{ \varvec{u} \in \tilde{\mathscr {X}} \mid u_0=u_{N}=0, \quad u_{1}=-u_{-1}, \quad u_{N-1}=-u_{N+1} \} \simeq {\mathbb {R}}^{N-1}. \end{aligned}$$

We denote this restricted discrete Hamiltonian with \(H_\mathrm {d}\). For \(\varvec{u}\), \(\varvec{v} \in \mathscr {X}\), \(H_\mathrm {d}\) is defined as

$$\begin{aligned}&H_\mathrm {d}(\varvec{u}, \varvec{v}, u_\mathrm {h}, p_\mathrm {h})\nonumber \\&\quad = \left[ \frac{c^2}{2} \left( \frac{u_1}{\varDelta x} \right) ^2 \right] \frac{\varDelta x}{2} \nonumber \\&\quad + \left[ \frac{c^2}{2} \frac{\left( \delta _x^+u_1 \right) ^2 + \left( u_1/\varDelta x \right) ^2}{2} + \frac{\kappa ^2}{2} \left( \frac{u_2-2u_1}{\varDelta x^2} \right) ^2 \right] \varDelta x \nonumber \\&\quad + \sum _{l=2}^{N-2} \left[ \frac{c^2}{2} \frac{\left( \delta _x^+u_l \right) ^2 + \left( \delta _x^-u_l \right) ^2}{2} +\frac{\kappa ^2}{2}\left( \delta _x^{\langle 2\rangle }u_l \right) ^2 \right] \varDelta x \nonumber \\&\quad + \left[ \frac{c^2}{2} \frac{\left( {u_{N-1}}/{\varDelta x} \right) ^2 + \left( \delta _x^-u_{N-1} \right) ^2}{2} +\frac{\kappa ^2}{2}\left( \frac{-2u_{N-1}+u_{N-2}}{\varDelta x^2} \right) ^2 \right] \varDelta x \nonumber \\&\quad + \left[ \frac{c^2}{2} \left( \frac{u_{N-1}}{\varDelta x} \right) ^2 \right] \frac{\varDelta x}{2} \nonumber \\&\quad +\sum _{l=1}^{N-1}\frac{1}{2}v_l^2\varDelta x+\frac{1}{2M_\mathrm {h}}p_\mathrm {h}^2 + \frac{K_\mathrm {h}}{\alpha +1}\left( [\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}} -u_\mathrm {h}]^{+} \right) ^{\alpha +1}, \end{aligned}$$
(9)

where \(\langle \cdot , \cdot \rangle _{\mathscr {X}}\) means

$$\begin{aligned} \langle \varvec{u},\varvec{v}\rangle _{\mathscr {X}} = \sum _{l=1}^{N-1} u_lv_l \varDelta x\quad \hbox { for all }\quad \varvec{u} , \varvec{v}\in \mathscr {X}. \end{aligned}$$
(10)

We note that this is equivalent to (7) under the boundary condition (8).

We calculate the partial derivatives of \(H_\mathrm {d}\) to obtain the gradient in the Hamilton equation. The partial derivatives of \(H_\mathrm {d}\) with respect to \(u_l\) and \(v_l\) \((l=1,\ldots ,N-1)\) are

$$\begin{aligned} \frac{\partial H_\mathrm {d}}{\partial u_1}&= \Bigl \{ - c^2\frac{-2u_1+u_2}{\varDelta x^2} + \kappa ^2\frac{5u_1-4u_2+u_3}{\varDelta x^4} + \varepsilon _1 K_\mathrm {h}([\langle \varvec{\varepsilon },\varvec{u}\rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha \Bigr \}\varDelta x, \\ \frac{\partial H_\mathrm {d}}{\partial u_2}&= \Bigl \{ - c^2\frac{u_1-2u_2+u_3}{\varDelta x^2}+\kappa ^2\frac{-4u_1+6u_2-4u_3+u_4}{\varDelta x^4}\\&\quad + \varepsilon _2 K_\mathrm {h}([\langle \varvec{\varepsilon },\varvec{u}\rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha \Bigr \}\varDelta x,\\ \frac{\partial H_\mathrm {d}}{\partial u_l}&= \Bigl \{ - c^2\frac{u_{l-1}-2u_l+u_{l+1}}{\varDelta x^2} + \kappa ^2\frac{u_{l-2}-4u_{l-1}+6u_l-4u_{l+1}+u_{l+2}}{\varDelta x^4}\\&\quad + \varepsilon _l K_\mathrm {h}([\langle \varvec{\varepsilon },\varvec{u}\rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha \Bigr \}\varDelta x, \quad (l=3,4,\dots ,N-3) \\ \frac{\partial H_\mathrm {d}}{\partial u_{N-2}}&= \Bigl \{ - c^2\frac{u_{N-3}-2u_{N-2}+u_{N-1}}{\varDelta x^2} + \kappa ^2\frac{u_{N-4}-4u_{N-3}+6u_{N-2}-4u_{N-1}}{\varDelta x^4} \\&\quad + \varepsilon _{N-2} K_\mathrm {h}([\langle \varvec{\varepsilon },\varvec{u}\rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha \Bigr \}\varDelta x,\\ \frac{\partial H_\mathrm {d}}{\partial u_{N-1}}&= \Bigl \{ - c^2\frac{u_{N-2}-2u_{N-1}}{\varDelta x^2} + \kappa ^2\frac{ u_{N-3}-4u_{N-2}+5u_{N-1}}{\varDelta x^4} \\&\quad + \varepsilon _{N-1} K_\mathrm {h}([\langle \varvec{\varepsilon },\varvec{u}\rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha \Bigr \}\varDelta x, \\ \frac{\partial H_\mathrm {d}}{\partial v_l}&=v_l\varDelta x, \quad (l=1,2,\dots ,N-1) \end{aligned}$$

and the partial derivatives with respect to \(u_\mathrm {h}\) and \(p_\mathrm {h}\) are

$$\begin{aligned} \frac{\partial H_\mathrm {d}}{\partial u_\mathrm {h}}&=-K_\mathrm {h}([\langle \varvec{\varepsilon },\varvec{u}\rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha , \\ \frac{\partial H_\mathrm {d}}{\partial p_\mathrm {h}}&=\frac{p_\mathrm {h}}{M_\mathrm {h}}. \end{aligned}$$

We denote the difference matrices of order \(N-1\) with

$$\begin{aligned} D_2= \frac{1}{\varDelta x^2} \begin{pmatrix} -2 &{} 1 &{} &{} &{} \\ 1 &{} -2 &{} 1 &{} &{} \\ &{} \ddots &{} \ddots &{} \ddots &{} \\ &{} &{} 1 &{} -2 &{} 1 \\ &{} &{} &{} 1 &{} -2 \\ \end{pmatrix},\quad \ D_4= \frac{1}{\varDelta x^4} \begin{pmatrix} 5 &{} -4 &{} 1 &{} &{} &{} &{} \\ -4 &{} 6 &{} -4 &{} 1 &{} &{} &{} \\ 1 &{} -4 &{} 6 &{} -4 &{} 1 &{} &{} \\ &{} \ddots &{} \ddots &{} \ddots &{} \ddots &{} \ddots &{} \\ &{} &{} 1 &{} -4 &{} 6 &{} -4 &{} 1 \\ &{} &{} &{} 1 &{} -4 &{} 6 &{} -4 \\ &{} &{} &{} &{} 1 &{} -4 &{} 5 \\ \end{pmatrix}. \end{aligned}$$

The semi-discrete scheme is defined by the following separable Hamiltonian system:

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t} \begin{pmatrix} \varvec{u} \\ u_\mathrm {h}\\ \varvec{v} \\ p_\mathrm {h}\end{pmatrix}&= \begin{pmatrix} &{}&{}I&{}\varvec{0}\\ &{}&{}{\varvec{0}}^\top &{}1\\ -I&{}\phantom {-}\varvec{0}&{}&{}\\ {\varvec{0}}^\top &{}-1&{}&{} \end{pmatrix} \,\, \begin{pmatrix} \nabla _{\varvec{u}} H_\mathrm {d}^\top&\nabla _{u_\mathrm {h}} H_\mathrm {d}&\nabla _{\varvec{v}} H_\mathrm {d}^\top&\nabla _{p_\mathrm {h}} H_\mathrm {d}\end{pmatrix}^\top \nonumber \\&=\begin{pmatrix} &{}&{}I&{}\varvec{0}\\ &{}&{}{\varvec{0}}^\top &{}1\\ -I&{}\phantom {-}\varvec{0}&{}&{}\\ {\varvec{0}}^\top &{}-1&{}&{} \end{pmatrix} \,\, \begin{pmatrix} -c^2 D_2 \varvec{u}+\kappa ^2D_4\varvec{u} + \varvec{\varepsilon }K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^{\alpha }\\ - K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha \\ \varvec{v}\\ \frac{p_\mathrm {h}}{M_\mathrm {h}} \end{pmatrix}, \end{aligned}$$
(11)

where \(\nabla _q H_\mathrm {d}\) means the gradient of \(H_\mathrm {d}\) in the q direction associated with the inner product (10). In the following theorem we address the energy behavior of this model.

Theorem 3

The semi-discretized Hamiltonian system (11) preserves the discretized energy function \(H_\mathrm {d}\) (9), i.e.

$$\begin{aligned} \frac{\mathrm {d}H_\mathrm {d}}{\mathrm {d}t}=0. \end{aligned}$$

Remark 3

This theorem is generalized to include the damping terms in Theorem 5, and a proof is given there.

4 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

$$\begin{aligned} \frac{\mathrm {d}y(t)}{\mathrm {d}t} =F(y(t), z(t)), \quad \frac{\mathrm {d}z(t)}{\mathrm {d}t}=G(y(t), z(t)). \end{aligned}$$
(12)

Definition 1

Let \(a_{ij}, \ b_i\), and \({\hat{a}}_{ij}, \ {\hat{b}}_i\) be the coefficients of two Runge–Kutta methods. Then, an s-stage PRK method for (12) with a step size \(\varDelta t\) is given by

$$\begin{aligned} k_i&=F\left( y^n+\varDelta t \sum _{j=1}^{s}a_{ij}k_j, z^n+\varDelta t \sum _{j=1}^{s}{\hat{a}}_{ij}l_j\right) , \nonumber \\ l_i&=G\left( y^n+\varDelta t \sum _{j=1}^{s}a_{ij}k_j, z^n+\varDelta t \sum _{j=1}^{s}{\hat{a}}_{ij}l_j\right) , \nonumber \\ y^{n+1}&=y^n+\varDelta t \sum _{i=1}^{s}b_ik_i, \quad z^{n+1}=z^n+\varDelta t \sum _{i=1}^{s}{\hat{b}}_il_i. \end{aligned}$$
(13)

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].

Theorem 4

(Symplectic Partitioned Runge–Kutta (SPRK) Method) An s-stage PRK method (13) is symplectic if it satisfies the conditions

$$\begin{aligned} {\left\{ \begin{array}{ll} b_i={\hat{b}}_i,\\ b_i {\hat{a}}_{ij}+{\hat{b}}_j a_{ji}-b_i{\hat{b}}_j=0, \end{array}\right. }\quad \hbox { for }\quad \ i=1, \dots , s, \quad j=1, \dots , s. \end{aligned}$$

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 these coefficients becomes symplectic and explicit for separable Hamiltonian systems. Following [50], we write the coefficients of Table 3 as

$$\begin{aligned} (b_1, b_2, \dots , b_s)[{\hat{b}}_1, {\hat{b}}_2, \dots , {\hat{b}}_s]. \end{aligned}$$

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.

Table 3 The Butcher tableau of an s-stage SPRK method
Full size table

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.

figure a

We tested three numerical schemes for the computation of piano sounds. These schemes are obtained by applying the SPRK methods to

$$\begin{aligned}&\frac{\mathrm {d}}{\mathrm {d}t} \begin{pmatrix} \varvec{u}\\ u_\mathrm {h}\\ \varvec{v}\\ p_\mathrm {h}\end{pmatrix}= \begin{pmatrix} \phantom {-}O&II&O \end{pmatrix} \begin{pmatrix} {F(\varvec{u}, u_\mathrm {h},\varvec{v}, p_\mathrm {h})}\\ {G(\varvec{v}, p_\mathrm {h})} \end{pmatrix}, \\&\quad F({\varvec{u}, u_\mathrm {h},\varvec{v}, p_\mathrm {h}}) \nonumber \\&\quad := \begin{pmatrix} -c^2 D_2 \varvec{u}+\kappa ^2D_4\varvec{u}+\varvec{\varepsilon }K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^{\alpha }(1+\mu v_{\mathrm {r}, \mathrm {d}})+d_1\varvec{v}-d_3D_2\varvec{v}\\ -K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha (1+\mu v_{\mathrm {r}, \mathrm {d}}) \end{pmatrix},\nonumber \\&G({\varvec{v}, p_\mathrm {h}}):= \begin{pmatrix} \varvec{v}^\top&\frac{p_\mathrm {h}}{M_\mathrm {h}} \end{pmatrix}^\top ,\quad v_{\mathrm {r}, \mathrm {d}} := \langle \varvec{\varepsilon }, \varvec{v} \rangle _{\mathscr {X}} - \frac{p_\mathrm {h}}{M_\mathrm {h}}, \nonumber \end{aligned}$$
(14)

or, equivalently,

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{u}_l = v_l, \quad (l=1,2,\dots ,N-1)\\ \dot{u}_\mathrm {h}= \frac{p_\mathrm {h}}{M_\mathrm {h}},\\ \dot{v}_1 = c^2\frac{-2u_1+u_2}{\varDelta x^2} - \kappa ^2\frac{5u_1-4u_2+u_3}{\varDelta x^4}\\ \qquad \quad -\,\, \varepsilon _1 K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}} -u_\mathrm {h}]^+)^\alpha (1+\mu v_{\mathrm {r},\mathrm {d}})-d_1v_1+d_3\frac{-2v_1+v_2}{\varDelta x^2},\\ \dot{v}_2 = c^2 \delta _x^{\langle 2\rangle }u_2 - \kappa ^2\frac{-4u_1+6u_2-4u_3+u_4}{\varDelta x^4}\\ \qquad \quad -\,\, \varepsilon _2 K_\mathrm {h}([\langle \varvec{\varepsilon },\varvec{u}\rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha (1+\mu v_{\mathrm {r},\mathrm {d}})-d_1v_2+d_3\delta _x^{\langle 2\rangle }v_2,\\ \dot{v}_l = c^2 \delta _x^{\langle 2\rangle }u_l - \kappa ^2 \delta _x^{\langle 4\rangle }u_l - \varepsilon _l K_\mathrm {h}([\langle \varvec{\varepsilon },\varvec{u}\rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha (1+\mu v_{\mathrm {r},\mathrm {d}})\\ \qquad \quad -\,\, d_1 v_l + d_3\delta _x^{\langle 2\rangle }v_l, \quad (l=3,\ldots ,N-3)\\ \dot{v}_{N-2} = c^2 \delta _x^{\langle 2\rangle }u_{N-2} - \kappa ^2\frac{u_{N-4}-4u_{N-3}+6u_{N-2}-4u_{N-1}}{\varDelta x^4} \\ \qquad \quad -\,\, \varepsilon _{N-2} K_\mathrm {h}([\langle \varvec{\varepsilon },\varvec{u}\rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha (1+\mu v_{\mathrm {r},\mathrm {d}}) - d_1 v_{N-2} + d_3 \delta _x^{\langle 2\rangle }v_{N-2},\\ \dot{v}_{N-1} = c^2\frac{u_{N-2}-2u_{N-1}}{\varDelta x^2} - \kappa ^2\frac{u_{N-3}-4u_{N-2}+5u_{N-1}}{\varDelta x^4} \\ \qquad \quad -\,\, \varepsilon _{N-1} K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}} -u_\mathrm {h}]^+)^\alpha (1+\mu v_{\mathrm {r},\mathrm {d}}) - d_1 v_{N-1} + d_3 \frac{v_{N-2}-2v_{N-1}}{\varDelta x^2},\\ \dot{p}_\mathrm {h}= K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}} -u_\mathrm {h}]^+)^\alpha (1+\mu v_{\mathrm {r},\mathrm {d}}), \end{array}\right. } \end{aligned}$$
(15)

where \(\delta _x^{\langle 4\rangle }:=\delta _x^{\langle 2\rangle }\delta _x^{\langle 2\rangle }\). 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_\mathrm {d}\) (9) is not increasing:

$$\begin{aligned} \frac{\mathrm {d}H_\mathrm {d}}{\mathrm {d}t}=&-\,\mu v_{\mathrm {r}, \mathrm {d}}^2K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha - d_1\sum _{l=1}^{N-1} v_l^2\varDelta x \\&-\, d_3\Biggl [ \frac{(\delta _x^+v_0)^2 + (\delta _x^-v_0)^2}{2}\frac{\varDelta x}{2} +\sum _{l=1}^{N-1} \frac{(\delta _x^+v_l)^2+(\delta _x^-v_l)^2}{2}\varDelta x \\&+\,\frac{(\delta _x^+v_N)^2 + (\delta _x^-v_N)^2}{2}\frac{\varDelta x}{2} \Biggr ]\le 0. \end{aligned}$$

Moreover, if \(\mu =d_1=d_3=0\), the system is conservative according to Theorem 3.

Proof

From the differentiation of \(H_\mathrm {d}\), we obtain

$$\begin{aligned} \frac{\mathrm {d}H_\mathrm {d}}{\mathrm {d}t}&=\frac{\mathrm {d}{{\tilde{H}}}_\mathrm {d}}{\mathrm {d}t}\\&=\left[ v_0\dot{v}_0 + c^2 \frac{(\delta _x^+u_0)(\delta _x^+\dot{u}_0) + (\delta _x^-u_0)(\delta _x^-\dot{u}_0)}{2} + \kappa ^2 (\delta _x^{\langle 2\rangle }u_0)(\delta _x^{\langle 2\rangle }\dot{u}_0) \right] \frac{\varDelta x}{2}\\&\quad + \sum _{l=1}^{N-1} \Biggl [ v_l\dot{v}_l + c^2 \frac{(\delta _x^+u_l)(\delta _x^+\dot{u}_l) + (\delta _x^-u_l)(\delta _x^-\dot{u}_l)}{2} + \kappa ^2(\delta _x^{\langle 2\rangle }u_l)(\delta _x^{\langle 2\rangle }\dot{u}_l) \Biggr ]\varDelta x\\&\quad +\, \left[ v_N\dot{v}_N + c^2 \frac{(\delta _x^+u_{N})(\delta _x^+\dot{u}_{N}) + (\delta _x^-u_{N})(\delta _x^-\dot{u}_{N})}{2} + \kappa ^2 (\delta _x^{\langle 2\rangle }u_N)(\delta _x^{\langle 2\rangle }\dot{u}_N) \right] \frac{\varDelta x}{2}\\&\quad +\, \frac{p_\mathrm {h}}{M_\mathrm {h}}\dot{p}_\mathrm {h}\,+\, K_\mathrm {h}([\langle \varvec{\varepsilon },\varvec{u}\rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha (\langle \varvec{\varepsilon },\dot{\varvec{u}}\rangle _{\mathscr {X}}-\dot{u}_\mathrm {h}). \end{aligned}$$

Instead of the integration by parts, we introduce the summation by parts [22]:

$$\begin{aligned}&a_0(\delta _x^+b_0) \frac{\varDelta x}{2} + \sum _{l=1}^{N-1} a_l(\delta _x^+b_l) \varDelta x + a_N(\delta _x^+b_N) \frac{\varDelta x}{2} \\&\quad + (\delta _x^-a_0)b_0 \frac{\varDelta x}{2} + \sum _{l=1}^{N-1} (\delta _x^-a_l)b_l \varDelta x + (\delta _x^-a_N)b_N \frac{\varDelta x}{2} \\&= -\frac{1}{2}\left( a_0b_1+a_{-1}b_0 \right) + \frac{1}{2}\left( a_Nb_{N+1}+a_{N-1}b_N \right) . \end{aligned}$$

Using the summation by parts, we obtain

$$\begin{aligned} \frac{\mathrm {d}H_\mathrm {d}}{\mathrm {d}t}&= -\Biggl [ \frac{c^2}{2}\left\{ \frac{(\delta _x^+u_0)\dot{u}_1 + (\delta _x^+u_{-1})\dot{u}_0}{2} + \frac{\dot{u}_0(\delta _x^-u_1) + \dot{u}_{-1}(\delta _x^-u_0)}{2} \right\} \\&\quad + \kappa ^2 \frac{(\delta _x^{\langle 2\rangle }u_0)(\delta _x^-\dot{u}_1) + (\delta _x^{\langle 2\rangle }u_{-1})(\delta _x^-\dot{u}_0)}{2} \Biggr ] \\&\quad + \kappa ^2 \frac{\dot{u}_0(\delta _x^-\delta _x^{\langle 2\rangle }u_1) + \dot{u}_{-1}(\delta _x^-\delta _x^{\langle 2\rangle }u_0)}{2} \\&\quad +\Biggl [ v_0 \dot{v}_0 - c^2(\delta _x^{\langle 2\rangle }u_0) \dot{u}_0 + \kappa ^2 (\delta _x^{\langle 4\rangle }u_0)\dot{u}_0 + K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}} - u_\mathrm {h}]^+)^\alpha \varepsilon _0\dot{u}_0 \Biggr ]\frac{\varDelta x}{2}\\&\quad + \sum _{l=1}^{N-1} \Biggl [ v_l \dot{v}_l - c^2(\delta _x^{\langle 2\rangle }u_l) \dot{u}_l + \kappa ^2 (\delta _x^{\langle 4\rangle }u_l)\dot{u}_l + K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}} - u_\mathrm {h}]^+)^\alpha \varepsilon _l\dot{u}_l \Biggr ]\varDelta x \\&\quad +\Biggl [ v_N \dot{v}_N - c^2(\delta _x^{\langle 2\rangle }u_N) \dot{u}_N + \kappa ^2 (\delta _x^{\langle 4\rangle }u_N)\dot{u}_N + K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}} - u_\mathrm {h}]^+)^\alpha \varepsilon _N\dot{u}_N \Biggr ] \frac{\varDelta x}{2} \\&\quad -\kappa ^2 \frac{\dot{u}_N(\delta _x^-\delta _x^{\langle 2\rangle }u_{N+1}) + \dot{u}_{N-1}(\delta _x^-\delta _x^{\langle 2\rangle }u_N)}{2} \\&\quad + \Biggl [ \frac{c^2}{2}\left\{ \frac{(\delta _x^+u_N)\dot{u}_{N+1} + (\delta _x^+u_{N-1})\dot{u}_N}{2} + \frac{\dot{u}_N(\delta _x^-u_{N+1}) + \dot{u}_{N-1}(\delta _x^-u_N)}{2} \right\} \\&\quad + \kappa ^2 \frac{(\delta _x^{\langle 2\rangle }u_N)(\delta _x^-\dot{u}_{N+1}) + (\delta _x^{\langle 2\rangle }u_{N-1})(\delta _x^-\dot{u}_N)}{2} \Biggr ] \\&\quad + \frac{p_\mathrm {h}}{M_\mathrm {h}} \dot{p}_\mathrm {h}- K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}} - u_\mathrm {h}]^+)^\alpha \dot{u}_\mathrm {h}\\&= \sum _{l=1}^{N-1} \Biggl [ v_l \dot{v}_l - c^2(\delta _x^{\langle 2\rangle }u_l) \dot{u}_l + \kappa ^2 (\delta _x^{\langle 4\rangle }u_l)\dot{u}_l + K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}} - u_\mathrm {h}]^+)^\alpha \varepsilon _l\dot{u}_l \Biggr ]\varDelta x \\&\quad + \frac{p_\mathrm {h}}{M_\mathrm {h}} \dot{p}_\mathrm {h}- K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}} - u_\mathrm {h}]^+)^\alpha \dot{u}_\mathrm {h}. \end{aligned}$$

We used the boundary conditions (8) at the last equality. Substituting the semi-discretized scheme (15), we get

$$\begin{aligned} \frac{\mathrm {d}H_\mathrm {d}}{\mathrm {d}t}&= \sum _{l=1}^{N-1} \left\{ \dot{v}_l - c^2 (\delta _x^{\langle 2\rangle }u_l) + \kappa ^2 (\delta _x^{\langle 4\rangle }u_l) + K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}} - u_\mathrm {h}]^+)^\alpha \varepsilon _l \right\} v_l \varDelta x \\&\quad + \frac{p_\mathrm {h}}{M_\mathrm {h}} K_\mathrm {h}([\langle \varvec{\varepsilon },\varvec{u}\rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha (1+\mu v_{\mathrm {r},\mathrm {d}}) - K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}} - u_\mathrm {h}]^+)^\alpha \frac{p_\mathrm {h}}{M_\mathrm {h}}\\&= \mu v_{\mathrm {r},\mathrm {d}}\frac{p_\mathrm {h}}{M_\mathrm {h}} K_\mathrm {h}([\langle \varvec{\varepsilon }, \varvec{u} \rangle _{\mathscr {X}} - u_\mathrm {h}]^+)^\alpha - \mu v_{\mathrm {r},\mathrm {d}} \langle \varvec{\varepsilon }, \varvec{v} \rangle _{\mathscr {X}}K_\mathrm {h}([\langle \varvec{\varepsilon },\varvec{u} \rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha \\&\quad - d_1 \sum _{l=1}^{N-1} v_l^2 \varDelta x + d_3 \sum _{l=1}^{N-1} (\delta _x^{\langle 2\rangle }v_l) v_l \varDelta x\\&= - \mu v_{\mathrm {r},\mathrm {d}}^2 K_\mathrm {h}([\langle \varvec{\varepsilon },\varvec{u} \rangle _{\mathscr {X}}-u_\mathrm {h}]^+)^\alpha - d_1 \sum _{l=1}^{N-1} v_l^2 \varDelta x \\&\quad - d_3 \Biggl [ \frac{(\delta _x^+v_0)^2+(\delta _x^-v_0)^2}{2} \frac{\varDelta x}{2} + \sum _{l=1}^{N-1} \frac{(\delta _x^+v_l)^2+ (\delta _x^-v_l)^2}{2}\varDelta x \\&\quad + \frac{(\delta _x^+v_N)^2+(\delta _x^-v_N)^2}{2} \frac{\varDelta x}{2} \Biggr ] \le 0. \end{aligned}$$

\(\square \)

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:

$$\begin{aligned}&(-1/24, 3/4, 7/24)[1, -2/3, 2/3],\\&(\omega , \nu , \omega , 0)[\omega /2, (\omega +\nu )/2, (\omega +\nu )/2, \omega /2],\\&(-1/48, 3/8, 7/24, 3/8, -1/48, 0)[1/2, -1/3, 1/3, 1/3, -1/3, 1/2], \end{aligned}$$

where \(\omega =(2+2^{1/3}+2^{-1/3})/3\) and \(\nu =1-2\omega \); see [51]. The initial conditions are

$$\begin{aligned} u(0,x)=v(0,x)=0, \quad u_\mathrm {h}(0)=0.0005, \quad p_\mathrm {h}=-4.0\,M_\mathrm {h}. \end{aligned}$$

We employed the approximated delta function \(\varvec{\varepsilon }(l_\mathrm {h})\) used in [9]:

$$\begin{aligned} (\varvec{\varepsilon }(l_\mathrm {h}))_l= {\left\{ \begin{array}{ll} -\beta _{\mathrm {h}}(\beta _{\mathrm {h}}-1)(\beta _{\mathrm {h}}-2)/(6\varDelta x), \quad &{}\text{( }l\,=\,l_{\mathrm {h}}-1\text{) }\\ (\beta _{\mathrm {h}}-1)(\beta _{\mathrm {h}}+1)(\beta _{\mathrm {h}}-2)/(2\varDelta x), \quad &{}\text{( }l\,=\,l_\mathrm {h}\text{) }\\ -\beta _{\mathrm {h}}(\beta _{\mathrm {h}}+1)(\beta _{\mathrm {h}}-2)/(2\varDelta x), \quad &{} \text{( }l\,=\,l_{\mathrm {h}}+1\text{) }\\ \beta _{\mathrm {h}}(\beta _{\mathrm {h}}-1)(\beta _{\mathrm {h}}+1)/(6\varDelta x), \quad &{} \text{( }l\,=\,l_{\mathrm {h}}+2\text{) }\\ 0. \quad &{} \text{(otherwise) } \end{array}\right. } \end{aligned}$$

According to [9], this is a 3rd-order approximation of \(\delta (x-(\beta _\mathrm {h}+l_\mathrm {h})\varDelta 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_\mathrm {h}=0.2N\) and \(\beta _\mathrm {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 shows the comparison of the numerical solutions by the 4-stage 4th-order SPRK under the various values of \(\varDelta t\) and \(\varDelta x\) with the following exact solution under the boundary condition (3):

$$\begin{aligned} u(t,x)&= \sum _{m=1}^{15} \sin \left( m\pi \frac{x}{L}\right) \left( A_m\cos (\sqrt{\lambda _m} t) + B_m\sin (\sqrt{\lambda _m} t) \right) , \end{aligned}$$
(16)
$$\begin{aligned} \lambda _m&=c^2\left( \frac{m\pi }{L}\right) ^2+\kappa ^2\left( \frac{m\pi }{L}\right) ^4. \end{aligned}$$
(17)

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 \(\varDelta t\rightarrow 0\) and \(\varDelta x\rightarrow 0\). This is in fact quantitatively confirmed in Table 4, where the \(L^2\)-norm and the \(L^\infty \)-norm are defined by

$$\begin{aligned}&||\varvec{e}(t)||_2=\left( \sum _{l=1}^{N-1} | e_l(t) |^2 \varDelta x\right) ^{1/2}, \quad ||\varvec{e}(t)||_\infty =\mathrm {max}_{l=1,\ldots ,N-1} | e_l(t) |,\\&\varvec{e}(t) = (e_1(t) \ldots \ e_{N-1}(t)), \quad e_l(t)={{\hat{u}}}_l^{t/\varDelta t}-u(t,l\varDelta x). \end{aligned}$$
Fig. 3
figure 3

Comparison of the numerical solutions by the 4-stage 4th-order SPRK with the exact solution by using the piano model without the terms regarding the hammer and the damping of the string. This graph shows the numerical solution \(\varvec{u}\) at \(t=0.1\) in the various conditions. The numerical solutions indeed converge to the exact solution as \(\varDelta t \rightarrow 0\) and \(\varDelta x \rightarrow 0\)

Full size image
Table 4 Errors of the numerical solution of (11) without the hammer and the damping terms by the 4-stage 4th-order SPRK
Full size table
Fig. 4
figure 4

The numerical solutions of the scheme derived by using \(\tilde{{\tilde{H}}}_\mathrm {d}\). Although \(\varDelta t\) and \(\varDelta x\) are the same or less, the waveform of that by \({{\tilde{H}}}_\mathrm {d}\) is closer to the exact solution than the other two results

Full size image
Fig. 5
figure 5

Comparison of the 4-stage 4th-order SPRK with the exact solution by using the piano model without the terms regarding the hammer and the damping of the string. These graphs show the gap of the discrete energy between the computed and the exact value with \(A_m=B_m=10^{-5}\) and \(N=100\). The discrete energy also converges to the exact value as \(\varDelta t\rightarrow 0\)

Full size image
Fig. 6
figure 6

The waveforms of the string computed with the three SPRK methods with \(N=1000\) for 1 s. The graph on the top is the numerical solution computed with the 4-stage 4th-order method, the one in the middle is computed with the 3-stage 3rd-order method and the bottom one with the 6-stage 4th-order method

Full size image

We also show the numerical solutions of the scheme derived by using \(\tilde{{\tilde{H}}}_\mathrm {d}\) in Fig. 4. Although we used the same or less \(\varDelta t\) and \(\varDelta x\), the waveform of that by \({{\tilde{H}}}_\mathrm {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 \(\varDelta t=1/(44{,}100\cdot 200)\) and \(N=1000\) hereinafter if it is not specifically noted.

Fig. 7
figure 7

The waveform of the string computed with the 4-stage 4th-order SPRK method with \(N=1000\). This is the waveform from 0.4 to 0.45 s. We find that some waves are regularly repeated and expect that this includes the integer multiple of the frequency of 261.63 Hz which corresponds to the note C4

Full size image
Fig. 8
figure 8

The spectrum of the numerical solution computed with the 4-stage 4th-order SPRK method with \(N=1000\). The spectrum under 10,000 Hz is shown at the top and the enlarged one under 2000 Hz at the bottom. The frequency of a note C4 is 261.63 Hz and the large peaks are observed near the integer multiples of 261.63 Hz

Full size image

Figure 6 shows the numerical solutions obtained by the 3-stage 3rd-order, the 4-stage 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_\mathrm {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\times 10^{-2}\). We excluded the damping terms in this numerical test so that the energy is preserved.

Fig. 9
figure 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 \(\varDelta t = 1/(44{,}100 \cdot 10)\), \(1/(44{,}100 \cdot 100)\). A similar result to Fig. 3 which shows the convergence of the energy as \(\varDelta t\rightarrow 0\) is observed

Full size image

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 behavior 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 \(\varDelta t\rightarrow 0\) is shown is again observed in this test.

Fig. 10
figure 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

Full size image
Fig. 11
figure 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

Full size image
Fig. 12
figure 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 \(\varDelta t\) tends to 0

Full size image

Figures  10, 11, 12 show the result when the number of points N is changed from 1000 to 50. We used \(l_{\mathrm {h}}=0.2N\) and \(\beta _\mathrm {h}=0.00945\) so that the hammer strikes the same position (\(x \approx 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; 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\times 10^{-2}\), is shown in Fig. 12. The values of \(H_\mathrm {d}\) are still within a certain fixed range and converge to the exact value by \(\varDelta t\rightarrow 0\) as well as in the case illustrated in Fig. 9.

5 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 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 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 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.