◁ The Fortnightly Tambourine

ideal_string

Ideal String Equation

\[\frac{\partial^{2}}{\partial t^{2}} y{\left(x,t \right)} = c^{2} \frac{\partial^{2}}{\partial x^{2}} y{\left(x,t \right)} \]

Where:

• \(y{\left(x,t \right)} \left(\text{m}\right)\): vertical displacement at position \(x\) and time \(t\)
• \(c \left(\frac{\text{m}}{\text{s}}\right)=\sqrt{\frac{\text{T}_0}{\mu}}=\sqrt{\frac{\text{T}_0}{\rho \text{A}}}=2 \text{L} f_{0}\): speed of the wave in the string
• \(\text{T}_0 \left(\text{N}\right)=\rho \text{A} \left(2 \text{L} f_{0}\right)^{2}\): Tension on the string
• \(\mu \left(\frac{\text{kg}}{\text{m}}\right)=\rho \text{A}\): linear density of the string
• \(\rho \left(\frac{\text{kg}}{\text{m}^{3}}\right)\): mass density of the string
• \(\text{A} \left(\text{m}^{2}\right)\): cross-sectional area of the string
• \(\text{E} \left(\frac{\text{N}}{\text{m}^{2}}\right)\): Elastic, Tensile or Young’s modulus of the string / bar material. Setting \(\text{E} = 0\) leads to no stiffness.
• \(\text{I} \left(\text{m}^{4}\right)=\frac{\pi r^{4}}{4}\): Second moment of area (area moment of inertia) of the string’s cross-section
• \(\kappa \left(\frac{\text{m}^{2}}{\text{s}}\right)=\sqrt{\frac{\text{E} \text{I}}{\mu}}=\sqrt{\frac{\text{E} \text{I}}{\rho \text{A}}}=\frac{\sqrt{\text{E}} r}{2 \sqrt{\rho}}\): bar stiffness
• \(\beta_{1} \left(\frac{1}{\text{s}}\right)\): damping coefficient 01
• \(\beta_{2} \left(\text{s}\right)\): damping coefficient 02
• Units:
  • \(\text{N} \left(\text{ newton }\right)\)
  • \(\text{kg} \left(\text{ kilogram }\right)\)
  • \(\text{m} \left(\text{ meter }\right)\)
  • \(\text{s} \left(\text{ second }\right)\)

Difference Equation

This equation approximates the continuous partial derivatives using finite differences:

\[- \frac{2 y{\left(x,t \right)}}{{\Delta}t^{2}} + \frac{y{\left(x,t - {\Delta}t \right)}}{{\Delta}t^{2}} + \frac{y{\left(x,t + {\Delta}t \right)}}{{\Delta}t^{2}} = c^{2} \left(- \frac{2 y{\left(x,t \right)}}{{\Delta}x^{2}} + \frac{y{\left(x - {\Delta}x,t \right)}}{{\Delta}x^{2}} + \frac{y{\left(x + {\Delta}x,t \right)}}{{\Delta}x^{2}}\right) \]

Discretizing with the following substitutions:

• Continuous time \(t\) becomes discrete time steps \(t=\text{n} {\Delta}t,\) where \(\text{n}=0,1,2,\cdots,\text{N}-1\)
• Continuous position \(x\) becomes discrete grid points \(x=\text{i} {\Delta}x,\) where \(\text{i}=0,1,2,\cdots,\text{M}-1\)

Applying these substitutions, the difference equation becomes a discrete equation relating values at grid points:

\[\frac{{y}_{\text{i},\text{n} + 1}}{{\Delta}t^{2}} + \frac{{y}_{\text{i},\text{n} - 1}}{{\Delta}t^{2}} - \frac{2 {y}_{\text{i},\text{n}}}{{\Delta}t^{2}} = c^{2} \left(\frac{{y}_{\text{i} + 1,\text{n}}}{{\Delta}x^{2}} + \frac{{y}_{\text{i} - 1,\text{n}}}{{\Delta}x^{2}} - \frac{2 {y}_{\text{i},\text{n}}}{{\Delta}x^{2}}\right) \]

Solving the discrete equation for the value at the next time step, \({y}_{\text{i},\text{n} + 1}\), yields the explicit finite difference scheme:

\[{y}_{\text{i},\text{n} + 1} = \frac{c^{2} {\Delta}t^{2} \left({y}_{\text{i} + 1,\text{n}} + {y}_{\text{i} - 1,\text{n}} - 2 {y}_{\text{i},\text{n}}\right)}{{\Delta}x^{2}} - {y}_{\text{i},\text{n} - 1} + 2 {y}_{\text{i},\text{n}} \]

Independent Variables

We want to control the simulation specifying the following variables:

• \(f_{0} \left(\text{Hz}\right)\): natural frequency of the string or bar
• \(f_{s} \left(\text{Hz}\right)\): sampling frequency
• \(\text{L} \left(\text{m}\right)\): length of the string or bar
• \(d \left(\text{s}\right)\): duration of the simulation
• Units:
  • \(\text{Hz} \left(\text{ hertz }\right)\)
  • \(\text{m} \left(\text{ meter }\right)\)
  • \(\text{s} \left(\text{ second }\right)\)

Dependent Variables

These are derived from the independent variables and stability requirements:

• \(\text{N} \left(\text{unitless}\right)=d f_{s}\): number of time steps
• \(\text{M} \left(\text{unitless}\right)=\frac{\text{L}}{{\Delta}x}\): number of space steps
• \({\Delta}t \left(\text{s}\right)=\frac{d}{\text{N}}=\frac{1}{f_{s}}\): Delta time (time step size)
• \({\Delta}x \left(\text{m}\right)=\frac{\text{L}}{\text{M}}\): Delta length (spatial step size)
• Units:
  • \(\text{m} \left(\text{ meter }\right)\)
  • \(\text{s} \left(\text{ second }\right)\)

\({\Delta}x\):

Determining the Spatial Step using Von Neumann Stability Analysis To ensure the numerical simulation doesn’t produce spurious growing oscillations (instability), we perform a Von Neumann stability analysis. We substitute a test solution (ansatz) representing a single Fourier mode into the finite difference scheme.

We use the wave-form ansatz: \({y}_{\text{i},\text{n}} = \text{Y} e^{i \left(\text{i} k {\Delta}x + \text{n} w {\Delta}t\right)}\)

where \(k\) is the wavenumber and \(w\) is the angular frequency.

Substituting the ansatz into the finite difference scheme gives:

\[\text{Y} e^{i w {\Delta}t} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t} = \frac{\text{Y} c^{2} {\Delta}t^{2} e^{i k {\Delta}x} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t}}{{\Delta}x^{2}} - \frac{2 \text{Y} c^{2} {\Delta}t^{2} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t}}{{\Delta}x^{2}} + \frac{\text{Y} c^{2} {\Delta}t^{2} e^{- i k {\Delta}x} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t}}{{\Delta}x^{2}} + 2 \text{Y} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t} - \text{Y} e^{- i w {\Delta}t} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t} \]

Dividing by the common factor \(\text{Y} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t}\), rearranging terms and simplifying leads to a relationship between the frequency \(w\) and the wavenumber \(k\) (the numerical dispersion relation):

\[e^{i w {\Delta}t} + e^{- i w {\Delta}t} = \frac{c^{2} {\Delta}t^{2} \left(e^{i k {\Delta}x} - 2 + e^{- i k {\Delta}x}\right)}{{\Delta}x^{2}} + 2 \]

\[\cos{\left(w {\Delta}t \right)} = \frac{c^{2} {\Delta}t^{2} \left(\cos{\left(k {\Delta}x \right)} - 1\right) + {\Delta}x^{2}}{{\Delta}x^{2}} \]

Solving for the angular frequency \(w\):

\[w = \left[ \frac{- \operatorname{acos}{\left(\frac{c^{2} {\Delta}t^{2} \cos{\left(k {\Delta}x \right)} - c^{2} {\Delta}t^{2} + {\Delta}x^{2}}{{\Delta}x^{2}} \right)} + 2 \pi}{{\Delta}t}, \ \frac{\operatorname{acos}{\left(\frac{c^{2} {\Delta}t^{2} \cos{\left(k {\Delta}x \right)} - c^{2} {\Delta}t^{2} + {\Delta}x^{2}}{{\Delta}x^{2}} \right)}}{{\Delta}t}\right] \]

For the frequency \(w\) to be real (ensuring non-growing solutions), the argument of the arccosine function must be within the range [-1, 1]:

\[-1 \le \frac{c^{2} {\Delta}t^{2} \cos{\left(k {\Delta}x \right)} - c^{2} {\Delta}t^{2} + {\Delta}x^{2}}{{\Delta}x^{2}} \le 1 \]

The right-hand side of the inequality is always true. The left-hand side gives the stability condition:

\[-1 \le \frac{c^{2} {\Delta}t^{2} \cos{\left(k {\Delta}x \right)} - c^{2} {\Delta}t^{2} + {\Delta}x^{2}}{{\Delta}x^{2}} \]

\(\cos{\left(k {\Delta}x \right)}\) has a minimum value of -1 (when \(k {\Delta}x = \pi, 3\pi, \dots\)). Substituting this minimum value gives the most restrictive condition:

\[-1 \leq \frac{- 2 c^{2} {\Delta}t^{2} + {\Delta}x^{2}}{{\Delta}x^{2}} \]

\[c {\Delta}t \leq {\Delta}x \]

This condition, often written as \(\frac{c {\Delta}t}{{\Delta}x} \le 1\), ensures numerical stability. It means the numerical domain of dependence must contain the physical domain of dependence; intuitively, the wave shouldn’t travel more than one spatial step (\({\Delta}x\)) in one time step (\({\Delta}t\)).

Fundamental frequency limits:

\[\frac{2 \text{L} f_{0}}{f_{s}} \leq \frac{\text{L}}{\text{M}} \]

\[f_{0} \leq \frac{f_{s}}{2 \text{M}} \]

When \(f_{s} = 44100\) and \(\text{M} = 3\):

\[f_{0} \leq 7350 \]