Generally, the equation of motion is a second-order differential equation with respect to the time of the position vector, so you have to solve this to simulate the motion. However, there are cases where a general solution cannot be easily obtained, in which case numerical calculations will be performed. In such a case, the method of finding the numerical solution of the differential equation is implemented in a module called ode in SciPy, so let's calculate using it.
scipy.integrate.ode http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html
For the sake of simplicity, let us consider the free fall problem from the height h for the mass points described by the mass m and the position x. Here, the air resistance is negligible, and the force applied to the mass point is only gravity.
From the above assumption, the equation of motion is
m \ddot{x} = - g
It will be. However, g is a gravitational constant. Also, in physics convention, the dots above the variable represent the time derivative. If there are two dots, it is the second derivative.
Well, this can clearly find an analytical solution,
x = - \frac{1}{2} g t^2 + h
However, if the time is inserted into this, there is no source or child, so it is used only for comparison with the numerical solution.
Well, this is the main subject. Since SciPy's ode module (probably like any other) can only solve first-order ODEs, it gives the following change of variables.
v := \dot{x}
Therefore, the equation of motion is
\dot{x} = v
\dot{v} = - \frac{g}{m}
It will be. Here, the derivative is moved to the left side, but this is not a matter of chance or appearance, but a constraint of the ode module. For linear differential equations, it should generally have the following form:
\frac{\mathrm{d}}{\mathrm{d}t}
\begin{bmatrix}
x_1 \\
\vdots \\
x_n
\end{bmatrix} =
\begin{bmatrix}
a_{11} & a_{12} & \dots & a_{1n} \\
\vdots & \ddots & & \vdots \\
a_{n1} & \dots & & a_{nn}
\end{bmatrix}
\begin{bmatrix}
x_1 \\
\vdots \\
x_n
\end{bmatrix} +
\begin{bmatrix}
b_1 \\
\vdots \\
b_n
\end{bmatrix}
In this issue,
\frac{\mathrm{d}}{\mathrm{d}t}
\begin{bmatrix}
x_1 \\
x_2
\end{bmatrix} =
\begin{bmatrix}
0 & 1 \\
0 & 0
\end{bmatrix}
\begin{bmatrix}
x_1 \\
x_2
\end{bmatrix} +
\begin{bmatrix}
0 \\
- \frac{g}{m}
\end{bmatrix}
is. ** In the code shown below, "first formula" and "second formula" appear in the expanded formula **.
Here, the following initial conditions are obtained from the descriptions "falling from height h "and" free fall ".
x_1(0) = h
x_2(0) = 0
Now that both are in the form of first-order ODEs, we can finally solve them using the ode module.
Now let's make this into python code. That said, this can be done almost mechanically. Here is the full code.
#-*- coding:utf-8 -*-
import numpy as np
from scipy.integrate import odeint
g = 9.8 #Gravitational constant
m = 1.0 #mass
h = 10 #Initial position
def f(x, t):
ret = [
x[1], #Right side of equation 1
-g / m #Right side of equation 2
]
return ret
def main():
#initial state
x0 = [
h, #Initial conditions of the first equation
0 #Initial condition of the second equation
]
#Interval to calculate
#The arguments are in the order of "start time", "end time", and "increment".
t = np.arange(0, 10, 0.1)
#Integrate
x = odeint(f, x0, t)
#Display the result (print as it is for the time being)
print(x)
if __name__ == '__main__':
main()
If you do this, you will get output similar to the following:
$ python free_fall_sample.py
[[ 1.00000000e+01 0.00000000e+00]
[ 9.95100000e+00 -9.80000000e-01]
[ 9.80400000e+00 -1.96000000e+00]
[ 9.55900000e+00 -2.94000000e+00]
...
[ -4.60596000e+02 -9.60400000e+01]
[ -4.70249000e+02 -9.70200000e+01]]
It is given in the form of a double array and stores a set of solutions for each time. The set of solutions is that the first element is the solution of the first equation and the second element is the solution of the second equation (in short, as you see it).
As a test, comparing the numerical solution at t = 10 with the value calculated from the analytical solution,
Numerical solution= -4.70249000e+02 = -470.2 ...
Analytical solution= -9.8 / 2 * (10 * 10) + 10 = -480.0 ...
The error is about 10 [m](* Koryor pointed out that it is currently being corrected. Please see the comment section for details). There are several ways to improve this accuracy, but the easiest way is to reduce the step size. In code, this is the part:
#Interval to calculate
#The arguments are in the order of "start time", "end time", and "increment".
t = np.arange(0, 10, 0.01) # <=Change
With this change, the numerical solution is 479 and the error is around 1.0.
Numerical solution= -4.79020490e+02 = -479.0 ...
Analytical solution= -9.8 / 2 * (10 * 10) + 10 = -480.0 ...
If reducing the step size does not provide sufficient accuracy or takes too long, consider another algorithm.
By the way, if you look at the solution that came out seriously, you can see that the position of the mass point is negative. This is a state in which a mass point is sunk into the ground, which a student who aspires to physics experiences once. Since the boundary conditions and the repulsion of the ground are ignored, the mass point pierces toward the ground and further speeds up to penetrate from the other side of the earth, but to prevent this from happening, air resistance and repulsion of the ground Need to be defined. This time it's a simple case, so you'll notice it right away, but if it's a complicated case, it's easy to overlook it, so don't forget to calmly look at the results once you've solved it.
I found the numerical solution of the ordinary differential equation with scipy. Once you learn how to use it, it's very easy and simple to use, so if you want to solve the equation of motion, remember it.
Recommended Posts