Modes of the Ball-Channel Pendulum Linear Model¶
import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as plt
from resonance.linear_systems import BallChannelPendulumSystem
%matplotlib inline
A (almost) premade system is available in resonance. The only thing missing is the function that calculates the canonical coefficients.
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
Once the system is completely defined the mass, damping, and stiffness matrices can be calculated and inspected:
M, C, K = sys.canonical_coefficients()
M
C
K
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.
L = la.cholesky(M)
L
The transpose can be computed with np.transpose(), L.transpose() or L.T for short:
np.transpose(L)
L.transpose()
L.T
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.
L @ L.T
inv() computes the inverse, giving $\left(\mathbf{L}^T\right)^{-1}$:
la.inv(L.T)
$\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.
la.inv(L) @ M @ la.inv(L.T)
$\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.
Ktilde = la.inv(L) @ K @ la.inv(L.T)
Ktilde
The entries of $\tilde{\mathbf{K}}$ can be accessed as so:
k11 = Ktilde[0, 0]
k12 = Ktilde[0, 1]
k21 = Ktilde[1, 0]
k22 = Ktilde[1, 1]
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.
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
Calculate the eigenfrequencies of the system¶
$\omega_i = \sqrt{\lambda_i}$
omega1 = np.sqrt(lam1)
omega1
omega2 = np.sqrt(lam2)
omega2
And in Hertz:
fn1 = omega1/2/np.pi
fn1
fn2 = omega2/2/np.pi
fn2
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$.
v1 = np.array([-k12 / (k11 - lam1), 1])
v2 = np.array([-k12 / (k11 - lam2), 1])
Check that they are orthogonal, i.e. the dot product should be zero.
np.dot(v1, v2)
The norm() function calculates the Euclidean norm, i.e. the vector's magnitude and the vectors can be normalized like so:
v1_hat = v1 / np.linalg.norm(v1)
v2_hat = v2 / np.linalg.norm(v2)
v1_hat
v2_hat
np.linalg.norm(v1_hat)
For any size $\tilde{\mathbf{K}}$ the eig() function can be used to calculate the eigenvalues and the normalized eigenvectors with one function call:
evals, evecs = np.linalg.eig(Ktilde)
evals
evecs
The columns of evecs correspond to the entries of evals.
P = evecs
P
If P contains columns that are orthnormal, then $\mathbf{P}^T \mathbf{P} = \mathbf{I}$. Check this with:
P.T @ P
$\mathbf{P}$ can be used to find the matrix $\Lambda$ that decouples the differential equations.
Lam = P.T @ Ktilde @ P
Lam
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$.
S = la.inv(L.T) @ P
S
u1 = S[:, 0]
u2 = S[:, 1]
u1
u2
Define the initial coordinates as a scalar factor of the second eigenvector, which sets these values to small angles.
c0 = S[:, 1] / 400
np.rad2deg(c0)
Set the initial speeds to zero:
s0 = np.zeros(2)
s0
The initial mass normalized coordinates and speeds are then:
q0 = L.T @ c0
q0
qd0 = L.T @ s0
qd0
Calculate the modal freqencies in radians per second.
ws = np.sqrt(evals)
ws
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.
phi1 = np.arctan2(ws * P[:, 0] @ q0, P[:, 0] @ qd0)
phi1
phi2 = np.arctan2(ws * P[:, 1] @ q0, P[:, 1] @ qd0)
phi2
All $\phi$'s can be calculated in one line using NumPy's broadcasting feature:
phis = np.arctan2(ws * P.T @ q0, P.T @ qd0)
phis
The phase shifts for this particular initial condition are $\pm90$ degrees.
np.rad2deg(phis)
Now calculate the modal participation factors.
$$ c_i = \frac{\hat{q}^T_{0i} \bar{q}(0)}{\sin\phi_i} $$cs = P.T @ q0 / np.sin(phis)
cs
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¶
t = np.linspace(0, 5, num=500)
cs[1] * np.sin(ws[1] * t)
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.html#broadcasting for more info. The tile() function is used to repeat u1 as many times as needed.
# 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
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:
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
Simulate and plot the first mode:
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$'])
Simulate and plot the second mode:
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$'])
Compare this to the free response from the system:
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)
Simulate with arbitrary initial conditions.
sys.coordinates['theta'] = np.deg2rad(12.0)
sys.coordinates['phi'] = np.deg2rad(3.0)
traj = sys.free_response(5.0)
traj[['theta', 'phi']].plot()
