differential_equations/dde/methods/
bs23.rs

1//! BS23 DDENumericalMethod for Delay Differential Equations.
2
3use crate::{
4    Error, Status,
5    alias::Evals,
6    dde::{DDE, DDENumericalMethod, methods::h_init::h_init},
7    interpolate::Interpolation,
8    traits::{CallBackData, Real, State},
9    utils::{constrain_step_size, validate_step_size_parameters},
10};
11use std::collections::VecDeque;
12
13/// Bogacki-Shampine 3(2) method adapted for DDEs.
14/// 3rd order method with embedded 2nd order error estimation and
15/// 3rd order dense output (cubic Hermite interpolation). FSAL property.
16///
17/// # Example
18///
19/// ```
20/// use differential_equations::prelude::*;
21/// use differential_equations::dde::methods::BS23;
22/// use nalgebra::{Vector2, vector};
23///
24/// let mut bs23 = BS23::new()
25///    .rtol(1e-5)
26///    .atol(1e-5);
27///
28/// let t0 = 0.0;
29/// let tf = 5.0;
30/// let y0 = vector![1.0, 0.0];
31/// let phi = |t| { // History function for t <= t0
32///     if t <= 0.0 { vector![1.0, 0.0] } else { panic!("phi called for t > t0") }
33/// };
34/// struct ExampleDDE;
35/// impl DDE<1, f64, Vector2<f64>> for ExampleDDE {
36///     fn diff(&self, t: f64, y: &Vector2<f64>, yd: &[Vector2<f64>; 1], dydt: &mut Vector2<f64>) {
37///        dydt[0] = yd[0][1];
38///        dydt[1] = -yd[0][0] - 1.0 * y[1];
39///     }
40///
41///     fn lags(&self, t: f64, y: &Vector2<f64>, lags: &mut [f64; 1]) {
42///         lags[0] = 1.0; // Constant delay tau = 1.0
43///     }
44/// }
45/// let problem = DDEProblem::new(ExampleDDE, t0, tf, y0, phi);
46/// let solution = problem.solve(&mut bs23).unwrap();
47///
48/// let (t, y) = solution.last().unwrap();
49/// println!("BS23 Solution at t={}: ({}, {})", t, y[0], y[1]);
50/// ```
51///
52/// # Settings (similar to BS23)
53/// * `rtol`, `atol`, `h0`, `h_max`, `h_min`, `max_steps`, `safe`, `fac1`, `fac2`, `beta`, `max_delay`.
54///
55/// # Default Settings (typical for BS23)
56/// * `rtol`   - 1e-3
57/// * `atol`   - 1e-6
58/// * `h0`     - None
59/// * `h_max`   - None
60/// * `h_min`   - 0.0
61/// * `max_steps` - 100_000
62/// * `safe`   - 0.9
63/// * `fac1`   - 0.2 (1/5 for order 3, can be tuned)
64/// * `fac2`   - 10.0
65/// * `beta`   - 0.0 (No PI stabilization by default for BS23, can be enabled)
66/// * `max_delay` - None
67pub struct BS23<const L: usize, T: Real, V: State<T>, H: Fn(T) -> V, D: CallBackData> {
68    pub h0: T,
69    t: T,
70    y: V,
71    h: T,
72    pub rtol: T,
73    pub atol: T,
74    pub h_max: T,
75    pub h_min: T,
76    pub max_steps: usize,
77    // BS23 is not typically for stiff problems, so n_stiff related fields are omitted for simplicity.
78    pub safe: T,
79    pub fac1: T,
80    pub fac2: T,
81    pub beta: T,
82    pub max_delay: Option<T>,
83    expo1: T,
84    facc1: T,
85    facc2: T,
86    facold: T,
87    fac11: T,
88    fac: T,
89    status: Status<T, V, D>,
90    steps: usize,
91    n_accepted: usize,
92    a: [[T; 4]; 3], // A has 3 rows for k2, k3, k4 based on k1, k2, k3
93    b: [T; 3],      // For y_new (k1, k2, k3)
94    c: [T; 3],      // For time points of k2, k3, k4
95    er: [T; 4],     // For error estimation (k1, k2, k3, k4)
96    k: [V; 4],      // k1, k2, k3, k4 (k4 is f(t+h, y_new))
97    y_old: V,
98    t_old: T,
99    h_old: T,
100    cont: [V; 4], // For cubic Hermite dense output
101    cont_buffer: VecDeque<(T, T, T, [V; 4])>,
102    phi: Option<H>,
103    t0: T,
104    tf: T,
105    posneg: T,
106    lags: [T; L], // Stores tau_i values
107    yd: [V; L],   // Stores y(t-tau_i)
108}
109
110impl<const L: usize, T: Real, V: State<T>, H: Fn(T) -> V, D: CallBackData>
111    DDENumericalMethod<L, T, V, H, D> for BS23<L, T, V, H, D>
112{
113    fn init<F>(&mut self, dde: &F, t0: T, tf: T, y0: &V, phi: H) -> Result<Evals, Error<T, V>>
114    where
115        F: DDE<L, T, V, D>,
116    {
117        let mut evals = Evals::new();
118
119        self.t = t0;
120        self.y = *y0;
121        self.t0 = t0;
122        self.tf = tf;
123        self.posneg = (tf - t0).signum();
124        self.phi = Some(phi);
125
126        if L > 0 {
127            dde.lags(self.t, &self.y, &mut self.lags);
128            for i in 0..L {
129                if self.lags[i] <= T::zero() {
130                    return Err(Error::BadInput {
131                        msg: "All lags must be positive.".to_string(),
132                    });
133                }
134                let t_delayed = self.t - self.lags[i];
135                // Assuming phi is defined for t <= t0.
136                // If t_delayed > t0 (for forward integration), it's an issue.
137                // (t_delayed - self.t0) * self.posneg > T::zero() indicates t_delayed is "beyond" t0.
138                if (t_delayed - self.t0) * self.posneg > T::default_epsilon() {
139                    // Allow for t_delayed == t0
140                    return Err(Error::BadInput {
141                        msg: format!(
142                            "Initial delayed time {} is out of history range (t <= {}).",
143                            t_delayed, self.t0
144                        ),
145                    });
146                }
147                self.yd[i] = (self.phi.as_ref().unwrap())(t_delayed);
148            }
149        }
150        dde.diff(self.t, &self.y, &self.yd, &mut self.k[0]); // k1 = f(t0, y0, yd(t0-tau))
151        evals.fcn += 1;
152
153        if self.h0 == T::zero() {
154            let h_est = h_init(
155                dde,
156                self.t,
157                self.tf,
158                &self.y,
159                3,
160                self.rtol,
161                self.atol,
162                self.h_min,
163                self.h_max,
164                self.phi.as_ref().unwrap(),
165                &self.k[0],
166                &mut evals,
167            );
168            self.h0 = h_est;
169        }
170
171        match validate_step_size_parameters::<T, V, D>(
172            self.h0, self.h_min, self.h_max, self.t, self.tf,
173        ) {
174            Ok(h0_validated) => self.h = h0_validated,
175            Err(status) => return Err(status),
176        }
177
178        self.t_old = self.t;
179        self.y_old = self.y;
180        self.h_old = self.h; // Can be zero if h0 was zero and h_init failed to produce non-zero.
181
182        self.steps = 0;
183        self.n_accepted = 0;
184        self.status = Status::Initialized;
185        Ok(evals)
186    }
187
188    fn step<F>(&mut self, dde: &F) -> Result<Evals, Error<T, V>>
189    where
190        F: DDE<L, T, V, D>,
191    {
192        let mut evals = Evals::new();
193
194        if self.steps >= self.max_steps {
195            self.status = Status::Error(Error::MaxSteps {
196                t: self.t,
197                y: self.y,
198            });
199            return Err(Error::MaxSteps {
200                t: self.t,
201                y: self.y,
202            });
203        }
204
205        let t_current_step_start = self.t;
206        let y_current_step_start = self.y;
207        let k0_current_step_start = self.k[0]; // k1 for this step, from previous step's k4 (FSAL) or init
208
209        let mut min_lag_abs = T::infinity();
210        if L > 0 {
211            let temp_y_for_lags = y_current_step_start + k0_current_step_start * self.h; // Estimate y at t+h
212            dde.lags(
213                t_current_step_start + self.h,
214                &temp_y_for_lags,
215                &mut self.lags,
216            );
217            for i in 0..L {
218                min_lag_abs = min_lag_abs.min(self.lags[i].abs());
219            }
220        }
221
222        let max_iter: usize = if L > 0 && min_lag_abs < self.h.abs() && min_lag_abs > T::zero() {
223            5
224        } else {
225            1
226        };
227
228        let mut y_new_from_iter = y_current_step_start;
229        let mut k_from_iter = self.k; // Stores k1, k2, k3, k4 of the converged iteration
230
231        let mut last_y_for_errit_calc = V::zeros();
232        let mut iteration_failed_to_converge = false;
233
234        for iter_idx in 0..max_iter {
235            if iter_idx > 0 {
236                last_y_for_errit_calc = y_new_from_iter;
237            }
238            self.k[0] = k0_current_step_start; // k1
239
240            // k2
241            let mut ti = t_current_step_start + self.c[0] * self.h; // c[0] is for k2 (0.5)
242            let mut yi = y_current_step_start + self.k[0] * (self.a[0][0] * self.h);
243            if L > 0 {
244                dde.lags(ti, &yi, &mut self.lags);
245                self.lagvals(ti, &yi);
246            }
247            dde.diff(ti, &yi, &self.yd, &mut self.k[1]); // k2 stored in self.k[1]
248
249            // k3
250            ti = t_current_step_start + self.c[1] * self.h; // c[1] is for k3 (0.75)
251            yi = y_current_step_start + self.k[1] * (self.a[1][1] * self.h); // BS23: y + h*A_32*k2
252            if L > 0 {
253                dde.lags(ti, &yi, &mut self.lags);
254                self.lagvals(ti, &yi);
255            }
256            dde.diff(ti, &yi, &self.yd, &mut self.k[2]); // k3 stored in self.k[2]
257
258            // y_new based on k1, k2, k3
259            y_new_from_iter = y_current_step_start
260                + (self.k[0] * self.b[0] + self.k[1] * self.b[1] + self.k[2] * self.b[2]) * self.h;
261
262            // k4 (FSAL property, f(t+h, y_new))
263            let t_new_val = t_current_step_start + self.h;
264            // For BS23, A[2][...] are for k4. c[2] is for k4.
265            // The k4 calculation in dde23 is f(tnew, ynew).
266            if L > 0 {
267                dde.lags(t_new_val, &y_new_from_iter, &mut self.lags);
268                self.lagvals(t_new_val, &y_new_from_iter);
269            }
270            dde.diff(t_new_val, &y_new_from_iter, &self.yd, &mut self.k[3]); // k4 stored in self.k[3]
271
272            evals.fcn += 3; // k2, k3, k4 evaluations
273            k_from_iter.copy_from_slice(&self.k);
274
275            if max_iter > 1 && iter_idx > 0 {
276                let mut errit_val = T::zero();
277                let n_dim = y_current_step_start.len();
278                for i_dim in 0..n_dim {
279                    let scale = self.atol
280                        + self.rtol
281                            * last_y_for_errit_calc
282                                .get(i_dim)
283                                .abs()
284                                .max(y_new_from_iter.get(i_dim).abs());
285                    if scale > T::zero() {
286                        let diff_val =
287                            y_new_from_iter.get(i_dim) - last_y_for_errit_calc.get(i_dim);
288                        errit_val += (diff_val / scale).powi(2);
289                    }
290                }
291                if n_dim > 0 {
292                    errit_val = (errit_val / T::from_usize(n_dim).unwrap()).sqrt();
293                }
294
295                if errit_val <= self.rtol * T::from_f64(0.1).unwrap() {
296                    break;
297                }
298            }
299            if iter_idx == max_iter - 1 && max_iter > 1 {
300                iteration_failed_to_converge = true;
301            }
302        }
303
304        if iteration_failed_to_converge {
305            self.h = (self.h.abs() * T::from_f64(0.5).unwrap()).max(self.h_min.abs()) * self.posneg;
306            if L > 0
307                && min_lag_abs > T::zero()
308                && self.h.abs() < T::from_f64(2.0).unwrap() * min_lag_abs
309            {
310                self.h = min_lag_abs * self.posneg;
311            }
312            self.h = constrain_step_size(self.h, self.h_min, self.h_max);
313            self.status = Status::RejectedStep;
314            return Ok(evals);
315        }
316
317        let mut err_final = T::zero();
318        let n = y_current_step_start.len();
319        for i in 0..n {
320            let sk = self.atol
321                + self.rtol
322                    * y_current_step_start
323                        .get(i)
324                        .abs()
325                        .max(y_new_from_iter.get(i).abs());
326            // Error = h * sum(E_i * k_i)
327            let err_comp = k_from_iter[0].get(i) * self.er[0]
328                + k_from_iter[1].get(i) * self.er[1]
329                + k_from_iter[2].get(i) * self.er[2]
330                + k_from_iter[3].get(i) * self.er[3];
331            let erri = self.h * err_comp;
332            if sk > T::zero() {
333                err_final += (erri / sk).powi(2);
334            }
335        }
336        if n > 0 {
337            err_final = (err_final / T::from_usize(n).unwrap()).sqrt();
338        }
339
340        self.fac11 = err_final.powf(self.expo1);
341        let fac_beta = if self.beta > T::zero() && self.facold > T::zero() {
342            self.facold.powf(self.beta)
343        } else {
344            T::one()
345        };
346        self.fac = self.fac11 / fac_beta;
347        self.fac = self.facc2.max(self.facc1.min(self.fac / self.safe));
348        let mut h_new_final = self.h / self.fac;
349
350        let t_new_val = t_current_step_start + self.h;
351
352        if err_final <= T::one() {
353            self.facold = err_final.max(T::from_f64(1.0e-4).unwrap());
354            self.n_accepted += 1;
355
356            // Dense output coefficients for cubic Hermite
357            // y(s) = c0 + s(c1 + s(c2 + s*c3)) where s = (t_int - t_old)/h_old
358            // c0 = y_old
359            // c1 = h_old * k_old
360            // c2 = 3*(y_new - y_old) - h_old*(2*k_old + k_new)
361            // c3 = -2*(y_new - y_old) + h_old*(k_old + k_new)
362            let k_old_for_cont = k_from_iter[0]; // k1 at t_current_step_start
363            let k_new_for_cont = k_from_iter[3]; // k4 at t_new_val
364            let y_diff_cont = y_new_from_iter - y_current_step_start;
365
366            self.cont[0] = y_current_step_start;
367            self.cont[1] = k_old_for_cont * self.h;
368            self.cont[2] = y_diff_cont * T::from_f64(3.0).unwrap()
369                - (k_old_for_cont * T::from_f64(2.0).unwrap() + k_new_for_cont) * self.h;
370            self.cont[3] = y_diff_cont * T::from_f64(-2.0).unwrap()
371                + (k_old_for_cont + k_new_for_cont) * self.h;
372
373            self.cont_buffer
374                .push_back((t_current_step_start, t_new_val, self.h, self.cont));
375
376            if let Some(max_delay_val) = self.max_delay {
377                let prune_time = if self.posneg > T::zero() {
378                    t_new_val - max_delay_val
379                } else {
380                    t_new_val + max_delay_val
381                };
382                while let Some((buf_t_start, buf_t_end, _, _)) = self.cont_buffer.front() {
383                    if (self.posneg > T::zero() && *buf_t_end < prune_time)
384                        || (self.posneg < T::zero() && *buf_t_start > prune_time)
385                    {
386                        self.cont_buffer.pop_front();
387                    } else {
388                        break;
389                    }
390                }
391            }
392
393            self.y_old = y_current_step_start;
394            self.t_old = t_current_step_start;
395            self.h_old = self.h;
396
397            self.k[0] = k_from_iter[3]; // FSAL: k4 becomes k1 for next step
398            self.y = y_new_from_iter;
399            self.t = t_new_val;
400
401            if let Status::RejectedStep = self.status {
402                h_new_final = self.h_old.min(h_new_final);
403                self.status = Status::Solving;
404            }
405        } else {
406            h_new_final = self.h / self.facc1.min(self.fac11 / self.safe);
407            self.status = Status::RejectedStep;
408        }
409
410        self.steps += 1;
411        self.h = constrain_step_size(h_new_final, self.h_min, self.h_max);
412        Ok(evals)
413    }
414
415    fn t(&self) -> T {
416        self.t
417    }
418    fn y(&self) -> &V {
419        &self.y
420    }
421    fn t_prev(&self) -> T {
422        self.t_old
423    }
424    fn y_prev(&self) -> &V {
425        &self.y_old
426    }
427    fn h(&self) -> T {
428        self.h
429    }
430    fn set_h(&mut self, h: T) {
431        self.h = h;
432    }
433    fn status(&self) -> &Status<T, V, D> {
434        &self.status
435    }
436    fn set_status(&mut self, status: Status<T, V, D>) {
437        self.status = status;
438    }
439}
440
441impl<const L: usize, T: Real, V: State<T>, H: Fn(T) -> V, D: CallBackData> Interpolation<T, V>
442    for BS23<L, T, V, H, D>
443{
444    fn interpolate(&mut self, t_interp: T) -> Result<V, Error<T, V>> {
445        // This interpolates within the last successfully completed step [t_old, t]
446        if (t_interp - self.t_old) * self.posneg < T::zero()
447            || (t_interp - self.t) * self.posneg > T::zero()
448        {
449            // If t_interp is exactly t_old or t, still proceed.
450            if (t_interp - self.t_old).abs() > T::default_epsilon()
451                && (t_interp - self.t).abs() > T::default_epsilon()
452            {
453                return Err(Error::OutOfBounds {
454                    t_interp,
455                    t_prev: self.t_old,
456                    t_curr: self.t,
457                });
458            }
459        }
460
461        let s = if self.h_old == T::zero() {
462            if (t_interp - self.t_old).abs() < T::default_epsilon() {
463                T::zero()
464            } else {
465                T::one()
466            }
467        } else {
468            (t_interp - self.t_old) / self.h_old
469        };
470
471        // y_interp = cont[0] + s * (cont[1] + s * (cont[2] + s * cont[3]))
472        let y_interp = self.cont[0] + (self.cont[1] + (self.cont[2] + self.cont[3] * s) * s) * s;
473        Ok(y_interp)
474    }
475}
476
477impl<const L: usize, T: Real, V: State<T>, H: Fn(T) -> V, D: CallBackData> BS23<L, T, V, H, D> {
478    pub fn new() -> Self {
479        Self::default()
480    }
481
482    // Builder methods
483    pub fn rtol(mut self, rtol: T) -> Self {
484        self.rtol = rtol;
485        self
486    }
487    pub fn atol(mut self, atol: T) -> Self {
488        self.atol = atol;
489        self
490    }
491    pub fn h0(mut self, h0: T) -> Self {
492        self.h0 = h0;
493        self
494    }
495    pub fn h_max(mut self, h_max: T) -> Self {
496        self.h_max = h_max;
497        self
498    }
499    pub fn h_min(mut self, h_min: T) -> Self {
500        self.h_min = h_min;
501        self
502    }
503    pub fn max_steps(mut self, max_steps: usize) -> Self {
504        self.max_steps = max_steps;
505        self
506    }
507    pub fn safe(mut self, safe: T) -> Self {
508        self.safe = safe;
509        self
510    }
511    pub fn fac1(mut self, fac1: T) -> Self {
512        self.fac1 = fac1;
513        self.facc1 = T::one() / fac1;
514        self
515    }
516    pub fn fac2(mut self, fac2: T) -> Self {
517        self.fac2 = fac2;
518        self.facc2 = T::one() / fac2;
519        self
520    }
521    pub fn beta(mut self, beta: T) -> Self {
522        self.beta = beta;
523        self
524    }
525    pub fn max_delay(mut self, max_delay: T) -> Self {
526        self.max_delay = Some(max_delay.abs());
527        self
528    }
529
530    fn lagvals(&mut self, t_stage: T, _y_stage: &V) {
531        // y_stage not used if interpolating from buffer
532        for i in 0..L {
533            let t_delayed = t_stage - self.lags[i];
534            if (t_delayed - self.t0) * self.posneg <= T::default_epsilon() {
535                // t_delayed <= t0 (forward) or t_delayed >= t0 (backward)
536                self.yd[i] = (self.phi.as_ref().unwrap())(t_delayed);
537            } else {
538                // Search cont_buffer (most recent first)
539                let mut found_in_buffer = false;
540                for (buf_t_start, buf_t_end, buf_h, buf_cont) in self.cont_buffer.iter().rev() {
541                    // Check if t_delayed is within [buf_t_start, buf_t_end]
542                    if (t_delayed - *buf_t_start) * self.posneg >= -T::default_epsilon()
543                        && (t_delayed - *buf_t_end) * self.posneg <= T::default_epsilon()
544                    {
545                        let s = if *buf_h == T::zero() {
546                            if (t_delayed - *buf_t_start).abs() < T::default_epsilon() {
547                                T::zero()
548                            } else {
549                                T::one()
550                            }
551                        } else {
552                            (t_delayed - *buf_t_start) / *buf_h
553                        };
554                        self.yd[i] =
555                            buf_cont[0] + (buf_cont[1] + (buf_cont[2] + buf_cont[3] * s) * s) * s;
556                        found_in_buffer = true;
557                        break;
558                    }
559                }
560                if !found_in_buffer {
561                    // Extrapolate using the most recent interval if t_delayed is beyond it,
562                    // or if buffer is empty (should ideally not happen if t_delayed > t0).
563                    // This part might need more robust handling (e.g., error or specific extrapolation strategy)
564                    if let Some((buf_t_start, _buf_t_end, buf_h, buf_cont)) =
565                        self.cont_buffer.back()
566                    {
567                        let s = if *buf_h == T::zero() {
568                            T::one()
569                        } else {
570                            (t_delayed - *buf_t_start) / *buf_h
571                        }; // Extrapolation
572                        self.yd[i] =
573                            buf_cont[0] + (buf_cont[1] + (buf_cont[2] + buf_cont[3] * s) * s) * s;
574                    } else {
575                        // Fallback: if cont_buffer is empty and t_delayed > t0, this is an issue.
576                        // This implies we are trying to access a future point not yet computed or history not covering it.
577                        // For now, use phi, though this might be incorrect if t_delayed > t0.
578                        // A proper error or specific handling for "future" lookups in early steps might be needed.
579                        self.yd[i] = (self.phi.as_ref().unwrap())(t_delayed);
580                        // Consider panic or error: panic!("Lag value lookup failed for t_delayed={}", t_delayed);
581                    }
582                }
583            }
584        }
585    }
586}
587
588// Bogacki-Shampine Coefficients (FSAL variant)
589// k1 = f(t,y)
590// k2 = f(t+0.5h, y + 0.5h*k1)
591// k3 = f(t+0.75h, y + 0.75h*k2)
592// y_new = y + h*(2/9 k1 + 1/3 k2 + 4/9 k3)
593// k4 = f(t+h, y_new)
594// Error: y_new - y_hat where y_hat = y + h*(7/24 k1 + 1/4 k2 + 1/3 k3 + 1/8 k4)
595// E = b_sol - b_err_coeffs = [2/9-7/24, 1/3-1/4, 4/9-1/3, 0-1/8]
596//   = [-5/72, 1/12, 1/9, -1/8]
597
598const BS23_C: [f64; 3] = [1.0 / 2.0, 3.0 / 4.0, 1.0]; // Time points for k2, k3, k4
599const BS23_A: [[f64; 4]; 3] = [
600    // A_ij * h * k_j
601    [1.0 / 2.0, 0.0, 0.0, 0.0],             // k2 depends on k1
602    [0.0, 3.0 / 4.0, 0.0, 0.0],             // k3 depends on k2
603    [2.0 / 9.0, 1.0 / 3.0, 4.0 / 9.0, 0.0], // k4 (f(t+h, y_new)) uses y_new which depends on k1,k2,k3. This row is not directly used for y_i stages.
604];
605const BS23_B_SOL: [f64; 3] = [2.0 / 9.0, 1.0 / 3.0, 4.0 / 9.0]; // For solution (uses k1, k2, k3)
606const BS23_E: [f64; 4] = [-5.0 / 72.0, 1.0 / 12.0, 1.0 / 9.0, -1.0 / 8.0]; // For error estimate (uses k1, k2, k3, k4)
607
608impl<const L: usize, T: Real, V: State<T>, H: Fn(T) -> V, D: CallBackData> Default
609    for BS23<L, T, V, H, D>
610{
611    fn default() -> Self {
612        let c_conv = BS23_C.map(|x| T::from_f64(x).unwrap());
613        let a_conv = BS23_A.map(|row| row.map(|x| T::from_f64(x).unwrap()));
614        let b_sol_conv = BS23_B_SOL.map(|x| T::from_f64(x).unwrap());
615        let er_conv = BS23_E.map(|x| T::from_f64(x).unwrap());
616
617        let expo1_final = T::one() / T::from_f64(3.0).unwrap();
618
619        let fac1_default = T::from_f64(0.2).unwrap(); // Can be tuned
620        let fac2_default = T::from_f64(10.0).unwrap();
621
622        BS23 {
623            t: T::zero(),
624            y: V::zeros(),
625            h: T::zero(),
626            h0: T::zero(),
627            rtol: T::from_f64(1e-3).unwrap(),
628            atol: T::from_f64(1e-6).unwrap(),
629            h_max: T::infinity(),
630            h_min: T::zero(),
631            max_steps: 100_000,
632            safe: T::from_f64(0.9).unwrap(),
633            fac1: fac1_default,
634            fac2: fac2_default,
635            beta: T::zero(), // No PI by default for BS23
636            max_delay: None,
637            expo1: expo1_final,
638            facc1: T::one() / fac1_default,
639            facc2: T::one() / fac2_default,
640            facold: T::from_f64(1.0e-4).unwrap(),
641            fac11: T::zero(),
642            fac: T::zero(),
643            status: Status::Uninitialized,
644            steps: 0,
645            n_accepted: 0,
646            a: a_conv,
647            b: b_sol_conv,
648            c: c_conv,
649            er: er_conv,
650            k: [V::zeros(); 4],
651            y_old: V::zeros(),
652            t_old: T::zero(),
653            h_old: T::zero(),
654            cont: [V::zeros(); 4],
655            cont_buffer: VecDeque::new(),
656            phi: None,
657            t0: T::zero(),
658            tf: T::zero(),
659            posneg: T::zero(),
660            lags: [T::zero(); L],
661            yd: [V::zeros(); L],
662        }
663    }
664}
differential_equations/dde/methods/bs23.rs

differential_equations/dde/methods/
bs23.rs
