単振り子の運動方程式から得られるエネルギー保存則から，Python の SymPyを使って（かつ SciPy と NumPy は使わずに）数値積分によって周期を求める。導出については以下を参照：
ライブラリの import
from sympy.abc import *
from sympy import *
# SymPy Plotting Backends (SPB) でグラフを描く
from spb import *
# ax 用
import matplotlib.pyplot as plt
# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']
SymPy の N() で数値積分する
振幅 $\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 = 80^{\circ}$ として計算してみます。
# ラジアンに変換
theta0 = rad(80)
# Tp の計算
sqrt(2)/pi*integrate(1/sqrt(cos(theta)-cos(theta0)), (theta, 0, theta0))
解析的に積分できない場合，SymPy は式をそのまま返します。
これを数値的に求めるには，以下のように N() で囲みます。
Jupyter Notebook では，以下のようにセルの1行目に %%time を入れると，計算にかかった時間を表示します。
%%time
theta0 = rad(80)
N(sqrt(2)/pi *
integrate(1/sqrt(cos(theta)-cos(theta0)), (theta, 0, theta0)))
# けっこう時間がかかる
SymPy の N() で Tp1(th0)
初めから N() で数値積分したを返す関数を Tp1(th0) として定義します。
def Tp1(th0):
"""度で与えられた振幅 th0 から
規格化された単振り子の規格化された周期を求める"""
# 度で与えられた th0 からラジアンへ変換
theta0 = rad(th0)
ans = N(sqrt(2)/pi *
integrate(1/sqrt(cos(theta)-cos(theta0)), (theta, 0, theta0)))
return ans
%%time
Tp1(80)
SymPy は解析的な微分や積分ができるので便利ですが，数値積分の場合だと妙に時間がかかります。
SymPy の積分そのものが遅くてつかいものにならないというわけではありません，念のため。積分範囲（端点）で被積分関数が発散するようなケースを N() で数値計算する場合，いろいろ考えて陰であの手この手を試しているためかと推察します。のちに示すように，この積分は適当な変数変換によって，積分範囲内で被積分関数が発散しない形にできます。その場合には，N() もそんなに時間がかからずに答えを出してくれますよ。
置換積分で被積分関数が発散しないようにする
このへんにまとめているように，変数変換によって
\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 の N() で Tp2(th0)
被積分関数が発散しないように変数変換を行ったあとに，数値積分する関数を Tp2(th0) として定義する。
def Tp2(th0):
theta0 = N(rad(th0))
m = sin(theta0/2)**2
ans = 2/N(pi) * N(integrate(1/sqrt(1-m*sin(t)**2), (t, 0, pi/2)))
return ans
%%time
Tp2(80)
第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() が使える。
elliptic_k(m) $\displaystyle \equiv \int_0^{\pi/2} \frac{dt}{\sqrt{1-m \sin^2 t}}$。 $m \equiv k^2$ であることに注意。
第1種完全楕円積分 elliptic_k() を使って周期を返す関数を Tp3(th0) として定義する。
def Tp3(th0):
theta0 = N(rad(th0))
m = sin(theta0/2)**2
ans = 2/N(pi) * elliptic_k(m)
return ans
%%time
Tp3(80)
Tp3(th0) のグラフ
# 振幅
th0s = [1] + [10*i for i in range(1,10)]
# 周期
Tp3s = []
# 各振幅ごとの周期
for th0 in th0s:
Tp3s += [Tp3(th0)]
print('θ_0 = %2d° のとき，' % th0, end = '')
print('周期 T = %.15f' % Tp3s[-1])
# SymPy Plotting Backends でグラフ
# 細かな設定はバックエンドの matplotlib ax で
# おまじない。これで ax が使える。
fig, ax = plt.subplots()
# x の目盛を 10 刻みに
ax.set_xticks([10*i for i in range(1, 10)])
# グリッドを点線に
ax.grid(ls = 'dotted')
graphics(
list_2d(th0s, Tp3s, scatter = True,
rendering_kw={"ms":3}),
xlabel = "振幅 θ (°)",
ylabel = "規格化された周期 T",
title = "単振り子の振幅と周期",
xlim = (0, 91),
ylim = (0.98, 1.2),
grid = False, # 一旦 False にして...
ax = ax # ax の設定を反映
);
まとめ
SymPy の N(integrate()) で数値積分した場合と第一種完全楕円積分 elliptic_k(m) を使った場合の周期を求めた。N() を使った場合，かなり時間がかかることがわかった。
しかし，計算時間には大差があるものの，数値計算の結果は同じ値を返すことも以下のことからわかる。
th0s = [1] + [10*i for i in range(1,10)]
%%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)):
print('Tp1(%2d)=%.15f' % (th0s[i], Tp1s[i]))
print('Tp2(%2d)=%.15f' % (th0s[i], Tp2s[i]))
print('Tp3(%2d)=%.15f' % (th0s[i], Tp3s[i]))
print('')