# Modes of the Ball-Channel Pendulum Linear Model¶

In [1]:
import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as plt
from resonance.linear_systems import BallChannelPendulumSystem

In [2]:
%matplotlib inline


A (almost) premade system is available in resonance. The only thing missing is the function that calculates the canonical coefficients.

In [3]:
sys = BallChannelPendulumSystem()

In [4]:
sys.constants

Out[4]:
{'mp': 0.012, 'mb': 0.0035, 'r': 0.1, 'l': 0.2, 'g': 9.81}
In [5]:
sys.states

Out[5]:
{'theta': 0.17453292519943295, 'phi': -0.17453292519943295, 'alpha': 0.0, 'beta': 0.0}
In [6]:
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 [7]:
sys.canonical_coeffs_func = can_coeffs


Once the system is completely defined the mass, damping, and stiffness matrices can be calculated and inspected:

In [8]:
M, C, K = sys.canonical_coefficients()

In [9]:
M

Out[9]:
array([[ 5.15e-04, -3.50e-05],
[-3.50e-05,  3.50e-05]])
In [10]:
C

Out[10]:
array([[0., 0.],
[0., 0.]])
In [11]:
K

Out[11]:
array([[0.023544 , 0.0034335],
[0.0034335, 0.0034335]])

## 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.

In [12]:
L = la.cholesky(M)
L

Out[12]:
array([[ 0.02269361,  0.        ],
[-0.00154228,  0.00571151]])

The transpose can be computed with np.transpose(), L.transpose() or L.T for short:

In [13]:
np.transpose(L)

Out[13]:
array([[ 0.02269361, -0.00154228],
[ 0.        ,  0.00571151]])
In [14]:
L.transpose()

Out[14]:
array([[ 0.02269361, -0.00154228],
[ 0.        ,  0.00571151]])
In [15]:
L.T

Out[15]:
array([[ 0.02269361, -0.00154228],
[ 0.        ,  0.00571151]])

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.

In [16]:
L @ L.T

Out[16]:
array([[ 5.15e-04, -3.50e-05],
[-3.50e-05,  3.50e-05]])

inv() computes the inverse, giving $\left(\mathbf{L}^T\right)^{-1}$:

In [17]:
la.inv(L.T)

Out[17]:
array([[ 44.06526492,  11.89898149],
[  0.        , 175.08501336]])

$\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.

In [18]:
la.inv(L) @ M @ la.inv(L.T)

Out[18]:
array([[ 1.00000000e+00, -1.82171447e-17],
[-1.16059609e-17,  1.00000000e+00]])

$\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.

In [19]:
Ktilde = la.inv(L) @ K @ la.inv(L.T)
Ktilde

Out[19]:
array([[ 45.71650485,  38.83489484],
[ 38.83489484, 122.89287015]])

The entries of $\tilde{\mathbf{K}}$ can be accessed as so:

In [20]:
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.

In [21]:
lam1 = (k11 + k22) / 2 + np.sqrt((k11 + k22)**2 - 4 * (k11 * k22 - k12*k21)) / 2
lam1

Out[21]:
139.05134855775452
In [22]:
lam2 = (k11 + k22) / 2 - np.sqrt((k11 + k22)**2 - 4 * (k11 * k22 - k12*k21)) / 2
lam2

Out[22]:
29.558026442245463

# Calculate the eigenfrequencies of the system¶

$\omega_i = \sqrt{\lambda_i}$

In [23]:
omega1 = np.sqrt(lam1)
omega1

Out[23]:
11.792003585385926
In [24]:
omega2 = np.sqrt(lam2)
omega2

Out[24]:
5.436729388358912

And in Hertz:

In [25]:
fn1 = omega1/2/np.pi
fn1

Out[25]:
1.876755659571523
In [26]:
fn2 = omega2/2/np.pi
fn2

Out[26]:
0.8652823564102976

# 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$.

In [27]:
v1 = np.array([-k12 / (k11 - lam1), 1])

In [28]:
v2 = np.array([-k12 / (k11 - lam2), 1])


Check that they are orthogonal, i.e. the dot product should be zero.

In [29]:
np.dot(v1, v2)

Out[29]:
2.220446049250313e-16

The norm() function calculates the Euclidean norm, i.e. the vector's magnitude and the vectors can be normalized like so:

In [30]:
v1_hat = v1 / np.linalg.norm(v1)
v2_hat = v2 / np.linalg.norm(v2)

In [31]:
v1_hat

Out[31]:
array([0.38415493, 0.92326864])
In [32]:
v2_hat

Out[32]:
array([-0.92326864,  0.38415493])
In [33]:
np.linalg.norm(v1_hat)

Out[33]:
1.0

For any size $\tilde{\mathbf{K}}$ the eig() function can be used to calculate the eigenvalues and the normalized eigenvectors with one function call:

In [34]:
evals, evecs = np.linalg.eig(Ktilde)

In [35]:
evals

Out[35]:
array([ 29.55802644, 139.05134856])
In [36]:
evecs

Out[36]:
array([[-0.92326864, -0.38415493],
[ 0.38415493, -0.92326864]])

The columns of evecs correspond to the entries of evals.

In [37]:
P = evecs
P

Out[37]:
array([[-0.92326864, -0.38415493],
[ 0.38415493, -0.92326864]])

If P contains columns that are orthnormal, then $\mathbf{P}^T \mathbf{P} = \mathbf{I}$. Check this with:

In [38]:
P.T @ P

Out[38]:
array([[1., 0.],
[0., 1.]])

$\mathbf{P}$ can be used to find the matrix $\Lambda$ that decouples the differential equations.

In [39]:
Lam = P.T @ Ktilde @ P
Lam

Out[39]:
array([[ 2.95580264e+01, -5.32907052e-15],
[ 0.00000000e+00,  1.39051349e+02]])

# 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$.

In [40]:
S = la.inv(L.T) @ P
S

Out[40]:
array([[ -36.11302491,  -27.9138454 ],
[  67.25977163, -161.6505027 ]])
In [41]:
u1 = S[:, 0]
u2 = S[:, 1]

In [42]:
u1

Out[42]:
array([-36.11302491,  67.25977163])
In [43]:
u2

Out[43]:
array([ -27.9138454, -161.6505027])

Define the initial coordinates as a scalar factor of the second eigenvector, which sets these values to small angles.

In [44]:
c0 = S[:, 1] / 400

Out[44]:
array([ -3.99836383, -23.1547289 ])

Set the initial speeds to zero:

In [45]:
s0 = np.zeros(2)
s0

Out[45]:
array([0., 0.])

The initial mass normalized coordinates and speeds are then:

In [46]:
q0 = L.T @ c0
q0

Out[46]:
array([-0.00096039, -0.00230817])
In [47]:
qd0 = L.T @ s0
qd0

Out[47]:
array([0., 0.])

Calculate the modal freqencies in radians per second.

In [48]:
ws = np.sqrt(evals)
ws

Out[48]:
array([ 5.43672939, 11.79200359])

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)}$$
In [49]:
phi1 = np.arctan2(ws * P[:, 0] @ q0, P[:, 0] @ qd0)
phi1

Out[49]:
-1.5707963267948966
In [50]:
phi2 =  np.arctan2(ws * P[:, 1] @ q0, P[:, 1] @ qd0)
phi2

Out[50]:
1.5707963267948966

All $\phi$'s can be calculated in one line using NumPy's broadcasting feature:

In [51]:
phis = np.arctan2(ws * P.T @ q0, P.T @ qd0)
phis

Out[51]:
array([-1.57079633,  1.57079633])

The phase shifts for this particular initial condition are $\pm90$ degrees.

In [52]:
np.rad2deg(phis)

Out[52]:
array([-90.,  90.])

Now calculate the modal participation factors.

$$c_i = \frac{\hat{q}^T_{0i} \bar{q}(0)}{\sin\phi_i}$$
In [53]:
cs = P.T @ q0 / np.sin(phis)
cs

Out[53]:
array([-0.    ,  0.0025])

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¶

In [54]:
t = np.linspace(0, 5, num=500)

In [55]:
cs[1] * np.sin(ws[1] * t)

Out[55]:
array([ 0.00000000e+00,  2.94704029e-04,  5.85298503e-04,  8.67731172e-04,
1.13806360e-03,  1.39252609e-03,  1.62757023e-03,  1.83991842e-03,
2.02660951e-03,  2.18504017e-03,  2.31300113e-03,  2.40870801e-03,
2.47082620e-03,  2.49848950e-03,  2.49131214e-03,  2.44939421e-03,
2.37332024e-03,  2.26415106e-03,  2.12340900e-03,  1.95305666e-03,
1.75546956e-03,  1.53340299e-03,  1.28995358e-03,  1.02851618e-03,
7.52736439e-04,  4.66460024e-04,  1.73678969e-04, -1.21523984e-04,
-4.15032324e-04, -7.02753167e-04, -9.80674337e-04, -1.24492031e-03,
-1.49180626e-03, -1.71788943e-03, -1.92001717e-03, -2.09537087e-03,
-2.24150528e-03, -2.35638260e-03, -2.43840089e-03, -2.48641645e-03,
-2.49975971e-03, -2.47824459e-03, -2.42217113e-03, -2.33232125e-03,
-2.20994788e-03, -2.05675747e-03, -1.87488622e-03, -1.66687027e-03,
-1.43561033e-03, -1.18433125e-03, -9.16537048e-04, -6.35962019e-04,
-3.46518698e-04, -5.22432818e-05,  2.42760651e-04,  5.34379361e-04,
8.18546319e-04,  1.09129890e-03,  1.34883366e-03,  1.58755936e-03,
1.80414704e-03,  1.99557645e-03,  2.15917817e-03,  2.29267082e-03,
2.39419289e-03,  2.46232869e-03,  2.49612809e-03,  2.49511976e-03,
2.45931776e-03,  2.38922134e-03,  2.28580798e-03,  2.15051974e-03,
1.98524317e-03,  1.79228300e-03,  1.57433001e-03,  1.33442347e-03,
1.07590882e-03,  8.02390957e-04,  5.17684000e-04,  2.25758101e-04,
-6.93159241e-05, -3.63423360e-04, -6.52462971e-04, -9.32404189e-04,
-1.19934332e-03, -1.44955799e-03, -1.67955902e-03, -1.88613912e-03,
-2.06641760e-03, -2.21788054e-03, -2.33841582e-03, -2.42634262e-03,
-2.48043483e-03, -2.49993816e-03, -2.48458062e-03, -2.43457638e-03,
-2.35062274e-03, -2.23389039e-03, -2.08600714e-03, -1.90903517e-03,
-1.70544231e-03, -1.47806758e-03, -1.23008166e-03, -9.64942630e-04,
-6.86347781e-04, -3.98182026e-04, -1.04463747e-04,  1.90711247e-04,
4.83226832e-04,  7.69003969e-04,  1.04405759e-03,  1.30455215e-03,
1.54685513e-03,  1.76758771e-03,  1.96367182e-03,  2.13237315e-03,
2.27133919e-03,  2.37863213e-03,  2.45275577e-03,  2.49267651e-03,
2.49783764e-03,  2.46816721e-03,  2.40407896e-03,  2.30646658e-03,
2.17669124e-03,  2.01656262e-03,  1.82831367e-03,  1.61456945e-03,
1.37831056e-03,  1.12283157e-03,  8.51695034e-04,  5.68681880e-04,
2.77738634e-04, -1.70775905e-05, -3.11655673e-04, -6.01887815e-04,
-8.83726818e-04, -1.15324253e-03, -1.40667663e-03, -1.64049506e-03,
-1.85143731e-03, -2.03656184e-03, -2.19328714e-03, -2.31942775e-03,
-2.41322465e-03, -2.47336990e-03, -2.49902477e-03, -2.48983152e-03,
-2.44591835e-03, -2.36789760e-03, -2.25685727e-03, -2.11434576e-03,
-1.94235037e-03, -1.74326950e-03, -1.51987928e-03, -1.27529483e-03,
-1.01292678e-03, -7.36433784e-04, -4.49671449e-04, -1.56638587e-04,
1.38578551e-04,  4.31863255e-04,  7.19125760e-04,  9.96360281e-04,
1.25970087e-03,  1.50547532e-03,  1.73025639e-03,  1.93090957e-03,
2.10463683e-03,  2.24901557e-03,  2.36203250e-03,  2.44211162e-03,
2.48813626e-03,  2.49946461e-03,  2.47593871e-03,  2.41788662e-03,
2.32611785e-03,  2.20191209e-03,  2.04700136e-03,  1.86354583e-03,
1.65410374e-03,  1.42159568e-03,  1.16926392e-03,  9.00627137e-04,
6.19431390e-04,  3.29597866e-04,  3.51682016e-05, -2.59751872e-04,
-5.51049787e-04, -8.34663484e-04, -1.10663806e-03, -1.36318091e-03,
-1.60071463e-03, -1.81592689e-03, -2.00581661e-03, -2.16773584e-03,
-2.29942668e-03, -2.39905272e-03, -2.46522473e-03, -2.49701994e-03,
-2.49399499e-03, -2.45619206e-03, -2.38413830e-03, -2.27883847e-03,
-2.14176095e-03, -1.97481725e-03, -1.78033533e-03, -1.56102719e-03,
-1.31995102e-03, -1.06046854e-03, -7.86198152e-04, -5.00964481e-04,
-2.08745017e-04,  8.63853318e-05,  3.80311064e-04,  6.68933477e-04,
9.48227821e-04,  1.21429942e-03,  1.46343800e-03,  1.69216939e-03,
1.89730401e-03,  2.07598131e-03,  2.22570970e-03,  2.34440127e-03,
2.43040089e-03,  2.48250933e-03,  2.49999995e-03,  2.48262885e-03,
2.43063827e-03,  2.34475319e-03,  2.22617126e-03,  2.07654607e-03,
1.89796410e-03,  1.69291560e-03,  1.46425993e-03,  1.21518560e-03,
9.49165896e-04,  6.69910367e-04,  3.81313147e-04,  8.73986342e-05,
-2.07734625e-04, -4.99971090e-04, -7.85235614e-04, -1.05955027e-03,
-1.31908983e-03, -1.56023510e-03, -1.77962337e-03, -1.97419535e-03,
-2.14123779e-03, -2.27842134e-03, -2.38383302e-03, -2.45600288e-03,
-2.49392456e-03, -2.49706923e-03, -2.46539305e-03, -2.39933774e-03,
-2.29982440e-03, -2.16824074e-03, -2.00642163e-03, -1.81662361e-03,
-1.60149333e-03, -1.36403072e-03, -1.10754714e-03, -8.35619153e-04,
-5.52038719e-04, -2.60760278e-04,  3.41543841e-05,  3.28592774e-04,
6.18449040e-04,  8.99681227e-04,  1.16836764e-03,  1.42076153e-03,
1.65334335e-03,  1.86286980e-03,  2.04641912e-03,  2.20143177e-03,
2.32574613e-03,  2.41762869e-03,  2.47579817e-03,  2.49944342e-03,
2.48823471e-03,  2.44232835e-03,  2.36236447e-03,  2.24945817e-03,
2.10518387e-03,  1.93155343e-03,  1.73098809e-03,  1.50628466e-03,
1.26057656e-03,  9.97290111e-04,  7.20096764e-04,  4.32861892e-04,
1.39590896e-04, -1.55626652e-04, -4.48674034e-04, -7.35464797e-04,
-1.01199973e-03, -1.27442265e-03, -1.51907414e-03, -1.74254261e-03,
-1.94171187e-03, -2.11380457e-03, -2.25642092e-03, -2.36757218e-03,
-2.44570839e-03, -2.48973996e-03, -2.49905288e-03, -2.47351729e-03,
-2.41348927e-03, -2.31980590e-03, -2.19377356e-03, -2.03714973e-03,
-1.85211848e-03, -1.64126002e-03, -1.40751470e-03, -1.15414202e-03,
-8.84675200e-04, -6.02871857e-04, -3.12661653e-04, -1.80914804e-05,
2.76730972e-04,  5.67694498e-04,  8.50741701e-04,  1.12192558e-03,
1.37746455e-03,  1.61379521e-03,  1.82762200e-03,  2.01596317e-03,
2.17619236e-03,  2.30607523e-03,  2.40380060e-03,  2.46800573e-03,
2.49779528e-03,  2.49275385e-03,  2.45295175e-03,  2.37894401e-03,
2.27176263e-03,  2.13290223e-03,  1.96429917e-03,  1.76830458e-03,
1.54765153e-03,  1.30541696e-03,  1.04497876e-03,  7.69968661e-04,
4.84221586e-04,  1.91722192e-04, -1.03450709e-04, -3.97181021e-04,
-6.85372768e-04, -9.64007206e-04, -1.22919886e-03, -1.47724973e-03,
-1.70470080e-03, -1.90838036e-03, -2.08544815e-03, -2.23343501e-03,
-2.35027732e-03, -2.43434575e-03, -2.48446798e-03, -2.49994508e-03,
-2.48056123e-03, -2.42658673e-03, -2.33877423e-03, -2.21834825e-03,
-2.06698810e-03, -1.88680445e-03, -1.68030990e-03, -1.45038395e-03,
-1.20023284e-03, -9.33344870e-04, -6.53441693e-04, -3.64426475e-04,
-7.03294435e-05,  2.24748310e-04,  5.16692019e-04,  8.01430618e-04,
1.07499352e-03,  1.33356597e-03,  1.57354226e-03,  1.79157599e-03,
1.98462676e-03,  2.15000252e-03,  2.28539717e-03,  2.38892267e-03,
2.45913539e-03,  2.49505623e-03,  2.49618429e-03,  2.46250384e-03,
2.39448455e-03,  2.29307491e-03,  2.15968906e-03,  1.99618702e-03,
1.80484877e-03,  1.58834247e-03,  1.34968723e-03,  1.09221103e-03,
8.19504279e-04,  5.35369799e-04,  2.43769754e-04, -5.12295839e-05,
-3.45514541e-04, -6.34981406e-04, -9.15593653e-04, -1.18343823e-03,
-1.43478014e-03, -1.66611448e-03, -1.87421537e-03, -2.05618092e-03,
-2.20947366e-03, -2.33195599e-03, -2.42191991e-03, -2.47811092e-03,
-2.49974544e-03, -2.48652180e-03, -2.43862438e-03, -2.35672111e-03,
-2.24195410e-03, -2.09592373e-03, -1.92066637e-03, -1.71862591e-03,
-1.49261975e-03, -1.24579947e-03, -9.81606906e-04, -7.03726141e-04,
-4.16032135e-04, -1.22536691e-04,  1.72667490e-04,  4.65463876e-04,
7.51769514e-04,  1.02759196e-03,  1.28908496e-03,  1.53260207e-03,
1.75474752e-03,  1.95242356e-03,  2.12287367e-03,  2.26372097e-03,
2.37300138e-03,  2.44919103e-03,  2.49122748e-03,  2.49852454e-03,
2.47098045e-03,  2.40897931e-03,  2.31338570e-03,  2.18553265e-03,
2.02720304e-03,  1.84060470e-03,  1.62833971e-03,  1.39336803e-03,
1.13896627e-03,  8.68681981e-04,  5.86284191e-04,  2.95710851e-04,
1.01391495e-06, -2.93697160e-04, -5.84312719e-04, -8.66780220e-04,
-1.13716074e-03, -1.39168391e-03, -1.62680048e-03, -1.83923183e-03,
-2.02601566e-03, -2.18454734e-03, -2.31261618e-03, -2.40843632e-03,
-2.47067156e-03, -2.49845405e-03, -2.49139639e-03, -2.44959698e-03,
-2.37363870e-03, -2.26458078e-03, -2.12394398e-03, -1.95368944e-03,
-1.75619132e-03, -1.53420365e-03, -1.29082200e-03, -1.02944023e-03,
-7.53703241e-04, -4.67456095e-04, -1.74690420e-04,  1.20511258e-04,
4.14032444e-04,  7.01780077e-04,  9.79741606e-04,  1.24404094e-03,
1.49099252e-03,  1.71715267e-03,  1.91936766e-03,  2.09481767e-03,
2.24105610e-03,  2.35604370e-03,  2.43817701e-03,  2.48631070e-03,
2.49977356e-03,  2.47837786e-03,  2.42242195e-03,  2.33268613e-03,
2.21042173e-03,  2.05733368e-03,  1.87555676e-03,  1.66762578e-03])

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.

In [56]:
# cs[1] * np.sin(ws[1] * t) * u1

In [57]:
c1 = cs[1] * np.sin(ws[1] * t) * np.tile(u1, (len(t), 1)).T
c1.shape

Out[57]:
(2, 500)

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:

In [58]:
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

-0.0 5.436729388358913 -1.5707963267948966 [-36.11302491  67.25977163]
0.0025000000000000005 11.792003585385926 1.5707963267948966 [ -27.9138454 -161.6505027]

In [59]:
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:

In [60]:
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.set_xlabel('Time [s]')
ax.set_ylabel('Angle [deg]')
ax.legend([r'$\theta$', r'$\phi$'])

Out[60]:
<matplotlib.legend.Legend at 0x7f5d38729910>

Simulate and plot the second mode:

In [61]:
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.set_xlabel('Time [s]')
ax.set_ylabel('Angle [deg]')
ax.legend([r'$\theta$', r'$\phi$'])

Out[61]:
<matplotlib.legend.Legend at 0x7f5d37ee3ee0>

Compare this to the free response from the system:

In [62]:
sys.coordinates['theta'] = c0[0]
sys.coordinates['phi'] = c0[1]

sys.speeds['alpha'] = 0
sys.speeds['beta'] = 0

In [63]:
traj = sys.free_response(5.0)

In [64]:
traj[['theta', 'phi']].plot()

Out[64]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f5d37eb4c40>
In [65]:
sys.animate_configuration(fps=30, repeat=False)

Out[65]:
<matplotlib.animation.FuncAnimation at 0x7f5d37d64700>

Simulate with arbitrary initial conditions.

In [66]:
sys.coordinates['theta'] = np.deg2rad(12.0)

traj[['theta', 'phi']].plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7f5d37d41e50>