DYNAMIC GRIDS FOR FINITE-DIFFERENCE SCHEMES IN MUSICAL INSTRUMENT SIMULATIONS

For physical modelling sound synthesis, many techniques are available; time-stepping methods (e.g., finite-difference time-domain (FDTD) methods) have an advantage of flexibility and generality in terms of the type of systems they can model. These methods do, however, lack the capability of easily handling smooth parameter changes while retaining optimal simulation quality and stability, something other techniques are better suited for. In this paper, we propose an efficient method to smoothly add and remove grid points from a FDTD simulation under sub-audio rate parameter variations. This allows for dynamic parameter changes in physical models of musical instruments. An instrument such as the trombone can now be modelled using FDTD methods, as well as physically impossible instruments where parameters such as e.g. material density or its geometry can be made time-varying. Re-sults show that the method does not produce (visible) artifacts and stability analysis is ongoing.


INTRODUCTION
The operation of most musical instruments can be subdivided into excitation and resonator components [1]. Examples of excitationresonator combinations are the bow and violin and the lips and trumpet. In most instruments, the parameters describing the excitation are continuously varied by the performer to play the instrument. As an example, the bow velocity, bow position and bow force for stringed instruments, and lip pressure and frequency for brass instruments. Naturally, the resonator is also altered by fingering the strings of the violin or pressing valves on the trumpet to change the instrument pitch. But, even under such variable playing conditions, physical properties of the resonators do not change: the string length and tension stay the same and the total tube length remains unchanged; it is only the portions that resonate that are shortened or lengthened.
There are several examples where the parameters of the resonator are also modified. A prime example of this is the trombone, where the tube length is dynamically changed in order to generate different pitches. The slide whistle is another example in this category. Guitar strings are another category where the tension can be smoothly modulated during performance using the fretting finger or a whammy bar to create smooth pitch glides. The same kind of tension modulation is used for the membranes of timpani or "hourglass drums" to change the pitch. It is these direct parameter mod- ifications of the resonators that we are interested in simulating. In addition to simulating existing instruments, one could potentially simulate instruments that can be manipulated in physically impossible ways. Examples of this could be to dynamically change material properties such as density or stiffness, or even the geometry and size of the instrument where this is physically impossible.
Finite-difference time-domain (FDTD) methods are flexible and generalisable techniques which have seen increased use in physical modelling sound synthesis applications [2]. The normal approach, for a given system such as a musical instrument, is to describe its motion by a set of partial differential equations (PDEs). The instrument is then represented over a spatial grid, and a timestepping method is developed, yielding a fully discrete approximation to the target PDE system. In many cases, the system itself is static, so that the defining parameters do not change over time. In others, such as the trombone and others mentioned above, this is not the case, and various technical challenges arise when trying to design a simulation using FDTD methods; all relate to the choice of the spatial grid. For example, the grid density is usually closely tied to the parameters themselves through a stability condition. Also, adding and removing points from the grid is nontrivial and can cause audible artifacts and new stability concerns. The default approach of defining a grid globally, according to a very conservative stability condition, as done in [3], is possible, but introduces numerical dispersion and bandlimiting effects. Full-grid interpolation [2,Ch. 5] could be used to change between grid configurations, but extremely high sample rates are necessary to avoid audible artifacts and low-passing effects, rendering any implementation offline.
In this paper, a new method is proposed, allowing the efficient and smooth insertion and deletion of grid points from 1D finitedifference grids to allow for dynamic parameter changes. We are interested in varying parameters 'slowly' (i.e., at sub-audio rate corresponding to human gestural control). In a companion paper we present a physical model of the trombone using the method proposed in this paper [4]. Notice that other techniques do allow for dynamic parameter changes but come with their own drawbacks [2]. Examples of dynamic parameters using modal synthesis [5] are shown in [6,7] and digital waveguides [8] are shown in [9]. This paper is structured as follows: Section 2 presents the 1D wave equation, to be used as an illustrative example for the proposed method. Section 3 gives an introduction to numerical methods, stability and simulation quality. The proposed method for dynamic grids is then presented in Section 4 and applied to the 1D wave equation. Section 5 shows the results of an analysis performed on the method, which are discussed in Section 6. Finally, concluding remarks and future perspectives are given in 7.

CONTINUOUS SYSTEMS
The wave equation is a useful starting point for investigations of time-varying behaviour in musical instruments. In 1D, the wave equation may be written as and is defined over spatial domain x ∈ [0, L], for length L (in m) and time t ≥ 0 (in s). c (in m/s) is the wave speed. The dependent variable q = q(x, t) in Eq. (1) may be interpreted as the transverse displacement of an ideal string, or the acoustic pressure in the case of a cylindrical tube. Two possible choices of boundary conditions are and describe 'fixed' or 'free' boundary respectively in the case of an ideal string, and 'open' or 'closed' conditions respectively in the case of a cylindrical acoustic tube.

Dynamic parameters
In the case of the 1D wave equation, only the wave speed c and length L can be altered (in the case of an acoustic tube, only L is variable, and for a string, c could exhibit variations through changes in tension). If the same boundary condition is used at both ends of the domain, and under static conditions, the fundamental frequency f0 of vibration can be calculated according to In the dynamic case, and under slow (sub-audio rate) variations of c or L, Eq. (3) still holds approximately. From Eq. (3), one can easily conclude that in terms of fundamental frequency, halving the length in Eq. (1) is identical to doubling the wave speed and vice versa. Looking at Eq. (1) in isolation, f0 is the only behaviour that can be changed. One can thus leave L fixed and allow time variation in c, so that c = c(t), which will prove easier to work with in the following sections. This fact can more easily be seen if Eq. (1) is scaled or non-dimensionalised as in [2], where scaled domain x ′ = x/L ⇒ x ′ ∈ [0, 1] and γ = c/L such that f0 = γ/2. For clarity, however, we will employ a fully dimensional representation here.

NUMERICAL METHODS
This section will provide a brief introduction to physical modelling using FDTD methods, including details on stability and quality of the simulations based on these methods. In this section, c is assumed constant.

Discretisation
In FDTD methods, the first step is the definition of a grid. The spatial variable can be discretised using x l = lh with integer l ∈ {0, . . . , N }. The grid spacing h (in m) is the distance between adjacent grid points, and the total number of points covering the domain, including endpoints, is N + 1. Here, integer N describes the total number of intervals between the grid points, and thus the total domain length is L = N h. The temporal variable can be discretised using tn = nk with positive integer n, time step k = 1/fs (in s) for sample rate fs (in Hz). The state variable q can then be approximated using q n l ≊ q(x = lh, t = nk).
The following operators can then be applied to q n l to get the following approximations to the derivatives in Eq. (1) Substituting these definitions into Eq. (1) yields the following finite-difference (FD) scheme δttq n l = c 2 δxxq n l .
Expanding the operators as in (4) and solving for q n+1 l yields the following update equation which is suitable for direct software implementation. Here, is referred to as the Courant number, constrained by numerical stability conditions, and also has an impact on the quality and behaviour of the simulation. This will be described in detail in Sections 3.2 and 3.3.
In the FD scheme described in Eq. (5), the boundary locations are at l = 0 and l = N . Substituting these locations into Eq. (6) seemingly introduces the need of grid points outside of the defined domain, namely q n −1 and q n N +1 . These can be referred to as virtual grid points and can be accounted for using the boundary conditions in Eq. (2). Discretising these yields is a second-order accurate approximation of the first-order spatial derivative. The Dirichlet condition in (8a) says that the displacements of q at the boundary locations are always 0. In practice, this means that these grid points do not need to be updated and the spatial range of calculation for Eq. (6) then becomes l ∈ {1, . . . , N − 1}. If the Neumann condition is used, the boundary points do need to be updated as these are not necessarily 0; rather, their 'slope' is 0. Eq. (8b) can then be expanded to yield defnitions for these virtual grid points q n −1 = q n 1 and q n N +1 = q n N −1 .
Now that the full system is described, audio output at sample rate fs can be drawn from the state q n l in Eq. (6) at 0 < l < N (when using Dirichlet boundary conditions).

Stability
Explicit FDTD methods for hyperbolic systems such as the 1D wave equation must necessarily satisfy a stability condition. In the case of the update in Eq. (6) it can be shown -using von Neumann analysis [10] -that the system is stable if which is referred to as the Courant-Friedrichs-Lewy (CFL) condition. The more closely λ approaches this condition with equality, the higher the quality of the simulation (see Section 3.3) and if λ = 1, Eq. (6) yields an exact solution to Eq. (1). If λ > 1 the system will become unstable. Recalling (7), Eq. (11) can be rewritten in terms of grid spacing h to get This shows that the CFL condition in (11) puts a lower bound on the grid spacing, determined by the sample rate and wave speed. Usually, the following steps are taken to calculate λ: where ⌊·⌋ denotes the flooring operation. In other words, condition (12) is first satisfied with equality and used to calculate number of intervals N . Thereafter, h is recalculated based on integer N and used to calculate λ. The calculation of λ in Eq. (13) can be compactly rewritten as The flooring operation causes the CFL condition in (11) to not always be satisfied with equality and results in a reduced simulation quality described in the following section.

Simulation Quality
Choosing λ < 1 in Eq. (6) will decrease the simulation quality in two ways. Firstly, it will decrease the maximum frequency that the simulation is able to produce, i.e., it will decrease the bandwidth of the output sound of the system. By analysing the scheme in Eq. (6), it can be shown that the maximum frequency produced by the system can be calculated using fmax = fs sin −1 (λ)/π [2, Chap. 6]. Note that only a small deviation of λ from condition (11) leads to a large reduction in output bandwidth. Secondly, choosing λ < 1 causes numerical dispersion. Harmonic partials become unnaturally closely spaced at higher frequencies (i.e. spurious inharmonicity increases) as λ decreases, which is generally undesirable.

THE DYNAMIC GRID
The time variation of the wave speed c leads to various complications in the simulation framework presented above. First of all, a change in c causes a change in λ according to Eq. (14), affecting the simulation quality and bandwidth. Secondly, and more importantly, a change in c could result in a change in N through Eq. (13). As N directly relates to the number of grid points, this raises questions as to where and especially how one would add and remove points to the grid according to the now-dynamic wave speed.
We propose a method that allows for a non-integer number of intervals to smoothly change between grid configurations, i.e, the number of grid points used. This removes the necessity of the flooring operation in Eqs. (13) and (14), and consequently satisfies the CFL condition in (11) with equality at all times. Introducing fractional number of intervals N , where N = ⌊N ⌋, Eq. (3) can be rewritten in terms of N by substituting the calculation of N from This shows that if λ = 1, N solely determines the fundamental frequency of the simulation. Ideally, a method that dynamically changes the grid size of a FD scheme should r1. generate an output with a fundamental frequency f0 which is proportional to the wave speed c (f0 ∝ c), r2. allow for a fractional number of intervals N to smoothly (without audible artifacts) transition between different grid configurations, r3. generate an output containing N − 1 modes which are integer multiples of the fundamental (fp = f0p with integer p), r4. work in real time to have a playable simulation.
These requirements will be used in Section 6 to evaluate the proposed method.

Proposed Method
In the following, the location of a grid point (in m from the left boundary) q l at time index n is denoted by x n q l . Furthermore, some variables are now time dependent as indicated by superscript n. These are c n , h n , N n , N n and f n 0 .

System Setup
Consider two grid functions, u n lu and w n lw defined over discrete domains lu ∈ {0, . . . , M n } and lw ∈ {0, . . . , M n w } respectively with integers M n = ⌈0.5N n ⌉ with ⌈·⌉ denoting the ceiling operation and M n w = ⌊0.5N n ⌋, i.e., half the number of points allowed by the stability condition, plus one for overlap. The two grid functions are assumed to lie adjacent to each other on the same domain x. For now, the grid locations lu = M n and lw = 0 are assumed to overlap so that x n u M n = x n w 0 = M n h n , and are referred to as the inner boundaries. The grid locations lu = 0 and lw = M n w are placed at x n u 0 = 0 and x n w M n w = L and will be referred to as the outer boundaries. See Figure 1a. The following boundary conditions are then imposed: In other words, grid points at the outer boundaries are fixed, according to the usual Dirichlet condition, and those at the inner boundaries are free. It is important to note that the Neumann condition is just used as a starting point for the method here, but will be modified in Section 4.1.2. The systems can then be connected at the inner boundaries using a rigid connection DAFx.3 Proceedings of the 24 th International Conference on Digital Audio Effects (DAFx20in21), Vienna, Austria, September 8-10, 2021 Notice that this condition only needs to be satisfied when the inner boundaries perfectly overlap, which is not always the case when c n is varied (see Section 4.1.2).
To sum up, a grid function with N intervals as per Eq. (13) is divided into two separate subsystems connected at their respective inner boundaries.
With the above boundary conditions imposed, the following state vectors can be defined: with T denoting the transpose operation, and have M n and M n w points respectively. Note that the grid points at the outer boundaries are excluded as they are 0 at all times due to the Dirichlet boundary condition in (8a). A vector concatenating (18) is then defined as Even though the new system has an extra (overlapping) grid point, the behaviour of the new system should be identical to that of the original system in Eq. (5) with (static) N n = N n . That this holds will be shown below.
Using u n lu and w n lw in the context of the 1D wave equation, a system of FD schemes can be defined as with spreading operators Ju(x n i ) = applying the effect of the connection F n (in m 2 /s 2 ) to grid points u n M n and w n 0 respectively. Expanding the spatial operators in system (20) at the inner boundaries, recalling the Neumann condition in (16b) and the definition for the virtual grid points needed for this condition in Eq. (10) yields It is important to note that the time index n in M n will not be affected by the δtt operator and all obtained terms after expansion (Eq. (4a)) will use the same value for M n . Because of the rigid connection in (17), it is also true that δttu n M n = δttw n 0 (if x n u M n = x n w 0 ), and F n can be calculated by setting the right side of the equations in (22) equal to each other: Substituting this into system (22) after expansion of the second-

Changing the Grid
The previous section describes the case in which N n is an integer. We now continue by varying c n such that this is not the case.
The locations of the outer boundaries x n u 0 and x n w Mw are fixed: If the wave speed c n is then decreased, and consequently the grid spacing h n according to Eq. (12) (with equality), all other points move towards their respective outer boundary (see Figure 1b). Calculating h n this way allows this method to always satisfy the CFL condition in Eq. (11) with equality, solving issues regarding simulation quality and numerical dispersion described in Section 3.3.

DAFx.4
Proceedings of the 24 th International Conference on Digital Audio Effects (DAFx20in21), Vienna, Austria, September 8-10, 2021 As mentioned in Section 4.1.1, the state of the virtual grid points at the inner boundaries are defined as u n M n +1 = w n 1 and w n −1 = u n M n −1 when the inner boundaries perfectly overlap (i.e., x n u M n = x n w 0 ). If this is not the case (x n u M n ̸ = x n w 0 ) a Lagrangian interpolator I(x n i ) at location x n i (in m from the left boundary) can be used to calculate the value of these virtual grid points (also see Figure 1c for reference). The interpolator I is a row-vector with the same length as U n (from Eq. (19)) and its values depend on the interpolation order. In the following, the fractional part of N n is defined as and for clarity, I and U n are indexed by m. Now, consider the following quadratic interpolator and its flipped version where the shift in the latter is necessary to transform the location x n i to the correct indices of U n . When applied to Eq. (19) this yields the definitions for the virtual grid points These definitions for the virtual grid points at the inner boundaries will replace the Neumann condition in Eq. (16b). One can show that when N n is an integer, and thus α = 0, Eqs. (26a) and (26b) can be substituted as w n 1 and u n M n −1 into Eqs. (23a) and (23b) respectively (as these acted as virtual grid points u n M n +1 and w n −1 ). Then recalling Eq. (17) it can be seen that the system reduces to (23) and exhibits the same exact behaviour as the usual case in Eq. (5). Now that the virtual grid points at the inner boundaries are not determined by the Neumann boundary condition in (16b), but rather by the definitions in Eqs. (26), system (20) can simply be re-written to δttu n lu = (c n ) 2 δxxu n lu , δttw n lw = (c n ) 2 δxxw n lw , where the Dirichlet condition in (16a) is (still) used for the outer boundaries and the Neumann condition at the inner boundaries in (16b) is replaced by the definitions in (26): u n M n +1 = I ← 2 (x n u M n +1 )U n and w n −1 = I2(x n w −1 )U n . (28) Figure 2: The moment when a point is added to u at location x n u M n +1 in Eq. (29). This figure shows an extreme case where this location is far from x n w 0 , i.e., α ̸ ≈ 0 in Eq. (30).

Adding and Removing Grid Points
When c n , and consequently h n , are decreased and the inner boundary points surpass the virtual points (i.e. x n u M n ≤ x n w −1 and x n w 0 ≥ x n u M n +1 ), this means that N n > N n−1 . A point is then added to the right boundary of u and the left boundary of w (for both time indices n and n − 1) in an alternating fashion: and cubic Lagrangian interpolator with I ← 3 being a flipped, not shifted (as I ← 2 in Eq. (25b)) version of (30). See Figure 2. Notice that N n is only going to be slightly bigger than an integer at the moment that a point is added and Eq. (24) will return α ≳ 0. This means that that I3 ≈ [0, 0, 1, 0] and the displacement of the newly added point is nearly fully based on the grid point at the inner boundary of the other system.
Removing grid points happens when c n , and consequently h n , are increased and x n u M n > x n w 0 (or N n < N n−1 ). Grid points are simply removed from u and w (again for both n and n − 1) in an alternating fashion according to In Eqs. (29) and (31), the even and odd conditions can be inverted. To keep the difference between u and w a maximum of one grid point, the ceiling and flooring operations when calculating M n and M n w will need to be inverted as well. Until now, only adding and removing points in the center of the original system has been considered. This location could be moved anywhere along the grid, the limit being one point from the boundary. In other words, both u n lu and w n lw need to have at least one point (excluding the grid points at the outer boundaries). Furthermore, one does not have to add and remove points from u and w in an alternating fashion as in (29), but can just add to and remove from, for example, u leaving w the same size throughout the simulation. In the extreme case where M n = N n − 1 and M n w = 1 (leaving w n lw with only one moving grid point, w n 0 ) the method still works.

Displacement correction
A problem that arises when increasing c n , is that it is possible that the displacements u n M n ̸ ≈ w n 0 at the time when a grid point needs to be removed. As the grid locations x n u M n ≈ x n w 0 at the time of removal, this violates the rigid connection in (17) and causes audible artifacts. A method is proposed that decreases the relative displacement of the inner boundaries the closer their grid-locations are together, i.e., the closer α in (24) is to 0. We thus extend system (27) with an artificial spring force as Using centred temporal averaging and difference operators the correction effect is defined as with the difference in displacement between the inner boundaries and damping coefficient σ0. Furthermore, β scales the effect of the displacement correction and is defined as where ε ≪ 1 prevents a division by 0. Despite the operators in (33) introducing states at n + 1, it is possible to calculate the force explicitly (such as in [2] or [11]). Furthermore, it can be shown that even when ε = 0 this calculation is always defined. In that case, as α → 0, β → ∞ which acts as a rigid connection such as Eq. (17). Essentially, the displacement correction attempts to have η n → 0 in Eq. (35) as α → 0 to satisfy the rigid connection in Eq. (17). Although the correction presented here is not based on some physical process, it can be justified by the fact that large differences in displacement between two spatially adjacent points is not physical. Notice that when c n is decreased, the rigid connection will not be violated as u n M n ≈ w n 0 when a point is added. This is due to the fact that I3 ≈ [0, 0, 1, 0] and either u n M n or w n 0 is the newly added point which almost solely based on the other.

Summary
Here, Section 4.1 is summarised and describes the final version of the proposed method.
The proposed method subdivides a grid function q n l with N intervals into two grid functions u n lu and w n lw with M n and M n w intervals respectively for a total of N n + 2 grid points. Knowing that λ = 1 ∀n, Eq. (6), written for both grid functions, becomes w n+1 lw = w n lw +1 + w n lw −1 − w n−1 lw .
Due to the Dirichlet boundary condition in (16a) imposed at the outer boundaries of the system, u n 0 and w n Mw are 0 at all times and do not have to be included in the calculation. The ranges of calculation for Eq. (37a) and (37b) then become lu ∈ {1, . . . , M n } and lw ∈ {0, . . . , M n w − 1} respectively. The grid points at the inner boundaries are calculated by expanding (27) (ignoring the displacement correction for now) where virtual grid points u n M n +1 and w n −1 can be calculated using Eq. (26).
Then, when N n > N n−1 , a point is added to u n and u n−1 (or w n and w n−1 ) using Eq. (29), and when N n < N n−1 , a point is removed from the same vectors using Eq. (31). In order to prevent audible artifacts when increasing c n (and thus decreasing N n ) due to a violation of the rigid connection in (17), a method is proposed in Eq. (32) to ensure that the grid points at the inner boundaries have a similar displacement when one of them is removed.

Implementation
A MATLAB implementation of the proposed method and audio examples can be found online 1 and Algorithm 1 shows the order of calculation of this implementation. Especially important to take into account, is to only retrieve a change in c n at time index n once before all other calculations. This is to ensure that u n lu and w n lw are calculated with the same α and β for all lu and lw.

Modes
Writing system (32) in matrix form, one can perform a modal analysis while changing c n to obtain the frequencies and damping coefficients for each mode over time. As a test case, the wave speed of a system running at fs = 44100 Hz is linearly varied from c 0 = 2940 (N 0 = 15) to c n end = 2205 (N n end = 20) where nend = tdurfs is the simulation length in samples and tdur = 10 s. Grid points are added and removed as close to the right boundary as possible, i.e., M n = N n − 1 and M n w = 1 (similar behaviour can be observed if M n = 1 and M n w = N n − 1). The result of the analysis is shown in Figure 3a where higher damping (induced by the displacement correction) is indicated using thinner and bluer lines. Figure 3b shows the resulting spectrogram, with the displacement correction deactivated, of the system excited with u 0 1 = 1 and the output retrieved at u n 1 , and Figure 3c shows a system with the same excitation but the change in c n inverted (N n = 20 → 15) and displacement correction activated.
In the following, the lowest mode generated by the analysis is referred to as f n 1 and should ideally be equal to f n 0 calculated using Eq. (3). The first thing one can observe from Figure 3a is that the frequencies of the modes decrease as c n decreases (as desired). The lower the mode, the more linear this decrease happens. Between N n = 15 and N n = 16, f n 1 maximally deviates by −0.15 cents. In this same interval f n 15 maximally deviates by −67 cents. This deviation gets less as N n increases. Experiments with higher even-ordered Lagrange interpolators show that these frequency deviations become smaller, but not by a substantial amount. The quadratic interpolator has thus been chosen for being simpler and more flexible while not being substantially worse than higher order interpolators.
Another observation from Figure 3a is that there are always N n modes present, corresponding to the number of moving points of the system. As can be seen in Figure 3b the highest mode is not excited. If the system is excited when N n is not an integer, the highest mode will also be excited. Comparing the implementation of the system using this method with integer N n (without changing c n ) to a normal implementation of the 1D wave equation (shown in Section 3) with (static) N n = N n , identical outputs are observed, even though the latter has N n − 1 moving points.

Displacement Correction
In the experiments, σ0 = 1 in Eq. (34). The displacement correction has a low-pass-comb-filtering effect on the system, where the position and amount of damped regions directly relates to the position of where grid points are added and removed. The best behaviour, i.e., least affecting lower frequencies, is when grid points are added and removed as close to the boundary as possible, i.e., M n = N n − 1, and only has one damped region as shown in Figures 3a and 3c.

Limit on Rate of Change of c
The current implementation of the proposed method can only add or remove a maximum one point per sample using Eqs. (29) and (31). The rate of change of f n 0 according to (15) is thus limited by |N n − N n−1 | ≤ 1. Though this is the maximum limitation on speed, a much lower limitation needs to be placed to keep the system well-behaved. The usual stability and energy analyses performed on FD schemes are not valid anymore in the time-varying case. Frozen coefficient analysis as in [10] could be applied here and hold for slowly varying coefficients, but is left for future work.  Figure 3a along the x-axis and applying this to a 10 s simulation).

DISCUSSION
To decide whether the proposed method is satisfactory, the results presented in the previous section are compared to the method requirements listed in Section 4.
It can be argued that the frequency deviations of f n 1 from f n 0 are sufficiently small to say that r1 is satisfied. As for r2, a fractional number of intervals N n has been introduced and smooth transitions are indeed observed from Figure 3b, in the case when c n is decreased and N n is increased. When c n is increased instead, the displacement correction prevents (visible) artifacts when DAFx.7 Proceedings of the 24 th International Conference on Digital Audio Effects (DAFx20in21), Vienna, Austria, September 8-10, 2021 grid points are removed as seen in Figure 3c. Despite this, the filtering effect that the displacement correction has on the system (mentioned in Section 5.2) is not ideal as it creates damped regions in the spectrum of the output sound. The least intrusive filtering happens when points are added and removed as close to the boundary as possible, i.e., when M n = 1 or M n w = 1 where the damping only occurs in the higher end of the spectrum. Although artifacts do not show in Figure 3c, to confirm the absence of audible artifacts, formal listening tests have to be carried out. Furthermore, higher speeds of parameter variation might cause artifacts anyway. The value of σ0 could therefore also be made dynamic and depending on the rate of change of c n to have a higher effect when c n is increased faster and vice versa. Either way, as this is still not ideal, another method for reducing artifacts that less affects the frequency content of the system should be devised, if possible. Furthermore, higher modes will be lost after decreasing N and will not return after increasing N again. They can, however, be activated again by re-exciting the system.
The modal analysis in Figure 3a shows that the method generates N n rather than N n − 1 modes as set by r3. However, the output does contain the correct number of modes as shown in Figure 3b due to the highest mode not being excited. This is a result of the rigid connection imposed on the inner boundaries, forcing them to have the same displacement and act as one point. The latter part of r3, however, is not satisfied. The modes deviate from integer multiples of f n 0 , moreso for higher modes. Other interpolation techniques could be investigated to improve the behaviour and decrease this deviation.
Finally, the method only adds a few extra calculations for the inner boundaries so r4 is also easily satisfied.
Although the results bring forward some drawbacks of the proposed method, such as modal frequency deviations, and filtering effects, most of these affect the higher frequencies of the output. First of all, human frequency sensitivity becomes very limited above 3000 Hz [12] making high-frequency deviations much less important perceptually. Secondly, the physical systems one usually tries to model contain high-frequency losses, causing higher modes to usually not have very high amplitudes to begin with. Finally, N n is usually much bigger in the systems one models, where frequency deviations happen to a much smaller degree.
As of now, some aspects of the proposed method still lack physical justification (such as the displacement correction), but are shown to yield the desired behaviour and fulfil the requirements to a satisfactory degree. Despite this, further work needs to be done to physically justify the choices made in this paper.

CONCLUSIONS AND PERSPECTIVES
This paper presents a method to change grid configurations of finite-difference schemes to allow for dynamic parameter changes. The method allows the stability condition that these schemes rely on can be satisfied with equality at all times, minimising numerical dispersion and bandlimiting issues. Grid points are shown to be added and removed smoothly and do not cause artifacts when switching between grid configurations. Listening tests will need to be performed to carried out to confirm the lack of audible artifacts.
The proposed method might not provide an exact solution to the problem of time-varying systems, and not all choices are physically justified, but it does circumvent the need for upsampling and higher orders of computations necessary to approximate this solution with, for example, full-grid interpolation.
Although this method has only been applied to the 1D wave equation it could be applied to many other 1D systems. Other parameters, such as material density or stiffness could also be made dynamic, going beyond what is physically possible. An application of the method that could be investigated is that of non-linear systems, such as the Kirchhoff-Carrier string model [13] where the tension is modulated based on the state of the system.
Other future work includes creating an adaptive version of the displacement correction that changes its effect depending on the speed at which the grid is changed. Finally, stability and energy analyses will have to be performed to show the limits on changes in parameters and grid configurations.

ACKNOWLEDGMENTS
This work is supported by NordForsk's Nordic University Hub Nordic Sound and Music Computing Network NordicSMC, project number 86892.