Tension modulated string

Step 1:

Start with the non-linear tension modulated string model. The model is defined by the following equations:

\[ \begin{align} \rho A \frac{\partial^2 y}{\partial t^2} + EI \frac{\partial^4 y}{\partial x^4} - T_{NL}(y) \frac{\partial^2 y}{\partial x^2} + d_{1} \frac{\partial y}{\partial t} + d_{3} \frac{\partial^3 y}{\partial t^3} = f(x,t) \end{align} \]

where \(T_{NL}\) is the non-linear tension function truncated at the third order:

\[ \begin{aligned} T_{N L}(y) & =T_0+T_1(y(x, t)) \\ & =T_0+\frac{E A}{2 l_0} \int_0^{l_0} y^2(x, t) d x . \end{aligned} \]

Step 2:

Perform a transformation of the PDE with respect to the spatial variable \(x\). We can do this using the transformation kernel \(K(x, \beta_{\mu})\).

\[ \begin{aligned} K\left(x, \beta_\mu\right) & =\sin \left(\frac{\mu \pi}{l_0} x\right) \\ \beta_\mu^4 & =E I\left(\frac{\mu \pi}{l_0}\right)^4+T_0\left(\frac{\mu \pi}{l_0}\right)^2 \end{aligned} \]

where \(\beta_\mu\) is the frequency variable. We see that the \(K(x, \beta_{\mu})\) is a function that generates the fourier sine series of the spatial variable \(x\).

The spatial tranformation of the non-linear term \(T_{1}\) is given as:

\[ \mathcal{T}\left\{T_1\right\} = \bar{b}(\mu,y,\bar{y}) \]

We now can write the transformed PDE as:

\[ \boxed{ \begin{aligned} \rho A \ddot{\bar{y}}(\mu, t) + \left(d_3 \eta_\mu^2+d_1\right) \dot{\bar{y}}(\mu, t) + \beta_\mu^4 \bar{y}(\mu, t)-\bar{b}(\mu, y, \bar{y}) = 0 \end{aligned} } \]

Where \(\eta_\mu=-\left(\mu \pi / l_0\right)^2\). The initial conditions are given as:

\[ \begin{aligned} \bar{y}_0(\mu) &= \mathcal{T}{\left\{y_0(x)\right\}} \\ \dot{\bar{y}}_0(\mu) &= \mathcal{T}{\left\{\dot{y}_0(x)\right\}} = 0 \end{aligned} \]

Step 3:

Discretize the non-linear term \(\bar{b}(\mu, y, \bar{y})\)

\[ \begin{aligned} \bar{b}^d(\mu) & =\eta_\mu^2 \frac{E A \bar{y}(\mu)}{2 l_0} \int_0^{l_0}\left[\sum_{\nu=1}^M \frac{1}{N_\nu} \bar{y}^d(\nu) K(\xi, \nu)\right]^{\prime 2} d x \\ & =\eta_\mu^2 \frac{E A \pi^2 \bar{y}(\mu)}{l_0^4} \sum_{\nu=1}^M \nu^2\left(\bar{y}^d(\nu)\right)^2 \end{aligned} \]

Step 4

Write as a system of first order ODEs:

\[ \begin{aligned} \dot{\bar{y}}(\mu, t) &= \bar{v}(\mu, t) \\ \dot{\bar{v}}(\mu, t) &= \frac{-\left(d_3 \eta_\mu^2+d_1\right)}{\rho A} \bar{v}(\mu, t) - \frac{\beta_\mu^4}{\rho A} \bar{y}(\mu, t) + \frac{\bar{b}(\mu, y, \bar{y})}{\rho A} \end{aligned} \]

This is a system of first-order ODEs. In matrix form, we can write this as:

\[ \boxed{ \begin{aligned} \mathbf{\dot{\bar{y}}} &= \mathbf{\bar{v}} \\ \mathbf{\dot{\bar{v}}} &= -\mathbf{M_v} \mathbf{\bar{v}} - \mathbf{M_y} \mathbf{\bar{y}} + \mathbf{\bar{b}} \end{aligned} } \]

The Matrices \(\mathbf{M_v}\), \(\mathbf{M_y}\) are diagonal matrices whose order is equal to the number of modes \(M\).

\[ \begin{aligned} \mathbf{M_v} &= \text{diag}\left(\frac{d_3 \eta_\mu^2+d_1}{\rho A}\right) \\ \mathbf{M_y} &= \text{diag}\left(\frac{\beta_\mu^4}{\rho A}\right) \end{aligned} \]

The vector \(\mathbf{b}\) contains the non-linear terms and is given as:

\[ \mathbf{\Eta} = \text{diag}\left(\eta_\mu^2\right) \\ T_{1} = \frac{E A \pi^2}{l_0^4} \mathbf{\bar{y}}(\mu)^T \cdot \Eta \cdot \mathbf{\bar{y}} \\ \mathbf{\bar{b}} = T_{1} \text{diag}(\frac{\eta_\mu^2}{\rho A}) \mathbf{\bar{y}} \]

from matplotlib import pyplot as plt

The parameters of the model are:

• \(N\): number of points in the string
• \(T\): total time of simulation
• \(f_s\): sampling rate
• \(l\): length of string
• \(A\): string cross-sectional area
• \(E\): Young’s modulus
• \(I\): string moment of inertia
• \(\rho\): density
• \(Ts_0\): initial tension
• \(d_{1}\): string frequency independent damping
• \(d_{3}\): string frequency dependent damping
• \(M\): number of modes to include in the modal expansion

---

source

Wave1dSolverTensionModulated

 Wave1dSolverTensionModulated (sampling_rate:int=48000,
                               final_time:float=0.5, n_gridpoints:int=101,
                               length:float=1.0, A:float=1.9634e-07,
                               I:float=2.454e-14, rho:float=7800,
                               E:int=190000000000.0, d1:float=0.004,
                               d3:float=6e-05, Ts0:float=150,
                               n_max_modes:int=50,
                               use_nonlinear:bool=True,
                               method:str='DOP853', rtol:float=1e-12,
                               atol:float=1e-14)

*Tension modulated wave equation solver.

This solver is based on the paper “Sound Synthesis With Tension Modulated Nonlinearities” and the accompanying code to the paper “Physical Modeling using Recurrent Neural Networks with Fast Convolutional Layers”.*

	Type	Default	Details
sampling_rate	int	48000	1/s Temporal sampling frequency
final_time	float	0.5	s Duration of the simulation
n_gridpoints	int	101	pts/m Spatial sampling grid
length	float	1.0	m Length at rest
A	float	1.9634e-07	m**2 Cross section area
I	float	2.454e-14	m**4 Moment of intertia
rho	float	7800	kg/m**3 Density
E	int	190000000000.0	Pa Young’s modulus
d1	float	0.004	kg/(ms) Frequency independent loss
d3	float	6e-05	kg m/s Frequency dependent loss
Ts0	float	150	N Tension
n_max_modes	int	50	Number of modal expansion terms
use_nonlinear	bool	True	Use nonlinear wave equation
method	str	DOP853	Integration method
rtol	float	1e-12	Relative tolerance
atol	float	1e-14	Absolute tolerance

---

source

nonlinear_string_func

 nonlinear_string_func (t, state, n_max_modes, use_nonlinear, pre, H, H_1,
                        M_v, M_y)

	Type	Details
t		time
state		state vector
n_max_modes		number of modes
use_nonlinear		use nonlinear wave equation
pre		precalculated constant
H		H matrix
H_1		H_1 matrix
M_v		M_v matrix
M_y		M_y matrix
Returns	ndarray	state at timestep t and position x (u(x, t))

Test the solver

from physmodjax.solver.generator import (
    Gaussian,
    generate_initial_condition,
)

length = 0.65
A = 0.5188e-6
I = 0.141e-12
rho = 1140
E = 5.4e9 
d1 = 8e-5
d3 = 1.4e-5
Ts0 = 60.97
n_gridpoints = 101
n_max_modes = 99
sr = 96000


solver = Wave1dSolverTensionModulated(
    sampling_rate=sr,
    final_time=0.1,
    A=A,
    I=I,
    rho=rho,
    E=E,
    d1=d1,
    d3=d3,
    Ts0=Ts0,
    length=length,
    n_gridpoints=n_gridpoints,
    n_max_modes=n_max_modes,
    use_nonlinear=True,
    # d3=1e-7,
)

# initial conditions
ic_max_amplitude = 0.02
u0, v0 = generate_initial_condition(
    rng=np.random.default_rng(354),
    generator=Gaussian(num_points=n_gridpoints),
    ic_max_amplitude=ic_max_amplitude,
    ic_amplitude_random=False,
)

t, u, v = solver.solve(u0, v0)
print(f"u shape: {u.shape}")

dx: 0.006500000000000001 in meters
dt: 1.0416666666666666e-05 in seconds
number of points (n_gridpoints): (101,)
time in samples (nt): (9600,)

n_samples_plot = 7000
fig, ax = plt.subplots(1, 2, figsize=(5, 6))
ax[0].pcolormesh(solver.grid, t[-n_samples_plot:], u[-n_samples_plot:, :])
ax[0].set_title("u")
ax[0].set_xlabel("x")
ax[0].set_ylabel("t")
ax[1].pcolormesh(solver.grid, t[-n_samples_plot:], v[-n_samples_plot:, :])
ax[1].set_title("v")
ax[1].set_xlabel("x")
ax[1].set_ylabel("t")
fig.tight_layout()
plt.show()
# plot the solution at the endpoints
plt.figure(figsize=(10, 5))
plt.plot(t, u[:, 0], label="left")
plt.plot(t, u[:, -1], label="right", alpha=0.5)
plt.legend()
plt.show()

# Plot FFT of the solution for a point on the string, in db scale
plt.figure(figsize=(10, 5))
freq = np.fft.rfftfreq(u.shape[-1], d=solver.dt)
plt.plot(freq, 20 * np.log10(np.abs(np.fft.rfft(u[17, :]))), label="fft")
# plt.figure(figsize=(10, 5))
# freq = np.fft.rfftfreq(u.shape[0], d=solver.dt)
# plt.plot(freq, np.abs(np.fft.rfft(u[:, 17])), label="fft")
# plt.xlim(3000, 8000)
plt.xlim(100, 5000)
plt.legend()
plt.show()

# Plot spectrogram of the solution for a point on the string
plt.figure(figsize=(10, 5))
plt.specgram(u[:, 17], Fs=solver.sampling_rate, NFFT=2048, noverlap=1024)
plt.ylim(0, 20000)
plt.show()

# Plot FFT of the initial condition
plt.figure(figsize=(10, 5))
freq = np.fft.rfftfreq(u.shape[-1], d=solver.dx)
# plt.plot(freq, np.abs(np.fft.rfft(u0)), label="fft")
plt.plot(np.abs(np.fft.rfft(u0)), label="fft")
# plt.xlim(0, 1)
plt.legend()

# plot the initial conditions
plt.figure(figsize=(10, 5))
plt.plot(solver.grid, u0, label="u0")
plt.plot(solver.grid, u[0, :], label="u0_sol", alpha=0.5)
plt.legend()
plt.show()

# Plot the difference between the initial condition given and the one obtained from the solution
plt.figure(figsize=(10, 5))
plt.plot(solver.grid, u0 - u[0], label="u0 - u0_sol")
plt.legend()
plt.show()

from IPython import display as ipd

# Play the solution as audio
ipd.Audio(u[:, 27], rate=solver.sampling_rate)