To find the moment of inertia and the principal axis of inertia from the inertia tensor, the operation of diagonalizing the matrix is performed. I think it's quite difficult to write a diagonalization program, but using NumPy makes it surprisingly easy. I'll read the theoretical story elsewhere, and only explain the procedure for calculations using NumPy.
[Moment of Inertia-Wikipedia](https://ja.wikipedia.org/wiki/%E6%85%A3%E6%80%A7%E3%83%A2%E3%83%BC%E3%83%A1% E3% 83% B3% E3% 83% 88 # .E6.85.A3.E6.80.A7.E4.B8.BB.E8.BB.B8.E3.81.A8.E4.B8.BB.E6 .85.A3.E6.80.A7.E3.83.A2.E3.83.BC.E3.83.A1.E3.83.B3.E3.83.88) Coordinate transformation of inertia tensor
Set the output of NumPy so that it is easy to read. Here, it is set to be fixed and displayed with 4 digits after the decimal point.
import numpy as np
np.set_printoptions(precision=4, suppress=True)
In particular, the readability of the minimum value differs greatly.
-2.86013404e-15
→ 0.
Define the matrix representing the inertia tensor as ndarray. However, since the inertia tensor is a symmetric matrix, make it symmetric.
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 is a function that finds eigenvalues and eigenvectors. This can be used to determine the moment of inertia and the principal axis of inertia.
I_p, P = np.linalg.eig(I)
print('Main moment of inertia: \n', I_p)
print('Moment of inertia spindle: \n', P.T)
output
Main moment of inertia:
[ 33.8923 18.5542 7.5536]
Moment of inertia spindle:
[[ 0.8716 0.4103 0.2683]
[ 0.4706 -0.8535 -0.2238]
[-0.1371 -0.3213 0.937 ]]
The first return value is an array of eigenvalues. This is the main moment of inertia. The second return value is an arrangement of three eigenvectors as column vectors, and each vector serves as the main axis of inertia. Since the row vector is easier to understand as it looks, it is transposed and output.
We have already obtained what we need, but let's confirm that diagonalization gives the principal moment of inertia. The principal moment of inertia can be obtained by diagonalizing the inertia tensor $ I $ using the matrix $ P $. Diagonal matrix of principal moment of inertia = $ P ^ {-1} I P $ Reference: Diagonalization -Wikipedia
numpy.linalg.inv is a function to find the inverse matrix.
The @
operator represents the product of matrices. (Note that in NumPy the *
operator is the product of the elements.)
I_p = np.linalg.inv(P) @ I @ P
print(I_p)
output
[[ 33.8923 -0. 0. ]
[ -0. 18.5542 0. ]
[ -0. -0. 7.5536]]
The diagonal component obtained by this is the main moment of inertia. Only diagonal components can be extracted using the diag function.
print(np.diag(I_p))
output
[ 33.8923 18.5542 7.5536]
Since the @
operator that calculates the matrix product can be used in Python 3.5 or later, the numpy.dot function is used instead in Python 3.4 and earlier versions.