radau5_step_audit/
radau5_step_audit.rs

1//! Step-count audit for the corrected Radau5 implementation.
2//!
3//! Prints `n_accept / n_reject / n_eval / n_jac / n_lu` for the five reference
4//! problems used to validate Radau5 against SciPy's port of Hairer's radau5.f.
5
6use numra_ode::problem::OdeProblem;
7use numra_ode::radau5::Radau5;
8use numra_ode::solver::{Solver, SolverOptions};
9
10fn run<F>(name: &str, rhs: F, t0: f64, tf: f64, y0: Vec<f64>, rtol: f64, atol: f64)
11where
12    F: Fn(f64, &[f64], &mut [f64]) + Copy,
13{
14    let problem = OdeProblem::new(rhs, t0, tf, y0.clone());
15    let opts = SolverOptions::default().rtol(rtol).atol(atol);
16    let res = Radau5::solve(&problem, t0, tf, &y0, &opts).expect("solve failed");
17    println!(
18        "{:<32} acc={:>4}  rej={:>3}  rhs={:>5}  jac={:>3}  lu={:>4}",
19        name,
20        res.stats.n_accept,
21        res.stats.n_reject,
22        res.stats.n_eval,
23        res.stats.n_jac,
24        res.stats.n_lu,
25    );
26}
27
28fn main() {
29    println!("Radau5 step-count audit (after corrections)");
30    println!("--------------------------------------------");
31
32    run(
33        "stiff decay (rtol=1e-2)",
34        |_, y, dy| dy[0] = -100.0 * y[0],
35        0.0,
36        0.1,
37        vec![1.0],
38        1e-2,
39        1e-4,
40    );
41
42    run(
43        "exponential (rtol=1e-6)",
44        |_, y, dy| dy[0] = y[0],
45        0.0,
46        1.0,
47        vec![1.0],
48        1e-6,
49        1e-8,
50    );
51
52    run(
53        "linear 2D (rtol=1e-4)",
54        |_, y, dy| {
55            dy[0] = -y[0] + y[1];
56            dy[1] = -y[0] - y[1];
57        },
58        0.0,
59        1.0,
60        vec![1.0, 0.0],
61        1e-4,
62        1e-6,
63    );
64
65    let mu_mild = 10.0;
66    run(
67        "van der Pol mu=10 (rtol=1e-4)",
68        move |_, y, dy| {
69            dy[0] = y[1];
70            dy[1] = mu_mild * (1.0 - y[0] * y[0]) * y[1] - y[0];
71        },
72        0.0,
73        2.0,
74        vec![2.0, 0.0],
75        1e-4,
76        1e-6,
77    );
78
79    let mu_stiff = 100.0;
80    run(
81        "van der Pol mu=100 (rtol=1e-3)",
82        move |_, y, dy| {
83            dy[0] = y[1];
84            dy[1] = mu_stiff * (1.0 - y[0] * y[0]) * y[1] - y[0];
85        },
86        0.0,
87        20.0,
88        vec![2.0, 0.0],
89        1e-3,
90        1e-5,
91    );
92}
radau5_step_audit/radau5_step_audit.rs

radau5_step_audit/
radau5_step_audit.rs
