単振り子の運動方程式から得られるエネルギー保存則から，Maxima を使って数値積分によって周期を求める。導出については以下を参照：
数値積分 quad_qags()
振幅 $\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}
積分部分は数値積分 quad_qags() を使う。（romberg() ではなく quad_qags() を使う理由はすぐわかります。）
この積分が一見やっかいなのは，積分範囲の上限 $\theta = \theta_0$ で被積分関数の分母がゼロになることです。
$\theta = \theta_0$ で被積分関数が発散してしまうので，何かエラーでも起こるのかと思いきや， quad_qags() は大変優秀で何事もなかったかのように数値積分してしまいます！
ためしに $\theta_0 = 60^{\circ}$ の場合の積分値を数値積分 quad_qags() で求めてみます。
/* 度をラジアンに変換する関数を準備しておく */
radians(deg):= deg * float(%pi)/180$
th: 60$
theta0: radians(th)$
quad_qags(1/sqrt(cos(theta)-cos(theta0)), theta, 0, theta0);
quad_qags() は４つの要素のリストを返します:
積分の近似値
見積もられた近似の絶対誤差
非積分関数の評価数
エラーコード
積分を使って定義された周期を Tp1(th) と定義する。
Tp1(th):= block(
/* 度で入力した th をラジアンになおす */
theta0: radians(th),
/* quad_gags() で数値積分し，tmp に格納する */
tmp: quad_qags(1/sqrt(cos(theta)-cos(theta0)), theta, 0, theta0),
/* 積分の近似値を tmp[1] として取り出し，係数をかけて周期をかえす */
return(sqrt(2.0)/float(%pi) * tmp[1])
)$
/* 振幅と周期のグラフを描くために x, y に追加していく */
x: []$ /* 振幅 */
y: []$ /* 周期 */
/* 参考として，ほぼ単振動となるであろう θ_0 = 1°の場合 */
th: 1$
x: append(x, [th])$
y: append(y, [Tp1(th)])$
/* リストに append したばかりの最後の値を表示 */
printf(true, "θ_0 = ~2d° のとき ", x[length(x)])$
printf(true, "周期 T = ~8,6f~%", y[length(x)])$
/* θ_0 を 10° から 90° まで 10° ごとにかえて計算 */
for th:10 thru 90 step 10 do(
x: append(x, [th]),
y: append(y, [Tp1(th)]),
printf(true, "θ_0 = ~2d° のとき ", x[length(x)]), /* 改行しない */
printf(true, "周期 T = ~8,6f~%", y[length(x)]) /* ~% で改行 */
)$
/* 振幅と周期のグラフを描く */
draw2d(
font = "Arial", font_size = 16,
title = "単振り子の振幅と周期の関係",
/* 表示範囲 */
xrange = [0, 92], yrange = [0.95, 1.2],
/* グリッド */
grid = true,
/* 座標軸のラベル */
xlabel = "振幅 θ (°)", ylabel = "規格化された周期 T",
/* 点を塗りつぶした丸に */
point_type = 7,
/* 点の大きさ */
point_size = 0.5,
points(x, y)
)$
第1種完全楕円積分 elliptic_kc(m)
ありがたいことに，Maxima では以下のように定義された第一種完全楕円積分が elliptic_kc(m) として使えます。
$$\mbox{elliptic_kc(m)} = K(m) \equiv \int_0^{\frac{\pi}{2}}\frac{1}{\sqrt{1 – m \sin^2 x}} dx$$
第1種完全楕円積分 elliptic_kc(m) を使って，振幅 $\theta_0$ の単振り子の周期 $T_p$ を以下のように定義する。（$m = k^2$ に注意。）
Tp(th):= block(
theta0: radians(th),
m: sin(theta0/2)**2,
return(2/float(%pi) * elliptic_kc(m))
)$
/* 参考として，ほぼ単振動となるであろう θ_0 = 1°の場合を表示 */
th: 1$
printf(true, "θ_0 = ~2d° のとき ", th)$
printf(true, "周期 T = ~14,12f~%", Tp(th))$
/* θ_0 を 10° から 90° まで 10° ごとにかえて表示 */
for th:10 thru 90 step 10 do(
printf(true, "θ_0 = ~2d° のとき ", th),
printf(true, "周期 T = ~14,12f~%", Tp(th))
)$
/* 振幅と周期のグラフを描く */
draw2d(
font = "Arial", font_size = 16,
title = "単振り子の振幅と周期の関係",
/* 表示範囲 */
xrange = [0, 90], yrange = [0.95, 1.2],
/* グリッド */
grid = true,
/* 座標軸のラベル */
xlabel = "振幅 θ (°)", ylabel = "規格化された周期 T",
/* 陽関数 y = f(x) としてグラフにする */
explicit(Tp(t), t, 1, 90), file_name = mfuriE02A, terminal = 'svg
)$
