被積分関数をかえてみる

簡単な例なので被積分関数をいじってみる。dy/dx = f(x, y) = 2x というのは y に依存しないので簡単すぎたのだ。ならば、 f(x,y) = x + sqrt(y) とやっても解析的には同じことである
が、 y がずれると積分結果がずれていってしまうという意味でよい感じになる。

次の結果で上が４次ルンゲクッタ、下が前進差分。もちろんバイアスがないルンゲクッタで定数倍改善していることも注目に値するが、ステップを 1/10 にしたときに１桁以上精度が改善しているこ
とに注目される。

0.3000000000000000 0.9717281348583426 -0.0282718651416574
0.1000000000000000 0.9945397310971597 -0.0054602689028403
0.0300000000000000 0.9991023660366382 -0.0008976339633618
0.0100000000000000 0.9998272390531332 -0.0001727609468668
0.0030000000000000 0.9999716121300442 -0.0000283878699558
0.0010000000000000 0.9999945367413657 -0.0000054632586343
0.0003000000000000 0.9999991022947341 -0.0000008977052659
0.0001000000000000 0.9999998272366565 -0.0000001727633435
0.0000300000000000 0.9999999716130413 -0.0000000283869587
0.0000100000000000 0.9999999945402596 -0.0000000054597404
0.0000030000000000 0.9999999991117637 -0.0000000008882363
0.0000010000000000 0.9999999998104413 -0.0000000001895587
0.0000003000000000 1.0000000000950195 0.0000000000950195
0.0000001000000000 1.0000000004349120 0.0000000004349121
0.0000000300000000 0.9999999989783790 -0.0000000010216210
0.0000000100000000 0.9999999959029633 -0.0000000040970367

0.3000000000000000 0.5100000000000001 -0.4899999999999999
0.1000000000000000 0.8100000000000004 -0.1899999999999996
0.0300000000000000 0.9410999999999989 -0.0589000000000011
0.0100000000000000 0.9800999999999989 -0.0199000000000011
0.0030000000000000 0.9940109999999988 -0.0059890000000012
0.0010000000000000 0.9980009999999985 -0.0019990000000015
0.0003000000000000 0.9994001100000743 -0.0005998899999257
0.0001000000000000 0.9998000100001390 -0.0001999899998610
0.0000300000000000 0.9999400011009805 -0.0000599988990195
0.0000100000000000 0.9999800001035221 -0.0000199998964779
0.0000030000000000 0.9999940000204693 -0.0000059999795307
0.0000010000000000 0.9999979999842054 -0.0000020000157946
0.0000003000000000 0.9999994001235168 -0.0000005998764832
0.0000001000000000 0.9999998004403844 -0.0000001995596156
0.0000000300000000 0.9999999389792776 -0.0000000610207224
0.0000000100000000 0.9999999759031382 -0.0000000240968618

import java.text.DecimalFormat;

public strictfp class Z {

static double gradient(double x, double y) {
return x + Math.sqrt(y);
}

public static double rk4(double xmin, double xmax, double ymin,
double step) {
double dx = step;
double xstop = xmax - dx;
double inv6 = 1.0/6.0;
double k1, k2, k3, k4;
double x = xmin;
double y = ymin;
boolean continuep = true;
while (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;
}
return y;
}

public static double euler(double xmin, double xmax, double ymin,
double step) {
double dx = step;
double xstop = xmax - dx;
double dydx;
double x = xmin;
double y = ymin;
boolean continuep = true;
while (continuep) {
if (x > xstop) {
dx = xmax - x;
continuep = false;
}
dydx = gradient(x, y);
y += dydx * dx;
x += dx;
}
return y;
}

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;
for (double step: dxs) {
double y = rk4(xmin, xmax, ymin, step);
System.out.println(
df.format(step)
+ " " + df.format(y)
+ " " + df.format(y - (xmax * xmax - xmin * xmin))
);
}
System.out.println("");
for (double step: dxs) {
double y = euler(xmin, xmax, ymin, step);
System.out.println(
df.format(step)
+ " " + df.format(y)
+ " " + df.format(y - (xmax * xmax - xmin * xmin))
);
}
}
}