from matplotlib import pyplot as plt
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}} \]
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
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 |
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, )
= 0.65
length = 0.5188e-6
A = 0.141e-12
I = 1140
rho = 5.4e9
E = 8e-5
d1 = 1.4e-5
d3 = 60.97
Ts0 = 101
n_gridpoints = 99
n_max_modes = 96000
sr
= Wave1dSolverTensionModulated(
solver =sr,
sampling_rate=0.1,
final_time=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=True,
use_nonlinear# d3=1e-7,
)
# initial conditions
= 0.02
ic_max_amplitude = generate_initial_condition(
u0, v0 =np.random.default_rng(354),
rng=Gaussian(num_points=n_gridpoints),
generator=ic_max_amplitude,
ic_max_amplitude=False,
ic_amplitude_random
)
= solver.solve(u0, v0)
t, u, v 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,)
= 7000
n_samples_plot = plt.subplots(1, 2, figsize=(5, 6))
fig, 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")
ax[
fig.tight_layout()
plt.show()# plot the solution at the endpoints
=(10, 5))
plt.figure(figsize0], label="left")
plt.plot(t, u[:, -1], label="right", alpha=0.5)
plt.plot(t, u[:,
plt.legend()
plt.show()
# Plot FFT of the solution for a point on the string, in db scale
=(10, 5))
plt.figure(figsize= np.fft.rfftfreq(u.shape[-1], d=solver.dt)
freq 20 * np.log10(np.abs(np.fft.rfft(u[17, :]))), label="fft")
plt.plot(freq, # 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)
100, 5000)
plt.xlim(
plt.legend()
plt.show()
# Plot spectrogram of the solution for a point on the string
=(10, 5))
plt.figure(figsize17], Fs=solver.sampling_rate, NFFT=2048, noverlap=1024)
plt.specgram(u[:, 0, 20000)
plt.ylim(
plt.show()
# Plot FFT of the initial condition
=(10, 5))
plt.figure(figsize= np.fft.rfftfreq(u.shape[-1], d=solver.dx)
freq # plt.plot(freq, np.abs(np.fft.rfft(u0)), label="fft")
abs(np.fft.rfft(u0)), label="fft")
plt.plot(np.# plt.xlim(0, 1)
plt.legend()
# plot the initial conditions
=(10, 5))
plt.figure(figsize="u0")
plt.plot(solver.grid, u0, label0, :], label="u0_sol", alpha=0.5)
plt.plot(solver.grid, u[
plt.legend()
plt.show()
# Plot the difference between the initial condition given and the one obtained from the solution
=(10, 5))
plt.figure(figsize- u[0], label="u0 - u0_sol")
plt.plot(solver.grid, u0
plt.legend() plt.show()
from IPython import display as ipd
# Play the solution as audio
27], rate=solver.sampling_rate) ipd.Audio(u[:,