mathhook_core/calculus/ode/numerical/
euler.rs

1//! Euler's method for numerical ODE solving
2//!
3//! Implements the forward Euler method for first-order ODEs: y' = f(x, y)
4//! Formula: y_{n+1} = y_n + h * f(x_n, y_n)
5
6use crate::calculus::ode::first_order::ODEError;
7
8/// Solves a first-order ODE using forward Euler method
9///
10/// # Arguments
11///
12/// * `f` - The derivative function f(x, y) where dy/dx = f(x, y)
13/// * `x0` - Initial x value
14/// * `y0` - Initial y value
15/// * `x_end` - Final x value
16/// * `step` - Step size h
17///
18/// # Returns
19///
20/// Vector of (x, y) solution points
21///
22/// # Examples
23///
24/// ```
25/// use mathhook_core::calculus::ode::numerical::euler::euler_method;
26///
27/// let solution = euler_method(
28///     |x, _y| x,
29///     0.0,
30///     0.0,
31///     1.0,
32///     0.1
33/// );
34///
35/// assert!(solution.len() > 0);
36/// assert_eq!(solution[0], (0.0, 0.0));
37/// ```
38pub fn euler_method<F>(f: F, x0: f64, y0: f64, x_end: f64, step: f64) -> Vec<(f64, f64)>
39where
40    F: Fn(f64, f64) -> f64,
41{
42    if step <= 0.0 {
43        return vec![(x0, y0)];
44    }
45
46    let num_steps = ((x_end - x0) / step).abs().ceil() as usize;
47    let mut solution = Vec::with_capacity(num_steps + 1);
48
49    let mut x = x0;
50    let mut y = y0;
51    solution.push((x, y));
52
53    let h = if x_end > x0 { step } else { -step };
54
55    for _ in 0..num_steps {
56        if (x_end > x0 && x >= x_end) || (x_end < x0 && x <= x_end) {
57            break;
58        }
59
60        let slope = f(x, y);
61        y += h * slope;
62        x += h;
63
64        solution.push((x, y));
65    }
66
67    solution
68}
69
70/// Solves a first-order ODE using Euler method with Result type
71///
72/// # Arguments
73///
74/// * `f` - The derivative function f(x, y)
75/// * `x0` - Initial x value
76/// * `y0` - Initial y value
77/// * `x_end` - Final x value
78/// * `step` - Step size
79///
80/// # Returns
81///
82/// Result containing vector of (x, y) solution points
83pub fn solve_euler<F>(
84    f: F,
85    x0: f64,
86    y0: f64,
87    x_end: f64,
88    step: f64,
89) -> Result<Vec<(f64, f64)>, ODEError>
90where
91    F: Fn(f64, f64) -> f64,
92{
93    if step <= 0.0 {
94        return Err(ODEError::InvalidInput {
95            message: "Step size must be positive".to_owned(),
96        });
97    }
98
99    if !x0.is_finite() || !y0.is_finite() || !x_end.is_finite() {
100        return Err(ODEError::InvalidInput {
101            message: "Initial values and endpoints must be finite".to_owned(),
102        });
103    }
104
105    Ok(euler_method(f, x0, y0, x_end, step))
106}
107
108#[cfg(test)]
109mod tests {
110    use super::*;
111
112    #[test]
113    fn test_euler_constant_derivative() {
114        let solution = euler_method(|_x, _y| 2.0, 0.0, 0.0, 1.0, 0.1);
115
116        assert_eq!(solution.len(), 11);
117        assert_eq!(solution[0], (0.0, 0.0));
118
119        let (x_final, y_final) = solution.last().unwrap();
120        assert!((x_final - 1.0).abs() < 1e-10);
121        assert!((y_final - 2.0).abs() < 0.1);
122    }
123
124    #[test]
125    fn test_euler_linear_ode() {
126        let solution = euler_method(|x, _y| x, 0.0, 0.0, 1.0, 0.1);
127
128        let (x_final, y_final) = solution.last().unwrap();
129        assert!((x_final - 1.0).abs() < 1e-10);
130        assert!((y_final - 0.5).abs() < 0.1);
131    }
132
133    #[test]
134    fn test_euler_exponential_growth() {
135        let solution = euler_method(|_x, y| y, 0.0, 1.0, 1.0, 0.1);
136
137        let (x_final, y_final) = solution.last().unwrap();
138        assert!((x_final - 1.0).abs() < 1e-10);
139
140        let expected = 1.0_f64.exp();
141        let relative_error = (y_final - expected).abs() / expected;
142        assert!(relative_error < 0.1);
143    }
144
145    #[test]
146    fn test_euler_backward_integration() {
147        let solution = euler_method(|x, _y| x, 1.0, 0.5, 0.0, 0.1);
148
149        assert!(solution.len() > 1);
150        let (x_first, y_first) = solution[0];
151        let (x_final, _y_final) = solution.last().unwrap();
152
153        assert_eq!((x_first, y_first), (1.0, 0.5));
154        assert!((x_final - 0.0).abs() < 1e-10);
155    }
156
157    #[test]
158    fn test_euler_zero_step_size() {
159        let solution = euler_method(|x, _y| x, 0.0, 0.0, 1.0, 0.0);
160
161        assert_eq!(solution.len(), 1);
162        assert_eq!(solution[0], (0.0, 0.0));
163    }
164
165    #[test]
166    fn test_solve_euler_invalid_input() {
167        let result = solve_euler(|x, _y| x, 0.0, 0.0, 1.0, -0.1);
168        assert!(result.is_err());
169
170        let result = solve_euler(|x, _y| x, f64::NAN, 0.0, 1.0, 0.1);
171        assert!(result.is_err());
172    }
173
174    #[test]
175    fn test_euler_variable_step() {
176        let solution1 = euler_method(|x, _y| x, 0.0, 0.0, 1.0, 0.1);
177        let solution2 = euler_method(|x, _y| x, 0.0, 0.0, 1.0, 0.05);
178
179        let (_, y1) = solution1.last().unwrap();
180        let (_, y2) = solution2.last().unwrap();
181
182        assert!(solution2.len() > solution1.len());
183        assert!((y2 - 0.5).abs() < (y1 - 0.5).abs());
184    }
185}
mathhook_core/calculus/ode/numerical/euler.rs

mathhook_core/calculus/ode/numerical/
euler.rs
