二重振り子カオス!

こんにちは、 t-kuhn です。今日は二重振り子を Unity でシミュレーションしてみた話について書きたいと思います。

 

上の動画の二重振り子の動きがこの記事の数学を Unity で実装した結果です。Unity の物理エンジンを使わずに微分方程式を数値計算でシミュレーションするアプローチを取っています。シミュレーションに使った微分方程式は以下の二つです。

$$
\begin{align}
(m_1 + m_2) l_1 \ddot{\theta}_1 + m_2 l_2 \ddot{\theta}_2 \cos(\theta_1 - \theta_2)
\\+m_2 l_2 \dot{\theta}^2_2 \sin(\theta_1 - \theta_2) + g (m_1 + m_2) \sin \theta_1 &= 0
\tag{0.1}
\\
\\
m_2 l_2 \ddot{\theta}_2 + m_2 l_1 \ddot{\theta}_1 cos(\theta_1 - \theta_2) \\
-m_2 l_1 \dot{\theta}^2_1 sin(\theta_1 - \theta_2) + m_2 g \sin \theta_2 &= 0
\tag{0.2}
\end{align}
$$

「何も定義せずに突然こんな方程式を出されても…!」と顔をしかめてほしいです。「そもそもどこから来た微分方程式だよ!どうやって導出したの?」と疑問を持ってくれればなおさらです。のちほど説明します。

ここでまず伝えたいのは
「ある物理システムを現わしている微分方程式があれば数値計算でシミュレーションできる」
ということです。

Unity で物理システムのシミュレーションを作るときは次のようなステップを踏んでいきます。

  1. システムを分析、図形を用いて変数を定義
  2. 微分方程式を導出
  3. 二階常微分方程式を2つの一階常微分方程式に分割(Euler 法を使うための準備)
  4. Euler 法でシミュレーションする

いうまでもなく、物理システムによっては少し異なってきますが、この記事で取り上げている二重振り子の場合はこの4ステップでした。

1. システムを分析、図形を用いて変数を定義

図1

手書きで少し読みにくいところはあるが、この図があるだけで前章の式(0.1)と(0.2)がだいぶ読めるようになったのではないでしょうか。

  • $m_1$ と $m_2$ は振り子のそれぞれ質量(単位は $\mathrm{kg}$ )
  • $l_1$ と $l_2$ は質量を持たないリンクの長さ(単位は $\mathrm{m}$ )
  • $\theta_1$ は座標系のy軸と $l_1$ がなす角度(単位は $\mathrm{rad}$ )
  • $\theta_2$ は座標系のy軸と並行な線と $l_2$ がなす角度(単位は $\mathrm{rad}$ )

ここで $\theta_2$ の定義に注目してもらいたいです。 $\theta_1$ が変化するとリンク $l_1$ が回りますが、リンク $l_2$ の指す方向は一切変わりません。Unityの二重振り子モデルもそれを反映しなければならないので、 $\theta_2$ を手動でいじると下のビデオのように動かなければなりません。

 

少し細かいところではありますが、Unity 内のモデルの定義と数学的モデルの定義が合っていないとシミュレーションが予想通りに動かないから要注意です。

2. 微分方程式を導出

図を描いて何がどう繋がっていて、角度はどう測っているかを定義したから、微分方程式の導出に進められます。二重振り子の微分方程式を導き出すのは容易なことではない。進むべき道を見せてくれるガイド的な何かがないと不安なのでこの記事をガイドに数学を進めることにした。

ガイド記事を見てみると前章の式(0.1)と(0.2)がガイド記事の式(14)と(19)に該当することがわかる。ガイド記事には導出方法の概要が書かれているが、細かいところは省かれているからガイド記事の議論を追いながらできるだけ分かりやすく、式変形の途中結果も載せて説明していきたいと思います。

まず最初に振り子のそれぞれの質点 $m_1$ と $m_2$ の位置をxy座標で表すと

$$
\begin{align}
x_1 &= l_1 \sin{\theta_1} \\
y_1 &= - l_1 \cos{\theta_1} \\
x_2 &= l_1 \sin{\theta_1} + l_2 \sin{\theta_2} \\
y_2 &= - l_1 \cos{\theta_1} - l_2 \cos{\theta_2}
\tag{2.1}
\end{align}
$$
となります。この結果は Fig. 1 にしたの図にあるように直角三角形を描いて三角関数(sin と cos)の定義を適用すると求まる。

図2

次はシステムの位置エネルギー $V$ と運動エネルギー $T$ を求めておく必要がある。後ほど微分方程式の導出に使います。

システムの位置エネルギー $V$ は

$$
V = m_1 g y_1 + m_2 g y_2
\tag{2.2}
$$

である。(2.2) の $y_1$ と $y_2$ は (2.1) の方程式にも出ていますね。(2.1) の以下の二つの式を

$$
\begin{align}
y_1 &= - l_1 \cos{\theta_1} \\
y_2 &= - l_1 \cos{\theta_1} - l_2 \cos{\theta_2}
\end{align}
$$

(2.2) に代入すると

$$
\begin{align}
V &= m_1 g (-l_1 \cos{\theta_1}) + m_2 g (-l_1 \cos{\theta_1} - l_2 \cos{\theta_2})
\\
&= -m_1 g l_1 \cos{\theta_1} - m_2 g l_1 \cos{\theta_1} - m_2 g l_2 \cos{\theta_2})
\\
&= -(m_1 + m_2)g l_1 \cos(\theta_1) - m_2 g l_2 \cos{\theta_2}
\tag{2.3}
\end{align}
$$
が得られます。

次はシステムの運動エネルギー $T$ です。古典力学の運動エネルギーの定義が
$$
E_k = \frac{1}{2}mv^2
$$
であることを思い出しておきます。ここで一度私たちが考察しているシステムの質点 $m_1$ と $m_2$ のそれぞれの速度 $v_1$ と $v_2$ がどのように求まるかについて考える必要がある。角度 $\theta$ の角速度が $\dot{\theta}$ になります。そして原点から距離 $r$ 離れた質点を角速度 $\dot{\theta}$ で原点を中心に回すと質点の速度 $v$ が
$$
v = \dot{\theta} r
\tag{2.4}
$$
になるので
$$
v_1 = l_1 \dot{\theta}_1
\tag{2.5}
$$
と 質点 $m_1$ の速度 $v_1$ が求まりました。次に $v_2$ を求めましょう。(2.1) の以下の二つの方程式からスタートします。
$$
\begin{align}
x_2 &= l_1 \sin{\theta_1} + l_2 \sin{\theta_2} \\
y_2 &= - l_1 \cos{\theta_1} - l_2 \cos{\theta_2}
\tag{2.6}
\end{align}
$$

$x_2$ と $y_2$ を時間 $t$ で微分することによってそれそれの成分の速度 $v_{x_2}$ と $v_{y_2}$ を求めておきます。

$$
\begin{align}
v_{x_2} = l_1 \dot{\theta}_1  \cos{\theta_1} + l_2 \dot{\theta}_2 \cos{\theta_2}\\
v_{y_2} = l_1 \dot{\theta}_1 \sin{\theta_1} + l_2 \dot{\theta}_2 \sin{\theta_2}
\tag{2.7}
\end{align}
$$

ここで「$t$ で微分するって言っても、(2.6) に $t$ なんかどこにもないぞ!」と思う人がいるかもしれない。これは $\theta(t)$ を $\theta$ と省略してるからです。この記事の $\theta$ はみんな本当は $\theta(t)$ と思ってください。数式が見やすくなるから省略しているわけです。

(2.6) を時間 $t$ で微分するとき連鎖律
$$
\frac{d}{dt}f(g(t)) = f'(g(t)) g'(t)
$$
が適用される。

さて (2.7) のそれぞれの速度成分から $v_2$ を求めよう。ピタゴラスの定理から
$$
v_2 = \sqrt{v_{x_2}^2 + v{y_2}^2}
\tag{2.8}
$$

そして (2.8) に (2.7) を代入して整理すると

$$
\begin{align}
v_2 &= \sqrt{(l_1 \dot{\theta}_1 \cos{\theta_1} + l_2 \dot{\theta}_2 \cos{\theta_2})^2(l_1 \dot{\theta}_1 \sin{\theta_1} + l_2 \dot{\theta}_2 \sin{\theta_2})^2} \\
&= \sqrt{l_1^2 \dot{\theta}^2_1 (\cos^2{\theta_1} + \sin^2{\theta_1}) + l_2^2 \dot{\theta}^2_2 (\cos^2{\theta_2} + \sin^2{\theta_2}) + l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \cos(\theta_1 - \theta_2)} \\
&= \sqrt{l_1^2 \dot{\theta}^2_1 + l_2^2 \dot{\theta}^2_2 + l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \cos(\theta_1 - \theta_2)}
\tag{2.9}
\end{align}
$$
が得られる。これで $v_1$ と $v_2$ が求まったのでやっと運動エネルギー $T$ に戻れます。運動エネルギーの定義
$$
E_k = \frac{1}{2}mv^2
$$
から私たちが考察しているシステムの運動エネルギー $T$ は
$$
T = \frac{1}{2}m_1 v_1^2 + \frac{1}{2}m_2 v_2^2
\tag{2.10}
$$
と求まります。

(2.10) に (2.5) で求まった $v_1$ と (2.9) で求まった $v_2$ を代入すると
$$
T = \frac{1}{2}m_1 l_1^2 \dot{\theta}_1^2 + \frac{1}{2}m_2 \left(l_1^2 \dot{\theta}^2_1 + l_2^2 \dot{\theta}^2_2 + 2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \cos(\theta_1 - \theta_2)\right)
\tag{2.11}
$$
となります。

これでシステムの位置エネルギー $V$ と運動エネルギー $T$ が求まった!微分方程式を導出するために Lagrangian Mechanics の数学を使います。一般的に Newtonian Mechanics を用いて力の釣り合いと $F = ma$ の関係から微分方程式を導出することができるが、二重振り子のように、質点の軌道が制約されているシステムの場合、 Lagrangian Mechanics を使った方が都合がいい。ということで Lagrangian Mechanics の世界に突入!

Lagrangian $\mathcal{L}$ が
$$
\mathcal{L} =  T - V
\tag{2.12}
$$
と定義される。

(2.12) に (2.3) と (2.11) を代入すると
$$
\begin{align}
\mathcal{L} =&  \frac{1}{2}m_1 l_1^2 \dot{\theta}_1^2 + \frac{1}{2}m_2 \left(l_1^2 \dot{\theta}^2_1 + l_2^2 \dot{\theta}^2_2 + 2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \cos(\theta_1 - \theta_2)\right)  \\
&- \left( -(m_1 + m_2)g l_1 \cos(\theta_1) - m_2 g l_2 \cos{\theta_2} \right) \\
=& \frac{1}{2}(m_1 + m_2) l_1^2 \dot{\theta}^2_1 + \frac{1}{2}m_2l_2^2\dot{\theta}^2_2 + m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \cos(\theta_1 - \theta_2) \\
&+(m_1 + m_2)g l_1 cos{\theta_1} + m_2 g l_2 \cos{\theta_2}
\tag{2.13}
\end{align}
$$
が得られます。

次は Euler-Lagrange 方程式の登場です。この方程式は Newtonian Mechanics の有名な運動の法則の Lagrangian Mechanics 版です。$\theta_1$ についての Euler-Lagrange 方程式は

$$
\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{\theta}_1}\right) - \frac{\partial \mathcal{L}}{\partial \theta_1} = 0
\tag{2.14}
$$

となります。

もう少しで二重振り子の微分方程式が求まります!

(2.14) の計算を三つに分けるとやりやすいので

$$
\begin{align}
&\frac{\partial \mathcal{L}}{\partial \dot{\theta}_1} \tag{2.15}\\
&\frac{d}{dt}\left( \frac{\partial \mathcal{L}}{\partial \dot{\theta}_1} \right)  \tag{2.16}\\
&\frac{\partial \mathcal{L}}{\partial \theta_1} \tag{2.17}\\
\end{align}
$$
と計算量を分割しよう。

まず (2.15) に (2.13) を代入して偏微分を計算すると
$$
\begin{align}
\frac{\partial \mathcal{L}}{\partial \dot{\theta}_1} =&
\frac{\partial}{\partial \dot{\theta}_1} \bigg(\frac{1}{2}(m_1 + m_2) l_1^2 \dot{\theta}^2_1 + \frac{1}{2}m_2l_2^2\dot{\theta}^2_2 + m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \cos(\theta_1 - \theta_2) \\
&+(m_1 + m_2)g l_1 cos{\theta_1} + m_2 g l_2 \cos{\theta_2}\bigg) \\
=& (m_1 + m_2)l_1^2 \dot{\theta}_1 + m_2 l_1 l_2 \dot{\theta}_2 \cos(\theta_1 - \theta_2) \tag{2.18}
\end{align}
$$
が得られます。

次に (2.16) に (2.18) を代入して微分を実行すると
$$
\begin{align}
\frac{d}{dt}\left( \frac{\partial \mathcal{L}}{\partial \dot{\theta}_1} \right)  =& \frac{d}{dt}\left( (m_1 + m_2)l_1^2 \dot{\theta}_1 + m_2 l_1 l_2 \dot{\theta}_2 \cos(\theta_1 - \theta_2)  \right) \tag{2.19}\\
=& (m_1 + m_2) l_1^2 \ddot{\theta}_1 + m_2 l_1 l_2 \ddot{\theta}_2 \cos(\theta_1 - \theta_2) \\
&- m_2 l_1 l_2 \dot{\theta}_2 \sin(\theta_1 - \theta_2) (\dot{\theta}_1 - \dot{\theta}_2) \tag{2.20}
\end{align}
$$
になります。

(2.19) の微分を実行するときに連鎖律
$$
\frac{d}{dt}f(g(t)) = f'(g(t)) g'(t)
$$
と積の法則
$$
\frac{d}{dt}f(t)g(t) = f'(t)g(t) + f(t)g'(t)
$$
が適用される。

最後に (2.17) に (2.13) を代入して偏微分を実行すると
$$
\begin{align}
\frac{\partial \mathcal{L}}{\partial \theta_1} =&  \frac{\partial}{\partial \theta_1}
\bigg(\frac{1}{2}(m_1 + m_2) l_1^2 \dot{\theta}^2_1 + \frac{1}{2}m_2l_2^2\dot{\theta}^2_2 + m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \cos(\theta_1 - \theta_2) \\
&+(m_1 + m_2)g l_1 cos{\theta_1} + m_2 g l_2 \cos{\theta_2}\bigg) \\
=& -m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \sin(\theta_1 - \theta_2) - (m_1 + m_2) g l_1 \sin{\theta_1} \tag{2.21}
\end{align}
$$
が得られます。

これで (2.14) の準備のための計算が終わりました。 (2.14) をもう一度書いておきます。
$$
\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{\theta}_1}\right) - \frac{\partial \mathcal{L}}{\partial \theta_1} = 0
\tag{2.14}
$$

(2.14) に (2.20) と (2.21) を代入して整理すると

$$
\begin{align}
&(m_1 + m_2) l_1^2 \ddot{\theta}_1 + m_2 l_1 l_2 \ddot{\theta}_2 \cos(\theta_1 - \theta_2)m_2 l_1 l_2 \dot{\theta}_2 \sin(\theta_1 - \theta_2) (\dot{\theta}_1 - \dot{\theta}_2) \\
&- \left(-m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \sin(\theta_1 - \theta_2) - (m_1 + m_2) g l_1 \sin{\theta_1}\right) \\
\\
=& (m_1 + m_2) l_1^2 \ddot{\theta}_1 + m_2 l_1 l_2 \ddot{\theta}_2 \cos(\theta_1 - \theta_2) - m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \sin(\theta_1 - \theta_2) \\
& m_2 l_1 l_2 \dot{\theta}^2_2\sin(\theta_1 - \theta_2) + m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \sin(\theta_1 - \theta_2) + (m_1 + m_2) g l_1 \sin{\theta_1} \\
\\
=& (m_1 + m_2) l_1^2 \ddot{\theta}_1 + m_2 l_1 l_2 \ddot{\theta}_2 \cos(\theta_1 - \theta_2) + m_2 l_1 l_2 \dot{\theta}^2_2 \sin(\theta_1 - \theta_2) \\
&+ (m_1 + m_2) g l_1 \sin{\theta_1} = 0 \tag{2.22}
\end{align}
$$
が得られます。

最後に (2.22) の両辺を $l_1$ で割ると
$$
\begin{align}
(m_1 + m_2) l_1 \ddot{\theta}_1 + m_2 l_2 \ddot{\theta}_2 \cos(\theta_1 - \theta_2)
\\+m_2 l_2 \dot{\theta}^2_2 \sin(\theta_1 - \theta_2) + g (m_1 + m_2) \sin \theta_1 &= 0
\tag{0.1}
\end{align}
$$
と、$\theta_1$ についての微分方程式が得られました!

同様に、$\theta_2$ における Euler-Lagrange 方程式
$$
\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{\theta}_2}\right) - \frac{\partial \mathcal{L}}{\partial \theta_2} = 0
\tag{2.23}
$$
の微分と偏微分を実行すると
$$
\begin{align}
m_2 l_2 \ddot{\theta}_2 + m_2 l_1 \ddot{\theta}_1 cos(\theta_1 - \theta_2) \\
-m_2 l_1 \dot{\theta}^2_1 sin(\theta_1 - \theta_2) + m_2 g \sin \theta_2 &= 0
\tag{0.2}
\end{align}
$$
と $\theta_2$ についての微分方程式が求まります。

3. Euler 法を使うための準備

Euler 法で前章で得られた二重振り子の微分方程式のシミュレーションを行いたいが、実は (0.1) と (0.2) をそのまま使えないのです。Euler 法で一階常微分方程式を数値計算的にシミュレーションできるが、(0.1) と (0.2) はどちらも二階常微分方程式。

Euler 法が使えるように (0.1) と (0.2) それぞれの二階常微分方程式を二つの一階常微分方程式に分ければいい。どうすれば分割できる?と疑問に思う人もいるかもしれません。新しい変数を次のように定義すればよい。
$$
\begin{align}
\lambda_1 &= \dot{\theta}_1 \tag{3.1} \\
\dot{\lambda}_1 &= \ddot{\theta}_1 \tag{3.2}
\\
\lambda_2 &= \dot{\theta}_2 \tag{3.3} \\
\dot{\lambda}_2 &= \ddot{\theta}_2 \tag{3.4}
\end{align}
$$

(3.1) ~ (3.4) を (0.1) に代入すると

$$
\begin{align}
(m_1 + m_2) l_1 \dot{\lambda}_1 + m_2 l_2 \dot{\lambda}_2 \cos(\theta_1 - \theta_2) \\
+ m_2 l_2 \lambda^2_2 \sin(\theta_1 - \theta_2) + g (m_1 + m_2) \sin \theta_1 &= 0
\tag{3.5}
\end{align}
$$

が得られます。

次に (3.5) を $\dot{\lambda}_1$ について解くと

$$
\begin{align}
\dot{\lambda}_1 &= \frac{- m_2 l_2 \dot{\lambda}_2 \cos(\theta_1 - \theta_2) - m_2 l_2 \lambda^2_2 \sin(\theta_1 - \theta_2) - g (m_1 + m_2) \sin \theta_1}{(m_1 + m_2) l_1}
\tag{3.6}
\end{align}
$$

となります。

同様に (3.1) ~ (3.4) を (0.2) に代入して

$$
\begin{align}
m_2 l_2 \dot{\lambda}_2 + m_2 l_1 \dot{\lambda}_1 cos(\theta_1 - \theta_2)
- m_2 l_1 \lambda^2_1 sin(\theta_1 - \theta_2) + m_2 g \sin \theta_2 &= 0
\tag{3.7}
\end{align}
$$

(3.7) を $\dot{\lambda}_2$ について解くと

$$
\begin{align}
\dot{\lambda}_2 &= \frac{m_2 l_1 \lambda^2_1 sin(\theta_1 - \theta_2) - m_2 l_1 \dot{\lambda}_1 cos(\theta_1 - \theta_2) - m_2 g \sin \theta_2}{m_2 l_2}
\tag{3.8}
\end{align}
$$
が得られます。

これでそれぞれの二階常微分方程式をそれぞれ二つの一階常微分方程式に分割することに成功しました。分割後の四つの一階常微分方程式をもう一度まとめて書いておくと
$$
\begin{align}
\dot{\theta}_1  &= \lambda_1 \tag{3.1} \\
\dot{\lambda}_1 &= \frac{- m_2 l_2 \dot{\lambda}_2 \cos(\theta_1 - \theta_2) - m_2 l_2 \lambda^2_2 \sin(\theta_1 - \theta_2) - g (m_1 + m_2) \sin \theta_1}{(m_1 + m_2) l_1}
\tag{3.6} \\
\dot{\theta}_2 &= \lambda_2 \tag{3.3} \\
\dot{\lambda}_2 &= \frac{m_2 l_1 \lambda^2_1 sin(\theta_1 - \theta_2) - m_2 l_1 \dot{\lambda}_1 cos(\theta_1 - \theta_2) - m_2 g \sin \theta_2}{m_2 l_2}
\tag{3.8} \\
\end{align}
$$
となります。

この四つの方程式に Newton 記法の微分点が二つ付いている変数が一つもないことに注意してほしい。みんなが一階常微分方程式になっているのでこれらに Euler 法を適用できます!

4. Euler 法でシミュレーションする

やっと Euler 法を適用するところまで来ました。まず Euler 法の仕組みを見てみましょう。
$$
\begin{align}
\dot{y} &= f(y, t) \tag{4.1}\\
y(t_0) &= y_0 \tag{4.2}
\end{align}
$$
(4.1) で表せるような一階常微分方程式があるとしよう。そして (4.2) でこの微分方程式の初期状態を定義しよう。このときに、Euler 法で時間を $\epsilon$ だけ進めると
$$
\begin{align}
t_1 &= t_0 + \epsilon \\
y_1 &= y_0 + \epsilon f(y_0, t_0) \tag{4.3}
\end{align}
$$
$y_1$ が定まります。ここで、 $\epsilon$ はシミュレーションのタイムステップで、小さければ小さいほどシミュレーションがより正確になる。今回の二重振り子シミュレーションでは「0.0001」と設定しました。(4.3) はシミュレーションを開始させて一番最初のステップを実行させたときになりますが、その後も同じパターンで進みます。

ここで重要なのが、Euler 法が $\dot{y}_n = f(y_n, t_n)$ を入力として与えられて $y_{n+1}$ を返してくれるところです!つまり、 $\dot{y}$ から $y$ が得られる。

もう一度準備した微分方程式を見てみましょう。

$$
\begin{align}
\dot{\theta}_1  &= \lambda_1 \tag{3.1} \\
\dot{\lambda}_1 &= \frac{- m_2 l_2 \dot{\lambda}_2 \cos(\theta_1 - \theta_2) - m_2 l_2 \lambda^2_2 \sin(\theta_1 - \theta_2) - g (m_1 + m_2) \sin \theta_1}{(m_1 + m_2) l_1}
\tag{3.6} \\
\dot{\theta}_2 &= \lambda_2 \tag{3.3} \\
\dot{\lambda}_2 &= \frac{m_2 l_1 \lambda^2_1 sin(\theta_1 - \theta_2) - m_2 l_1 \dot{\lambda}_1 cos(\theta_1 - \theta_2) - m_2 g \sin \theta_2}{m_2 l_2}
\tag{3.8} \\
\end{align}
$$

$\dot{\lambda}_1$ を Euler 法で処理して $\lambda_1$ が得られます。さらに $\dot{\theta}_1  = \lambda_1 $ の関係を用いて Euler 法で $\dot{\theta}_1$ を処理して $\theta_1$ が得られます。同様に $\theta_2$ も得られるので $\theta_1$ と $\theta_2$ を毎フレーム Unity の中で作っておいた二重振り子に適用するとシミュレーション完成です!

 

ちなみに、「そもそもシミュレーションなんかしないで微分方程式を解けばいいんじゃない?」と疑問に思う人もいるかもしれませんが、二重振り子の微分方程式は解けません。上のビデオを見るとその原因を垣間見ることができます。二重振り子にはカオスの種が埋もれていて、初期状態がほぼ同じな二重振り子でも時間が経つとふるまいが大きく異なってきます。

5. 最後に

二重振り子プロジェクトを Github で公開しています: link

Q:「Euler 法より効率的な数値計算法いっぱいあるんじゃない?何で Euler 法にしたんだろう?」
A: 「ブログ記事で一番説明が簡単になるのが Euler 法だったからです。上記のリンクの振り子シミュレーションプロジェクトには Euler 法以外にも Runge-Kutta (RK4) を実装したので、いろいろ比較できます。」

リンク集: