はじめに
微分方程式は数学だけでなく、物理学、化学、生物学など多岐にわたる分野においてとても重要な役割を担っている。例えば物体の単振動について考えてみる。この場合、現在のような特定の時刻における加速度の情報を物体に力より求めることは運動方程式という微分方程式を用いれば極めて単純に計算することができる。一方で時刻が変化したときの物体の運動の様子を運動方程式から推定するのはかなり面倒になる。というのは時間というものは実数のように連続的なものであると仮定しているため、理想上では運動方程式を各時刻ごとに立式し、それらの和を取ることで初めて、位置や速度を求めることができる。なので、微分方程式は立式自体容易だが数学的に解くのは難しい場合がほとんどである。そこで、まず今回は比較的容易にイメージできる単振り子の安定性をエネルギー保存の法則から理論的に推定する。そしてその特性が正しいと言えるのか計算機による数値計算を用いて視覚的に確認する。そのときに安定性判別などで用いられる相平面という概念を簡単に説明する。
単振り子の微分方程式
以下のような、シンプルな単振り子を考える。
この系の質点に対して運動方程式を立式してみる。ただし、動径方向の角度を$\theta$とする。
厳密式
m(l \ddot{\theta})=-mg \sin{\theta}
つまり、
\ddot{\theta}=-\frac{g }{l}\sin{\theta}
構造的には極めてシンプルだが、$\sin\theta$が極めて複雑に作用するため、この微分方程式の解は楕円関数という特殊な関数を用いなければならない。
近似式
そこで、一般的には単振り子の使用範囲は、$|\theta|<<1$つまり単振り子の振れ幅は長さに対して極めて小さいと近似的にみなせるとし、以下の式が成立すると仮定する。
\ddot{\theta}=-\frac{g }{l}\theta
これなら、計算機や大学以上の高度な数学の知識を用いずとも一般解を求めることができる。
今回は、$\omega =\sqrt{\frac{g }{l}}$としたとき、
\theta=A\sin(\omega t+\alpha),\dot{\theta}=A\omega\cos(\omega t+\alpha)
と表すことができる。ただし、$A,\alpha$は定数である。
（成立する理由は、実際に解を微分方程式に代入してみるとよい。）
相平面について
横軸に、物体の変位(今回は$\theta$)を、縦軸にはその時間微分である$\dot{\theta}$を対応させる、
近似式
まず、簡単な近似式の場合から説明する。
三角関数の相互関係の式より、
\theta=A\sin(\omega t+\alpha),\dot{\theta}=A\omega\cos(\omega t+\alpha)
は、
(\frac{\theta}{A})^2+(\frac{\dot{\theta}}{A\omega})^2=1
となる。つまり相平面は楕円になる。
ちなみにこれは、わざわざ微分方程式を解かなくても時間項を削除するエネルギー保存の導出方法を使っても導出することができる。
\ddot{\theta}=-\frac{g }{l}\theta
の両辺に$\frac{1}{2}\dot{\theta}$を掛けて積分をする。
\int \frac{1}{2}\dot{\theta}\ddot{\theta}=-\int \frac{1}{2}\frac{g }{l}\theta\dot{\theta}
\dot{\theta}^{2}=-\frac{g}{l}\theta^2+C
ただし、$C$は定数とする。このように、確かに楕円の方程式を示すようである。
厳密式
それでは、以下に示す厳密式について考えてみる。
\ddot{\theta}=-\frac{g }{l}\sin{\theta}
近似式と同様に、両辺に$\frac{1}{2}\dot{\theta}$を掛け積分する。
\int \frac{1}{2}\dot{\theta}\ddot{\theta}=-\int \frac{1}{2}\frac{g }{l}\sin{\theta}\dot{\theta}
したがって、
\dot{\theta}^2=-\frac{g}{2l}\int \sin \theta d\theta
ここで、初期条件を$\theta=\theta_0,\dot{\theta}=0$とすると、
\dot{\theta}^2=-\frac{g}{2l}\int^{\theta}_{\theta_0} \sin \theta d\theta
\dot{\theta}^2=\frac{g}{l}(\cos \theta -\cos\theta_0)
プログラム
さて、以上の議論をもとに相平面を計算機によて微分方程式を数値計算することによって描写することを試みる。
そこで、以下の記事のプログラムを改造することを考える。
\begin{equation}
\left\{ \,
\begin{aligned}
& \frac{dx}{dt} =y \\
& \frac{dy}{dt} =f(x)\\
\end{aligned}
\right.
\end{equation}
と表される連立微分方程式は、$t$を消去することで以下のように表すことができる。
\frac{dy}{dx}=\frac{f(x)}{y}
したがって、それらを反映させたプログラムを以下に示す。
ただし、厳密式と近似式2つ共用なので、上手くコメント文を切り分けてほしい。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
# ドット数
n=1000
#xy座標系における描写範囲
L1=5*math.pi
L2=10
# x軸の刻み幅
dx= 2*L1/n
g=9.8
l=1
x=np.linspace(-L1 ,L1, n)
y=np.linspace(-L2, L2, n)
X, Y = np.meshgrid(x, y)
U = np.zeros((n, n))
V = np.zeros((n, n))
for p in range(n):
for q in range(n):
xx=X[p][q]
yy=Y[p][q]
##微分方程式
##一般式の場合
f_y=-(g/l)*math.sin(xx)/yy
##近似式の場合
#f_y=-(g/l)*xx/yy
U[p][q]=dx
V[p][q]=f_y*dx
plt.streamplot(X, Y, U, V, color='blue', density=1.5, linewidth=0.5, arrowsize=1.5)
plt.xlabel("変位")
plt.ylabel("速度")
plt.savefig('単振り子の相平面_厳密式.png')
#plt.savefig('単振り子の相平面_近似式.png')
plt.show()
これを実行すると以下のようになる。
結果
近似式
近似式での相平面は以下のようになった。
たしかに、理論的な考察と同様に楕円形を示すことが分かった。
厳密式
厳密式での相平面は以下のようになった。
一方で、厳密式であるがこのように、$\theta=\pi$などの、最高点に質点が存在するとき、質点の変位は、少しずれただけで速度が正の方向もしくは負の方向に急変する。つまり最高点に質点が存在しているときが単振り子は不安定であるといえる。
まとめ
今回は、単振り子という基本的な系においての安定性を議論するために相平面という概念を用いた。また、微分方程式が複雑で解くのが難しい場合でもエネルギー保存則の導出のように積分を計算を行うことで時間成分を消去し、相平面の方程式自体は理論的に求めることができた。最後に、計算機を用いた数値計算で相平面を実際に描写することで、単振り子の安定性についても議論した。
参考文献