ライブラリの import
# NumPy も使います
import numpy as np
# 数値積分 scipy.integrate.quad()
from scipy.integrate import quad
# 第1種完全楕円積分 K(m)
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):
"""度で与えられた振幅 th から
規格化された単振り子の規格化された周期を求める"""
# 度で与えられた th からラジアンへ変換
theta0 = np.radians(th0)
ans = quad(f1, 0, theta0, args=theta0)
#
return ans[0]
セルの先頭に %%time を書くと，セルの実行（計算）にかかった時間を表示します。
%%time
theta0 = np.radians(80)
quad(f1, 0, theta0, args=theta0)
（SymPy の関数で定義した場合とくらべて）圧倒的に短い時間で計算してくれます。ただし，このままだと，（被積分関数が端点で発散するという問題もあり）誤差が $10^{-11}$ 程度あります。
quad() の精度に関するパラメータは epsabs と epsrel です。以下のように設定してみると…
%%time
theta0 = np.radians(80)
quad(f1, 0, theta0, args=theta0, epsabs = 1e-12, epsrel = 1e-12)
絶対誤差は小さくなりますが，theta0 の値によっては以下のような警告が出るようなので，あまりいじらないことにします。
%%time
theta0 = np.radians(1)
quad(f1, 0, theta0, args=theta0, epsabs = 1e-12, epsrel = 1e-12)
Tp1(th0) のグラフ
Python の Matplotlib によるグラフ作成には，plt.*** という関数のみを使った plt (pyplot) 流（pyplot インターフェースとも）と，ax.*** という関数のみを使った ax 流（オブジェクト指向インターフェースとも）の2つの方法があります。
ここでは，以下のページの説明に従い，ax.*** のみを使ってグラフを作成します。
plt.*** のみを使ってグラフを作成する方法については，以下のページにまとめています。
# 振幅
th0s = [1] + [10*i for i in range(1,10)]
# 周期
Tp1s = []
# 各振幅ごとの周期
for th0 in th0s:
Tp1s += [Tp1(th0)]
print('θ_0 = %2d° のとき，' % th0, end = '')
print('周期 T = %10.8f' % Tp1s[-1])
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
# 横軸に振幅 th0s, 縦軸に周期 Tp1s
ax.scatter(th0s, Tp1s, s=10)
ax.set_title("単振り子の振幅と周期");
ax.set_xlabel("振幅 θ (°)")
ax.set_ylabel("規格化された周期 T")
ax.set_xlim(0, 91)
ax.set_ylim(0.98, 1.2)
ax.set_xticks(list(range(10,91,10)))
ax.grid(ls='dotted');
被積分関数が発散しないように変数変換して置換積分
このへんにまとめているように，変数変換によって
\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]
%%time
theta0 = np.radians(80)
quad(f2, 0, np.pi/2, args=theta0)
被積分関数が発散しないように変数変換してから 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種完全楕円積分が使える。 マニュアルは以下を参照：
ellipk(m) $\displaystyle = 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(80)
○練習：誤差と計算時間の評価
Tp1(th0), Tp2(th0), Tp3(th0)の誤差および計算時間を調べる。
th0 = 80 の場合にはすでに出している。他の角度についても計算してみる。
%%time
Tp1s = []
for th0 in th0s:
Tp1s += [Tp1(th0)]
%%time
Tp2s = []
for th0 in th0s:
Tp2s += [Tp2(th0)]
%%time
Tp3s = []
for th0 in th0s:
Tp3s += [Tp3(th0)]
for i in range(len(th0s)):
th0 = th0s[i]
print('Tp1(%2d)=%.15f' % (th0, Tp1s[i]))
print('Tp2(%2d)=%.15f' % (th0, Tp2s[i]))
print('Tp3(%2d)=%.15f' % (th0, Tp3s[i]))
print('')