準シンプレクティック法による２階線形微分方程式の解法のスプレッドシート

今日、ブログで紹介した「準シンプレクティック法による２階線形常微分方程式の解法」で使用されてたスプレッドシートを、
ddt³さんからいただいたので、公開したにゃ。

エクセル版
https://docs.google.com/spreadsheets/d/e/2PACX-1vRqvxRTz9uW0sDkyxhJb4rzDr6zzJqUlatgu8jDGBNqXeKWU3qBwfk7j9l-aT6j2BSL19VMX8FmZJkS/pub?output=xlsx

シートが９つもある、巨大なファイル（2MB）なので注意するにゃ。
スプレッドシート内で使われている式そのものは単純なものなのだけれど、何しろ計算量が多いのでこうなってしまうにゃ。
近似計算による（打ち切り誤差や丸め）誤差の蓄積、累積を調べるためにには、計算する時間の範囲を長くしないといけないから、これはしょうがないにゃ。
たとえば、惑星の軌道計算は、数年とか、場合によっては数億年先まで計算しないといけない場合がある。
常微分方程式の数値計算でよく使われる（４次の）ルンゲ・クッタ法は、局所的な計算精度はいい（微分方程式の解が不連続でない場合、その計算結果は厳密解とほとんど一致する）けれど、その計算結果は力学的なエネルギー保存則を満たさないので、エネルギー保存則からの僅かなズレ、誤差が計算を進めるたびに新たな誤差を産み、その誤差が新たな誤差を生むという悪循環に陥り、雪だるま式に誤差が増大し、最終的にとんでもない計算結果になってしまうことがある。
こうした誤差の時間進展を調べるためには、計算を何度も何度も繰り返さないといけない――特に、４次のルンゲ・クッタ法は打ち切り誤差が非常に小さいから、短い計算では、誤差の蓄積（テイラー展開の打ち切りによる伝播誤差というべきか）がまったないように見える――ので、どうしてもこうなってしまう。

　　【参考記事】
　　誤差の蓄積
　　https://nekodamashi-math.blog.so-net.ne.jp/2017-11-15-1
　　オイラー法の誤差の内訳
　　https://nekodamashi-math.blog.so-net.ne.jp/2017-11-16-1

ルンゲ・クッタとは違って、シンプレクティック法はエネルギー保存則を満足するので、何度も何度も計算を繰り返しても誤差の累積は小さい。だから、局所的な誤差ではルンゲ・クッタ法に軍配が上がるけれど、大域的な誤差ではシンプレクティック法に軍配が挙がる。
こんなことを書いてもわかりづらいと思うので、例えば、ネムネコが昔に書いた次の記事などを見て欲しい。
　オイラー法とシンプレクティック法を用いた１次元調和振動子の解法
　https://nekodamashi-math.blog.so-net.ne.jp/2017-06-29-4
１次精度のオイラー法を用いた計算結果

シンプレクティック法による計算結果

局所的な（打ち切り）誤差は同レベルなのに、オイラー法は計算を繰り返す（時間を進展させる）と力学的なエネルギーがドンドン増えていくのに対し、ある範囲内ではあるが、シンプレクティック法ではエネルギー保存則が保存され、どれだけ計算を繰り返しても（長いタイムスケールで計算しても）逸脱は常にその範囲内に収まっていることがわかる。

で、（微分方程式などの）数値計算で、「精度がいい」といったら、普通、「局所的な（打ち切り）誤差が小さい」という意味なんだケロ。だけど、この他に、「大域的な誤差（数値計算の専門家、特に数学系の専門家のいう”大域的な誤差”と意味が違う異なる可能性大(^^ゞ）」というものがあって、「局所的な誤差が小さい」近似計算法が必ずしも「大域的な誤差が小さい」とは限らないんだ。
普通の数値計算をする場合、「大域的な誤差」はそれほど気にする必要はないのだけれど、惑星や宇宙探査機の軌道計算などのように数年といった長い時間スケールの計算をする場合、「大域的な誤差」を気にする必要がある。
木星探査のために打ち上げたのに、最初のうちは予定していた軌道通りに進んでいたが、次第にその予想軌道から外れ、探査機は木星ではなく土星に向かっていたなんてことが起きるかもしれないから。