さて、差分スキーム自体の話にいく前に、さっきの精度打ち止めの理由が情報落ちであることを確認しておく。
(用語訂正:さっきは桁落ちといったがふつうの用語ではなかった)。
情報落ちとは、指数部が大きく異なる浮動小数点数を加減算すると、増分の精度が数形式の精度より大きく劣ることである。
桁落ちとは、絶対値がほとんど同じで、なおかつ符号の反対な浮動小数点数を加算(符号が同じ数を減算といっても同じこと)すると、結果の相対精度が大きく劣ることで、ちょっと違う。
ステップが非常に小さい積分は、この情報落ちにぴったり当てはまる。
回避法は簡単で、積算時に中間積算を行って、あるていどまとまったところで積算変数への足しこみを行えばよい。次の例では1024回づつ中間積算することによって、ステップが 1e-9 でも精度打ち止めにならないよう
になった。(おそらくは精度を3桁改善できるだろうが、ものすごく時間がかかるので確認しない)
0.0001000000000000 0.9998999999999811 -0.0001000000000189
0.0000300000000000 0.9999700002000160 -0.0000299997999840
0.0000100000000000 0.9999900000000238 -0.0000099999999762
0.0000030000000000 0.9999970000020161 -0.0000029999979839
0.0000010000000000 0.9999989999999600 -0.0000010000000400
0.0000003000000000 0.9999996999999685 -0.0000003000000315
0.0000001000000000 0.9999998999999479 -0.0000001000000521
0.0000000300000000 0.9999999700006352 -0.0000000299993648
0.0000000100000000 0.9999999900010098 -0.0000000099989902
0.0000000030000000 0.9999999970078517 -0.0000000029921483
0.0000000010000000 0.9999999990116513 -0.0000000009883487
import java.text.DecimalFormat;
public strictfp class Z {
public static void main(String[] args) {
DecimalFormat df = new DecimalFormat(" 0.0;-0.0");
df.setMinimumIntegerDigits(1);
df.setMaximumIntegerDigits(1);
df.setMinimumFractionDigits(16);
df.setMaximumFractionDigits(16);
double dxs[] = {1.0e-4, 3.0e-5, 1.0e-5, 3.0e-6, 1.0e-6, 3.0e-7,
1.0e-7, 3.0e-8, 1.0e-8, 3.0e-9, 1.0e-9};
double xmin = 0.0;
double xmax = 1.0;
double ymin = 0.0;
int midlimit = 1024;
for (double step: dxs) {
double dx = step;
double xstop = xmax - dx;
double x, y, xmid, ymid, dydx;
boolean continuep;
int count;
for (x = xmin, y = ymin, xmid = 0.0, ymid = 0.0, continuep = true,
count = 0; continuep; ) {
dydx = (x + xmid) * 2.0;
if (x + xmid > xstop) {
dx = xmax - (x + xmid);
count = midlimit;
continuep = false;
}
ymid += dx * dydx;
xmid += dx;
if (++count > midlimit) {
x += xmid;
y += ymid;
xmid = ymid = 0.0;
count = 0;
}
}
System.out.println(
df.format(step)
+ " " + df.format(y)
+ " " + df.format(y - (xmax * xmax - xmin * xmin))
);
}
}
}
0 件のコメント:
コメントを投稿