Pythonで常微分方程式を解く(ルンゲクッタ法、オイラー法、解析的手法)

今日はPython常微分方程式を解いてみます。
空気抵抗がある状態でのボールの自由落下における速度と軌跡を求めましょう。

数値計算

数値計算では4次のルンゲクッタ法を使用して計算します。
ルンゲクッタ法は数値解析において微分方程式の初期値問題に対して比較的良い近似解を与える方法ですが、それに関しては詳しくはWikipediaなどをご覧ください。

import numpy as np
import matplotlib.pyplot as plt

g = 9.81 #重力加速度(m/s^2)
k = 0.24 #空気抵抗係数(kg/m)
m = 10.0 #ボール重さ(kg)

def f(v): #関数の定義
    return (m*g-k*v**2)/m

vtrj = [0.0] #速度のtrajectory配列
xtrj = [100.0] #位置のtrajectory配列
array_t = [0.0] #時間発展配列

v = vtrj[0] #初期速度
x = xtrj[0] #初期位置
t = array_t[0] #初期時間
h = 0.1e0 #ルンゲクッタ法の刻み幅

#ルンゲクッタ法のメイン計算
while x >= 0.0e0:
    k1 = h*f(v)
    k2 = h*f(v+0.5*k1)
    k3 = h*f(v+0.5*k2)
    k4 = h*f(v+k3)
    v = v + (k1+2*k2+2*k3+k4)/6
    x = x - v*h
    t = t + h
    vtrj = np.append(vtrj,v)
    xtrj = np.append(xtrj,x)
    array_t = np.append(array_t,t)

#グラフ描写
plt.plot(array_t,vtrj)
plt.plot(array_t,xtrj)

f:id:mochirixi:20180326164745p:plain 終端速度が20.2m/s、100m地点から落とすと6.4秒後に地面に到達することが求められた。
位置の時間発展にはオイラー法を使用している。(ただしこの場合はルンゲクッタ法でも同じ結果が得られる。)
オイラー法の詳しい解説はWikipediaなどを参照してください。

解析的手法

Sympyを使用して解析的に常微分方程式を解くことも可能である。

import numpy as np
import matplotlib.pyplot as plt
from sympy import *

v = Function('v')
t = Symbol('t')

g = 9.81 #重力加速度(m/s^2)
k = 0.24 #空気抵抗係数(kg/m)
m = 10.0 #ボール重さ(kg)

eq1 = dsolve(v(t).diff(t,1)-(m*g-k*v(t)**2)/m)
print(eq1)

dsolveを使用して微分方程式を解いている。

Eq(t + 1.03045701422812*log(1.0*v(t) - 20.2175666191557) - 1.03045701422812*log(1.0*v(t) + 20.2175666191557), C1)

しかし解析解は上のように直感的に理解しにくいのと、場合によっては求められないときもあるので基本的には数値的解法をするのが良いと思う。

今日はここまで。