Return to コンピュータ演習

単振り子のエネルギー保存則から周期を求める準備

単振り子の運動方程式から得られるエネルギー保存則から,数値積分によって周期を求めるための準備。

問題(3択)

まずは予想を立ててみよう。

単振り子の振幅が大きくなると(「十分に小さい」という近似が成り立たなくなると),単振り子の周期は,単振動の場合の周期にくらべて…

  1. 「振り子の等時性」があるから変化しない,一定のままである。
  2. 振幅が大きくなればなるほど,周期は長くなる。
  3. 振幅が大きくなればなるほど,周期は短くなる。

 

保存量(エネルギー保存則)

無次元化された運動方程式

$$ \frac{d^2\theta}{dT^2} = – 4\pi^2 \sin\theta$$

の両辺に $\displaystyle \frac{d\theta}{dT}$ をかけて整理すると,

\begin{eqnarray}
\frac{d\theta}{dT}\frac{d^2\theta}{dT^2}+4\pi^2 \sin\theta  \frac{d\theta}{dT}&=& 0 \\
\frac{d}{dT}\left\{
\frac{1}{2}\left( \frac{d\theta}{dT} \right)^2 – 4\pi^2 \cos\theta
\right\} &=& 0 \\
\therefore\ \ \frac{1}{2}\left( \frac{d\theta}{dT} \right)^2 – 4\pi^2 \cos\theta
&=& \mbox{const.} \equiv \varepsilon
\end{eqnarray}

つまり,

$$
\varepsilon = \frac{1}{2}\left( \frac{d\theta}{dT} \right)^2 – 4\pi^2 \cos\theta
$$

という保存量が存在することがわかる。これはエネルギー保存則に相当する。

$T = 0$ で $\displaystyle \theta = \theta_0, \frac{d\theta}{dT} = 0$ という初期条件を課すと,定数 $\varepsilon$ の値は以下のように定まる。
$$ \varepsilon = – 4\pi^2 \cos\theta_0$$

つまり,$\theta = \theta_0$ から $\theta = 0$ まで $\frac{d\theta}{dT} < 0$ の状況を考えることにして

\begin{eqnarray}
\frac{1}{2}\left( \frac{d\theta}{dT} \right)^2 &=& 4\pi^2 \left( \cos\theta-\cos\theta_0\right)\\
\frac{d\theta}{dT} &=& – 2\pi \sqrt{2(\cos\theta-\cos\theta_0)} \\
-\frac{1}{2\pi} \frac{d\theta}{\sqrt{2(\cos\theta-\cos\theta_0)}} &=& dT
\end{eqnarray}

これは変数分離形の微分方程式であり,$\theta$ が $\theta_0$ から $0$ までの時間の4倍が周期 $T_p$ になることから,

\begin{eqnarray}
-\frac{1}{2\pi} \int^0_{\theta_0} \frac{d\theta}{\sqrt{2(\cos\theta-\cos\theta_0)}} &=& \int_0^{\frac{T_p}{4}} dT
\end{eqnarray}

すなわち

$$
T_p = \frac{2}{\pi} \int_0^{\theta_0} \frac{1}{\sqrt{2(\cos\theta-\cos\theta_0)}}d\theta
$$

被積分関数が発散すると困る?

この積分は解析的には解けないので数値積分すればいいと思いますが,このままの形だと,積分区間の上限で被積分関数の分母がゼロになってしまい,つまりは被積分関数が発散してしまうので,大変なことがおこりそう…

Maxima の quad_qags() や Python の scipy.integrate.quad() なら無問題

… などと,いろいろ心配しましたが,被積分関数が発散するからといって積分そのものが発散するわけではありません。数値積分をする関数である Maxima の quad_qags()romberg() はダメ)や,Python の scipy.integrate.quad() は,$\theta = \theta_0$ で被積分関数が発散してしまうような場合でも,何事もなかったように数値積分してしまいます。

第1種完全楕円積分

ここでは,上記の積分が変数変換によって,第1種完全楕円積分 $K(k)$ を使って
\begin{eqnarray}
\int_0^{\theta_0} \frac{1}{\sqrt{2(\cos\theta-\cos\theta_0)}} d\theta
&=& K(k) \\&\equiv&
\int_0^{\pi/2} \frac{dt}{\sqrt{1 – k^2 \sin^2 t}}, \quad k = \sin\frac{\theta_0}{2}
\end{eqnarray}

で書けることを示す。この形だと,$k^2 < 1$ なら積分区間内で被積分関数が発散することもなく,数値積分に困ることもない。

まず,分母は半角の公式を使って

\begin{eqnarray}
\sqrt{2 (\cos\theta-\cos\theta_0)} &=&
\sqrt{2 \left\{ \left(1-2\sin^2\frac{\theta}{2}\right)
-\left(1-2\sin^2\frac{\theta_0}{2}\right)
\right\} }\\
&=& 2\sqrt{\sin^2\frac{\theta_0}{2} – \sin^2\frac{\theta}{2}}
\end{eqnarray}

ここで,

$$\displaystyle \sin\frac{\theta_0}{2} \equiv k, \quad
\sin\frac{\theta}{2} \equiv k \sin t
$$

とおくと,

\begin{eqnarray}
2\sqrt{\sin^2\frac{\theta_0}{2} – \sin^2\frac{\theta}{2}}
&=& 2 k \sqrt{1 – \sin^2 t} \\
&=& 2 k \cos t
\end{eqnarray}

一方,分子は

$$ \sin\frac{\theta}{2} \equiv k \sin t $$

の微分をとることにより,

$$ \cos\frac{\theta}{2} \frac{d\theta}{2} = k \cos t \, dt$$

$$
\therefore \ d\theta = \frac{2 k \cos t}{\cos\frac{\theta}{2}} dt
= \frac{2 k \cos t}{\sqrt{1-\sin^2\frac{\theta}{2}}} dt
= \frac{2 k \cos t}{\sqrt{1-k^2 \sin^2 t}} dt
$$

したがって,

\begin{eqnarray}
\int_0^{\theta_0} \frac{1}{\sqrt{2(\cos\theta-\cos\theta_0)}} d\theta
&=&\int_0^{\pi/2}
\frac{1}{2 k \cos t} \frac{2 k \cos t}{\sqrt{1-k^2 \sin^2 t}} dt\\
&=&\int_0^{\pi/2} \frac{dt}{\sqrt{1-k^2 \sin^2 t}}\\
&\equiv& K(k)
\end{eqnarray}

第1種完全楕円積分 $K(k)$ を使うと,振幅 $\theta_0$ の単振り子の周期 $T_p$ は

\begin{eqnarray}
T_p(\theta_0) &=&
\frac{2}{\pi} \int_{0}^{\theta_0} \frac{1}{\sqrt{2(\cos\theta-\cos\theta_0)}}d\theta \\
&=& \frac{2}{\pi} K(k) \\
&=& \frac{2}{\pi}\int_0^{\pi/2} \frac{dt}{\sqrt{1-k^2 \sin^2 t}}, \quad k \equiv \sin\frac{\theta_0}{2}
\end{eqnarray}

SciPy で単振り子のエネルギー保存則から周期を求める

単振り子の運動方程式から得られるエネルギー保存則から,Python の SciPy と NumPy を使って(かつ SymPy は使わずに)数値積分によって周期を求める。導出については以下を参照:

  • 単振り子のエネルギー保存則から周期を求める準備

SymPy で単振り子のエネルギー保存則から周期を求める

単振り子の運動方程式から得られるエネルギー保存則から,Python の SymPyを使って(かつ SciPy と NumPy は使わずに)数値積分によって周期を求める。導出については以下を参照:

  • 単振り子のエネルギー保存則から周期を求める準備

Maxima で単振り子のエネルギー保存則から周期を求める

単振り子の運動方程式から得られるエネルギー保存則から,Maxima を使って数値積分によって周期を求める。導出については以下を参照:

  • 単振り子のエネルギー保存則から周期を求める準備