<jupyter_start><jupyter_text>Modes of the Ball-Channel Pendulum Linear Model<jupyter_code>import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as plt
from resonance.linear_systems import BallChannelPendulumSystem
%matplotlib widget<jupyter_output><empty_output><jupyter_text>A (almost) premade system is available in `resonance`. The only thing missing is the function that calculates the canonical coefficients.<jupyter_code>sys = BallChannelPendulumSystem()
sys.constants
sys.states
def can_coeffs(mp, mb, l, g, r):
    M = np.array([[mp * l**2 + mb * r**2, -mb * r**2],
                  [-mb * r**2, mb * r**2]])
    C = np.zeros((2, 2))
    K = np.array([[g * l * mp, g * mb * r],
                  [g * mb * r, g * mb * r]])
    return M, C, K
sys.canonical_coeffs_func = can_coeffs<jupyter_output><empty_output><jupyter_text>Once the system is completely defined the mass, damping, and stiffness matrices can be calculated and inspected:<jupyter_code>M, C, K = sys.canonical_coefficients()
M
C
K<jupyter_output><empty_output><jupyter_text>Convert to mass normalized form (calculate $\tilde{\mathbf{K}}$)First calculate the Cholesky lower triangular decomposition matrix of $\mathbf{M}$, which is symmetric and postive definite.<jupyter_code>L = la.cholesky(M)
L<jupyter_output><empty_output><jupyter_text>The transpose can be computed with `np.transpose()`, `L.transpose()` or `L.T` for short:<jupyter_code>np.transpose(L)
L.transpose()
L.T<jupyter_output><empty_output><jupyter_text>Check that $\mathbf{L}\mathbf{L}^T$ returns $M$. Note that in Python the `@` operator is used for matrix multiplication. The `*` operator will do elementwise multiplication.<jupyter_code>L @ L.T<jupyter_output><empty_output><jupyter_text>`inv()` computes the inverse, giving $\left(\mathbf{L}^T\right)^{-1}$:<jupyter_code>la.inv(L.T)<jupyter_output><empty_output><jupyter_text>$\mathbf{L}^{-1}\mathbf{M}\left(\mathbf{L}^T\right)^{-1} = \mathbf{I}$. Note that the off-diagonal terms are very small numbers. The reason these are not precisely zero is due to floating point arithmetic and the associated truncation errors.<jupyter_code>la.inv(L) @ M @ la.inv(L.T)<jupyter_output><empty_output><jupyter_text>$\tilde{\mathbf{K}} = \mathbf{L}^{-1}\mathbf{K}\left(\mathbf{L}^T\right)^{-1}$. Note that this matrix is symmetric. It is guaranteed to be symmetric if $\mathbf{K}$ is symmetric.<jupyter_code>Ktilde = la.inv(L) @ K @ la.inv(L.T)
Ktilde<jupyter_output><empty_output><jupyter_text>The entries of $\tilde{\mathbf{K}}$ can be accessed as so:<jupyter_code>k11 = Ktilde[0, 0]
k12 = Ktilde[0, 1]
k21 = Ktilde[1, 0]
k22 = Ktilde[1, 1]<jupyter_output><empty_output><jupyter_text>Calculate the eigenvalues of $\tilde{\mathbf{K}}$The eigenvalues of this 2 x 2 matrix are found by forming the characteristic equation from:$$\textrm{det}\left(\tilde{\mathbf{K}} - \lambda \mathbf{I}\right) = 0$$and solving the resulting quadratic polynomial for its roots, which are the eigenvalues.<jupyter_code>lam1 = (k11 + k22) / 2 + np.sqrt((k11 + k22)**2 - 4 * (k11 * k22 - k12*k21)) / 2
lam1
lam2 = (k11 + k22) / 2 - np.sqrt((k11 + k22)**2 - 4 * (k11 * k22 - k12*k21)) / 2
lam2<jupyter_output><empty_output><jupyter_text>Calculate the eigenfrequencies of the system$\omega_i = \sqrt{\lambda_i}$<jupyter_code>omega1 = np.sqrt(lam1)
omega1
omega2 = np.sqrt(lam2)
omega2<jupyter_output><empty_output><jupyter_text>And in Hertz:<jupyter_code>fn1 = omega1/2/np.pi
fn1
fn2 = omega2/2/np.pi
fn2<jupyter_output><empty_output><jupyter_text>Calculate the eigenvectors of $\tilde{\mathbf{K}}$The eigenvectors can be found by substituting the value for $\lambda$ into:$$\tilde{\mathbf{K}}\hat{q}_0 = \lambda \hat{q}_0$$and solving for $\hat{q}_0$.<jupyter_code>v1 = np.array([-k12 / (k11 - lam1), 1])
v2 = np.array([-k12 / (k11 - lam2), 1])<jupyter_output><empty_output><jupyter_text>Check that they are orthogonal, i.e. the dot product should be zero.<jupyter_code>np.dot(v1, v2)<jupyter_output><empty_output><jupyter_text>The `norm()` function calculates the Euclidean norm, i.e. the vector's magnitude and the vectors can be normalized like so:<jupyter_code>v1_hat = v1 / np.linalg.norm(v1)
v2_hat = v2 / np.linalg.norm(v2)
v1_hat
v2_hat
np.linalg.norm(v1_hat)<jupyter_output><empty_output><jupyter_text>For any size $\tilde{\mathbf{K}}$ the `eig()` function can be used to calculate the eigenvalues and the normalized eigenvectors with one function call:<jupyter_code>evals, evecs = np.linalg.eig(Ktilde)
evals
evecs<jupyter_output><empty_output><jupyter_text>The columns of `evecs` correspond to the entries of `evals`.<jupyter_code>P = evecs
P<jupyter_output><empty_output><jupyter_text>If P contains columns that are orthnormal, then $\mathbf{P}^T \mathbf{P} = \mathbf{I}$. Check this with:<jupyter_code>P.T @ P<jupyter_output><empty_output><jupyter_text>$\mathbf{P}$ can be used to find the matrix $\Lambda$ that decouples the differential equations.<jupyter_code>Lam = P.T @ Ktilde @ P
Lam<jupyter_output><empty_output><jupyter_text>Formulate solution to ODEs (simulation) The trajectory of the coordinates can be found with:$$\bar{c}(t) = \sum_{i=1}^n c_i \sin(\omega_i t + \phi_i) \bar{u}_i$$where$$\phi_i = \arctan \frac{\omega_i \hat{q}_{0i}^T \bar{q}(0)}{\hat{q}_{0i}^T \dot{\bar{q}}(0)}$$and$$c_i = \frac{\hat{q}^T_{0i} \bar{q}(0)}{\sin\phi_i}$$$c_i$ are the modal participation factors and reflect what propotional of each mode is excited given specific initial conditions. If the initial conditions are the eigenmode, $\bar{u}_i$, the all but the $i$th $c_i$ will be zero.A matrix $\mathbf{S} = \left(\mathbf{L}^T\right)^{-1} = \begin{bmatrix}\bar{u}_1 \quad \bar{u}_2\end{bmatrix}$ can be computed such that the columns are $\bar{u}_i$.<jupyter_code>S = la.inv(L.T) @ P
S
u1 = S[:, 0]
u2 = S[:, 1]
u1
u2<jupyter_output><empty_output><jupyter_text>Define the initial coordinates as a scalar factor of the second eigenvector, which sets these values to small angles.<jupyter_code>c0 = S[:, 1] / 400
np.rad2deg(c0)<jupyter_output><empty_output><jupyter_text>Set the initial speeds to zero:<jupyter_code>s0 = np.zeros(2)
s0<jupyter_output><empty_output><jupyter_text>The initial mass normalized coordinates and speeds are then:<jupyter_code>q0 = L.T @ c0
q0
qd0 = L.T @ s0
qd0<jupyter_output><empty_output><jupyter_text>Calculate the modal freqencies in radians per second.<jupyter_code>ws = np.sqrt(evals)
ws<jupyter_output><empty_output><jupyter_text>The phase shifts for each mode can be found. Note that it is important to use `arctan2()` so that the quadrant and thus sign of the arc tangent is properly handled.$$\phi_i = \arctan \frac{\omega_i \hat{q}_{0i}^T \bar{q}(0)}{\hat{q}_{0i}^T \dot{\bar{q}}(0)}$$<jupyter_code>phi1 = np.arctan2(ws * P[:, 0] @ q0, P[:, 0] @ qd0)
phi1
phi2 =  np.arctan2(ws * P[:, 1] @ q0, P[:, 1] @ qd0)
phi2<jupyter_output><empty_output><jupyter_text>All $\phi$'s can be calculated in one line using NumPy's broadcasting feature:<jupyter_code>phis = np.arctan2(ws * P.T @ q0, P.T @ qd0)
phis<jupyter_output><empty_output><jupyter_text>The phase shifts for this particular initial condition are $\pm90$ degrees.<jupyter_code>np.rad2deg(phis)<jupyter_output><empty_output><jupyter_text>Now calculate the modal participation factors.$$c_i = \frac{\hat{q}^T_{0i} \bar{q}(0)}{\sin\phi_i}$$<jupyter_code>cs = P.T @ q0 / np.sin(phis)
cs<jupyter_output><empty_output><jupyter_text>Note that the first participation factor is zero. This is because we've set the initial coordinate to be a scalar function of the second eigenvector.  Simulate<jupyter_code>t = np.linspace(0, 5, num=500)
cs[1] * np.sin(ws[1] * t)<jupyter_output><empty_output><jupyter_text>The following line will give an error because the dimensions of `u1` are not compatible with the dimensions of the preceding portion. It is possible for a single line to work like this if you take advatnage of NumPy's broadcasting rules. See https://scipy-lectures.org/intro/numpy/operations.htmlbroadcasting for more info. The `tile()` function is used to repeat `u1` as many times as needed.<jupyter_code># cs[1] * np.sin(ws[1] * t) * u1
c1 = cs[1] * np.sin(ws[1] * t) * np.tile(u1, (len(t), 1)).T
c1.shape<jupyter_output><empty_output><jupyter_text>`tile()` can be used to create a 2 x 1000 vector that repeats the vector $\hat{u}_i$ allowing a single line to calculate the mode contribution. Now use a loop to calculate the contribution of each mode and build the summation of contributions from each mode:<jupyter_code>ct = np.zeros((2, len(t)))  # 2 x m array to hold coordinates as a function of time
for ci, wi, phii, ui in zip(cs, ws, phis, S.T):
    print(ci, wi, phii, ui)
    ct += ci * np.sin(wi * t + phii) * np.tile(ui, (len(t), 1)).T
def sim(c0, s0, t):
    """Returns the time history of the coordinate vector, c(t) given the initial state and time.
    
    Parameters
    ==========
    c0 : ndarray, shape(n,)
    s0 : ndarray, shape(n,)
    t : ndarray, shape(m,)
    
    Returns
    =======
    c(t) : ndarray, shape(n, m)
    
    """
    q0 = L.T @ c0
    qd0 = L.T @ s0
    ws = np.sqrt(evals)
    phis = np.arctan2(ws * P.T @ q0, P.T @ qd0)
    cs = P.T @ q0 / np.sin(phis)
    c = np.zeros((2, 1000))
    for ci, wi, phii, ui in zip(cs, ws, phis, S.T):
        c += ci * np.sin(wi * t + phii) * np.tile(ui, (len(t), 1)).T
    return c<jupyter_output><empty_output><jupyter_text>Simulate and plot the first mode:<jupyter_code>t = np.linspace(0, 5, num=1000)

c0 = S[:, 0] / np.max(S[:, 0]) * np.deg2rad(10)
s0 = np.zeros(2)

fig, ax = plt.subplots()
ax.plot(t, np.rad2deg(sim(c0, s0, t).T))
ax.set_xlabel('Time [s]')
ax.set_ylabel('Angle [deg]')
ax.legend([r'$\theta$', r'$\phi$'])<jupyter_output><empty_output><jupyter_text>Simulate and plot the second mode:<jupyter_code>t = np.linspace(0, 5, num=1000)

c0 = S[:, 1] / np.max(S[:, 1]) * np.deg2rad(10)
s0 = np.zeros(2)

fig, ax = plt.subplots()
ax.plot(t, np.rad2deg(sim(c0, s0, t).T))
ax.set_xlabel('Time [s]')
ax.set_ylabel('Angle [deg]')
ax.legend([r'$\theta$', r'$\phi$'])<jupyter_output><empty_output><jupyter_text>Compare this to the free response from the system:<jupyter_code>sys.coordinates['theta'] = c0[0]
sys.coordinates['phi'] = c0[1]

sys.speeds['alpha'] = 0
sys.speeds['beta'] = 0
traj = sys.free_response(5.0)
traj[['theta', 'phi']].plot()
sys.animate_configuration(fps=30, repeat=False)<jupyter_output><empty_output><jupyter_text>Simulate with arbitrary initial conditions.<jupyter_code>sys.coordinates['theta'] = np.deg2rad(12.0)
sys.coordinates['phi'] = np.deg2rad(3.0)
traj = sys.free_response(5.0)
traj[['theta', 'phi']].plot()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1546857357922469, 'fraction_numerical': 0.025780955965374483, 'mean_word_length': 5.334326579261025, 'pattern_counts': {'":': 0, '<': 145, '<?xml version=': 0, '>': 145, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '36', 'dir': 'jupyter-structured-clean-dedup', 'id': '4716', 'license': 'MIT', 'n_tokens_mistral': 4208, 'n_tokens_neox': 4016, 'n_words': 1025, 'path': 'notebooks/12-2020/modes_of_a_ball_channel_pendulum.ipynb', 'repo_name': 'gbrault/resonance', 'size': '10628'}