Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
4
Help us understand the problem. What is going on with this article?
@alchemist

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

More than 3 years have passed since last update.

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

参考:

慣性モーメント - Wikipedia
慣性テンソルの座標変換

計算手順

  1. 慣性テンソルを表す行列$I$を作る。
  2. numpy.linalg.eigで固有値と固有ベクトルを束ねた行列$P$を求める。
  3. 行列$P$を使って慣性テンソル$I$を対角化する。(省略可)

0. 準備

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

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

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

1. 慣性テンソルを表す行列Iを作る。

慣性テンソルを表す行列を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]]

2. numpy.linalg.eigで固有値と固有ベクトルを束ねた行列Pを求める。

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

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

3. 行列Pを使って慣性テンソルIを対角化する。(省略可)

すでに必要なものは得られているが、対角化によって主慣性モーメントが得られることを確認してみる。
慣性テンソル$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関数を使う。

4
Help us understand the problem. What is going on with this article?
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
alchemist
CAEの技術サポートとコンサルティングの仕事をしています。 おもにNXの構造解析と機構解析を担当しています。 最近はCAD/CAEのカスタマイズの仕事を始めました。 仕事ではC#とPythonを使っています。

Comments

No comments
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account Login
4
Help us understand the problem. What is going on with this article?