Maxima を使って運動方程式を数値的に解き，単振り子の振幅と周期の関係を調べる。導出については，以下を参照。
一様重力場中の単振り子
単振り子のひもの長さを $\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 倍の時間まで，という意味ですよ。）
/* 連立微分方程式の右辺 */
F1(V):= V;
F2(theta):= - 4 * %pi**2 * sin(theta);
/* 単振動の場合
F2(theta):= - 4 * %pi**2 * theta; */
まず，$\theta_0 = 80^{\circ}$ の場合に，刻み幅 $h$ と数値計算の精度について確認しておく。
Maxima の Runge-Kutta 法の関数 rk() では，数値計算の精度に関係するのは刻み幅 $h$ だけであるからね。
/* 度をラジアンに変換する関数を準備しておく */
radians(deg):= deg * float(%pi)/180$
/* 初期条件と計算範囲 */
T0: 0$
th: 80$
theta0: radians(th)$ /* ラジアンになおす */
V0: 0$
T1: 2$
/* N を変えて数値計算... */
for N in [200, 2000, 20000] do(
/* 刻み幅 */
h: float((T1 - T0)/N),
/* Runge-Kutta 法で数値計算 */
ansN: rk([F1(V), F2(theta)],
[theta, V], [theta0, V0],
[T, T0, T1, h]),
/* 最後の値を表示し，誤差を確認。 */
printf(true, "N = ~6d, h = ~,8f, V(T=2) = ~15,13f~%",
N, h, ansN[length(ansN)][3])
)$
上の結果から $N=2000$ で小数点以下10桁程度の精度があるとわかる。ここでは $N=20000$ で計算を行うことにしてみる。グラフを描くだけなら $N$ をそんなに大きくとる必要はありませんが，あとで理由がわかります。
/* 初期条件と計算範囲 */
T0: 0$
V0: 0$
T1: 2$
/* N を大きめにとる理由は，あとでわかる */
N: 20000$
/* 刻み幅 */
h: float((T1 - T0)/N)$
/* 振幅の初期値を変えて数値計算... */
for th in [10, 60, 80] do(
theta0: radians(th),
ansdeg[th]: rk([F1(V), F2(theta)],
[theta, V], [theta0, V0],
[T, T0, T1, h])
)$
数値計算できたら，結果をグラフにしてみます。
/* あらかじめグラフにする points のリストを作成しておく */
i: 0$
for th in [10, 60, 80] do(
/* 線の色 */
i: i + 1,
col[th]: i,
/* 凡例に振幅の初期値をいれる */
legend[th]: concat("θ_0 = ", th, "°"),
/* たくさん点を描くと時間がかかるので
プロットするデータを少し間引く */
pts[th]: makelist([ansdeg[th][i][1], ansdeg[th][i][2]], i, 1, N+1, 100)
)$
tmp: makelist([key =legend[th],
color = col[th],
points(pts[th])], th, [80, 60, 10])$
draw2d(
font = "Arial", font_size = 16,
title = "単振り子の振幅と周期",
user_preamble = "set key samplen 1",
/* 表示範囲 */
xrange = [0, 2], yrange = [-1.5, 1.5],
/* グリッド */
grid = true,
xaxis = true,
/* 座標軸のラベル */
xlabel = "規格化された時間 T", ylabel = "振幅 θ (rad)",
point_type = dot,
points_joined = true,
line_width = 1.5,
tmp
)$
縦軸の振幅をラジアンから度に変えたグラフを描きます。
なお，以下では直接使いませんが，参考までに Maxima でリスト内の全ての要素を2倍する例は以下の通りです。
a: [1, 2, 3];
a*2;
/* ラジアンから度へ変換する関数を準備しておく */
degrees(rad):= rad * 180/float(%pi)$
/* あらかじめグラフにする points のリストを作成しておく */
i: 0$
for th in [10, 60, 80] do(
/* 線の色 */
i: i + 1,
col[th]: i,
/* 凡例 */
legend[th]: concat("θ_0 = ", th, "°"),
/* プロットするデータを少し間引き，縦軸をラジアンから度に。 */
pts[th]: makelist([ansdeg[th][i][1],
degrees(ansdeg[th][i][2])], i, 1, N+1, 100)
)$
tmp: makelist([key =legend[th],
color = col[th],
points(pts[th])], th, [80,60,10])$
draw2d(
font = "Arial", font_size = 16,
title = "単振り子の振幅と周期",
user_preamble = "set key samplen 1",
/* 表示範囲 */
xrange = [0, 2], yrange = [-90, 90],
/* グリッド */
grid = true,
xaxis = true,
/* 座標軸のラベル */
xlabel = "規格化された時間 T", ylabel = "振幅 θ (°)",
point_type = dot,
points_joined = true,
line_width = 1.5,
tmp
)$
周期は速度の向きが変わる時刻から…
初速度 $V(0) = 0$ ではじまる振り子の運動は，半周期ごとに $V = 0$ となるはずである。実際の数値計算では，刻み幅 $h$ ごとの値を求めているので，厳密に $V=0$ となる瞬間が得られるとは限らない。したがって，速度の向きが変わる時刻，つまり $V_i \times V_{i+1} \leq 0$ となる時刻を探すことにする。
ただグラフを描くためだけなら，N: 20000 などという大きな値をとる必要はないのだが，このように速度の向きが変わる瞬間をとらえたいときには，なるべく刻み幅を小さくすると良いはず。
このような理由で N の値を大きくしている。
for th in [10, 60, 80] do(
printf(true, "θ_0 = ~2d° のとき，速度の向きが変わる時刻~%", th),
for i:2 thru N do(
/* if then の後に複数行の実行文を書く場合は () で囲む */
if ansdeg[th][i][3] * ansdeg[th][i+1][3] <= 0 then (
t: (ansdeg[th][i][1] + ansdeg[th][i+1][1])/2,
printf(true, " T = ~,5f", t)
)
),
printf(true, "~%") /* 改行 */
)$
以上の結果から，例えば振幅の初期条件を $\theta_0 = 60^{\circ}$ とした時の周期は $T = 1.07315$，すなわち単振動の場合の周期
$$T_0 = 2\pi \sqrt{\frac{\ell}{g}}$$
の $1.07315$ 倍だということになる。
問題
では本題。$\theta_0 = 10^{\circ}$ から $10^{\circ}$ きざみで $90^{\circ}$ まで， $T = 0$ から $ T = 2$ まで解き，振幅 $\theta_0$ に対する周期 $T$ を求めてみよ。
ヒント：
for th in [10, 60, 80]
などとなっているところを
for th:10 thru 90 step 10
となどすればいいかも。また，上記の例では半周期ごとに T = と表示しているが，周期は2回目に速度の向きが変わる時間だから，その時間だけを表示するようにしてください。
解答例
/* 振幅と周期のグラフを描くために x, y に追加していく */
x: []$ /* 振幅 */
y: []$ /* 周期 */
/* 参考として，ほぼ単振動となるであろう θ = 1°の場合 */
th: 1$
theta0: radians(th)$
ansdeg[th]: rk([F1(V), F2(theta)],
[theta, V], [theta0, V0],
[T, T0, T1, h])$
printf(true, "θ_0 = ~2d° のとき ", th)$
count: 0$
for i:2 thru N do(
/* if then の後に複数の実行文を書く場合は () で囲む */
if ansdeg[th][i][3] * ansdeg[th][i+1][3] <= 0 then (
count: count + 1,
/* 周期は2回目に速度の向きが変わる時間だから */
if count = 2 then (
t: (ansdeg[th][i][1] + ansdeg[th][i+1][1])/2,
x: append(x, [th]),
y: append(y, [t]),
printf(true, "周期 T = ~,5f~%", t),
return() /* 周期を表示したら for ループから抜ける */
)
)
)$
/* 振幅の初期値を変えて数値計算... */
for th:10 thru 90 step 10 do(
theta0: radians(th),
ansdeg[th]: rk([F1(V), F2(theta)],
[theta, V], [theta0, V0],
[T, T0, T1, h])
)$
for th:10 thru 90 step 10 do(
printf(true, "θ_0 = ~2d° のとき ", th),
count: 0,
for i:2 thru N do(
/* if then の後に複数の実行文を書く場合は () で囲む */
if ansdeg[th][i][3] * ansdeg[th][i+1][3] <= 0 then (
count: count + 1,
/* 周期は2回目に速度の向きが変わる時間だから */
if count = 2 then (
t: (ansdeg[th][i][1] + ansdeg[th][i+1][1])/2,
x: append(x, [th]),
y: append(y, [t]),
printf(true, "周期 T = ~,5f~%", t),
return() /* 周期を表示したら for ループから抜ける */
)
)
)
)$
/* オプションを設定してグラフを描く */
draw2d(
font = "Arial", font_size = 16,
title = "単振り子の振幅と周期の関係",
/* 表示範囲 */
xrange = [-2, 92], yrange = [0.95, 1.2],
/* グリッド */
grid = true,
/* 座標軸のラベル */
xlabel = "振幅 θ (°)", ylabel = "規格化された周期 T",
/* 点を塗りつぶした丸に */
point_type = 7,
/* 点の大きさ */
point_size = 0.5,
points(x, y)
)$
