differential_equations/methods/erk/
mod.rs

1//! Explicit Runge-Kutta (ERK) methods
2
3mod adaptive;
4mod dormandprince;
5mod fixed;
6
7use std::{collections::VecDeque, marker::PhantomData};
8
9use crate::{
10    methods::Delay,
11    status::Status,
12    traits::{CallBackData, Real, State},
13};
14
15/// Runge-Kutta solver that can handle:
16/// - Fixed-step methods with cubic Hermite interpolation
17/// - Adaptive step methods with embedded error estimation and cubic Hermite interpolation
18/// - Advanced methods with dense output interpolation using Butcher tableau coefficients
19///
20/// # Type Parameters
21///
22/// * `E`: Equation type (e.g., Ordinary, Delay, Stochastic)
23/// * `F`: Family type (e.g., Adaptive, Fixed, DormandPrince, etc.)
24/// * `T`: Real number type (f32, f64)
25/// * `Y`: State vector type
26/// * `D`: Callback data type
27/// * `const O`: Order of the method
28/// * `const S`: Number of stages in the method
29/// * `const I`: Total number of stages including interpolation (equal to S for methods without dense output)
30pub struct ExplicitRungeKutta<
31    E,
32    F,
33    T: Real,
34    Y: State<T>,
35    D: CallBackData,
36    const O: usize,
37    const S: usize,
38    const I: usize,
39> {
40    // Domain of problem
41    t0: T,
42
43    // Initial Step Size
44    pub h0: T,
45
46    // Current Step Size
47    h: T,
48
49    // Current State
50    t: T,
51    y: Y,
52    dydt: Y,
53
54    // Previous State
55    h_prev: T,
56    t_prev: T,
57    y_prev: Y,
58    dydt_prev: Y,
59
60    // Stage values
61    k: [Y; I],
62
63    // Constants from Butcher tableau
64    c: [T; I],
65    a: [[T; I]; I],
66    b: [T; S],
67    bh: Option<[T; S]>, // Lower order coefficients for embedded methods
68    er: Option<[T; S]>, // Error estimation coefficients
69
70    // Interpolation coefficients
71    bi: Option<[[T; I]; I]>, // Optional for methods without dense output
72    cont: [Y; O],
73
74    // Settings
75    pub rtol: T,
76    pub atol: T,
77    pub h_max: T,
78    pub h_min: T,
79    pub max_steps: usize,
80    pub max_rejects: usize,
81    pub safety_factor: T,
82    pub min_scale: T,
83    pub max_scale: T,
84
85    // Iteration tracking
86    stiffness_counter: usize,
87    non_stiffness_counter: usize,
88    steps: usize,
89
90    // Status
91    status: Status<T, Y, D>,
92
93    // Method info
94    order: usize,
95    stages: usize,
96    dense_stages: usize,
97    fsal: bool, // First Same As Last (FSAL) property
98
99    // Family classification
100    family: PhantomData<F>,
101
102    // Equation type
103    equation: PhantomData<E>,
104
105    // DDE (Delay) support
106    history: VecDeque<(T, Y, Y)>, // (t, y, dydt)
107    max_delay: Option<T>,         // Minimum delay for DDEs so that buffer can be emptied
108}
109
110impl<E, F, T: Real, Y: State<T>, D: CallBackData, const O: usize, const S: usize, const I: usize>
111    Default for ExplicitRungeKutta<E, F, T, Y, D, O, S, I>
112{
113    fn default() -> Self {
114        Self {
115            t0: T::zero(),
116            h0: T::zero(),
117            h: T::zero(),
118            t: T::zero(),
119            y: Y::zeros(),
120            dydt: Y::zeros(),
121            h_prev: T::zero(),
122            t_prev: T::zero(),
123            y_prev: Y::zeros(),
124            dydt_prev: Y::zeros(),
125            k: [Y::zeros(); I],
126            c: [T::zero(); I],
127            a: [[T::zero(); I]; I],
128            b: [T::zero(); S],
129            bh: None,
130            er: None,
131            bi: None,
132            cont: [Y::zeros(); O],
133            rtol: T::from_f64(1.0e-6).unwrap(),
134            atol: T::from_f64(1.0e-6).unwrap(),
135            h_max: T::infinity(),
136            h_min: T::zero(),
137            max_steps: 10_000,
138            max_rejects: 100,
139            safety_factor: T::from_f64(0.9).unwrap(),
140            min_scale: T::from_f64(0.2).unwrap(),
141            max_scale: T::from_f64(10.0).unwrap(),
142            stiffness_counter: 0,
143            non_stiffness_counter: 0,
144            steps: 0,
145            status: Status::Uninitialized,
146            order: O,
147            stages: S,
148            dense_stages: I,
149            fsal: false,
150            family: PhantomData,
151            equation: PhantomData,
152            history: VecDeque::new(),
153            max_delay: None,
154        }
155    }
156}
157
158impl<E, F, T: Real, Y: State<T>, D: CallBackData, const O: usize, const S: usize, const I: usize>
159    ExplicitRungeKutta<E, F, T, Y, D, O, S, I>
160{
161    /// Set the relative tolerance for error control
162    pub fn rtol(mut self, rtol: T) -> Self {
163        self.rtol = rtol;
164        self
165    }
166
167    /// Set the absolute tolerance for error control
168    pub fn atol(mut self, atol: T) -> Self {
169        self.atol = atol;
170        self
171    }
172
173    /// Set the initial step size
174    pub fn h0(mut self, h0: T) -> Self {
175        self.h0 = h0;
176        self
177    }
178
179    /// Set the minimum allowed step size
180    pub fn h_min(mut self, h_min: T) -> Self {
181        self.h_min = h_min;
182        self
183    }
184
185    /// Set the maximum allowed step size
186    pub fn h_max(mut self, h_max: T) -> Self {
187        self.h_max = h_max;
188        self
189    }
190
191    /// Set the maximum number of steps allowed
192    pub fn max_steps(mut self, max_steps: usize) -> Self {
193        self.max_steps = max_steps;
194        self
195    }
196
197    /// Set the maximum number of consecutive rejected steps before declaring stiffness
198    pub fn max_rejects(mut self, max_rejects: usize) -> Self {
199        self.max_rejects = max_rejects;
200        self
201    }
202
203    /// Set the safety factor for step size control (default: 0.9)
204    pub fn safety_factor(mut self, safety_factor: T) -> Self {
205        self.safety_factor = safety_factor;
206        self
207    }
208
209    /// Set the minimum scale factor for step size changes (default: 0.2)
210    pub fn min_scale(mut self, min_scale: T) -> Self {
211        self.min_scale = min_scale;
212        self
213    }
214
215    /// Set the maximum scale factor for step size changes (default: 10.0)
216    pub fn max_scale(mut self, max_scale: T) -> Self {
217        self.max_scale = max_scale;
218        self
219    }
220
221    /// Get the order of the method
222    pub fn order(&self) -> usize {
223        self.order
224    }
225
226    /// Get the number of stages in the method
227    pub fn stages(&self) -> usize {
228        self.stages
229    }
230
231    /// Get the number of terms in the dense output interpolation polynomial
232    pub fn dense_stages(&self) -> usize {
233        self.dense_stages
234    }
235}
236
237impl<F, T: Real, Y: State<T>, D: CallBackData, const O: usize, const S: usize, const I: usize>
238    ExplicitRungeKutta<Delay, F, T, Y, D, O, S, I>
239{
240    /// Set the maximum delay for DDEs
241    pub fn max_delay(mut self, max_delay: T) -> Self {
242        self.max_delay = Some(max_delay);
243        self
244    }
245}
differential_equations/methods/erk/mod.rs

differential_equations/methods/erk/
mod.rs
