Non-iterative, conservative schemes for geometrically exact nonlinear string
vibration
Michele DUCCESCHI(1) , Stefan BILBAO(1)(1) Acoutics and Audio Group, University of Edinburgh, UK, michele.ducceschi@ed.ac.uk

Abstract
This work presents novel finite difference schemes for the solution of the geometrically exact nonlinear model ofwave propagation in strings. Stability is enforced by ensuring conservation in the discrete setting. The schemes
are here derived from a suitable quadratisation of the nonlinear potential, yielding ultimately a non-iterative
update equation requiring the solution of a sparse linear system per time step. Comparison with previous iterative
schemes is carried out, showing convergence to a unique solution. The new schemes are extremely efficient,yielding speedups of about half an order of magnitude over previously available schemes.
Keywords: Finite Difference Schemes, Energy Methods, Nonlinear Wave Equation, Piano Acoustics

1

INTRODUCTION

Computer simulation of nonlinear phenomena is a subject of growing interest in acoustics and musical acoustics.
Though linear models are adequate in many situations, exceptions are many: in one dimension, a nonlinear wave
equation is required to fully capture the behaviour of piano strings, for instance. When it comes to simulation,
the algorithm designer is faced with two main concerns: one is ensuring convergence of the schemes; the otheris efficiency. Among the usual approaches in numerical simulation for nonlinear systems are iterative schemes,
obtained via discretisation methods such as the bilinear transform (i.e. the trapezoid rule.) These techniqueslead to stable and convergent schemes, but can be quite inefficient, especially when solving for distributed
nonlinearities such as those occurring in a piano string. Recently, a new class of methods has been proposed,
relying on a suitable quadratisation of the nonlinear potential energy. Such methods appeared first in virtualanalog simulation (1, 2, 3), and have been extended to include the case of non-invertible potentials (4, 5, 6), andhere to the case of a fully distributed nonlinearity. Applications in fluid dynamics have also been devised (7).
The proposed schemes ultimately lead to update equations which are resolvable without recourse to iterative
methods, and which maintain the notion of an energy balance leading to a numerical stability guarantee. In
Section 2 the continuous models of nonlinear string vibration are presented, along with considerations regarding
the quadratisation of the potential, and a linear expansion. In Section 3 two conservative finite difference
schemes are given: the first, iterative in nature, follows from a discretisation obtained via the trapezoid rule;the second is the newly proposed non-iterative scheme. Finally, in Section 4 some numerical experiments aredetailed.

2

MODEL EQUATIONS

In musical acoustics, it is often assumed that a uniform string vibrates according to the simple wave equation.
There are two cases when this assumption is no longer valid: when the string has intrinsic stiffness, and when
it vibrates at higher amplitudes. Stiffness is indeed a prominent effect across all musical strings, and it canbe incorporated by means of an appropriate linear beam model. From a numerical standpoint, it has beenstudied in numerous other works [see, e.g. (8, 9)], so here the focus will be exclusively on the second effect,
nonlinearities. In the case of a string vibrating in isolation, these are of geometric type and arise when thevibration amplitude is large. It is in fact possible to derive a model which is geometrically exact, by calculating
the extension of an infinitesimal string element under tension (10). Ultimately, the stretching of the string leads

5507

to tension modulation (as opposed to the linear case, where tension is constant.) In vector form, this model is

q


EA 2
(1 + p)2 + q2 − 1 − p ≥ 0 ∀ (q, p) (1)ρA∂t2 w = ∂x ∇(q,p) Ξ , with Ξ(q, p) =q + p2 + (EA − T0 )2
where Ξ is interpreted as a positive-definite potential function. Above, ρ is the volumetric density, and A is thearea of the cross section, E is Young’s modulus and T0 is the applied tension. The string occupies a domainx ∈ D = [0, L]. The symbol ∂sn denotes the nth partial derivative along s, and ∇(q,p) is the gradient taken along
q and p. Assuming the motion to be constrained in a single plane, one can identify w(t, x) = T [u(t, x), ζ (t, x)],where u is the transverse (or vertical) displacement, and ζ is the longitudinal (or horizontal) displacement.
Hence, q , ∂x u and p , ∂x ζ . It is convenient to rewrite (1) in a different, yet equivalent form, as
ρA∂t2 u = T0 ∂x2 u + ∂x (∂q φ )(2a)ρA∂t2 ζ = T0 ∂x2 ζ + ∂x (∂ p φ )ζ

(2b)

where

q
2
EA − T022
φ=(1 + p) + q − 1
2
Notice that φ is positive definite, if EA > T0 (certainly the case for all musical strings.)(3)

2.1 Inner product and energy analysis in the continuous case
An energy analysis of (2) is readily available upon the introduction of an inner product and associated norm.
For two well-behaved functions c, d, one has
Z L
q
(4)
Inner Product: hc, diD =c d dx,Norm: kckD = hc, ciD
0Taking an inner product of (2a) with ∂t u, and of (2b) with ∂t ζ gives (after integration by parts)
dH
=B
dt

(5)

with H being the total energy, and with B being the boundary conditions. The analytic expressions are
 T 

ρA 
0
k∂t uk2D + k∂t ζ k2D +kqk2D + kpk2D + < φ , 1 >DH=
2
2
B = [∂t u (T0 q + ∂q φ ) + ∂t ζ (T0 p + ∂ p φ )]L0

(6a)
(6b)

A convenient set of boundary conditions which make B vanish, used in this work, is thenu = ζ = 0 at x = 0, L

(7)

2.2 Quadratisation of the nonlinear potential
A convenient alternate form of the equations of motion can be derived by quadratising the potential φ . Hence,takepψ = 2φ(8)
One hasρA∂t2 u = T0 ∂x2 u + ∂x (ψ∂q ψ)(9a)ρA∂t2 ζ = T0 ∂x2 ζ + ∂x (ψ∂ p ψ)

(9b)

The energy analysis in this case yields an energy balance analogous to (5), where now
 T 
 kψk2ρA 
0Dk∂t uk2D + k∂t ζ k2D +kqk2D + kpk2D +222
B = [∂t u (T0 q + ψ∂q ψ) + ∂t ζ (T0 p + ψ∂ p ψ)]L0
H=

5508

(10a)
(10b)

Though (9) and (2) are completely equivalent (so long as ψ is well defined, i.e. so long as, from (8), φ is nonnegative), the two forms lead to very different discrete methods: system (2) will require a nonlinear iterativeroot finder, whereas system (9) will be calculated by means of a simple matrix inversion.
2.3 Small oscillations and natural wave speeds
One may wonder what the behaviour of the nonlinear wave equation (2) is under small amplitudes. An expansion of φ to second order yields
EA − T0 2
φ≈
p
(11)
2
Inserting this value in (2) gives
ρA∂t2 u = T0 ∂x2 u

ρA∂t2 ζ = EA∂x2 ζ

(12a)

(12b)

Hence, for small amplitudes, the system is completely uncoupled, andpboth transverse andp longitudinal motions
are given by the linear wave equation with speeds, respectively, vu = T0 /ρA and vζ = E/ρ. Notice that formusical strings, normally vζ  vu , and this will have consequences in terms of the choice of the grid spacingin the finite difference schemes.

3

DISCRETE MODEL

Solutions to the nonlinear wave equation are sought by means of the finite difference method. The continuousfunctions u(t, x), ζ (t, x) are approximated by grid functions unm , u(nk, mh), ζmn , ζ (nk, mh), where the indicesm, n are positive integers, and where k is the time step, and h is the grid spacing. It is assumed that n ≥ 0
and m ∈ M = [0, M], and thus the grid spacing h divides the domain length L in M equal intervals. The
continuous functions p(t, x), q(t, x) are approximated by pnm−1/2 , qnm−1/2 , defined on an interleaved spatial gridm ∈ M = [1, M]. In the following, in order to keep the index notation to a minimum, it will be assumed that thenotation u stands for unm , and similarly for ζ , p and q (these last two defined on an interleaved grid.) Indiceswill be specified only when needed.
Finite difference operators are now introduced. The identity and time shifting operators are
1wnm = wnm ,n+1
et+ wnm = wm

et− wnm = wn−1
m

(13)

δt· ,

et+ − et−
,
2k

(14)

1 + et−,2µt· ,Time difference operators may be introduced as
δt+ ,et+ − 1,kδt− ,1 − et−,k

1 + et+,2µt− ,δtt , δt+ δt−

Time averaging operators areµt+ ,et+ + et−2

(15)Space shifting operators are defined as
ex+ wnm = wnm+1 ,ex− wnm = wnm−1(16)

Similar definitions hold for spacial operators acting on an interleaved grid function. From those, space difference
operators are given as
ex+ − 11 − ex−
δx+ ,, δx− ,
, δxx , δx+ δx−(17)hhPartial differentiation of the potential φ with respect to p, q is given as
δq· φ (qn , pn ) ,φ (qn+1 , pn ) − φ (qn−1 , pn ),qn+1 − qn−1

δ p· φ (qn , pn ) ,φ (qn , pn+1 ) − φ (qn , pn−1 )pn+1 − pn−1

(18)

5509

A particular averaging of φ will also be needed. This isµq− |p− φ (qn , pn ) ,φ (qn−1 , pn ) + φ (qn , pn−1 )2(19)

For two grid functions cm , dm defined over D = d− ≤ m ≤ d+ one can define an inner product and associatednorm, in analogy to (4), asd+

Inner Product: hc, diD =∑ h cm dm ,qhc, ciD(20)

ρAδtt ζ = T0 δxx ζ + δx+ (δ p· φ )(21b)

Norm: kckD =m=d−

3.1 Iterative, conservative scheme
A discretisation of (2) can be written as
ρAδtt u = T0 δxx u + δx+ (δq· φ )(21a)where q , δx− u, p , δx− ζ . This discretisation follows from applying the trapezoid rule to the nonlinear potential.
An energy conservation property can be found by taking an inner product of (21a) with δt· u, and of (21b) with
δt· ζ , and by summing the resulting equations. Applying summation by parts gives
δt+ hn−1/2 = b

(22)

where
hn−1/2 =

 T 

ρA 
0
kδt− un k2M + kδt− ζ n k2M +hqn , et− qn iM + hpn , et− pn iM + 1, µq− |p− φ n M
2
2

(23)

and where the boundary term b vanishes under a choice of fixed end conditions, i.e. u0 = uM = ζ0 = ζM = 0 ∀n.
Notice that (22) is discrete version of (5). However, as opposed to the continuous energy H which is always
non-negative, the discrete energy h in (23) will be non-negative if and only if (11)sT0h≥k , h(t)(24)ρA
Non-negativity of the energy function assures stability and convergence for scheme (21).
Remark. Though condition (24) leads to stable simulations, choosing h to be close to h(t) would lead in
this case to some aliasing problems. pRemembering the linear expansion (12b), the natural grid spacing forthe longitudinal grid is indeed h(l) = E/ρk, which is much larger than h(t) . In order to avoid such spacialsampling problems, in the following the choice h = h(l) will be enforced, along with frequency oversampling.
More efficient schemes, avoiding oversampling, may be devised by interpolating between the longitudinal andtransverse grids. However, this possibility will not be explored here, as the focus is only on the convergenceand stability properties of the schemes 
Both (21a) and (21b) can conveniently be multiplied on the left by δx− . This yields the following two equations
ρAδtt q = T0 δxx q + δxx (δq· φ )(25a)ρAδtt p = T0 δxx p + δxx (δ p· φ )(25b)Then, definer(q) = qn+1 − qn−1 ,5510

r(p) = pn+1 − pn−1 ,a(q) = qn−1 ,a(p) = pn−1(26)

Owing to these definitions, one may recast system (25) in the following form!

φ r(q) + a(q) , pn − φ a(q) , pn
ρA
n
n−1
n
r − 2q + 2q
= T0 δxx q + δxx
k2 (q)
r(q)
!

n, a
n, r

−
φ
q
φ
q
+
a
ρA
(p)
(p)
(p)
r − 2pn + 2pn−1 = T0 δxx pn + δxx
k2 (p)r(p)This form of the equations allows to solve for the unknowns r(q) , r(p) , using a convenient iterative root finder,such as Newton-Raphson. Existence and uniqueness of the solution may be proven by using convexity of φ ,however the proof will be omitted here. The displacements u, ζ can be recovered from r(q) , r(p) by solving thefollowing linear systems
δx− un+1 = r(q) + δx− un−1 ,δx− ζ n+1 = r(p) + δx− ζ n−1(27)

3.2 Novel non-iterative, conservative scheme
A discretisation of (9) follows as
ρAδtt u = T0 δxx u + δx+ ((µt+ ψ) g)

ρAδtt ζ = T0 δxx ζ + δx+ ((µt+ ψ) f )

(28a)

(28b)

n−1/2

n
. The system is completed by an extra
Above, it is understood that ψ = ψm−1/2 , g = gnm−1/2 and f = fm−1/2
equation, relating the time derivative of ψ to f , g, asδt+ ψ = gδt· q + f δt· p(29)

Finally, g and f are calculated as the analytic derivatives of ψ with respect to p, q, at the current time step, i.e.
gnm−1/2 = ∂q ψnfm−1/2= ∂pψ

,p=pnm−1/2 ,q=qnm−1/2(30)p=pnm−1/2 ,q=qnm−1/2

Energy conservation can be shown once more by taking taking an inner product of (28a) with δt· u, and of (28b)with δt· ζ . Summing the two resulting equations, and making use of (29) gives again an energy balance of theform (22), where in this case
2
 T 

ψ n−1/2 M
ρA 
0
n
n
n 2
n 2
n
n
kδt− u kM + kδt− ζ kM +hq , et− q iM + hp , et− p iM +
h
=
(31)
2
2
2
and where the boundary term vanishes under the same fixed end conditions as before. Inspection of the energyfunction allows to write the same stability condition as (24), but as pointed out in the remark in section 3.1,
the natural longitudinal grid spacing h(l) will be enforced here, to avoid aliasing.
Using the identityk(32)µt+ ψ n−1/2 = δt+ ψ n−1/2 + ψ n−1/22
into (28a), (28b), and owing to (29), scheme (28) can be cast in a convenient matrix-vector update of the form
n−1/2

Awn+1 = b

(33)

n+1
where wn+1 , T [un+1 ,ζζ] and where−− 41 D+x Λ f g Dx
ρA−
I − 41 D+x Λ f 2 Dx
k2

!

A=

ρA
−
I − 14 D+x Λ g 2 Dx
k2−
− 14 D+x Λ f g Dxb=ρA
− n−1 − 1 D+Λ D−ζ n−1
ψ n−1/2 ◦ gn ) − 14 D+(2un − un−1 ) + T0 Dxx un + D+x (ψx Λ g2 Dx u
4 x fg x
k2
n
n−1
n
ρA
+ (ψ
n−1/2 ◦ fn ) − 1 D+Λ D−ζ n−1 − 1 D+Λ D− un−1
ζ
ζ
ψ
(2ζ
−ζ
)
+
T
D
ζ
+
D
2
0 xx
x
x
4 x f
4 x fg x
k2

(34)
!
(35)

5511

+
T −
In the above, D−
x is the M × M − 1 space derivative matrix; Dx = − Dx ; Dxx is the second space derivative
matrix; Λ ν is the M ×M diagonal matrix where the main diagonal is the vector ν; I is the M −1×M −1 identity
matrix. Also, the symbol ◦ denotes the element-wise product. Once wn+1 is known, one may go on and updateψ using (29). One appealing feature of both schemes (21) and (28) is that they are virtually insensitive tothe choice of the potential function: in other words, the schemes leads to energy conservation regardless ofthe particular choice for φ , so long, of course, that it is positive definite, and its gradient exists and is well
behaved. One attractive possibility arises as one may use the skeleton of the non-iterative scheme to test out
various nonlinear potentials, effortlessly. One particularly useful potential, here, is given by the Morse andIngard nonlinear string model, obtained as truncated series expansion of the geometrically exact case. Thoughconservative finite difference schemes for this particular model have been derived in the past (12), the current
schemes represent a reliable tool for the analyst wishing to experiment with different models.

4

NUMERICAL EXPERIMENTS

Figure 1. Low amplitude initial conditions, U0 = 1 mm. The grey line is the solution obtained with the linear
wave equation, the black line is the iterative nonlinear scheme, and the blue line is the proposed nonlinear
non-iterative scheme. The energy error is defined as (hn−1/2 − h1/2 )/h1/2 . All simulations are obtained with asample rate fs = 48 · 104 Hz. The number of iterations for Newton-Raphson is 10.

Figure 2. Medium amplitude initial conditions, U0 = 1 cm. Colour scheme and sample rate as per Fig. 1.
In the following, the schemes will be compared. The simulations consider the lossless case, and initial conditions in the form of a gaussian distribution with zero initial velocity, i.e.
−u(0, x) = U0 e

(x/L+1/2)2
2σ 2,∂t u(0, x) = 0(36)

The free parameter U0 (measured in metres) allows to control the level of nonlinearity in the string. In thefollowing, σ = 0.24 will be used. The readout point is placed at rp = 0.72L. String parameters are selected
here as so to resemble a low pitched piano string, though only for illustrative purposes. They are given as
L = 1 m,ρ = 8000 kg/m3 ,A = 2.827 · 10−7 m2 ,T0 = 100 N,

E = 2 · 1011 Pa(37)

As a first experiment, take a look at Fig. 1 and Fig. 2. As expected, under low amplitude initial conditions theoutput of both the iterative and of the non-iterative schemes is virtually indistinguishable from the simulationof the simple linear wave equation. As the amplitude is increased, some nonlinear effects come into existence5512

and as the longitudinal component becomes more prominent. Notice that energy is conserved to machine accuracy for both schemes, although some jumps are observed in the Newton-Raphson (this could be alleviated byincreasing the number of iterations.)
A second experiment is considered in Fig. 3. There, high amplitude initial conditions are enforced, resulting instrongly nonlinear dynamics. Increasing the sample rate shows a gradual regularisation of the waveforms, which
begin to converge to a common solution, as expected.

Figure 3. High amplitude initial conditions, U0 = 10 cm. Top row: fs = 48 · 104 Hz, Middle Row: fs = 3 · 48 · 104
Hz, Bottom Row: fs = 5 · 48 · 104 Hz. Black: iterative scheme. Blue: non-iterative scheme. The number of
iterations for Newton-Raphson is 10.
Finally, Fig. 4 shows the output spectra of the simulations obtained with low, medium and high amplitude
initial conditions. Notice the progressive shift of the spectra away from the harmonic series predicted by linear
theory.

Figure 4. Output spectra of the non-iterative scheme. Left to right: U0 = 1 mm, U0 = 1 cm, U0 = 10 cm. For
all simulations fs = 48 · 104 Hz.

4.1 Computational Testing
A quick computational test was performed, on a 3.6GHz Intel core i5 CPU found in a 2011 iMac, and running
Matlab R2018b. The test was run by simply recording the compute times for 0.1s of simulation. The times
reported in table 1 are an average over 10 repetitions of the same simulation. It is seen that the new schemes
perform much more efficiently than the iterative schemes, with speedups of about half an order of magnitude.

5

CONCLUSIONSIn this work, a novel non-iterative, conservative finite difference scheme was offered for the solution of thegeometrically exact nonlinear string vibration. The new schemes are comparable to iterative schemes obtained

5513

Table 1. Compute times for 0.1 s, under high amplitude conditions. The oversampling factor is given for a basesample rate fs = 48 · 103 Hz. The iterative scheme was run with 10 iterations per time step.
Oversampling
1
5
10

Iterative (s)Non-Iterative (s)Speedup

3.1
22.7
69.3

0.65
5.2
15.5

4.80×
4.36×
4.47×

via the trapezoid rule in terms of resolving power, but they are much more efficient. Also, the skeleton of thescheme can be used to simulate an arbitrary nonlinear potential, provided that it is positive definite and that its
gradient is well defined. These schemes follow as a particular application of emerging methods in numericalanalysis, based on a quadratisation of the nonlinear potential.

ACKNOWLEDGMENTS
The first author wishes to thank the Leverhulme Trust, who is supporting his research with an Early Career
Fellowship.

References[1] A. Falaize. Modélisation, simulation, génération de code et correction de systèmes multi-physiques audios: Approche parréseau de composants et formulation hamiltonienne à ports. PhD thesis, Université Pierre et Marie Curie, Paris, July 2016.
[2] N. Lopes, T. Hélie, and A. Falaize. Explicit second-order accurate method for the passive guaranteed simulation of porthamiltonian systems. In Proc. 5th IFAC 2015, Lyon, France, July 2015.
[3] A. Falaize and T. Hélie. Passive guaranteed simulation of analog audio circuits: A port-Hamiltonian approach. Appl. Sci., 6:
273 – 273, 2016.
[4] M. Ducceschi and S. Bilbao. Non-iterative solvers for nonlinear problems: The case of collisions (under review). In 22ndInt. Conf. on Dig. Audio Eff. (DAFx 2019), Birmingham, UK, September 2019.
[5] S. Bilbao, M. Ducceschi, and C. Webb. A large-scale real-time modular physical modeling sound synthesis system (underreview). In 22nd Int. Conf. on Dig. Audio Eff. (DAFx 2019), Birmingham, UK, September 2019.
[6] M. Ducceschi and S. Bilbao. A physical model for the prepared piano. In Proc. 26th Int. Cong. of Sound Vib. (ICSV 2019),
Montreal, Canada, July 2019.
[7] X. Yang. Linear and unconditionally energy stable schemes for the binary fluid-surfactant phase field model. Comp. Methods
Appl. Mech. Eng., 318:1005–1029, 2017.
[8] J. Chabassier and S. Imperiale. Stability and dispersion analysis of improved time discretisation for prestressed Timoshenko
systems. Application to the stiff piano string. Wave Motion, 50(3):456 – 480, 4 2013.
[9] M. Ducceschi and S. Bilbao. Conservative finite difference time domain schemes for the prestressed Timoshenko, shear andEuler-Bernoulli beam equations. Wave Motion, 89:142–165, 2019.
[10] A. Chaigne and J. Kergomard. Acoustics of musical instruments. Springer, New York, USA, 2016.
[11] S. Bilbao. Numerical sound synthesis: Finite difference schemes and simulation in musical acoustics. Wiley, Chichester,UK, 2009.
[12] S. Bilbao. Conservative numerical methods for nonlinear strings. J. Acoust. Soc. Am., 118:3316–3327, 2005.

5514