減衰振動【Pythonと物理#4】

今回のターゲット：減衰振動
減衰振動
Pythonで計算する
実行結果
最後に

今回のターゲット：減衰振動

　前回は調和振動子（単振動）の微分方程式を計算してみました。

今回のターゲット：調和振動子（単振動）の運動　前回は物体の放物運動をアニメーションにしてみました。今回は調和振動子（単振動）の微分方程式を解いてみます。調和振動子の微分方程式調和振動子とは、ある点からの距離 $x$ 比例した力を受けて運動す...

今回は減衰振動の微分方程式を計算してみます。

減衰振動

減衰振動は前回の調和振動子（単振動）の運動に速度に比例する抵抗力が加わる形になります。微分方程式は以下のようになります。

$$
\begin{eqnarray}
\begin{array}{l}
m\cfrac{d^{2}x}{dt^{2}} = -kx-c\cfrac{dx}{dt}\\
\end{array}
\end{eqnarray}
$$

調和振動子のときと同様にして微分方程式を１階の微分方程式に書き直します。

$$
\begin{eqnarray}
\begin{array}{l}
m\cfrac{dv}{dt} = -kx-cv
\end{array}
\end{eqnarray}
$$

Pythonで計算する

コードは下のようにしました。前回とほぼ変わりません。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint 
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML

#諸々の条件設定
m = 1.0 #物体の重さを m = 1.0 kgとしておきます
k = 100 #ばね定数は k = 100 N/mとしておきます
c = 1 #減衰係数は、c = 1 N・s/m
ini = [0.1, 0.0]#バネを 0.1 m伸ばした状態をt=0として、初速v0 = 0.0 m/sとしておきます
t = np.arange(0, 10, 0.001) #時間は2.0秒程度に設定しておきます

#計算に必要な関数(微分方程式)の定義
#vとdvdtの計算（戻り値は速度vと加速度dvdt）
def func(cond, t, m, k, c):
    x, v = cond
    dvdt = (-k/m) * x - (c/m) * v
    return [v, dvdt]

#x座標とy座標の計算　1階微分方程式を解く
#注意：odeintは微分方程式(おそらく1階微分方程式のみ?)を解くライブラリ
#引数は(解きたい微分方程式, 初期値, 変数, args = (パラメータ1, パラメータ2))
#なお、パラメータはタプルで渡して、1つしかない場合は、args = (パラメータ1, )とする
condition = odeint(func, ini, t, args = (m, k, c))

#以下は可視化のためにグラフを書くコード
fig, ax = plt.subplots()
#今回はx座標の時間変化を見たいので、x[:, 0]とする
#vを見たい場合はx[:, 0]とする
ax.plot(t, condition[:,0], color='crimson', label='x')

ax.set_xlabel('t [sec]')
ax.set_ylabel('x [m]')

ax.tick_params(direction='in', axis='x')
ax.tick_params(direction='in', axis='y')

ax.set_xlim(0, 10)
#ax.set_ylim(-0.001, 0.001)

plt.show()