単振り子の運動方程式から得られるエネルギー保存則から，Python を使って数値積分によって周期を求める。導出については以下を参照：
数値積分 scipy.integrate.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}
これを Tp1(theta0) と定義する。積分部分は数値積分 scipy.integrate.quad() を使う。
Python で必要なパッケージを import します。
# NumPy も使います
import numpy as np
# Matplotlib でグラフを描きます
import matplotlib.pyplot as plt
# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
平方根や三角関数などは NumPy の関数 np.sqrt()，np.cos()，np.sin() などを使います。scipy.integrate.quad() については授業で簡単な説明をしています。
# 数値積分 scipy.integrate.quad() を import
from scipy.integrate import quad
# 被積分関数の定義
def f(x, theta0):
return 1/np.sqrt(np.cos(x)-np.cos(theta0))
$\theta = \theta_0$ で被積分関数の分母はゼロになってしまうので，何かエラーでも起こるのかと思いきや，scipy.integrate.quad() は大変優秀で何事もなかったかのように数値積分してしまいます！
ためしに $\theta_0 = 60^{\circ}$ の場合の積分値を数値積分 scipy.integrate.quad() で求めてみます。
# 振幅を°で
th = 60
theta0 = np.radians(th)
quad(f, 0, theta0, args=theta0)
quad() は (積分の近似, 見積もられた近似の絶対誤差) を返します。
# 周期 Tp1(theta0) の定義
def Tp1(th):
"""度で与えられた振幅 th から
規格化された単振り子の周期を求める"""
# 度で与えられた th からラジアンへ変換
theta0 = np.radians(th)
tmp = quad(f, 0, theta0, args=theta0)
# 積分値のみを tmp[0] としてとる
return np.sqrt(2)/np.pi * tmp[0]
help(Tp1)
# 振幅と周期のグラフを描くために x, y に追加していく
x = []
y = []
# 参考として，ほぼ単振動となるであろう θ_0 = 1°の場合
th = 1
x.append(th)
y.append(Tp1(th))
# リストに append したばかりの最後の値を表示
print('θ_0 = %2d° のとき ' % x[-1], end = '')
print('周期 T = %10.8f' % y[-1])
# θ_0 を 10° から 90° まで 10° ごとにかえて計算
for th in range(10, 91, 10):
x.append(th)
y.append(Tp1(th))
print('θ_0 = %2d° のとき ' % x[-1], end = '')
print('周期 T = %10.8f' % y[-1])
# 横軸に振幅 x, 縦軸に周期 y
plt.scatter(x, y, 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("単振り子の振幅と周期");
第1種完全楕円積分 elliptic_kc(m)
ありがたいことに，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$ に注意。）
from scipy.special import ellipk
def Tp(th):
theta0 = np.radians(th)
m = np.sin(theta0/2)**2
return(2/np.pi*ellipk(m))
# 振幅と周期のグラフを描くために x1, y1 に追加していく
x1 = []
y1 = []
th = 1
x1.append(th)
y1.append(Tp(th))
print('θ_0 = %2d° のとき ' % x1[-1], end = '')
print('周期 T = %14.12f' % y1[-1])
for th in range(2, 91):
x1.append(th)
y1.append(Tp(th))
# 10°ごとに表示
if th % 10 == 0:
print('θ_0 = %2d° のとき ' % x1[-1], end = '')
print('周期 T = %14.12f' % y1[-1])
# 横軸に振幅 x1, 縦軸に周期 y1
# 1°ごとのデータがあるので，plot で線に
plt.plot(x1, y1)
plt.grid(linestyle='dotted')
plt.xlabel("振幅 θ (°)")
plt.xlim(0, 90)
plt.ylim(0.95, 1.2)
plt.ylabel("規格化された周期 T")
plt.title("単振り子の振幅と周期");
plt.savefig('pfuriE02A.svg');
