My Blog Logo
Blog
数値シミュレーション －二重振り子－

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

数値シミュレーション解説 第7回

カテゴリ:

シミュレーション

タグ:

シミュレーション数学

公開日: 8/19/2018

更新日: 12/23/2025

読了目安:

5分

数値シミュレーション解説の前回の記事

数値シミュレーション －各解法の誤差比較－

数値シミュレーション解説の次回の記事

数値シミュレーション －ローレンツ・アトラクタ－

このシリーズの記事一覧（全8回）
  1. 第1回 数値シミュレーション －数理モデル－
  2. 第2回 数値シミュレーション －微分方程式－
  3. 第3回 数値シミュレーション －オイラー法－
  4. 第4回 数値シミュレーション －ホイン法－
  5. 第5回 数値シミュレーション －ルンゲ＝クッタ法－
  6. 第6回 数値シミュレーション －各解法の誤差比較－
  7. 第7回 数値シミュレーション －二重振り子－
  8. 第8回 数値シミュレーション －ローレンツ・アトラクタ－
更新履歴
  • 2025/12/23: プログラムを C# 14.0 に修正
こんにちは、Rafka です。
今回は手計算では厳密解が求められないような微分方程式で数理モデルが表される物理現象を、 ルンゲ＝クッタ法を用いてシミュレーションを行っていきます。

二重振り子

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

2重振り子

www.aihara.co.jp のファビコン
www.aihara.co.jp
（新しいタブで開く）
θ˙1=ϑ1θ˙2=ϑ2ϑ˙1=g(sin⁡θ2C−m1,2m2sin⁡θ1)−(l1ϑ12C+l2ϑ22)Sl1(m1,2m2−C2)ϑ˙2=gm1,2m2(Csin⁡θ1−sin⁡θ2)+(m1,2m2l1ϑ12+l2ϑ22C)Sl2(m1,2m2−C2)\begin{split} \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{split}
各記号は以下の通りです。
  • m1m_1 : 1 つ目の振り子の質量
  • l1l_1 : 1 つ目の振り子の長さ
  • θ1\theta_1 : 1 つ目の振り子の角度
  • ϑ1\vartheta_1 : 1 つ目の振り子の角速度
  • m2m_2 : 2 つ目の振り子の質量
  • l2l_2 : 2 つ目の振り子の長さ
  • θ2\theta_2 : 2 つ目の振り子の角度
  • ϑ2\vartheta_2 : 2 つ目の振り子の角速度
  • m1,2m_{1, 2} : m1m_1 と m2m_2 の和
また、ドットはそれぞれ記号の時間微分を表しています。 最後に、SS と CC はそれぞれ以下の値です。
S=sin⁡(θ1−θ2)C=cos⁡(θ1−θ2)\begin{split} S &= \sin(\theta_1 - \theta_2) \\ C &= \cos(\theta_1 - \theta_2) \end{split}
数理モデルが分かったので、実際にプログラムしてシミュレーションしてみましょう。

シミュレーション結果

各振り子の重さをそれぞれ 1.01.0、長さも 1.01.0 とし、 初期の角度を θ1=90∘,θ2=90∘\theta_1 = 90^\circ, \theta_2 = 90^\circ、 角速度はどちらも 00 として、 時間刻み 1/29.971 / 29.97s で約 6060 秒間のシミュレーションを行いました。 結果は各振り子の角度で記録されるので、 その値から各質点の位置に変換し、 さらに角度が 0∘0^\circ の時に真下を向くように −90∘-90^\circ 回転させた時の位置の変化を グラフで描画してみました。
前回の斜方投射のグラフと同じように、 時間による変化をアニメーションさせることができるので、 ボタンやスライドバーを操作して遊んでみてください。
因みに二重振り子のシミュレーションは、 カオスと呼ばれる系の代表例でして、 初期値鋭敏性という性質から予測できない軌道を示すことが知られています。
二重振り子のシミュレーションも、 初期値をほんの少し変更するだけで、 上のグラフの軌道とは全く異なる軌道が生成されるので、 是非今まで示してきたプログラムを組み合わせて、 自分の手元で初期値を変えながらシミュレーションを行ってみてください。

まとめ

二重振り子のシミュレーションを通して、 厳密解が求められないような微分方程式についても数値解を計算し、 その答えがもっともらしい事を確認できたと思います。
次回からはついに、脳のシミュレーションの話に入っていきます。 手始めに次回は、 脳を構成するパーツである神経細胞（ニューロン）の性質と、 数理モデルについて解説していきます。
Rafka

Rafka

Software Engineer

Rafka って言います。品川でソフトウェアエンジニアやってます。

Limbus Company の攻略と感想、中央競馬 G1 レースの予想、参加したイベントの感想などを書いています。

大学では計算論的神経科学や高性能計算を研究していたので、数値シミュレーションの解説も書いています。

Ace Combat 8 を買いました！シリーズファンなのでめっちゃ楽しんでいます！

数値シミュレーション解説の前回の記事

数値シミュレーション －各解法の誤差比較－

数値シミュレーション解説の次回の記事

数値シミュレーション －ローレンツ・アトラクタ－

前回の記事

数値シミュレーション －各解法の誤差比較－

次回の記事

新しい記事の事とか

関連記事

数値シミュレーション －ローレンツ・アトラクタ－
数値シミュレーション －ローレンツ・アトラクタ－

シミュレーション

約4分

数値シミュレーション解説 第8回

シミュレーション
数学
カオス

カオスの代表例であるローレンツ・アトラクタについて解説し、3D ビジュアライゼーションで初期値鋭敏性を体験します。

数値シミュレーション －各解法の誤差比較－
数値シミュレーション －各解法の誤差比較－

シミュレーション

約10分

数値シミュレーション解説 第6回

シミュレーション
数学

今まで紹介してきた数値解法の誤差比較を、斜方投射を題材に行っていきます。

数値シミュレーション －ルンゲ＝クッタ法－
数値シミュレーション －ルンゲ＝クッタ法－

シミュレーション

約6分

数値シミュレーション解説 第5回

シミュレーション
数学

常微分方程式の数値解法の 1 つであるルンゲ＝クッタ法の導出と実装を行います。

数値シミュレーション －ホイン法－
数値シミュレーション －ホイン法－

シミュレーション

約11分

数値シミュレーション解説 第4回

シミュレーション
数学

常微分方程式の数値解法の 1 つであるホイン法の導出と実装を行います。

数値シミュレーション －オイラー法－
数値シミュレーション －オイラー法－

シミュレーション

約8分

数値シミュレーション解説 第3回

シミュレーション
数学

常微分方程式の数値解法の 1 つであるオイラー法の導出と実装を行います。