In [1]:
from IPython.display import Image
import numpy as np
from scipy import linalg as LA
import sympy as sp
import math
In [2]:
Image(filename='two_springs_two_masses.png')
Out[2]:
In [3]:
# Construct mass matrix
m = 20
m1 = m
m2 = 2 * m
mMatrix = np.array([[m1, 0], [0, m2]])
mMatrix
Out[3]:
array([[20,  0],
       [ 0, 40]])
In [4]:
# Construct spring stiffness matrix
k = 1000
k1 = k
k2 = 2 * k
kMatrix = np.array([[k1 + k2, -k2], [-k2, k2]])
kMatrix
Out[4]:
array([[ 3000, -2000],
       [-2000,  2000]])
In [5]:
# Find the eigenvalues and eigenvactors
values, vectors = LA.eig(a = kMatrix, b = mMatrix)
vectors = vectors.T
vectors = np.array([[vectors[0,0] / vectors[0, 1] , 1.0], [vectors[1,0] / vectors[1, 1], 1.0]])
values, vectors
Out[5]:
(array([ 186.60254038+0.j,   13.39745962+0.j]),
 array([[-2.73205081,  1.        ],
        [ 0.73205081,  1.        ]]))
In [6]:
# Ordering by the magnitude of frequency
arg = np.argsort(values)
values = np.array(values)[arg]
vectors = np.array(vectors)[arg]
values, vectors
Out[6]:
(array([  13.39745962+0.j,  186.60254038+0.j]),
 array([[ 0.73205081,  1.        ],
        [-2.73205081,  1.        ]]))
In [7]:
# Natural frequency
natural_frequencies = np.sqrt(values)
natural_frequencies
Out[7]:
array([  3.66025404+0.j,  13.66025404+0.j])
In [8]:
# Check for orthogonality with respect to M
a = np.dot(np.dot(vectors, mMatrix), vectors.T)
a
Out[8]:
array([[  5.07179677e+01,  -7.10542736e-15],
       [ -7.10542736e-15,   1.89282032e+02]])
In [9]:
# Drop small off-diagonal terms to make the mode normalization easier
indexes = a < 10e-8
a[indexes] = 0
a
Out[9]:
array([[  50.7179677,    0.       ],
       [   0.       ,  189.2820323]])
In [10]:
# Construct Normal Modal Matrix
phi = np.dot(vectors.T, LA.inv(np.sqrt(a)))
phi
Out[10]:
array([[ 0.10279223, -0.19857935],
       [ 0.1404168 ,  0.07268509]])
In [11]:
# Initial condition in principal corrdinates
y0 = np.dot(phi.T, np.dot(mMatrix, np.array([1, -1])))
y0
Out[11]:
array([-3.56082742, -6.87899034])
In [12]:
# Response in principal coordinates
t = sp.Symbol('t')
y = y0 * np.array([sp.cos(natural_frequencies[0] * t), sp.cos(natural_frequencies[1] * t)])
y
Out[12]:
array([-3.5608274213041*cos(3.66025403784439*t),
       -6.87899033839188*cos(13.6602540378444*t)], dtype=object)
In [13]:
# Response in x corrdinates
x = np.dot(phi, y)
x
Out[13]:
array([ -0.366025403784439*cos(3.66025403784439*t) + 1.36602540378444*cos(13.6602540378444*t),
       -0.5*cos(3.66025403784439*t) - 0.5*cos(13.6602540378444*t)], dtype=object)
In [13]:
In [ ]: