proof_engine/solver/
stability.rs

1//! Stability analysis — automatic step size selection and stiff equation detection.
2
3use super::ode::{OdeSystem, OdeState, OdeSolver, OdeMethod};
4
5/// Results of a stability analysis.
6#[derive(Debug, Clone)]
7pub struct StabilityAnalysis {
8    pub is_stiff: bool,
9    pub recommended_method: OdeMethod,
10    pub recommended_dt: f64,
11    pub max_eigenvalue_estimate: f64,
12    pub stability_ratio: f64,
13}
14
15impl StabilityAnalysis {
16    /// Analyze an ODE system at the given state to determine stiffness and recommend solver settings.
17    pub fn analyze(system: &dyn OdeSystem, state: &OdeState, dt: f64) -> Self {
18        let n = system.dimension();
19        let h = 1e-6;
20
21        // Estimate Jacobian eigenvalues via power iteration on finite differences
22        let mut max_eig = 0.0f64;
23        let mut min_eig = f64::MAX;
24
25        let mut f0 = vec![0.0; n];
26        system.evaluate(state.t, &state.y, &mut f0);
27
28        for j in 0..n {
29            let mut y_perturbed = state.y.clone();
30            y_perturbed[j] += h;
31            let mut f1 = vec![0.0; n];
32            system.evaluate(state.t, &y_perturbed, &mut f1);
33
34            // Approximate column j of the Jacobian
35            let col_norm: f64 = (0..n).map(|i| ((f1[i] - f0[i]) / h).powi(2)).sum::<f64>().sqrt();
36            if col_norm > max_eig { max_eig = col_norm; }
37            if col_norm < min_eig { min_eig = col_norm; }
38        }
39
40        let stiffness_ratio = if min_eig > 1e-15 { max_eig / min_eig } else { 1.0 };
41        let is_stiff = stiffness_ratio > 1000.0;
42
43        let stability_limit = if max_eig > 1e-15 { 2.8 / max_eig } else { dt };
44        let recommended_dt = (stability_limit * 0.8).min(dt);
45
46        let recommended_method = if is_stiff {
47            OdeMethod::ImplicitEuler
48        } else if stiffness_ratio > 100.0 {
49            OdeMethod::RungeKutta45
50        } else {
51            OdeMethod::RungeKutta4
52        };
53
54        Self {
55            is_stiff,
56            recommended_method,
57            recommended_dt,
58            max_eigenvalue_estimate: max_eig,
59            stability_ratio: stiffness_ratio,
60        }
61    }
62}
63
64#[cfg(test)]
65mod tests {
66    use super::*;
67    use crate::solver::ode::HarmonicOscillator;
68
69    #[test]
70    fn harmonic_oscillator_not_stiff() {
71        let sys = HarmonicOscillator { omega: 1.0 };
72        let state = OdeState { t: 0.0, y: vec![1.0, 0.0] };
73        let analysis = StabilityAnalysis::analyze(&sys, &state, 0.01);
74        assert!(!analysis.is_stiff);
75        assert!(analysis.recommended_dt > 0.0);
76    }
77}
proof_engine/solver/stability.rs

proof_engine/solver/
stability.rs
