(21) 振り子の周期 - 読み書きプログラミング

振り子の運動について考えます。

振り子のHamiltonianは以下のように書けます。

$H(l, \theta)=\frac{m}{2}(\frac{l}{mL})^2-mgL(1-\cos\theta)$

振り子の正準運動方程式はHamiltonianから以下の通りに求められます。

$\frac{d\theta}{dt}=\frac{\partial H}{\partial l}=\frac{l}{mL^2}$

$\frac{dl}{dt}=-\frac{\partial H}{\partial\theta}=-mgL\sin\theta$

ここで、lは角運動量、θは角度、mは質量、gは重力加速度、Lは振り子の長さです。

Euler法でMaximaに周期を計算させてみましょう。
以下の二階微分方程式に直してからスタートします。

$\frac{d^2\theta}{{dt}^2}=-\frac{g}{L}\sin\theta$

初期条件 $\theta=\frac{\pi}{2},\frac{d\theta}{dt}=0$ とします。また簡単のため、L=gと考えます。

/* 5.1.13m */
theta : %pi/2$
dtheta : 0$
dt : 0.02$
t : 0$
while theta > 0 do (
  ddtheta : -sin(theta),
  dtheta : dtheta + ddtheta*dt,
  theta : theta + dtheta * dt,
  t : t + dt)$
period : 4 * t;

$7.440000000000005$

正確な周期は以下の積分で表されます。

$P=4\sqrt{\frac{L}{2g}}\int_{0}^{\frac{\pi}{2}}{\frac{d\theta}{\sqrt{\cos\theta}}}$

L=gとしてこれをMaximaで数値積分してみましょう。

/* 5.1.17m */
p : sqrt(8) * quad_qags(1 / sqrt(cos(theta)), theta, 0, %pi / 2)$
p[1], numer;

$7.416298709205079$

使用した機能

quad_qags

Euler法で求めた値が近かったことがわかりました。

次に数値積分をせずに計算してましょう。

sqrt(8)*integrate(1/sqrt(cos(theta)), theta, 0, %pi/2);
%, numer;

$\sqrt{2}\,\beta\left( \frac{1}{4},\frac{1}{2}\right)$

$7.416298709205502$

使用した機能

βはベータ函数です。こちらの値がより精確な値となります。

θの初期条件を $\frac{\pi}{2}$ ではなく、一般に $\theta_0$ とすると、周期は以下の積分で表されます。

$P=4\sqrt{\frac{L}{g}}\int_{0}^{\theta_0}{\frac{d\theta}{\sqrt{2(\cos\theta-\cos\theta_0)}}}$

これは、式変形により、第一種完全楕円関数K(k)で表すことができます。

$P=4\sqrt{\frac{L}{g}}\int_{0}^{\frac{\pi}{2}}{\frac{d\phi}{\sqrt{1-\sin^2\frac{\theta_0}{2}\sin^2\phi}}=4\sqrt{\frac{L}{g}}K(\sin\frac{\theta_0}{2})$

Maximaの楕円積分ライブラリでも、 $\theta=\frac{\pi}{2}$ の値を確認しておきましょう。

4*elliptic_kc(sin(%pi/4)^2);
%, numer;

$\frac{2\,{\pi }^{\frac{3}{2}}}{{\gamma\left( \frac{3}{4}\right) }^{2}}$

$7.416298709205486$

使用した機能

elliptic_kc

$K(\sin\frac{\pi}{4})$ を求めるのに、elliptic_kc(sin(%pi/4)^2)を使っていることに注意してください。Maximaの楕円積分ライブラリはAbramowitz&Stegunを参照して設計されているため、通常の楕円積分と定義が異なっています。

Hamiltonianまで戻って、Runge-Kutta法を使って、数値的に解いてみましょう。

/* 5.1.25m */
load(dynamics)$
h(v, x) := v ^ 2 / 2 - cos(x)$
dt : 0.02$
phaseplot : rk([diff(h(v, x), v), -diff(h(v, x), x)], [x, v], [float(%pi / 2), 0], [t, 0, 8, dt])$
plot2d([discrete, makelist([x[2], x[3]], x, phaseplot)])$

使用した機能

この計算結果から $\theta$ の符号が変わる時刻を内挿し、周期を計算してみましょう。

/* 5.1.27m */
j : 1$
while j <= length(phaseplot) do (
  if phaseplot[j][2] < 0 then return(),
  j : j + 1)$
before : abs(phaseplot[j - 1][2])$
after : abs(phaseplot[j][2])$
period : 4 * (before*phaseplot[j][1] + after*phaseplot[j - 1][1]) / (after + before)$