Time integrators
A collection of efficient time integration methods for the damped oscillator with optional nonlinearities and initial conditions or excitation inputs.
Stormer-verlet
The stormer-verlet is a centered difference scheme used to approximate derivatives. In the present case we use it for approximating the second and first derivatives of an oscillator defined by the following differential equation:
\[ \ddot{q} + c \dot{q} + k q = f(t) \]
The finite difference operators are defined as:
\[ \begin{align} \delta_t q &= \frac{q^{n+1} - q^{n-1}}{2 h} \\ \delta_{tt} q &= \frac{q^{n+1} - 2 q^n + q^{n-1}}{h^2} \end{align} \]
for the first derivative and the second derivative respectively. The difference equation for the oscillator is then:
\[ \delta_{tt} q + c \delta_t q + k q = f(t) \]
after expanding and some algebraic manipulation to isolate \(q^{n+1}\) we get:
\[ \bigl(\tfrac{1}{h^2} + \tfrac{c}{2h}\bigr)\,q^{n+1} + \bigl(-\tfrac{2}{h^2} + k\bigr) q^n + \bigl(\tfrac{1}{h^2} - \tfrac{c}{2h}\bigr)\,q^{n-1} = f(t^n). \]
\[ \boxed{ q^{n+1} = \frac{2 h^2}{2 + c h} \Bigl[ f(t^n) + \Bigl(\tfrac{2}{h^2} - k\Bigr) q^n + \Bigl(-\tfrac{1}{h^2} + \tfrac{c}{2h}\Bigr) q^{n-1} \Bigr]. } \]
To make it more readable we can define the following constants:
\[ \begin{align} a &= \frac{2 h^2}{2 + c h} \\ b &= \frac{2}{h^2} - k \\ c &= -\frac{1}{h^2} + \frac{c}{2h} \end{align} \]
make_tm_nl_fn
make_tm_nl_fn (lambda_mu, factors)
Returns a function that computes nl given q.
make_vk_nl_fn
make_vk_nl_fn (H)
Returns a function that computes nl given q.
make_identity_nl_fn
make_identity_nl_fn ()
string_tau_with_density
string_tau_with_density (string_params)
plate_tau_with_density
plate_tau_with_density (plate_params)
second_order_step
second_order_step (u0:jax.Array, v0:jax.Array, dt:float, gamma2_mu:jax.Array, omega_mu_squared:jax.Array)
One step of 2nd-order Taylor expansion for damped oscillator. Returns: (u1, v1)
| Type | Details | |
|---|---|---|
| u0 | Array | initial conditions (n_modes,) |
| v0 | Array | initial conditions (n_modes,) |
| dt | float | time step |
| gamma2_mu | Array | damping (n_modes,) |
| omega_mu_squared | Array | frequency (n_modes,) |
rk4_step
rk4_step (u0:jax.Array, v0:jax.Array, dt:float, gamma2_mu:jax.Array, omega_mu_squared:jax.Array)
One step of RK4 for the second-order damped oscillator. Returns: (u1, v1)
| Type | Details | |
|---|---|---|
| u0 | Array | initial conditions (n_modes,) |
| v0 | Array | initial conditions (n_modes,) |
| dt | float | time step |
| gamma2_mu | Array | damping (n_modes,) |
| omega_mu_squared | Array | frequency (n_modes,) |
solve_sv_initial_conditions
solve_sv_initial_conditions (gamma2_mu, omega_mu_squared, u0:jax.Array, v0:jax.Array, dt:float, n_steps:int, nl_fn:collections.abc.Callable)
| Type | Details | |
|---|---|---|
| gamma2_mu | (n_modes,) | |
| omega_mu_squared | (n_modes,) | |
| u0 | Array | initial conditions (n_modes,) |
| v0 | Array | initial conditions (n_modes,) |
| dt | float | |
| n_steps | int | |
| nl_fn | Callable |
solve_sv_excitation
solve_sv_excitation (gamma2_mu, omega_mu_squared, modal_excitation:jax.Array, dt:float, nl_fn:collections.abc.Callable)
| Type | Details | |
|---|---|---|
| gamma2_mu | (n_modes,) | |
| omega_mu_squared | (n_modes,) | |
| modal_excitation | Array | (T, n_modes) |
| dt | float | |
| nl_fn | Callable |
solve_sv_vk_jax_scan
solve_sv_vk_jax_scan (A_inv:jax.Array, B:jax.Array, C:jax.Array, modal_excitation:jax.Array, nl_fn:collections.abc.Callable)
| Type | Details | |
|---|---|---|
| A_inv | Array | |
| B | Array | |
| C | Array | |
| modal_excitation | Array | (T, n_modes) |
| nl_fn | Callable |
Discretised transfer function method
solve_tf_initial_conditions
solve_tf_initial_conditions (gamma2_mu, omega_mu_squared, u0:jax.Array, v0:jax.Array, dt:float, n_steps:int, nl_fn:collections.abc.Callable)
Solve using transfer-function (TF) based recurrence.
| Type | Details | |
|---|---|---|
| gamma2_mu | ||
| omega_mu_squared | ||
| u0 | Array | initial conditions (n_modes,) |
| v0 | Array | initial conditions (n_modes,) |
| dt | float | |
| n_steps | int | |
| nl_fn | Callable |
solve_tf_excitation
solve_tf_excitation (gamma2_mu, omega_mu_squared, modal_excitation:jax.Array, dt:float, nl_fn:collections.abc.Callable)
Solve using transfer-function (TF) based recurrence.
| Type | Details | |
|---|---|---|
| gamma2_mu | ||
| omega_mu_squared | ||
| modal_excitation | Array | (T, n_modes) |
| dt | float | |
| nl_fn | Callable |
Sinusoidal solve
Solve the system of ODEs using complex exponentials
The system of ODEs is given by:
\[ \ddot{q} + 2 \gamma \dot{q} + \omega^2 q = 0 \]
where \(\gamma\) is the damping coefficient and \(\omega\) is the frequency. The damped frequencies are given by:
\[ \tilde{\omega} = \omega \sqrt{1 - \gamma^2} = \sqrt{\omega^2 - \gamma^2} \] The eigenvalues after the dispersion relation are then \[ s_\pm = -\gamma \pm i \tilde{\omega} \]
and discrete time eigenvalues (poles) are \[ z_\pm = e^{s_\pm \Delta t} = e^{-\gamma \Delta t} e^{\pm i \tilde{\omega} \Delta t} \]
solve_sinusoidal
solve_sinusoidal (gamma2_mu, omega_mu_squared, ic, n_steps, dt)
Solve the system of ODEs using complex exponentials NB: this assumes the ic is only for positions and that the initial velocities are 0
| Type | Details | |
|---|---|---|
| gamma2_mu | jnp.ndarray | Damping coefficients |
| omega_mu_squared | jnp.ndarray | Squared frequencies |
| ic | jnp.ndarray | Initial conditions |
| n_steps | int | Number of steps |
| dt | float | Time step |
| Returns | jnp.ndarray | Modal solution |
solve_sinusoidal_excitation
solve_sinusoidal_excitation (gamma2_mu, omega_mu_squared, modal_excitation:jax.Array, dt:float)
Solve the modal system with sinusoidal response for external excitation using parallel scan.
| Type | Details | |
|---|---|---|
| gamma2_mu | jnp.ndarray | Damping coefficients (n_modes,) |
| omega_mu_squared | jnp.ndarray | Squared frequencies (n_modes,) |
| modal_excitation | Array | (T, n_modes) |
| dt | float | Time step |
| Returns | jnp.ndarray | Modal solution (T, n_modes) |
solve_tf_ic
solve_tf_ic (gamma2_mu, omega_mu_squared, ic, n_steps, dt, nl_fn)