Return to コンピュータ演習

単振り子の運動方程式を数値的に解く準備

単振り子の振幅と周期の関係を数値的に調べる。振幅が「十分に小さい」時には単振動になり,単振り子の周期は振幅に依存しないことが知られている。これを

振り子の等時性

という。

では,振幅が大きい場合(「十分に小さい」という近似が成り立たない場合),単振り子の周期は振幅にどのように依存するか。

問題(3択)

まずは予想を立ててみよう。

単振り子の振幅が大きくなると(「十分に小さい」という近似が成り立たなくなると),単振り子の周期は,単振動の場合の周期にくらべて…

  1. 「振り子の等時性」があるから変化しない,一定のままである。
  2. 振幅が大きくなればなるほど,周期は長くなる。
  3. 振幅が大きくなればなるほど,周期は短くなる。

 

単振り子の運動方程式の導出

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

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

となることを示す。

運動方程式

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

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

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

保存量(角運動量の鉛直成分)

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

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

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

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

この式は,角運動量 $\displaystyle \boldsymbol{L} \equiv m \boldsymbol{r}\times\frac{d\boldsymbol{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
$$

極座標
$$\boldsymbol{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$ 平面内に限られることがわかる。

最終的な運動方程式

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

すると,$\boldsymbol{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} = – \frac{g}{\ell} \theta$$

となり,これは単振動となる。

単振動の周期

単振動の周期 $\tau_0$ は

$$ \tau_0 = 2\pi \sqrt{\frac{\ell}{g}} $$

となり,これは初期条件で与えられる振れ幅(振幅)$\theta_0$ に依存しない。

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

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

運動方程式

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

は解析的には解けない。

運動方程式の無次元化

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

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

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

となる。また,規格化された時間 $T$ を

$$ T \equiv \frac{t}{\tau_0}$$

のように定義すると

$$ \tau_0^2 \frac{d^2}{dt^2} = \frac{d^2}{dT^2}$$

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

$$ \frac{d^2\theta}{dT^2} = – 4\pi^2 \sin\theta$$

Runge-Kutta 法で数値的に解くための準備

この2階常微分方程式を数値的に解くわけだが,常微分方程式の数値解法である Runge-Kutta 法(ルンゲ=クッタ法)を使うために,2階常微分方程式を連立の1階常微分方程式系にしておく。(昔,若かりし頃に必要に迫られて2階常微分方程式の数値解法を調べることになったが,教科書には1階常微分方程式の数値解法として Runge-Kutta 法が紹介されているのみで,「2階」常微分方程式の数値解法については一切記載がなく,「私はいったいどうすればいいのだ… 」と呆然としてしまったことがあるので,くどいようですが念のため。)

2階常微分方程式を Runge-Kutta 法で数値的に解くためには,連立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}

とし,$T = T_0$ から $ T_1$ まで解く。

Python で単振り子の運動方程式を解いて振幅と周期の関係を調べる

Python を使って運動方程式を数値的に解き,単振り子の振幅と周期の関係を調べる。導出については,以下を参照。

  • 「単振り子の運動方程式を数値的に解く準備」