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

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

 

 

　以前このプログで、時間に関する1階偏微分方程式を時間に関する2階偏微分方程式に直して、強引にシンプレクティック法を適用してさんざんな目に合ったおぼえがあるので、とりあえず今回はどうなんだろぉ～とは思いました。また数値解法上で時間反転を考えると、ちょっとわかりにくい事態になるんですよね(^^;)。

 

　ちなみにシンプレクティック法は基本的に、2階微分方程式に関する専用解法です。ここで2階微分方程式を一般に運動方程式と呼ぶとすると、運動方程式に従う運動はシンプレクティック変換を満たさねばならない事がわかります。具体的には位置と運動量の組を(q(t)，p(t))として、時間がτ進んだ状態(q(t＋τ)，p(t＋τ))への変換、

 

　　(q(t)，p(t)) → (q(t＋τ)，p(t＋τ))

 

はシンプレクティック変換になってるという訳です。

　数値解法は一般に時間ステップ幅τによって、点列、

 

　　(q(t)，p(t)) → (q(t＋τ)，p(t＋τ)) → (q(t＋2τ)，p(t＋2τ)) → ・・・

 

を与えます。この(q(t＋nτ)，p(t＋nτ))から(q(t＋(n＋1)τ)，p(t＋(n＋1)τ))を数値的に計算する過程がシンプレクティック変換になっているものを、シンプレクティック法と言います。

 

　シンプレクティック法がリープ・フロッグ法（蛙飛び法）とも言われるのは、シンプレクティック変換の標準的構成法を使うと、変位場と速度（運動量）場を時間に従って交互に更新する計算になるからです。電磁場方程式の場合は、電場と磁場を交互に更新する蛙飛びです(^^)。

 

 

１．やってみますか

 

　ネコ先生から出された宿題は、

 

　　

 

でした。qを変位と呼びます。シンプレクティック法は2階微分方程式の専用解法なので、(1)を強引に2階微分方程式に持ち込みます。(1)を時間tで微分し、

 

　　

 

すなわち、

 

　　

 

になります。

　(2)はエネルギー曲線を持ちます。両辺にq'をかけ、合成関数の微分公式の逆を使い、時間tで積分すれば、

 

　　

 

が得られます。Ｃは初期条件で決まる積分定数。

　ここでpは、

 

　　

 

で定義される運動量（速度）です。エネルギー曲線(3)とオイラー法との関係を調べるために、(2)を連立1階微分方程式へ引き直します。

 

　　

 

　(5)の1段目は(4)そのものです。2段目は1段目の意味を考慮すれば(2)そのもので、じっさいに1段目を微分して2段目に代入すれば、(2)が得られます。(5)に対するオイラー法は時間ステップ幅をτとして、

 

　　

 

でしょう。これはネコ先生が与えたオイラー法と本質的に同等な計算手順です。

 

　ネコ先生の宿題は1階微分方程式(1)に対するものだったので、初期条件はt＝0でq(0)＝1だけでした。(5)は2階微分方程式(2)と同等なので初期条件は2個必要ですが、宿題に合わせれば(1)より、t＝0でq(0)＝p(0)＝1です。しかしこの解ではオイラー法の特徴が見えにくいので、初期条件をq(0)＝1，p(0)＝0とします。

　(6)でτ＝0.1とした点列を(p，q)空間にプロットします。(p，q)空間は位相空間と言われ、微分方程式の解の定性的性質を見定める便利手段として、けっこう多用されます。図-1の計算はExcelで行いました。残念ながら、タイガー計算機（← 炊飯器じゃないよ(^^)）などは用いておりません(^^;)。Ｃは－0.5となります。

 

　時間範囲は0≦t≦1だったので、τ＝0.1から(6)を10回繰り返します（赤点列）。実際に計算すると、(p(1)，q(1))＝(1.12253，1.47121)になります。さらにこれを初期条件とし、τ＝－0.1として(6)を10回繰り返した時間反転結果が緑の点列です。最終結果は、(p(0)，q(0))＝(0.00000，0.90438)になります。という訳で、時間反転しても最初の初期条件には戻りません。

 

 

２．オイラー法の定性的（大域的）性質を調べる

 

　どうしてこうなったかを考えるために、定量的な局所精度を調べるのではなく、運動方程式の不変性に注目して定性的（大域的）にオイラー法の性質を調べてみます。そこで基準になるのは、厳密解を与える解法の条件とはどんなものかという都合の良い考えです。あるとすれば理想的な解法の条件です。

 

　　(S-1) 変換：(q(t)，p(t)) → (q(t＋τ)，p(t＋τ))は、シンプレクティック変換である。

　　(S-2) 変換：(q(t)，p(t)) → (q(t＋τ)，p(t＋τ))は、運動方程式を不変に保つ。

 

　(S-1)は冒頭で述べた、運動方程式に従う運動の必要条件です。(S-2)は、数値解の点列が同じ運動方程式を満たし続けるという条件です。(q(t＋τ)，p(t＋τ))＝(Ｑ，Ｐ)と書くことにすれば、厳密解を与える限り、(q，p)は厳密な運動方程式、

 

　　

 

を満たすはずです。同じ理由から、変換後も(Ｑ，Ｐ)に対して厳密な運動方程式が成り立たねばなりません。

 

 

　条件(S-2)さえあれば数値解法は、厳密解を与えそうな気がします。しかし同じ運動方程式を満たし続けたとしても、同じ初期条件から出発した点列であるという保証はありません。そこを制約するのが条件(S-1)です。じつはオイラー法は、条件(S-2)を満たすのです。

 

　(6)は線形変換、

 

　　

 

ですから逆行列をとる事により、

 

　　

 

が得られ、(7)を時間tで微分すれば、

 

　　

 

になるのは明らかです。

　(7)，(8)を運動方程式(5)、すなわち、

 

　　

 

に代入すれば、

　　

　ここで、

　　

ゆえに、

 

　　

 

と(q(t＋τ)，p(t＋τ))＝(Ｑ，Ｐ)に対して(9)と同じ形(10)が得られます。

 

　よってオイラー法(6)がシンプレクティック変換になっていれば、オイラー法は厳密解を与える事になります！。

 

　安心してください。ちゃんとシンプレクティック変換ではありませんよ(^^)。(S-2)を満たすオイラー法が、シンプレクティック変換なら(S-1)も満たすので、それは厳密解を与えるはずですがオイラー法の点列は図-1から、厳密解のエネルギー曲線から外れた点列になっています。よってオイラー法は自明にシンプレクティック変換ではありません。

 

　つまりオイラー法は、運動方程式は満たし続けるけれど、異なる初期条件から出発した解にジャンプし続けているのです。実際t＝0，0.1，0.2，・・・に対するオイラー法の点列(q(t)，p(t))について、エネルギーの表式である(3)を計算してみると図-2となり、オイラー法においてはエネルギーが保存しません。そしてその理由は明白です。何故ならオイラー法は計算手順(6)より、エネルギー曲線の接線方向へ飛び出し続けるというくせを持つからです。そしてその理由も明白です。

 

　(6)の時間更新手続きはqもpも同時に更新するので、等速運動させながら同時に等加速度運動もさせるというものです。そのような運動は原理的に不可能であり、当然それを許容するような運動方程式もあり得ないので、オイラー法がエネルギー保存則を満たす訳はないのです。

 

　だったら図-1の時間反転部分も同じじゃないですか！。オイラー法は定性的（大域的）に考えると、数値解がどこにジャンプして飛んで行ってしまうかわからない解法なのです。時間反転しても、最初の初期条件にもどる保証はありません。