differential_equations/methods/erk/
mod.rs

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

differential_equations/methods/erk/
mod.rs
