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

differential_equations/methods/erk/
mod.rs
