オイラー法とシンプレクティック法 その２

オイラー法とシンプレクティック法　その２

 

 

３．しかし逆変換はある

 

　しかし(6)の逆変換はあります。それが(7)です。(7)を繰り返し用いれば、(p(1)，q(1))＝(1.12253，1.47121)から初期条件(p(0)，q(0))＝(0，1)に戻れるのは明らかです。時間反転とは、時間に関して解の進行を逆向きに進める事なので、これは時間反転ではないのでしょうか？。この定義からは、時間反転ではあるでしょう。

　ただこれは、オイラー法による時間反転ではないのです。それは(6)と(7)でτを－τにかえただけの同じ計算手続きになっていないからです。ここから、オイラー法は時間反転を満たさない、と言われます。

 

４．シンプレクティック法

 

　オイラー法の時間更新手続き(6)を、以下のようにほんのわずか変更します。

 

　　（１１）

 

　明らかに(11)の局所精度は、(6)とほとんど変わらないはずです。しかし(11)は等速運動させてから等加速度運動させる計算手順になっています。このような運動を許容する運動方程式はあるはずです。じっさいシンプレクティック変換の標準的構成法によって、(11)を与える事はたやすいです。たやすいのでその詳細は省略します(^^;)。そういう訳で(11)は、シンプレクティック変換です(^^)。

 

　シンプレクティック変換ならば、何らかの運動方程式に従った運動を表します。その運動方程式は、(11)の局所精度が(6)と同等程度という事を考えれば、もとの系に対する近似運動方程式になるはずです。しかも近似エネルギー曲線も持つはずです。もとの系のエネルギー曲線の表式と近似エネルギーの表式を比較すれば、大域的精度までいっきょに明らかになるかも知れません。こうしてなんとなく性質の良さそうな近似解法への展望が開けます(^^)。

 

　要するに問題の理論的背景は、(11)の計算手順が厳密解を与える運動方程式をみつける問題に帰着されます。シンプレクティック法は、少なくとも(2)のような自励系に対しては既に確立された方法です。オーソライズされた方法は非常にスタイリッシュで、少なくとも初見では、とても読みたくなるようなものではありません(^^;)。こんな具合です・・・。

 

　まず(11)の1段目を、運動エネルギー1/2×p²に適当な偏微分作用素Ｄ_p()が作用したものとみなし、同様に2段目をポテンシャルエネルギー－1/2×q²にＤ_q()が作用したものとみなします。Ｄ_p()とＤ_q()は非可換作用素です。(11)を繰り返す計算手続き全体は、Ｄ_p()とＤ_q()から構成される指数演算子を用いて、

 

　　　　

 

と表されます。右辺のＨ(p，q)が求める近似運動方程式のエネルギーです。エネルギーがわかればその形から、運動方程式も導けます。普通の指数法則から、

 

　　

 

と思いたくなりますが、Ｄ_p()とＤ_q()は非可換なのでそうはならず、非可換作用素代数（リー代数）のBCH公式（ベイカー・キャンベル・ハウスドルフ公式）を使ってＨ(p，q)を求めます。BCH公式は、クソ力任せな計算で証明されますが、Ｈ(p，q)を求める計算もクソ力任せです(^^;)。

 

　オーソライズされた方法は、(11)が厳密解を与える系を力任せに構成するという方法です。もう一つのやり方としては、(11)が厳密解を与える系を探すという方向もあり得ます。こっちの方が多少初等的（？）と思えるのですが、非常に地道な大量の計算をしなければならないので、やっぱり詳細は省略します(^^;)。こんな具合です・・・。

　まず(3)のポテンシャルエネルギー：－1/2×q²を、一般にＵ(q)で表します。要はもとの系として、

 

　　

 

の形のエネルギーを持つ自励系を想定する、という意味です。そうするとシンプレクティック変換(11)は、

 

 

　　

 

と書けます。シンプレクティック変換(13)が厳密解を与える系を探すために、近似系のエネルギーに、

　　

 

の形を仮定し、これに変換に対して運動方程式が不変になるという条件(S-2)を要求します。要するに(13)で(14)を(Ｐ，Ｑ)に書き換えた時に、同じ形になるという要求です。

　で、・・・初等的なムチャ長ぁ～い計算を実行すると・・・、(14)の未知関数f₀，f₁，f₂，・・・に対して、次の偏微分方程式の系列が得られます。Ｄ₁，Ｄ₂を偏微分作用素、

 

　　

 

として、

 

　　

 

　ここにn＝0，1，2，・・・で、n＝0のとき右辺は0とします。また0≦k，mに対し、

 

　　

 

は、Ｄ₁，Ｄ₂に含まれる偏微分演算子の係数を、偏微分演算子に対して定数と考えてとった積です。(15)もなかなか力任せの計算になりそうですが、Ｄ₁作用素が線形演算子なので、なんとか処理できます。

　n＝0の時は、

 

　　

 

なので、f₀として、

 

　　

 

の形の解を選べるのは、けっこうすぐわかります。すなわちf₀として想定した系のエネルギーの表式をとれます。n＝1では、

 

　　

 

ですが、右辺に(17)を代入し少し頑張ると、

 

　　

 

を得ます。n＝2で、

 

　　

 

　(17)と(18)を代入し、かなり頑張ると、

 

 

　　

 

が出ます。以下同様に(15)を解き続ける事になりますが、3≦nでは右辺の項数も増えて、非常に非常に頑張らないといけない事になり、一般項を与えるためには、やはりクソ力任せ計算となります(^^;)。しかし(2)のような線形系においては、ここまでで十分です。

　(17)，(18)，(19)でＵ(q)＝－1/2×q²とします。そうすると、

　　

 

　(20)と(22)を比べると、(22)は(20)の定数倍になっているのがわかりますが、(16)，(17)から明らかなようにこれは、斉次方程式Ｄ₁ｇ＝0の特解で、1≦nで非斉次になる(15)の解の不定性を表すものだと解釈できます。従って(22)を採用するか否かは、人間が判断して良いものです。目的は、(12)に出来るだけ近い(14)を定める事でした。よって(22)の不定性を無視し、f₂＝0を採用します。

 

　3≦nで(15)右辺のf₀，f₁に関する和は、pまたはqの3階微分以上の項の和になるので、それらは0です。n＝3で右辺に残るのはf₂に関する和ですが、f₂＝0と決めました。従ってＤ₁f₃＝0となり、f₃＝0を解に選べます。以下帰納的に3≦nでＤ₁f_n＝0が得られ、f₂＝f₃＝f₄＝・・・＝0を選べます。

　以上まとめれば、

 

　　

 

が求める近似系のエネルギーです。近似運動方程式は正準方程式の手順によって、

 

 

　　

 

と得られますが、(24)の厳密解が(11)で与えられるとわかってる以上、(24)を使うことはほとんどありません。それよりも、重要なのは(12)と(23)の比較です。(23)と(12)の差は、|pq|のτの1次オーダーです。最も荒い見当で|p²|～|q²|～|pq|程度とすれば、(11)では相対誤差をτの1次オーダーに保ったまま、どこまでも計算できる事になります。こうして、数値解の大域的性質がいっきょに手に入ります。

 

　相対誤差なので今回のような単調増加する解では、絶対誤差は増え続け、局所精度ではルンゲクッタ法やオイラー法にさえ及ばない事態も起きますが、厳密解が大きくなりすぎ他の解法では計算不能になっても、シンプレクティック法は相対誤差を保ちながら計算可能な場合も良くあります。

 

　現実には無限大に逝っちゃう解にはあまり興味がありません（少なくとも工学では）。重要なのはどこまでいっても有界な解です。このときポアンカレの再帰定理より、解は1度通過した点の近傍に必ず帰ってきます。という事は相対誤差一定なので、いったん増えた絶対誤差が減少もします。このような性質は、他の解法にはないものです。ここから特に周期解（振動系）において、抜群の性能を発揮します。こういう事もあり、シンプレクティック法に数値誤差の蓄積はないといわれます。

 

 

５．シンプレクティック法の宣伝

 

　ちょっとシンプレクティック法の宣伝をしちゃいますね(^^)。1次の陽的シンプレクティック法の一般形は以下です。系のエネルギーを、

 

　　

 

として、

 

　　

 

　ちなみにネコ先生のあげたリープ・フロッグ法は、2次の陽的シンプレクティック法の一種です。系のエネルギーが(25)のように、運動量だけの関数と変位だけの関数に分離されるとき、8次の陽的シンプレクティック法まで知られています（本当はどこまでもいける）。大域的精度はn次解法ならばτのn次精度です。

 

　(26)は本当は陰解法です。時間更新手続きの両辺に(t＋τ)が出てくるからです。にもかかわらず陽的と言われるのは、あたかもオイラー法のごとく順番に代入計算してけばOKだからです。これはＶ(p)やＵ(q)が非線形であっても同じで、非線形陰解法に憑いてまわる、時間更新ステップごとの繰り返し計算はいっさい不要になります。つまり非線形現象に対しても計算時間は、陽解法とほとんど同じです（高速）。かつプログラム簡単（Excelでもできちゃう！）。

 

　しかも無条件安定の精度保証付きです！。陽解法の宿命である、解の発散や縮退は起こりません！。使ってみたくなりませんか？(^^)。

 

 

６．シンプレクティック法における時間反転

 

　1次の陽的シンプレクティック法(11)を、初期条件(p(0)，q(0))＝(0，1)，τ＝0.1で計算したものが図-3です。シンプレクティック法の方がオイラー法より若干はずれた解を与えますが、近似エネルギー曲線にはしっかり載ってます。

(p(1)，q(1))＝(1.17309，1.48394)です。さらにτ＝－0.1で時間反転を行うと・・・。

 

　あれ～(^^;)もとに戻らないじゃん。しかもオイラー法の時より、はずれが大きい。・・・そりゃそうなんですよ。

 

　τを－τにしたという事は、近似エネルギー曲線、

 

　　

 

を、

 

　　

 

に変えたって事なんですから。Ｃ＝－0.32592は初期条件(p(1)，q(1))＝(1.17309，1.48394)から計算しました。オイラー法は異なる初期条件の解曲線に数値解の点列が飛んでく解法でしたが、シンプレクティック法では逆に、人間の方で異なる初期条件の解曲線へ数値解を乗せ換えてやったようなものです。実際にやってみるとわかりますが、アクア色の点列は、後の方のエネルギー曲線にしっかり載ってるのです(^^;)。

 

　　（２７）

 

 

 

　つまりシンプレクティック法でも解を逆進させるためには、逆変換を使うしかないという事です。逆変換は(27)になります。これは(11)でq(t)，p(t)が左辺に来るように、移項しただけのものです。(27)を見ると、なんだか(11)と似ています。違いは、pの更新を先にやってからqの更新を行う点です。qとpの立場が入れ替わっています。じつはシンプレクティック変換では、変位qと運動量pの区別はないのです。qとpの事を人間が勝手にそう名付けただけ、という事になります。qとpを入れ替えるシンプレクティック変換が存在するからです。だから(27)もシンプレクティック変換なのです(^^)。従って（同じではないけれど）シンプレクティック法による、正確な解の逆進が可能です。この事態をさして、シンプレクティック法には時間反転性があるといわれます。

 

　数値解法で時間反転を考えると、このようにちょっとわかりにくい事態になりますね(^^;)。

 

 

７．q'＝qの場合

 

 

　最後に、懸案の(p(0)，q(0))＝(1，1)のケースです（q'＝qに相当）。結果は図-4になります。

 

　q'＝qの場合はエネルギー曲線が直線で、オイラー法の接線飛び出しが直線の傾きに一致して点列は正確に厳密解q＝e^tのエネルギー曲線に載りますが、シンプレクテッィク変換ではないため、誤差は蓄積して行きます。シンプレクテッィク法では最初からエネルギー曲線がずれてますが、t＝1ではオイラー法より厳密解に近いです。この場合の評価は微妙ですなぁ～(^^;)。