41 lines
1 KiB
Text
41 lines
1 KiB
Text
use fmt;
|
|
use math;
|
|
|
|
export fn main() void = {
|
|
rk4_driver(&f, 0.0, 10.0, 1.0, 0.1);
|
|
};
|
|
|
|
fn rk4_driver(func: *fn(_: f64, _: f64) f64, t_init: f64, t_final: f64, y_init: f64, h: f64) void = {
|
|
let n = ((t_final - t_init) / h): int;
|
|
let tn: f64 = t_init;
|
|
let yn: f64 = y_init;
|
|
let i: int = 1;
|
|
|
|
fmt::printfln("{: 2} {: 18} {: 21}", "t", "y(t)", "absolute error")!;
|
|
fmt::printfln("{: 2} {: 18} {: 21}", tn, yn, math::absf64(exact(tn) - yn))!;
|
|
|
|
for (i <= n; i += 1) {
|
|
yn = rk4(func, tn, yn, h);
|
|
tn = t_init + (i: f64)*h;
|
|
|
|
if (i % 10 == 0) {
|
|
fmt::printfln("{: 2} {: 18} {: 21}\t", tn, yn, math::absf64(exact(tn) - yn))!;
|
|
};
|
|
};
|
|
};
|
|
|
|
fn rk4(func: *fn(_: f64, _: f64) f64, t: f64, y: f64, h: f64) f64 = {
|
|
const k1 = func(t, y);
|
|
const k2 = func(t + 0.5*h, y + 0.5*h*k1);
|
|
const k3 = func(t + 0.5*h, y + 0.5*h*k2);
|
|
const k4 = func(t + h, y + h*k3);
|
|
return y + h/6.0 * (k1 + 2.0*k2 + 2.0*k3 + k4);
|
|
};
|
|
|
|
fn f(t: f64, y: f64) f64 = {
|
|
return t * math::sqrtf64(y);
|
|
};
|
|
|
|
fn exact(t: f64) f64 = {
|
|
return 1.0/16.0 * math::powf64(t*t + 4.0, 2.0);
|
|
};
|