Scipyでモード解析における一般固有値問題を解く

機械振動学におけるモード解析では、剛性行列Kと質量行列Mから一般固有値問題

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

を解く必要があり、これは両辺に行列が2つある一般固有値問題と呼ばれる。なお、

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

は標準固有値問題と呼ばれる。一般固有値問題を解くには、両辺に左から$M^{-1}$を掛けて標準固有値問題にして解く方法もあるが、scipy.linalg.eighを使えば、いきなり解くことができる。なお、K,Mは実対称行列なのでscipy.linalg.eigでなくscipy.linalg.eighの方を使う（なぜだかわからないがscipy.linalg.eigだとうまくいかない）

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

    omega = np.sqrt(w)
    # 昇順にソートする
    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)

正規化は、$\phi_i^T M\phi_i=1$となるようにするが、scipy.linalg.eighでは既にそうなるように出力されるようである。もし自分で計算する場合は$M=diag(m_1,m_2,\cdots,m_n)$なので$\phi_i = \phi_i / \sqrt{m_i}$のように新しい固有ベクトル（固有振動モード）を計算すればよい。

上記により、剛性行列Kも、質量行列Mも固有ベクトルを横に並べた行列[$\Phi$]により対角化され、以下のようになる

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

計算結果：

[[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]]