Python で連立微分方程式を解き描画する方法

python で微分方程式を特には

from scipy.integrate import odeint

を使用する．

以下では下記コマンドで適宜ライブラリを追加してほしい．

pip install scipy
pip install numpy
pip install matplotlib japanize-matplotlib

微分方程式

微分方程式を解く．

例として自由落下運動を使用する．

重力 $g$ と空気抵抗 $k$ が働くとすると，運動方程式は以下のようになる．

m \frac{dv}{dt} = mg - kv \\　\\ \therefore \frac{dv}{dt} = \frac{kv}{m} - g

これを $t = 0$ の時に $v = 0$ として微分方程式を python で解いて描画する．

以下コードを示す．

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import japanize_matplotlib

# 微分方程式の定義
def func(v, t, m, k, g):
    dvdt = - k * ( v - m * g / k ) / m
    return dvdt


# 質量と係数と重力加速度を指定
m = 100.0
k = 10.0
g= 9.8

# 速度 v の初期値を指定
v = 0.0

# t を 0 から 10 で 1000 個取る
t = np.linspace(0.0, 200.0, 1000)

# 連立常微分方程式をとく
var_list = odeint(func, v, t, args=(m, k, g))

# 以下描画
fig = plt.figure(figsize=(6.4 * 1, 4.8))

ax1 = fig.add_subplot(1, 1, 1, xlabel='$t$', ylabel='速度')

ax1.plot(t, var_list[:, 0], label="$v$")
ax1.legend()
ax1.grid()

fig.savefig(f'AirResistance.png')

終端速度は $\frac{mg}{k} = 98$ となるのでうまく描画されている．

参考

https://manabitimes.jp/physics/1931

連立微分方程式

続いて連立微分方程式を解く．

例として振動現象を上げる．

$u$ と v の濃度が

\epsilon \frac{du}{dt} = -u (u-a) (u-1) - v + I \\\\ \frac{dv}{dt} = bu - cv

で与えられる．

初期値を $t = 0$ のとき $(u, v) = (0, 0)$ として描画する．

以下コードを示す．

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import japanize_matplotlib

def func(var, t, epsilon, a, b, c, I):
    dudt = ( - var[0] * ( var[0] - a ) * ( var[0] - 1 ) - var[1] + I ) / epsilon
    dvdt = b * var[0] - c * var[1]
    return [dudt, dvdt]

epsilon = 0.05
a = 2.5
b = 2
c = 1

I = 2

# t を 0 から 10 で 1000 個取る
t = np.linspace(0.0, 10.0, 1000)

# u,v の初期値を 0 とする
var_init = [0, 0]

# 連立常微分方程式をとく
var_list = odeint(func, var_init, t, args=(epsilon, a, b, c, I))

# 以下描画
fig = plt.figure(figsize=(6.4 * 2, 4.8))

ax1 = fig.add_subplot(1, 2, 1, xlabel='$t$', ylabel='濃度')
ax2 = fig.add_subplot(1, 2, 2, xlabel='$u$', ylabel='$v$')

ax1.plot(t, var_list[:, 0], label="$u$")
ax1.plot(t, var_list[:, 1], label="$v$")
ax1.legend()
ax1.grid()

ax2.plot(var_list[:, 0], var_list[:, 1], label="")
ax2.grid()

fig.savefig(f'VibrationPhenomena.png')