私はシミュレーションではシンプレティック数値積分法が大好きでよく使います。
何がいいって、実装が簡単で、長く計算しても発散しない（誤差がある一定を超えないこと）。
アメリカPCで見ようと思ったら、文字化けしてみえねぇー

申し訳ないと思いつつ、自分のブログに載せて見させてもらいます。
http://www-cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/

ちなみに、この特性をKalman Filterに応用できれば長期に使っても誤差が発散しないセンシングができるのでは？って考えて、調べ中です。

第1章 Symplectic Integratorの紹介

Euler法やRunge-Kutta法などの差分法は、ある物理量の現在の時間微分から未来を予測する方法で数値積分を行なっていました。しかし、この方法では原理的に誤差が生じ、それが蓄積され、系のエネルギーの発散などの実際にはあり得ない結果となります。
Δt を小さくすると、その誤差は少なくなりますが、逆に計算回数が増えるので、かけ算などでの計算誤差が入ってきてしまいます。
数値積分ではこの誤差は仕方のないことだと思われがちなのですが、この問題を解決する新しい数値積分法が1980年代の終りごろから登場しました。その名も Symplectic Integrator(シンプレクティック数値積分法)です。そして 1990年代前半に鈴木増雄先生(現東京大学名誉教授) によって演算子の指数の指数積展開の理論として体系化されました。
この数値積分法はいままでの差分法とは全く異なる理論により作られています。その詳細は以降の章で解説することにして、その結果は系のエネルギーを原理的に保存する数値積分法となるのです。しかも使い方がEulerによく似ていて誰でも簡単に利用できるのです。
この数値積分法の利用分野は古典力学系に限らず量子多体系、乱流解析にも利用可能と言われています。まさに次世代のコンピュータ物理学を支える重要な理論かつ技術なのです。
この数値積分法がどのくらい凄いものなのかをお見せしましょう。例として、鈴木増雄先生が論文で使われたHamiltonianを用います。そのHamiltonianは、

(00)
と与えられます。初期値としてt=0でｑx=2,ｑy=1,ｐx=0,ｐy=0を設定します。

この系のエネルギーは理論的には初期値 E=0で一定のはずです。では数値積分ではどうなるでしょう。この系のエネルギーの時間発展を1次Euler法、 4次Runge-Kutta法、1次Symplectic法、4次Symplectic法でそれぞれ dt=0.001でt=2000まで計算した結果を以下に載せます。

1次Euler法：あっというまにエネルギーは発散してしまいt=20までしか表示できませんでした。もし発散しなければt=2000までの計算にかかる時間は40秒です。

4次Runge-Kutta法：t=300付近から段階的にエネルギーが減少しています。エネルギーの発散が起こらないだけEuler法よりはましでしょう。計算にかかった時間は300秒です。時間がかかります。

1次Symplectic法：エネルギーが奇妙な振舞いをしています。いたる所々で非常に激しい振動が起こっています。まるでカオスの分岐図がフラクタルの様にならんでいる様です。振動の後にはエネルギーは元の値に戻ります。計算にかかった時間は40秒です。1次Euler法と同程度です。

4次Symplectic法：エネルギーが上向きのみに所々で振動しています。振動の振幅は1次Symplectic法のものより1000分の1程度に小さくなっています。これも振動の後にはエネルギーは元の値に戻ります。計算にかかった時間は130秒です。1次Symplectic法の約3倍です。
このようにSymplectic法はエネルギーが振動するにもかかわらずまた元のエネルギーに落ち着くという性質がある数値積分法なのです。また計算時間がそれほどかからないことも長所です。このエネルギーの振動は次数を上げることで劇的に小さくすることができるのです。
Symplectic法のすばらしさがおわかりいただけたことでしょう。ではその詳細な解説を次章から始めます。

第2章微分作用素による運動方程式の記述

1粒子の運動をHamilton形式で見ることから始めましょう。粒子の位置を q 、運動量を p とし、この粒子の Hamitonianを H とします。運動方程式は以下のようになります。

(01)
ここで (p,q) をベクトルとして見てみましょう。すると、このベクトルの時間変化は、Hamitonianの形と、そのベクトル自身によって決まることになります。
つまりベクトルに対して、Hamitonianから決まる特殊な作用を施すと、その時間変化のベクトルになると見ることができます。その作用素を D と表すことにしましょう。

(02)
このように書くと D が行列の様ですが、行列はなく、ベクトルからその時間変化を表すベクトルに変換する作用素です。この作用素は時刻tによって変化しません。 Hのみによってきまるのです。
ところで、普通のスカラー量のx,kに対しての微分方程式とその解

(03)
と(02)式の作用素の方程式とを比較するとその解は

(04)
となるような気がしませんか？

ここで、作用素の指数が登場しましたが、これの意味は次の通りです。

(05)
つまり exp のテイラー展開のそのままです。テイラー展開ですので、t があまり大きいと具合がわるいので、以後は有限の短い時間 Δt に置き換えます。 D を (p,q) に作用させると、(02)式よりその1階微分が求まります。
D を (p,q) に2回作用させると、

(06)
となるので2階微分が求まります。

よって D を n 回作用させると (p,q) の n 階微分が求まることになります。

(07)
ここで(04)式にこの指数関数の定義式を代入し、作用の処理を施すと次の結果になります。

(08)
これは (p(Δt),q(Δt)) の (p(0),q(0)) まわりのテイラー展開になるではないですか。よって(04)式には正当性があるのです。

次にこの Exp{Δt D} という謎の作用素が具体的にどのようなものなのかを調べてみましょう。
D はH から決まる作用素です。H は運動エネルギー K とポテンシャルエネルギー Uとの和でした。よって作用素D も K による要素と U による要素の和になるのです。
つまり、それぞれを、Du Dk とすれば、
D = Du + Dk
と分けられるということです。

スカラーの指数関数の法則 exp{A+B}=exp{A}exp{B} がこの作用素の指数に対しても成り立つとしたらどうなるでしょう。

(09)
これが成り立てば、(04)式は

(10)
となります。

つまり、(p(0),q(0)) を Exp(Δt Dk)で変換したものをさらに Exp(Δt Du) で変換することになるのです。
では、この Exp(Δt Dk) (p(0),q(0)) とは、具体的に何なのでしょう？
Dk は、U が無い系での H が作る作用素 D です。つまり自由粒子というわけです。自由粒子ならばその運動の計算は簡単かつ厳密です。運動量は変化しないで、位置が一定方向に変化するだけですから。
なので、Exp(Δt Dk)(p(0),q(0)) が表すベクトルは粒子を自由粒子としてΔt後の状態まで積分しているのです。自由粒子の運動なので、積分は厳密に行われます。この厳密な積分が Symplectic法の鍵となっています。
結局、作用素 Exp(Δt Dk) の効果は次式で表されます。

(11)
次に、Exp(Δt Du)(p(0),q(0)) について同様に考えます。
D_u は、K がない H の Dです。運動エネルギーがないというのは意味がわかりませんが、Uだけなので、運動方程式より位置 q が変化しないで、運動量 p だけが変化していることになります。q が変化しないとU が変わらないので p の変化は一定になります。よってこれも自由粒子と同じ様に積分が厳密にできます。その効果は次式で表されます。

(12)
このように、作用素 Exp(Δt D) が Exp(Δt Dk) と Exp(Δt Du) の積に分離できると、個々の作用素は粒子の運動を厳密に計算できるタイプの作用素となります。その作用素を組み合わせた作用素もまた粒子の運動を厳密に計算できるものになるのです。
先の2つの作用素を (p(0),q(0)) に作用させて得られる (p(Δt),q(Δt))は

(13)

(14)
となります。これを1次のSymplectic解法と呼びます。
あまりEuler法とかわりないように思われるでしょうが実は大きな違いがあるのです。

第4章非可換微分作用素の指数法則

(13)(14)式の1次のSymplectic法がEuler法とあまり変わり無いように見えるのは、実は仮定した指数法則が成り立っていないからなのです。
でも、なんらかの形で二つの作用素が積の形で分離できると、運動の厳密な予測ができるのです。実はそれが無限積の形で実現されるのです。
A,Bを一般の非可換の作用素とし、tをスカラーとすると、

(15)
計算で無限積にするのは無理なので有限積となります。この係数 u_i,k_iを求めるには、上式をExpの定義にしたがってテイラー展開し、作用素の同じ項でまとめ、その係数を比較し、 u_i,k_iの連立方程式を作り、それを解くことによって得られます。

試しに一つ例として求めてみましょう。有限積で、さぼって r=2 とします。さらにt の3乗のオーダーを無視することにします。 (15)式は

(16)
となります。この左辺は

(17)
と展開されます。右辺は少し込み入ってますが、t^3 を無視して計算をすすめます。

(18)
これから係数を比較して、u_i,k_iの連立方程式を作ります。

(19)
独立でない方程式があるのでこれだけでは係数は一意には定まりませんが、これは我々に係数を選ぶ権利があることになります。なので展開が最も簡単になるように k_2=0 としましょう。こうすると解は次のようになります。

(20)
が求まりました。

よって(15)式の2次近似式は

(21)
となります。指数の並び方が左右対称になっていますね。実はこの対称性が後で重要になります。そのため先程、 k_2=0 を選んだのです。

この結果を(04)式に用いると、

(22)
作用素が分解されたので具体的な計算にはいれます。

まず最初の作用素による変換で (p(0),q(0)) が (p_1,q_1) に写ります。

(23)
次の作用素による変換で (p_1,q_1) が (p_2,q_2) に写ります。

(24)
最後の作用素による変換で (p_2,q_2) が (p(Δt),q(Δt)) に写ります。

(25)

こうして1つの時間ステップ \Delta t を進めることができました。以上のプロセスで運動方程式を解いていく方法を 2次のSymplectic解法と呼びます。

第5章影のHamiltonian

Symplectic法の作り方と使い方はわかりましたが、この方法がなぜそれほど注目されているかがまるでわかりません。Symplectic Integratorの最大の特徴のエネルギー保存について解説します。
作用素を無限積で表示して、それを有限積で近似しました。このときの近似で数値積分に誤差が起こります。しかし、その誤差の振る舞に、とてつもない性質があるのです。
指数関数の有限積は、一つの指数関数に必ずまとめらることが知られています。すなわち、

(26)
この D' はもちろん D とは違います。その違いは、 u_i,k_i を求めるときに無視した Δt^nのオーダーです。さらにこのD' に対して、これを生成するハミルトニアンH' が存在します。
有限積による積分では、この影のハミルトニアン H' の系での粒子の運動を完全に予測するものなのです。エネルギーもこの系では厳密に保存されます。
H' の H との差は Δt^n のオーダーなので、 H'から解かれた(p,q)でHの値を求めると、振動はするものの、その振幅はΔt^nのオーダーであって決してエネルギーの単調増大減少は起こり得ないのです。

これがSymplectic数値積分法がエネルギーを保存する原理です。近似度を高めることは、エネルギーの振動の振幅を小さくするためなのです。

第6章指数作用素のフラクタル分割(前編)

より高次のSymplectic Integratorを構築するには(15)式の指数積展開をより細かく分割しなくてはなりません。例として指数4項(実質3項)からなる2次の分割を(16)-(22)式に示しましたがたった2次でも計算がこれほど面倒です。この方法で高次の分割を求めるのは非常に困難です。
1990年にこの分割を容易に無限次まで行う画期的な方法が鈴木増雄先生(当時東京大学理学部物理学科教授、現同名誉教授東京理科大学教授)によって考案されました。この分割により細分化された時間ステップの様子にフラクタルの性質があることからFractal Decomposition (フラクタル分割)と命名されました。
この章ではこのフラクタル分割の理論について解説します。

指数作用素の分割のおさらい

非可換な作用素 A,B の和の指数にはスカラーの指数分解公式が成り立ちません。そのためこの指数を個々の作用素の指数の積として表すためには(15)式の無限積展開とならざるをえません。 tのm次までの展開で分割します。

(27)
この指数の積の数 rはm次分割が得られるように調整されます。各指数の係数 a_i,b_iを決定することがこの展開作業の主題です。 Πの部分を以後 S_m(t)と表すことにします。

(28)
S_1(t)は、明らかに

(29)
です。S_2(t)は、(22)式の通り

(30)
です。

2次分割による3次分割の構成

S_3(t)を求めるには少々工夫が必要です。この工夫が鈴木先生が思い付かれた方法です。
まずExp(t(A+B))を次のようにある定数sを用いて書き改めます。

(31)
作用素(A+B)と(A+B)は可換なのですから、この等式は近似ではなく厳密な恒等式です。

(31)式の右辺の2つの指数のそれぞれをtの3次の項まで Tylar展開します。そしてこの3次の項が全体の指数積に及ぼす影響を考えると、それはこの2つの3次の項の和であることに気が付きます。つまり

(32)
の項のことです。もしこの項が0になるのなら、それは(31)式の 2つの指数のそれぞれのtの3次の項は全体の指数積になんら影響を及ぼさないことになります。つまりこれらに誤差が含まれていても全体ではその誤差はなくなるのです。全体が 3次となり4次の項で始めて誤差がでるのです。

(32)式の項が0になる条件は明らかに

(33)
です。このとき(31)式の2つの指数をそれぞれ2次のS_2(t) で表すと、全体は3次になることはもうおわかりでしょう。

結局、

(34)
となります。これで3次の指数積分割はできあがりです。

m-1次分割によるm次分割の構成

一般にm次の分割を考えます。m-1次のそれが求まっているとしてこれを用いて再帰的にm次を構成します。
(31)式の2つの指数をそれぞれm-1次のS_{m-1}(t) で表して、S_{m-1}(t)が持つm次の不正確さが全体に現れないように、sを調整します。すなわち次の通りです。

(35)
こうして再帰的に無限の精度まで高次の指数積分解が可能なのです。各指数の係数a_i,b_iはこの方法によりコンピュータで簡単に求まります。

第7章指数作用素のフラクタル分割(後編)

より安定な分解法

先の方法により無限次のSymplectic Integratorが構築されるのですが同じ次数のSymplectic Integratorによる計算でも、エネルギーの振動をより押えた安定なSymplectic Integratorを構築することができます。
(31)式でのExp(t(A+B))の恒等分解を更に複雑な形で恒等分解することにより、その安定なSymplectic Integratorが構築されます。
変形の仕方にはいろいろあるのですが、例えば(35)式を複雑化して

(36)
と再帰定義する分解方法があります。この方法の明瞭な利点は s_mが実数であるため物理計算に即実用可能なことです。
また分割の形が前後でSymmetric(対称)です。時間推進作用素であるこの作用素がSymmetricであると、時間反転に対して対称となります。つまりS_m(t)S_m(-t)=1 となります。この種の分割をSymmetric Decompsitionと呼びます。

さらに別な形式の変形として鈴木先生が推奨されている次の分割を紹介します。

(37)
これもSymmetric Decompsitionです。1次あげることで指数の項の数が約5倍になるので計算量が著しく多くなるのが困りますが、数値積分の精度が遥かにあがります。

下図に(36)式型と(37)式型の分解によるSymmetric数値積分の結果を載せます。Hamitonianなどの系の条件は (00)式で例示したものと同じです。

3項分解型3次Symplectic法：計算にかかった時間は 130秒です。

5項分解型3次Symplectic法：計算にかかった時間は 180秒です。振動の振幅が3項型のものよりさらに20分の1になっています。

2m-1奇数次精度と2m偶数次精度の関係

Symmetric Decompsitionのうち奇数次の分解S_{2m-1}(t) についてさらに考えます。元の指数とこの分解との関係は次の通りです。

(38)
ここでR_2m(A,B)は、S_{2m-1}(t)の不正確さのうちt^{2m}の項のA,Bによる多項式の係数を表します。このR_2m(A,B)の意外な性質を調べます。

(38)式のtを-tに反転した式を用意してそれと(38)式の辺辺を掛け合わせます。S_{2m-1}(t)が Symmetric Decompsitionであることを考慮すると

(39)
となります。t=0を代入してもこの式は成り立つはずなので結局 R_2m(A,B)=0となってしまいます。

(38)式を見直せば、S_{2m-1}(t)にはt^{2m}の不正確さがないことがわかります。つまり2m-1次の分割 S_{2m-1}(t)は実は2m次の分割も兼ねるのです。
これにより(37)式のタイプの分解法は次のように拡張されます。

(40)
このタイプの分解法ではS_1(t),S_2(t)は共に(30)式を用います。偶数次分割の2次おきに高次の分解を作るので、指数の総数が格段に少なくなります。

フラクタル分解

(40)式で示される各次数の指数積展開での各項の展開係数をコンピュータで計算します。
Symplectic Integratorでのこの展開係数の意味は、短い時間ステップ Δtうちの各係数の分を作用素 A(=Du),B(=Dk)のどちらか一方だけで時間発展するという意味です。2種の作用素があるので時間はのべ2Δt進めることになります。
各係数の指数ひとつを作用する作業を1stepとして、各stepで時間がどのように進んでいくかを以下の図に示します。
4次分解S_4(t)

6次分解S_6(t)

8次分解S_8(t)

次数があがるとこのように折れ線グラフがより複雑になります。その複雑さの入り組み方はまさにfractal的です。高次分解が低次分解を再帰的に組みあわせていることからもその構造がfractalになることは察しがつきます。

このような性質からこの指数作用素の分解はfractal decomposition と命名されました。

第8章 Symplectic Integratorのプログラム例

Symplectic Integratorのための分解係数の求め方が完全に理解できたところで実際にこれを計算してみましょう。前の節で紹介したFractal分解は、再帰的な定義式でしたので、プログラムでも再帰呼び出しで係数を計算するのが簡単です。
ここに簡単な分解係数の計算プログラムを用意しましたのでこれを眺めて使って下さい。
Fractal分解の係数を計算するプログラム fracdeco.cc のダウンロード
係数が求まったら、これを使って実際にSymplectic数値積分を行ってみましょう。計算する物理系は本文章中でも例として使われた (00)式の系です。
example.ccのダウンロード
少々不完全な所もありますが、C言語とSymplectic法を理解していれば対処できるでしょう。

//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
//  Program example.cc
//  Fractal分解を利用した Symplectic integrator による
//  物理計算のプログラム例
//  Copyright (C) by 渡辺尚貴 1998年 8月14日
//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=

#include <stdio.h>
#include <math.h>
#include "vector2.h"
#include "nxgraph.h"

// fracdeco.cc の出力結果をこのファイルにでも入れておいて下さい。
#include "symplectic.h"

//---- 諸定数の設定(好み合わせて変更して下さい)
#define dT             (1.0/1024)  // 基本 時間刻幅 
#define FINAL_T        (2048.0)    // 最終時刻

#define POINTS_PER_DOT (16)        // 縦1dot当たりのプロット数

// プロットする際のエネルギー軸の拡大率(かなり大きな値が必要です)
#define MAG            (128*256.0) // 2次の分割向け
//#define MAG            (128*128*256.0) // 4次の分割向け
//#define MAG            (128*128*128*128*128.0) // 6次の分割向け
//#define MAG            (128*128*128*128*256.0) // 8次の分割向け


//---- 諸定数の定義(手を触れないで下さい)
#define SCALE_WIDTH    (512)       // 横軸の長さ(dot)
#define SCALE_HEIGHT   (256)       // 縦軸の長さ(dot)
#define XMARGIN        (40)        // 画面の枠と軸との隙間の横幅(dot)
#define YMARGIN        (16)        // 画面の枠と軸との隙間の縦幅(dot)

#define WIN_WIDTH    (SCALE_WIDTH +2*XMARGIN)
#define WIN_HEIGHT   (SCALE_HEIGHT+2*YMARGIN)
#define XO           (XMARGIN)
#define YO           (YMARGIN+SCALE_HEIGHT/2)
#define STEP_PER_POINT (int(FINAL_T/SCALE_WIDTH/POINTS_PER_DOT/dT))


//---- 座標軸や目盛などを描く関数
void MakeSheet( void )
{
  double ymax, yrange, y_value;
  int y_height;

    // 縦軸と横軸を描く
  NXDrawLine( XO, YO-SCALE_HEIGHT/2, XO, YO+SCALE_HEIGHT/2 );
  NXDrawLine( XO, YO+SCALE_HEIGHT/2, XO+SCALE_WIDTH, YO+SCALE_HEIGHT/2 );

    // 拡大率(MAG)に合わせて縦軸の目盛を選んで描く
  ymax   = (SCALE_HEIGHT/2)/MAG;
  yrange = pow( 10, (int)floor( log10( ymax ) ) );

  if( 5*yrange < ymax ){
    y_value  = 5*yrange;
  }else if( 2*yrange < ymax ){
    y_value  = 2*yrange;
  }else{
    y_value  = yrange;
  }
  y_height = int(y_value*MAG);

    // 縦軸の＋側に目盛を入れる
  NXDrawMoveto( XO-3, YO-y_height );
  NXDrawLinerel( 6, 0 );
  NXDrawPrintf( 0, YO-y_height+4,"%+5.0e", +y_value );

    // 縦軸の原点に目盛を入れる
  NXDrawPrintf( 0, YO+4, "%0" );

    // 縦軸のー側に目盛を入れる
  NXDrawMoveto( XO-3, YO+y_height );
  NXDrawLinerel( 6, 0 );
  NXDrawPrintf( 0, YO+y_height+4,"%+5.0e", -y_value );

    // 縦軸の変数名を描く
  NXDrawPrintf( XO-16, YO-SCALE_HEIGHT/2, "Energy" );

    // 横軸の原点に目盛を入れる
  NXDrawPrintf( XO-6, YO+SCALE_HEIGHT/2+16, "0" );

    // 横軸の中間に目盛を入れる
  NXDrawMoveto( XO+SCALE_WIDTH/2, YO+SCALE_HEIGHT/2-3 );
  NXDrawLinerel( 0, 6 );
  NXDrawPrintf( XO+SCALE_WIDTH/2-12, YO+SCALE_HEIGHT/2+16, "1000" );

    // 横軸の右端に目盛を入れる
  NXDrawMoveto( XO+SCALE_WIDTH, YO+SCALE_HEIGHT/2-3 );
  NXDrawLinerel( 0, 6 );
  NXDrawPrintf( XO+SCALE_WIDTH-12, YO+SCALE_HEIGHT/2+16, "2000" );

    // 横軸の変数名を描く
  NXDrawPrintf( XO+SCALE_WIDTH-16, YO+SCALE_HEIGHT/2-8, "Time" );

    // Symplectic法の種類を表す数値を描く
  NXDrawPrintf( WIN_WIDTH/2, YMARGIN, "m=%d  r=%d", symp_m, symp_r );

  NXFlush();
}

//---- 時間を横軸、エネルギーを縦軸にプロットする関数
void Plot( double t, double E )
{
  NXDrawPoint( XO+int(SCALE_WIDTH*t/FINAL_T), YO - int(MAG*E) );
}

//---- 質点に働く力を計算する関数
Vector2 Force( const Vector2& q )
{
  return( Vector2(-q.x*q.y*q.y, -q.x*q.x*q.y ) );
}

//---- 質点のエネルギーを計算する関数
double Energy( const Vector2& q, const Vector2& p )
{
  return( 0.5*(p.x*p.x + p.y*p.y + q.x*q.x*q.y*q.y) - 2.0 );
}


int main( void )
{
    // 質点の位置 q と運動量 p の変数の定義と初期化
  Vector2 q(2,1), p(0,0);

  NXOpenWindow( "example of symplectic integrator", WIN_WIDTH, WIN_HEIGHT );

    // 座標軸などの描画
  MakeSheet();

    // 時間発展のループ
  for( double t=0.0 ; t<FINAL_T ; ){
      // いくらかstepを踏んでからプロットするためのループ
    for( int step=0 ; step<STEP_PER_POINT ; step++ ){
        // Symplectic integratorの心臓部
      for( int r=0 ; r<symp_r ; r++ ){
 Vector2 f =  Force( q );   // 力の計算
 p += dT * symp_u[r] * f;   // Duの効果による運動量の更新
 q += dT * symp_k[r] * p;   // Dkの効果による位置の更新
      }
      t+=dT;
    }

      // エネルギーの時間変化を描く
    Plot( t, Energy( q, p ) );
  }

  printf("プログラムは終了です。クリックしてください。\n");
  XEvent ev;
  NXCheckEvent( NX_WAIT, ev );
  NXCloseWindow();
}

アメリカ博士号＋αへの道（In Georgia Tech）

2015年5月1日金曜日

シンプレティック数値積分法

第1章 Symplectic Integratorの紹介

第2章微分作用素による運動方程式の記述

第4章非可換微分作用素の指数法則

第5章影のHamiltonian

第6章指数作用素のフラクタル分割(前編)

指数作用素の分割のおさらい

2次分割による3次分割の構成

m-1次分割によるm次分割の構成

第7章指数作用素のフラクタル分割(後編)

より安定な分解法

2m-1奇数次精度と2m偶数次精度の関係

フラクタル分解

第8章 Symplectic Integratorのプログラム例

0 件のコメント:

コメントを投稿

2015年5月1日金曜日

シンプレティック数値積分法

第1章 Symplectic Integratorの紹介

第2章 微分作用素による運動方程式の記述

第4章 非可換微分作用素の指数法則

第5章 影のHamiltonian

第6章 指数作用素のフラクタル分割(前編)

指数作用素の分割のおさらい

2次分割による3次分割の構成

m-1次分割によるm次分割の構成

第7章 指数作用素のフラクタル分割(後編)

より安定な分解法

2m-1奇数次精度と2m偶数次精度の関係

フラクタル分解

第8章 Symplectic Integratorのプログラム例

0 件のコメント:

コメントを投稿

第2章微分作用素による運動方程式の記述

第4章非可換微分作用素の指数法則

第5章影のHamiltonian

第6章指数作用素のフラクタル分割(前編)

第7章指数作用素のフラクタル分割(後編)