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° きざみで大きくしていった場合，ふりこの周期はどうなるか？ というのが問題。
ルンゲ・クッタ法で解くための準備
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
Python で必要なパッケージを import します。
# NumPy も使います
import numpy as np
# Matplotlib でグラフを描きます
import matplotlib.pyplot as plt
# グラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
運動方程式に出てくる円周率 $\pi$ や三角関数 $\sin\theta$ は NumPy で定義されているものを使います。
import numpy as np
として import しているので，NumPy で定義されているものを使う場合には np. をつけます。
例えば，円周率 $\pi = $ np.pi，三角関数 $\sin\theta = $ np.sin(theta) などです。
連立常微分方程式の右辺を定義
# 連立微分方程式の右辺
# 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]
数値解の精度を確認
まず，$\theta_0 = 80^{\circ}$ の場合に，数値計算の精度について確認しておく。
SciPy の scipy.integrate.solve_ivp() では，数値計算の精度に関係するのは rtol と atol である。
# SciPy による2階常微分方程式の数値解法
# solve_ivp
from scipy.integrate import solve_ivp
# 初期条件
T0 = 0
# 80°をラジアンに
theta0 = np.radians(80)
V0 = 0
# 計算範囲
T1 = 2
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) = %15.13f' %(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)
# 計算結果を保存しておくリスト
ansdeg = []
# 振幅の初期値を変えて数値計算...
for th in [10, 60, 80]:
theta0 = np.radians(th)
tmp = solve_ivp(F, [T0, T1], [theta0, V0],
t_eval = t_list,
rtol = 1.e-13,
atol = 1.e-13)
# ansdeg に [th, theta, V] の順に保存しておく
ansdeg.append([th, tmp.y[0], tmp.y[1]])
計算結果のグラフ
数値計算できたら，結果をグラフにしてみます。
なお今回は直接使いませんが，点の個数が非常に多い場合，少し間引いてグラフにしてみようかな？と思うかも知れませんので，参考までにリストの要素を一定間隔ごとにする例を示しておきます。
a = list(range(11))
print(a)
# リストの要素を5ごとに表示
print(a[::5])
# [::-1] とすると逆順に表示
print(a[::-1])
# 逆順にする別解
list(reversed(a))
# ansdeg に保存された計算結果を
# 逆順に使うときの for の使用例
# なぜ逆順に plt.plot させるのか，わかりますか？
for tmp in reversed(ansdeg):
# 凡例に振幅の初期値を表示
key = 'θ0 = '+str(tmp[0])+'°'
# 横軸に t_list，縦軸に tmp[1] (theta) をとって plot
plt.plot(t_list, tmp[1], label = key)
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベル
plt.xlabel("規格化された時間 T")
plt.ylabel("振幅 θ (rad)")
# 横軸縦軸の表示範囲
plt.xlim(0, 2)
plt.ylim(-1.5, 1.5)
# グラフのタイトル
plt.title("単振り子の振幅と周期")
# 判例を表示
plt.legend();
縦軸の単位を度に変えたグラフ
縦軸の振幅をラジアンから度に変えたグラフを描きます。
以下では直接使いませんが，参考までに，Python でリスト内の全ての要素を2倍する例を示しておきます。
例：
# Python の標準リストでは...
a = [1, 2, 3]
# こう書くと，想定？とは別の意味になる
2*a
# Python 的な方法で全ての要素を2倍する
[2*n for n in a]
# NumPy の配列に変換してから2倍すると...
2 * np.array(a)
また，角度を度からラジアンへ変換するのは np.radians() でしたが，逆にラジアンから度への変換は np.degrees() です。
例：
angd = [30, 45, 60]
angr = np.radians(angd)
# 度をラジアンに変換した値を表示
print(angr)
# ラジアンを度に変換した値を表示
print(np.degrees(angr))
for tmp in reversed(ansdeg):
key = 'θ0 = '+str(tmp[0])+'°'
# リスト tmp[1] の全ての要素をラジアンから度に変換
plt.plot(t_list, np.degrees(tmp[1]), label = key)
plt.grid(linestyle='dotted')
plt.xlabel("規格化された時間 T")
plt.xlim(0, 2)
plt.ylim(-90, 90)
plt.ylabel("振幅 θ (°)")
plt.title("単振り子の振幅と周期")
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 tmp in ansdeg:
# tmp の繰り返しはこのインデント
print('θ_0 = %2d° のとき' % tmp[0])
for i in range(1, N-1):
# i の繰り返しはこのインデント
# 速度の向きが変わるかどうかを判定
if tmp[2][i] * tmp[2][i+1] < 0:
# if が成り立つ時の実行文はこのインデント
t = (t_list[i] + t_list[i+1])/2
# 速度の向きが変わる時刻を表示
# end = '' で改行なしに
print(' T = %8.6f' % t, end = '')
# 次の tmp にいくまえに改行
print()
以上の結果から，例えば振幅の初期条件を $\theta_0 = 80^{\circ}$ とした時の周期は $T = 1.137495$，すなわち単振動の場合の周期
$$T_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 th in [10, 60, 80]:
などとなっているところを
for th in range(10, 91, 10):
などとすればいいかも。また，上記の例では半周期ごとに T = と表示しているが，周期は2回目に速度の向きが変わる時間だから，その時間だけを表示するようにしてください。
解答例
ansdeg = []
# ほぼ単振動となるであろう θ = 1°の場合
th = 1
theta0 = np.radians(th)
tmp = solve_ivp(F, [T0, T1], [theta0, V0],
t_eval = t_list,
rtol = 1.e-13,
atol = 1.e-13)
# ansdeg に [th, theta, V] の順に保存しておく
ansdeg.append([th, tmp.y[0], tmp.y[1]])
# 振幅の初期値を変えて数値計算...
for th in range(10, 91, 10):
theta0 = np.radians(th)
tmp = solve_ivp(F, [T0, T1], [theta0, V0],
t_eval = t_list,
rtol = 1.e-13,
atol = 1.e-13)
# ansdeg に [th, theta, V] の順に保存しておく
ansdeg.append([th, tmp.y[0], tmp.y[1]])
# 振幅と周期のグラフを描くために x, y に追加していく
x = []
y = []
for tmp in ansdeg:
print('θ_0 = %2d° のとき ' % tmp[0], end = '')
x.append(tmp[0])
count = 0
for i in range(1, N-1):
if tmp[2][i] * tmp[2][i+1] < 0:
count += 1
if count == 2:
t = (t_list[i] + t_list[i+1])/2
print('周期 T = %7.5f' % t)
y.append(t)
break
# 横軸に振幅 x, 縦軸に周期 y
plt.scatter(x, y, s=10)
plt.grid(linestyle='dotted')
plt.xlabel("振幅 θ (°)")
# 90°のときの点をちゃんと表示させるために
plt.xlim(0, 91)
plt.xticks(list(range(0,91,10)))
plt.ylim(0.95, 1.2)
plt.ylabel("規格化された周期 T")
plt.title("単振り子の振幅と周期");
