単振り子の運動方程式から得られるエネルギー保存則から，Python の SciPy と NumPy を使って（かつ SymPy は使わずに）数値積分によって周期を求める。導出については以下を参照：
必要なモジュールの import
# NumPy も使います
import numpy as np
# 数値積分 scipy.integrate.quad()
from scipy.integrate import quad
# 第1種完全楕円積分 K(k)
from scipy.special import ellipk
# Matplotlib でグラフを描きます
import matplotlib.pyplot as plt
# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
平方根や三角関数などは NumPy の関数 np.sqrt()，np.cos()，np.sin() などを使います。scipy.integrate.quad() については授業で簡単な説明をしています。
数値計算に限ることにし，SymPy の関数を混在させないことにして，NumPy で定義された関数や定数のみを使うようにします。
SciPy の quad() で数値積分する
振幅 $\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}
積分部分は数値積分 scipy.integrate.quad() を使う。
$\theta = \theta_0$ で被積分関数の分母はゼロになってしまうので，何かエラーでも起こるのかと思いきや，scipy.integrate.quad() は優秀で何事もなかったかのように数値積分してしまいます！
quad() は (積分の近似値, 見積もられた近似の絶対誤差) の2つの数値を返します。
SciPy の quad() で Tp1(th0)
# 被積分関数の定義
# quad() を使う場合は積分変数名は x 決め打ち
def f1(x, theta0):
return np.sqrt(2)/(np.pi*np.sqrt(np.cos(x)-np.cos(theta0)))
# 周期 Tp1(theta0) の定義
def Tp1(th0):
"""度で与えられた振幅 th0 から
規格化された単振り子の規格化された周期を求める"""
# 度で与えられた th0 からラジアンへ変換
theta0 = np.radians(th0)
ans = quad(f1, 0, theta0, args=theta0)
# 積分値のみを ans[0] としてとり，周期をかえす
# 参考までに誤差 ans[1] も返す
return ans[0], ans[1]
セルの先頭に %%time を書くと，セルの実行（計算）にかかった時間を表示します。
%%time
Tp1(60)
（SymPy の関数で定義した場合とくらべて）圧倒的に短い時間で計算してくれます。ただし，このままだと，（被積分関数が端点で発散するという問題もあり）誤差が小さくならないようです。
Tp1(th0) のグラフ
# 振幅と周期のグラフを描くために x1, y1 に追加していく
x1 = []
y1 = []
# 参考として，ほぼ単振動となるであろう θ_0 = 1°の場合
th0 = 1
x1 += [th0]
# 数値積分値のみを取り出すには Tp1(th0)[0]
y1 += [Tp1(th0)[0]]
# リストに append したばかりの最後の値を表示
print('θ_0 = %2d° のとき ' % x1[-1], end = '')
print('周期 T = %10.8f' % y1[-1])
# θ_0 を 10° から 90° まで 10° ごとにかえて計算
for th0 in range(10, 91, 10):
x1 += [th0]
y1 += [Tp1(th0)[0]]
print('θ_0 = %2d° のとき ' % x1[-1], end = '')
print('周期 T = %10.8f' % y1[-1])
# 横軸に振幅 x1, 縦軸に周期 y1
plt.scatter(x1, y1, s=10)
plt.grid(linestyle='dotted')
plt.xlabel("振幅 θ (°)")
plt.xlim(0, 91)
plt.xticks(list(range(10,91,10)))
plt.ylim(0.95, 1.2)
plt.ylabel("規格化された周期 T")
plt.title("単振り子の振幅と周期");
被積分関数が発散しないように変数変換して置換積分
このへんにまとめているように，変数変換によって
\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}
と書ける。この形にすると，積分範囲内で被積分関数が発散することもない。
SciPy の quad() で Tp2(th0)
# 被積分関数の定義
# quad() を使う場合は積分変数名は x 決め打ち
def f2(x, theta0):
k = np.sin(theta0/2)
return 2/(np.pi * np.sqrt(1-k**2*np.sin(x)**2))
def Tp2(th0):
theta0 = np.radians(th0)
ans = quad(f2, 0, np.pi/2, args=theta0)
# 積分値のみを ans[0] としてとり，係数をかけて周期をかえす
return ans[0], ans[1]
%%time
Tp2(60)
被積分関数が発散しないように変数変換してから quad() で数値積分すると，誤差が小さく，より短時間で実行しているのがわかります。
第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) \\
&\equiv& \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 の ellipk(m) で Tp3(th0)
ありがたいことに，SciPy では第1種完全楕円積分が使える。 マニュアルは以下を参照：
$$\mbox{ellipk(m)} = K(m) \equiv \int_0^{\frac{\pi}{2}}\frac{1}{\sqrt{1 – m \sin^2 x}} dx$$
第1種完全楕円積分 scipy.special.ellipk(m) を使って，振幅 $\theta_0$ の単振り子の周期 $T_p$ を以下のように定義する。（$m = k^2$ に注意。）
def Tp3(th0):
theta0 = np.radians(th0)
m = np.sin(theta0/2)**2
return 2/np.pi*ellipk(m)
%%time
Tp3(60)
練習
Tp1(th0), Tp2(th0), Tp3(th0)の誤差および計算時間を調べる。
th0 = 60 の場合にはすでに出している。他の角度についても計算してみる。
th0_list = np.append(1, np.arange(10, 91, 10))
th0_list
%%time
TP1 = []
for th0 in th0_list:
TP1 += [Tp1(th0)[0]]
%%time
TP2 = []
for th0 in th0_list:
TP2 += [Tp2(th0)[0]]
%%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('')