Conservative Numerical Methods for Nonlinear String Dynamics: Non-planar Vibration

TL;DR

This paper introduces numerical schemes for simulating nonlinear nonplanar string vibrations that preserve energy and angular momentum, enabling stable and accurate simulations even under strong nonlinearities.

Contribution

It presents new numerical methods emphasizing conservation properties for nonlinear string dynamics, with detailed implementation and stability analysis.

Findings

Numerical schemes effectively conserve energy and angular momentum.

Stable simulations achieved under strong nonlinear conditions.

Implementation details facilitate practical application.

Abstract

Nonlinear string vibration, in particular the case of nonplanar motion, has been an area of intense study for many years. Numerical simulation methods, essential for the comparison between measured data and theory, have received somewhat less attention. In this article, various numerical schemes for nonlinear nonplanar string dynamics are presented, with an emphasis on discrete conservation of energy and angular momentum. Simple numerical stability conditions may be arrived at, even under strongly nonlinear conditions, by employing these conservation properties. Full implementation details and various numerical examples are presented, and several topics, including a discussion of numerical loss models and spectral methods, are dealt with in the Appendix.

Tables

Table 1. Table 1: Nomenclature Table.

Symbol	definition
$A$	string cross-sectional area
${𝐀}_{{𝐬}^{(\bullet )}}$	update matrix, scheme ${𝐬}^{(\bullet )}$
${𝒜}_{𝐒}$, ${𝒜}_{𝐊}$, ${𝒜}_{{𝐬}^{(\bullet )}}$, ${𝒜}_{{𝐤}^{(\bullet )}}$	angular momentum, model $𝐒$, $𝐊$, scheme ${𝐬}^{(\bullet )}$, ${𝐤}^{(\bullet )}$
${𝐁}_{{𝐬}^{(\bullet )}}$	update matrix, scheme ${𝐬}^{(\bullet )}$
${ℬ}_{𝐒,\mathscr{H}}$, ${ℬ}_{𝐊,\mathscr{H}}$	boundary term, energy conservation of system $𝐒$, $𝐊$
${ℬ}_{𝐒,𝒜}$, ${ℬ}_{𝐊,𝒜}$	boundary term, angular momentum conservation of system $𝐒$, $𝐊$
${ℬ}_{{𝐬}^{(\bullet )},\mathscr{H}}$, ${ℬ}_{{𝐤}^{(\bullet )},\mathscr{H}}$	boundary term, energy conservation of system ${𝐬}^{(\bullet )}$, ${𝐤}^{(\bullet )}$
${ℬ}_{{𝐬}^{(\bullet )},𝒜}$, ${ℬ}_{{𝐤}^{(\bullet )},𝒜}$	boundary term, angular momentum conservation of system ${𝐬}^{(\bullet )}$, ${𝐤}^{(\bullet )}$
$𝐃$	differentiation matrix
$𝒟$	discrete spatial domain
$e_{\bullet }$	shift operator
$E$	Young’s modulus
$h_t$	time step
$h_x$	grid spacing
$𝒢$, ${𝒢}_{{𝐤}^{(\bullet )}}$	nonlinear scaling, system $𝐊$, scheme ${𝐤}^{(\bullet )}$
${\mathscr{H}}_{𝐒}$, ${\mathscr{H}}_{𝐊}$, ${\mathscr{H}}_{{𝐬}^{(\bullet )}}$, ${\mathscr{H}}_{{𝐤}^{(\bullet )}}$	total energy, model $𝐒$, $𝐊$, scheme ${𝐬}^{(\bullet )}$, ${𝐤}^{(\bullet )}$
$L$	string length
$n$	time index
$p$	$={\xi }_x$ until Section 4, $={\delta }_{x-}\xi$ thereafter
$𝐏$	auxiliary scheme matrix
$𝐪$	$={𝜼}_x$ until Section 4, $={\delta }_{x-}𝜼$ thereafter
${𝐐}_{\left(1\right)},{𝐐}_{\left(2\right)}$	auxiliary scheme matrices
$t$	time variable (non-dimensionalized from Section 2.1 onwards)
$T_0$	nominal string tension
${𝒯}_{𝐒}$, ${𝒯}_{𝐊}$, ${𝒯}_{{𝐬}^{(\bullet )}}$, ${𝒯}_{{𝐤}^{(\bullet )}}$	kinetic energy, model $𝐒$, $𝐊$, scheme ${𝐬}^{(\bullet )}$, ${𝐤}^{(\bullet )}$
$𝐮$	vector of discrete longitudinal displacements
${𝐯}_{\left(1\right)},{𝐯}_{\left(2\right)}$	vectors of discrete transverse displacements
${𝒱}_{𝐒}$, ${𝒱}_{𝐊}$, ${𝒱}_{{𝐬}^{(\bullet )}}$, ${𝒱}_{{𝐤}^{(\bullet )}}$	potential energy, model $𝐒$, $𝐊$, scheme ${𝐬}^{(\bullet )}$, ${𝐤}^{(\bullet )}$
$𝐰$	discrete state update vector
$x$	spatial independent variable (non-dimensionalized from Section 2.1 onwards)
$\alpha$	dimensionless string parameter $(=T_0/EA$)
$\beta$	= $\frac{\alpha -1}{2}$
${\delta }_{\bullet }$	difference operator
${\eta }_{\left(1\right)},{\eta }_{\left(2\right)}$	transverse displacements
$𝜼$	vector transverse displacement (non-dimensionalized from Section 2.1 onwards)
$\widehat{𝜼}$	Fourier expansion coefficients of vector transverse displacement
$\lambda$	Courant number $\left(=h_t/h_x\right)$
${\mu }_{\bullet }$	discrete averaging operator
$\nu$	free parameter, implicit scheme
$\xi$	longitudinal displacement (non-dimensionalized from Section 2.1 onwards)
$\rho$	linear mass density
${\sigma }_{\xi }$, ${\sigma }_{𝜼}$	loss parameter (longitudinal, transverse)
$\tau$	free parameter, implicit scheme
$\varphi$	$=p+\frac{1-\alpha }{2}{𝐪}^T𝐪$
$𝝍$	$=\left(\alpha +\frac{1-\alpha }{2}\left({𝐪}^T𝐪+2p\right)\right)𝐪$

Table 2. Table 2: Finite difference schemes for systems S , and K . All instances of a grid function ξ 𝜉 \xi or 𝜼 𝜼 \eta   refer to that function at grid location i 𝑖 i and time step n 𝑛 n , i.e., ξ i n superscript subscript 𝜉 𝑖 𝑛 \xi_{i}^{n} and 𝜼 𝜼 \eta n i superscript subscript absent 𝑖 𝑛 {}_{i}^{n} . The quantities p 𝑝 p and 𝐪 𝐪 {\bf q} are defined in Eq. ( 35 ). All schemes are defined over the spatial interval i ∈ 𝒟 𝑖 𝒟 i\in{\mathcal{D}} , and for n ≥ 0 𝑛 0 n\geq 0 . The set 𝒟 + superscript 𝒟 {\mathcal{D}}^{+} is defined in the third paragraph of Section 4 .

	Defining Equations
${𝐬}^{\left(a\right)}$		${\varphi }_{{𝐬}^{\left(a\right)}}=p+\frac{1-\alpha }{2}{𝐪}^T𝐪$
		${𝝍}_{{𝐬}^{\left(a\right)}}=\alpha 𝐪+\frac{1-\alpha }{2}\left({𝐪}^T𝐪+2p\right)𝐪$
${𝐬}^{\left(b\right)}$		${\varphi }_{{𝐬}^{\left(b\right)}}=p+\frac{1-\alpha }{2}{𝐪}^T𝐪$
		${𝝍}_{{𝐬}^{\left(b\right)}}=\alpha 𝐪+\frac{1-\alpha }{2}\left({𝐪}^T𝐪+2p\right)\left({\mu }_{to}𝐪\right)$
${𝐬}^{\left(c\right)}$	${\delta }_{t+}{\delta }_{t-}\xi ={\delta }_{x+}{\varphi }_{\bullet }$	${\varphi }_{{𝐬}^{\left(c\right)}}=p+\frac{1-\alpha }{2}{𝐪}^T{\mu }_{to}𝐪$
	${\delta }_{t+}{\delta }_{t-}𝜼={\delta }_{x+}{𝝍}_{\bullet }$	${𝝍}_{{𝐬}^{\left(c\right)}}=\alpha 𝐪+\frac{1-\alpha }{2}{𝐪}^T𝐪{\mu }_{to}𝐪+2\left({\mu }_{t+}{\mu }_{t-}p\right)\left(𝐪\right)$
${𝐬}^{\left(d\right)}$		${\varphi }_{{𝐬}^{\left(d\right)}}=p+\frac{1-\alpha }{2}{𝐪}^T{\mu }_{to}𝐪$
		${𝝍}_{{𝐬}^{\left(d\right)}}=\alpha 𝐪+\frac{1-\alpha }{2}{𝐪}^T\left({\mu }_{to}𝐪\right)𝐪+2\left({\mu }_{t+}{\mu }_{t-}p\right)\left(𝐪\right)$
${𝐬}^{\left(e\right)}$		${\varphi }_{{𝐬}^{\left(e\right)}}=p+\frac{1-\alpha }{2}{\mu }_{to}\left({𝐪}^T𝐪\right)$
		${𝝍}_{{𝐬}^{\left(e\right)}}=\alpha 𝐪+\frac{1-\alpha }{2}\left({\mu }_{to}\left({𝐪}^T𝐪+2p\right)\right)\left({\mu }_{to}𝐪\right)$
${𝐤}^{\left(a\right)}$	${\delta }_{t+}{\delta }_{t-}𝜼=\alpha {𝒢}_{\bullet }{\delta }_{x+}𝐪$	${𝒢}_{{𝐤}^{\left(a\right)}}=1+\frac{1}{2\alpha }{\Vert 𝐪\Vert }_{D^{+}}^2$
${𝐤}^{\left(b\right)}$		${𝒢}_{{𝐤}^{\left(b\right)}}=1+\frac{1}{2\alpha }{\mu }_{t+}{⟨𝐪,e_{t-}𝐪⟩}_{D^{+}}$

Table 3. Table 3: Conserved angular momentum, and boundary terms for the schemes given in Table 2 . In all cases for which expressions are provided, it is true that δ t + ​ 𝒜 ∙ = ℬ ∙ subscript 𝛿 limit-from 𝑡 subscript 𝒜 ∙ subscript ℬ ∙ \delta_{t+}{\mathcal{A}}_{\bullet}={\mathcal{B}}_{\bullet} . The lack of a conserved angular momentum is indicated by the symbol “ − ⁣ − -- ”.

	Conserved angular momentum	Boundary term
${𝐬}^{\left(a\right)}$	${𝒜}_{{𝐬}^{\left(a\right)}}={⟨{\mu }_{t-}\widetilde{𝜼},{\delta }_{t-}𝜼⟩}_{𝒟}$	${ℬ}_{𝒜,{𝐬}^{\left(a\right)}}={\widetilde{𝜼}}_N^T{𝝍}_{N+1}-{\widetilde{𝜼}}_0^T{𝝍}_0$
${𝐬}^{\left(b\right)}$	${𝒜}_{{𝐬}^{\left(b\right)}}={⟨{\mu }_{t-}\widetilde{𝜼},{\delta }_{t-}𝜼⟩}_{𝒟}+\frac{\alpha h_t}{2}{⟨\widetilde{𝐪},e_{t-}𝐪⟩}_{{𝒟}^{+}}$	${ℬ}_{𝒜,{𝐬}^{\left(b\right)}}=\left({\mu }_{to}{\widetilde{𝜼}}_N^T\right){𝝍}_{N+1}-\left({\mu }_{to}{\widetilde{𝜼}}_0^T\right){𝝍}_0$
${𝐬}^{\left(c\right)}$	$--$	$--$
${𝐬}^{\left(d\right)}$	${𝒜}_{{𝐬}^{\left(d\right)}}={⟨{\mu }_{t-}\widetilde{𝜼},{\delta }_{t-}𝜼⟩}_{𝒟}$	${ℬ}_{𝒜,{𝐬}^{\left(d\right)}}={\widetilde{𝜼}}_N^T{𝝍}_{N+1}-{\widetilde{𝜼}}_0^T{𝝍}_0$
${𝐬}^{\left(e\right)}$	${𝒜}_{{𝐬}^{\left(e\right)}}={⟨{\mu }_{t-}\widetilde{𝜼},{\delta }_{t-}𝜼⟩}_{𝒟}+\frac{\alpha h_t}{2}{⟨\widetilde{𝐪},e_{t-}𝐪⟩}_{{𝒟}^{+}}$	${ℬ}_{𝒜,{𝐬}^{\left(e\right)}}=\left({\mu }_{to}{\widetilde{𝜼}}_N^T\right){𝝍}_{N+1}-\left({\mu }_{to}{\widetilde{𝜼}}_0^T\right){𝝍}_0$
${𝐤}^{\left(a\right)}$	${𝒜}_{{𝐤}^{\left(a\right)}}={⟨{\mu }_{t-}\widetilde{𝜼},{\delta }_{t-}𝜼⟩}_{𝒟}$	${ℬ}_{𝒜,{𝐤}^{\left(a\right)}}={𝒢}_{{𝐤}^{\left(a\right)}}\left({\widetilde{𝜼}}_N^T{𝐪}_{N+1}-{\widetilde{𝜼}}_0^T{𝐪}_0\right)$
${𝐤}^{\left(b\right)}$	${𝒜}_{{𝐤}^{\left(b\right)}}={⟨{\mu }_{t-}\widetilde{𝜼},{\delta }_{t-}𝜼⟩}_{𝒟}$	${ℬ}_{𝒜,{𝐤}^{\left(b\right)}}={𝒢}_{{𝐤}^{\left(b\right)}}\left({\widetilde{𝜼}}_N^T{𝐪}_{N+1}-{\widetilde{𝜼}}_0^T{𝐪}_0\right)$

Table 4. Table 4: Discrete kinetic and potential energies for the schemes given in Table 2 ; their sum will be conserved. The symbol − ⁣ − -- indicates that the scheme is not conservative. When such quantities are indicated, the discrete energy conservation property δ t + ​ ( 𝒯 𝐬 ( ∙ ) + 𝒱 𝐬 ( ∙ ) ) = ℬ ℋ , 𝐬 ( ∙ ) subscript 𝛿 limit-from 𝑡 subscript 𝒯 superscript 𝐬 ∙ subscript 𝒱 superscript 𝐬 ∙ subscript ℬ ℋ superscript 𝐬 ∙ \delta_{t+}\left({\mathcal{T}}_{{\bf s}^{(\bullet)}}+{\mathcal{V}}_{{\bf s}^{(\bullet)}}\right)={\mathcal{B}}_{{\mathcal{H}},{\bf s}^{(\bullet)}} holds, where ℬ ℋ , 𝐬 ( ∙ ) subscript ℬ ℋ superscript 𝐬 ∙ {\mathcal{B}}_{{\mathcal{H}},{\bf s}^{(\bullet)}} is given in Eq. ( 37 ).

	Kinetic energy	Potential energy
${𝐬}^{\left(a\right)}$	$--$	$--$
${𝐬}^{\left(b\right)}$		${𝒱}_{{𝐬}^{\left(b\right)}}=\frac{\alpha }{2}{⟨p,e_{t-}p⟩}_{{𝒟}^{+}}+\frac{\alpha }{2}{⟨𝐪,e_{t-}𝐪⟩}_{{𝒟}^{+}}$
		$+\frac{1-\alpha }{2}{⟨p+\frac{1}{2}{𝐪}^T𝐪,e_{t-}\left(p+\frac{1}{2}{𝐪}^T𝐪\right)⟩}_{{𝒟}^{+}}$
${𝐬}^{\left(c\right)}$		${𝒱}_{{𝐬}^{\left(c\right)}}=\frac{1}{2}{⟨p,e_{t-}p⟩}_{{𝒟}^{+}}+\frac{\alpha }{2}{⟨𝐪,e_{t-}𝐪⟩}_{{𝒟}^{+}}$
		$+\frac{1-\alpha }{2}\left({\Vert {\mu }_{t-}p+\frac{1}{2}{𝐪}^T{e}_{t-}𝐪\Vert }_{{𝒟}^{+}}^2-{\Vert {\mu }_{t-}p\Vert }_{{𝒟}^{+}}^2\right)$
	${𝒯}_{{𝐬}^{(\bullet )}}=\frac{1}{2}\left({\Vert {\delta }_{t-}\xi \Vert }_{𝒟}^2+{\Vert {\delta }_{t-}𝜼\Vert }_{𝒟}^2\right)$	$+\frac{1-\alpha }{8}{\Vert {\widetilde{𝐪}}^T{e}_{t-}𝐪\Vert }_{{𝒟}^{+}}^2$
${𝐬}^{\left(d\right)}$		${𝒱}_{{𝐬}^{\left(d\right)}}=\frac{1}{2}{⟨p,e_{t-}p⟩}_{{𝒟}^{+}}+\frac{\alpha }{2}{⟨𝐪,e_{t-}𝐪⟩}_{{𝒟}^{+}}$
		$+\frac{1-\alpha }{2}\left({\Vert {\mu }_{t-}p+\frac{1}{2}{𝐪}^T{e}_{t-}𝐪\Vert }_{{𝒟}^{+}}^2-{\Vert {\mu }_{t-}p\Vert }_{{𝒟}^{+}}^2\right)$
${𝐬}^{\left(e\right)}$		${𝒱}_{{𝐬}^{\left(e\right)}}=\frac{1}{2}{⟨p,e_{t-}p⟩}_{{𝒟}^{+}}+\frac{\alpha }{2}{⟨𝐪,e_{t-}𝐪⟩}_{{𝒟}^{+}}$
		$+\frac{1-\alpha }{2}{\mu }_{t-}\left({\Vert p+\frac{1}{2}{𝐪}^T𝐪\Vert }_{{𝒟}^{+}}^2-{\Vert p\Vert }_{{𝒟}^{+}}^2\right)$
${𝐤}^{\left(a\right)}$	$--$	$--$
${𝐤}^{\left(b\right)}$	${𝒯}_{{𝐤}^{\left(b\right)}}=\frac{1}{2}{\Vert {\delta }_{t-}𝜼\Vert }_{𝒟}^2$	${𝒱}_{{𝐤}^{\left(b\right)}}=\frac{\alpha }{2}{⟨𝐪,e_{t-}𝐪⟩}_{{𝒟}^{+}}\left(1+\frac{1}{4\alpha }{⟨𝐪,e_{t-}𝐪⟩}_{{𝒟}^{+}}\right)$

Table 5. Table 5: Stability conditions and bounds on solution size for the schemes given in Table 2 . Schemes not possessing such conditions are indicated by the symbol − ⁣ − -- in the accompanying rows. Two types of bounds are given: in the third column, general bounds on the size of the solution, which hold under any conservative boundary conditions, and in the fourth column, better bounds available, when one of the ends of the string is fixed.

	Stab. conditions	General bounds	Bounds under fixed conditions
${𝐬}^{\left(a\right)}$	$--$	$--$	$--$
${𝐬}^{\left(b\right)}$	$--$	$--$	$--$
${𝐬}^{\left(c\right)}$	$\lambda \le \sqrt{1/\alpha }$	${\Vert {\xi }^n\Vert }_{𝒟}\le {\Vert {\xi }^0\Vert }_{𝒟}+h_{t}n\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(c\right)}}^0}{1-{\lambda }^2}}$	${\Vert \xi \Vert }_{𝒟}\le N{h}_x\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(c\right)}}^0}{\alpha }}+\frac{h_t}{2}\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(c\right)}}^0}{1-{\lambda }^2}}$
	$\lambda \le 1$	${\Vert {𝜼}^n\Vert }_{𝒟}\le {\Vert {𝜼}^0\Vert }_{𝒟}+h_{t}n\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(c\right)}}^0}{1-\alpha {\lambda }^2}}$	${\Vert 𝜼\Vert }_{𝒟}\le N{h}_x\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(c\right)}}^0}{\alpha }}+\frac{h_t}{2}\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(c\right)}}^0}{1-\alpha {\lambda }^2}}$
${𝐬}^{\left(d\right)}$	$\lambda \le \sqrt{1/\alpha }$	${\Vert {\xi }^n\Vert }_{𝒟}\le {\Vert {\xi }^0\Vert }_{𝒟}+h_{t}n\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(d\right)}}^0}{1-{\lambda }^2}}$	${\Vert \xi \Vert }_{𝒟}\le N{h}_x\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(c\right)}}^0}{\alpha }}+\frac{h_t}{2}\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(d\right)}}^0}{1-{\lambda }^2}}$
	$\lambda \le 1$	${\Vert {𝜼}^n\Vert }_{𝒟}\le {\Vert {𝜼}^0\Vert }_{𝒟}+h_{t}n\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(d\right)}}^0}{1-\alpha {\lambda }^2}}$	${\Vert 𝜼\Vert }_{𝒟}\le N{h}_x\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(d\right)}}^0}{\alpha }}+\frac{h_t}{2}\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(d\right)}}^0}{1-\alpha {\lambda }^2}}$
${𝐬}^{\left(e\right)}$	$\lambda \le \sqrt{1/\alpha }$	${\Vert {\xi }^n\Vert }_{𝒟}\le {\Vert {\xi }^0\Vert }_{𝒟}+h_{t}n\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(e\right)}}^0}{1-\left(2-\alpha \right){\lambda }^2}}$	${\Vert \xi \Vert }_{𝒟}\le N{h}_x\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(e\right)}}^0}{\alpha }}+\frac{h_t}{2}\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(e\right)}}^0}{1-\left(2-\alpha \right){\lambda }^2}}$
	$\lambda \le \sqrt{\frac{1}{2-\alpha }}$	${\Vert {𝜼}^n\Vert }_{𝒟}\le {\Vert {𝜼}^0\Vert }_{𝒟}+h_{t}n\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(e\right)}}^0}{1-\alpha {\lambda }^2}}$	${\Vert 𝜼\Vert }_{𝒟}\le N{h}_x\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(e\right)}}^0}{\alpha }}+\frac{h_t}{2}\sqrt{\frac{2{\mathscr{H}}_{{𝐬}^{\left(e\right)}}^0}{1-\alpha {\lambda }^2}}$
${𝐤}^{\left(a\right)}$	$--$	$--$	$--$
${𝐤}^{\left(b\right)}$	$\lambda \le \sqrt{1/\alpha }$	${\Vert {𝜼}^n\Vert }_{𝒟}\le {\Vert {𝜼}^0\Vert }_{𝒟}+h_{t}n\sqrt{\frac{2{\mathscr{H}}_{{𝐤}^{\left(b\right)}}^0}{1-\alpha {\lambda }^2}}$	${\Vert 𝜼\Vert }_{𝒟}\le N{h}_x\sqrt{\frac{2{\mathscr{H}}_{{𝐤}^{\left(b\right)}}^0}{\alpha }}+\frac{h_t}{2}\sqrt{\frac{2{\mathscr{H}}_{{𝐤}^{\left(b\right)}}^0}{1-\alpha {\lambda }^2}}$

Table 6. Table 6: State update matrices for schemes 𝐬 ( ∙ ) superscript 𝐬 ∙ {\bf s}^{(\bullet)} , under fixed boundary conditions. The block matrices 𝐀 𝐬 ( ∙ ) subscript 𝐀 superscript 𝐬 ∙ {\bf A}_{{\bf s}^{(\bullet)}} and 𝐁 𝐬 ( ∙ ) subscript 𝐁 superscript 𝐬 ∙ {\bf B}_{{\bf s}^{(\bullet)}} , given below, are used in the update of Eq. ( 56 ). The symbol ⋅ ⋅ \cdot indicates a an N − 1 𝑁 1 N-1 by N − 1 𝑁 1 N-1 zero matrix. System 𝐬 ( e ) superscript 𝐬 𝑒 {\bf s}^{(e)} does not have a state update form.

${𝐬}^{\left(a\right)}$
${𝐬}^{\left(b\right)}$
${𝐬}^{\left(c\right)}$
${𝐬}^{\left(d\right)}$
${𝐬}^{\left(e\right)}$	No matrix update form

Table 7. Table 7: Conserved angular momentum (nondimensional, × 10 − 7 absent superscript 10 7 \times 10^{-7} ) for various algorithms for the string described in the caption of Fig. 2 , against time step n 𝑛 n .

n	${𝒜}_{{𝐬}^{\left(a\right)}}$	${𝒜}_{{𝐬}^{\left(b\right)}}$	${𝒜}_{{𝐬}^{\left(d\right)}}$	${𝒜}_{{𝐤}^{\left(a\right)}}$	${𝒜}_{{𝐤}^{\left(b\right)}}$
1	2.0000000000000	1.9999950753362	2.0000000000000	2.0000000000000	2.0000000000000
2	2.0000000000000	1.9999950753362	2.0000000000000	2.0000000000000	2.0000000000000
3	2.0000000000000	1.9999950753362	2.0000000000000	2.0000000000000	2.0000000000000
4	2.0000000000000	1.9999950753362	2.0000000000000	2.0000000000000	2.0000000000000
5	2.0000000000000	1.9999950753362	2.0000000000000	2.0000000000000	2.0000000000000
100	2.0000000000000	1.9999950753362	2.0000000000007	2.0000000000000	2.0000000000001

Table 8. Table 8: Conserved energy (nondimensional, × 10 − 7 absent superscript 10 7 \times 10^{-7} ) for various algorithms for the string described in the caption of Fig. 2 , against time step n 𝑛 n .

n	${\mathscr{H}}_{{𝐬}^{\left(b\right)}}$	${\mathscr{H}}_{{𝐬}^{\left(c\right)}}$	${\mathscr{H}}_{{𝐬}^{\left(d\right)}}$	${\mathscr{H}}_{{𝐤}^{\left(b\right)}}$
1	9.245104334637	9.245104334637	9.245104316451	6.821328138420
2	9.245104334637	9.245104334637	9.245104316451	6.821328138420
3	9.245104334637	9.245104334637	9.245104316451	6.821328138420
4	9.245104334637	9.245104334637	9.245104316451	6.821328138420
5	9.245104334637	9.245104334637	9.245104316451	6.821328138420
100	9.245104334635	9.245104334635	9.245104316452	6.821328138419

Table 9. Table 9: Conserved angular momentum 𝒜 𝐬 ( d ) subscript 𝒜 superscript 𝐬 𝑑 {\mathcal{A}}_{{\bf s}^{(d)}} (nondimensional) for scheme 𝐬 ( d ) superscript 𝐬 𝑑 {\bf s}^{(d)} , against time step n 𝑛 n , for the string of parameters defined in the caption to Fig. 2 , under initial conditions ( 58 ) of increasing magnitude. Values of γ 1 subscript 𝛾 1 \gamma_{1} and γ 2 subscript 𝛾 2 \gamma_{2} are given in the table.

n	${\gamma }_1=0.001{\gamma }_2=0.00001$	${\gamma }_1=0.01{\gamma }_2=0.0001$	${\gamma }_1=0.02{\gamma }_2=0.0002$	${\gamma }_1=0.04{\gamma }_2=0.0004$
1	5.000000000000$\times {10}^{-9}$	5.000000000000$\times {10}^{-7}$	2.000000000000$\times {10}^{-6}$	8.000000000000$\times {10}^{-6}$
2	5.000000000000$\times {10}^{-9}$	5.000000000000$\times {10}^{-7}$	2.000000000000$\times {10}^{-6}$	8.000000000000$\times {10}^{-6}$
3	5.000000000000$\times {10}^{-9}$	5.000000000000$\times {10}^{-7}$	2.000000000000$\times {10}^{-6}$	8.000000000000$\times {10}^{-6}$
4	5.000000000000$\times {10}^{-9}$	5.000000000000$\times {10}^{-7}$	2.000000000000$\times {10}^{-6}$	8.000000000000$\times {10}^{-6}$
5	5.000000000000$\times {10}^{-9}$	5.000000000000$\times {10}^{-7}$	2.000000000000$\times {10}^{-6}$	8.000000000000$\times {10}^{-6}$
100	5.000000000002$\times {10}^{-9}$	5.000000000002$\times {10}^{-7}$	1.999999999999$\times {10}^{-6}$	8.000000000000$\times {10}^{-6}$

Table 10. Table 10: Conserved energy ℋ 𝐬 ( d ) subscript ℋ superscript 𝐬 𝑑 {\mathcal{H}}_{{\bf s}^{(d)}} (nondimensional) for scheme 𝐬 ( d ) superscript 𝐬 𝑑 {\bf s}^{(d)} , against time step n 𝑛 n , for the string of parameters defined in the caption to Fig. 2 , under initial conditions ( 58 ) of increasing magnitude. Values of γ 1 subscript 𝛾 1 \gamma_{1} and γ 2 subscript 𝛾 2 \gamma_{2} are given in the table.

n	${\gamma }_1=0.001{\gamma }_2=0.00001$	${\gamma }_1=0.01{\gamma }_2=0.0001$	${\gamma }_1=0.02{\gamma }_2=0.0002$	${\gamma }_1=0.04{\gamma }_2=0.0004$
1	5.22012775452$\times {10}^{-10}$	9.72106301924$\times {10}^{-8}$	9.34410431645$\times {10}^{-7}$	1.24667283005$\times {10}^{-5}$
2	5.22012775452$\times {10}^{-10}$	9.72106301924$\times {10}^{-8}$	9.34410431645$\times {10}^{-7}$	1.24667283005$\times {10}^{-5}$
3	5.22012775452$\times {10}^{-10}$	9.72106301924$\times {10}^{-8}$	9.34410431645$\times {10}^{-7}$	1.24667283005$\times {10}^{-5}$
4	5.22012775452$\times {10}^{-10}$	9.72106301924$\times {10}^{-8}$	9.34410431645$\times {10}^{-7}$	1.24667283005$\times {10}^{-5}$
5	5.22012775452$\times {10}^{-10}$	9.72106301924$\times {10}^{-8}$	9.34410431645$\times {10}^{-7}$	1.24667283005$\times {10}^{-5}$
100	5.22012775452$\times {10}^{-10}$	9.72106301925$\times {10}^{-8}$	9.34410431645$\times {10}^{-7}$	1.24667283005$\times {10}^{-5}$

Table 11. Table 11: Summary of properties of schemes 𝐬 ( ∙ ) superscript 𝐬 ∙ {\bf s}^{(\bullet)} and 𝐤 ( ∙ ) superscript 𝐤 ∙ {\bf k}^{(\bullet)} .

	angular momentum conserving	energy conserving	stability conditions	existence/uniqueness
${𝐬}^{\left(a\right)}$	Yes	No	No	Yes
${𝐬}^{\left(b\right)}$	Yes	Yes	No	Yes
${𝐬}^{\left(c\right)}$	No	Yes	Yes	Yes
${𝐬}^{\left(d\right)}$	Yes	Yes	Yes	Yes
${𝐬}^{\left(e\right)}$	Yes	Yes	Yes	No
${𝐤}^{\left(a\right)}$	Yes	No	No	Yes
${𝐤}^{\left(b\right)}$	Yes	Yes	Yes	Yes

Equations

$$\rho {\xi }_{tt}$$

$$\displaystyle\rho\mbox{\boldmath$\eta$}_{tt}\!$$

$$x^{\prime}=x/L\qquad\xi=\xi/L\qquad\mbox{\boldmath$\eta$}=\mbox{\boldmath$\eta$}/L\qquad t^{\prime}=t\sqrt{\frac{EA}{\rho L^{2}}}$$

$${\xi }_{tt}$$

$$\displaystyle\mbox{\boldmath$\eta$}_{tt}$$

$${\xi }_{tt}$$

$$\displaystyle\mbox{\boldmath$\eta$}_{tt}\!$$

$$p=\xi_{x}\qquad{\bf q}=\mbox{\boldmath$\eta$}_{x}$$

$$\mbox{\boldmath$\eta$}_{tt}=\alpha{\mathcal{G}}{\bf q}_{x}$$

$$\mathcal{G}=1+\frac{1}{2\alpha }{\int }_0^1{\mathbf{q}}^T\mathbf{q}dx$$

$$\int_{0}^{1}\xi_{t}\xi_{tt}-\xi_{t}\phi_{x}dx=0\qquad\qquad\int_{0}^{1}\mbox{\boldmath$\eta$}_{t}^{T}\mbox{\boldmath$\eta$}_{tt}-\mbox{\boldmath$\eta$}_{t}^{T}\mbox{\boldmath$\psi$}_{x}dx=0$$

$$\displaystyle\int_{0}^{1}\xi_{t}\xi_{tt}+\mbox{\boldmath$\eta$}_{t}^{T}\mbox{\boldmath$\eta$}_{tt}+p_{t}\phi+{\bf q}^{T}_{t}\mbox{\boldmath$\psi$}dx={\mathcal{B}}_{{\bf S},{\mathcal{H}}}$$

$${\mathcal{B}}_{{\bf S},{\mathcal{H}}}=\left(\xi_{t}\phi+\mbox{\boldmath$\eta$}^{T}_{t}\mbox{\boldmath$\psi$}\right)\Big{|}_{0}^{1}$$

$$\frac{d}{dt}{\mathcal{H}}_{\mathbf{S}}={\mathcal{B}}_{\mathbf{S},\mathcal{H}}\text{mbox}\mathrm{for}{\mathcal{H}}_{\mathbf{S}}={\mathcal{T}}_{\mathbf{S}}+{\mathcal{V}}_{\mathbf{S}}$$

$${\mathcal{T}}_{\mathbf{S}}$$

$${\mathcal{V}}_{\mathbf{S}}$$

$$\frac{d}{dt}{\mathcal{H}}_{\mathbf{K}}={\mathcal{B}}_{\mathbf{K},\mathcal{H}}\text{mbox}\mathrm{for}{\mathcal{H}}_{\mathbf{K}}={\mathcal{T}}_{\mathbf{K}}+{\mathcal{V}}_{\mathbf{K}}$$

$${\mathcal{T}}_{\mathbf{K}}$$

$${\mathcal{V}}_{\mathbf{K}}$$

$${\mathcal{B}}_{\mathbf{K},\mathcal{H}}$$

$${\tilde{\mathbf{f}}}^T\mathbf{f}=0$$

$$\int_{0}^{1}\tilde{\mbox{\boldmath$\eta$}}^{T}\mbox{\boldmath$\eta$}_{tt}-\tilde{\mbox{\boldmath$\eta$}}^{T}\mbox{\boldmath$\psi$}_{x}dx=0$$

$$\int_{0}^{1}\left((\tilde{\mbox{\boldmath$\eta$}}^{T}\mbox{\boldmath$\eta$}_{t})_{t}-\tilde{\mbox{\boldmath$\eta$}}^{T}_{t}\mbox{\boldmath$\eta$}_{t}\right)+\tilde{{\bf q}}^{T}\mbox{\boldmath$\psi$}dx={\mathcal{B}}_{{\bf S},{\mathcal{A}}}$$

$${\mathcal{B}}_{{\bf S},{\mathcal{A}}}=\tilde{\mbox{\boldmath$\eta$}}^{T}\mbox{\boldmath$\psi$}\Big{|}_{0}^{1}$$

$$\frac{d}{dt}{\mathcal{A}}_{\mathbf{S}}={\mathcal{B}}_{\mathbf{S},\mathcal{A}}$$

$${\mathcal{A}}_{{\bf S}}=\int_{0}^{1}\tilde{\mbox{\boldmath$\eta$}}^{T}\mbox{\boldmath$\eta$}_{t}dx$$

$$\frac{d}{dt}{\mathcal{A}}_{\mathbf{K}}={\mathcal{B}}_{\mathbf{K},\mathcal{A}}$$

$${\mathcal{A}}_{\mathbf{K}}$$

$${\mathcal{B}}_{\mathbf{K},\mathcal{A}}$$

$$\xi=0\qquad\mbox{{\rm and}}\qquad\mbox{\boldmath$\eta$}={\bf 0}\\$$

$$p=0\text{mbox}\mathrm{and}\mathbf{q}=0$$

$${\mathcal{H}}_{\mathbf{S}}(t)={\mathcal{H}}_{\mathbf{S}}(0){\mathcal{A}}_{\mathbf{S}}(t)={\mathcal{A}}_{\mathbf{S}}(0)$$

$${\mathcal{H}}_{\mathbf{K}}(t)={\mathcal{H}}_{\mathbf{K}}(0){\mathcal{A}}_{\mathbf{K}}(t)={\mathcal{A}}_{\mathbf{K}}(0)$$

$$⟨\mathbf{f},\mathbf{g}⟩={\int }_0^1{\mathbf{f}}^T\mathbf{g}dx\Vert \mathbf{f}\Vert =⟨\mathbf{f},\mathbf{f}{⟩}^{1/2}$$

Full text

Conservative Numerical Methods for Nonlinear String Dynamics: Non-planar Vibration

Technical Report, Acoustics and Audio Group, University of Edinburgh

Stefan Bilbao

May 9 2006

Nonlinear string vibration, in particular the case of nonplanar motion, has been an area of intense study for many years. Numerical simulation methods, essential for the comparison between measured data and theory, have received somewhat less attention. In this article, various numerical schemes for nonlinear nonplanar string dynamics are presented, with an emphasis on discrete conservation of energy and angular momentum. Simple numerical stability conditions may be arrived at, even under strongly nonlinear conditions, by employing these conservation properties. Full implementation details and various numerical examples are presented, and several topics, including a discussion of numerical loss models and spectral methods, are dealt with in the Appendix.

1 Introduction

The dynamics of strings under nonlinear conditions has been under study at least since the early work of Kirchhoff [1] and then Carrier [2]; later, Anand [3] and Narasimha [4] extended this work considerably. This easily defined, but nonetheless extremely complex system exhibits a wide variety of behaviour characteristic of nonlinear systems; of particular interest is the “whirling” phenomenon, which is peculiar to the case of non-planar motion [5], [6]. Nonlinear string dynamics in three dimensions serves as an excellent test problem not only for a variety of analysis techniques, but also for the construction of numerical methods, which are the focus of this article.

The study of numerical methods which inherit discrete conservation laws from a continuous model system has been ongoing for some time [7], [8], and relates to early work on the so-called energy-method [9]. Most often, the systems under study are of a general form, e.g., the nonlinear Klein-Gordon equation, which was approached by Vu-quoc and Li [10] and Li and Vu-quoc [11], and various systems including a single polarization transverse-only model of string dynamics which were discussed by Furihata [12]. In all these cases, the form of the nonlinearity is left unspecified. In various commonly-encountered models of string dynamics, the nonlinearity is often simplified using a series approximation. Such a simplification allows obvious benefits in the analysis of the string, particularly when the nonlinearity is approximated using quadratic or cubic terms [13], as it often is. At the same time, the possibility of exploiting various algebraic symmetries in a nonlinear difference scheme also appears; as will be shown in this article, there are many distinct ways of designing numerical schemes for the same nonlinear string system, which vary considerably in terms of their conservation properties, ease of use, and, most importantly, their stability properties. Conservative difference schemes for the Kirchhoff-Carrier string, and for a more general string undergoing planar motion have been discussed by this author in [14] and [15], respectively.

It is perhaps worth mentioning here that one extremely interesting recent application of numerical simulation techniques for nonlinear mechanical systems such as the string is in the area of musical sound synthesis. Such physical modeling synthesis, as it is often called, has been in existence for some time now, and for some systems, real-time performance is now possible using personal computers. There has been some important work on using standard numerical methods such as finite difference schemes [16], [17], [18], but the dominant techniques have been based around efficient structures with their roots in digital filter design; probably the best known are digital waveguides [19]. Recently, there have been some attempts at sound synthesis based on nonlinear string vibration, in order to model effects such as pitch glides under high amplitude plucking conditions [20], as well as the phantom partial phenomenon which occurs in piano strings [21]. Though this article is intended for a general audience, some commentary on this topic will appear at various points throughout this article.

In Section 2, a general model of nonlinear string dynamics is presented, followed by two simpler forms, one employing the series approximations to the nonlinearity mentioned above (system S), and a further simplified form of the Kirchhoff-Carrier variety (system K). The section concludes with a brief presentation of the energy and angular momentum conservation properties of these models, and in particular, the bounds on the growth of the solution which result from the former property. Section 3 is a short recap of the properties of finite difference operators and inner product spaces, with a view toward applications in the construction of conservative schemes. In Section 4, various difference approximations to systems S and K are presented, followed by an analysis of their discrete energy and angular momentum conservation properties. The discrete conserved energy for each of the schemes is further examined, first to determine conditions under which, when it indeed exists, it is positive, and if so, what bounds may be placed on the size of the solution. The section concludes with a brief look at the schemes in the forms in which they will be implemented. Numerical results are presented in Section 5, with a special emphasis on phenomena which are inherent to motion in three dimensions, and in particular the so-called “whilrling” behaviour [5], [6]. Several other topics are briefly addressed in the Appendices, namely the generalization to the case of linear damping in Appendix A, a loosening of the stability condition on the time step in Appendix B, and finally, in Appendix C, a brief look at a spectral-type method for the integration of system K.

2 Nonlinear String Models

A general model of nonlinear string dynamics, discussed by many authors, and summarized succinctly by Morse and Ingard [22] and which can be related to the geometrically-exact theory of beams [23], is given by the following set of equations:

[TABLE]

Here, $\xi(x,t)$ and the two-element column vector \mbox{\boldmath\eta}(x,t)=[\eta_{(1)},\eta_{(2)}]^{T} describe, respectively, the longitudinal and transverse deviation of a point on the string as a function of time $t\geq 0$ and distance along the string $x\in[0,L]$. (The superscript T indicates a vector transpose.) Such a point, located at Cartesian coordinates $(x,0,0)$ when the string is at rest, will have dynamic coordinates $(x+\xi,\eta_{(1)},\eta_{(2)})$. See Figure 1. $E$, $A$, $\rho$ and $T_{0}$ are Young’s modulus, cross-sectional area, linear mass density, and nominal tension for the string, all assumed constant here. Subscripts indicate differentiation with respect to the named independent variable. System (1) is by no means the most general model of string dynamics; higher order effects may be modelled as well, as per the work of Narasimha [4] and Kurmyshev [24]; such improved models may fall outside the range of the techniques presented here.

System (1) requires the specification of the initial conditions $\xi(x,0)$, $\xi_{t}(x,0)$, \mbox{\boldmath\eta}(x,0), and \mbox{\boldmath\eta}_{t}(x,0), in order that the solution exist and be uniquely defined. A discussion of boundary conditions is postponed briefly until Section 2.4.

2.1 A Nondimensionalized Form

System (1) may be nondimensionalized by introducing the variables

[TABLE]

which, when substituted in (1) and primes removed, leads to

[TABLE]

which depends on a single parameter $\alpha=T_{0}/EA$, and which is defined over $x\in[0,1]$

2.2 Approximate Forms

There are various levels of approximation to system (3); the most direct forms result from Taylor series approximations to the nonlinearity (i.e., the term in large parentheses). An approximation to first order uncouples the longitudinal and transverse motion (and the two transverse polarizations from one another), yielding linear wave equations, with wave speeds 1 (longitudinal) and $\sqrt{\alpha}$ (both transverse polarizations). An approximation to second order is sometimes employed [21], [25], but most common in the study of nonlinear string vibration is an approximation to third order, following Anand [3] and Morse [22] in making use of the observation that \xi=O(|\mbox{\boldmath\eta}|^{2}), which is true for metallic strings:

System S

[TABLE]

For notational simplicity, the symbols $p$ and ${\bf q}$ have been introduced; they are defined as

[TABLE]

Numerical methods for system ${\bf S}$, in its planar form, have been applied to the problem of piano string vibration at high amplitudes [21], in order to generate perceptually important “phantom partials,” [26] which result from coupling between longitudinal and transverse motion.

Under certain conditions, namely that (1) $1\gg\alpha$, (2) the term $\xi_{tt}$ in Eq. (4a) may be neglected, and (3) the longitudinal displacement $\xi(x,t)$ is zero at $x=0$ and $x=1$, system S may be reduced to a single equation in the vector transverse displacement $\eta$  alone:

System K

[TABLE]

where ${\mathcal{G}}$ is defined by

[TABLE]

System K above, often referred to as being of Kirchhoff-Carrier type [1], [2], is far simpler to deal with that system S, both analytically and numerically, for the simple reason that the nonlinearity, characterized by ${\mathcal{G}}$, is averaged over the string length, and does not have any spatial dependence (i.e., it is a scalar quantity). Eq. (6) can be viewed, in a crude sense, as a wave equation with a wave speed which is dependent on variations in total string length. This type of system has been employed, in the context of digital waveguides [19] in order to produce synthetic sound for strings under high-amplitude plucking conditions [20]. This type of nonlinearity is often referred to as “tension modulation” in the sound synthesis community.

2.3 Conserved Quantities

Given systems S and K, which will serve as the models to be solved numerically in this article, it is worthwhile to spend some time examining their conservation properties. In the present case of continuously variable systems, this is quite straightforward.

2.3.1 Energy Conservation

For system S, multiplying Eq. (4a) by $\xi_{t}$ and left-multiplying Eq. (4b) by \eta$${}^{T}_{t} and then integrating over the interval $[0,1]$ gives

[TABLE]

Integrating by parts for the latter term under the integral in each equation and adding the results gives, employing definitions Eqs. (5),

[TABLE]

where the boundary term ${\mathcal{B}}_{{\bf S},{\mathcal{H}}}$ is given by

[TABLE]

Eq. (2.3.1) can be written as

[TABLE]

where $d/dt$ represents a total derivative with respect to time $t$, and the scalar quantities ${\mathcal{T}}_{{\bf S}}$ and ${\mathcal{V}}_{{\bf S}}$ are given by

[TABLE]

Clearly, ${\mathcal{T}}_{{\bf S}}$ and ${\mathcal{V}}_{{\bf S}}$ represent the kinetic and potential energy, respectively, of the string as described by system S, and ${\mathcal{H}}_{{\bf S}}$ the total energy, whose rate of change is dependent only on the boundary term ${\mathcal{B}}_{{\bf S},{\mathcal{H}}}$.

For system K, a similar property may be derived, i.e.,

[TABLE]

with

[TABLE]

2.3.2 Angular Momentum Conservation

In order to examine the conservation of angular momentum, it is useful to define, for any two-vector ${\bf f}=[f_{(1)},f_{(2)}]^{T}$, the operation $\tilde{}$ by $\tilde{{\bf f}}=[-f_{(2)},f_{(1)}]^{T}$. It should be clear that for any such vector ${\bf f}$, it must be true that

[TABLE]

Considering only the angular momentum of system S about the string axis (i.e., the $x$-component), left-multiply Eq. (4b) by \tilde{\mbox{\boldmath\eta}}^{T} and integrate over the interval $[0,1]$ to get

[TABLE]

Expanding the first term and integrating the second by parts gives

[TABLE]

where

[TABLE]

Finally, by applying identity (9) to the second and third terms under the integral above (note that $\psi$  is proportional to ${\bf q}$), Eq. (10) can be reduced to

[TABLE]

with

[TABLE]

which is the total angular momentum of system S, in the $x$ direction.

For system K, the analysis is similar, and gives

[TABLE]

with

[TABLE]

${\mathcal{A}}_{{\bf K}}$ is the total angular momentum of system K.

2.4 Boundary Conditions

In the interest of simplifying the analysis somewhat, and of rendering systems S and K conservative, it is useful to specify several boundary conditions. Fixed conditions, at an end of the string, are defined by

[TABLE]

and free conditions by

[TABLE]

Note that the free conditions above also imply that $\phi$ and $\psi$  vanish.

If either of these conditions holds at each of $x=0$ and $x=1$, then both ${\mathcal{B}}_{{\bf S},{\mathcal{H}}}$ and ${\mathcal{B}}_{{\bf S},{\mathcal{A}}}$ vanish, and as a result

[TABLE]

In other words, the energy and angular momentum of system S remain constant for all $t$ and equal to their initial values. Many other conditions, including mixtures of the above conditions, obviously lead to conservative behaviour as well.

For system K, of the two conditions given above, only the conditions (11) above is allowable (recall that $\xi=0$ at either end of the string is a starting point in the derivation of system K), and again leads, when enforced at both ends of the string, to

[TABLE]

2.5 Bounds on Solution Size

In the interest of simplifying notation, it is useful to introduce the spatial $L_{2}$ inner product and norm of two column vector functions of $x\in[0,1]$ and perhaps $t$, ${\bf f}$ and ${\bf g}$, containing the same number of elements. These are given by

[TABLE]

Obviously, such norms and inner products remain functions of the time variable $t$; when necessary, this will be indicated, e.g., $\|{\bf f}\|(t)$.

Returning to the expressions for kinetic and potential energy, and angular momentum of system S, it is then possible to rewrite them as

[TABLE]

${\mathcal{T}}_{{\bf S}}$ is non-negative, and ${\mathcal{V}}_{{\bf S}}$ is as well, under the condition

[TABLE]

which is the case for moderately stiff strings, and which will be assumed henceforth in this article.

Similarly, for system K, the kinetic and potential energies and angular momentum may be written as

[TABLE]

2.5.1 General Bounds

For system S, conservative boundary conditions lead to Eq. (13), and it is clearly true then that ${\mathcal{T}}_{{\bf S}}(t)\leq{\mathcal{H}}_{{\bf S}}(t)={\mathcal{H}}_{{\bf S}}(0)$, further implying that

[TABLE]

It is simple to translate the above bounds on the time derivatives of the displacements to bounds on the displacements themselves. For the longitudinal displacement, for example, one may write, employing the Cauchy-Schwartz inequality [27] and the first of the bounds (16),

[TABLE]

implying that

[TABLE]

In other words, growth of the $L_{2}$ norm of the longitudinal displacement is at most linear. An identical bound follows for the transverse displacement, i.e.,

[TABLE]

For system K, a similar bound on the transverse displacement holds, under conservative boundary conditions:

[TABLE]

2.5.2 Bounds under Fixed Conditions

For a given type of motion (i.e., longitudinal or transverse), if at least one of the ends is fixed, then better bounds are possible. For instance, consider condition (11) applied at $x=0$, paying special attention to the first condition, i.e., $\xi(0,t)=0$. The following bound is immediate:

[TABLE]

which implies, furthermore, that

[TABLE]

(If $\xi=0$ at both $x=0$ and $x=1$, then the above bound may be tightened to $\|\xi\|\leq\|p\|/2$.) Finally, using the fact that ${\mathcal{V}}_{{\bf S}}(t)\leq{\mathcal{H}}_{{\bf S}}(t)={\mathcal{H}}_{{\bf S}}(0)$, one may then conclude that

[TABLE]

Similar bounds may be found for $\eta$, again under fixed transverse conditions at at least one end of the string, for both system S and system K.

3 Grid Functions and Finite Difference Operators

In this section, a short review of grid functions as employed in finite difference schemes and the related difference operators is provided, with an eye toward applications in the study of schemes for conservative systems. Due to the vector nature of the differential equations to be studied here, definitions are framed here in terms of vector-valued grid functions, but it should be kept in mind that in most cases, such definitions reduce simply to the scalar case.

3.1 Grid Functions

A grid function ${\bf f}_{i}^{n}$ is defined as a column vector taking on values at the collection of points indexed by integers $i$ and $n\geq 0$; $({\bf f}_{i}^{n})^{T}$ is its transpose. It is intended as an approximation to a continuously-variable function ${\bf f}(x,t)$ at the location $x=ih_{x}$, $t=nh_{t}$, where $h_{x}$ and $h_{t}$ are the grid spacing and time step, respectively. In order that the connection with the underlying model problem be maintained, in this article, a grid function will be described using the same variable name as the continuous function it approximates. In addition, if a grid function is presented without one or both of its indices, it is assumed to have general indices $i$ and $n$.

As for the continuous case, for any column two-vector grid function ${\bf f}=[f_{(1)},f_{(2)}]$, the grid function $\tilde{{\bf f}}$ is defined by $\tilde{{\bf f}}=[-f_{(2)},f_{(1)}]$. It is then always true that, for any two two-vector grid functions ${\bf f}$ and ${\bf g}$,

[TABLE]

and, furthermore,

[TABLE]

where here, “0” is interpreted as a grid function taking on the value zero for all $i$ and $n$.

The important parameter $\lambda$, defined by

[TABLE]

is crucial in that numerical stability conditions are framed in terms of it, as per the linear case [28].

3.2 Temporal Operators

The basic temporal operators are the unit forward and backwards shifts, defined in terms of their effect on a grid function ${\bf f}$ by

[TABLE]

The forward, backward and central difference operators may be defined simply in terms of these shifts as

[TABLE]

(Here, the symbol 1 corresponds to the identity operation.) All of these serve as approximations to a first time derivative; a centered approximation to a second derivative is given by

[TABLE]

Forward, backward and central time-averaging operators, defined by

[TABLE]

are approximations to the identity operator. Another averaging operator, useful in the context of the construction of conservative schemes, is given by

[TABLE]

Note also that

[TABLE]

In the energetic analysis of difference schemes, the following identities are indispensable: for any grid function ${\bf f}$,

[TABLE]

where “$\star$” stands for any of “+”, “-” or “$\cdot$”, and

[TABLE]

The following identities are more useful when examining the conservation of angular momentum. For any two grid functions ${\bf f}$ and ${\bf g}$ of the same number of elements,

[TABLE]

Both are analogous to the product rule of differentiation. In particular, if for some two-vector grid function ${\bf q}$, it is true that ${\bf f}=\tilde{{\bf q}}$, and ${\bf g}=\delta_{t-}{\bf q}$, it then follows immediately that

[TABLE]

Another useful identity is the following:

[TABLE]

3.3 Spatial Operators

The only spatial difference operators which will be employed here are the forward and backward difference, defined here in terms of their action on the grid function ${\bf f}_{i}^{n}$ and by

[TABLE]

and are both approximations to a first spatial derivative; a centered approximation to the second derivative is then given by $\delta_{x+}\delta_{x-}$.

The operators $\delta_{t+}$, $\delta_{t-}$, $\delta_{to}$, $\mu_{t+}$, $\mu_{t-}$, $\mu_{to}$, $\delta_{x+}$, and $\delta_{x-}$ all commute.

3.4 Inner Products and Norms

The spatial inner product which will be of use in the present article is defined in terms of two grid functions ${\bf f}$ and ${\bf g}$, again of the same number of elements, and over the finite range of indices $i=r,\ldots,s$, for $r$, $s$ integer such that $s\geq r$:

[TABLE]

It is worth noting that for any vector grid functions ${\bf f}$ and ${\bf g}$, and any scalar grid function $l$,

[TABLE]

The definition of the norm follows in the usual way as

[TABLE]

The two standard inequalities which follow from the above definitions are the triangle inequality

[TABLE]

and the Cauchy-Schwartz inequality

[TABLE]

Summation by parts follows from the definition of the inner product as

[TABLE]

It is important to note that the identities defined in Section 3.2 involving the product of two grid functions extend immediately when an inner product is formed. For instance, identity (27) becomes

[TABLE]

when an inner product is taken over some spatial domain ${\mathcal{D}}$.

3.5 Bounds

It is straightforward to relate a bound on the norm of a time difference of a grid function to the norm of the grid function itself. For instance, consider a vector-valued grid function ${\bf f}_{i}^{n}$, and suppose it is true that

[TABLE]

over some spatial interval ${\mathcal{D}}$, for all $n$, and for some constant $K$. It then follows, from the triangle inequality 29, that

[TABLE]

Thus, given a bound such as (31), growth of the norm of a grid function is at most linear; this is independent of any boundary condition considerations.

The following bound on the norm of a grid function in terms of its spatial difference follows directly from the triangle inequality (29):

[TABLE]

If, for any scalar grid function $f_{i}$ (which could be a component of a vector grid function), it is true that $f_{r}=0$, then for any $r+1\leq m\leq s$, it is true that

[TABLE]

where “1” refers to a scalar grid function consisting of a sequence of ones. Then, by the Cauchy-Schwartz inequality,

[TABLE]

which implies, furthermore, that

[TABLE]

4 Finite Difference Schemes for Systems S and K

In Table 2, several finite difference schemes for systems S and K are presented; these are indicated by s and k, with a distinguishing superscript. The variety of schemes, in particular for system S, is intended to illustrate the many subtle differences among the schemes, with respect to conservation properties, numerical stability, as well as ease of implementation. The number of possible schemes, even of limited stencil, is of course much larger than that indicated here. Typically, in sound synthesis applications, explicit schemes similar to ${\bf s}^{(a)}$ are used [21].

All the schemes given in the table for system S have the same form, indicated in the second column of the table, which is a direct discretization of Eqs. (4), in $\phi$ and $\psi$. The shorthand forms

[TABLE]

are used throughout the rest of this article. Distinctions among the various schemes are due to variations in the discretization of $\phi$ and $\psi$, given explicitly in the third column of the table. Similarly, for system K, the two schemes given have the form of a direct discretization of Eq. (6), and variations are due to the way in which the quantity ${\mathcal{G}}$ is discretized.

All the schemes below are consistent with systems S or K and accurate to second order in both time and space (it is simplest to see this by rewriting the schemes in a first-order transmission-line form, as per [15]). They are two-step schemes, and for initialization, values of the grid functions $\xi^{n}$ and \mbox{\boldmath\eta}^{n} are required at the first two time steps, i.e., for $n=0$ and $n=1$. The spatial domain of the problem will be limited to $i\in{\mathcal{D}}=[0,\ldots,N]$. In dealing with boundary conditions, the set ${\mathcal{D}}^{+}=[1,\ldots,N]$ is also of great utility.

As mentioned above, each scheme can be considered from various points of view. Does it possess conserved analogues of angular momentum and energy? Are there simple conditions under which numerical stability can be ensured? Is the scheme explicit or implicit, and if implicit, do existence and uniqueness conditions for the numerical solution follow? All of these points will be dealt with in turn in the following sections.

4.1 Conservation of Angular Momentum

The conservation of angular momentum (in the $x$-direction) is perhaps the simplest property to examine. Consider first the simple scheme s(a)**, as given in Table 2. Taking the inner product of the second equation with \tilde{\mbox{\boldmath\eta}} over the domain ${\mathcal{D}}$ gives

[TABLE]

where the second and third equalities above follow from summation by parts (Eq. (30)) and identity (18) (note that \mbox{\boldmath\psi}_{{\bf s}^{(a)}} is a scalar multiple of ${\bf q}$). The boundary term $\mathcal{B}_{{\mathcal{A}},{\bf s}^{(a)}}$ is given in Table 3. From identity (26), it then follows that

[TABLE]

where

[TABLE]

can be identified easily with the angular momentum. Clearly, if the boundary term $\mathcal{B}_{{\mathcal{A}},{\bf s}^{(a)}}$ vanishes, then angular momentum is conserved by scheme ${\bf s}^{(a)}$.

The treatment of system ${\bf s}^{(b)}$ is similar, except that it is now necessary to take the inner product of the second equation with \mu_{to}\tilde{\mbox{\boldmath\eta}}, instead of \tilde{\mbox{\boldmath\eta}}, giving

[TABLE]

where the boundary term $\mathcal{B}_{{\mathcal{A}},{\bf s}^{(b)}}$ is given in Table 3. The expression for \mbox{\boldmath\psi}_{{\bf s}^{(b)}}, from Table 2, is made up of two terms; the first is a scalar multiple of ${\bf q}$, and the second a multiple of $\mu_{to}{\bf q}$. Thus, applying identity (18),

[TABLE]

where, in the second and third equalities above, commutativity of the operators $\mu_{to}$ and $\delta_{x-}$ and identity (27) have been used, respectively. Finally, applying identity (18), this can be rewritten as

[TABLE]

where

[TABLE]

This form of the angular momentum is distinct from the quantity conserved under scheme ${\bf s}^{(a)}$, but note that in the limit as $h_{t}$ becomes small, the two definitions approach one another. Scheme ${\bf s}^{(c)}$ does not possess a simple conserved quantity analogous to an angular momentum, but schemes ${\bf s}^{(d)}$ and ${\bf s}^{(e)}$ do—their analysis is nearly identical to that of ${\bf s}^{(a)}$ and ${\bf s}^{(b)}$, respectively, and the conserved quantities and boundary terms are given in Table 3.

For schemes ${\bf k}^{(a)}$ and ${\bf k}^{(b)}$, the analysis is very similar, and conservation of angular momentum holds in either case, with conserved quantities as given in Table 3.

4.2 Conservation of Energy

Conservation of energy for the schemes ${\bf s}^{(\bullet)}$ is always arrived at in the following way: given the general form of the scheme, shown in the second column of Table 2, take the inner product over the domain ${\mathcal{D}}$ of the first equation with $\delta_{to}\xi$, and the second with \delta_{to}\mbox{\boldmath\eta}. After using summation by parts (Eq. (30)) and adding the resulting equations one arrives at

[TABLE]

where $\phi_{{\bf s}^{(\bullet)}}$ and \mbox{\boldmath\psi}_{{\bf s}^{(\bullet)}} depend on the choice of scheme, and are given in Table 2. The boundary term ${\mathcal{B}}_{{\mathcal{H}},{\bf s}^{(\bullet)}}$ is of the same form for all the schemes given for system S:

[TABLE]

Clearly, the first term on the left hand side of Eq. (36) behaves as a difference approximation to the first derivative of the kinetic energy. The second two terms will vary from scheme to scheme, depending on the forms of $\phi_{{\bf s}^{(\bullet)}}$ and \mbox{\boldmath\psi}_{{\bf s}^{(\bullet)}}.

Scheme ${\bf s}^{(a)}$ does not possess a conserved energy, but the other four do. For example, for scheme ${\bf s}^{(b)}$, using the forms of $\phi_{{\bf s}^{(b)}}$ and \mbox{\boldmath\psi}_{{\bf s}^{(b)}} given in Table 2, one may write

[TABLE]

and

[TABLE]

and thus, combining Eqs. (4.2) and (4.2),

[TABLE]

where ${\mathcal{V}}_{{\bf s}^{(b)}}$ is as given in Table 4. Thus, one has, immediately,

[TABLE]

One may proceed in a similar vein for schemes ${\bf s}^{(c)}$, ${\bf s}^{(d)}$ and ${\bf s}^{(e)}$; all these algorithms are energy-conserving, with expressions for discrete kinetic and potential energy given in Table 4.

The two schemes for system K, ${\bf k}^{(a)}$ and ${\bf k}^{(b)}$, given in Table 2 differ only in the treatment of the quantity ${\mathcal{G}}$. For either scheme, one may take the inner product with \delta_{to}\mbox{\boldmath\eta} to get

[TABLE]

and, again using summation by parts (Eq. (30)),

[TABLE]

Here, the boundary term ${\mathcal{B}}_{{\mathcal{H}},{\bf k}^{(\bullet)}}$ is given by

[TABLE]

For scheme ${\bf k}^{(a)}$, substitution of the expression ${\mathcal{G}}_{{\bf k}^{(a)}}$ (given in Table 2) does not lead to an energy-conservation property. But for scheme ${\bf k}^{(b)}$, using ${\mathcal{G}}_{{\bf k}^{(b)}}$ (also given in Table 2) yields

[TABLE]

and, using identity (23),

[TABLE]

Finally, one arrives at

[TABLE]

where the expressions ${\mathcal{T}}_{{\bf k}^{(b)}}$, and ${\mathcal{V}}_{{\bf k}^{(b)}}$ are as given in Table 4. This is the desired discrete energy conservation property.

4.3 Conservative Boundary Conditions

Fixed boundary conditions, defined by

[TABLE]

are a direct counterpart to the continuous conditions (11), and are assumed to hold at an endpoint of the interval ${\mathcal{D}}$, i.e., for $i=0$ or $i=N$.

Free boundary conditions (12) at, e.g., the left end of the discrete domain, can be approximated by

[TABLE]

where it is recalled that for all the schemes discussed here, the definitions (35) hold. At the right end, such free conditions are given by

[TABLE]

For systems ${\bf k}^{(\bullet)}$, which depend only on $\eta$, the second of each pair of conditions above suffices to characterize a boundary condition as fixed or free.

Given the forms of ${\mathcal{B}}_{{\mathcal{A}},{\bf s}^{(\bullet)}}$, ${\mathcal{B}}_{{\mathcal{A}},{\bf k}^{(\bullet)}}$ given in Table 3, and ${\mathcal{B}}_{{\mathcal{H}},{\bf s}^{(\bullet)}}$ and ${\mathcal{B}}_{{\mathcal{H}},{\bf k}^{(\bullet)}}$ given by definitions (37) and (40), respectively, and recalling the various forms of $\phi_{{\bf s}^{(\bullet)}}$ and \mbox{\boldmath\psi}_{{\bf s}^{(\bullet)}} given in Table 2, it should be clear that a choice of a fixed discrete boundary condition (41), or a free boundary condition such as (42) or (43) at each of the endpoints of the domain ${\mathcal{D}}$ leads to a vanishing of the boundary terms. Under such conditions, any such scheme is fully conservative, i.e., one has

[TABLE]

and thus

[TABLE]

Such conservative boundary conditions will be assumed for the remainder of this article.

4.4 Numerical Stability

Under further conditions, numerical stability of the difference schemes may follow immediately from the discrete conserved energy quantities given in Table 4. (The schemes which do not possess an energy conservation property, namely ${\bf s}^{(a)}$ and ${\bf k}^{(a)}$, will not be examined in this section.) The main goal, in this section, is to find conditions, if they exist, under which the discrete conserved energy is positive for all possible choices of the state variables $\xi$ and $\eta$. If such conditions exist, then bounds on the solution size in terms of initial conditions (and hence a numerical stability guarantee) can be obtained.

4.4.1 Positivity of Discrete Conserved Energy

For the schemes ${\bf s}^{(b)}$, ${\bf s}^{(c)}$, ${\bf s}^{(d)}$, and ${\bf s}^{(e)}$, a discrete energy conservation property exists, and in particular, under conservative boundary conditions such as those discussed in Section 4.3, one has

[TABLE]

It is useful to begin with scheme ${\bf s}^{(d)}$, which has a conserved energy of a particularly simple form. Considering the form for ${\mathcal{V}}_{{\bf s}^{(d)}}$ given in Table 4, which is

[TABLE]

one may then write, employing identity (24),

[TABLE]

Applying definitions (35) and the bound (33) yields the following inequality

[TABLE]

Note that the parameter $\lambda$, defined in Eq. (19), has been introduced here. One further has

[TABLE]

Keeping in mind the condition (15), which is assumed a priori, then the discrete conserved energy ${\mathcal{H}}_{{\bf s}^{(d)}}$ will be non-negative under the conditions

[TABLE]

which have the form of Courant-Friedrichs-Lewy type conditions [28], [29], which often result from a Fourier or von Neumann type analysis of difference schemes in the linear case. These conditions, and similar conditions for schemes to be discussed shortly, are given in the second column of Table 5.

Referring to Table 4, the expression ${\mathcal{V}}_{{\bf s}^{(c)}}$ for the discrete potential energy for scheme ${\bf s}^{(c)}$ differs from ${\mathcal{V}}_{{\bf s}^{(d)}}$ only by a single term, which is non-negative; following steps similar to the above, one may easily derive the following inequality:

[TABLE]

Thus the positivity conditions (48) derived above for scheme ${\bf s}^{(d)}$ hold for scheme ${\bf s}^{(c)}$ as well.

The analysis for scheme ${\bf s}^{(e)}$ is similar to that performed above for schemes ${\bf s}^{(d)}$, and ${\bf s}^{(c)}$, but the positivity conditions are slightly different. One may derive the following inequality:

[TABLE]

For positivity, condition (48b) holds as before, but condition (48a) must be modified to

[TABLE]

For schemes ${\bf s}^{(c)}$, ${\bf s}^{(d)}$, and ${\bf s}^{(e)}$, the analysis above is simple, because the contributions of the nonlinearity to the expressions for the discrete conserved energy (i.e., the terms which are not simply quadratic in the expressions for ${\mathcal{V}}_{{\bf s}^{(c)}}$, ${\mathcal{V}}_{{\bf s}^{(d)}}$ and ${\mathcal{V}}_{{\bf s}^{(e)}}$) are themselves non-negative. Thus the determination of positivity conditions reduces, essentially, to analysis of a linear problem. For scheme ${\bf s}^{(b)}$, however, this is not the case; the term in the expression for ${\mathcal{V}}_{{\bf s}^{(b)}}$ resulting from the nonlinearity is not necessarily positive, and indeed, can be negative and unbounded (the choices $p=0$ and $e_{t-}p=-e_{t-}{\bf q}^{T}{\bf q}$ illustrate this point simply). No global positivity condition is available for scheme ${\bf s}^{(b)}$.

For the energy-conserving scheme ${\bf k}^{(b)}$, similar positivity conditions are also immediate. Beginning from the expression for ${\mathcal{V}}_{{\bf k}^{(b)}}$, one may write, employing identity (24) and inequality (33),

[TABLE]

implying

[TABLE]

The total energy is then positive under condition (48a).

4.4.2 Bounds on Solution Size

In this section, only schemes with a conditionally positive discrete conserved energy are considered, namely schemes ${\bf s}^{(c)}$, ${\bf s}^{(d)}$, ${\bf s}^{(e)}$ and ${\bf k}^{(b)}$, under conservative boundary conditions.

Consider first the inequalities (49) and (47), which give lower bounds on the conserved energy for schemes ${\bf s}^{(c)}$ and ${\bf s}^{(d)}$. Under stability conditions (48), it is then clear that, for either scheme,

[TABLE]

For scheme ${\bf s}^{(e)}$, the bounds, under stability conditions (48b) and (51), are slightly different; now, one has

[TABLE]

For scheme ${\bf k}^{(b)}$, the bound on $\eta$  follows immediately as

[TABLE]

Any of these bounds on the norm of a time difference of a grid function may be converted to a general bound on the norm of the grid function itself, through an application of inequality (32). For example, for the bounds given above on $\delta_{t-}\xi$ in (53) for schemes ${\bf s}^{(c)}$ or ${\bf s}^{(d)}$, it then follows that

[TABLE]

Thus growth of the longitudinal displacement is at most linear in time, for any conservative boundary conditions. A similar bound can be found for the transverse displacement, using the second of bounds (53); bounds for the schemes ${\bf s}^{(e)}$ and ${\bf k}^{(b)}$, arrived at in a nearly identical way, are given in Table 5.

Under fixed boundary conditions, tighter bounds may be obtained. For schemes ${\bf s}^{(c)}$, ${\bf s}^{(d)}$, and ${\bf s}^{(e)}$, the following bounds on the quantities $p$ and ${\bf q}$ follow from the appropriate stability conditions ((48) for ${\bf s}^{(c)}$ and ${\bf s}^{(d)}$, and (48b) and (51) for ${\bf s}^{(e)}$) and inequalities (49), (47) and (50):

[TABLE]

For scheme ${\bf k}^{(b)}$, the bound, from stability condition (48b) and inequality (52), is simply

[TABLE]

For any of the above bounds on $\|\mu_{t-}p\|_{{\mathcal{D}}^{+}}$ or $\|\mu_{t-}{\bf q}\|_{{\mathcal{D}}^{+}}$, bounds on $\|\xi\|_{{\mathcal{D}}}$ or \|\mbox{\boldmath\eta}\|_{{\mathcal{D}}} may be obtained in the following way. Considering, for example, $\eta$ , one may immediately write, using identity (22) and the triangle inequality (29),

[TABLE]

and, furthermore, under scheme ${\bf s}^{(c)}$ for example,

[TABLE]

where in the first inequality above, the inequality (34) and the second of the bounds (53) have been employed, and in the second, the second of the bounds (55) has been used.

This bound, and similar forms for the other schemes of interest in this section are given in Table 5.

4.5 Implementation Details

It is useful to write the schemes defined in Table 2 in forms suitable for computer implementation. In order to avoid presenting a multiplicity of slightly different cases, it is assumed in this section that boundary conditions are of the fixed type (i.e., (41) for system s) at both ends of the string, which are the conditions of interest in many cases. In this case, the values of \mbox{\boldmath\eta}_{i}^{n} and $\xi_{i}^{n}$ at the endpoints of the spatial domain (i.e., at $i=0$ and $i=N$) are fixed at zero, and need not be considered. It is helpful to introduce the column vectors ${\bf u}^{n}$, ${\bf v}_{(1)}^{n}$ and ${\bf v}_{(2)}^{n}$ defined by

[TABLE]

and the state vector ${\bf w}^{n}$ containing all the displacements at time step $n$, defined by

[TABLE]

Schemes ${\bf s}^{(a)}$, ${\bf s}^{(b)}$, ${\bf s}^{(c)}$, and ${\bf s}^{(d)}$ may all be written in the form of two-step matrix recursions of the form

[TABLE]

where ${\bf A}_{{\bf s}^{(\bullet)}}^{n}$ and ${\bf B}_{{\bf s}^{(\bullet)}}^{n}$ are square matrices which depend on the values of the state at time step $n$, i.e., ${\bf w}^{n}$. The fact that ${\bf A}_{{\bf s}^{(\bullet)}}^{n}$ must be inverted (or, rather, a linear system solved) reflects the implicit nature of the schemes. For system ${\bf s}^{(a)}$, the matrix ${\bf A}_{{\bf s}^{(\bullet)}}^{n}$ is diagonal, and the algorithm is explicit. The complete forms of the state update matrices are given in Table 6, where the scaled differentiation matrix ${\bf D}$ is defined as

[TABLE]

${\bf D}$ is an $N$ by $N-1$ matrix. In addition, it is useful to define the diagonal $N$ by $N$ matrices ${\bf P}$, ${\bf Q}_{(1)}$ and ${\bf Q}_{(2)}$ by

[TABLE]

In Table 6, ${\bf I}$ refers to the $N-1$ by $N-1$ identity matrix, and the constant $\beta$ is defined as

[TABLE]

If a scheme may be written in a state update form as per Eq. (56), one has an immediate proof of existence and uniqueness of solutions, and this is indeed the case for schemes ${\bf s}^{(a)}$, ${\bf s}^{(b)}$, ${\bf s}^{(c)}$, and ${\bf s}^{(d)}$. In all these cases, although the difference scheme is, as a whole, strongly nonlinear and implicit, the coupling among the state variables at the current time step is linear in character. It is important to note, however, that this property is independent of the existence of discrete conservation laws—for scheme ${\bf s}^{(e)}$, which is energy conserving, there is no state update form. Thus a proof of conditions for existence and uniqueness for solutions to ${\bf s}^{(e)}$ will be considerably more difficult, if even possible to obtain. It remains true, however, that if a solution does exist, it will be stable under condition (48).

For scheme ${\bf k}^{(a)}$, which is explicit, the implementation is immediate. For scheme ${\bf k}^{(b)}$, however, the situation is slightly different. Though it might appear, from the form of ${\mathcal{G}}_{{\bf k}^{(b)}}$, that the scheme would be implicit, it is in fact possible to write it in an explicit form in the following way. For the quantity ${\mathcal{G}}_{{\bf k}^{(b)}}$, one may proceed as follows:

[TABLE]

where in the second, third and fourth equalities, identity (25), the definition of system ${\bf k}^{(b)}$ from Table 2, and identity (21) have been used, respectively. and thus

[TABLE]

which renders scheme ${\bf k}^{(b)}$ fully explicit. If the boundary conditions are fixed or free, one may go further and write

[TABLE]

5 Numerical Examples

5.1 Energy and Angular Momentum Conservation

As a basic illustration of the conservative properties of the schemes presented here, consider the case of a steel string string, of parameters given in the caption to Fig. 2. The initial conditions of the string are

[TABLE]

and are discretized as

[TABLE]

In other words, the string is subjected to an initial displacement in the form of the first linear modal configuration in the polarization corresponding to $\eta_{(1)}$, and an initial velocity of a similar form in the other polarization. The initial longitudinal displacement and velocity is assumed to be zero. Shown in Fig. 2 are the results of simulations using scheme ${\bf s}^{(a)}$, where the transverse string displacements $\eta_{(1)}$ and $\eta_{(2)}$ at the string center are plotted as a function of time in each case. In this, and all examples in this section, the parameter $\lambda$ is chosen as close to the bounds (48b) as possible.

Tables 7 and 8 illustrate the numerical conservation properties of those algorithms which possess discrete conserved angular momentum and energy, respectively. The fluctuations observed in some cases (typically in the 12th place) are due to numerical round-off error.

As an example of the effect of the nonlinearity, consider a representative scheme, such as, e.g., ${\bf s}^{(d)}$, which is both angular momentum and energy conserving, applied to the string of parameters as given in the caption to Fig. 2, and again with initial conditions as given by Eqs. (58). In Fig. 3 are shown transverse displacements at the string center, for different values of $\gamma_{1}$ and $\gamma_{2}$. Notice in particular the change in the oscillating frequency. Values of the conserved quantities ${\mathcal{A}}_{{\bf s}^{(d)}}$ and ${\mathcal{H}}_{{\bf s}^{(d)}}$ are given in tables 9 and 10, respectively.

5.2 Instability of Planar Motion

The phenomenon of the instability of purely planar motion is inherent to nonlinear strings, and has been examined in depth by various authors [30], [31], [32], [33]. (By “instability,” one refers here to the tendency for motion which is confined to a single plane to be transferred to the perpendicular polarization, and not instability in the sense of explosive growth of the solution.)

Consider a string under fixed conditions, with numerical initial conditions of the form

[TABLE]

where $\theta_{i}$ is a uniformly distributed random variable taking values over $(-1,1)$. In this case, the initial conditions correspond to a stationary modal distribution for $\eta_{(1)}$, of amplitude $\gamma_{1}$ and to a random perturbation for $\eta_{(2)}$, of amplitude $\gamma_{(2)}$. When $\gamma_{2}=0$, the state of the string remains in the $\eta_{(1)}$ polarization for all future time. But when $\gamma_{2}$ is not identically zero, even if very small, energy will be transferred from the $\eta_{(1)}$ polarization to the $\eta_{(2)}$ polarization, and whirling will occur. In Figs. 4 and 5 are shown simulation results, again for the string of parameters given in the caption to Fig. 2, under the initial conditions (60), using a very small value of $\gamma_{2}$ relative to $\gamma_{1}$ (values given in the captions). Both polarizations are shown in the figures. In the case of schemes ${\bf s}^{(\bullet)}$, energy is transferred nearly entirely from the $\eta_{(1)}$ polarization to the $\eta_{(2)}$, over approximately 100 000 time steps, and continues to oscillate henceforth. For schemes ${\bf k}^{(\bullet)}$, the transfer is much faster, and at the same time, limited, in the sense that the total energy transferred is much smaller. This planar instability, while physical, does not lead to explosive growth (indeed, angular momentum and energy conservation properties, where possessed by a given scheme, are not violated, as are stability conditions when implied by the latter form of conservation), but to a great deal of variation in the results of different schemes, particularly for schemes ${\bf s}^{(\bullet)}$.

5.3 Numerical Instability of Scheme ${\bf s}^{(b)}$

Scheme ${\bf s}^{(b)}$ is perhaps the most interesting discussed here, in that it does indeed possess a discrete conserved energy, but this property does not lead to simple stability conditions. Instability which develops is inherently tied to ill-conditioning of the matrices${\bf A}_{{\bf s}^{(b)}}$ and ${\bf B}_{{\bf s}^{(b)}}$, from the matrix update form Eq. (56), which leads to large numerical round off error. (Note in particular that of the matrices shown in Table 6, only ${\bf B}_{{\bf s}^{(b)}}$ is asymmetric.).

Under low amplitude initial conditions, scheme ${\bf s}^{(b)}$ performs similarly to the other schemes for system ${\bf S}$, as discussed in Section 5.1. Consider, though, planar initial conditions (58) with $\gamma_{1}=0.1$ and $\gamma_{2}=0$. Though it is difficult to see in Figure 6(b), the energy remains conserved and constant for several time steps, but then suddenly begins to fluctuate, and quickly begins taking on negative values. The displacement itself, shown in the left panel, experiences quantized jumps in amplitude, the largest coinciding with the moment at which the energy begins taking on negative values. If the simulation is allowed to proceed, similar discrete jumps in amplitude continue to occur; this type of instability is obviously of a very different nature from the exponential growth often seen in numerical schemes. Note in particular the interesting quantization effects in the fluctuations in energy.

The schemes ${\bf s}^{(c)}$ and ${\bf s}^{(d)}$ are well-behaved under these conditions.

6 Conclusions

The main topic of this article has been the construction of numerical schemes for nonlinear strings which possess conservation properties; several different schemes for two distinct string models (systems ${\bf S}$ and ${\bf K}$) have been presented, which serve to highlight various distinctions among the schemes, which are summarized in Table 11. The angular momentum and energy conservation properties are generally independent; a scheme may possess one or the other, or both, or (though an example has not been presented here) neither. If energy conservation is implied by a scheme, stability conditions may follow, if the expression for energy can be shown to be generally positive for at least some range of choices of the material parameters and $\lambda$, the Courant number as defined by Eq. (19). But, as in the case of scheme ${\bf s}^{(b)}$, this does not necessarily follow from energy conservation. Finally, there is the issue of computability, which is also independent of the other properties. Certain schemes are fully explicit (such as ${\bf s}^{(a)}$, ${\bf k}^{(a)}$ and ${\bf k}^{(b)}$), and others (namely ${\bf s}^{(b)}$, ${\bf s}^{(c)}$, ${\bf s}^{(d)}$ and ${\bf s}^{(e)}$) are formally implicit. There is a distinction to be made here, however, between schemes ${\bf s}^{(b)}$, ${\bf s}^{(c)}$, and ${\bf s}^{(d)}$, for which the implicit nature of the scheme is manifested as a linear system inversion performed at each time step, and ${\bf s}^{(e)}$, for which this is not possible. Thus, for this former set of schemes, there is at hand a simple proof of existence and uniqueness of solutions (i.e., the matrix update form Eq. (56)), whereas for the latter, there is not, and an implementation will require iterative methods for solution. It is important to reiterate that the simple implementation property of the former set of schemes is dependent on the form of the nonlinearity, which contains terms up to order three; for more general nonlinear forms (such as the exact nonlinear string system; see, e.g., [22]), this is not the case. It is also tempting to conclude that, at least for a sufficiently complex nonlinear system such as ${\bf S}$, there is no fully explicit scheme which preserves energy; the same is clearly not true for the simpler system ${\bf K}$, for which the conservative scheme ${\bf k}^{(b)}$ may be written in an explicit form.

Several interesting features of these methods were examined in numerical examples in Section 5. Planar instability of string motion was illustrated in Section 5.2; under perturbed planar initial conditions, it is difficult to obtain consistent results across various different schemes for the same system; it may be possible to relate the resulting rate of oscillation to conservation of angular momentum, though this is a large topic, to be left to a future work. The largest issue, however, is that of discrete conservation itself, under finite-precision machine arithmetic. As was illustrated in Section 5.1, energy and angular momentum may be conserved in a scheme exactly, in infinite precision, but when numerical round-off errors occur, conservation is lost. Under moderate vibration amplitudes, the resulting fluctuations will be small, often on the order of “machine epsilon,” but as the nonlinearity becomes stronger, these fluctuations can become large, even in floating point arithmetic, and numerical stability may be violated. Perhaps the best way of examining the effect of round-off will be through an examination of the matrix update form Eq. (56), though this is a large separate topic which cannot be entered into in any detail here. The treatment of simulation techniques under finite-precision arithmetic is given scant attention in the literature, which is surprising; it is perhaps worth mentioning that this issue has been dealt with in the simulation of electrical networks (and also digital filtering), in particular by Fettweis (in the context of wave digital filters [34]) and Smith (for transmission lines and digital waveguides [19]), essentially through the decomposition of conservative numerical methods into unitary matrix transformations.

Finally, though the focus of this article has been on conservative methods for nonlinear string vibration, it should be clear that there are various more complex systems with nonlinearities of a similar simple form (i.e., resulting from a low-order series expansion); chief among these are the Berger and von Karman models of plate vibration [13], [35], which are, in a sense, direct analogues to systems ${\bf K}$ and ${\bf S}$ for the string, respectively, with added fourth-order terms modelling stiffness. The extension of the methods discussed here to these systems is immediate, and investigation is currently under way.

Appendix A Linear Damping

The introduction of linear damping terms to the systems under study here affects the analysis of difference schemes in only a very minor way. If such terms are included in system ${\bf S}$, Eqs. (4) will be modified as

[TABLE]

where $\sigma_{\xi}$ and \sigma_{{\scriptsize\mbox{\boldmath\eta}}} are non-negative constants, and where $\phi$ and $\psi$ are as defined in Eq. (4). The energetic analysis is similar to before, except that one will have

[TABLE]

and if the boundary term ${\mathcal{B}}_{{\bf S}}$ vanishes, one has

[TABLE]

Thus all bounds on solution size, discussed in Section 2.5, remain unchanged.

The conservation of angular momentum ${\mathcal{A}}_{{\bf S}}$, defined as before, is generalized to

[TABLE]

In this case, under conservative boundary conditions ${\mathcal{B}}_{{\bf S},{\mathcal{A}}}=0$, one then has

[TABLE]

In other words, the angular momentum decays exponentially.

Consider one of the difference schemes ${\bf s}^{(c)}$, ${\bf s}^{(d)}$ or ${\bf s}^{(e)}$, which is both energy conserving and for which CFL-like conditions exist for the positivity of the discrete energy. For any of these schemes, an extension to the lossy case may be discretized as

[TABLE]

Discrete energy conservation generalizes simply to

[TABLE]

and if the boundary terms vanish, one then has

[TABLE]

and all derived bounds on solution size again hold as before.

For schemes ${\bf s}^{(a)}$ and ${\bf s}^{(d)}$, it is simple to show that

[TABLE]

and if the boundary terms vanish, one has a geometric decay in the angular momentum of the form

[TABLE]

(Notice, however, that if

[TABLE]

the decay will be highly oscillatory, and unphysical.) For schemes ${\bf s}^{(b)}$ and ${\bf s}^{(e)}$, such exact geometric decay does not follow.

When a linear damping term is added to system ${\bf K}$, it is easy to show that a similar monotonic decrease in total energy, as well as an exponential decay in angular momentum also follows. For both schemes ${\bf k}^{(a)}$ and ${\bf k}^{(b)}$, when the linear damping term is added, similarly to for schemes ${\bf s}^{(\bullet)}$ above, angular momentum will decrease geometrically, and for scheme ${\bf k}^{(b)}$, discrete energy will decrease monotonically.

Appendix B Improved Stability Conditions for Schemes ${\bf s}^{(\bullet)}$

For any of the schemes for system ${\bf S}$ which are conservative and for which a global stability condition may be derived (namely schemes ${\bf s}^{(c)}$, ${\bf s}^{(d)}$, and ${\bf s}^{(e)}$), it is true that the condition on the time step $h_{t}$, for a given grid spacing $h_{x}$ can be exceedingly small, from conditions such as (48) and (51). On the other hand, these schemes are already implicit, and there is thus no loss in efficiency in generalizing the treatment of the linear part of system ${\bf S}$ to improve on these strict bounds on $h_{t}$.

Turning attention to Table 2, any of schemes ${\bf s}^{(\bullet)}$ may be generalized directly as

[TABLE]

where $\phi_{{\bf s}^{(\bullet)},\tau}$ and \mbox{\boldmath\psi}_{{\bf s}^{(\bullet)},\nu} are defined in terms of the free parameters $\tau$ and $\nu$ as

[TABLE]

Schemes (66) reduce to the forms shown in Table 2 when $\tau=\nu=0$.

Consider first the energetic analysis of the generalized systems (66); this is essentially the same as that carried out in Section 4.2, except for the new linear terms parameterized by $\tau$ and $\nu$. For those schemes with an energy conservation property, namely ${\bf s}^{(b)}$, ${\bf s}^{(c)}$, ${\bf s}^{(d)}$ and ${\bf s}^{(e)}$, the potential energy under this generalization is modified to

[TABLE]

In other words, for positive $\tau$ and $\nu$, the additional contribution to the potential energy is positive.

Consider now, as an example, the effect on scheme ${\bf s}^{(d)}$. The potential energy is generalized from the expression (44) to

[TABLE]

The expression for the kinetic energy ${\mathcal{T}}_{{\bf s}^{(\bullet)}}$ as given in Table 4 remains unchanged. Clearly, if $\tau\geq 1$ and $\nu\geq 1$, then the potential energy is always positive, as will be the energy ${\mathcal{H}}_{{\bf s}^{(d)},\tau,\nu}$. Under these conditions, the generalized algorithm is stable for any choice of time step $h_{t}$. If $\tau\leq 1$, or $\nu\leq 1$, the analysis is very similar to that carried out previously. One may obtain similar improvements on the bound on the time step for schemes ${\bf s}^{(d)}$ and ${\bf s}^{(e)}$, as well as ${\bf k}^{(a)}$ and ${\bf k}^{(b)}$; scheme ${\bf s}^{(b)}$ can still be shown to have a conserved energy which is non-positive for at least some choices of state variables.

As mentioned earlier in this section, the generalized schemes will not require more operations per time step; some modifications to the update matrices given in Section 4.5 will be necessary, but the sparsity remains nearly the same.

Appendix C Energy-conserving Spectral Method for System K

Due to the special form of the nonlinearity in system ${\bf K}$, an alternative analysis is possible using spatial Fourier series expansion techniques [2], [36], [30]; such analysis leads naturally to the construction of highly accurate spectral-type numerical solution methods [37], [38], which like the simpler difference schemes discussed in the main body of this article, are conservative. The same is not true for system ${\bf S}$.

Consider system ${\bf K}$ under fixed boundary conditions. An expansion for $\eta$  of the form

[TABLE]

where the time-dependent vector expansion coefficients are given by \hat{\mbox{\boldmath\eta}}_{k}(t), for $k=1,\ldots,\infty$, thus satisfies the boundary conditions automatically. One may then rewrite system ${\bf K}$ as the infinite system of ordinary differential equations

[TABLE]

This may be time-discretized immediately, and the infinite system of equations truncated to $M$ terms to form system ${\bf k}^{(s)}$, defined as

[TABLE]

where the form of ${\mathcal{G}}_{{\bf k}^{(s)}}$ under discretization is left unspecified for the moment. Notice that the approximation to the spatial derivative operators above is spectrally accurate, and is exact in the limit as $M$ becomes large.

Introducing the inner product of two sets of vector expansion coefficients $\hat{\bf f}_{k}$ and $\hat{\bf g}_{k}$ of dimension $M$ by

[TABLE]

and the associated norm by

[TABLE]

one may then take the inner product of Eq. (67) with \delta_{to}\hat{\mbox{\boldmath\eta}} to get

[TABLE]

or

[TABLE]

Now, in analogy with scheme ${\bf k}^{(b)}$, one may define ${\mathcal{G}}_{{\bf k}^{(s)}}$ as

[TABLE]

and one again arrives at an expression for conserved energy as

[TABLE]

with

[TABLE]

Conditions for positivity may be arrived at by rewriting the expression for potential energy as

[TABLE]

which gives a lower bound for the total energy as

[TABLE]

and the positivity condition is easily read off as

[TABLE]

Given this positivity condition, bounds on the solution size (i.e., bounds on the norm of \hat{\mbox{\boldmath\eta}}) may be derived exactly as in Section 2.5. Bounds on \hat{\mbox{\boldmath\eta}} may be simply related to bounds on $\eta$  itself through an application of Parseval’s Theorem [27]. Conservation of angular momentum also holds for system ${\bf k}^{(s)}$, and is trivial to show, through an inner product of system (67) with \tilde{\hat{\mbox{\boldmath\eta}}}.