4次ルンゲクッタの威力（？）

さて、4次ルンゲクッタにしたところ、困ったことに誤差がなくなってしまった。
前進差分と違って、バイアスがないし、高次なものだから２次関数のシミュレートなんて厳密にできてしまう。刻みが 0.3 ですら解析解と全ビット同じ値をたたき出すのである。結果
として、ステップが小さくなるにつれて大きくなる情報落ちだけが誤差となる。

教材としては適当ではないが、情報落ちの可視化としては適当なのでメモしておく。

結果

ステップ 積分結果 誤差（対解析解）
0.3000000000000000 1.0000000000000000 0.0000000000000000
0.1000000000000000 1.0000000000000002 0.0000000000000002
0.0300000000000000 0.9999999999999991 -0.0000000000000009
0.0100000000000000 0.9999999999999990 -0.0000000000000010
0.0030000000000000 0.9999999999999990 -0.0000000000000010
0.0010000000000000 0.9999999999999988 -0.0000000000000012
0.0003000000000000 1.0000000000000695 0.0000000000000695
0.0001000000000000 1.0000000000001080 0.0000000000001079
0.0000300000000000 1.0000000000007920 0.0000000000007920
0.0000100000000000 1.0000000000031826 0.0000000000031826
0.0000030000000000 1.0000000000076747 0.0000000000076747
0.0000010000000000 0.9999999999833878 -0.0000000000166122
0.0000003000000000 1.0000000001106562 0.0000000001106562
0.0000001000000000 1.0000000003876117 0.0000000003876117
0.0000000300000000 0.9999999990878823 -0.0000000009121177
0.0000000100000000 0.9999999963603514 -0.0000000036396486

コード

import java.text.DecimalFormat;

public strictfp class Z {

static double gradient(double x, double y) {
return x * 2.0;
}

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[] = {3.0e-1, 1.0e-1, 3.0e-2, 1.0e-2, 3.0e-3, 1.0e-3, 3.0e-4,
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};
double xmin = 0.0;
double xmax = 1.0;
double ymin = 0.0;
int midlimit = 1;
for (double step: dxs) {
double dx = step;
double xstop = xmax - dx;
double x, y, k1, k2, k3, k4;
double inv6 = 1.0/6.0;
boolean continuep;
for (x = xmin, y = ymin, continuep = true; continuep; ) {
if (x > xstop) {
dx = xmax - x;
continuep = false;
}
k1 = gradient(x, y);
k2 = gradient(x + 0.5 * dx, y + 0.5 * dx * k1);
k3 = gradient(x + 0.5 * dx, y + 0.5 * dx * k2);
k4 = gradient(x + dx, y + dx * k3);
y += dx * (k1 + 2.0 * (k2 + k3) + k4) * inv6;
x += dx;
}
System.out.println(
df.format(step)
+ " " + df.format(y)
+ " " + df.format(y - (xmax * xmax - xmin * xmin))
);
}
}
}