はじめに
一般的に運動方程式は位置ベクトルの時間に関する 2階微分方程式になるので、運動をシミュレートするためにはこれを解かないといけないですよね。しかし一般解が簡単に求められない場合もあり、その場合は数値計算を行うことになります。
そのような場合に、微分方程式の数値解を求める方法が SciPy では ode というモジュールで実装されているので、それを用いて計算してみます。
scipy.integrate.ode
http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html
問題設定
今回は簡単のため、質量 m、位置 x で記述される質点に関して、高さ h からの自由落下問題を考えます。ここで、空気抵抗は無視できるとし、質点に加わる力は重力のみとします。
運動方程式
以上の仮定より運動方程式は、
m \ddot{x} = - g
となります。但、g は重力定数です。また、物理学の慣習で、変数の上のドットは時間微分を表します。ドットが2つついている場合は２階微分です。
さて、これは明らかに解析解を求めることができて、
x = - \frac{1}{2} g t^2 + h
となるのですが、これに時刻を突っ込んでしまうと元も子もないので、数値解との比較にのみ使用します。
scipy で解けるように式変形
さて、本題です。
SciPy の odeモジュールは（おそらく他のものと同じように）、１階の常微分方程式しか解けないので、以下の変数変換を与えます。
v := \dot{x}
従って運動方程式は、
\dot{x} = v
\dot{v} = - \frac{g}{m}
となります。ここで、微分を左辺に寄せていますが、これは偶然とか見栄えの問題ではなくて、odeモジュールの制約です。線形微分方程式に関しては、一般的に次のような形にならないといけません。
\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}
今回の取り扱う問題では、
\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}
です。 以下に示すコード中で「第一式」とか「第二式」とか出てくるのはこれを展開した式のことです 。
ここで「高さ h からの落下である」、「自由落下である」という記述から、以下の初期条件を得ます。
x_1(0) = h
x_2(0) = 0
これでどちらも1階常微分方程式の形になったので、やっと odeモジュールを使用して解けます。
scipy で常微分方程式を解く
さて、これを python のコードに起こしていきます。
と言っても、これはほとんど機械的にできます。コードの全容を示します。
#-*- coding:utf-8 -*-
import numpy as np
from scipy.integrate import odeint
g = 9.8 # 重力定数
m = 1.0 # 質量
h = 10 # 初期位置
def f(x, t):
ret = [
x[1], # 第一式の右辺
-g / m # 第二式の右辺
]
return ret
def main():
# 初期状態
x0 = [
h, # 第一式の初期条件
0 # 第二式の初期条件
]
# 計算するインターバル
# 引数は、「開始時刻」、「終了時刻」、「刻み」の順
t = np.arange(0, 10, 0.1)
# 積分する
x = odeint(f, x0, t)
# 結果を表示する（とりあえずそのまま print）
print(x)
if __name__ == '__main__':
main()
これを実行すると、以下の様な出力が得られます。
$ 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]]
2重配列の形で与えられていて、時刻ごとに解の組が格納されています。
解の組は、第一要素が第一式の解、第二要素が第二式の解です（要するに見たままです）。
試しに t = 10 における数値解を解析解から計算される値と比較してみると、
数値解 = -4.70249000e+02 = -470.2・・・
解析解 = -9.8 / 2 * (10 * 10) + 10 = -480.0・・・
となり、誤差は 10[m] 程度出ています（※ koryorさんよりご指摘いただきまして、現在修正中です。詳細はコメント欄をご覧ください）。この精度をあげるためにはいくつか方法がありますが、簡単なものとしては、刻み幅を小さくするという方法があります。コードで言うと、この部分です：
# 計算するインターバル
# 引数は、「開始時刻」、「終了時刻」、「刻み」の順
t = np.arange(0, 10, 0.01) # <= 変更
この変更により、数値解は 479 となり、誤差は 1.0 程度になります。
数値解 = -4.79020490e+02 = -479.0・・・
解析解 = -9.8 / 2 * (10 * 10) + 10 = -480.0・・・
もし刻み幅を小さくしても十分な精度が出なかったり、時間がかかりすぎてしまう場合には、別のアルゴリズムを検討してください。
考察
ちなみに、真面目に出てきた解を見つめると、質点の位置が負になっていることがわかります。
これは物理を志した学生が一度は体験する、地面に質点がめり込んで突き進んでいる状態です。
境界条件や地面の反発を無視しているため、質点は地面に向かって突き刺さりさらに速度を上げて地球の反対側から突き抜けますが、こうなってしまうことを防ぐためには空気抵抗や地面の反発力を定義してあげる必要があります。今回は簡単なケースなので直ぐに気づきますが、複雑なケースだと見落としがちなので解き終わったら冷静に結果を見つめることを忘れないようにしてください。
まとめ
scipy で常微分方程式の数値解を求めました。
いったん使い方を覚えてしまえば非常に簡単で、また使い方も単純なので、運動方程式を解きたくなったら思い出してください。