38 lines
1.1 KiB
D
38 lines
1.1 KiB
D
import std.stdio, std.math, std.typecons;
|
|
|
|
alias FP = real;
|
|
alias FPs = Typedef!(FP[101]);
|
|
|
|
void runge(in FP function(in FP, in FP)
|
|
pure nothrow @safe @nogc yp_func,
|
|
ref FPs t, ref FPs y, in FP dt) pure nothrow @safe @nogc {
|
|
foreach (immutable n; 0 .. t.length - 1) {
|
|
immutable FP
|
|
dy1 = dt * yp_func(t[n], y[n]),
|
|
dy2 = dt * yp_func(t[n] + dt / 2.0, y[n] + dy1 / 2.0),
|
|
dy3 = dt * yp_func(t[n] + dt / 2.0, y[n] + dy2 / 2.0),
|
|
dy4 = dt * yp_func(t[n] + dt, y[n] + dy3);
|
|
t[n + 1] = t[n] + dt;
|
|
y[n + 1] = y[n] + (dy1 + 2.0 * (dy2 + dy3) + dy4) / 6.0;
|
|
}
|
|
}
|
|
|
|
FP calc_err(in FP t, in FP calc) pure nothrow @safe @nogc {
|
|
immutable FP actual = (t ^^ 2 + 4.0) ^^ 2 / 16.0;
|
|
return abs(actual - calc);
|
|
}
|
|
|
|
void main() {
|
|
enum FP dt = 0.10;
|
|
FPs t_arr, y_arr;
|
|
|
|
t_arr[0] = 0.0;
|
|
y_arr[0] = 1.0;
|
|
runge((t, y) => t * y.sqrt, t_arr, y_arr, dt);
|
|
|
|
foreach (immutable i; 0 .. t_arr.length)
|
|
if (i % 10 == 0)
|
|
writefln("y(%.1f) = %.8f Error: %.6g",
|
|
t_arr[i], y_arr[i],
|
|
calc_err(t_arr[i], y_arr[i]));
|
|
}
|