[PYTHON] Find the inertial spindle and moment of inertia from the inertial tensor with NumPy

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.

reference:

[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

Calculation procedure

1. Create a matrix $ I $ that represents the inertia tensor.
2. Use numpy.linalg.eig to find the matrix $ P $ that bundles the eigenvalues and eigenvectors.
3. Diagonalize the inertia tensor $ I $ using the matrix $ P $. (Optional)

0. Preparation

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.

1. Create a matrix I that represents the inertia tensor.

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

2. Use numpy.linalg.eig to find the matrix P that bundles the eigenvalues and eigenvectors.

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.

3. Diagonalize the inertia tensor_I_ using the matrix P. (Optional)

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]

Notes

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.