In [1]:
import numpy as np

from resonance.linear_systems import BallChannelPendulumSystem
In [2]:
import matplotlib.pyplot as plt
In [3]:
%matplotlib inline
In [4]:
sys = BallChannelPendulumSystem()
In [5]:
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
In [6]:
sys.canonical_coeffs_func = can_coeffs
In [7]:
M, C, K = sys.canonical_coefficients()
In [8]:
M
Out[8]:
array([[ 5.15e-04, -3.50e-05],
       [-3.50e-05,  3.50e-05]])
In [9]:
C
Out[9]:
array([[0., 0.],
       [0., 0.]])
In [10]:
K
Out[10]:
array([[0.023544 , 0.0034335],
       [0.0034335, 0.0034335]])
In [11]:
L = np.linalg.cholesky(M)
L
Out[11]:
array([[ 0.02269361,  0.        ],
       [-0.00154228,  0.00571151]])
In [12]:
L.T
Out[12]:
array([[ 0.02269361, -0.00154228],
       [ 0.        ,  0.00571151]])
In [13]:
L @ L.T
Out[13]:
array([[ 5.15e-04, -3.50e-05],
       [-3.50e-05,  3.50e-05]])
In [14]:
np.linalg.inv(L.T)
Out[14]:
array([[ 44.06526492,  11.89898149],
       [  0.        , 175.08501336]])
In [15]:
np.linalg.inv(L) @ M @ np.linalg.inv(L.T)
Out[15]:
array([[ 1.00000000e+00, -1.82171447e-17],
       [-1.16059609e-17,  1.00000000e+00]])
In [16]:
Ktilde = np.linalg.inv(L) @ K @ np.linalg.inv(L.T)
Ktilde
Out[16]:
array([[ 45.71650485,  38.83489484],
       [ 38.83489484, 122.89287015]])
In [17]:
k11 = Ktilde[0, 0]
k12 = Ktilde[0, 1]
k21 = Ktilde[1, 0]
k22 = Ktilde[1, 1]

Eigenvalues

In [18]:
lam1 = (k11 + k22) / 2 + np.sqrt((k11 + k22)**2 - 4 * (k11 * k22 - k12*k21)) / 2
lam1
Out[18]:
139.05134855775452
In [19]:
lam2 = (k11 + k22) / 2 - np.sqrt((k11 + k22)**2 - 4 * (k11 * k22 - k12*k21)) / 2
lam2
Out[19]:
29.558026442245463

Eigenfrequencies

In [20]:
omega1 = np.sqrt(lam1)
omega1
Out[20]:
11.792003585385926
In [21]:
omega2 = np.sqrt(lam2)
omega2
Out[21]:
5.436729388358912

Eigenvectors

In [22]:
v1 = np.array([-k12 / (k11 - lam1), 1])
v1 / np.linalg.norm(v1)
Out[22]:
array([0.38415493, 0.92326864])
In [23]:
v2 = np.array([-k12 / (k11 - lam2), 1])
v2 / np.linalg.norm(v2)
Out[23]:
array([-0.92326864,  0.38415493])

Using NumPy

In [24]:
evals, evecs = np.linalg.eig(Ktilde)
In [25]:
evals
Out[25]:
array([ 29.55802644, 139.05134856])
In [26]:
evecs
Out[26]:
array([[-0.92326864, -0.38415493],
       [ 0.38415493, -0.92326864]])

Modal coordinates

In [27]:
P = evecs
In [28]:
P.T @ P
Out[28]:
array([[1., 0.],
       [0., 1.]])
In [29]:
Lam = P.T @ Ktilde @ P
Lam
Out[29]:
array([[ 2.95580264e+01, -5.32907052e-15],
       [ 0.00000000e+00,  1.39051349e+02]])

The eigenvectors can be put back into units of x with:

In [30]:
S = np.linalg.inv(L.T) @ P
S
Out[30]:
array([[ -36.11302491,  -27.9138454 ],
       [  67.25977163, -161.6505027 ]])

The trajectory of building's coordinates can be found with:

$$ \mathbf{x}(t) = \sum_{i=1}^n c_i \sin(\omega_i t + \phi_i) \mathbf{u}_i $$

where

$$ \phi_i = \arctan \frac{\omega_i \mathbf{v}_i^T \mathbf{q}_0}{\mathbf{v}_i^T \dot{\mathbf{q}}_0} $$

and

$$ c_i = \frac{\mathbf{v}^T_i \mathbf{q}_0}{\sin\phi_i} $$

$d_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, $\mathbf{u}_i$, the all but the $i$th $d_i$ will be zero.

In [31]:
x0 = S[:, 1] / 300
np.rad2deg(x0)
Out[31]:
array([ -5.33115177, -30.87297187])
In [32]:
x0[0]/x0[1]
Out[32]:
0.17268022639045402
In [33]:
xd0 = np.zeros(2)
xd0
Out[33]:
array([0., 0.])
In [34]:
q0 = L.T @ x0
q0
Out[34]:
array([-0.00128052, -0.00307756])
In [35]:
q0[0]/q0[1]
Out[35]:
0.4160814257260209
In [36]:
ws = np.sqrt(evals)
ws
Out[36]:
array([ 5.43672939, 11.79200359])
In [37]:
phis = np.arctan2(ws * P.T @ q0, P.T @ xd0)
phis
Out[37]:
array([-1.57079633,  1.57079633])
In [38]:
cs = P.T @ q0 / np.sin(phis)
cs
Out[38]:
array([2.16840434e-19, 3.33333333e-03])
In [39]:
t = np.linspace(0, 5, num=1000)
In [40]:
x = np.zeros((2, 1000))
for ci, wi, phii, ui in zip(cs, ws, phis, S.T):
    x += ci * np.sin(wi * t + phii) * np.tile(ui, (len(t), 1)).T
In [41]:
def sim(x0, xd0, t):
    q0 = L.T @ x0
    ws = np.sqrt(evals)
    phis = np.arctan2(ws * P.T @ q0, P.T @ xd0)
    cs = P.T @ q0 / np.sin(phis)
    x = np.zeros((2, 1000))
    for ci, wi, phii, ui in zip(cs, ws, phis, S.T):
        x += ci * np.sin(wi * t + phii) * np.tile(ui, (len(t), 1)).T
    return x
In [42]:
t = np.linspace(0, 5, num=1000)
x0 = S[:, 0] / np.max(S[:, 0]) * np.deg2rad(10)
xd0 = np.zeros(2)
plt.figure()
plt.plot(t, sim(x0, xd0, t).T)
Out[42]:
[<matplotlib.lines.Line2D at 0x7f418317f190>,
 <matplotlib.lines.Line2D at 0x7f418317f2b0>]
In [43]:
t = np.linspace(0, 5, num=1000)
x0 = S[:, 1] / np.max(S[:, 1]) * np.deg2rad(10)
xd0 = np.zeros(2)
plt.figure()
plt.plot(t, sim(x0, xd0, t).T)
Out[43]:
[<matplotlib.lines.Line2D at 0x7f4182942e50>,
 <matplotlib.lines.Line2D at 0x7f41829500a0>]
In [44]:
sys.coordinates['theta'] = x0[0]
sys.coordinates['phi'] = x0[1]

sys.speeds['alpha'] = 0
sys.speeds['beta'] = 0
In [45]:
traj = sys.free_response(2 * np.pi / ws[1] * 10)
In [46]:
traj[['theta', 'phi']].plot()
Out[46]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f4182859cd0>
In [47]:
sys.animate_configuration()
Out[47]:
<matplotlib.animation.FuncAnimation at 0x7f41827bc670>