differential_equations/methods/bdf/
mod.rs

1//! Adaptive variable-order Backward Differentiation Formula solver.
2//!
3//! Implements Backward Differentiation Formula orders 1 through 5 with automatic order selection and
4//! adaptive step size control. The solver starts at order 1 and adjusts
5//! both order and step size to balance accuracy and efficiency.
6//!
7//! Based on the approach described in:
8//! - Hairer, E., & Wanner, G. (1996). Solving Ordinary Differential
9//!   Equations II, Section IV.2 (Backward Differentiation Formula methods).
10
11mod adaptive;
12
13use std::marker::PhantomData;
14
15use crate::{
16    linalg::Matrix,
17    status::Status,
18    tolerance::Tolerance,
19    traits::{Real, State},
20};
21
22/// Maximum Backward Differentiation Formula order (order 5 is the highest practically useful order).
23pub const MAX_ORDER: usize = 5;
24const BACKWARD_DIFFERENTIATION_FORMULA_ROWS: usize = MAX_ORDER + 3;
25
26/// Adaptive variable-order Backward Differentiation Formula solver.
27///
28/// Uses Backward Differentiation Formula orders 1-5 with automatic order and step size selection.
29/// Newton iteration resolves the implicit Backward Differentiation Formula equations. Dense output
30/// via polynomial interpolation is supported.
31///
32/// The solver starts at order 1 and considers increasing the order after
33/// enough successful steps. On step rejection, the order is decreased.
34pub struct BackwardDifferentiationFormula<E, T: Real, Y: State<T>> {
35    pub h0: T,
36    pub rtol: Tolerance<T>,
37    pub atol: Tolerance<T>,
38    pub h_max: T,
39    pub h_min: T,
40    pub max_steps: usize,
41    pub safety_factor: T,
42    pub min_scale: T,
43    pub max_scale: T,
44    pub filter: fn(T) -> T,
45    pub newton_tol: T,
46    pub max_newton_iter: usize,
47    pub max_order: usize,
48
49    h: T,
50    t: T,
51    y: Y,
52    dydt: Y,
53
54    t_prev: T,
55    y_prev: Y,
56    h_prev: T,
57    d: [Y; BACKWARD_DIFFERENTIATION_FORMULA_ROWS],
58
59    order: usize,
60    n_equal_steps: usize,
61
62    gamma: [T; MAX_ORDER + 1],
63    alpha: [T; MAX_ORDER + 1],
64    error_const: [T; MAX_ORDER + 2],
65
66    jacobian: Matrix<T>,
67    newton_matrix: Matrix<T>,
68    ip: Vec<usize>,
69    lu_valid: bool,
70
71    status: Status<T, Y>,
72    steps: usize,
73
74    _phantom: PhantomData<E>,
75}
76
77impl<E, T: Real, Y: State<T>> Default for BackwardDifferentiationFormula<E, T, Y> {
78    fn default() -> Self {
79        Self {
80            h0: T::zero(),
81            rtol: Tolerance::Scalar(T::from_f64(1e-8).unwrap()),
82            atol: Tolerance::Scalar(T::from_f64(1e-10).unwrap()),
83            h_max: T::infinity(),
84            h_min: T::zero(),
85            max_steps: 10_000,
86            safety_factor: T::from_f64(0.9).unwrap(),
87            min_scale: T::from_f64(0.2).unwrap(),
88            max_scale: T::from_f64(10.0).unwrap(),
89            filter: |h| h,
90            newton_tol: T::zero(),
91            max_newton_iter: 4,
92            max_order: MAX_ORDER,
93
94            h: T::zero(),
95            t: T::zero(),
96            y: Y::zeros(),
97            dydt: Y::zeros(),
98
99            t_prev: T::zero(),
100            y_prev: Y::zeros(),
101            h_prev: T::zero(),
102            d: core::array::from_fn(|_| Y::zeros()),
103
104            order: 1,
105            n_equal_steps: 0,
106
107            gamma: [T::zero(); MAX_ORDER + 1],
108            alpha: [T::zero(); MAX_ORDER + 1],
109            error_const: [T::zero(); MAX_ORDER + 2],
110
111            jacobian: Matrix::zeros(0, 0),
112            newton_matrix: Matrix::zeros(0, 0),
113            ip: Vec::new(),
114            lu_valid: false,
115
116            status: Status::Uninitialized,
117            steps: 0,
118
119            _phantom: PhantomData,
120        }
121    }
122}
123
124impl<E, T: Real, Y: State<T>> BackwardDifferentiationFormula<E, T, Y> {
125    pub fn adaptive() -> Self {
126        Self::default()
127    }
128
129    pub fn rtol<V: Into<Tolerance<T>>>(mut self, rtol: V) -> Self {
130        self.rtol = rtol.into();
131        self
132    }
133
134    pub fn atol<V: Into<Tolerance<T>>>(mut self, atol: V) -> Self {
135        self.atol = atol.into();
136        self
137    }
138
139    pub fn h0(mut self, h0: T) -> Self {
140        self.h0 = h0;
141        self
142    }
143
144    pub fn h_min(mut self, h_min: T) -> Self {
145        self.h_min = h_min;
146        self
147    }
148
149    pub fn h_max(mut self, h_max: T) -> Self {
150        self.h_max = h_max;
151        self
152    }
153
154    pub fn max_steps(mut self, max_steps: usize) -> Self {
155        self.max_steps = max_steps;
156        self
157    }
158
159    pub fn safety_factor(mut self, safety_factor: T) -> Self {
160        self.safety_factor = safety_factor;
161        self
162    }
163
164    pub fn newton_tol(mut self, newton_tol: T) -> Self {
165        self.newton_tol = newton_tol;
166        self
167    }
168
169    pub fn max_newton_iter(mut self, max_newton_iter: usize) -> Self {
170        self.max_newton_iter = max_newton_iter;
171        self
172    }
173
174    pub fn max_order(mut self, max_order: usize) -> Self {
175        self.max_order = max_order.clamp(1, MAX_ORDER);
176        self
177    }
178
179    pub fn filter(mut self, filter: fn(T) -> T) -> Self {
180        self.filter = filter;
181        self
182    }
183
184    pub fn order(&self) -> usize {
185        self.order
186    }
187}
differential_equations/methods/bdf/mod.rs

differential_equations/methods/bdf/
mod.rs
