next up previous
Next: 4 多体問題 Up: 物理シミュレーションの手法と結果の検証 Previous: 2 シンプレクティック法

3 ケプラー問題

(とても長い時間のあと)

一応作ってみました。積分法は、折角なのでシンプレクティックというののほうです。こっちのほうがプログラム簡単ですし。

pt 

//
// kepler.C
//
// J. Makino 2002/06/03
#include  <stdlib.h>
#include  <math.h>
#include  <iostream>

using namespace std;

const int ndim=2;

int main()
{
    double x[ndim],v[ndim],t;
    double k;
    int norbits, nsteps, ndiag;
    cerr << "Enter norbits, nsteps, ndiag :";
    cin >> norbits >> nsteps >> ndiag;
    double dt = 2*M_PI/nsteps;
    cerr << "Enter vx,vy (normalized to 1):";
    cin >> v[0]>>v[1];
    double vabs = 0;
    for (int k = 0;k< ndim;k++) vabs +=v[k]*v[k];
    vabs = sqrt(vabs);
    for(int k = 0;k<ndim;k++) v[k]/= vabs;
    cerr << "v = " << v[0] << " " << v[1] << endl;
    x[0] = 1 ; x[1] = 0;
    t = 0;
    for(int i = 0 ; i<norbits*nsteps;i++){
        double a[ndim];
        double r2 = 0;
        for (int k = 0;k< ndim;k++) r2 += x[k]*x[k];
        double r = sqrt(r2);
        double r3inv = 1.0/(r*r2);
        t += dt;
        for (int k = 0;k< ndim;k++) a[k] = -x[k]*r3inv;
        for (int k = 0;k< ndim;k++) v[k] += a[k]*dt;
        for (int k = 0;k< ndim;k++) x[k] += v[k]*dt;
        //
        // calculate energy etc...
        //
        if ( ((i+1) % ndiag) == 0){
            double v2 =0;
            for (int k = 0;k< ndim;k++) v2 += v[k]*v[k];
            double etot = 0.5*v2 - 1.0/r;
            cout.precision(12);
            cout << "T, E,x = " << t <<" " << etot << " " << -0.5 -etot
                 << " " << x[0] << " " << x[1]<<endl;
        }
    }
}

とりあえず周期は$2\pi$ になるような初期条件から始めることにしてます。円軌道の場合で時間刻みをいろいろ変えてみたんですが、なんか調和振動の場合と変わんないみたいなんですけど、、、

赤木： 変わらないっていうと？

学生B：ええと、 $(1,0)$からの円軌道で、$x$と$y$の1周期後の誤差をみると やっぱり時間刻みの3乗や2乗になってます。

赤木： 円軌道でない場合はやってみた？

学生B：やってないですけど、おんなじじゃないんですか？

赤木： まあ、折角プログラム書いたんだから、確かめてみたら。

学生B：はあ、、、

(調べる)

あれ、ほんの少しでも円軌道からずれると、いきなり全然駄目っていうか、誤 差が時間刻みの1乗になっちゃいます。え、どうしてなんですか？

赤木： これもちゃんと説明すると長いんだけど、要するに円軌道ってい うのは特別な場合になってるのよ。と、これでは全然説明になってないわね、 まあちゃんと計算すればすぐわかるはずよ。つまり、円軌道の場合には本当の 解自体は調和振動と同じでしょ、で、誤差が小さくとすれば、誤差の従う差分 方程式を書き下せるわけね。これは線型の方程式になって、前の調和振動の場合と同じ振舞いをするわけ。

学生B：え、でもそれなんか変じゃないですか。誤差が小さければ線型化した 方程式を書けるのは円軌道でなくても同じじゃないですか。それに方程式自体 同じようなものにならないですか？

赤木： それはねえ、本当はそうなの。ではなにが違うかってのは、まあ 計算してみればすぐにわかるわ。こういうところの理屈がわかって、ちゃんと 計算できるのは大事なことよ。そこまでいちいち私が説明したんではしょうが ないでしょ。

学生B：ええと、それはそうなんですが、でも自分でやってできなかったら困 るじゃないですか。僕は後で赤木さんに聞けるけど、読者はどうするんですか？

赤木： あーもう、細かいわね。問い合わせが殺到するようなら作者が Web3 に説明をだすでしょ。それはいいんだけど、誤差が一周のあいだにどう変 わるかって調べてみた？例えばエネルギーとか。

学生B：みてないですけど、一次の公式なんだから単調に増えるか減るかするだけじゃないんですか？

赤木：ちょっとやってみてよ。

学生B：やってみました、、、え、一周期でエネルギーがもとに戻る、、、

赤木： そういうことね。これは調和振動の場合もエネルギーを書いてみ れば同じだったはずね。円軌道になるはずのが楕円になったんだから、周期的 にエネルギーが変わるわけ。もとが楕円のケプラー軌道なら、数値解もまあ楕 円みたいなものがでてきてその周期でエネルギーも変わるわけ。ただ、数値解 は厳密な楕円じゃないし、そもそももとのところに戻らないから、エネルギー が変わらなければいいというものでもないけど。

こういうふうになるというのがシンプレクティック法の特徴の1つなの。真の 解が周期解なら、数値解も周期解になるのね。もちろん数値誤差があるから同 じ周期解にはならないわけで、ずれの大きさが誤差で決まるわけ。

学生B：ということは、シンプレクティック法でも1次ではなくてもっと高次の 方法を使えば誤差を簡単に小さくできるわけですよね。そういう高次の方法っ ていうのはあるんですか？

赤木： もちろんいくらでもあるの。ここ15年くらいでいろいろ研究が進 んだし。 まず2次からいきましょう。2次はこういうの


$\displaystyle v_{i+1/2}$ $\textstyle =$ $\displaystyle v_i + \Delta t a(x_i)/2$ (15)
$\displaystyle x_{i+1}$ $\textstyle =$ $\displaystyle x_{i} + \Delta t v_{i+1/2}$ (16)
$\displaystyle v_{i+1}$ $\textstyle =$ $\displaystyle v_{i+1/2} + \Delta ta(x_{i+1})/2$ (17)

学生B：なんか変な格好ですね。 $v_{i+1/2}$ ってのは中間で定義されるってことですか？でも、それ消去して

$\displaystyle x_{i+1}$ $\textstyle =$ $\displaystyle x_{i} + \Delta t v_{i} + \Delta t^2 a(x_i)/2$ (18)
$\displaystyle v_{i+1}$ $\textstyle =$ $\displaystyle v_{i} + \Delta t [a(x_i)+a(x_{i+1})]/2$ (19)

としてもおんなじじゃないですか？

赤木： あ、これはもちろんおんなじよ。私の書いた形のほうが足し算一 つ少ないから速いかもしれないけど、普通は加速度計算するほうに時間かかる から関係ないわね。いじましいことをいえば、計算量を減らすには逆に $v_{i+1/2}$ だけを残して

$\displaystyle v_{i+1/2}$ $\textstyle =$ $\displaystyle v_{i-1/2} + \Delta t a(x_i)$ (20)
$\displaystyle x_{i+1}$ $\textstyle =$ $\displaystyle x_{i} + \Delta t v_{i+1/2}$ (21)

とすればいいわね。最初は上の式で $v_{i}$ から計算するわけ。この形のこ とをリープフロッグ公式ということが多いわね。ここで書いたのはみんな同じ だから、どれもリープフロッグと呼んで間違いではないけど。

学生B：あれ、これだと一次の公式とおんなじじゃないですか。どうしてこれで2次になるんですか？

赤木： 2番目の形をみれば2次でしょ、$x$ はテイラー展開の2次まで入っ てるし、 $v$ は台形公式だから。台形公式って知ってるわよね？ 17 式が定義。普通は陰公式で、 $i+1$ での値を求めるのに代 数方程式を解かないといけないけど、この場合は加速度が速度に依存しないか ら陽的に計算できるわけ。 最後の形だと $x$ と $v$ がずれて定義さ れてるのが仕掛けということになるわね。

学生B：えーと、式がそうだってのはわかるんですが、なんか変な気がするん ですが、、、 大体 $x$ と $v$ で式が違うのが変じゃないですか。

$\displaystyle x_{i+1}$ $\textstyle =$ $\displaystyle x_{i} + \Delta t (v_i+v_{i+1})/2$ (22)
$\displaystyle v_{i+1}$ $\textstyle =$ $\displaystyle v_{i} + \Delta t [a(x_i)+a(x_{i+1})]/2$ (23)

とは、、、あ、これだと陰公式になっちゃうんですね。 $x_{i+1}$ を計算す るのに $v_{i+1}$ が必要で、それには$a(x_{i+1})$ がいるのか。

赤木： そうね。台形公式はいろいろいい性質があるけど、普通に使うと さっきいったように陰公式で代数方程式を解く必要がでてくるわけ。リープフロッグは、形式的に は

$\displaystyle x_{i+1}$ $\textstyle =$ $\displaystyle x_{i} + \Delta t (v_i+v_{i+1})/2 - \Delta t^2 [a(x_{i+1})-a(x_{i})]/4$ (24)
$\displaystyle v_{i+1}$ $\textstyle =$ $\displaystyle v_{i} + \Delta t [a(x_i)+a(x_{i+1})]/2$ (25)

と書けるでしょ。$x$の最後の項は $O(\Delta t^3)$ だから、公式の次数には関係ないのよ。

学生B：ええと、それって変じゃないですか？ $[a(x_{i+1})-a(x_{i})]/\Delta t$ は加速度の微分の近似になってるから

\begin{displaymath} x_{i+1} = x_{i} + \Delta t v_{i} + \Delta t^2 a(x_i)/2 + \Delta t^2 [a(x_{i+1})-a(x_{i})]/6 \end{displaymath} (26)

としたら 3次になりませんか？陰公式だから計算面倒ですが、、、

赤木： あら、なかなか冴えてるじゃない。まあ、そのへんも、自分で調 べてみて納得しないと、本とか読んだだけじゃわからないのよね。というわけ でそれは宿題ということにして、リープフロッグでプログラム書いてみて。で、 メインプログラムだけってのでなくて、もうちょっと違う公式とかそ ういうのを使いやすいように書いてみてくれないかしら。

学生B：やってみます。
(無限とも思える時間の後) こんなのでどうでしょう？

pt 

//
// kepler-leap.C
//
// J. Makino 2002/06/05
#include  <stdlib.h>
#include  <math.h>
#include  <iostream>

using namespace std;

const int ndim=2;


void accel(double x[ndim],
           double v[ndim],
           double a[ndim])
{
    double r2 = 0;
    for (int k = 0;k< ndim;k++) r2 += x[k]*x[k];
    double r = sqrt(r2);
    double r3inv = 1.0/(r*r2);
    for (int k = 0;k< ndim;k++) a[k] = -x[k]*r3inv;
}    
void update(double y[ndim],
            double ydot[ndim],
            double dt)
{
    for (int k = 0;k< ndim;k++) y[k] += ydot[k]*dt;
}

void leapfrog(double x[ndim],
              double v[ndim],
              double a[ndim],
              double dt)
{
    double dthalf = dt/2;
    update(v,a,dthalf);
    update(x,v,dt);
    accel(x,v,a);
    update(v,a,dthalf);
}

void diag(double x[ndim],
          double v[ndim],
          double t)
{
    double v2 =0;
    double r = 0;
    for (int k = 0;k< ndim;k++) v2 += v[k]*v[k];
    for (int k = 0;k< ndim;k++) r += x[k]*x[k];
    r = sqrt(r);
    double etot = 0.5*v2 - 1.0/r;
    cout.precision(12);
    cout << "T, E,x = " << t <<" " << etot << " " << -0.5 -etot
         << " " << x[0] << " " << x[1]<<endl;
}    

int main()
{
    double x[ndim],v[ndim],t;
    double k;
    int norbits, nsteps, ndiag;
    cerr << "Enter norbits, nsteps, ndiag :";
    cin >> norbits >> nsteps >> ndiag;
    double dt = 2*M_PI/nsteps;
    cerr << "Enter vx,vy (normalized to 1):";
    cin >> v[0]>>v[1];
    double vabs = 0;
    for (int k = 0;k< ndim;k++) vabs +=v[k]*v[k];
    vabs = sqrt(vabs);
    for(int k = 0;k<ndim;k++) v[k]/= vabs;
    cerr << "v = " << v[0] << " " << v[1] << endl;
    x[0] = 1 ; x[1] = 0;
    t = 0;
    double a[ndim];
    accel(x,v,a);
    for(int i = 0 ; i<norbits*nsteps;i++){
        leapfrog(x,v,a,dt);
        t+= dt;
        //
        // calculate energy etc...
        //
        if ( ((i+1) % ndiag) == 0){
            diag(x,v,t);
        }
    }
}

赤木： 初期設定がいま一つな感じもするけど、まあいいんじゃないかし ら。精度とかは見た？

学生B：ええ、ちゃんとエネルギーとかが 2 次になってるのは確認しました。これよりもっと精度をあげようと思うとどんな方法があるんですか？

赤木： そうね、一番簡単なのは、リープフロッグの組合せで次数をあげ ていく方法ね。図2 みたいに、最初行き過ぎて、それからもどってまた行き過ぎ るというふうに3回リープフロッグを使って、誤差がうちけしあうようにする わけ。 進む量をどうすればいいかっていうと、1ステップでの誤差は $\Delta t^3$ に比例するわけだから、進む量と戻る量を $d_1\Delta t, d_2\Delta t$ とすれば、

$\displaystyle 2d_1 -d_2$ $\textstyle =$ $\displaystyle 1$ (27)
$\displaystyle 2d_1^3$ $\textstyle =$ $\displaystyle d_2^3$ (28)

を満たすように決めればいいわけ。

Figure 2: 4次シンプレクティック法
\begin{figure}\begin{center} \leavevmode \epsfxsize 6 cm \epsffile{yoshida4.eps}\end{center}\end{figure}

学生B：あ、なるほど、そうすると、プログラムは例えば 

void symp4(double x[ndim],
              double v[ndim],
              double a[ndim],
              double dt)
{
    const double d1 = 1.0/(2.0-pow(2.0,1.0/3));
    const double d2 = 1-2*d1;
    leapfrog(x,v,a,d1*dt);
    leapfrog(x,v,a,d2*dt);
    leapfrog(x,v,a,d1*dt);
}
というのでいいわけですね。

赤木： そういうことね。これは割合簡単でしょ？

学生B：そうですね。3体とか、粒子数が多いのをやってみたいんですが、、、 解析解がある場合の話だけではつまんないですし。

赤木： そうね、まあ、なんか書いてみたら？

学生B：はい、、、

next up previous
Next: 4 多体問題 Up: 物理シミュレーションの手法と結果の検証 Previous: 2 シンプレクティック法
Jun Makino
平成14年10月4日