numerics_ode/
adams_bashforth.rs

1//! Adams-Bashforth 2-step explicit multistep method.
2//!
3//! Requires a bootstrap step (computed with RK4) to start.
4//!
5//! ```text
6//! y_{n+1} = y_n + (h/2)(3·f(x_n, y_n) - f(x_{n-1}, y_{n-1}))
7//! ```
8//!
9//! Local truncation error: O(h³), global: O(h²).
10
11/// Solve an IVP using the 2-step Adams-Bashforth method.
12///
13/// The first step is bootstrapped with a single RK4 step.
14///
15/// # Arguments
16///
17/// * `f`     — Right-hand side.
18/// * `x0`    — Initial x.
19/// * `y0`    — Initial y.
20/// * `x_end` — Terminal x.
21/// * `n`     — Number of steps (≥ 2 for AB2 to be used; if `n == 1`, falls back to Euler).
22///
23/// # Returns
24///
25/// `(xs, ys)` — all computed points.
26pub fn solve(f: &dyn Fn(f64, f64) -> f64, x0: f64, y0: f64, x_end: f64, n: usize) -> (Vec<f64>, Vec<f64>) {
27    assert!(n >= 1, "number of steps must be ≥ 1");
28    let h = (x_end - x0) / n as f64;
29    let mut xs = Vec::with_capacity(n + 1);
30    let mut ys = Vec::with_capacity(n + 1);
31    xs.push(x0);
32    ys.push(y0);
33
34    if n == 1 {
35        // fallback: single Euler step
36        let y1 = y0 + h * f(x0, y0);
37        xs.push(x0 + h);
38        ys.push(y1);
39        return (xs, ys);
40    }
41
42    // Bootstrap: use RK4 for the first step
43    let x0_val = x0;
44    let y0_val = y0;
45    let k1 = f(x0_val, y0_val);
46    let k2 = f(x0_val + h / 2.0, y0_val + h * k1 / 2.0);
47    let k3 = f(x0_val + h / 2.0, y0_val + h * k2 / 2.0);
48    let k4 = f(x0_val + h, y0_val + h * k3);
49    let y1 = y0_val + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
50    let x1 = x0_val + h;
51    xs.push(x1);
52    ys.push(y1);
53
54    let mut f_prev = f(x0_val, y0_val);
55    let mut x_curr = x1;
56    let mut y_curr = y1;
57
58    for _ in 1..n {
59        let f_curr = f(x_curr, y_curr);
60        let y_next = y_curr + (h / 2.0) * (3.0 * f_curr - f_prev);
61        x_curr += h;
62        xs.push(x_curr);
63        ys.push(y_next);
64        f_prev = f_curr;
65        y_curr = y_next;
66    }
67
68    (xs, ys)
69}
70
71#[cfg(test)]
72mod tests {
73    use super::*;
74
75    #[test]
76    fn exponential_growth() {
77        let f = |_x: f64, y: f64| y;
78        let (_, ys) = solve(&f, 0.0, 1.0, 1.0, 100);
79        let exact = 1.0_f64.exp();
80        assert!((ys[100] - exact).abs() < 5e-4);
81    }
82
83    #[test]
84    fn convergence_order() {
85        let f = |_x: f64, y: f64| y;
86        let exact = 1.0_f64.exp();
87        let e1 = (solve(&f, 0.0, 1.0, 1.0, 50).1[50] - exact).abs();
88        let e2 = (solve(&f, 0.0, 1.0, 1.0, 100).1[100] - exact).abs();
89        let ratio = e1 / e2;
90        // AB2 is second order: ratio should be ~4 (2²)
91        assert!(ratio > 2.5 && ratio < 7.0, "ratio = {ratio}");
92    }
93
94    #[test]
95    fn constant_rhs() {
96        let f = |_x, _y| 0.0;
97        let (_, ys) = solve(&f, 0.0, 10.0, 1.0, 10);
98        for &y in &ys {
99            assert!((y - 10.0).abs() < 1e-12);
100        }
101    }
102
103    #[test]
104    fn linear_rhs() {
105        let f = |_x, _y| 3.0;
106        let (_, ys) = solve(&f, 0.0, 0.0, 2.0, 100);
107        assert!((ys[100] - 6.0).abs() < 1e-8);
108    }
109
110    #[test]
111    fn sinusoidal_rhs() {
112        let f = |x: f64, _y: f64| x.cos();
113        let (_, ys) = solve(&f, 0.0, 0.0, std::f64::consts::PI / 2.0, 200);
114        assert!((ys[200] - 1.0).abs() < 1e-4);
115    }
116
117    #[test]
118    fn exponential_decay() {
119        let f = |_x: f64, y: f64| -y;
120        let (_, ys) = solve(&f, 0.0, 1.0, 1.0, 200);
121        assert!((ys[200] - (-1.0_f64).exp()).abs() < 1e-5);
122    }
123
124    #[test]
125    fn single_step_fallback() {
126        let f = |_x, y| y;
127        let (xs, ys) = solve(&f, 0.0, 1.0, 1.0, 1);
128        assert_eq!(xs.len(), 2);
129        assert_eq!(ys.len(), 2);
130    }
131
132    #[test]
133    #[should_panic(expected = "number of steps must be ≥ 1")]
134    fn panics_on_zero_steps() {
135        let f = |_x, _y| 0.0;
136        solve(&f, 0.0, 1.0, 1.0, 0);
137    }
138
139    #[test]
140    fn quadratic_rhs() {
141        let f = |x: f64, _y: f64| 2.0 * x;
142        let (_, ys) = solve(&f, 0.0, 0.0, 3.0, 200);
143        assert!((ys[200] - 9.0).abs() < 1e-4);
144    }
145
146    #[test]
147    fn negative_direction() {
148        let f = |_x, y| y;
149        let (_, ys) = solve(&f, 1.0, 1.0_f64.exp(), 0.0, 200);
150        assert!((ys[200] - 1.0).abs() < 1e-3);
151    }
152}
numerics_ode/adams_bashforth.rs

numerics_ode/
adams_bashforth.rs
