何がいいって、実装が簡単で、長く計算しても発散しない（誤差がある一定を超えないこと）。
アメリカ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)と与えられます。初期値として
この系のエネルギーは理論的には初期値
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形式で見ることから始めましょう。 粒子の位置を
(01)ここで
つまりベクトルに対して、Hamitonianから決まる特殊な作用を施すと、 その時間変化のベクトルになると見ることができます。 その作用素を
(02)このように書くと
ところで、普通のスカラー量のx,kに対しての微分方程式とその解
(03)と(02)式の作用素の方程式とを比較するとその解は
(04)となるような気がしませんか？
ここで、作用素の指数が登場しましたが、これの意味は次の通りです。
(05)つまり
(06)となるので2階微分が求まります。
よって
(07)ここで(04)式にこの指数関数の定義式を代入し、作用の処理を施すと 次の結果になります。
(08)これは
次にこの
つまり、それぞれを、
D = Du + Dk
と分けられるということです。
スカラーの指数関数の法則
(09)これが成り立てば、(04)式は
(10)となります。
つまり、
では、この
なので、
結局、作用素
(11)次に、
(12)このように、作用素
先の2つの作用素を
(13)
(14)となります。これを1次のSymplectic解法と呼びます。
あまりEuler法とかわりないように思われるでしょうが 実は大きな違いがあるのです。
第4章 非可換微分作用素の指数法則
(13)(14)式の1次のSymplectic法がEuler法とあまり変わり無いように 見えるのは、実は仮定した指数法則が成り立っていないからなのです。
でも、なんらかの形で二つの作用素が積の形で分離できると、運動の厳密な 予測ができるのです。実はそれが無限積の形で実現されるのです。
(15)計算で無限積にするのは無理なので有限積となります。この係数
試しに一つ例として求めてみましょう。有限積で、さぼって
(16)
(17)
(18)
(19)
(20)よって(15)式の2次近似式は
(21)この結果を(04)式に用いると、
(22)まず最初の作用素による変換で
(23)次の作用素による変換で
(24)最後の作用素による変換で
(25)こうして1つの時間ステップ
第5章 影のHamiltonian
Symplectic法の作り方と使い方はわかりましたが、この方法が なぜそれほど注目されているかがまるでわかりません。Symplectic Integratorの最大の特徴のエネルギー保存について解説します。
作用素を無限積で表示して、それを有限積で近似しました。このときの 近似で数値積分に誤差が起こります。しかし、その誤差の振る舞に、 とてつもない性質があるのです。
指数関数の有限積は、一つの指数関数に必ずまとめらることが 知られています。すなわち、
(26)この
有限積による積分では、この影のハミルトニアン
これがSymplectic数値積分法がエネルギーを保存する原理です。 近似度を高めることは、エネルギーの振動の振幅を小さくするためなのです。
第6章 指数作用素のフラクタル分割(前編)
より高次のSymplectic Integratorを構築するには(15)式 の指数積展開をより細かく分割しなくてはなりません。例として 指数4項(実質3項)からなる2次の分割を(16)-(22)式に示しましたが たった2次でも計算がこれほど面倒です。この方法で 高次の分割を求めるのは非常に困難です。
1990年にこの分割を容易に無限次まで行う画期的な方法が 鈴木増雄先生(当時 東京大学理学部物理学科教授、現 同名誉教授 東京理科 大学教授)によって考案されました。この分割により細分化された時間 ステップの様子にフラクタルの性質があることからFractal Decomposition (フラクタル分割)と命名されました。
この章ではこのフラクタル分割の理論について解説します。
指数作用素の分割のおさらい
非可換な作用素
(27)この指数の積の数
(28)
(29)です。
(30)です。
2次分割による3次分割の構成
まず
(31)作用素
(31)式の右辺の2つの指数のそれぞれを
(32)の項のことです。もしこの項が0になるのなら、それは(31)式の 2つの指数のそれぞれの
(32)式の項が0になる条件は明らかに
(33)です。このとき(31)式の2つの指数をそれぞれ2次の
結局、
(34)となります。これで3次の指数積分割はできあがりです。
m-1次分割によるm次分割の構成
一般に(31)式の2つの指数をそれぞれ
(35)こうして再帰的に無限の精度まで高次の指数積分解が可能なのです。 各指数の係数
第7章 指数作用素のフラクタル分割(後編)
より安定な分解法
先の方法により無限次のSymplectic Integratorが構築されるのですが 同じ次数のSymplectic Integratorによる計算でも、エネルギーの振動を より押えた安定なSymplectic Integratorを構築することができます。(31)式での
変形の仕方にはいろいろあるのですが、例えば(35)式を複雑化して
(36)と再帰定義する分解方法があります。この方法の明瞭な利点は
また分割の形が前後でSymmetric(対称)です。時間推進作用素であるこの 作用素がSymmetricであると、時間反転に対して対称となります。 つまり
さらに別な形式の変形として鈴木先生が推奨されている次の分割を 紹介します。
(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のうち奇数次の分解
(38)ここで
(38)式の
(39)(38)式を見直せば、
これにより(37)式のタイプの分解法は次のように拡張されます。
(40)フラクタル分解
(40)式で示される各次数の指数積展開での各項の展開係数をコンピュータで計算します。Symplectic Integratorでのこの展開係数の意味は、短い時間ステップ Δtうちの各係数の分を作用素
各係数の指数ひとつを作用する作業を1stepとして、 各stepで時間がどのように進んでいくかを以下の図に示します。
次数があがるとこのように折れ線グラフがより複雑になります。 その複雑さの入り組み方はまさに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();
}
