久々に数値解析をやってみたくなりました。題材としては多くの人がやっている二重振り子が面白そうです。では言語は何にしようかと迷うところですが、コードを書いて計算して結果を描画するまで可能な「Octave」が良さそうです。
いざやってみようとすると結構苦労しました。難しい点をいくつかの関門に見立てて以下の記事にまとめました。本記事には特に有益な情報はありませんが、結果をYouTube動画にアップして説明欄にソースコードを載せております。興味がありましたらご参照ください。(描画の部分は割愛しております。)
第一関門：ラグランジュ方程式
二重振り子の解析のソースコードはあらゆる言語でWeb上に情報があります。しかし大まかな流れを理解していないとコピー＆ペーストだけでは解析実行は難しいです。ラグランジュ方程式による解法を学んでおくとさらなるN重振り子への展開もスムーズにできそうです。ラグランジュ方程式の導出についてもWeb上に多くの情報があるのでそちらに任せます。
とりあえず私はWeb上の資料を読んで、最終式の添字等が間違っていないかを確認しました。経験上、他人の情報でコードを書く際、i,j,kなどの添字の間違いで正しい結果が出ずにハマる場合が多々あります。なので最低限の確認をしておいた方が良いです。
第二関門：オイラー法（ホイン法）
ステップごとに計算していく方法としてオイラー法が一番簡単なようです。しかし誤差が大きい？というのが一般的な理解のようです。またアニメーションは長い時間の方が面白いので、できるだけ時間ステップ幅は大きく取りたいところです。そうすると２次的な近似のホイン法が良さそうです。
（さらに４次のルンゲクッタもありますが私はホイン法でやる事にしました。）
第三関門：高速描画
アニメーションを作る上で描画をある程度高速に行えないと作業効率が悪く、凝った動画を作成する事ができません。
Octaveでの具体例を言うと、単にplot関数を使用して描画をすると大量のフレーム(10,000コマなど)のアニメーションにとんでもない時間がかかってしまいます。Matlabでの情報になりますが、plotではなくlineなどの関数を使用すると良いとのことです。Octaveでもこの方法は有効でした。
第四関門：計算の高速化
二重振り子ならばそれほど高速化を意識してコードを書かなくても良いかもしれません。しかし三重、四重となってくると多くのステップの計算時間は膨大なものとなってくるのである程度高速化を意識してコードを書く必要があります。
高速化の方法としては、基本的な事項になりますが、計算時間の大半が連立１次方程式(Ax=b)を解く箇所なので「inv(A)」は使用ぜずに「A\b」を使うなどがあります。
さらにホイン法では2回、A\bを解くのですが行列Aについては不変となると思われます。この場合、予め行列AをLU分解しておくと高速に解けるというのが一般的な見解と思われますが、現状のOctaveではその恩恵は受けられないようです。
(Matlabならば高速化が実現できるようです。)
以上、Octaveによる二重振り子解析に関する記事でした。
おまけ：スマホアプリにしました
余力があったのでスマホアプリにしてみました。スマホ版はOctave(インタープリタ言語)に比べて爆速で計算が行えます。なので計算しながら描画をする仕様が可能です。良かったらダウンロードして遊んでみてください。無料です。
Android: 魔改造振り子 – Double Pendulum
https://play.google.com/store/apps/details?id=jp.co.g_llc.pendulum
iPhone:魔改造振り子 – Double Pendulum i
https://apps.apple.com/jp/app/double-pendulum-i/id1563599973
