非減衰単振り子の運動（Part 2 -シミュレーション編-）

非減衰単振り子の運動（Part 2 -シミュレーション編-）

はじめに

前回は、非減衰単振り子の運動方程式を導出しました。
本稿ではその続きとして、非線形モデルと線形化モデルを 4 次のルンゲ–クッタ法（RK4）を用いて数値解析し、その挙動を比較していきます。
前回記事はこちら: 非減衰単振り子の運動（Part 1 - シミュレーション準備編 -）

---

単振り子のモデル

1. 概要

Fig.1 Schematic of a simple pendulum

Fig.1に示すように、単振り子の構成要素は以下です。

• 質点 $m$（糸の先端）
• 糸の長さ $l$
• 重力加速度 $g$

仮定：

1. 糸は質量がなく、伸縮しない。
2. 空気抵抗や軸の摩擦はない。
3. 振れ角は十分小さい場合、小振幅近似が可能。

2. 運動方程式

非線形モデル：

$$\ddot{\theta} + \frac{g}{l} \sin\theta = 0 \tag{1}$$

線形化モデル：

$$\ddot{\theta} + \frac{g}{l} \theta = 0 \tag{2}$$

---

4 次のルンゲ・クッタ法（RK4）の公式

初期値問題の解法として4 次のルンゲ・クッタ法(単にルンゲ・クッタとも呼ぶ)がよく用いられている。この方法は$f(x,y)$を 4 回計算する必要があるが、計算精度が良いことで知られている。次の公式(3)で計算する。

$$\begin{aligned} y_{n+1} = y_n + \frac{h}{6}(k_1+2k_2+2k_3+k_4) \end{aligned} \tag{3}$$

ただし、

$$\begin{aligned} k_1 &= f\left(x_n, y_n \right) \\ k_2 &= f\left( x_n + \frac{h}{2}, y_n + \frac{h}{2}k_1 \right) \\ k_3 &= f\left( x_n + \frac{h}{2}, y_n + \frac{h}{2}k_2 \right) \\ k_4 &= f\left( x_n + h, y_n + hk_3 \right) \\ \end{aligned} \tag{4}$$

である。

ルンゲ–クッタ法の詳しい解説は下のサイトが参考になります。
Python での常微分方程式の数値解法（Runge–Kutta 法）

---

シミュレーション１

まずは、初期角度 $\theta_0$ を変化させたときに、非線形モデルと線形モデルで振り子の角度 $\theta(t)$ がどのように変化するかを比較してみます。
計算には Octave を用い、パラメータは次の通りです。

simulation_1.m

% simulation parameters
Tsim = 10;   % シミュレーション時間 [s]
dt   = 0.01; % 刻み幅 [s]

% initial condition
theta_o  = [0.03, 0.5, 0.9];  % [rad] 初期角度
dtheta_o = 0.0; % [rad/s] 初期角速度

% parameters
g = 9.81;  % [m/s^2] 重力加速度
l = 0.5;   % [m] 糸の長さ

計算結果

Fig.2 Simulation results for three different initial angles.

Fig.2に計算結果を示します。初期角度が大きくなるにつれて、非線形モデルの周期が長くなっているのが確認できます。

$\theta_0=0.03$ [rad]では両方のモデルが良く一致していますが、 $\theta_0=0.5, 0.9$[rad]で周期の差がはっきり確認できますね。これは、線形モデルの周期は一定であるのに対し、非線形モデルでは振幅に依存しているためと考えられます。

振幅依存の周期（非線形モデル）

非線形モデルの振動周期は次式(5)で計算できます。

$$T_{\text{non-linear}} = \frac{4K(k)}{\omega} \quad \left( k = \sin \frac{\alpha}{2}, \; \omega = \sqrt{\frac{g}{l}} \right) \tag{5}$$

ここで $\alpha$ は振り子の最大角を表します。摩擦がない今回の条件では、この最大角 $\alpha$ は初期条件で与えた角度 $\theta_0$ に対応します。右辺に現れる積分は第 1 種完全楕円積分 $K(k)$ です。

したがって、非線形モデルの周期は $\omega$ のみならず、最大角 $\alpha (=\theta_0)$にも依存することが分かります。そのためシミュレーション結果では、初期角度 $\theta_0$ が大きいほど振動周期が長くなる様子が確認できます。

一方、線形モデルの周期は次式 (6) のように表され、$\omega$ のみに依存するため、非線形モデルとは異なる結果を示します。

$$T_{\text{linear}} = 2\pi \sqrt{\frac{l}{g}} \tag{6}$$

式(5)の導出は以下のページにわかりやすく解説されています。

🔗 新潟大学 工学部「単振り子の運動とヤコビの楕円関数」

運動シミュレーション

次に、それぞれの初期角度について運動をシミュレーションします。
青い球が線形モデル、グレーの球が非線形モデルを表しており、両者の運動の違いを確認できます。
Fig.5 ($\theta_0 = 0.03$ rad) では、線形モデルだけ確認できますが、実際は両モデルの運動がほぼ重なっています。

Fig.5 Pendulum motion (initial angle $\theta_0 = 0.03$ rad). The blue ball represents the linear model and the gray ball represents the nonlinear model.

Fig.6 Pendulum motion (initial angle $\theta_0 = 0.50$ rad). The blue ball represents the linear model and the gray ball represents the nonlinear model.

Fig.7 Pendulum motion (initial angle $\theta_0 = 0.90$ rad). The blue ball represents the linear model and the gray ball represents the nonlinear model.

ポイント

線形モデルは本来、振れ角が微小であることを前提としています。
そのため、初期角度が大きい場合には非線形モデルの挙動に追従できなくなります。
したがって「近似を行う」ということは、その適用範囲や意味を十分に理解しておく必要があります。

---

シミュレーション ２ : ブランコをイメージしよう!

次に、振り子の長さを 2 m に設定してみます。
単なる興味本位です。
このときは初期角度 $\theta_0$ に加えて、初速度 $\dot{\theta}_0$ も与えます。
これはちょうど ブランコに乗って漕ぎ始める状況をイメージすると分かりやすいでしょう。

simulation_2.m

% initial condition
theta_o  = 0.5;  % [rad] 初期角度
dtheta_o = -3.0; % [rad/s] 初期角速度

% parameters
l = 2.0;   % [m] ブランコの長さ

計算結果

Fig.8 Pendulum angle response with initial velocity ($\theta_0 = 0.5$ rad, $\dot{\theta}_0 = -3.0$ rad/s, $l = 2.0$ m).

Fig.8 に計算結果を示します。シミュレーション 1 と同様に、非線形モデルでは周期が長くなる様子が確認できます。さらに今回は初速度を与えているため、振幅が線形モデルよりも大きくなっている点にも注目してください。

運動シミュレーション

Fig.9 Pendulum motion ($\theta_0 = 0.5$ rad, $\dot{\theta}_0 = -3.0$ rad/s, $l = 2.0$ m).

長さを 2 m にしたことで、振り子はよりゆったりと揺れるようになります。
直感的にも、ブランコをこいでいると上の位置で滞空している時間が長く感じられるとおもいます。

※ 実際には摩擦などの影響によって振動は減衰していく点に注意してください。（あとくつ飛ばしもできます）

---

シミュレーション ３ : ブランコで一回転!?（おまけ）

さらに初速度を大きくするとどうなるでしょうか。
実際のブランコであれば、先ほどの結果でも十分にスリルのある体験になると思いますが、
シミュレーションの良いところは、その先を行けるところです。

simulation_3.m

% initial condition
theta_o  = 0.5;  % [rad] 初期角度
dtheta_o = -5.0; % [rad/s] 初期角速度（後ろに蹴り出すイメージ）

% parameters
l = 2.0;   % [m] ブランコの長さ

計算結果

Fig.10に示すように、非線形モデルが振り切っていますね。これはFig.11のアニメーションを見ると何が起きているか確認できます。

Fig.10 Pendulum angle response with initial velocity ($\theta_0 = 0.5$ rad, $\dot{\theta}_0 = -5.0$ rad/s, $l = 2.0$ m).

Fig.11 Pendulum motion ($\theta_0 = 0.5$ rad, $\dot{\theta}_0 = -5.0$ rad/s, $l = 2.0$ m).

そうです。初速を大きくしたので、ブランコ（非線形）が何回転もしています。実際は摩擦などの影響で減衰していきますが、乗る人が頑張ればできるかもしれません...

まとめ

線形化するといろいろ便利なことが多いですが、その意味と適用範囲についてよく知っておかなければなりません。
また、多くの力学問題では解析解を得ることは難しいため、今回のように 数値計算（RK4 など）を用いたシミュレーションが強力な手段となります。