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

今日、ブログで紹介した「準シンプレクティック法による２階線形常微分方程式の解法」で使用されてたスプレッドシートを、
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
１次精度のオイラー法を用いた計算結果

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

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

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

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

［仲間に入れて欲しいなぁ～(^^)］

　いまさら次郎なんですが、

 

　　https://nekodamashi-math.blog.so-net.ne.jp/archive/c2306081239-1

 

の2018-09-16に対するコメントです。

　　

　2階微分方程式じゃないっすか！。2階微分方程式なら準シンプレクティック法でしょうと、個人的に思うだけですが(^^;)、だから仲間に入れて下さいな(^^)。

 

　シンプレクティック法は2階微分方程式専用の方法で、2階微分方程式を運動方程式と解釈します。特に(1)でy'の項がない場合を保存系といい、保存系に対するシンプレクティック法は確立されています。

　シンプレクティック法の最大の特徴は、近似計算に対応する近似運動方程式が存在する点で、保存系の場合、そのために厳密なエネルギー保存則を大域的に一様な精度で再現する、近似エネルギー保存則がある事になり、理論的には誤差が蓄積しないという結論になります。

　シンプレクティック・オイラー法の計算手順は極めて簡単です。等速運動させてから等加速度運動させるです。そうすれば近似運動方程式があるはずだからです。(1)でy'＝pとおけば(1)は、

　　

と連立1階微分方程式に書き直せます。(2)から導かれるシンプレクティック・オイラー法の計算手続きは、

　　

となります。Δtは時間分割幅です。

　(3)の上段が等速運動です。下段は、その結果を用いた等加速度運動となります。ただし(3)のような定式化は公式には認められていないので、ここでは(3)を準シンプレクティック法と呼びます。もし(2)で－3pの項がなければ、公に認められたシンプレクティック法（保存系）です。

 

　自分はシンプレクティック法の本質は、近似計算に対応する近似運動方程式が存在する事だと思っています。(1)はエネルギーが減少する減衰系で保存系ではありませんが、近似計算に対応する近似運動方程式が存在するというシンプレクティック法の本質は、保存系であろうとそうでなかろうと同じと思う訳です。証明出来ませんけど・・・(^^;)。

　なのでエネルギー減衰系にシンプレクティック法をそのまま適用しても、減衰するエネルギー曲線を大域的に一様な精度で再現する解が得られるだろうと考えます。もう一回言いますけど、証明できていません(^^;)。

 

　(3)の下段をp(t＋Δt)について解きます。結果は(4)です。

　　

　(4)を単純にExcelで計算した結果が図-1です。Δt＝0.1。

 

　青丸のS-1が準シンプレクティック法，黒丸のEuがオイラー法。青点線と黒点線が、対応する相対誤差。

　・・・あれ？、準シンプレクティック法がオイラー法に負けてるじゃないの。シンプレクティック法はオイラー法より、はるかにはるかに精度が良いはずじゃないのか。どうした、準シンプレクティック！(^^;

 

　ところで(3)を見た時、右辺のpはp(t＋Δt)じゃなくp(t)じゃ駄目なのかい？って気はしませんか。

　－3p(t＋Δt)Δtは、時間間隔Δtの間のpに比例した運動量減衰の近似です。シンプレクティック法の特徴は、運動量pの更新に、更新された変位y(t＋Δt)を使うところなので、減衰項3pについてもなんとなくp(t＋Δt)を使った訳ですが、3p(t)じゃ駄目なんですか？（← と自分で言うな(^^;)）。そこで単純に、

　　

を採用してみると図-2となり、突如精度が改善されます。図-2のS-2です。

　しかしいずれにしろ、3p(t)Δtであろうと3p(t＋Δt)Δtであろうと、この項は時間間隔Δtの間のpに比例した運動量減衰の近似には違いありません。だったらt～t＋Δtに台形公式を適用して、

　　

とする手だってあるはずです。

　(6)から(4)に相当するものを作ると、

　　

 

です。結果は図-3のS-3。

 

 

 

　予想通り、S-3はS-1とS-2の中間に来ました。精度はS-1についでいいくらい。

　・・・でもシンプレクティック法が真価を発揮するのは、長時間積分において。短時間積分では、ルンゲクッタ方なんかには絶対に後れを取ります。今回はオイラー法にさえ後れを取った。という訳で、時間ステップ幅をΔt＝0.01と1/10にし、t＝0～30としてステップ数を100倍にしてみよう

 

 

 

　ステップ数を100倍にした時のt＝0～3の挙動は図-4です。Δtを1/10にすると、どの解法もそれなりに厳密解を近似するようになり、精度は10倍になって、精度の良さの順番はさっきと同じ。がしかし、t＝0～30での相対誤差を比較してみると、S-3だけが安定している。他の解法は、誤差の蓄積が如実に見える！。

　いちおうのシンプレクティック法と思えるS-1は、局所的精度でもオイラー法に負けてるので、長時間積分もオイラー法より成績が悪い。S-2は半オイラー法なので局所的精度では上回ったが、後半はハイブリッドの悪い処取りが災いして、一番成績が悪くなった気がする。S-3が最もシンプレクティック法らしい挙動。何故こうなのか？。

　(3)，(5)，(6)の2段目の減衰項を移項すると、

　　

となる。(3')，(5')，(6')の右辺はどれも、減衰がない場合の本来のシンプレクティック法の運動量更新になっている。一方左辺は、減衰による運動量減少を逆向きに加えたもの。

　ここでダランベールの原理のように考える。真の運動量p(t＋Δt)に減衰による運動量変化を逆向きに加えたら、無減衰運動（保存系）と同じだと。そういう仮想系には、本来のシンプレクティック法を適用出来るはず。

　よって減衰による運動量変化を最も正確に評価した(6')すなわちS-3が、最もシンプレクティックぽくなったと考えられる。つまり(3')すなわちS-1では、中途半端だった訳だ。

 

　・・・という解釈はいちおう成り立つんだけれど、準シンプレクティック法にはまだ色々考えなきゃならん事があるようで・・・。しかたないので、時間推進演算子を本気で勉強しようかな(^^;)。