単振り子の運動方程式から得られるエネルギー保存則から，Python の SymPyを使って（かつ SciPy と NumPy は使わずに）数値積分によって周期を求める。導出については以下を参照：
必要なモジュールの import
from sympy.abc import *
from sympy import *
# SymPy Plotting Backends (SPB) でグラフを描く
from spb import *
# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']
SymPy の .evalf() で数値積分する
振幅 $\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}
これを Tp(th0) と定義する。
ためしに，$\theta_0 = 60^{\circ}$ として計算してみます。
# ラジアンに変換
theta0 = 60 * pi/180
# Tp の計算
sqrt(2)/pi * integrate(1/sqrt(cos(theta)-cos(theta0)), (theta, 0, theta0))
解析的に積分できない場合，SymPy は式をそのまま返します。
これを数値的に求めるには，以下のように .evalf() をつけます。
_.evalf()
# けっこう時間がかかる
SymPy の .evalf() で Tp1(th0)
初めから .evalf() で数値積分したを返す関数を Tp1(th0) として定義します。
def Tp1(th0):
"""度で与えられた振幅 th0 から
規格化された単振り子の規格化された周期を求める"""
# 度で与えられた th0 からラジアンへ変換
theta0 = th0 * pi/180
ans = (sqrt(2)/pi * integrate(1/sqrt(cos(theta)-cos(theta0)),
(theta, 0, theta0))).evalf()
return ans
Jupyter Notebook では，以下のようにセルの1行目に %%time を入れると，計算にかかった時間を表示します。
%%time
Tp1(60)
# けっこう時間がかかり，使いものにならない。
SymPy は解析的な微分や積分ができるので便利ですが，数値積分の場合だと妙に時間がかかります。
SymPy の積分そのものが遅くてつかいものにならないというわけではありません，念のため。積分範囲（端点）で被積分関数が発散するようなケースを .evalf() で数値計算する場合，いろいろ考えて陰であの手この手を試しているためかと推察します。のちに示すように，この積分は適当な変数変換によって，積分範囲内で被積分関数が発散しない形にできます。その場合には，.evalf() もそんなに時間がかからずに答えを出してくれますよ。
置換積分で被積分関数が発散しないようにする
このへんにまとめているように，変数変換によって
\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}
と書ける。この形にすると，積分範囲内で被積分関数が発散することもない。
SymPy の .evalf() で Tp2(th0)
def Tp2(th0):
theta0 = th0 * pi/180
m = float(sin(theta0/2)**2)
ans = (2/pi * integrate(1/sqrt(1-m*sin(t)**2),
(t, 0, pi/2))).evalf()
return ans
%%time
Tp2(60)
被積分関数が積分範囲内で発散しないように変数変換したあとの置換積分による周期 Tp2() は Tp1() と同じ数値を返しますが，所要時間は Tp1() よりも短くなります。
第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}
SymPy の elliptic_k() で Tp3(th0)
ありがたいことに，SymPy では第1種完全楕円積分 elliptic_k() が使える。
def Tp3(th0):
theta0 = th0 * float(pi)/180
m = sin(theta0/2)**2
ans = 2/float(pi) * elliptic_k(m)
return ans
%%time
Tp3(60)
Tp3(th0) のグラフ
p = plot(Tp3(x), (x, 1, 90),
xlabel = "振幅 θ (°)",
ylabel = "規格化された周期 $T$",
title = "単振り子の振幅と周期",
xlim = (0, 91),
ylim = (0.95, 1.2),
show = False);
ax = p.ax
# x の主目盛を 10 刻みに
ax.set_xticks([i for i in range(10, 91, 10)]);
SPB の plot() で陽関数 $y = Tp3(x)$ をグラフに描く際，最初，上記のように Warning が出るかもしれない。Warning が嫌なら，plot_list() で描くのもいいかもしれない。
# x 軸（横軸）のリスト
x_list = [i for i in range(1,91)]
# y 軸（縦軸）のリスト
y_list = []
for x in x_list:
y_list += [Tp3(x)]
# SymPy Plotting Backends (SPB)
p = plot_list(x_list, y_list,
xlabel = "振幅 θ (°)",
ylabel = "規格化された周期 $T$",
title = "単振り子の振幅と周期",
xlim = (0, 91),
ylim = (0.95, 1.2),
show = False);
ax = p.ax
# x の主目盛を 10 刻みに
ax.set_xticks([i for i in range(10, 91, 10)]);
p.save("SfuriTp3B.svg");
まとめ
SymPy の integrate().evalf() で数値積分した場合と第一種完全楕円積分 elliptic_k(m) を使った場合の周期を求めた。.evalf() を使った場合，かなり時間がかかることがわかった。
しかし，計算時間には大差があるものの，数値計算の結果は同じ値を返すことも以下のことからわかる。
th0_list = [10*i for i in range(1,10,2)]
%%time
TP1 = []
for th0 in th0_list:
TP1 += [Tp1(th0)]
%%time
TP2 = []
for th0 in th0_list:
TP2 += [Tp2(th0)]
%%time
TP3 = []
for th0 in th0_list:
TP3 += [Tp3(th0)]
for i in range(len(th0_list)):
th0 = th0_list[i]
print('Tp1(%2d)=%.15f' % (th0, TP1[i]))
print('Tp2(%2d)=%.15f' % (th0, TP2[i]))
print('Tp3(%2d)=%.15f' % (th0, TP3[i]))
print('')