シミュレーション 数値計算

数値シミュレーション -二重振り子-

こんにちは、Rafkaです。

前回 → 数値シミュレーション -各解法の誤差比較-

今回は手計算では厳密解が求められないような微分方程式で数理モデルが表される物理現象を、ルンゲ=クッタ法を用いてシミュレーションを行っていきます。

 

二重振り子


厳密解が求められないような数理モデルというわけで、題材には二重振り子を選びました。二重振り子とは、振り子の先にもう一つ振り子が付いている振り子のことです。

二重振り子の数理モデルは以下のサイトを参考にさせていただきました。

二重振り子|MECHANICS OF CHAOS MAN — カオス人形の仕組み —

$$\begin{align}\dot{\theta}_1 &= \vartheta_1\\ \dot{\theta}_2 &= \vartheta_2\\ \dot{\vartheta}_1 &= \frac{g(\sin\theta_2C – \frac{m_{1,2}}{m_2}\sin\theta_1) – (l_1\vartheta_1^2C + l_2\vartheta_2^2)S}{l_1(\frac{m_{1,2}}{m_2} – C^2)}\\ \dot{\vartheta}_2 &= \frac{g\frac{m_{1,2}}{m_2}(C\sin\theta_1 – \sin\theta_2) + (\frac{m_{1,2}}{m_2}l_1\vartheta_1^2 + l_2\vartheta_2^2C)S}{l_2(\frac{m_{1,2}}{m_2} – C^2)}\end{align}$$

各記号は、\(m_1\) は1つ目の振り子の質量、\(l_1\) は1つ目の振り子の長さ、\(\theta_1\) は1つ目の振り子の角度、\(\vartheta_1\) は1つ目の振り子の角速度、\(m_2\) は2つ目の振り子の質量、\(l_2\) は2つ目の振り子の長さ、\(\theta_2\) は2つ目の振り子の角度、\(\vartheta_2\) は2つ目の振り子の角速度、\(m_{1, 2}\) は \(m_1\) と \(m_2\) の和、ドットはそれぞれ記号の時間微分を表しています。また、\(S\) と \(C\) はそれぞれ以下の値です。

$$\begin{align}S &= \sin(\theta_1 – \theta_2)\\ C &= \cos(\theta_1 – \theta_2)\end{align}$$

数理モデルが分かったので、実際にプログラムしてシミュレーションしてみましょう。

 

シミュレーション結果


本記事シリーズの公開リポジトリに、二重振り子のシミュレーションを行うプロジェクトを、DoublePendulumという名前で追加・実装を行いました。

二重振り子のシミュレーションのために作成したプログラムが気になる方はこちらから確認できます。

各振り子の重さをそれぞれ \(1.0\)、長さも \(1.0\) とし、初期の角度を \(\theta_1 = 90^\circ,\;\theta_2 = 90^\circ\)、角速度はどちらも \(0\) として、時間刻み \(0.1\) で \(60\) 秒間のシミュレーションを行いました。結果は各振り子の角度で記録されるので、その値から各質点の位置に変換し、さらに角度が \(0^\circ\) の時に真下を向くように \(-90^\circ\) 回転させた時の位置の変化を以下のグラフで表します。

前回の斜方投射のグラフと同じように、時間による変化をアニメーションさせることができるので、ボタンやスライドバーを操作して遊んでみてください。

因みに二重振り子のシミュレーションは、カオスと呼ばれる系の代表例でして、初期値鋭敏性という性質から予測できない軌道を示すことが知られています。二重振り子のシミュレーションも、初期値をほんの少し変更するだけで上のグラフの軌道とは全く異なる軌道が生成されるので、是非今まで示してきたプログラムを組み合わせて、自分の手元で初期値を変えながらシミュレーションを行ってみてください。

もしかしたら後日、初期値を変更した時のプロットを数パターン追加するかもしれません。

 

まとめ


二重振り子のシミュレーションを通して、厳密解が求められないような微分方程式についても数値解を計算し、その答えがもっともらしい事を確認できたと思います。

次回からはついに脳のシミュレーションの話に入っていきます。手始めに次回は、脳を構成するパーツである神経細胞(ニューロン)の性質と、数理モデルについて解説していきます。

次回 →