Python を使って運動方程式を数値的に解き，単振り子の振幅と周期の関係を調べる。導出については，以下を参照。
一様重力場中の単振り子
単振り子のひもの長さを $\ell$，重力加速度の大きさを $g$，鉛直下向きからの振れ角を $\theta$ とすると，単振り子の運動方程式は
$$\frac{d^2\theta}{dt^2} = – \frac{g}{\ell} \sin\theta$$
となる。導出については，以下を参照。
振れ角が小さい場合には単振動となること
$|\theta| \ll 1$ の場合には，$\sin\theta \simeq \theta$ と近似することで，運動方程式は
$$\frac{d^2\theta}{dt^2} = – \frac{g}{\ell} \theta$$
となり，これは単振動となる。単振動の周期 $\tau_0$ は
$$\tau_0 = 2 \pi \sqrt{\frac{\ell}{g}}$$
となり，単振り子のひもの長さと重力加速度だけで決まり，初期条件で与えられる振れ角 $\theta_0$ に依存しない。これを「振り子の等時性」と言ったりする。
無次元化
単振動の周期 $\tau_0$ で無次元化した時間 $T$ を
$$T \equiv \frac{t}{\tau_0}$$
と定義する。この無次元化された時間座標で書いた運動方程式は
$$\frac{d^2 \theta}{dT^2} = – 4 \pi^2 \sin\theta$$
となる。
問題：振れ角が大きい場合の周期は？
振れ角が大きい，つまり $\sin\theta \simeq \theta$
と近似できない場合は，振り子の周期は一般に振れ角に依存するのではないかと考えられる。
そこで，$T = 0$ での振れ角 $\theta_0$
をたとえば 10° から 90° まで 10° きざみで大きくしていった場合，ふりこの周期はどうなるか？ というのが問題。
ルンゲ・クッタ法による数値解法：solve_ivp
2階常微分方程式をルンゲ・クッタ法で数値的に解くためには，連立1階微分方程式系になおす。
\begin{eqnarray}
\frac{d\theta}{dT} &=& F_1(V) = V \\
\frac{dV}{dT} &=& F_2(\theta) = – 4 \pi^2 \sin\theta
\end{eqnarray}
これを初期条件 $T = 0$ で
\begin{eqnarray}
\theta(0) &=& \theta_0 \\
V(0) &=& 0
\end{eqnarray}
とし，$\theta_0 = 10^{\circ}, \ 60^{\circ}, \ 80^{\circ}$ の場合に $T_0 = 0$ から $ T_1 = 2$ まで解く。
（くどいようですが，規格化された時間で $T_1 = 2$ というのは，単振動の場合の周期 $\tau_0$ の 2 倍の時間まで，という意味ですよ。）
必要なモジュールの import
# NumPy も使います
import numpy as np
# Matplotlib でグラフを描きます
import matplotlib.pyplot as plt
# グラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
# SciPy による2階常微分方程式の数値解法
# solve_ivp
from scipy.integrate import solve_ivp
運動方程式に出てくる円周率 $\pi$ や三角関数 $\sin\theta$ は NumPy で定義されているものを使います。
import numpy as np
として import しているので，NumPy で定義されているものを使う場合には np. をつけます。
例えば，円周率 $\pi = $ np.pi，三角関数 $\sin\theta = $ np.sin(theta) などです。
連立常微分方程式の右辺を定義
\begin{eqnarray}
\frac{d\theta}{dT} &=& F_1(V) = V \\
\frac{dV}{dT} &=& F_2(\theta) = – 4 \pi^2 \sin\theta
\end{eqnarray}
# 連立微分方程式の右辺
# scipy.integrate.solve_ivp() は
# dy/dt = F(t, y) の形を解く
def F(t, y):
# solve.ivp() に渡す関数の引数の順番に注意。t が先
# たとえ t に依存しなくても，必ずこう書く。
# 以下のようなインデント（行頭の字下げ）部分が定義
# Python ではインデントには意味があるのでしたね。
theta = y[0]
V = y[1]
F1 = V
F2 = -4 * np.pi**2 * np.sin(theta)
# 単振動の場合
# F2 = -4 * np.pi**2 * theta
return [F1, F2]
数値解の精度を確認
初期条件 $T = 0$ で
\begin{eqnarray}
\theta(0) &=& \theta_0 \\
V(0) &=& 0
\end{eqnarray}
まず，$\theta_0 = 80^{\circ}$ の場合に，数値計算の精度について確認しておく。
SciPy の scipy.integrate.solve_ivp() では，数値計算の精度に関係するのは rtol と atol である。
# 初期条件
T0 = 0
# 80°をラジアンに
theta0 = np.radians(80)
V0 = 0
# 計算範囲
T1 = 2
# ちょっと N が大きいがあとで意味がわかる
N = 2000
# 計算範囲：T0 から T1 までを N 等分したリストを作成
# このリストの時刻ごとに計算結果を出力する
t_list = np.linspace(T0, T1, N+1)
# 精度を変えて結果を確認。
for TOL in [1.e-12, 1.e-13, 4.e-14]:
# 以下のようなインデント（行頭の字下げ）部分が繰り返し範囲
# Python ではインデントには意味があるのでしたね。
ansN = solve_ivp(F, [T0, T1], [theta0, V0],
t_eval = t_list,
rtol = TOL,
atol = TOL)
# 最後の値を表示し，誤差を確認。単振動なら V(T=2)=0
# ansN.y[0] は theta, ansN.y[1] は V
# Python ではリストの最後（最後尾）の値は [-1] で参照する
print('TOL = %.2e, V(T=2) = %.12f' %(TOL, ansN.y[1][-1]))
上の結果から TOL = 1.e-13 で小数点以下 10 桁程度の精度があるとわかる。以後の計算は，これくらいで行う。
# 初期条件と計算範囲
T0 = 0
V0 = 0
T1 = 2
# ちょっと N が大きいがあとで意味がわかる
N = 200000
t_list = np.linspace(T0, T1, N+1)
# 初期条件 θ0 = 80° をラジアンに
th0 = 80
theta0 = np.radians(th0)
sol80 = solve_ivp(F, [T0, T1], [theta0, V0],
t_eval = t_list,
rtol = 1.e-13,
atol = 1.e-13)
solve_ivp() による数値解は sol80 に格納されています。値を参照するには以下のように。
# t_eval すなわち t_list は
sol80.t
# theta すなわち y[0] は
sol80.y[0]
# V すなわち y[1] は
sol80.y[1]
計算結果のグラフ
数値計算できたら，結果をグラフにしてみます。横軸に（規格化された）時間 $T$，縦軸に角度 $\theta$ をとってグラフにします。
plt.plot(sol80.t, sol80.y[0]);
初期条件を凡例に記し，サイエンティフィックなグラフらしく grid（格子線）等をつけてみます。
# 横軸縦軸のラベル。LaTeX 表記は $ で囲む。\ は \\ に。
plt.xlabel("規格化された時間 $T$")
plt.ylabel("振れ角 $\\theta$（ラジアン）")
# 凡例 label の設定。LaTeX 表記は $ で囲む。\ は \\ に。
key = "$\\theta_0 = 80°$"
# グリッド（格子線）の表示。私は点線が好み。
plt.grid(linestyle='dotted')
plt.plot(sol80.t, sol80.y[0], label = key)
# 凡例の表示
plt.legend();
縦軸の単位を度に変えたグラフ
縦軸の振れ角をラジアンから度に変えたグラフを描きます。度からラジアンへの変換は np.radians() でした。ラジアンから度へは np.degrees() です。
# 横軸縦軸のラベル。LaTeX 表記は $ で囲む。\ は \\ に。
plt.xlabel("規格化された時間 $T$")
plt.ylabel("角度 $\\theta$ (°)")
# 凡例 label の設定。LaTeX 表記は $ で囲む。\ は \\ に。
key = "$\\theta_0 = 80$ (°)"
# グリッド（格子線）の表示。私は点線が好み。
plt.grid(linestyle='dotted')
plt.plot(sol80.t, np.degrees(sol80.y[0]), label = key)
# 凡例の表示
plt.legend();
振れ角の初期値を変えて計算
$\theta_0 = 80^{\circ}$ の場合は数値計算できたので，$\theta_0 = 60^{\circ}, 10^{\circ}$，そして比較のために，ほぼ単振動と思われる $\theta_0 = 1^{\circ}$ の場合も計算します。
$\theta_0 = 80^{\circ}$ の場合のコードを4回分コピペして… というのもなんですので，for 文を使って繰り返し処理してみます。
# 振れ角の初期値以外は上と同じにするので省略
# 計算結果を保存しておくリスト
# sols に [th0, theta, V] の順に保存しておく
sols = []
for th0 in [80, 60, 10, 1]:
theta0 = np.radians(th0)
sol = solve_ivp(F, [T0, T1], [theta0, V0],
t_eval = t_list,
rtol = 1.e-13,
atol = 1.e-13)
sols.append([th0, sol.y[0], sol.y[1]])
計算結果をあわせてグラフに。（$\theta_0 = 1^{\circ}$ の場合を省略する場合は，for sol in sols[:3] に）
# 横軸縦軸のラベル。LaTeX 表記は $ で囲む。\ は \\ に。
plt.xlabel("規格化された時間 $T$")
plt.ylabel("振れ角 $\\theta$°")
for sol in sols:
# 凡例 label の設定。LaTeX 表記は $ で囲む。\ は \\ に。
key = "$\\theta_0 = %2d°$" % sol[0]
plt.plot(t_list, np.degrees(sol[1]), label = key)
# グリッド（格子線）の表示。私は点線が好み。
plt.grid(linestyle='dotted')
plt.xlim(0,2)
# 凡例の表示
plt.legend();
周期は速度の向きが変わる時刻から…
初速度 $V(0) = 0$ ではじまる振り子の運動は，半周期ごとに $V = 0$ となるはずである。実際の数値計算では，刻み幅
$$h \equiv \frac{T_1 – T_0}{N}$$
ごとの値を求めているので，厳密に $V=0$ となる瞬間が得られるとは限らない。したがって，速度の向きが変わる時刻，つまり $V_i \times V_{i+1} \leq 0$ となる時刻を探すことにする。
ただグラフを描くためだけなら，N = 200000 などという大きな値をとる必要はないのだが，このように速度の向きが変わる瞬間をとらえたいときには，なるべく刻み幅を小さくすると良いはず。
このような理由で N の値を大きくしている。
for sol in sols:
# インデント（字下げ）が大事！
print('$θ_0$ = %2d° のとき' % sol[0])
for i in range(1, N-1):
# i の繰り返しはこのインデント
# 速度の向きが変わるかどうかを判定
if sol[2][i] * sol[2][i+1] <= 0:
# if が成り立つ時の実行文はこのインデント
t = (t_list[i] + t_list[i+1])/2
# 速度の向きが変わる時刻を表示
# end = '' で改行なしに
print(' T = %8.6f' % t, end = '')
# 次にいくまえに改行
print()
以上の結果から，例えば振れ角の初期条件を $\theta_0 = 80^{\circ}$ とした時の周期は $T = 1.137495$，すなわち単振動の場合の周期
$$\tau_0 = 2\pi \sqrt{\frac{\ell}{g}}$$
の $1.137495$ 倍だということになる。
問題
では本題。$\theta_0 = 10^{\circ}$ から $10^{\circ}$ きざみで $90^{\circ}$ まで， $T = 0$ から $ T = 2$ まで解き，その結果から振れ角 $\theta_0$ に対する周期 $T$ を求めてみよ。
また，横軸に振れ角 $\theta_0$，縦軸に周期 $T$ をとったグラフを描け。
ヒント：
for th0 in [80, 60, 10, 1]:
などとなっているところを
for th0 in range(10, 91, 10):
などとすればいいかも。また，上記の例では半周期ごとに T = と表示しているが，周期は2回目に速度の向きが変わる時間だから，そのときの時間を周期として保存するように。
大学の力学等でも，単振り子の運動方程式は，振れ角が非常に小さくて $\sin\theta \simeq \theta$ と近似できるときの話しかしないかもしれない。
近似しないときには，単振り子の運動方程式は解析的に解けないからである。
ここでは，2階常微分方程式を数値的に解くことによって，単振り子の周期は振れ角に依存すること，振れ角が大きくなると周期が長くなることがわかった。「振り子の等時性は成り立たない」ことを理解できた！という意味で，これは重要な経験だと思うが，いかがでしょうか。