[PYTHON] Solve general eigenvalue problems in mode analysis with Scipy

In modal analysis in mechanical vibrations, the general eigenvalue problem from the stiffness matrix K and the mass matrix M

K\{\phi_i\}=\omega_i^2 M\{\phi_i\}

Needs to be solved, which is called the general eigenvalue problem with two matrices on each side. In addition, it should be noted

A\{\phi_i\}=\lambda \{\phi_i\}

Is called the standard eigenvalue problem. To solve the general eigenvalue problem, you can multiply both sides by $ M ^ {-1} $ from the left to make it a standard eigenvalue problem, but you can solve it suddenly by using scipy.linalg.eigh. Since K and M are real symmetric matrices, use scipy.linalg.eigh instead of scipy.linalg.eig (I don't know why, but scipy.linalg.eig doesn't work).

    w, v = scipy.linalg.eigh(K, M)
    print("eigen value:\n", w)

    omega = np.sqrt(w)
    #Sort in ascending order
    omega_v = [[row,list(v[i])] for i, row in enumerate(omega)]
    omega_v.sort(key=lambda x: x[0])
    omega = np.array([row[0] for row in omega_v])
    v = np.array([row[1] for row in omega_v])

    print("eigen vector:\n", v)
    print("omega[rad/s]=",omega)
    f = omega / (2*np.pi)
    print("f[Hz]=",f)
     
    print(v.T @ M @ v)
    print(v.T @ K @ v)

Normalization makes $ \ phi_i ^ T M \ phi_i = 1 $, but scipy.linalg.eigh seems to output so already. If you calculate it yourself, $ M = diag (m_1, m_2, \ cdots, m_n) $, so if you calculate a new eigenvector (natural vibration mode) like $ \ phi_i = \ phi_i / \ sqrt {m_i} $ Good.

From the above, both the stiffness matrix K and the mass matrix M are diagonalized by the matrix [$ \ Phi $] in which the eigenvectors are arranged side by side, and become as follows.

[\Phi]^T M [\Phi] = I \\
[\Phi]^T K [\Phi] = diag(\omega_1^2,\omega_2^2,\cdots,\omega_n^2)

Calculation result:

[[1.00000000e+00 0.00000000e+00 3.16285753e-17 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 3.16285753e-17]
 [2.94753911e-17 0.00000000e+00 1.00000000e+00 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  1.00000000e+00 0.00000000e+00]
 [0.00000000e+00 2.94753911e-17 0.00000000e+00 0.00000000e+00
  0.00000000e+00 1.00000000e+00]]


[[3.86469572e-01 0.00000000e+00 1.29304523e-15 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 9.66173929e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 1.80376128e-14]
 [1.13928941e-15 0.00000000e+00 1.10551030e+02 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 4.80769231e+02
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  1.25000000e+03 0.00000000e+00]
 [0.00000000e+00 2.90139143e-14 0.00000000e+00 0.00000000e+00
  0.00000000e+00 2.76377576e+03]]