調和振動子（単振動）【Pythonと物理#3】

今回のターゲット：調和振動子（単振動）の運動
調和振動子の微分方程式
Pythonで計算する
実行結果
最後に

今回のターゲット：調和振動子（単振動）の運動

　前回は物体の放物運動をアニメーションにしてみました。

物体の放物運動とアニメーション【Pythonと物理#2】

本記事はPythonでアニメーションを作成する方法を紹介しています。

今回は調和振動子（単振動）の微分方程式を解いてみます。

調和振動子の微分方程式

調和振動子とは、ある点からの距離 $x$ 比例した力を受けて運動する質点のことを言います。
２つの意味で古典的な調和振動子は、理想的な”ばね”に接続された質点の運動＝単振動があります。

今回は、下のようにばね定数$k$のばねに接続された、質量$m$の物体について考えます。
今、物体を $x_{0}$ mほど$x$軸の正方向に引っ張り、手を離すことを考えます。ただし、ばねが自然長の時の物体の位置を$x = 0.0$とします。

このとき、物体の運動方程式は以下のようになります。

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

この微分方程式は２階の常微分方程式です。しかし、pythonのodeintライブラリでは１階の微分方程式しか計算できません。そこで、以下のように１階の常微分方程式に書き直す必要があります。ただし、$v = dx/dt$です。

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

ここまで来たら、この微分方程式をPythonで計算をしてみます。

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としておきます
ini = [0.1, 0.0]#バネを 0.1 m伸ばした状態をt=0として、初速v0 = 0.0 m/sとしておきます
t = np.arange(0, 2, 0.001) #時間は2.0秒程度に設定しておきます

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

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

#以下は可視化のためにグラフを書くコード
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, 2)
#ax.set_ylim(-0.001, 0.001)

plt.show()