単振り子の運動方程式から得られるエネルギー保存則から，数値積分によって周期を求めるための準備。
保存量（エネルギー保存則）
無次元化された運動方程式
$$ \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
$$
この積分は解析的には解けないので数値積分すればいいと思いますが，このままの形だと，積分区間の上限で被積分関数の分母がゼロになってしまい，つまりは被積分関数が発散してしまうのでちょっと困る。
この積分が「第1種完全楕円積分」となることは別途説明することにして，ここではなんとか数値積分が可能な形にします。
被積分関数が発散する付近でのテイラー展開
この積分が一見やっかいなのは，積分の上限 $\theta = \theta_0$ で被積分関数の分母がゼロとなることです。
$$\lim_{\theta \rightarrow \theta_0-0} \sqrt{2(\cos\theta-\cos\theta_0)} = 0$$
そこで，$\cos\theta-\cos\theta_0$ を $\theta = \theta_0$ のまわりでテイラー展開をすると
\begin{eqnarray}
\cos\theta-\cos\theta_0 &\simeq& (\theta_0 – \theta) \sin\theta_0
– \frac{\cos\theta_0}{2}(\theta_0 – \theta)^2 \\
&=& (\theta_0 – \theta) \sin\theta_0 \left\{1-\frac{\cos\theta_0}{2\sin\theta_0} (\theta_0 – \theta) \right\}
\end{eqnarray}
ということは
\begin{eqnarray}
\frac{1}{\sqrt{2(\cos\theta-\cos\theta_0)}} &\simeq&
\frac{1}{\sqrt{2 (\theta_0 – \theta) \sin\theta_0}}\left\{1-\frac{\cos\theta_0}{2\sin\theta_0} (\theta_0 – \theta) \right\}^{-\frac{1}{2}} \\
&=& \frac{1}{\sqrt{2\sin\theta_0}}\left\{(\theta_0 – \theta)^{-\frac{1}{2}}+\frac{\cos\theta_0}{4\sin\theta_0} (\theta_0 – \theta)^{\frac{1}{2}} \right\}
\end{eqnarray}
そこで，このテイラー展開の表現を使うと
\begin{eqnarray}
\int_{\theta_0-\epsilon}^{\theta_0} \frac{1}{\sqrt{2(\cos\theta-\cos\theta_0)}}d\theta &\simeq&
\frac{1}{\sqrt{2\sin\theta_0}}\int_{\theta_0-\epsilon}^{\theta_0}
\left\{(\theta_0 – \theta)^{-\frac{1}{2}}+\frac{\cos\theta_0}{4\sin\theta_0} (\theta_0 – \theta)^{\frac{1}{2}}
\right\}d\theta \\
&=& \frac{1}{\sqrt{2\sin\theta_0}}
\Biggl[-2(\theta_0 – \theta)^{\frac{1}{2}}
-\frac{\cos\theta_0}{6\sin\theta_0} (\theta_0 – \theta)^{\frac{3}{2}}
\Biggr]_{\theta_0-\epsilon}^{\theta_0}\\
&=&
\frac{1}{\sqrt{2\sin\theta_0}}
\left\{ 2\sqrt{\epsilon} + \frac{\cos\theta_0}{6 \sin\theta_0} \epsilon\sqrt{\epsilon}
\right\}
\end{eqnarray}
\begin{eqnarray}
\int_{0}^{\theta_0} \frac{1}{\sqrt{2(\cos\theta-\cos\theta_0)}}d\theta &=&
\int_{0}^{\theta_0-\epsilon} \frac{1}{\sqrt{2(\cos\theta-\cos\theta_0)}}d\theta \\ && +
\int_{\theta_0-\epsilon}^{\theta_0} \frac{1}{\sqrt{2(\cos\theta-\cos\theta_0)}}d\theta \\
&\simeq&
\int_{0}^{\theta_0-\epsilon} \frac{1}{\sqrt{2(\cos\theta-\cos\theta_0)}}d\theta \\ && +
\frac{1}{\sqrt{2\sin\theta_0}}
\left\{ 2\sqrt{\epsilon} + \frac{\cos\theta_0}{6 \sin\theta_0} \epsilon\sqrt{\epsilon}
\right\}
\end{eqnarray}
として，$\epsilon$ に小さい値を設定して計算してみるという手もあるだろう。
テイラー展開を利用したこの近似を使うと，振幅 $\theta_0$ の単振り子の周期 $T_p$ は
\begin{eqnarray}
T_p(\theta_0, \epsilon) &=&
\frac{2}{\pi} \int_{0}^{\theta_0} \frac{1}{\sqrt{2(\cos\theta-\cos\theta_0)}}d\theta \\
&\simeq& \frac{\sqrt{2}}{\pi} \int_{0}^{\theta_0-\epsilon} \frac{1}{\sqrt{\cos\theta-\cos\theta_0}}d\theta \\
&& + \frac{\sqrt{2}}{\pi}\frac{1}{\sqrt{\sin\theta_0}}
\left\{ 2\sqrt{\epsilon} + \frac{\cos\theta_0}{6 \sin\theta_0} \epsilon\sqrt{\epsilon}
\right\}
\end{eqnarray}
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}