葛西 真寿　弘前大学大学院理工学研究科

一様重力場中の単ふりこ

gnuplot を使って一様重力場中の単ふりこの問題を数値的に解きます。

振幅が小さい時には単振動になり，振り子の周期は振幅に依存しないことが知られていますが，振幅が大きいときの（「十分小さい」という近似が成り立たない場合の）ふりこの周期は，振幅にどう依存するでしょうか？

単ふりこの運動方程式

単ふりこのひもの長さを $\ell$，重力加速度の大きさを $g$，鉛直下向きからの振れ角を $\theta$ とすると，単ふりこの運動方程式は

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

となることを示す。

単ふりこ

運動方程式

（以下，Jupyter Notebook の数式表示では，3次元ベクトルをboldsymbol $\boldsymbol{r}$ のように表すと普通の $r$ と区別がつきにくい懸念があるので，$\vec{r}$ のように表すことにする。）

単ふりこの先端のおもりは，ひもの張力 $\vec{T}$ と重力 $ m\vec{g}$ を受ける。鉛直下向きを $+z$ 方向とすると $\vec{g} = (0, 0, g)$ は重力加速度ベクトル, $g$ は重力加速度の大きさである。

したがって，運動方程式は

$$m \frac{d^2\vec{r}}{dt^2} = m\vec{g} + \vec{T}$$

保存量（角運動量の鉛直成分）

この式の両辺に位置ベクトル $\vec{r}$ の外積をかける。$\vec{T} \propto - \vec{r}$ であることに注意すると，$\vec{r}\times \vec{T} = \vec{0}$ であるから，運動方程式は以下のようになる。

$$\frac{d}{dt}\left(m \vec{r}\times\frac{d\vec{r}}{dt}\right) = m \vec{r}\times\vec{g}$$

この両辺に一定のベクトルである重力加速度ベクトル $\vec{g}$ の内積をかけてやると，

$$ \frac{d}{dt}\left\{\vec{g}\cdot\left(m\vec{r}\times\frac{d\vec{r}}{dt}\right)\right\} = 0 $$

この式は，角運動量 $\displaystyle \vec{L} \equiv m \vec{r}\times\frac{d\vec{r}}{dt}$ の鉛直成分，つまり重力加速度ベクトルにそった成分が保存することを意味する。

運動が平面内に限られる理由

この式を成分で具体的に書き下すと（両辺を $m g$ で割ったのち）

$$\frac{d}{dt}\left( x \frac{dy}{dt} - y \frac{dx}{dt}\right) = \frac{d}{dt}\left\{ x^2 \frac{d}{dt}\left(\frac{y}{x}\right) \right\} = 0 $$

極座標 $$\vec{r} = (x, y, z) = (\ell \sin\theta \cos\varphi, \ell \sin\theta \sin\varphi, \ell \cos\theta) $$

を使って書くと，

$$\frac{d}{dt}\left\{ x^2 \frac{d}{dt}\left(\frac{y}{x}\right) \right\} = \ell^2 \frac{d}{dt} \left(\sin^2\theta \frac{d\varphi}{dt} \right) = 0 $$

微分してゼロということは，中身が一定であることだから，

$$ \sin^2\theta \frac{d\varphi}{dt} = \mbox{const.} \equiv l_z $$

ここで，初期条件として $t = 0$ で $\displaystyle\frac{d\varphi}{dt} = 0$ とすると，$l_z = 0$ となり常に $\displaystyle\frac{d\varphi}{dt} = 0$，つまり $\varphi = \mbox{const.}$ でありつづけることが示される。

したがって，単ふりこの運動が（慣性系では）一定の平面内に限られるのは，以下の理由による：

  • 角運動量の鉛直成分（重力加速度ベクトルに平行な成分）は保存される。
    （上記の場合は，重力加速度ベクトルの向きは $z$ 方向なので，角運動量の $z$ 成分 $L_z = m \ell^2 l_z$ が一定である。）

  • 初期条件として $L_z = 0$ とすると，ずーっと $L_z = 0$ でありつづける。
    （上記の場合は初期条件として $\displaystyle \frac{d\varphi}{dt} = 0$ とすると，ずーっと $\displaystyle \frac{d\varphi}{dt} = 0$でありつづける，つまり $\varphi = \mbox{const.}$）

  • したがって，運動は $\varphi$ が一定の平面内に限られる。

すなわち，初期条件として，$\displaystyle y = 0, \frac{dy}{dt} = 0$ とすると，以後の 単ふりこの運動は $y = 0$ の $xz$ 平面内に限られることがわかる。

最終的な運動方程式

そこでこれ以後は，$\vec{r} = (x, y, z) = (\ell\sin\theta, 0, \ell\cos\theta)$ とおく。

すると，$\vec{r}$ と外積をとった運動方程式の残りの成分のうち，$x$ 成分は $ 0 = 0$ となり， $y$ 成分は以下のようになる。

\begin{eqnarray} \frac{d}{dt} \left( z \frac{dx}{dt} - x \frac{dz}{dt} \right) &=& -x g \\ \frac{d}{dt} \left\{z^2 \frac{d}{dt}\left(\frac{x}{z}\right)\right\} &=& -x g \\ \ell^2 \frac{d}{dt} \left\{ \cos^2\theta \frac{d}{dt}\left( \frac{\sin\theta}{\cos\theta}\right) \right\} &=& - \ell g \sin\theta \\ \ell^2 \frac{d^2\theta}{dt^2} =&=& - \ell g \sin\theta \end{eqnarray}

であるから，運動方程式は最終的に

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

となる。

振幅が小さい場合

$|\theta| \ll 1$ の場合には，$\sin\theta \simeq \theta$ と近似することにより，運動方程式は

$$ \frac{d^2\theta}{dt^2} \simeq - \frac{g}{\ell} \theta$$

となり，これは単振動となり，一般解はただちに以下のように求まる。

$$\theta = A \cos \omega t + B \sin \omega t, \quad \omega \equiv \sqrt{\frac{g}{\ell}}$$

$t = 0$ で $\theta = \theta_0, \ \dot{\theta} = 0$ という初期条件を課すと積分定数 $A, \ B$ が決まり，解は以下のようになる。

$$\theta = \theta_0 \cos \omega t$$

単振動の周期

単振動の周期 $T_0$ は，角振動数 $\omega$ から

$$ T_0 = \frac{2\pi}{\omega} = 2\pi \sqrt{\frac{\ell}{g}} $$

のように求まり，これは初期条件で与えられる振れ幅（振幅）$\theta_0$ に依存しない。

振幅が大きい場合の周期は？

振幅が大きい，つまり $\sin\theta \simeq \theta$ と近似できない場合は，単振動ではない運動をするはずだから，ふりこの周期は一般に振幅に依存するのではないかと考えられる。

そこで，振幅 $\theta_0$ をたとえば $10$° から $90$° まで $10$° きざみで大きくしていった場合，ふりこの周期はどうなるか？　というのが問題。

無次元化された運動方程式

解析的に解けない場合は，数値的に解くことになる。そんなとき，解くべき方程式を無次元化しておくとよい。

ここでは，単振動の場合の周期 $T_0$ の2乗を運動方程式の両辺にかける。 すると，右辺は

$$- T_0^2 \times \frac{g}{\ell} \sin\theta = - 4\pi^2 \sin\theta$$

となる。また，規格化された時間 $\tau$ を

$$ \tau \equiv \frac{t}{T_0}$$

のように定義すると

$$ T_0^2 \frac{d^2}{dt^2} = \frac{d^2}{d\tau^2}$$

となるので，無次元化された運動方程式は以下のようになる。

$$ \frac{d^2\theta}{d\tau^2} = - 4\pi^2 \sin\theta$$

まず，この2階常微分方程式を数値的に解いてみる。

運動方程式の数値解

ルンゲ・クッタ法

gnuplot で2階常微分方程式系をルンゲ・クッタ法で数値的に解く方法については，以下のリンクを参照。

無次元化された運動方程式

$$ \frac{d^2\theta}{d\tau^2} = - 4\pi^2 \sin\theta$$

を以下のような連立1階常微分方程式の系にしてルンゲ・クッタ法で数値的に解いています。

\begin{eqnarray} \theta &\Rightarrow& x \\ \frac{d\theta}{d\tau} &\Rightarrow& v\\ \tau &\Rightarrow& t \\ \ \\ \frac{dx}{dt} &=& v \equiv F_1(x,v,t) \\ \frac{dv}{dt} &=& - 4\pi^2 \sin x \equiv F_2(x,v,t) \end{eqnarray}

まず，$\theta_0 = 60$° の場合に，刻み幅を変えてみて精度を確認してみます。

In [1]:
F1(x, v, t) = v
F2(x, v, t) = -4*pi**2*sin(x)

t0 = 0.0
t1 = 2.0
# 振幅
x0 = 60*pi/180
v0 = 0.0

N = 2000
h = (t1 - t0)/N

# 計算結果を配列に入れます。
array T[N]
array X[N]
array V[N]

t = t0; x = x0; v = v0

do for [i=1: N]{
  k1 = h*F1(x, v, t)
  m1 = h*F2(x, v, t)
  k2 = h*F1(x+0.5*k1, v+0.5*m1, t+0.5*h)
  m2 = h*F2(x+0.5*k1, v+0.5*m1, t+0.5*h)
  k3 = h*F1(x+0.5*k2, v+0.5*m2, t+0.5*h)
  m3 = h*F2(x+0.5*k2, v+0.5*m2, t+0.5*h)
  k4 = h*F1(x+k3, v+m3, t+h)
  m4 = h*F2(x+k3, v+m3, t+h)
  X[i] = x + (k1 + 2.*k2 + 2.*k3 + k4)/6
  V[i] = v + (m1 + 2.*m2 + 2.*m3 + m4)/6
  T[i] = t + h 
  t = T[i]
  x = X[i]
  v = V[i]
}
print sprintf("t =%5.2f, theta =%18.15f, h = %f", T[N], X[N], h)

# N を変えてもう一度。精度確認のため。
N = 10000
h = (t1 - t0)/N

# 計算結果を配列に入れます。
array T[N]
array X[N]
array V[N]

t = t0; x = x0; v = v0

do for [i=1: N]{
  k1 = h*F1(x, v, t)
  m1 = h*F2(x, v, t)
  k2 = h*F1(x+0.5*k1, v+0.5*m1, t+0.5*h)
  m2 = h*F2(x+0.5*k1, v+0.5*m1, t+0.5*h)
  k3 = h*F1(x+0.5*k2, v+0.5*m2, t+0.5*h)
  m3 = h*F2(x+0.5*k2, v+0.5*m2, t+0.5*h)
  k4 = h*F1(x+k3, v+m3, t+h)
  m4 = h*F2(x+k3, v+m3, t+h)
  X[i] = x + (k1 + 2.*k2 + 2.*k3 + k4)/6
  V[i] = v + (m1 + 2.*m2 + 2.*m3 + m4)/6
  T[i] = t + h 
  t = T[i]
  x = X[i]
  v = V[i]
}
print sprintf("t =%5.2f, theta =%18.15f, h = %f", T[N], X[N], h)
Out[1]:
t = 2.00, theta = 0.695210056794716, h = 0.001000
t = 2.00, theta = 0.695210056875677, h = 0.000200

刻み幅 0.001 で小数点以下9桁ほどの精度が出ていることがわかります。以後は，刻み幅 h = 0.0002 で続けます。

In [2]:
# 初期条件
t0 = 0.0
t1 = 2.0
# 振幅が 10°の場合
x0 = 10*pi/180
v0 = 0.0

N = 10000
h = (t1 - t0)/N

# 計算結果を配列に入れます。
# 初期条件の1個分を余分に確保。
array T10[N+1]
array X10[N+1]
array V10[N+1]

t = t0; x = x0; v = v0

do for [i=1: N+1]{
  T10[i] = t
  X10[i] = x
  V10[i] = v
  k1 = h*F1(x, v, t)
  m1 = h*F2(x, v, t)
  k2 = h*F1(x+0.5*k1, v+0.5*m1, t+0.5*h)
  m2 = h*F2(x+0.5*k1, v+0.5*m1, t+0.5*h)
  k3 = h*F1(x+0.5*k2, v+0.5*m2, t+0.5*h)
  m3 = h*F2(x+0.5*k2, v+0.5*m2, t+0.5*h)
  k4 = h*F1(x+k3, v+m3, t+h)
  m4 = h*F2(x+k3, v+m3, t+h)
  x = x + (k1 + 2.*k2 + 2.*k3 + k4)/6
  v = v + (m1 + 2.*m2 + 2.*m3 + m4)/6
  t = t + h 
}

# 初期条件
t0 = 0.0
t1 = 2.0
# 振幅が 60°の場合
x0 = 60*pi/180
v0 = 0.0

N = 10000
h = (t1 - t0)/N

# 計算結果を配列に入れます。
# 初期条件の1個分を余分に確保。
array T60[N+1]
array X60[N+1]
array V60[N+1]

t = t0; x = x0; v = v0

do for [i=1: N+1]{
  T60[i] = t
  X60[i] = x
  V60[i] = v
  k1 = h*F1(x, v, t)
  m1 = h*F2(x, v, t)
  k2 = h*F1(x+0.5*k1, v+0.5*m1, t+0.5*h)
  m2 = h*F2(x+0.5*k1, v+0.5*m1, t+0.5*h)
  k3 = h*F1(x+0.5*k2, v+0.5*m2, t+0.5*h)
  m3 = h*F2(x+0.5*k2, v+0.5*m2, t+0.5*h)
  k4 = h*F1(x+k3, v+m3, t+h)
  m4 = h*F2(x+k3, v+m3, t+h)
  x = x + (k1 + 2.*k2 + 2.*k3 + k4)/6
  v = v + (m1 + 2.*m2 + 2.*m3 + m4)/6
  t = t + h 
}

plot

振幅 10° の場合と 60° の場合の数値解をプロットします。

In [3]:
set zeroaxis
set grid
set key sample 2
plot X10 using (T10[$1]):2 w l title "θ=10°", \
     X60 using (T60[$1]):2 w l title "θ=60°"
Gnuplot Produced by GNUPLOT 5.4 patchlevel 2 -1.5 -1 -0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 θ=10° θ=10° θ=60° θ=60°

判断しやすいように，初期条件の振幅 X10[1] および X60[1] で規格化した数値解をプロットします。

In [4]:
plot X10 using (T10[$1]):($2)/(X10[1]) w l title "θ=10°", \
     X60 using (T60[$1]):($2)/(X60[1]) w l title "θ=60°"
Gnuplot Produced by GNUPLOT 5.4 patchlevel 2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 θ=10° θ=10° θ=60° θ=60°

上のグラフを見る限りでは，振幅が60°と大きくなると，周期がのびているように見えます。

もう少し定量的に調べてみます。

数値計算の答えのうち，$\dot{\theta}$ の数値解は V10[i] や V60[i] に格納されています。この数値解をみて，$\theta$ の値が増加から減少に転じる時刻，つまり $\dot{\theta}$ の値が正から負に転じる時刻を以下のようにして調べてみます。

In [5]:
print "振幅 10°の場合："
do for [i=3:N]{
  if ((V10[i]>=0) & (V10[i+1]<=0)) {print "周期 =", (T10[i]+T10[i+1])/2}
}

print "振幅 60°の場合："
do for [i=3:N]{
  if ((V60[i]>=0) & (V60[i+1]<=0)) {print "周期 =", (T60[i]+T60[i+1])/2}
}
Out[5]:
振幅 10°の場合：
周期 =1.00189999999992
振幅 60°の場合：
周期 =1.07309999999991

問題

無次元化された運動方程式を数値的に解き，振幅 $\theta_0$ を $10$° から $90$° まで $10$° きざみで大きくしていった場合，ふりこの周期はどのように変化するかを示せ。

保存量（エネルギー保存則）

無次元化された運動方程式

$$ \frac{d^2\theta}{d\tau^2} = - 4\pi^2 \sin\theta$$

の両辺に $\displaystyle \frac{d\theta}{d\tau}$ をかけて整理すると，

$$\frac{d}{d\tau}\left\{ \frac{1}{2}\left( \frac{d\theta}{d\tau} \right)^2 - 4\pi^2 \cos\theta \right\} = 0 $$

つまり，

$$ \frac{1}{2}\left( \frac{d\theta}{d\tau} \right)^2 - 4\pi^2 \cos\theta = \mbox{const.} \equiv \varepsilon $$

という保存量が存在することがわかる。（これはエネルギー保存則です。）

$t = 0$ で $\displaystyle \theta = \theta_0, \frac{d\theta}{d\tau} = 0$ という初期条件を課すと，定数 $\varepsilon$ の値は以下のように定まる。 $$ \varepsilon = - 4\pi^2 \cos\theta_0$$

つまり，

$$\frac{d\theta}{d\tau} = \pm 2\pi \sqrt{2(\cos\theta-\cos\theta_0)}$$

これは変数分離形の微分方程式であり，$\theta$ が $0$ から $\theta_0$ までの時間の4倍が周期 $\bar{T}$ になることから，

$$ \frac{1}{2\pi} \int_0^{\theta_0} \frac{d\theta}{\sqrt{2(\cos\theta-\cos\theta_0)}} = \int_0^{\bar{T}/4} d\tau $$

すなわち

$$ \bar{T} = \frac{2}{\pi} \int_0^{\theta_0} \frac{d\theta}{\sqrt{2(\cos\theta-\cos\theta_0)}} $$

この積分は解析的には解けないので数値積分すればいいと思いますが，このままの形だと，積分区間の上限で被積分関数の分母がゼロになってしまい，つまりは被積分関数が発散してしまうのでちょっと困る。なので以下のようにうまく変数変換することで，数値積分が可能な形にします。

第1種完全楕円積分

ここでは，上記の積分が変数変換によって，第1種完全楕円積分

$$K(k) \equiv \int_0^{\pi/2} \frac{dt}{\sqrt{1 - k^2 \sin^2 t}}$$

で書けることを示す。この形だと，$k^2 < 1$ なら積分区間内で被積分関数が発散することもなく，数値積分に困ることもない。

まず，分母は半角の公式を使って

\begin{eqnarray} \sqrt{2 (\cos\theta-\cos\theta_0)} &=& \sqrt{2 \left\{ \left(1-2\sin^2\frac{\theta}{2}\right) -\left(1-2\sin^2\frac{\theta_0}{2}\right) \right\} }\\ &=& 2\sqrt{\sin^2\frac{\theta_0}{2} - \sin^2\frac{\theta}{2}} \end{eqnarray}

ここで，

$$\displaystyle \sin\frac{\theta_0}{2} \equiv k, \quad \sin\frac{\theta}{2} \equiv k \sin t $$

とおくと，

\begin{eqnarray} 2\sqrt{\sin^2\frac{\theta_0}{2} - \sin^2\frac{\theta}{2}} &=& 2 k \sqrt{1 - \sin^2 t} \\ &=& 2 k \cos t \end{eqnarray}

一方，分子は

$$ \sin\frac{\theta}{2} \equiv k \sin t $$

の微分をとることにより，

$$ \cos\frac{\theta}{2} \frac{d\theta}{2} = k \cos t \, dt$$$$ \therefore \ d\theta = \frac{2 k \cos t}{\cos\frac{\theta}{2}} dt = \frac{2 k \cos t}{\sqrt{1-\sin^2\frac{\theta}{2}}} dt = \frac{2 k \cos t}{\sqrt{1-k^2 \sin^2 t}} dt $$

したがって，

\begin{eqnarray} \int_0^{\theta_0} \frac{1}{\sqrt{2(\cos\theta-\cos\theta_0)}} d\theta &=&\int_0^{\pi/2} \frac{1}{2 k \cos t} \frac{2 k \cos t}{\sqrt{1-k^2 \sin^2 t}} dt\\ &=&\int_0^{\pi/2} \frac{dt}{\sqrt{1-k^2 \sin^2 t}}\\ &\equiv& K(k) \end{eqnarray}

EllipticK(k)

ありがたいことに，gnuplot では第一種完全楕円積分 が $K(k) \equiv$ EllipticK(k) として使えます。

これを使うと，振幅 $\theta_0$ のふりこの（規格化された）周期 $\bar{T}$は

$$\bar{T} = \frac{2}{\pi} K(k), \quad k \equiv \sin\frac{\theta_0}{2}$$

となります。

In [6]:
# 第一種完全楕円積分を使って
# 振幅 60°の場合の周期を求める
x0 = 60*pi/180
k = sin(x0/2)
print 2/pi *EllipticK(k)
Out[6]:
1.07318200714936

問題

第1種完全楕円積分を使って書かれたふりこの周期の式を使い，振幅 $\theta_0$ を $10$° から $90$° まで $10$° きざみで大きくしていった場合，ふりこの周期はどのように変化するかを示せ。

また，運動方程式を数値的に解いて得られた結果と比較せよ。

参考：第1種完全楕円積分の数値積分

第1種完全楕円積分

$$K(k) \equiv \int_0^{\pi/2} \left[1 - k^2 \sin^2 t \right]^{-1/2}\, dt$$

を数値積分で求めてみる。gnuplot で数値積分する方法については，以下のリンクを参照。

シンプソン法

In [7]:
f(k,t) = 1/sqrt(1 - k**2 * sin(t)**2) # 被積分関数

# a = 0 積分の下端。実数で入力。
# b =pi/2 積分の上端。実数で入力。
# N  分割数 N は偶数とすること（シンプソン法の決め事）

# a = 0
b = pi/2
Ksimp_f(k, N, i) = (i == 2 ? \
  (b)/(N*3.0)*(f(k,(i-2)*(b)/N) + 4.0*f(k,(i-1)*(b)/N) + f(k,i*(b)/N)) : \
  Ksimp_f(k, N, i-2) + \
  (b)/(N*3.0)*(f(k,(i-2)*(b)/N) + 4.0*f(k,(i-1)*(b)/N) + f(k,i*(b)/N)))

# 実際の計算の際は，Ksimp_f(k, N, N) と同じ N を2回書く
N = 10
print N," 分割の数値解は ", Ksimp_f(0.5, N, N)
N = 20
print N," 分割の数値解は ", Ksimp_f(0.5, N, N)
N = 40
print N," 分割の数値解は ", Ksimp_f(0.5, N, N)
# 以下の結果表示から，N = 40 で十分
Out[7]:
10 分割の数値解は 1.68575035481187
20 分割の数値解は 1.6857503548126
40 分割の数値解は 1.6857503548126
In [8]:
N = 40
# 数値積分を使って
# 振幅 60°の場合の周期を求める
x0 = 60*pi/180
k = sin(x0/2)
print 2/pi *Ksimp_f(k, N, N), "   数値積分"
print 2/pi *EllipticK(k),     "   EllipticK(k)を使った場合"
Out[8]:
1.07318200714936   数値積分
1.07318200714936   EllipticK(k)を使った場合