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.

Recommended Posts

Find the inertial spindle and moment of inertia from the inertial tensor with NumPy
See the power of speeding up with NumPy and SciPy
Find the distance from latitude and longitude (considering the roundness of the earth).
Find the waypoint from latitude and longitude (considering the roundness of the earth).
Find the general terms of the Tribonacci sequence with linear algebra and Python
Find a position above the threshold with NumPy
Deep Learning from scratch The theory and implementation of deep learning learned with Python Chapter 3
Find out the day of the week with datetime
Extract images and tables from pdf with python to reduce the burden of reporting
Let's transpose the matrix with numpy and multiply the matrices.
Learn Nim with Python (from the beginning of the year).
Find the sum of unique values with pandas crosstab
Visualize the range of interpolation and extrapolation with python
Find out the location of packages installed with pip
Find the "minimum passing score" from the "average score of examinees", "average score of successful applicants", and "magnification" of the entrance examination
I tried to find the entropy of the image with python
I tried to find the average of the sequence with TensorFlow
Play with the password mechanism of GitHub Webhook and Python
A story about predicting prefectures from the names of cities, wards, towns and villages with Jubatus
Find the white Christmas rate by prefecture with Python and map it to a map of Japan
Find the critical path of PERT using breadth-first search and depth-first search
I compared the speed of Hash with Topaz, Ruby and Python
Find out the age and number of winnings of prefectural governors nationwide
Find the optimal value of a function with a genetic algorithm (Part 2)
Find the transfer function of one degree of freedom system with PythonControl.
Capture GeneratorExit and detect the end of iteration from the generator side
Find the numerical solution of the second-order ordinary differential equation with scipy
[Required subject DI] Implement and understand the mechanism of DI with Go
Find the smallest index that meets the cumulative sum threshold with numpy
To improve the reusability and maintainability of workflows created with Luigi
From the introduction of JUMAN ++ to morphological analysis of Japanese with Python
Operate Firefox with Selenium from python and save the screen capture
Find the intersection of a circle and a straight line (sympy matrix)
What I did when I couldn't find the feature point with the optical flow of opencv and when I lost it
[Completed version] Try to find out the number of residents in the town from the address list with Python