diff --git a/Cargo.toml b/Cargo.toml index 94bf882..1e9f68b 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -24,6 +24,7 @@ roots = "0.0.8" approx = "0.5" num-dual = "0.15" criterion = "0.7.0" +serde_json = "1.0" [[bench]] name = "simple_1d" diff --git a/readme.md b/readme.md index 76e6b58..ffb7480 100644 --- a/readme.md +++ b/readme.md @@ -29,6 +29,7 @@ See [benchmark report](VERN7_BENCHMARK_REPORT.md) for detailed performance analy - **Lazy computation** of extra interpolation stages (Vern7) - **PID step size control** as an alternative to the PI controller - **Automatic Differentiation**, the same problem can be solved with floats or also differentiated with dual numbers +- **Runtime integrator selection**, the method can be a value, named in a configuration file rather than fixed as a type ### Future Improvements @@ -83,6 +84,36 @@ let solution = problem.solve(); let _interpolated_answer = solution.interpolate(4.4); ``` +## Choosing the integrator at runtime + +The integrator is normally a type parameter, fixed when the code is compiled. `Method` is an enum +over the integrators here that itself implements `Integrator`, so the choice can be a value instead: + +```rust +use nalgebra::Vector1; +use ordinary_diffeq::prelude::*; + +fn derivative(_t: f64, y: Vector1, _p: &()) -> Vector1 { + -y +} + +// `MethodName` derives `Serialize`/`Deserialize`, so this can come out of a +// configuration file instead of being written here. +let name: MethodName = serde_json::from_str(r#""Vern7""#).unwrap(); +let method: Method<1> = name.instantiate().a_tol(1e-12).r_tol(1e-12); + +let ode = ODE::new(&derivative, 0.0, 2.0, Vector1::new(1.0), ()); +let mut problem = Problem::new(ode, method, PIController::default()); +let solution = problem.solve(); + +// And the method that ran can be recorded in a form that selects it again +assert_eq!(method.name(), MethodName::Vern7); +``` + +Dispatch is a `match` on the variant, not a trait object, so nothing is erased: the same `Method` +works at a dual scalar and differentiates exactly as the integrator it wraps does. A `Vec>` +holds several methods for comparison over the same problem. + ## Automatic differentiation State, time and step size are generic over the `Scalar` trait. `f64` implements it by default, diff --git a/roadmap/README.md b/roadmap/README.md index a551874..46c0a05 100644 --- a/roadmap/README.md +++ b/roadmap/README.md @@ -68,12 +68,14 @@ Each feature below links to a detailed implementation plan in the `features/` di - **Effort**: Small - **Status**: All four integrators implement `Tolerances`; no call site changed -- [ ] **[Runtime Integrator Selection](features/41-runtime-integrator-selection.md)** +- [x] **[Runtime Integrator Selection](features/41-runtime-integrator-selection.md)** ✅ COMPLETED - Structural facts move from associated constants to `&self` methods, so an enum can dispatch and an integrator can be chosen at runtime or named in a config file - **Dependencies**: None (pairs with Tolerances Trait) - **Effort**: Medium + - **Status**: `Method` and `MethodName` added; benches measured at -0.7% to + +0.4%, within criterion's noise threshold ### Callbacks @@ -354,11 +356,11 @@ Each algorithm implementation should include: ## Progress Tracking Total Features: 41 -- Tier 1: 10 features (5/10 complete) ✅ +- Tier 1: 10 features (6/10 complete) ✅ - Tier 2: 11 features (0/11 complete) - Tier 3: 20 features (0/20 complete) -**Overall Progress: 12.2% (5/41 features complete)** +**Overall Progress: 14.6% (6/41 features complete)** ### Completed Features 1. ✅ BS3 (Bogacki-Shampine 3/2) - Tier 1 (2025-10-23) @@ -366,5 +368,6 @@ Total Features: 41 3. ✅ PID Controller - Tier 1 (2025-10-24) 4. ✅ Rosenbrock23 - Tier 1 (2025-10-24) 5. ✅ Tolerances Trait - Tier 1 (2026-09-30) +6. ✅ Runtime Integrator Selection - Tier 1 (2026-09-30) Last updated: 2026-09-30 diff --git a/roadmap/features/41-runtime-integrator-selection.md b/roadmap/features/41-runtime-integrator-selection.md index aa19463..aaf72eb 100644 --- a/roadmap/features/41-runtime-integrator-selection.md +++ b/roadmap/features/41-runtime-integrator-selection.md @@ -1,6 +1,10 @@ # Feature: Runtime Integrator Selection -**Status**: proposed (2026-09-30) +**Status**: completed (2026-09-30) + +Measured on the existing criterion benches, constants to methods: `bench_orbit` +-0.73%, `bench_simple_1d` +0.42%, `bench_interpolation_1d` -0.09%. Criterion +calls all three within its noise threshold, so the inlining claim below holds. ## Overview diff --git a/src/integrator/bs3.rs b/src/integrator/bs3.rs index 1a1da76..abce274 100644 --- a/src/integrator/bs3.rs +++ b/src/integrator/bs3.rs @@ -4,6 +4,9 @@ use super::super::ode::ODE; use super::{Integrator, StepResult, Tolerances}; use crate::scalar::Scalar; +/// The number of stages BS3 evaluates per step. +const STAGES: usize = 4; + /// Bogacki-Shampine 3/2 integrator trait for tableau coefficients pub trait BS3Integrator<'a> { const A: &'a [f64]; @@ -155,14 +158,25 @@ impl<'a, const D: usize> Integrator for BS3 where BS3: BS3Integrator<'a>, { - const ORDER: usize = 3; - const STAGES: usize = 4; - const ADAPTIVE: bool = true; - const DENSE: bool = true; + fn order(&self) -> usize { + 3 + } + + fn stages(&self) -> usize { + STAGES + } + + fn adaptive(&self) -> bool { + true + } + + fn dense(&self) -> bool { + true + } fn step(&self, ode: &ODE, h: T) -> StepResult { // Allocate storage for the 4 stages - let mut k: Vec> = vec![SVector::::zeros(); Self::STAGES]; + let mut k: Vec> = vec![SVector::::zeros(); STAGES]; // Stage 1: k1 = f(t, y) k[0] = (ode.f)(ode.t, ode.y, &ode.params); @@ -271,13 +285,12 @@ mod tests { use nalgebra::Vector1; #[test] - #[allow(clippy::assertions_on_constants)] fn test_bs3_creation() { - let _bs3: BS3<1> = BS3::new(); - assert_eq!(BS3::<1>::ORDER, 3); - assert_eq!(BS3::<1>::STAGES, 4); - assert!(BS3::<1>::ADAPTIVE); - assert!(BS3::<1>::DENSE); + let bs3: BS3<1> = BS3::new(); + assert_eq!(bs3.order(), 3); + assert_eq!(bs3.stages(), 4); + assert!(bs3.adaptive()); + assert!(bs3.dense()); } #[test] diff --git a/src/integrator/dormand_prince.rs b/src/integrator/dormand_prince.rs index ec2e799..000893b 100644 --- a/src/integrator/dormand_prince.rs +++ b/src/integrator/dormand_prince.rs @@ -4,6 +4,10 @@ use super::{Integrator, StepResult, Tolerances}; use crate::ode::ODE; use crate::scalar::Scalar; +/// The number of stages Dormand-Prince 5(4) evaluates per step. Also the offset +/// between the 5th and 4th order weights in `B`. +const STAGES: usize = 7; + /// Integrator Trait pub trait DormandPrinceIntegrator<'a> { const A: &'a [f64]; @@ -174,24 +178,35 @@ impl<'a, const D: usize> Integrator for DormandPrince45 where DormandPrince45: DormandPrinceIntegrator<'a>, { - const ORDER: usize = 5; - const STAGES: usize = 7; - const ADAPTIVE: bool = true; - const DENSE: bool = true; + fn order(&self) -> usize { + 5 + } + + fn stages(&self) -> usize { + STAGES + } + + fn adaptive(&self) -> bool { + true + } + + fn dense(&self) -> bool { + true + } fn step(&self, ode: &ODE, h: T) -> StepResult { - let mut k: Vec> = vec![SVector::::zeros(); Self::STAGES]; + let mut k: Vec> = vec![SVector::::zeros(); STAGES]; let mut next_y = ode.y; let mut err = SVector::::zeros(); let mut rcont5 = SVector::::zeros(); // Do the first of the summations k[0] = (ode.f)(ode.t, ode.y, &ode.params); next_y += k[0] * T::lift(Self::B[0]) * h; - err += k[0] * T::lift(Self::B[0] - Self::B[Self::STAGES]) * h; + err += k[0] * T::lift(Self::B[0] - Self::B[STAGES]) * h; let rcont1 = ode.y; rcont5 += k[0] * h * T::lift(Self::D[0]); // Then the rest - for i in 1..Self::STAGES { + for i in 1..STAGES { // Compute the ks let mut y_term = SVector::::zeros(); for (j, item) in k.iter().enumerate().take(i) { @@ -201,12 +216,12 @@ where // Use that and bis to calculate the y and error terms next_y += k[i] * h * T::lift(Self::B[i]); - err += k[i] * T::lift(Self::B[i] - Self::B[i + Self::STAGES]) * h; + err += k[i] * T::lift(Self::B[i] - Self::B[i + STAGES]) * h; rcont5 += k[i] * h * T::lift(Self::D[i]); } let rcont2 = next_y - ode.y; let rcont3 = k[0] * h - rcont2; - let rcont4 = rcont2 - k[Self::STAGES - 1] * h - rcont3; + let rcont4 = rcont2 - k[STAGES - 1] * h - rcont3; let rcont = vec![rcont1, rcont2, rcont3, rcont4, rcont5]; // Just take the real component of the error let mut error_norm_sq = 0.0; diff --git a/src/integrator/method.rs b/src/integrator/method.rs new file mode 100644 index 0000000..c9b03e4 --- /dev/null +++ b/src/integrator/method.rs @@ -0,0 +1,444 @@ +//! An integrator chosen at runtime. +//! +//! Which integrator solves a problem is normally a type parameter, fixed when +//! the code is compiled. [`Method`] is an enum over the integrators this crate +//! provides that itself implements [`Integrator`], so a value can stand where a +//! type would otherwise have to: an accuracy preset can carry a method as well +//! as a tolerance, a method can be named in a configuration file through +//! [`MethodName`], and a `Vec>` can hold several for comparison. +//! +//! Dispatch is a `match` on the variant, which the compiler resolves per call +//! site. There is no trait object and no scalar erasure, so autodiff keeps +//! working through a `Method` exactly as through the integrator it wraps. + +use nalgebra::SVector; +use serde::{Deserialize, Serialize}; + +use super::bs3::BS3; +use super::dormand_prince::DormandPrince45; +use super::rosenbrock::Rosenbrock23; +use super::vern7::Vern7; +use super::{Integrator, StepResult, Tolerances}; +use crate::ode::ODE; +use crate::scalar::Scalar; + +/// An integrator chosen at runtime. +/// +/// Every variant wraps a configured integrator, so tolerances are set on a +/// `Method` the same way they are set on the integrator itself. +/// +/// # Example +/// +/// ```rust +/// use ordinary_diffeq::prelude::*; +/// use nalgebra::Vector1; +/// +/// fn derivative(_t: f64, y: Vector1, _p: &()) -> Vector1 { +/// -y +/// } +/// +/// // The method is a value here, so it can come from a configuration file, a +/// // command line flag, or an accuracy preset. +/// let method: Method<1> = MethodName::Vern7.instantiate().a_tol(1e-12).r_tol(1e-12); +/// +/// let ode = ODE::new(&derivative, 0.0, 2.0, Vector1::new(1.0), ()); +/// let mut problem = Problem::new(ode, method, PIController::default()); +/// let solution = problem.solve(); +/// +/// let y = solution.interpolate(1.3)[0]; +/// assert!((y - (-1.3f64).exp()).abs() < 1e-10); +/// ``` +#[non_exhaustive] +#[derive(Debug, Clone, Copy)] +pub enum Method { + /// Bogacki-Shampine 3(2). + BS3(BS3), + /// Dormand-Prince 5(4). + DormandPrince45(DormandPrince45), + /// Verner's "most efficient" 7(6). + Vern7(Vern7), + /// Rosenbrock-W 2(3), for stiff problems. + Rosenbrock23(Rosenbrock23), +} + +/// Forwards a `&self` call to whichever integrator the method holds. +macro_rules! dispatch { + ($method:expr, $inner:ident => $call:expr) => { + match $method { + Method::BS3($inner) => $call, + Method::DormandPrince45($inner) => $call, + Method::Vern7($inner) => $call, + Method::Rosenbrock23($inner) => $call, + } + }; +} + +/// Forwards a consuming call and puts the result back in the same variant. +macro_rules! dispatch_rewrap { + ($method:expr, $inner:ident => $call:expr) => { + match $method { + Method::BS3($inner) => Method::BS3($call), + Method::DormandPrince45($inner) => Method::DormandPrince45($call), + Method::Vern7($inner) => Method::Vern7($call), + Method::Rosenbrock23($inner) => Method::Rosenbrock23($call), + } + }; +} + +impl Method { + /// Which method this is, without its configuration. + /// + /// The counterpart to [`MethodName::instantiate`], so a run can record the + /// method it used in a form that selects the same one again. + pub fn name(&self) -> MethodName { + match self { + Method::BS3(_) => MethodName::BS3, + Method::DormandPrince45(_) => MethodName::DormandPrince45, + Method::Vern7(_) => MethodName::Vern7, + Method::Rosenbrock23(_) => MethodName::Rosenbrock23, + } + } +} + +impl Integrator for Method { + fn order(&self) -> usize { + dispatch!(self, m => m.order()) + } + + fn stages(&self) -> usize { + dispatch!(self, m => m.stages()) + } + + fn adaptive(&self) -> bool { + dispatch!(self, m => m.adaptive()) + } + + fn dense(&self) -> bool { + dispatch!(self, m => m.dense()) + } + + fn main_stages(&self) -> usize { + dispatch!(self, m => m.main_stages()) + } + + fn extra_stages(&self) -> usize { + dispatch!(self, m => m.extra_stages()) + } + + fn total_dense_stages(&self) -> usize { + dispatch!(self, m => m.total_dense_stages()) + } + + fn step(&self, ode: &ODE, h: T) -> StepResult { + dispatch!(self, m => m.step(ode, h)) + } + + fn interpolate( + &self, + t_start: T, + t_end: T, + dense: &[SVector], + t: T, + ) -> SVector { + dispatch!(self, m => m.interpolate(t_start, t_end, dense, t)) + } + + fn compute_extra_stages( + &self, + ode: &ODE, + t_start: T, + y_start: SVector, + h: T, + main_stages: &[SVector], + ) -> Vec> { + dispatch!(self, m => m.compute_extra_stages(ode, t_start, y_start, h, main_stages)) + } +} + +impl Tolerances for Method { + fn absolute_tolerance(&self) -> SVector { + dispatch!(self, m => m.absolute_tolerance()) + } + + fn relative_tolerance(&self) -> f64 { + dispatch!(self, m => m.relative_tolerance()) + } + + fn a_tol(self, a_tol: f64) -> Self { + dispatch_rewrap!(self, m => m.a_tol(a_tol)) + } + + fn a_tol_full(self, a_tol: SVector) -> Self { + dispatch_rewrap!(self, m => m.a_tol_full(a_tol)) + } + + fn r_tol(self, r_tol: f64) -> Self { + dispatch_rewrap!(self, m => m.r_tol(r_tol)) + } +} + +impl From> for Method { + fn from(integrator: BS3) -> Self { + Method::BS3(integrator) + } +} + +impl From> for Method { + fn from(integrator: DormandPrince45) -> Self { + Method::DormandPrince45(integrator) + } +} + +impl From> for Method { + fn from(integrator: Vern7) -> Self { + Method::Vern7(integrator) + } +} + +impl From> for Method { + fn from(integrator: Rosenbrock23) -> Self { + Method::Rosenbrock23(integrator) + } +} + +/// The name of an integrator, without its configuration. +/// +/// Serializes as the variant name, so `integrator = "Vern7"` in a configuration +/// file deserializes to [`MethodName::Vern7`] and +/// [`instantiate`](MethodName::instantiate) turns it into a [`Method`] at +/// whatever dimension the problem needs. +/// +/// # Example +/// +/// ```rust +/// use ordinary_diffeq::prelude::*; +/// +/// let name: MethodName = serde_json::from_str(r#""Vern7""#).unwrap(); +/// let method: Method<6> = name.instantiate(); +/// assert_eq!(method.name(), MethodName::Vern7); +/// ``` +#[non_exhaustive] +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)] +pub enum MethodName { + /// Bogacki-Shampine 3(2). + BS3, + /// Dormand-Prince 5(4). + DormandPrince45, + /// Verner's "most efficient" 7(6). + Vern7, + /// Rosenbrock-W 2(3), for stiff problems. + Rosenbrock23, +} + +impl MethodName { + /// Builds the named integrator, at its default tolerances. + pub fn instantiate(&self) -> Method { + match self { + MethodName::BS3 => Method::BS3(BS3::new()), + MethodName::DormandPrince45 => Method::DormandPrince45(DormandPrince45::new()), + MethodName::Vern7 => Method::Vern7(Vern7::new()), + MethodName::Rosenbrock23 => Method::Rosenbrock23(Rosenbrock23::new()), + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::controller::PIController; + use crate::problem::Problem; + use nalgebra::Vector2; + + /// The harmonic oscillator, which every method solves and which has enough + /// curvature that steps land in different places at different tolerances. + fn derivative(_t: f64, y: Vector2, _p: &()) -> Vector2 { + Vector2::new(y[1], -y[0]) + } + + /// Solves to `t = 4` and reports the step times, the states, and the value + /// interpolated at ten points that fall strictly inside steps. + fn solve + Copy>(integrator: S) -> (Vec, Vec>, Vec) { + let ode = crate::ode::ODE::new(&derivative, 0.0, 4.0, Vector2::new(1.0, 0.0), ()); + let mut problem = Problem::new(ode, integrator, PIController::default()); + let solution = problem.solve(); + + let interpolated = (1..=10) + .map(|i| solution.interpolate(i as f64 * 0.37)[0]) + .collect(); + ( + solution.times.clone(), + solution.states.clone(), + interpolated, + ) + } + + /// Wrapping an integrator in `Method` must not change a single bit of the + /// answer, interpolation included — which for Vern7 exercises the lazy + /// extra stages, the path `extra_stages` governs. + #[test] + fn test_method_matches_wrapped_integrator() { + macro_rules! assert_identical { + ($integrator:expr, $variant:ident) => {{ + let integrator = $integrator; + let (direct_t, direct_y, direct_interp) = solve(integrator); + let (wrapped_t, wrapped_y, wrapped_interp) = solve(Method::$variant(integrator)); + + assert_eq!( + direct_t, + wrapped_t, + "{} step times differ", + stringify!($variant) + ); + assert_eq!( + direct_y, + wrapped_y, + "{} states differ", + stringify!($variant) + ); + assert_eq!( + direct_interp, + wrapped_interp, + "{} interpolated values differ", + stringify!($variant) + ); + // And the step count is nontrivial, so the comparison means something. + assert!( + direct_t.len() > 5, + "{} took too few steps", + stringify!($variant) + ); + }}; + } + + assert_identical!(BS3::<2>::new().a_tol(1e-8).r_tol(1e-8), BS3); + assert_identical!( + DormandPrince45::<2>::new().a_tol(1e-10).r_tol(1e-10), + DormandPrince45 + ); + assert_identical!(Vern7::<2>::new().a_tol(1e-12).r_tol(1e-12), Vern7); + assert_identical!( + Rosenbrock23::<2>::new().a_tol(1e-6).r_tol(1e-6), + Rosenbrock23 + ); + } + + /// The structural facts a `Method` reports are the wrapped integrator's. + #[test] + fn test_method_reports_wrapped_structure() { + let vern7 = Vern7::<2>::new(); + let wrapped = Method::Vern7(vern7); + assert_eq!(wrapped.order(), vern7.order()); + assert_eq!(wrapped.stages(), vern7.stages()); + assert_eq!(wrapped.main_stages(), vern7.main_stages()); + assert_eq!(wrapped.extra_stages(), vern7.extra_stages()); + assert_eq!(wrapped.total_dense_stages(), vern7.total_dense_stages()); + assert!(wrapped.adaptive() && wrapped.dense()); + + // Dormand-Prince has no extra stages, and the default applies through + // the enum as it does directly. + let dp45 = Method::DormandPrince45(DormandPrince45::<2>::new()); + assert_eq!(dp45.extra_stages(), 0); + assert_eq!(dp45.total_dense_stages(), dp45.main_stages()); + } + + /// A name round-trips through serde and instantiates the variant it names. + #[test] + fn test_method_name_round_trip() { + for (name, json) in [ + (MethodName::BS3, r#""BS3""#), + (MethodName::DormandPrince45, r#""DormandPrince45""#), + (MethodName::Vern7, r#""Vern7""#), + (MethodName::Rosenbrock23, r#""Rosenbrock23""#), + ] { + assert_eq!(serde_json::to_string(&name).unwrap(), json); + + let parsed: MethodName = serde_json::from_str(json).unwrap(); + assert_eq!(parsed, name); + + let method: Method<3> = parsed.instantiate(); + assert_eq!(method.name(), name); + } + } + + /// Tolerance setters forward through the enum and keep the variant. + #[test] + fn test_tolerances_through_the_enum() { + let method: Method<2> = MethodName::BS3.instantiate().a_tol(1e-11).r_tol(1e-9); + assert_eq!(method.name(), MethodName::BS3); + assert_eq!( + method.absolute_tolerance(), + SVector::::repeat(1e-11) + ); + assert_eq!(method.relative_tolerance(), 1e-9); + + let per = SVector::::new(1e-6, 1e-12); + let method = method.a_tol_full(per); + assert_eq!(method.absolute_tolerance(), per); + } + + /// The heterogeneous collection the type parameter used to rule out. + #[test] + fn test_collection_of_methods() { + let methods: Vec> = [ + MethodName::BS3, + MethodName::DormandPrince45, + MethodName::Vern7, + MethodName::Rosenbrock23, + ] + .iter() + .map(|name| name.instantiate().a_tol(1e-8).r_tol(1e-8)) + .collect(); + + for method in methods { + let (times, states, _) = solve(method); + // cos(4) at the end of the interval, to whatever each method manages. + let final_y = states.last().unwrap()[0]; + assert!( + (final_y - 4.0f64.cos()).abs() < 1e-3, + "{:?} ended at {final_y} after {} steps", + method.name(), + times.len() + ); + } + } + + /// The module header claims a `Method` differentiates exactly as the + /// integrator it wraps does, which is only true if nothing erases the + /// scalar. Vern7 is the demanding case: interpolating inside a step calls + /// the right-hand side again, in dual arithmetic, after the solve. + #[cfg(feature = "autodiff")] + #[test] + fn test_method_preserves_autodiff() { + use num_dual::Dual64; + + fn decay(_t: T, y: nalgebra::Vector1) -> nalgebra::Vector1 { + -y + } + + let y0 = Dual64::from_re(1.0).derivative(); + let ode = crate::ode::ODE::new( + &|t, y, _p: &()| decay(t, y), + Dual64::lift(0.0), + Dual64::lift(2.0), + nalgebra::Vector1::new(y0), + (), + ); + let method: Method<1> = MethodName::Vern7.instantiate().a_tol(1e-12).r_tol(1e-12); + let mut problem = Problem::new(ode, method, PIController::default()); + let solution = problem.solve(); + + let y = solution.interpolate(Dual64::lift(1.3))[0]; + + // y(t) = y0 e^-t, so both the value and dy/dy0 are e^-1.3. + assert!((y.re - (-1.3f64).exp()).abs() < 1e-10); + assert!((y.eps - (-1.3f64).exp()).abs() < 1e-10); + } + + /// `From` gets an already-configured integrator into the enum. + #[test] + fn test_from_integrator() { + let method: Method<2> = Vern7::<2>::new().a_tol(1e-9).into(); + assert_eq!(method.name(), MethodName::Vern7); + assert_eq!(method.absolute_tolerance(), SVector::::repeat(1e-9)); + } +} diff --git a/src/integrator/mod.rs b/src/integrator/mod.rs index 3cfba63..485ba37 100644 --- a/src/integrator/mod.rs +++ b/src/integrator/mod.rs @@ -5,6 +5,7 @@ use crate::scalar::Scalar; pub mod bs3; pub mod dormand_prince; +pub mod method; pub mod rosenbrock; pub mod vern7; @@ -23,27 +24,48 @@ pub struct StepResult { /// [`Integrator::interpolate`], or `None` from a method without dense /// output. /// - /// For a method with lazy dense output (`EXTRA_STAGES > 0`) the layout is - /// `[y_start, y_end, k_1, ..., k_MAIN_STAGES]`, and the extra stages are + /// For a method with lazy dense output (`extra_stages() > 0`) the layout is + /// `[y_start, y_end, k_1, ..., k_main_stages]`, and the extra stages are /// appended by [`Integrator::compute_extra_stages`] on first use. pub dense: Option>>, } /// Integrator Trait +/// +/// The structural facts about a method are `&self` methods rather than +/// associated constants, so that a value can answer them. That is what lets +/// [`Method`](method::Method) wrap the integrators in an enum and dispatch to +/// whichever one it holds, and so lets an integrator be chosen at runtime. None +/// of these appears in type position, and for a concrete integrator each call +/// inlines to the constant it returns. pub trait Integrator { - const ORDER: usize; - const STAGES: usize; - const ADAPTIVE: bool; - const DENSE: bool; + /// The order of the method's solution. + fn order(&self) -> usize; - /// Number of main stages stored in dense output (default: same as STAGES) - const MAIN_STAGES: usize = Self::STAGES; + /// The number of stages the method evaluates per step. + fn stages(&self) -> usize; + + /// Whether the method estimates its local error, and so can be stepped + /// adaptively. + fn adaptive(&self) -> bool; + + /// Whether the method produces dense output. + fn dense(&self) -> bool; + + /// Number of main stages stored in dense output (default: same as `stages`) + fn main_stages(&self) -> usize { + self.stages() + } /// Number of extra stages for full-order dense output (default: 0, no extra stages) - const EXTRA_STAGES: usize = 0; + fn extra_stages(&self) -> usize { + 0 + } /// Total stages when full dense output is computed - const TOTAL_DENSE_STAGES: usize = Self::MAIN_STAGES + Self::EXTRA_STAGES; + fn total_dense_stages(&self) -> usize { + self.main_stages() + self.extra_stages() + } /// Advances `ode` by `h` from its current state, without mutating it. fn step(&self, ode: &ODE, h: T) -> StepResult; diff --git a/src/integrator/rosenbrock.rs b/src/integrator/rosenbrock.rs index c692fc7..67a8910 100644 --- a/src/integrator/rosenbrock.rs +++ b/src/integrator/rosenbrock.rs @@ -229,10 +229,21 @@ impl Default for Rosenbrock23 { } impl Integrator for Rosenbrock23 { - const ORDER: usize = 2; - const STAGES: usize = 2; - const ADAPTIVE: bool = true; - const DENSE: bool = true; + fn order(&self) -> usize { + 2 + } + + fn stages(&self) -> usize { + 2 + } + + fn adaptive(&self) -> bool { + true + } + + fn dense(&self) -> bool { + true + } fn step(&self, ode: &ODE, h: T) -> StepResult { let t = ode.t; diff --git a/src/integrator/vern7.rs b/src/integrator/vern7.rs index bc377e6..11451c0 100644 --- a/src/integrator/vern7.rs +++ b/src/integrator/vern7.rs @@ -4,6 +4,13 @@ use super::super::ode::ODE; use super::{Integrator, StepResult, Tolerances}; use crate::scalar::Scalar; +/// The number of main stages Vern7 evaluates per step. +const STAGES: usize = 10; + +/// The number of further stages Vern7 computes, lazily, for full-order dense +/// output. +const EXTRA_STAGES: usize = 6; + /// Verner 7 integrator trait for tableau coefficients pub trait Vern7Integrator<'a> { const A: &'a [f64]; // Lower triangular A matrix (flattened) @@ -289,24 +296,40 @@ impl<'a, const D: usize> Integrator for Vern7 where Vern7: Vern7Integrator<'a> + Vern7ExtraStages<'a>, { - const ORDER: usize = 7; - const STAGES: usize = 10; - const ADAPTIVE: bool = true; - const DENSE: bool = true; + fn order(&self) -> usize { + 7 + } + + fn stages(&self) -> usize { + STAGES + } + + fn adaptive(&self) -> bool { + true + } + + fn dense(&self) -> bool { + true + } // Lazy dense output configuration - const MAIN_STAGES: usize = 10; - const EXTRA_STAGES: usize = 6; + fn main_stages(&self) -> usize { + STAGES + } + + fn extra_stages(&self) -> usize { + EXTRA_STAGES + } fn step(&self, ode: &ODE, h: T) -> StepResult { // Allocate storage for the 10 stages - let mut k: Vec> = vec![SVector::::zeros(); Self::STAGES]; + let mut k: Vec> = vec![SVector::::zeros(); STAGES]; // Stage 1: k[0] = f(t, y) k[0] = (ode.f)(ode.t, ode.y, &ode.params); // Compute remaining stages using the A matrix - for i in 1..Self::STAGES { + for i in 1..STAGES { let mut y_temp = ode.y; // A matrix is stored in lower triangular form, row by row // Row i has i elements (0-indexed), starting at position i*(i-1)/2 @@ -469,7 +492,7 @@ where k9 * b9_theta) * h; // If extra stages are available, add their contribution for full 7th order accuracy - if dense.len() >= 2 + Self::TOTAL_DENSE_STAGES { + if dense.len() >= 2 + STAGES + EXTRA_STAGES { // Extra stages are at indices 12-17 let k11 = &dense[12]; let k12 = &dense[13]; @@ -563,7 +586,7 @@ where let k8 = &main_stages[7]; let k9 = &main_stages[8]; - let mut extra_stages = Vec::with_capacity(Self::EXTRA_STAGES); + let mut extra_stages = Vec::with_capacity(EXTRA_STAGES); // Stage 11: uses k1, k4-k9 (7 coefficients) let mut y11 = y_start; diff --git a/src/lib.rs b/src/lib.rs index b8b065a..36e7b3e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -16,6 +16,7 @@ pub mod prelude { pub use super::controller::{PIController, PIDController}; pub use super::integrator::bs3::BS3; pub use super::integrator::dormand_prince::DormandPrince45; + pub use super::integrator::method::{Method, MethodName}; pub use super::integrator::rosenbrock::Rosenbrock23; pub use super::integrator::vern7::Vern7; pub use super::ode::ODE; diff --git a/src/problem.rs b/src/problem.rs index 9d58f43..6d903a1 100644 --- a/src/problem.rs +++ b/src/problem.rs @@ -52,7 +52,7 @@ where self.controller .set_next_step_guess(TryStep::NotYetAccepted(self.ode.t_end - self.ode.t)); } - let (mut new_y, mut curr_step, mut dense_option) = if S::ADAPTIVE { + let (mut new_y, mut curr_step, mut dense_option) = if self.integrator.adaptive() { // First, we try stepping with the "next step guess" to get the error let mut result = self .integrator @@ -175,13 +175,13 @@ where // Check if we need to compute extra stages for lazy dense output let dense_cell = &self.dense[end_index - 1]; - if S::EXTRA_STAGES > 0 { + if self.integrator.extra_stages() > 0 { let needs_extra = { let borrowed = dense_cell.borrow(); // Dense array format: [y0, y1, k1, k2, ..., k_main] - // If we have main stages only: 2 + MAIN_STAGES elements - // If we have all stages: 2 + MAIN_STAGES + EXTRA_STAGES elements - borrowed.len() < 2 + S::TOTAL_DENSE_STAGES + // If we have main stages only: 2 + main_stages() elements + // If we have all stages: 2 + total_dense_stages() elements + borrowed.len() < 2 + self.integrator.total_dense_stages() }; if needs_extra { @@ -189,7 +189,7 @@ where let mut dense = dense_cell.borrow_mut(); // Extract main stages (skip y0 and y1 at indices 0 and 1) - let main_stages = &dense[2..2 + S::MAIN_STAGES]; + let main_stages = &dense[2..2 + self.integrator.main_stages()]; // Compute extra stages lazily let extra_stages = self.integrator.compute_extra_stages(