Pythonを使って自由落下及びのシュミレーションをしてみようということでまとめました。

運動方程式の差分化

運動方程式の差分化
$$
F=ma
$$
高校物理でもおなじみの運動方程式である.
大学においては,運動方程式は微分を使った式で書かれることが多く,
$$
F=m\frac{d^2x}{dx^2}
$$
となる.上式より,以下の式が導き出される.
$$
\begin{eqnarray}
\frac{dv(t)}{dt} &=& a(t) \
\frac{dx(t)}{dt} &=& v(t)
\end{eqnarray}
$$
よって,この連立微分方程式をプログラムに落とし込めばシュミレーションが可能となる.
ここで微分の定義を思い出してほしい
$$
\frac{df(x)}{dx}=\frac{f(x+\Delta x)-f(x)}{\Delta x}
$$
但し,$\Delta x → 0$なる極限を取る.この関係より
$$
\begin{eqnarray}
a(t) &=& \frac{v(t+\Delta t)-v(t)}{\Delta t} \
v(t) &=& \frac{x(t+\Delta t)-x(t)}{\Delta t}
\end{eqnarray}
$$
が成り立つ.変形すると
$$
\begin{eqnarray}
v(t+\Delta t) &=& v(t)+a(t)\Delta t\tag{1}\
x(t+\Delta t) &=& x(t)+v(t)\Delta t\tag{2}\
\end{eqnarray}
$$
となる.この式変形で得られた結果を微分方程式の差分化という.
差分化により,微分している変数のごく小さな変化に関する方程式で近似できるというものである.

自由落下の運動のモデル化

放物運動に拡張する前に最もかんたんなモデルである,空気抵抗を考慮しない自由落下について考えていく.空気抵抗がない自由落下の運動方程式は,鉛直上向きを正にとると
x軸方向
$$
F=0
$$
y軸方向
$$
F=-mg
$$
となり,プログラム化すると

sample1.py

def m_step(x,y,x_vec,y_vec,stride):
    g = 9.8
    y += y_vec*stride
    x += x_vec*stride
    y_vec = y_vec - g*stride

    return [x,y,x_vec,y_vec]

def motion(x,y,x_vec,y_vec):
    stride = 0.3
    cur = [x,y,x_vec,y_vec]
    for i in range(1,100):
        cur = m_step(cur[0],cur[1],cur[2],cur[3],stride) 
        print(cur)
    return cur

motion(10,10,0,0)

となる.初期値(10,10)における自由落下の微分方程式である.ここで$\Delta x=0.3$としている.このプログラムの場合だとyが負である時でも処理を実行していく.つまり実際には地面があるにも関わらず地中深くに物体が潜り続けるということである.実際にはそんなことはありえないため,yに対して条件を設け,更にグラフ描写も行う.

sample2.py

import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
mpl.rcParams['font.family'] = 'Hiragino Sans'

def m_step(x,y,x_vec,y_vec,stride):
    g = 9.8
    y += y_vec*stride
    x += x_vec*stride
    y_vec = y_vec - g*stride

    return [x,y,x_vec,y_vec]

def motion(x,y,x_vec,y_vec):
    stride = 0.01
    cur = [x,y,x_vec,y_vec]
    xset = []
    yset_vec = []
    yset_y = []
    alltime = 0
    while (cur[1] > 0 ):
        cur = m_step(cur[0],cur[1],cur[2],cur[3],stride) 
        alltime += stride
        xset.append(alltime)
        yset_vec.append(abs(cur[3]))
        yset_y.append(cur[1])
    plt.xlabel("時刻:t")
    plt.ylabel("変位:y(t)")
    plt.plot(xset,yset_y)
    plt.show()
    plt.xlabel("時刻:t")
    plt.ylabel("y方向の速度の大きさ:v(t)")
    plt.plot(xset,yset_vec)
    plt.show()
    return cur

motion(10,10,0,0)

計算する際に条件としてstr[1] > 0とすることでyが負にならないようにしている.日本語フォントの設定は,以下の部分で行っている.'AppleGothic'を指定したところ"変"のみが文字化けしたため,'Hiragino Sans'としている.なぜだろう・・・。

sample2-日本語フォントの設定部.py

import matplotlib as mpl
mpl.rcParams['font.family'] = 'Hiragino Sans'

ともかく,得られたグラフが下記の2つである.

自由落下(空気抵抗あり)のモデル化

実際に,自由落下する際には空気抵抗を考慮する必要があり空気抵抗係数を$k$とすると
y軸方向
$$
ma = -mg + kv
$$
となる.これまでは,空気抵抗を考慮しなかったため加速度$a=-g$なり質量$m$は考慮しなくてよかった.しかし,空気抵抗を考えることで,
$$
a = -g + \frac{k}{m}v \tag{3}
$$
となり,質量$m$に応じて加速度も変化する.
(1)式に(3)式を代入すると

$$
\begin{eqnarray}
v(t+\Delta t)&=&v(t)+(-g+\frac{k}{m}v)\Delta t \
&=& v(t)-g\Delta t + \frac{k}{m}v\Delta t
\end{eqnarray}
$$
となる.この関係を用いてコード化すると

import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
mpl.rcParams['font.family'] = 'Hiragino Sans'

def m_step(x,y,x_vec,y_vec,stride):
    g = 9.8
    k = 0.1
    m = 1.0
    y = y + y_vec*stride
    x = x + x_vec*stride
    y_vec = y_vec - g*stride - k/m*y_vec*stride
    return [x,y,x_vec,y_vec]

def motion(x,y,x_vec,y_vec):
    stride = 0.01
    cur = [x,y,x_vec,y_vec]
    xset = []
    yset_vec = []
    yset_y = []
    alltime = 0
    while (cur[1] > 0 ):
        cur = m_step(cur[0],cur[1],cur[2],cur[3],stride) 
        alltime += stride
        xset.append(alltime)
        yset_vec.append(abs(cur[3]))
        yset_y.append(cur[1])
    plt.xlabel("時刻:t")
    plt.ylabel("変位:y(t)")
    plt.plot(xset,yset_y)
    plt.show()
    plt.xlabel("時刻:t")
    plt.ylabel("y方向の速度の大きさ:v(t)")
    plt.plot(xset,yset_vec)
    plt.show()
    return cur

motion(0,300,0,0)

となる.速度$v_y$は

となる.この時k=0.5,m=1.0である.

最後に

つづきを書いていこうかね