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

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

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

quad_qags() で数値積分する

振幅 $\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{\sqrt{2}}{\pi} \int_{0}^{\theta_0} \frac{1}{\sqrt{\cos\theta-\cos\theta_0}}d\theta
\end{eqnarray}

積分部分は数値積分 quad_qags() を使う。（romberg() ではなく quad_qags() を使う！）

この積分が一見やっかいなのは，積分範囲の上限 $\theta = \theta_0$ で被積分関数の分母がゼロになることです。$\theta = \theta_0$ で被積分関数が発散してしまうので，何かエラーでも起こるのかと思いきや，quad_qags() は大変優秀で何事もなかったかのように数値積分してしまいます！

ためしに $\theta_0 = 60^{\circ}$ の場合の積分値を数値積分 quad_qags() で求めてみます。

quad_qags() は４つの要素のリストを返します:

積分の近似値
見積もられた近似の絶対誤差
被積分関数の評価数
エラーコード

quad_qags() で Tp1(th0)

積分を使って定義された周期を Tp1(th0) と定義する。

置換積分で被積分関数が発散しないようにする

このへんにまとめているように，変数変換によって
\begin{eqnarray}
\int_0^{\theta_0} \frac{1}{\sqrt{2(\cos\theta-\cos\theta_0)}} d\theta
&=&
\int_0^{\pi/2} \frac{dt}{\sqrt{1 – k^2 \sin^2 t}}, \\ \quad k &\equiv& \sin\frac{\theta_0}{2}
\end{eqnarray}

と書ける。この形にすると，積分範囲内で被積分関数が発散することもない。

quad_qags() で Tp2(th0)

置換積分を使って定義された周期を Tp2(th0) と定義する。

第1種完全楕円積分を使う

このへんにまとめているように，
第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}

elliptic_kc(m) で Tp3(th0)

ありがたいことに，Maxima では以下のように定義された第一種完全楕円積分が elliptic_kc(m) として使えます。

$$\mbox{elliptic_kc(m)} = K(m) \equiv \int_0^{\frac{\pi}{2}}\frac{1}{\sqrt{1 – m \sin^2 x}} dx$$

第1種完全楕円積分 elliptic_kc(m) を使って，Tp3(th0) を以下のように定義する。（$m = k^2$ に注意。）

Tp3(th0) のグラフ

練習

Tp1(th0)，Tp2(th0)，Tp3(th0) の誤差を調べる。