数値積分の誤差の問題

質問

数値解析は初心者です。
惑星の運動を調べるために（３つの星です）、
微分方程式を数値的に解こうとしています。
エネルギーが保存する系では、
シンプレクティック法というもの方が良いと
教えてもらいいろいろな人のものを参考にして
使っています。

問題点は、シンプレクティック法の結果は、
ルンゲクッタよりエネルギーの保存の
精度はいいのですが
１～１０万ステップくらいの挙動が
計算するコンピューターで違うことが
分かりました。誤差が出ているのだと思います。

ルンゲクッタは保存は悪い代わりに
１０万ステップくらいまでマシンに依存した
結果はでません。

エネルギーは高い精度で保存しているので
バグではないと思います。

シンプレクティック法は６次と８次を使っています。

どなたかシンプレクティック法を使って
不安定な系を扱っている方で
シンプレクティック法の精度を上げる方法を
ご存知の方がいらっしゃったらお教えください。

回答No.1

tomtom_

情報が少ないので推測で書きます．
まず，計算機によって結果が異なるというのが，非常に小さな変化なら気にすることはないと思います．シンプクレクティック数値解法や多様体法，シンメトリック法はそんなもんです．

次に，6次や8次という高い精度をどうやって作っているかが気になります．もしガウス・ルジャンドル型なら，反復回数を打ち切るマシンイプシロンなんかがどうなってるかチェックしたいところです．
もし，対称分解法や対称合成法，フラクタル分解法を使っているなら，分解の仕方によって安定度が大きく変わります．例えば2次の解法から4次の解法を作る場合，

S4=S2(a τ)S2(a τ)

と分解するよりも，

S4=S2(a τ)S2(b τ)S2(a τ)

とか

S4=S2(a τ)S2(b τ)S2(c τ)S2(a τ)S2(b τ)

と分解するほうが，aやbなどの定数の値が小さくなって安定します．もちろん計算時間は増えます．
もしこの辺りのせいなら，刻み幅の値を変えて計算し，比べてみれば分かります．
刻み幅を半分にして，誤差が6次の方法で0.5^6に減っていなかったら，プログラムのバグです．
誤差がちゃんと減っているば，丸め誤差や打切り誤差のせいですから，安定した解法を使うことは意味があります．

お礼

初心者にも分かりやすく、
詳しいご説明大変ありがとうございます。
誤差について詳しくないのでだいぶわかりました。
ありがとうございます。

そんなものと考えればいいとのことで
少し安心しました。
ただルンゲクッタより計算機の依存性が大きかったのが
意外でしたので。

私がもらったシンプレクティック法のやり方は、
６次と８次はすでに２次で分解した状態に解いた
時の各２時の係数の表が入っていてそれを
読み込んで使っています。これは、吉田という人が
論文で発表してあるものだそうです。

やはり時間ステップを小さくすれば誤差は減るんですか。今ステップを０．０１単位時間にして、１万単位時間くらい計算してるのですが、ステップを０．００１にした場合、１万時間まで見るためには、１０倍のステップがいるために結局同じなのかなと勘違いしていました。

幅を変えるのをやってみます。
また報告します。

ありがとうございます。

補足

つかっているアルゴリズムは以下の通りです。
もっと良い分割の仕方や係数の表をご存じでしたらお教え頂ければと思います。

void symplectic_integrator_2(double *x, double dt)
{
double dxdt[2*DOF];
int i;

/* Step 1 */
time_derivatives(x, dxdt);
for (i = DOF; i < 2*DOF; i++)
x[i] += 0.5*dxdt[i]*dt;

/* Step 2 */
time_derivatives(x, dxdt);
for (i = 0; i < DOF; i++)
x[i] += dxdt[i]*dt;

/* Step 3 */
time_derivatives(x, dxdt);
for (i = DOF; i < 2*DOF; i++)
x[i] += 0.5*dxdt[i]*dt;
}

void symplectic_integrator_6(double *x, double dt)
{
int i;
double w0, w1, w2, w3;

w1 = -1.17767998417887;
w2 = 0.235573213359357;
w3 = 0.784513610477560;
w0 = 1.0 - 2.0 * (w1 + w2 + w3);

/* Step 1 */
symplectic_integrator_2(x, w3 * dt);

/* Step 2 */
symplectic_integrator_2(x, w2 * dt);

/* Step 3 */
symplectic_integrator_2(x, w1 * dt);

/* Step 4 */
symplectic_integrator_2(x, w0 * dt);

/* Step 5 */
symplectic_integrator_2(x, w1 * dt);

/* Step 6 */
symplectic_integrator_2(x, w2 * dt);

/* Step 7 */
symplectic_integrator_2(x, w3 * dt);

}
void symplectic_integrator_8(double *x, double dt)
{
int i;
double w0, w1, w2, w3, w4, w5, w6, w7;
w1 = 0.311790812418427;
w2 = -0.155946803821447*10.0;
w3 = -0.167896928259640*10.0;
w4 = 0.166335809963315*10.0;
w5 = -0.106458714789183*10.0;
w6 = 0.136934946416871*10.0;
w7 = 0.629030650210433;
w0 = 1. - 2. * (w1 + w2 + w3 + w4 + w5 + w6 + w7);

symplectic_integrator_2(x, w7 * dt);
symplectic_integrator_2(x, w6 * dt);
symplectic_integrator_2(x, w5 * dt);
symplectic_integrator_2(x, w4 * dt);
symplectic_integrator_2(x, w3 * dt);
symplectic_integrator_2(x, w2 * dt);
symplectic_integrator_2(x, w1 * dt);
symplectic_integrator_2(x, w0 * dt);
symplectic_integrator_2(x, w1 * dt);
symplectic_integrator_2(x, w2 * dt);
symplectic_integrator_2(x, w3 * dt);
symplectic_integrator_2(x, w4 * dt);
symplectic_integrator_2(x, w5 * dt);
symplectic_integrator_2(x, w6 * dt);
symplectic_integrator_2(x, w7 * dt);
}