はじめに
前回の記事からの続きです。前回の記事はこちら↓↓↓
https://qiita.com/nefu_chem/items/1aeea4d9e58539c23e0d
振り子の周期の厳密解を求める。
では,いよいよ本記事のメインに入っていきます。前回までに運動方程式からの導出では非線形微分方程式が出てしまい,厳密解を求めることが難しいというところを見ました。ですのでここからは方針を変えて,エネルギー保存則を使っていきます。
振り子のエネルギーは
E = \frac{1}{2}mv^2-mgl \cos \theta
とかけます。ただし振り子の位置エネルギーは最下点で$-mgl$となるように原点を選びます。摩擦力や空気抵抗などは無視しているので$E$は一定で,振り子の速さが最大値をとるときに振幅は最大になるので$E=-mgl\cos \theta_m$となります。また,振り子の速度は$v=l \dfrac{d\theta }{dt}$だから
\frac{d\theta }{dt} = \pm \omega_0 \sqrt{2(\cos \theta - \cos \theta_m)}
となります。
$t=0$のとき$\theta = 0$で$\dfrac{d\theta }{dt} > 0 $であるとして,$\theta = \theta_m$になるとき(振り子の振幅が最大になるとき)の時刻を$t=t_m$とします。周期はこの4倍になるはずなので
T = 4\int_0^{t_m}dt = 4\int_0^{\theta_m} (\frac{d\theta }{dt}) ^{-1} d\theta = 4\int_0^{\theta_m} \frac{d\theta}{\omega_0 \sqrt{2(\cos \theta - \cos \theta_m)}}
となります。ここで現れた積分
I = \int_0^{\theta_m} \frac{d\theta}{\sqrt{2(\cos \theta - \cos \theta_m)}}
は,変数を以下のように
\sin \frac{\theta }{2} = \sin \frac{\theta_m}{2}\sin \phi
$\phi $に変換してみると
I = \int_0^{\frac{\pi }{2}}\frac{d\phi }{\sqrt{1-\sin ^2 (\frac{\theta _m}{2})\sin^2 \phi }} = K(\sin \frac{\theta _m}{2})
となります。ここの$K(k)$は第一種完全楕円積分と呼ばれるもので,以下のように定義されます。
K(k)= \int_0^{\frac{\pi }{2}}\frac{d\phi }{\sqrt{1- k^2\sin^2 \phi }} = \int_0^1 \frac{du}{\sqrt{(1-u^2)(1-k^2u^2)}}
ここで$k$は0以上1以下の定数です。これが本記事で初めて出てきた特殊関数です。この定義を用いると,周期は
T = \frac{4}{\omega _0} K(\sin \frac{\theta _m}{2})
で表されます。
求めた周期を最大振幅$\theta_m$の関数として表したグラフを作りましょう。作ったものが以下です...とやりたかったんですがグラフ作成でエラーを吐かれて無理だったので,計算サイト(https://keisan.casio.jp/ )に手を貸してもらいます。このサイトのこのページ(https://keisan.casio.jp/exec/system/1166751749?lang=&charset=utf-8&var_k=1&ketasu=22#! )にあるグラフを見てみます。
このグラフの$k$が今は$\sin \frac{\theta _m}{2}$になっているので$k=1$が$\theta_m = 180 $に対応します。このグラフより最大振幅が大きいときは周期が正の無限大に発散し,小さいときは"ある値"に収束することが分かります。実はこの"ある値"に$\dfrac{4}{\omega_0}$をかけたものが$T_0$になります。これは振幅が小さいことを仮定して近似をしたことを考えれば納得できます。
これで周期については厳密解を求めることができました。でも実はもう少し計算を進めることができます,がそれは最後の記事に続きます。次の記事は振り子の運動についてもう少し深く考えます。
参考
高精度計算サイト https://keisan.casio.jp/