メモ(あるいは丸コピー)です。当方はMolecular Biology of the Cellのシステム生物学の記載をきっかけに微分方程式に興味を持つようになっている。
https://twitter.com/math153arclight/status/1059706184876228610より
https://twitter.com/math153arclight/status/1059706688532410369より
↑今後自分で書いてみたいときにリンクが切れていると辛いので画像はコピーさせていただいた。 ↓以下は別の方の解説
ルンゲクッタ法は微分方程式の数値的解法としてよく用いられる方法ですが，近似が悪くなってしまう場合があります．そのような例を一つ見てみましょう．
— ふぉく (@D_A_works) 2018年11月5日
こちらは，(0, 0)を中心とした-1/rのポテンシャル中で，初期位置(-1, -1), 初期速度(0.2, 0.0)の物体の運動を，4次のルンゲクッタ法と2次のシンプレクティック法で計算したものです．時間刻み幅はdt=1.0e-3としてあります．
— ふぉく (@D_A_works) 2018年11月5日
ご覧の通り，楕円軌道からずれてしまっています． pic.twitter.com/IQW3uCz7HC
ではどこで計算がずれてしまっているのか調べてみましょう．
— ふぉく (@D_A_works) 2018年11月5日
今，エネルギーが保存するので，エネルギーの時間変化を見るのが計算精度の確認方法の王道ですね．
横軸時間，縦軸エネルギーとした図がこちらになります．
明らかに，周期的に計算精度が悪くなってるのが分かります． pic.twitter.com/6RqEwtra3E
今，周期的に変化するものとして真っ先に思い浮かぶのは距離です．
— ふぉく (@D_A_works) 2018年11月5日
そこで，今度は縦軸は同じくエネルギーにしたまま，横軸を原点からの距離にしてグラフを書いてみます．
見て分かる通り，原点付近でエネルギーが全く安定していません．
これが，誤差の原因ですね． pic.twitter.com/6L8UbZ4SpC
さて，原点付近での誤差が大きいことがわかりました．
— ふぉく (@D_A_works) 2018年11月5日
これは考えてみれば当然で，原点付近では大きな力が働くため，細かく計算しないと誤差が大きくなり過ぎてしまいます．
簡単な解決策としては，（当たり前ですが）時間刻み幅を十分小さく取ればいいわけですね．
ところが，今はそんなことはしたくありません．
— ふぉく (@D_A_works) 2018年11月5日
なぜなら，先ほどの図を見ると分かるように，時間刻み幅を小さくしたいのは原点付近だけで，それ以外の大部分では今のままでもそんなに悪くないわけです．
それなのに無闇に時間刻み幅を小さくしてしまうと，無駄に計算時間がかかってしまいます．
時間刻みを欲しい精度に応じて分けることができれば，この問題は解決です．
— ふぉく (@D_A_works) 2018年11月5日
そこで今回ご紹介するのは，Runge–Kutta–Fehlberg法です．
この方法は，計算次数の違うもの同士を組み合わせて，計算精度が常に同じくらいになるように時間刻み幅を変えながら計算する方法です．
この方法は，基本的には，時間tにおける位置xと速度vから4次と5次のルンゲクッタ法を用いて，時間t+dtにおける位置の計算値x_4, x_5を得ます．この時，x_4とx_5の差が大きい場合はdtを小さくして計算をやり直す，とすることで，計算精度を安定させる方法です．https://t.co/duTWszbDxG
— ふぉく (@D_A_works) 2018年11月5日
さて，この方法を用いて先ほどの系を計算し直しましょう．
— ふぉく (@D_A_works) 2018年11月5日
今回はこちらのサイトを参考にさせていただきました．https://t.co/fa5ppKdFwq
精度の差の上限を1.0e-2とすると，このような図が得られます．
きちんと楕円っぽくなっていますね． pic.twitter.com/q8jFb3serc
また，エネルギーの時間変化について，先ほどのルンゲクッタ法と比較して見ても，エネルギーの誤差が少なく，安定した方法であることが見てわかります． pic.twitter.com/CWFnH3K5Xn
— ふぉく (@D_A_works) 2018年11月5日
ところが，これだけではRunge–Kutta–Fehlberg法が優れているということは言えません．なぜなら，時間刻み幅を調整する機構がついているため，先ほどのルンゲクッタ法に比べて計算時間が長くなってしまっているからです．
— ふぉく (@D_A_works) 2018年11月5日
手元のパソコンではRunge–Kutta–Fehlberg法の計算に，先ほどのルンゲクッタ法の5倍程度の時間がかかってしまっていたため，先ほどのルンゲクッタ法の時間刻み幅を0.2倍して計算した結果がこの図です．
— ふぉく (@D_A_works) 2018年11月5日
Runge–Kutta–Fehlberg法と同じように楕円軌道を描いているように見えますね． pic.twitter.com/7jfPtxVMgG
ところが，エネルギーの時間変化を見て見ると，時間刻み幅を短くしたルンゲクッタ法よりも今回のRunge–Kutta–Fehlberg法の方が，精度が良いことがわかります．
— ふぉく (@D_A_works) 2018年11月5日
これはおそらく，原点付近での時間の刻み幅が更に小さいために，エネルギーのずれば少なくなっているのだと考えられます． pic.twitter.com/rLLGxUAP1c
以上より，Runge–Kutta–Fehlberg法は，いわゆる普通の4次のルンゲクッタ法に比べて，精度が高い方法であることが確認されました．
— ふぉく (@D_A_works) 2018年11月5日
実際には，今回のように力の大きさや向きが急に変わるような系でなければ必要はないでしょう．
しかし，アルゴリズムも難しくないので，一度組んでみては如何でしょうか.
ところでツイッターも便利であるが、友人と連絡を密に取らずとも友人は友人であるはずなのでもうやめても良いように感じている。
