scivex_optim/ode/
euler.rs

1//! Forward Euler method for ODE initial value problems.
2//!
3//! The simplest first-order explicit method: `y_{n+1} = y_n + h * f(t_n, y_n)`.
4//! First-order accurate (`O(h)`), not recommended for stiff problems.
5
6use scivex_core::Float;
7
8use super::{OdeOptions, OdeResult};
9use crate::error::Result;
10
11/// Solve an ODE system using the forward Euler method.
12///
13/// `f(t, y) -> dy/dt` defines the system. `y0` is the initial state vector.
14///
15/// # Examples
16///
17/// ```
18/// # use scivex_optim::ode::{euler, OdeOptions};
19/// // dy/dt = -y, y(0) = 1 → y(t) = e^(-t)
20/// let result = euler(
21///     |_t: f64, y: &[f64]| vec![-y[0]],
22///     [0.0, 1.0],
23///     &[1.0],
24///     &OdeOptions::default(),
25/// ).unwrap();
26/// let y_final = result.y.last().unwrap()[0];
27/// assert!((y_final - (-1.0_f64).exp()).abs() < 0.02);
28/// ```
29pub fn euler<T, F>(f: F, t_span: [T; 2], y0: &[T], options: &OdeOptions<T>) -> Result<OdeResult<T>>
30where
31    T: Float,
32    F: Fn(T, &[T]) -> Vec<T>,
33{
34    let t0 = t_span[0];
35    let tf = t_span[1];
36    let h = options
37        .first_step
38        .unwrap_or_else(|| (tf - t0) / T::from_f64(100.0));
39    let n = y0.len();
40
41    let mut t = t0;
42    let mut y = y0.to_vec();
43    let mut t_values = vec![t];
44    let mut y_values = vec![y.clone()];
45    let mut n_evals: usize = 0;
46
47    while t < tf {
48        // Don't step past tf
49        let step = if t + h > tf { tf - t } else { h };
50
51        let dy = f(t, &y);
52        n_evals += 1;
53
54        for i in 0..n {
55            y[i] += step * dy[i];
56        }
57        t += step;
58
59        t_values.push(t);
60        y_values.push(y.clone());
61
62        // Event detection
63        if let Some(ref event_fn) = options.event_fn {
64            let val = event_fn(t, &y);
65            if val.abs() < T::from_f64(1e-12)
66                || (t_values.len() > 1 && {
67                    let prev_y = &y_values[y_values.len() - 2];
68                    let prev_t = t_values[t_values.len() - 2];
69                    let prev_val = event_fn(prev_t, prev_y);
70                    (prev_val > T::zero()) != (val > T::zero())
71                })
72            {
73                break;
74            }
75        }
76    }
77
78    let n_steps = y_values.len() - 1;
79    Ok(OdeResult {
80        t: t_values,
81        y: y_values,
82        n_evals,
83        n_steps,
84        success: true,
85    })
86}
87
88#[cfg(test)]
89mod tests {
90    use super::*;
91
92    #[test]
93    fn test_euler_exponential_decay() {
94        // dy/dt = -y, y(0) = 1 => y(t) = e^(-t)
95        let result = euler(
96            |_t: f64, y: &[f64]| vec![-y[0]],
97            [0.0, 1.0],
98            &[1.0],
99            &OdeOptions::default(),
100        )
101        .unwrap();
102
103        let y_final = result.y.last().unwrap()[0];
104        let expected = (-1.0_f64).exp();
105        // Euler is only first-order, so allow some error
106        assert!(
107            (y_final - expected).abs() < 0.02,
108            "y_final={y_final}, expected={expected}"
109        );
110    }
111
112    #[test]
113    fn test_euler_linear() {
114        // dy/dt = 1, y(0) = 0 => y(t) = t
115        let result = euler(
116            |_t: f64, _y: &[f64]| vec![1.0],
117            [0.0, 2.0],
118            &[0.0],
119            &OdeOptions::default(),
120        )
121        .unwrap();
122
123        let y_final = result.y.last().unwrap()[0];
124        assert!(
125            (y_final - 2.0).abs() < 1e-10,
126            "y_final={y_final}, expected=2.0"
127        );
128    }
129
130    #[test]
131    fn test_euler_system() {
132        // dy0/dt = y1, dy1/dt = -y0 (harmonic oscillator)
133        // y(0) = [1, 0] => y(t) = [cos(t), -sin(t)]
134        let opts = OdeOptions {
135            first_step: Some(0.001),
136            ..OdeOptions::default()
137        };
138        let result = euler(
139            |_t: f64, y: &[f64]| vec![y[1], -y[0]],
140            [0.0, 1.0],
141            &[1.0, 0.0],
142            &opts,
143        )
144        .unwrap();
145
146        let y_final = &result.y.last().unwrap();
147        let expected_y0 = 1.0_f64.cos();
148        let expected_y1 = -(1.0_f64.sin());
149        assert!(
150            (y_final[0] - expected_y0).abs() < 0.01,
151            "y0={}, expected={}",
152            y_final[0],
153            expected_y0
154        );
155        assert!(
156            (y_final[1] - expected_y1).abs() < 0.01,
157            "y1={}, expected={}",
158            y_final[1],
159            expected_y1
160        );
161    }
162
163    #[test]
164    fn test_euler_stores_trajectory() {
165        let result = euler(
166            |_t: f64, _y: &[f64]| vec![1.0],
167            [0.0, 1.0],
168            &[0.0],
169            &OdeOptions::default(),
170        )
171        .unwrap();
172
173        assert!(result.t.len() > 2);
174        assert_eq!(result.t.len(), result.y.len());
175        assert!((result.t[0] - 0.0).abs() < 1e-12);
176        assert!((*result.t.last().unwrap() - 1.0).abs() < 1e-12);
177    }
178}
scivex_optim/ode/euler.rs

scivex_optim/ode/
euler.rs
