[Pythonで物理シミュレーション] ばね振動の運動方程式をルンゲクッタ法で解いて可視化。連立微分方程式

2021/12/9　投稿
#0 メニュー
1 問題設定
2 2階微分方程式の数値解法
3 計算コード
4 位置、速度の時間変化: xt、vtグラフ
5 運動アニメーション
6 位置と速度の関係 xv空間中の解軌跡
7 減衰振動

※数式は読み飛ばしてグラフ、アニメーションを先に見ても構いません

#1 問題設定

ばねから復元力$-kx$を受けている状況を考える
ばねの自然長からの伸びを$x(t)$とする
この時の運動方程式は

\frac{d^2 x(t)}{dt^2} = -\omega^2 x(t)
\space\space\space\space
\omega = \sqrt{\frac{k}{m}}　\tag{1.1}

$m:$おもりの質量、$k:$ばね定数

#2 2階微分方程式の数値解法
4次のルンゲクッタ法を用いる

2階微分方程式を1階連立微分方程式に書き直す

式$(1.1)$のような2階微分方程式にルンゲクッタ法を適用するには1階の連立微分方程式に書き直す必要がある。
したがって以下のように書き換える

f_v(x(t)) := -\omega^2 x(t) \space\space\space\space
f_x(v(t)) := v(t)
\tag{2.1, 2.2}

\frac{d v(t)}{dt} = f_v(x(t))  \space\space\space\space
\frac{d x(t)}{dt} = f_x(v(t)) 
\tag{2.3, 2.4}

###連立微分方程式に対するルンゲクッタ法
上の式$(2.1)、(2.2)$にルンゲクッタ法を適用すると

k_{x1} = f_x(v(t)) \space\space\space\space
k_{v1} = f_v(x(t))   \tag{2.5, 2.6}

k_{x2} = f_x(v(t) + \frac{\Delta t}{2} k_{v1}) \space\space\space\space
k_{v2} = f_v(x(t) + \frac{\Delta t}{2} k_{x1})   \tag{2.7, 2.8}

k_{x3} = f_x(v(t) + \frac{\Delta t}{2} k_{v2}) \space\space\space\space
k_{v3} = f_v(x(t) + \frac{\Delta t}{2} k_{x2})   \tag{2.9, 2.10}

k_{x4} = f_x(v(t) + \Delta t k_{v3}) \space\space\space\space
k_{v4} = f_v(x(t) + \Delta t k_{x3})   \tag{2.11, 2.12}

####漸化式

x(t+\Delta t) = x(t) + \frac{\Delta t}{6} (k_{x1}+2k_{x2}+2k_{x3}+k_{x4}) \tag{2.13}

v(t+\Delta t) = v(t) + \frac{\Delta t}{6} (k_{v1}+2k_{v2}+2k_{v3}+k_{v4}) \tag{2.14}

###ベクトル形式にまとめる
　これをもとにプログラムを書けばいいが、かなり式が長くなってしまった。式$(2.3)$から式$(2.14)$をベクトル形式にまとめる
以下のものを定義する

\textbf{x}(t) := \begin{pmatrix}
x(t) \\
v(t)
\end{pmatrix} 
\tag{2.15}

連立微分方程式、式$(2.1)$から式$(2.4)$をベクトル形式にすると
連立微分方程式の右辺

\textbf{f}(\textbf{x}(t)) := \begin{pmatrix}
f_x(v) \\
f_v(x)
\end{pmatrix} 

= \begin{pmatrix}
v(t) \\
-\omega^2 x(t) 
\end{pmatrix} 

\tag{2.16}

\frac{d\textbf{x}(t)}{dt} = \textbf{f}(\textbf{x}(t))
\tag{2.17}

###ベクトル形式のルンゲクッタ法
式$(2.17)$にルンゲクッタ法を適応するには式$(2.5)$から式$(2.14)$をベクトル形式に書き換えたものを使う。

\textbf{k}_{i} := \begin{pmatrix}
k_{xi} \\
k_{vi}
\end{pmatrix}
\tag{2.18}

\textbf{k}_{1} = \textbf{f}(\textbf{x}(t)) 
\tag{2.19}

\textbf{k}_{2} = \textbf{f}(\textbf{x}(t) + \frac{\Delta t}{2} \textbf{k}_{1}) 
\tag{2.20}

\textbf{k}_{3} = \textbf{f}(\textbf{x}(t) + \frac{\Delta t}{2} \textbf{k}_{2}) 
\tag{2.21}

\textbf{k}_{4} = \textbf{f}(\textbf{x}(t) + \Delta t \textbf{k}_{3}) 
\tag{2.22}

\textbf{x}(t+\Delta t) = \textbf{x}(t) + \frac{\Delta t}{6} (\textbf{k}_{1}+2\textbf{k}_{2}+2\textbf{k}_{3}+\textbf{k}_{4}) \tag{2.19}

　コンピューターで複数本の方程式を扱う際はベクトル、行列形式のほうが配列にまとめられるから扱いやすい。今後、連立微分方程式と数値的に解くときはこのやり方で計算する。なおこのやり方は一般的に未知関数の数が増えても形式的に同じように書ける。

#3 計算コード

import numpy as np

# 物理定数
m_pa = 1.0 
k_pa = 1.0
omega_pa = (k_pa/m_pa)**0.5

# 数値解用変数
x_vec = np.zeros(2)

# 曲線plot用
x_line = []
v_line = []

# 初期条件
xmax = 10.0
x_vec[0] = xmax
x_vec[1] = -1e-6

# 時間変数
tmin = 0.0 
tmax = 10.0 
dt = 0.001
t = tmin

# 連立微分方程式の右辺
def f(omega, x_vec):
    x, v = x_vec[0], x_vec[1]
    re = np.array([
        v,
        -omega**2 * x - 0.3*v
    ])
    return re 


# アニメーション用データカット変数
sb_dt = 0.1
Nt = 1
count = 1
cat_val = int(sb_dt/dt)
x_lis = [x_vec.copy()]

# 曲線plot用
x_line.append(x_vec[0])
v_line.append(x_vec[1])

# ルンゲクッタ変数
beta = [0.0, 1.0, 2.0, 2.0, 1.0]
delta = [0.0, 0.0, dt/2, dt/2, dt]
k_rk = np.zeros((5, 2))

# 時間積分
while t < tmax:
    sum_x = np.zeros(2)
    for k in [1, 2, 3, 4]:
        k_rk[k, :] = f(omega_pa, x_vec + delta[k]*k_rk[k-1, :])
        sum_x += beta[k]*k_rk[k, :]
    #
    # 更新
    x_vec += (dt/6)*sum_x
    # 曲線plot用
    x_line.append(x_vec[0])
    v_line.append(x_vec[1])
    #
    # アニメーション用データ
    if count%cat_val == 0:
        x_lis.append(x_vec.copy())
        Nt += 1
    #
    count += 1
    t += dt

t_line = np.linspace(tmin, tmax, len(x_line))

#4 位置、速度の時間変化: xt、vtグラフ

###コード

import matplotlib.pyplot as plt
plt.xlabel('t[s]')
plt.plot(t_line, x_line, color='blue')
plt.plot(t_line, v_line, color='red')
plt.show()

###グラフ

青い曲線が位置$x(t)$の変化で、赤い曲線が速度$v(t)$の変化(横軸: 時刻$t$)

#5 運動アニメーション
###コード

import matplotlib.animation as animation

fig, ax = plt.subplots()

def animate_move(i):
    plt.cla()
    plt.xlabel('x')
    plt.xlim(-1.2*xmax, 1.2*xmax)
    plt.ylim(-1.0, 1.0)
    # 変数
    v = x_lis[i][1]
    x = x_lis[i][0]
    a = f(omega_pa, x_lis[i])[1]
    # 質点
    plt.scatter(x, 0.0, s=1000)
    # 速度
    plt.quiver(x, 0.2, v, 0.0, color='green', scale=50.0)
    # 加速度
    plt.quiver(x, 0.0, a, 0.0, color='red', scale=50.0)

animate = animate_move
kaisu = Nt - 1
anim = animation.FuncAnimation(fig,animate,frames=kaisu,interval=1000)
anim.save("anime_bool.gif", writer="imagemagick")
plt.show()

###アニメーション

赤い矢印は力$f=-kx(t)$で、緑の矢印は速度$v(t)$を表している
$x$:ばねの自然長からの伸び

#6 位置と速度の関係 xv空間中の解軌跡
###コード

fig, ax = plt.subplots()

vmin = min(v_line)
vmax = max(v_line)
def animate_xv(i):
    plt.cla()
    plt.xlabel('x')
    plt.ylabel('v')
    plt.xlim(-1.2*xmax, 1.2*xmax)
    plt.ylim(1.2*vmin, 1.2*vmax)
    # 変数
    v = x_lis[i][1]
    x = x_lis[i][0]
    plt.plot(x, v, 'bo')
    # ベクトル
    plt.quiver(x, v, v, -omega_pa**2 * x, scale=50.0)
    plt.plot(x_line, v_line)
    plt.grid()


animate = animate_xv
kaisu = Nt - 1
anim = animation.FuncAnimation(fig,animate,frames=kaisu,interval=1)
anim.save("anime_xv.gif", writer="imagemagick")
plt.show()

###アニメーション

青い曲線が解軌跡$\textbf{x} = (x,v)$
黒いベクトルが「原点から軌跡までの位置ベクトル」の微分値$\textbf{f}=(v,-\omega^2 x)$
軌跡のパラメーター$t$微分がの接線ベクトルになっている。

#7: 減衰振動
復元力に加えて空気抵抗$-\gamma v$が加わる場合でも計算してみた。この時の連立微分方程式の右辺, 式$(2.16)$は以下のように書き換えられる

\textbf{f}(\textbf{x}(t)) := \begin{pmatrix}
f_x(v) \\
f_v(x)
\end{pmatrix} 

= \begin{pmatrix}
v(t) \\
-\omega^2 x(t) - \gamma v(t)
\end{pmatrix} 

\tag{7.1}

コードを書き換える

# 連立微分方程式の右辺
def f(omega, x_vec):
    x, v = x_vec[0], x_vec[1]
    gamma = 0.3
    re = np.array([
        v,
        -omega**2 * x - gamma*v
    ])
    return re 

位置、速度の時間変化: xt、vtグラフ

青い曲線が位置$x(t)$の変化で、赤い曲線が速度$v(t)$の変化(横軸: 時刻$t$)

運動アニメーション

赤い矢印は力$f=-kx(t)-\gamma v(t)$で、緑の矢印は速度$v(t)$を表している
$x$:ばねの自然長からの伸び

位置と速度の関係 xv空間中の解軌跡

青い曲線が解軌跡$\textbf{x} = (x,v)$
黒いベクトルが「原点から軌跡までの位置ベクトル」の微分値$\textbf{f}=(v,-\omega^2 x - \gamma v)$
軌跡のパラメーター$t$微分がの接線ベクトルになっている。

End