scirs2_integrate/ode/methods/
adaptive.rs

1//! Adaptive ODE solver methods
2//!
3//! This module implements adaptive methods for solving ODEs,
4//! including Dormand-Prince (RK45), Bogacki-Shampine (RK23),
5//! and Dormand-Prince 8th order (DOP853) methods.
6
7use crate::common::IntegrateFloat;
8use crate::error::IntegrateResult;
9use crate::ode::types::{ODEMethod, ODEOptions, ODEResult};
10use scirs2_core::ndarray::{Array1, ArrayView1};
11
12/// Solve ODE using the Dormand-Prince method (RK45)
13///
14/// This is an adaptive step size method based on embedded Runge-Kutta formulas.
15/// It uses a 5th-order method with a 4th-order error estimate.
16///
17/// # Arguments
18///
19/// * `f` - ODE function dy/dt = f(t, y)
20/// * `t_span` - Time span [t_start, t_end]
21/// * `y0` - Initial condition
22/// * `opts` - Solver options
23///
24/// # Returns
25///
26/// The solution as an ODEResult or an error
27#[allow(dead_code)]
28pub fn rk45_method<F, Func>(
29    f: Func,
30    t_span: [F; 2],
31    y0: Array1<F>,
32    opts: ODEOptions<F>,
33) -> IntegrateResult<ODEResult<F>>
34where
35    F: IntegrateFloat,
36    Func: Fn(F, ArrayView1<F>) -> Array1<F>,
37{
38    // Initialize
39    let [t_start, t_end] = t_span;
40    let n_dim = y0.len();
41
42    // Determine initial step size if not provided
43    let h0 = opts.h0.unwrap_or_else(|| {
44        // Simple heuristic for initial step size
45        let _span = t_end - t_start;
46        _span / F::from_usize(100).expect("Operation failed")
47    });
48
49    // Determine minimum and maximum step sizes
50    let min_step = opts.min_step.unwrap_or_else(|| {
51        let _span = t_end - t_start;
52        _span * F::from_f64(1e-8).expect("Operation failed") // Minimal step size
53    });
54
55    let max_step = opts.max_step.unwrap_or_else(|| {
56        t_end - t_start // Maximum step can be the whole interval
57    });
58
59    // Current state
60    let mut t = t_start;
61    let mut y = y0.clone();
62    let mut h = h0;
63
64    // Storage for results
65    let mut t_values = vec![t_start];
66    let mut y_values = vec![y0.clone()];
67
68    // Statistics
69    let mut func_evals = 0;
70    let mut step_count = 0;
71    let mut accepted_steps = 0;
72    let mut rejected_steps = 0;
73
74    // Dormand-Prince coefficients
75    // Time steps
76    let c2 = F::from_f64(1.0 / 5.0).expect("Operation failed");
77    let c3 = F::from_f64(3.0 / 10.0).expect("Operation failed");
78    let c4 = F::from_f64(4.0 / 5.0).expect("Operation failed");
79    let c5 = F::from_f64(8.0 / 9.0).expect("Operation failed");
80    let c6 = F::one();
81
82    // Main integration loop
83    while t < t_end && step_count < opts.max_steps {
84        // Adjust step size for the last step if needed
85        if t + h > t_end {
86            h = t_end - t;
87        }
88
89        // Limit step size to bounds
90        h = h.min(max_step).max(min_step);
91
92        // Runge-Kutta stages
93        let k1 = f(t, y.view());
94
95        // Manually compute the stages to avoid type mismatches
96        let mut y_stage = y.clone();
97        for i in 0..n_dim {
98            y_stage[i] = y[i] + h * F::from_f64(1.0 / 5.0).expect("Operation failed") * k1[i];
99        }
100        let k2 = f(t + c2 * h, y_stage.view());
101
102        let mut y_stage = y.clone();
103        for i in 0..n_dim {
104            y_stage[i] = y[i]
105                + h * (F::from_f64(3.0 / 40.0).expect("Operation failed") * k1[i]
106                    + F::from_f64(9.0 / 40.0).expect("Operation failed") * k2[i]);
107        }
108        let k3 = f(t + c3 * h, y_stage.view());
109
110        let mut y_stage = y.clone();
111        for i in 0..n_dim {
112            y_stage[i] = y[i]
113                + h * (F::from_f64(44.0 / 45.0).expect("Operation failed") * k1[i]
114                    + F::from_f64(-56.0 / 15.0).expect("Operation failed") * k2[i]
115                    + F::from_f64(32.0 / 9.0).expect("Operation failed") * k3[i]);
116        }
117        let k4 = f(t + c4 * h, y_stage.view());
118
119        let mut y_stage = y.clone();
120        for i in 0..n_dim {
121            y_stage[i] = y[i]
122                + h * (F::from_f64(19372.0 / 6561.0).expect("Operation failed") * k1[i]
123                    + F::from_f64(-25360.0 / 2187.0).expect("Operation failed") * k2[i]
124                    + F::from_f64(64448.0 / 6561.0).expect("Operation failed") * k3[i]
125                    + F::from_f64(-212.0 / 729.0).expect("Operation failed") * k4[i]);
126        }
127        let k5 = f(t + c5 * h, y_stage.view());
128
129        let mut y_stage = y.clone();
130        for i in 0..n_dim {
131            y_stage[i] = y[i]
132                + h * (F::from_f64(9017.0 / 3168.0).expect("Operation failed") * k1[i]
133                    + F::from_f64(-355.0 / 33.0).expect("Operation failed") * k2[i]
134                    + F::from_f64(46732.0 / 5247.0).expect("Operation failed") * k3[i]
135                    + F::from_f64(49.0 / 176.0).expect("Operation failed") * k4[i]
136                    + F::from_f64(-5103.0 / 18656.0).expect("Operation failed") * k5[i]);
137        }
138        let k6 = f(t + c6 * h, y_stage.view());
139
140        let mut y_stage = y.clone();
141        for i in 0..n_dim {
142            y_stage[i] = y[i]
143                + h * (F::from_f64(35.0 / 384.0).expect("Operation failed") * k1[i]
144                    + F::zero() * k2[i]
145                    + F::from_f64(500.0 / 1113.0).expect("Operation failed") * k3[i]
146                    + F::from_f64(125.0 / 192.0).expect("Operation failed") * k4[i]
147                    + F::from_f64(-2187.0 / 6784.0).expect("Operation failed") * k5[i]
148                    + F::from_f64(11.0 / 84.0).expect("Operation failed") * k6[i]);
149        }
150        let k7 = f(t + h, y_stage.view());
151
152        func_evals += 7;
153
154        // 5th order solution
155        let mut y5 = y.clone();
156        for i in 0..n_dim {
157            y5[i] = y[i]
158                + h * (F::from_f64(35.0 / 384.0).expect("Operation failed") * k1[i]
159                    + F::zero() * k2[i]
160                    + F::from_f64(500.0 / 1113.0).expect("Operation failed") * k3[i]
161                    + F::from_f64(125.0 / 192.0).expect("Operation failed") * k4[i]
162                    + F::from_f64(-2187.0 / 6784.0).expect("Operation failed") * k5[i]
163                    + F::from_f64(11.0 / 84.0).expect("Operation failed") * k6[i]
164                    + F::zero() * k7[i]);
165        }
166
167        // 4th order solution
168        let mut y4 = y.clone();
169        for i in 0..n_dim {
170            y4[i] = y[i]
171                + h * (F::from_f64(5179.0 / 57600.0).expect("Operation failed") * k1[i]
172                    + F::zero() * k2[i]
173                    + F::from_f64(7571.0 / 16695.0).expect("Operation failed") * k3[i]
174                    + F::from_f64(393.0 / 640.0).expect("Operation failed") * k4[i]
175                    + F::from_f64(-92097.0 / 339200.0).expect("Operation failed") * k5[i]
176                    + F::from_f64(187.0 / 2100.0).expect("Operation failed") * k6[i]
177                    + F::from_f64(1.0 / 40.0).expect("Operation failed") * k7[i]);
178        }
179
180        // Error estimation
181        let mut err_norm = F::zero();
182        for i in 0..n_dim {
183            let sc = opts.atol + opts.rtol * y5[i].abs();
184            let err = (y5[i] - y4[i]).abs() / sc;
185            err_norm = err_norm.max(err);
186        }
187
188        // Step size control
189        let order = F::from_f64(5.0).expect("Operation failed"); // 5th order method
190        let exponent = F::one() / (order + F::one());
191        let safety = F::from_f64(0.9).expect("Operation failed");
192        let factor = safety * (F::one() / err_norm).powf(exponent);
193        let factor_min = F::from_f64(0.2).expect("Operation failed");
194        let factor_max = F::from_f64(5.0).expect("Operation failed");
195        let factor = factor.min(factor_max).max(factor_min);
196
197        if err_norm <= F::one() {
198            // Step accepted
199            t += h;
200            y = y5; // Use higher order solution
201
202            // Store results
203            t_values.push(t);
204            y_values.push(y.clone());
205
206            // Increase step size for next step
207            if err_norm <= F::from_f64(0.1).expect("Operation failed") {
208                // For very accurate steps, try a larger increase
209                h *= factor.max(F::from_f64(2.0).expect("Operation failed"));
210            } else {
211                h *= factor;
212            }
213
214            step_count += 1;
215            accepted_steps += 1;
216        } else {
217            // Step rejected
218            h *= factor.min(F::one());
219            rejected_steps += 1;
220
221            // If step size is too small, return error
222            if h < min_step {
223                return Err(crate::error::IntegrateError::StepSizeTooSmall(format!(
224                    "Step size {h} too small at t {t}"
225                )));
226            }
227        }
228    }
229
230    // Check if integration was successful
231    let success = t >= t_end;
232    let message = if !success {
233        Some(format!(
234            "Maximum number of steps ({}) reached",
235            opts.max_steps
236        ))
237    } else {
238        None
239    };
240
241    // Return the solution
242    Ok(ODEResult {
243        t: t_values,
244        y: y_values,
245        success,
246        message,
247        n_eval: func_evals,
248        n_steps: step_count,
249        n_accepted: accepted_steps,
250        n_rejected: rejected_steps,
251        n_lu: 0,  // No LU decompositions in explicit methods
252        n_jac: 0, // No Jacobian evaluations in explicit methods
253        method: ODEMethod::RK45,
254    })
255}
256
257/// Shared PI (proportional-integral) step-size controller in the style of
258/// Hairer, Nørsett & Wanner's "Lund stabilization" (see e.g. the reference
259/// `dopri5.f`/`dop853.f` codes and Gustafsson (1991), "Control theoretic
260/// techniques for stepsize selection in explicit Runge-Kutta methods",
261/// ACM TOMS 17(4)).
262///
263/// The controller combines the current step's tolerance-normalized local
264/// error norm (`err_norm`, where `<= 1` means "accept") with the *previous
265/// accepted* step's error norm (`err_prev`) to damp the step-size
266/// oscillations that a purely elementary (integral-only) controller can
267/// exhibit. `beta_gain = 0` degenerates it to a plain elementary
268/// controller (used for DOP853, matching the reference implementation's
269/// own default of no Lund stabilization).
270///
271/// Returns `(step_size_factor, updated_err_prev)`, where `step_size_factor`
272/// is the multiplier to apply to `h`. `err_prev` is only advanced when
273/// `accepted` is true, mirroring the reference behavior of only updating
274/// the controller's memory on accepted steps.
275#[allow(clippy::too_many_arguments)]
276fn pi_step_factor<F: IntegrateFloat>(
277    err_norm: F,
278    err_prev: F,
279    accepted: bool,
280    just_recovered_from_rejection: bool,
281    alpha: F,
282    beta_gain: F,
283    safety: F,
284    growth_min: F,
285    growth_max: F,
286) -> (F, F) {
287    let tiny = F::from_f64(1e-300).expect("Operation failed");
288    let err_eff = err_norm.max(tiny);
289    let fac_elementary = err_eff.powf(alpha);
290
291    if accepted {
292        let history = err_prev.max(tiny).powf(beta_gain);
293        let mut factor = safety / (fac_elementary * history);
294        factor = factor.min(growth_max).max(growth_min);
295        if just_recovered_from_rejection {
296            // Never grow immediately after recovering from a rejection.
297            factor = factor.min(F::one());
298        }
299        let updated_err_prev = err_norm.max(F::from_f64(1e-10).expect("Operation failed"));
300        (factor, updated_err_prev)
301    } else {
302        // On rejection, ignore step-size history entirely (matches the
303        // reference codes) and guarantee an actual shrink.
304        let factor = (safety / fac_elementary).min(F::one()).max(growth_min);
305        (factor, err_prev)
306    }
307}
308
309/// Solve ODE using the Bogacki-Shampine method (RK23)
310///
311/// This is an adaptive step size method based on embedded Runge-Kutta
312/// formulas. It uses the Bogacki & Shampine (1989) 3(2) pair: a 3rd-order
313/// solution is advanced at every step, with a 2nd-order embedded solution
314/// used purely for error estimation (local extrapolation). The method has
315/// the "first same as last" (FSAL) property: the 4th stage is evaluated at
316/// the accepted 3rd-order solution, so the embedded error estimate comes at
317/// the cost of only one extra function evaluation per step. Step sizes are
318/// chosen with a PI (proportional-integral) controller.
319///
320/// # Arguments
321///
322/// * `f` - ODE function dy/dt = f(t, y)
323/// * `t_span` - Time span [t_start, t_end]
324/// * `y0` - Initial condition
325/// * `opts` - Solver options
326///
327/// # Returns
328///
329/// The solution as an ODEResult or an error
330///
331/// # References
332///
333/// P. Bogacki, L.F. Shampine, "A 3(2) Pair of Runge-Kutta Formulas",
334/// Appl. Math. Lett. Vol. 2, No. 4, pp. 321-325, 1989.
335#[allow(dead_code)]
336pub fn rk23_method<F, Func>(
337    f: Func,
338    t_span: [F; 2],
339    y0: Array1<F>,
340    opts: ODEOptions<F>,
341) -> IntegrateResult<ODEResult<F>>
342where
343    F: IntegrateFloat,
344    Func: Fn(F, ArrayView1<F>) -> Array1<F>,
345{
346    // Initialize
347    let [t_start, t_end] = t_span;
348    let n_dim = y0.len();
349
350    // Determine initial step size if not provided
351    let h0 = opts.h0.unwrap_or_else(|| {
352        // Simple heuristic for initial step size
353        let _span = t_end - t_start;
354        _span / F::from_usize(100).expect("Operation failed")
355    });
356
357    // Determine minimum and maximum step sizes
358    let min_step = opts.min_step.unwrap_or_else(|| {
359        let _span = t_end - t_start;
360        _span * F::from_f64(1e-8).expect("Operation failed") // Minimal step size
361    });
362
363    let max_step = opts.max_step.unwrap_or_else(|| {
364        t_end - t_start // Maximum step can be the whole interval
365    });
366
367    // Current state
368    let mut t = t_start;
369    let mut y = y0.clone();
370    let mut h = h0;
371
372    // Storage for results
373    let mut t_values = vec![t_start];
374    let mut y_values = vec![y0.clone()];
375
376    // Statistics
377    let mut func_evals = 0;
378    let mut step_count = 0;
379    let mut accepted_steps = 0;
380    let mut rejected_steps = 0;
381
382    // Bogacki-Shampine (1989) coefficients.
383    let c2 = F::from_f64(0.5).expect("Operation failed");
384    let c3 = F::from_f64(0.75).expect("Operation failed");
385    let a21 = F::from_f64(0.5).expect("Operation failed");
386    let a32 = F::from_f64(0.75).expect("Operation failed");
387    let b1 = F::from_f64(2.0 / 9.0).expect("Operation failed");
388    let b2 = F::from_f64(1.0 / 3.0).expect("Operation failed");
389    let b3 = F::from_f64(4.0 / 9.0).expect("Operation failed");
390    // Embedded error-estimator weights: err = h * (e1*k1+e2*k2+e3*k3+e4*k4)
391    let e1 = F::from_f64(5.0 / 72.0).expect("Operation failed");
392    let e2 = F::from_f64(-1.0 / 12.0).expect("Operation failed");
393    let e3 = F::from_f64(-1.0 / 9.0).expect("Operation failed");
394    let e4 = F::from_f64(0.125).expect("Operation failed");
395
396    // PI step-size controller parameters. `alpha` is the elementary
397    // exponent 1/(error_order+1) with error_order=2 (the embedded
398    // 2nd-order solution); `beta_gain` is the PI (history feedback) gain.
399    let alpha = F::one() / F::from_f64(3.0).expect("Operation failed");
400    let beta_gain = F::from_f64(0.08).expect("Operation failed");
401    let safety = F::from_f64(0.9).expect("Operation failed");
402    let growth_min = F::from_f64(0.2).expect("Operation failed");
403    let growth_max = F::from_f64(10.0).expect("Operation failed");
404    let mut err_prev = F::one();
405    let mut just_rejected = false;
406
407    // Main integration loop
408    while t < t_end && step_count < opts.max_steps {
409        // Adjust step size for the last step if needed
410        if t + h > t_end {
411            h = t_end - t;
412        }
413
414        // Limit step size to bounds
415        h = h.min(max_step).max(min_step);
416
417        // Stage 1
418        let k1 = f(t, y.view());
419
420        // Stage 2
421        let mut y_stage = y.clone();
422        for i in 0..n_dim {
423            y_stage[i] = y[i] + h * a21 * k1[i];
424        }
425        let k2 = f(t + c2 * h, y_stage.view());
426
427        // Stage 3
428        let mut y_stage = y.clone();
429        for i in 0..n_dim {
430            y_stage[i] = y[i] + h * a32 * k2[i];
431        }
432        let k3 = f(t + c3 * h, y_stage.view());
433
434        // 3rd order solution (the point at which the FSAL stage is
435        // evaluated).
436        let mut y3 = y.clone();
437        for i in 0..n_dim {
438            y3[i] = y[i] + h * (b1 * k1[i] + b2 * k2[i] + b3 * k3[i]);
439        }
440
441        // FSAL stage: derivative at the accepted 3rd-order point.
442        let k4 = f(t + h, y3.view());
443        func_evals += 4;
444
445        // Embedded error estimate (2nd order vs. 3rd order), WRMS-style.
446        let mut err_norm = F::zero();
447        for i in 0..n_dim {
448            let sc = opts.atol + opts.rtol * y3[i].abs().max(y[i].abs());
449            let err_i = e1 * k1[i] + e2 * k2[i] + e3 * k3[i] + e4 * k4[i];
450            err_norm = err_norm.max((h * err_i / sc).abs());
451        }
452
453        let accepted = err_norm <= F::one();
454        let (factor, new_err_prev) = pi_step_factor(
455            err_norm,
456            err_prev,
457            accepted,
458            just_rejected,
459            alpha,
460            beta_gain,
461            safety,
462            growth_min,
463            growth_max,
464        );
465
466        if accepted {
467            // Step accepted
468            t += h;
469            y = y3;
470
471            // Store results
472            t_values.push(t);
473            y_values.push(y.clone());
474
475            err_prev = new_err_prev;
476            just_rejected = false;
477            h *= factor;
478
479            step_count += 1;
480            accepted_steps += 1;
481        } else {
482            // Step rejected
483            h *= factor;
484            just_rejected = true;
485            rejected_steps += 1;
486
487            // If step size is too small, return error
488            if h < min_step {
489                return Err(crate::error::IntegrateError::StepSizeTooSmall(format!(
490                    "Step size {h} too small at t {t}"
491                )));
492            }
493        }
494    }
495
496    // Check if integration was successful
497    let success = t >= t_end;
498    let message = if !success {
499        Some(format!(
500            "Maximum number of steps ({}) reached",
501            opts.max_steps
502        ))
503    } else {
504        None
505    };
506
507    // Return the solution
508    Ok(ODEResult {
509        t: t_values,
510        y: y_values,
511        success,
512        message,
513        n_eval: func_evals,
514        n_steps: step_count,
515        n_accepted: accepted_steps,
516        n_rejected: rejected_steps,
517        n_lu: 0,  // No LU decompositions in explicit methods
518        n_jac: 0, // No Jacobian evaluations in explicit methods
519        method: ODEMethod::RK23,
520    })
521}
522
523/// Number of stages used to advance the DOP853 solution.
524///
525/// The reference `dop853.f`/`dop853_coefficients.py` tableau has 16 rows:
526/// 12 stages for the main 8th/5th/3rd order solution and error estimate,
527/// a 13th "first same as last" (FSAL) stage, and 3 more stages used only to
528/// build a dense-output (continuous extension) interpolant. Neither error
529/// estimator weight vector (`DOP853_E5`, the `E3` derived from
530/// `DOP853_B`) has a nonzero coefficient on the FSAL stage, and this
531/// crate's ODE solver contract reconstructs dense/continuous output
532/// externally (see `ode::utils::dense_output::DenseSolution`, built via
533/// cubic Hermite interpolation over the returned samples) rather than
534/// through a per-step interpolant -- exactly the same contract already
535/// used by `rk45_method` in this file. The FSAL and dense-output-only
536/// stages are therefore omitted here as genuinely unneeded evaluations.
537const DOP853_STAGES: usize = 12;
538
539/// c_i: stage time fractions for the DOP853 tableau (Hairer, Nørsett &
540/// Wanner, "Solving Ordinary Differential Equations I", 2nd ed.).
541const DOP853_C: [f64; 12] = [
542    0.0,
543    0.05260015195876773,
544    0.0789002279381516,
545    0.1183503419072274,
546    0.2816496580927726,
547    0.3333333333333333,
548    0.25,
549    0.3076923076923077,
550    0.6512820512820513,
551    0.6,
552    0.8571428571428571,
553    1.0,
554];
555
556/// a_{i,j} strictly-lower-triangular stage coefficients. Row `i` (0-indexed)
557/// holds the weights on stages `1..=i` used to build stage `i+1`'s
558/// argument; row 0 (stage 1, always `f(t, y)`) needs no combination.
559const DOP853_A: [[f64; 12]; 12] = [
560    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
561    [
562        0.05260015195876773,
563        0.0,
564        0.0,
565        0.0,
566        0.0,
567        0.0,
568        0.0,
569        0.0,
570        0.0,
571        0.0,
572        0.0,
573        0.0,
574    ],
575    [
576        0.0197250569845379,
577        0.0591751709536137,
578        0.0,
579        0.0,
580        0.0,
581        0.0,
582        0.0,
583        0.0,
584        0.0,
585        0.0,
586        0.0,
587        0.0,
588    ],
589    [
590        0.02958758547680685,
591        0.0,
592        0.08876275643042054,
593        0.0,
594        0.0,
595        0.0,
596        0.0,
597        0.0,
598        0.0,
599        0.0,
600        0.0,
601        0.0,
602    ],
603    [
604        0.2413651341592667,
605        0.0,
606        -0.8845494793282861,
607        0.924834003261792,
608        0.0,
609        0.0,
610        0.0,
611        0.0,
612        0.0,
613        0.0,
614        0.0,
615        0.0,
616    ],
617    [
618        0.037037037037037035,
619        0.0,
620        0.0,
621        0.17082860872947386,
622        0.12546768756682242,
623        0.0,
624        0.0,
625        0.0,
626        0.0,
627        0.0,
628        0.0,
629        0.0,
630    ],
631    [
632        0.037109375,
633        0.0,
634        0.0,
635        0.17025221101954405,
636        0.06021653898045596,
637        -0.017578125,
638        0.0,
639        0.0,
640        0.0,
641        0.0,
642        0.0,
643        0.0,
644    ],
645    [
646        0.03709200011850479,
647        0.0,
648        0.0,
649        0.17038392571223998,
650        0.10726203044637328,
651        -0.015319437748624402,
652        0.008273789163814023,
653        0.0,
654        0.0,
655        0.0,
656        0.0,
657        0.0,
658    ],
659    [
660        0.6241109587160757,
661        0.0,
662        0.0,
663        -3.3608926294469414,
664        -0.868219346841726,
665        27.59209969944671,
666        20.154067550477894,
667        -43.48988418106996,
668        0.0,
669        0.0,
670        0.0,
671        0.0,
672    ],
673    [
674        0.47766253643826434,
675        0.0,
676        0.0,
677        -2.4881146199716677,
678        -0.590290826836843,
679        21.230051448181193,
680        15.279233632882423,
681        -33.28821096898486,
682        -0.020331201708508627,
683        0.0,
684        0.0,
685        0.0,
686    ],
687    [
688        -0.9371424300859873,
689        0.0,
690        0.0,
691        5.186372428844064,
692        1.0914373489967295,
693        -8.149787010746927,
694        -18.52006565999696,
695        22.739487099350505,
696        2.4936055526796523,
697        -3.0467644718982196,
698        0.0,
699        0.0,
700    ],
701    [
702        2.273310147516538,
703        0.0,
704        0.0,
705        -10.53449546673725,
706        -2.0008720582248625,
707        -17.9589318631188,
708        27.94888452941996,
709        -2.8589982771350235,
710        -8.87285693353063,
711        12.360567175794303,
712        0.6433927460157636,
713        0.0,
714    ],
715];
716
717/// b_i: 8th-order solution weights.
718const DOP853_B: [f64; 12] = [
719    0.054293734116568765,
720    0.0,
721    0.0,
722    0.0,
723    0.0,
724    4.450312892752409,
725    1.8915178993145003,
726    -5.801203960010585,
727    0.3111643669578199,
728    -0.1521609496625161,
729    0.20136540080403034,
730    0.04471061572777259,
731];
732
733/// 5th-order error-estimator weights (`err5 = dot(K, E5)`).
734const DOP853_E5: [f64; 12] = [
735    0.01312004499419488,
736    0.0,
737    0.0,
738    0.0,
739    0.0,
740    -1.2251564463762044,
741    -0.4957589496572502,
742    1.6643771824549864,
743    -0.35032884874997366,
744    0.3341791187130175,
745    0.08192320648511571,
746    -0.022355307863886294,
747];
748
749/// Tabulated adjustments turning the 8th-order weight vector `B` into the
750/// secondary 3rd-order error-check weight vector `E3`: only 3 of the 12
751/// stage weights differ between the two (Hairer, Nørsett & Wanner).
752const DOP853_E3_ADJUST: [(usize, f64); 3] = [
753    (0, 0.2440944881889764),
754    (8, 0.7338466882816118),
755    (11, 0.022058823529411766),
756];
757
758/// Solve ODE using the Dormand-Prince 8th order method (DOP853)
759///
760/// This is a high-accuracy adaptive step size method based on embedded
761/// Runge-Kutta formulas: the full Hairer-Nørsett-Wanner 12-stage 8th-order
762/// method, with a blended 5th/3rd-order embedded error estimate. The
763/// blend combines a primary 5th-order error estimator with a secondary
764/// 3rd-order check so that the estimate stays reliable even in the rare
765/// case where the 5th-order estimator alone would (spuriously) nearly
766/// vanish. Step sizes are chosen with the same PI-capable controller used
767/// by [`rk23_method`], with the PI (history) term disabled by default to
768/// match the reference implementation's own default (no Lund
769/// stabilization for DOP853 unless explicitly requested).
770///
771/// # Arguments
772///
773/// * `f` - ODE function dy/dt = f(t, y)
774/// * `t_span` - Time span [t_start, t_end]
775/// * `y0` - Initial condition
776/// * `opts` - Solver options
777///
778/// # Returns
779///
780/// The solution as an ODEResult or an error
781///
782/// # References
783///
784/// E. Hairer, S.P. Nørsett, G. Wanner, "Solving Ordinary Differential
785/// Equations I: Nonstiff Problems", 2nd ed., Springer, 1993, Section II.10.
786#[allow(dead_code)]
787pub fn dop853_method<F, Func>(
788    f: Func,
789    t_span: [F; 2],
790    y0: Array1<F>,
791    opts: ODEOptions<F>,
792) -> IntegrateResult<ODEResult<F>>
793where
794    F: IntegrateFloat,
795    Func: Fn(F, ArrayView1<F>) -> Array1<F>,
796{
797    // Initialize
798    let [t_start, t_end] = t_span;
799    let n_dim = y0.len();
800
801    // Determine initial step size if not provided
802    let h0 = opts.h0.unwrap_or_else(|| {
803        // Simple heuristic for initial step size
804        let _span = t_end - t_start;
805        _span / F::from_usize(100).expect("Operation failed")
806    });
807
808    // Determine minimum and maximum step sizes
809    let min_step = opts.min_step.unwrap_or_else(|| {
810        let _span = t_end - t_start;
811        _span * F::from_f64(1e-8).expect("Operation failed") // Minimal step size
812    });
813
814    let max_step = opts.max_step.unwrap_or_else(|| {
815        t_end - t_start // Maximum step can be the whole interval
816    });
817
818    // Current state
819    let mut t = t_start;
820    let mut y = y0.clone();
821    let mut h = h0;
822
823    // Storage for results
824    let mut t_values = vec![t_start];
825    let mut y_values = vec![y0.clone()];
826
827    // Statistics
828    let mut func_evals = 0;
829    let mut step_count = 0;
830    let mut accepted_steps = 0;
831    let mut rejected_steps = 0;
832
833    // Secondary 3rd-order check weight vector: E3 = B with 3 tabulated
834    // adjustments (Hairer, Nørsett & Wanner).
835    let mut e3 = DOP853_B;
836    for &(idx, adjust) in DOP853_E3_ADJUST.iter() {
837        e3[idx] -= adjust;
838    }
839
840    // Step-size controller parameters. `alpha` is the reference
841    // implementation's `expo1` with its default `beta = 0` (i.e. no
842    // Lund/PI stabilization for DOP853 by default).
843    let alpha = F::from_f64(0.125).expect("Operation failed");
844    let beta_gain = F::zero();
845    let safety = F::from_f64(0.9).expect("Operation failed");
846    let growth_min = F::from_f64(1.0 / 3.0).expect("Operation failed");
847    let growth_max = F::from_f64(6.0).expect("Operation failed");
848    let mut err_prev = F::one();
849    let mut just_rejected = false;
850
851    let n_dim_f = F::from_usize(n_dim).expect("Operation failed");
852    let err3_weight = F::from_f64(0.01).expect("Operation failed");
853    let non_finite_penalty = F::from_f64(1e30).expect("Operation failed");
854
855    // Main integration loop
856    while t < t_end && step_count < opts.max_steps {
857        // Adjust step size for the last step if needed
858        if t + h > t_end {
859            h = t_end - t;
860        }
861
862        // Limit step size to bounds
863        h = h.min(max_step).max(min_step);
864
865        // Compute the 12 stages of the DOP853 tableau.
866        let mut k_stages: Vec<Array1<F>> = Vec::with_capacity(DOP853_STAGES);
867        k_stages.push(f(t, y.view()));
868        for s in 1..DOP853_STAGES {
869            let mut y_stage = y.clone();
870            for (j, k_j) in k_stages.iter().enumerate().take(s) {
871                let a_sj = DOP853_A[s][j];
872                if a_sj != 0.0 {
873                    let a_sj = F::from_f64(a_sj).expect("Operation failed");
874                    for d in 0..n_dim {
875                        y_stage[d] += h * a_sj * k_j[d];
876                    }
877                }
878            }
879            let c_s = F::from_f64(DOP853_C[s]).expect("Operation failed");
880            k_stages.push(f(t + c_s * h, y_stage.view()));
881        }
882        func_evals += DOP853_STAGES;
883
884        // 8th-order solution.
885        let mut y8 = y.clone();
886        for (j, k_j) in k_stages.iter().enumerate() {
887            let b_j = DOP853_B[j];
888            if b_j != 0.0 {
889                let b_j = F::from_f64(b_j).expect("Operation failed");
890                for d in 0..n_dim {
891                    y8[d] += h * b_j * k_j[d];
892                }
893            }
894        }
895
896        // Blended 5th/3rd-order error estimate.
897        let mut err5_sq = F::zero();
898        let mut err3_sq = F::zero();
899        for d in 0..n_dim {
900            let scale = opts.atol + opts.rtol * y[d].abs().max(y8[d].abs());
901            let mut e5_d = F::zero();
902            let mut e3_d = F::zero();
903            for (j, k_j) in k_stages.iter().enumerate() {
904                let kjd = k_j[d];
905                let e5j = DOP853_E5[j];
906                if e5j != 0.0 {
907                    e5_d += F::from_f64(e5j).expect("Operation failed") * kjd;
908                }
909                let e3j = e3[j];
910                if e3j != 0.0 {
911                    e3_d += F::from_f64(e3j).expect("Operation failed") * kjd;
912                }
913            }
914            let e5s = e5_d / scale;
915            let e3s = e3_d / scale;
916            err5_sq += e5s * e5s;
917            err3_sq += e3s * e3s;
918        }
919        let denom = err5_sq + err3_weight * err3_sq;
920        let err_norm = if denom <= F::zero() {
921            F::zero()
922        } else {
923            let raw = h.abs() * err5_sq / (denom * n_dim_f).sqrt();
924            if raw.is_finite() {
925                raw
926            } else {
927                // A diverging/blown-up derivative produced a non-finite
928                // error estimate; treat it as maximally bad so the
929                // controller shrinks the step aggressively instead of
930                // stalling.
931                non_finite_penalty
932            }
933        };
934
935        let accepted = err_norm <= F::one();
936        let (factor, new_err_prev) = pi_step_factor(
937            err_norm,
938            err_prev,
939            accepted,
940            just_rejected,
941            alpha,
942            beta_gain,
943            safety,
944            growth_min,
945            growth_max,
946        );
947
948        if accepted {
949            // Step accepted
950            t += h;
951            y = y8;
952
953            // Store results
954            t_values.push(t);
955            y_values.push(y.clone());
956
957            err_prev = new_err_prev;
958            just_rejected = false;
959            h *= factor;
960
961            step_count += 1;
962            accepted_steps += 1;
963        } else {
964            // Step rejected
965            h *= factor;
966            just_rejected = true;
967            rejected_steps += 1;
968
969            // If step size is too small, return error
970            if h < min_step {
971                return Err(crate::error::IntegrateError::StepSizeTooSmall(format!(
972                    "Step size {h} too small at t {t}"
973                )));
974            }
975        }
976    }
977
978    // Check if integration was successful
979    let success = t >= t_end;
980    let message = if !success {
981        Some(format!(
982            "Maximum number of steps ({}) reached",
983            opts.max_steps
984        ))
985    } else {
986        None
987    };
988
989    // Return the solution
990    Ok(ODEResult {
991        t: t_values,
992        y: y_values,
993        success,
994        message,
995        n_eval: func_evals,
996        n_steps: step_count,
997        n_accepted: accepted_steps,
998        n_rejected: rejected_steps,
999        n_lu: 0,  // No LU decompositions in explicit methods
1000        n_jac: 0, // No Jacobian evaluations in explicit methods
1001        method: ODEMethod::DOP853,
1002    })
1003}
1004
1005#[cfg(test)]
1006mod tests {
1007    use super::*;
1008
1009    /// Classical fixed-step convergence-order check: setting
1010    /// `min_step == max_step == h` pins every attempted step to exactly
1011    /// `h`, since the `h.min(max_step).max(min_step)` clamp at the top of
1012    /// the main loop wins regardless of what the adaptive controller's
1013    /// accept/reject decision would otherwise have chosen. A very loose
1014    /// tolerance guarantees every fixed-size step is accepted, so this
1015    /// isolates the raw Butcher tableau's truncation-error order: halving
1016    /// `h` should shrink the global error by ~2^order.
1017    #[test]
1018    fn rk23_empirical_convergence_order_is_three() {
1019        let k = 3.0_f64;
1020        let f =
1021            move |_t: f64, y: ArrayView1<f64>| -> Array1<f64> { Array1::from_vec(vec![-k * y[0]]) };
1022        let t_end = 1.0_f64;
1023        let y0 = 2.0_f64;
1024        let y_exact = y0 * (-k * t_end).exp();
1025
1026        let mut prev_err: Option<f64> = None;
1027        let mut min_ratio = f64::INFINITY;
1028        let mut h = t_end / 8.0;
1029        for _ in 0..3 {
1030            let opts = ODEOptions {
1031                method: ODEMethod::RK23,
1032                rtol: 1.0,
1033                atol: 1.0,
1034                h0: Some(h),
1035                min_step: Some(h),
1036                max_step: Some(h),
1037                max_steps: 1_000_000,
1038                ..Default::default()
1039            };
1040            let result = rk23_method(f, [0.0, t_end], Array1::from_vec(vec![y0]), opts)
1041                .expect("RK23 fixed-step solve failed");
1042            assert!(
1043                result.success,
1044                "RK23 fixed-step integration did not complete"
1045            );
1046            let y_final = result.y.last().expect("empty result")[0];
1047            let err = (y_final - y_exact).abs();
1048
1049            if let Some(prev) = prev_err {
1050                min_ratio = min_ratio.min(prev / err);
1051            }
1052            prev_err = Some(err);
1053            h /= 2.0;
1054        }
1055
1056        // A genuine 3rd-order method shrinks the error by ~2^3=8x each
1057        // time h halves. Use a generous lower bound (5x) that a
1058        // mislabeled 1st-order Euler stand-in (ratio ~2x) cannot meet.
1059        assert!(
1060            min_ratio > 5.0,
1061            "RK23 empirical convergence order too low: min ratio={min_ratio} (expected ~8)"
1062        );
1063    }
1064
1065    #[test]
1066    fn dop853_empirical_convergence_order_is_eight() {
1067        let k = 3.0_f64;
1068        let f =
1069            move |_t: f64, y: ArrayView1<f64>| -> Array1<f64> { Array1::from_vec(vec![-k * y[0]]) };
1070        let t_end = 1.0_f64;
1071        let y0 = 2.0_f64;
1072        let y_exact = y0 * (-k * t_end).exp();
1073
1074        let mut prev_err: Option<f64> = None;
1075        let mut min_ratio = f64::INFINITY;
1076        let mut h = t_end / 2.0;
1077        for _ in 0..3 {
1078            let opts = ODEOptions {
1079                method: ODEMethod::DOP853,
1080                rtol: 1.0,
1081                atol: 1.0,
1082                h0: Some(h),
1083                min_step: Some(h),
1084                max_step: Some(h),
1085                max_steps: 1_000_000,
1086                ..Default::default()
1087            };
1088            let result = dop853_method(f, [0.0, t_end], Array1::from_vec(vec![y0]), opts)
1089                .expect("DOP853 fixed-step solve failed");
1090            assert!(
1091                result.success,
1092                "DOP853 fixed-step integration did not complete"
1093            );
1094            let y_final = result.y.last().expect("empty result")[0];
1095            let err = (y_final - y_exact).abs();
1096
1097            if let Some(prev) = prev_err {
1098                min_ratio = min_ratio.min(prev / err);
1099            }
1100            prev_err = Some(err);
1101            h /= 2.0;
1102        }
1103
1104        // A genuine 8th-order method shrinks the error by ~2^8=256x each
1105        // time h halves. Use a generous lower bound (50x).
1106        assert!(
1107            min_ratio > 50.0,
1108            "DOP853 empirical convergence order too low: min ratio={min_ratio} (expected ~256)"
1109        );
1110    }
1111
1112    /// Cross-method agreement + sanity check directly at the method-function
1113    /// level (bypassing the public `solve_ivp` dispatcher): RK23 and DOP853
1114    /// must agree closely with the crate's existing RK45 on a nonlinear,
1115    /// non-constant-trajectory problem (Van der Pol oscillator).
1116    #[test]
1117    fn rk23_dop853_agree_with_rk45_van_der_pol() {
1118        let mu = 1.0_f64;
1119        let f = move |_t: f64, y: ArrayView1<f64>| -> Array1<f64> {
1120            Array1::from_vec(vec![y[1], mu * (1.0 - y[0] * y[0]) * y[1] - y[0]])
1121        };
1122        let y0 = Array1::from_vec(vec![0.5_f64, 0.0]);
1123        let t_span = [0.0_f64, 3.0];
1124
1125        let opts_of = |method: ODEMethod| ODEOptions {
1126            method,
1127            rtol: 1e-10,
1128            atol: 1e-12,
1129            max_steps: 500_000,
1130            ..Default::default()
1131        };
1132
1133        let rk45 = rk45_method(f, t_span, y0.clone(), opts_of(ODEMethod::RK45))
1134            .expect("RK45 Van der Pol solve failed");
1135        let rk23 = rk23_method(f, t_span, y0.clone(), opts_of(ODEMethod::RK23))
1136            .expect("RK23 Van der Pol solve failed");
1137        let dop853 = dop853_method(f, t_span, y0, opts_of(ODEMethod::DOP853))
1138            .expect("DOP853 Van der Pol solve failed");
1139
1140        assert!(rk45.success && rk23.success && dop853.success);
1141
1142        let y_rk45 = rk45.y.last().expect("empty result");
1143        let y_rk23 = rk23.y.last().expect("empty result");
1144        let y_dop853 = dop853.y.last().expect("empty result");
1145
1146        for i in 0..2 {
1147            assert!(
1148                (y_rk23[i] - y_rk45[i]).abs() < 1e-5,
1149                "RK23 vs RK45 disagreement at index {i}: {} vs {}",
1150                y_rk23[i],
1151                y_rk45[i]
1152            );
1153            assert!(
1154                (y_dop853[i] - y_rk45[i]).abs() < 1e-5,
1155                "DOP853 vs RK45 disagreement at index {i}: {} vs {}",
1156                y_dop853[i],
1157                y_rk45[i]
1158            );
1159        }
1160    }
1161}
scirs2_integrate/ode/methods/adaptive.rs

scirs2_integrate/ode/methods/
adaptive.rs
