自由落下運動のシミュレーションと可視化

最初は高校生の物理で一番最初に習うであろう、自由落下運動を例にシミュレーションしてみよう。
簡単すぎると思うかもしれないが、初めての人にもとっつきやすいように誰でもわかるようなものからシミュレーションしていく。
そのうち、複雑なシミュレートもしていくから、とりあえずこれで我慢してほしい。

ということで、以下の問題について考えよう。

初期位置 $(x_0,y_0)$ から初期速度 $(v_{x0},v_{y0})$ で質点を投げ出したとき、
その質点はどのような動くだろうか？
ただし、空気抵抗は考えないこととする。

自由落下運動では、ある質点に下向きに $g$ 加速度が働いている。
上向きを正とすることによって、時刻 $t$ における変位は

\begin{equation}
y=v_0t-\frac{1}{2}gt^2
\end{equation}

となる。
上の式は、 $t=0$ のときの位置を考慮していないため、それを $y_0$ とすることで、正しくは次のようにあらわされる。

\begin{equation}
y=y_0+v_0t-\frac{1}{2}gt^2
\end{equation}

大学では、ベクトル表記を使うことの方が多い。早めになれるために、今のうちから使っておこう。このベクトルはプログラムにおける配列の考え方と一致してる。
上記の加速度は $x$ 方向に $0$ , $y$ 方向に $-g$ と考えることができるため、初期位置、初期速度、加速度をベクトルで表現すると次のようにあらわせる。

なお、大学風に、ベクトルは $\vec{x}$ のように上付き矢印ではなく、 $\boldsymbol{x}$ のように、太字で表記する。
(筆者も最初はこの表記に違和感を覚えていたが、高校で学習するベクトル＝矢印というイメージが必ずしも本質的ではないことを理解して、解消されてきた)。
\begin{equation}
\boldsymbol{a}=\left( \begin{array}{cc} 0\\ -g\\ \end{array} \right),\ \
\boldsymbol{x}_0=\left( \begin{array}{cc} x_0\\ y_0\\ \end{array} \right),\ \
\boldsymbol{v}_0=\left( \begin{array}{cc} v_{x0}\\ v_{y0}\\ \end{array} \right),\ \
\end{equation}

これにより、等加速度運動の式をベクトルで表記すると次のようにあらわすことができる。

\begin{equation}
\boldsymbol{x}=\boldsymbol{x}_0+\boldsymbol{v}_0t+\frac{1}{2}\boldsymbol{a}t^2
\end{equation}

さて、長々と書いてしまったが、上式をそのままプラグラムにぶち込めば、自由落下のシミュレーションは完成だ。めちゃくちゃ簡単！！
念のため、コードを全文書いておく

# -*- coding utf-8 -*-

import numpy as np

First_pos=np.array([0.0,0.0])   #初期位置(m)
First_vel=np.array([4.0,5.0])   #初期速度(m/s)

G=-9.80665      #重力加速度(m/s^2)

def main():
    print("|*時刻|*位置[x,y]|*速度[vx,vy]|")
    PM1=PointMass(First_pos,First_vel)
    for t in np.arange(0,5,0.2):
        xy,v=PM1.fall(t)
        xy=np.array2string(xy,precision=5,separator=',')
        v =np.array2string(v,precision=5,separator=',')
        print("| %.2f |" % (t),xy,"|",v,"|")

class PointMass:
    def __init__(self,xy0,v0):
        self.xy0=xy0
        self.v0=v0
        self.a=np.array([0.0,G])
    def fall(self,t):
        xy=self.xy0+self.v0*t+0.5*self.a*t**2
        v=self.v0+self.a*t
        return xy,v

if __name__=='__main__':
    main()

この結果は次のようになる。

時刻	位置[x,y]	速度[vx,vy]
0.00	[0.,0.]	[4.,5.]
0.20	[0.8 ,0.80387]	[4. ,3.03867]
0.40	[1.6 ,1.21547]	[4. ,1.07734]
0.60	[2.4 ,1.2348]	[ 4. ,-0.88399]
0.80	[3.2 ,0.86187]	[ 4. ,-2.84532]
1.00	[4. ,0.09668]	[ 4. ,-4.80665]

それぞれのデータを"|"でわけていたり(ブログ掲載用)、無駄にclassを使っていたり(自分の練習のため)、
学習するにあたって無駄なところが多いが、重要なのは25,26行目だ。

念のため、簡単に説明すると、5,6行目で初期位置、初期速度を与えている。
初期位置は原点、初期速度は $(4,5)$ を与えている。
そして、25行目で速度を、26行目で位置を計算している。
上の式をそのままぶち込んでいるだけだ。

上の解析結果の数値による情報は重要な時も多いが、どう言う振る舞いであるのかをイメージしにくい思う。
数値を見ただけでは質点がどのようなふるまいをするかは(この問題は簡単すぎてできると思うけど)イメージするのは難しい。

そこで、上のデータを何らかの図によってあらわそう。このことを可視化(visualization)という。
実はこの記事のテーマとして、可視化の重要性について論じようとも思っていた。

上のプログラムをmatplotlibパッケージを使ってアニメーションにしてみた。

f:id:atily529:20181005024657g:plain — 自由落下の結果1

おお～アニメーションを作ったのは初めてなので、なんかうれしい。
このように図で表わすことで、質点の動きがとてもイメージしやすくなったと思う(こんなことしなくてもわかるとか言わないで)。

ここに、更なる情報を加えることもできる。例えば、軌跡を表示してみよう。

f:id:atily529:20181005025408g:plain — 自由落下の結果2

このように軌跡を表示するだけで、なんとなく質点が放物線上に動きそうって予想できそう！

さらに、軌跡に速いところでは赤色に、遅くなるにつれて青色にするという情報も付け加えてみよう。

f:id:atily529:20181005042207g:plain — 自由落下の結果3

ここまでくれば、頂点で速度が最も遅くなるだろうと予測できそう！

ここまで説明すれば、可視化の重要性について説明できたと思う。そして、単に可視化といっても色々な工夫の仕方があるということも。
上のプログラムを一応アップロードした。
free-fall_animation1.py - Google ドライブ
色々稚拙で恥ずかしいが、いろいろ試してみてもいいかも？

pythonと、numpy，matplotlibのパッケージ、pillowがあれば動かせると思う。
暇ならやってみてね(修正案とかももしあったら教えてください。上のプログラムに反映するかはわからないけど、参考にさせていただきます！)