NumPyで慣性テンソルから慣性主軸と主慣性モーメントを求める #Python

慣性テンソルから慣性主軸と主慣性モーメント求めるには行列の対角化という操作を行う。
対角化のプログラムを書くのはかなり大変だと思うけど、NumPyを使うことで案外簡単に求めることができる。
理論的な話は他で読んでもらうことにして、NumPyを使った計算の手順だけを説明する。

計算手順

NumPyの出力を読みやすいように設定しておく。
ここでは小数点以下４桁で固定して表示するように設定する。

import numpy as np
np.set_printoptions(precision=4, suppress=True)

特に極小の値の読みやすさが大きく違ってくる。
-2.86013404e-15　→　0.

慣性テンソルを表す行列をndarrayとして定義する。
ただし慣性テンソルは対称行列なので、対称になるように。

I = np.array([[30, 5, 5],
              [5, 20, 5],
              [5, 5, 10]])
print(I)

output

[[30  5  5]
 [ 5 20  5]
 [ 5  5 10]]

numpy.linalg.eigは固有値と固有ベクトルを求める関数。
これを使って主慣性モーメントと慣性主軸を求めることができる。

I_p, P = np.linalg.eig(I)
print('主慣性モーメント: \n', I_p)
print('慣性主軸: \n', P.T)

output

主慣性モーメント: 
 [ 33.8923  18.5542   7.5536]
慣性主軸: 
 [[ 0.8716  0.4103  0.2683]
 [ 0.4706 -0.8535 -0.2238]
 [-0.1371 -0.3213  0.937 ]]

１つ目の戻り値は固有値の配列。これが主慣性モーメントとなる。
２つ目の戻り値は３つの固有ベクトルを列ベクトルとして並べたもので、それぞれのベクトルが慣性主軸となる。
行ベクトルの方が見た目としてわかりやすいので転置して出力している。

すでに必要なものは得られているが、対角化によって主慣性モーメントが得られることを確認してみる。
慣性テンソル$I$を行列$P$を使って対角化すると主慣性モーメントを得られる。
主慣性モーメントの対角行列 = $P^{-1} I P$
参考：対角化 - Wikipedia

numpy.linalg.invは逆行列を求める関数。
@演算子は行列の積を表す。（NumPyでは*演算子は要素どうしの積になることに注意。）

I_p = np.linalg.inv(P) @ I @ P
print(I_p)

output

[[ 33.8923  -0.       0.    ]
 [ -0.      18.5542   0.    ]
 [ -0.      -0.       7.5536]]

これで得られた対角成分が主慣性モーメント。
diag関数を使うと対角成分だけを取り出せる。

print(np.diag(I_p))

output

[ 33.8923  18.5542   7.5536]

行列の積を計算する@演算子が使えるのはPython 3.5以降なので、Python 3.4以前のバージョンでは代わりにnumpy.dot関数を使う。