Release 0.5.0

Runtime integrator selection. The structural facts on Integrator move from
associated constants to &self methods, which is a breaking change for anyone
implementing the trait or reading ORDER, STAGES, ADAPTIVE or DENSE. In exchange
Method dispatches to any of the four integrators as a value, and MethodName
names one in a configuration file.

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "ordinary-diffeq"
-version = "0.3.1"
+version = "0.5.0"
 edition = "2021"
 authors = ["Connor Johnstone"]
 description = "A library for solving differential equations based on the DifferentialEquations.jl julia library."
@@ -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"

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<f64>, _p: &()) -> Vector1<f64> {
+    -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<Method<D>>`
+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,

roadmap/README.md

@@ -58,6 +58,25 @@ Each feature below links to a detailed implementation plan in the `features/` di
   - **Dependencies**: None
   - **Effort**: Small
 
+### API Ergonomics
+
+- [x] **[Tolerances Trait](features/40-tolerances-trait.md)** ✅ COMPLETED
+  - Move `a_tol` / `a_tol_full` / `r_tol` off the inherent impls into a trait, and
+    add getters, so a generic `S` can set and read them
+  - Raised by `astrokit`, which declared its own copy as a workaround
+  - **Dependencies**: None
+  - **Effort**: Small
+  - **Status**: All four integrators implement `Tolerances`; no call site changed
+
+- [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
 
 - [ ] **[Discrete Callbacks](features/05-discrete-callbacks.md)**
@@ -336,17 +355,19 @@ Each algorithm implementation should include:
 
 ## Progress Tracking
 
-Total Features: 38
-- Tier 1: 8 features (4/8 complete) ✅
-- Tier 2: 12 features (0/12 complete)
-- Tier 3: 18 features (0/18 complete)
+Total Features: 41
+- Tier 1: 10 features (6/10 complete) ✅
+- Tier 2: 11 features (0/11 complete)
+- Tier 3: 20 features (0/20 complete)
 
-**Overall Progress: 10.5% (4/38 features complete)**
+**Overall Progress: 14.6% (6/41 features complete)**
 
 ### Completed Features
 1. ✅ BS3 (Bogacki-Shampine 3/2) - Tier 1 (2025-10-23)
 2. ✅ Vern7 (Verner 7th order) - Tier 1 (2025-10-24)
 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: 2025-10-24
+Last updated: 2026-09-30

roadmap/features/40-tolerances-trait.md

@@ -0,0 +1,129 @@
+# Feature: Tolerances Trait
+
+**Status**: completed (2026-09-30)
+
+## Overview
+
+Error tolerances are set through inherent builder methods on each integrator
+rather than through a trait:
+
+```rust
+impl<const D: usize> Vern7<D> where Vern7<D>: Integrator<D> {
+    pub fn a_tol(mut self, a_tol: f64) -> Self { .. }
+    pub fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self { .. }
+    pub fn r_tol(mut self, r_tol: f64) -> Self { .. }
+}
+```
+
+`DormandPrince45`, `BS3` and `Rosenbrock23` each have the same four methods,
+written out separately. There is no trait bound that says "this integrator has
+tolerances", so a downstream crate generic over `S: Integrator<D>` cannot set
+them, and cannot read back what they were.
+
+## Why This Feature Matters
+
+A library wrapping this one wants to offer accuracy presets:
+
+```rust
+Propagation::new(..).accuracy(Accuracy::Medium)   // sets a_tol and r_tol
+Propagation::new(..).a_tol(1e-10)
+```
+
+That is impossible over a generic `S`. `astrokit` worked around it by declaring
+its own trait and implementing it for three of the four integrators:
+
+```rust
+pub trait Tolerances<const D: usize>: Integrator<D> {
+    fn with_a_tol(self, a_tol: f64) -> Self;
+    fn with_a_tol_per_component(self, a_tol: SVector<f64, D>) -> Self;
+    fn with_r_tol(self, r_tol: f64) -> Self;
+}
+impl<const D: usize> Tolerances<D> for Vern7<D> { .. }
+```
+
+That works, but the trait belongs here. Downstream it has three problems: any
+integrator added to this crate is invisible to it until the downstream crate is
+updated, a second dependent would declare the same trait again, and the orphan
+rule means only crates that own the trait can extend it.
+
+Reading tolerances back matters separately. A run record that says which
+integrator ran but not at what tolerance does not describe the run.
+
+## Dependencies
+
+None.
+
+## Implementation Approach
+
+A separate trait rather than more methods on `Integrator`. `Integrator` already
+carries `const ADAPTIVE: bool` because not every method is adaptive, and a
+fixed-step method has no tolerances to set. Splitting them keeps that honest and
+leaves room for a fixed-step integrator that implements `Integrator` and not
+this.
+
+```rust
+/// An adaptive integrator's error tolerances.
+pub trait Tolerances<const D: usize>: Integrator<D> {
+    /// The absolute tolerance, per component.
+    fn absolute_tolerance(&self) -> SVector<f64, D>;
+
+    /// The relative tolerance.
+    fn relative_tolerance(&self) -> f64;
+
+    /// Sets the absolute tolerance for every component.
+    fn a_tol(self, a_tol: f64) -> Self;
+
+    /// Sets the absolute tolerance per component.
+    fn a_tol_full(self, a_tol: SVector<f64, D>) -> Self;
+
+    /// Sets the relative tolerance.
+    fn r_tol(self, r_tol: f64) -> Self;
+}
+```
+
+The setters keep the names they already have, so **no existing call site
+changes**. `Vern7::new().a_tol(1e-12).r_tol(1e-12)` continues to compile,
+resolving through the trait instead of the inherent impl, because `Integrator` is
+already in the prelude and `Tolerances` would be too.
+
+The getters need different names because `a_tol` is taken by the setter. Naming
+them `absolute_tolerance`/`relative_tolerance` avoids the collision and reads
+better than `a_tol_value`.
+
+### Steps
+
+1. Declare `Tolerances<D>` in `src/integrator/mod.rs`.
+2. For each of `Vern7`, `DormandPrince45`, `BS3`, `Rosenbrock23`: move the three
+   setter bodies from the inherent impl into an `impl Tolerances<D>`, and add the
+   two getters returning the stored fields. Keep `new()` inherent.
+3. Export `Tolerances` from the prelude.
+4. Delete `astrokit`'s local copy and import this one.
+
+### What not to do
+
+Do not add the tolerance methods to `Integrator` itself. That would oblige a
+fixed-step method to implement setters that mean nothing, and `const ADAPTIVE`
+exists precisely to admit such methods.
+
+Do not add `Default` or `Copy` as supertraits of `Integrator` while doing this,
+tempting as it is for downstream bounds. All four current integrators are both,
+but a multistep method holding Nordsieck history would reasonably be neither, and
+a downstream crate can bound on them explicitly where it needs them.
+
+## Testing Requirements
+
+- Every integrator round-trips: set a tolerance, read it back.
+- `a_tol` sets all components; `a_tol_full` sets them independently.
+- A generic function over `S: Tolerances<D>` compiles and changes step counts as
+  the tolerance changes, for each integrator.
+- Existing tests continue to pass unchanged, which is the point: the setter names
+  do not move.
+
+## Success Criteria
+
+- A downstream crate generic over `S: Tolerances<D>` can set and read tolerances
+  for every adaptive integrator in this crate.
+- No call site in this crate or its tests changes.
+- `astrokit` deletes its local `Tolerances`.
+
+**Effort**: Small

roadmap/features/41-runtime-integrator-selection.md

@@ -0,0 +1,178 @@
+# Feature: Runtime Integrator Selection
+
+**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
+
+Which integrator runs a problem is fixed at compile time, because `Integrator<D>`
+keeps its structural facts in associated constants:
+
+```rust
+pub trait Integrator<const D: usize> {
+    const ORDER: usize;
+    const STAGES: usize;
+    const ADAPTIVE: bool;
+    const DENSE: bool;
+    const MAIN_STAGES: usize = Self::STAGES;
+    const EXTRA_STAGES: usize = 0;
+    const TOTAL_DENSE_STAGES: usize = Self::MAIN_STAGES + Self::EXTRA_STAGES;
+}
+```
+
+`Problem` is generic over `S`, and `Problem::solve` reads `S::ADAPTIVE` to choose
+its stepping branch and `S::EXTRA_STAGES` / `S::TOTAL_DENSE_STAGES` to decide
+whether lazy stages still need computing. A value cannot answer those questions —
+only a type can.
+
+So the obvious workaround does not work. An enum wrapping the four integrators
+cannot implement `Integrator<D>`, because it would have to give one
+`EXTRA_STAGES` for all of them, and Vern7 has six where Dormand-Prince has none.
+
+`Box<dyn Integrator<D>>` is not available either, and will not become available:
+`step<T: Scalar, P>` and `interpolate<T: Scalar>` are generic methods, so the
+trait is not object-safe. Erasing the scalar would defeat the genericity that
+makes autodiff work. The enum is the practical route, not `dyn`.
+
+## Why This Feature Matters
+
+This surfaced building `astrokit` on top of this crate. Three things a dependent
+wants and cannot have:
+
+1. **Accuracy presets that include a method.** `Accuracy::High` should mean "tight
+   tolerances *and* a high-order method", since the two belong together — a 7th
+   order method at 1e-6 is waste, and a 3rd order method at 1e-12 is hopeless. As
+   it stands a preset can only set tolerances, and the method has to be named
+   separately as a type parameter.
+
+2. **An integrator named in a configuration file.** `astrokit` is built for a
+   study where every run parameter must be declared in a committed file and every
+   propagation must record enough to recreate itself. `integrator: "Vern7"`
+   deserialized from that file cannot select anything, because the choice is a
+   type. The integrator ends up recorded in the provenance but not settable from
+   it, which is backwards.
+
+3. **A heterogeneous collection.** Benchmarking or comparing methods over the same
+   problem means either a macro over method types or duplicated code per method,
+   where a `Vec<Method>` would do.
+
+None of these is fatal. `astrokit` has settled on Vern7 as a type-level default,
+which is the right method for its accuracy range anyway. But the constraint is
+structural rather than deliberate, and it is worth fixing before more dependents
+grow around it.
+
+## Dependencies
+
+None, though it interacts with [Tolerances Trait](40-tolerances-trait.md): an
+enum wrapper needs to forward tolerance setters too, which is much easier once
+those are trait methods.
+
+## Implementation Approach
+
+### Step 1: constants become methods
+
+The structural facts move from associated constants to methods taking `&self`:
+
+```rust
+pub trait Integrator<const D: usize> {
+    fn order(&self) -> usize;
+    fn stages(&self) -> usize;
+    fn adaptive(&self) -> bool;
+    fn dense(&self) -> bool;
+
+    fn main_stages(&self) -> usize { self.stages() }
+    fn extra_stages(&self) -> usize { 0 }
+    fn total_dense_stages(&self) -> usize { self.main_stages() + self.extra_stages() }
+
+    // step, interpolate, compute_extra_stages unchanged
+}
+```
+
+**This costs nothing measurable.** Every use of these constants today is already a
+runtime comparison, not a compile-time size:
+
+- `if S::ADAPTIVE` in `Problem::solve` becomes `if self.integrator.adaptive()`,
+  which monomorphizes to the same constant and folds the same way for a concrete
+  integrator.
+- `borrowed.len() < 2 + S::TOTAL_DENSE_STAGES` compares against a `Vec` length.
+- `&dense[2..2 + S::MAIN_STAGES]` slices a `Vec`.
+
+None of them appears in type position, so nothing needs the values to be consts.
+For a concrete `S` the method calls inline to the same constants they are now.
+
+This is a breaking change to anyone implementing `Integrator` by hand — but a
+mechanical one, and this crate has no external implementors.
+
+### Step 2: an enum in this crate
+
+With step 1 done, an enum can dispatch:
+
+```rust
+/// An integrator chosen at runtime.
+#[non_exhaustive]
+#[derive(Debug, Clone, Copy)]
+pub enum Method<const D: usize> {
+    BS3(BS3<D>),
+    DormandPrince45(DormandPrince45<D>),
+    Vern7(Vern7<D>),
+    Rosenbrock23(Rosenbrock23<D>),
+}
+
+impl<const D: usize> Integrator<D> for Method<D> {
+    fn extra_stages(&self) -> usize {
+        match self {
+            Method::Vern7(m) => m.extra_stages(),
+            Method::BS3(m) => m.extra_stages(),
+            ..
+        }
+    }
+    // and so on, including step and interpolate
+}
+```
+
+Each method body is a `match` forwarding to the wrapped integrator. Verbose, and
+a good candidate for a macro over the variant list.
+
+Worth deriving `Serialize`/`Deserialize` on a companion `MethodName` enum so a
+method can be named in a configuration file, and adding
+`MethodName::instantiate(&self) -> Method<D>`. That is what makes point 2 above
+work.
+
+### Step 3: presets
+
+Once a method is a value, an accuracy preset can carry one, pairing a tolerance
+with a method whose order suits it. That belongs in whichever crate owns the
+notion of a preset — probably this one, given
+[Default Algorithm Selection](13-default-algorithm-selection.md) is already on
+the roadmap and needs the same machinery.
+
+### The cost of not doing step 1
+
+If the constants must stay, the fallback is a parallel object-ish layer: a second
+trait with `&self` methods, blanket-implemented for everything that implements
+`Integrator`, and `Problem` generic over that instead. It works, it duplicates the
+trait surface, and every implementor carries both. Step 1 is cheaper and cleaner.
+
+## Testing Requirements
+
+- Every existing test passes after step 1, with method calls substituted for
+  constant reads.
+- A benchmark confirms no regression from the constants becoming methods — the
+  claim is that they inline, and that should be measured rather than assumed, on
+  the existing `simple_1d` and `orbit` benches.
+- `Method` gives bit-identical results to the integrator it wraps, for each
+  variant, on a problem with dense output interpolation (which exercises the
+  lazy-stage path where `extra_stages` matters).
+- A method name round-trips through serde and instantiates the right variant.
+
+## Success Criteria
+
+- `Problem` can be constructed with an integrator chosen at runtime.
+- No measurable performance change on the existing criterion benches.
+- A dependent crate can offer an accuracy preset that selects both a tolerance and
+  a method, and can read a method name out of a configuration file.
+
+**Effort**: Medium — step 1 is small and mechanical, step 2 is repetitive

src/integrator/bs3.rs

@@ -1,9 +1,12 @@
 use nalgebra::SVector;
 
 use super::super::ode::ODE;
-use super::{Integrator, StepResult};
+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];
@@ -80,21 +83,31 @@ where
             r_tol: 1e-8,
         }
     }
+}
 
-    /// Set absolute tolerance (same value for all components)
-    pub fn a_tol(mut self, a_tol: f64) -> Self {
+impl<const D: usize> Tolerances<D> for BS3<D>
+where
+    BS3<D>: Integrator<D>,
+{
+    fn absolute_tolerance(&self) -> SVector<f64, D> {
+        self.a_tol
+    }
+
+    fn relative_tolerance(&self) -> f64 {
+        self.r_tol
+    }
+
+    fn a_tol(mut self, a_tol: f64) -> Self {
         self.a_tol = SVector::<f64, D>::from_element(a_tol);
         self
     }
 
-    /// Set absolute tolerance (different value per component)
-    pub fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
+    fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
         self.a_tol = a_tol;
         self
     }
 
-    /// Set relative tolerance
-    pub fn r_tol(mut self, r_tol: f64) -> Self {
+    fn r_tol(mut self, r_tol: f64) -> Self {
         self.r_tol = r_tol;
         self
     }
@@ -145,14 +158,25 @@ impl<'a, const D: usize> Integrator<D> for BS3<D>
 where
     BS3<D>: 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<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D> {
         // Allocate storage for the 4 stages
-        let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::zeros(); Self::STAGES];
+        let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::zeros(); STAGES];
 
         // Stage 1: k1 = f(t, y)
         k[0] = (ode.f)(ode.t, ode.y, &ode.params);
@@ -261,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]

src/integrator/dormand_prince.rs

@@ -1,9 +1,13 @@
 use nalgebra::SVector;
 
-use super::{Integrator, StepResult};
+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];
@@ -73,25 +77,35 @@ where
             r_tol: 1e-8,
         }
     }
+}
 
-    /// The same absolute tolerance on every component.
-    pub fn a_tol(self, a_tol: f64) -> Self {
+impl<const D: usize> Tolerances<D> for DormandPrince45<D>
+where
+    DormandPrince45<D>: Integrator<D>,
+{
+    fn absolute_tolerance(&self) -> SVector<f64, D> {
+        self.a_tol
+    }
+
+    fn relative_tolerance(&self) -> f64 {
+        self.r_tol
+    }
+
+    fn a_tol(self, a_tol: f64) -> Self {
         Self {
             a_tol: SVector::<f64, D>::from_element(a_tol),
             r_tol: self.r_tol,
         }
     }
 
-    /// A per-component absolute tolerance.
-    pub fn a_tol_full(self, a_tol: SVector<f64, D>) -> Self {
+    fn a_tol_full(self, a_tol: SVector<f64, D>) -> Self {
         Self {
             a_tol,
             r_tol: self.r_tol,
         }
     }
 
-    /// The same relative tolerance on every component
-    pub fn r_tol(self, r_tol: f64) -> Self {
+    fn r_tol(self, r_tol: f64) -> Self {
         Self {
             a_tol: self.a_tol,
             r_tol,
@@ -164,24 +178,35 @@ impl<'a, const D: usize> Integrator<D> for DormandPrince45<D>
 where
     DormandPrince45<D>: 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<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D> {
-        let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::zeros(); Self::STAGES];
+        let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::zeros(); STAGES];
         let mut next_y = ode.y;
         let mut err = SVector::<T, D>::zeros();
         let mut rcont5 = SVector::<T, D>::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::<T, D>::zeros();
             for (j, item) in k.iter().enumerate().take(i) {
@@ -191,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;

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<Method<D>>` 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<f64>, _p: &()) -> Vector1<f64> {
+///     -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<const D: usize> {
+    /// Bogacki-Shampine 3(2).
+    BS3(BS3<D>),
+    /// Dormand-Prince 5(4).
+    DormandPrince45(DormandPrince45<D>),
+    /// Verner's "most efficient" 7(6).
+    Vern7(Vern7<D>),
+    /// Rosenbrock-W 2(3), for stiff problems.
+    Rosenbrock23(Rosenbrock23<D>),
+}
+
+/// 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<const D: usize> Method<D> {
+    /// 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<const D: usize> Integrator<D> for Method<D> {
+    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<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D> {
+        dispatch!(self, m => m.step(ode, h))
+    }
+
+    fn interpolate<T: Scalar>(
+        &self,
+        t_start: T,
+        t_end: T,
+        dense: &[SVector<T, D>],
+        t: T,
+    ) -> SVector<T, D> {
+        dispatch!(self, m => m.interpolate(t_start, t_end, dense, t))
+    }
+
+    fn compute_extra_stages<T: Scalar, P>(
+        &self,
+        ode: &ODE<T, D, P>,
+        t_start: T,
+        y_start: SVector<T, D>,
+        h: T,
+        main_stages: &[SVector<T, D>],
+    ) -> Vec<SVector<T, D>> {
+        dispatch!(self, m => m.compute_extra_stages(ode, t_start, y_start, h, main_stages))
+    }
+}
+
+impl<const D: usize> Tolerances<D> for Method<D> {
+    fn absolute_tolerance(&self) -> SVector<f64, D> {
+        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<f64, D>) -> 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<const D: usize> From<BS3<D>> for Method<D> {
+    fn from(integrator: BS3<D>) -> Self {
+        Method::BS3(integrator)
+    }
+}
+
+impl<const D: usize> From<DormandPrince45<D>> for Method<D> {
+    fn from(integrator: DormandPrince45<D>) -> Self {
+        Method::DormandPrince45(integrator)
+    }
+}
+
+impl<const D: usize> From<Vern7<D>> for Method<D> {
+    fn from(integrator: Vern7<D>) -> Self {
+        Method::Vern7(integrator)
+    }
+}
+
+impl<const D: usize> From<Rosenbrock23<D>> for Method<D> {
+    fn from(integrator: Rosenbrock23<D>) -> 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<const D: usize>(&self) -> Method<D> {
+        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<f64>, _p: &()) -> Vector2<f64> {
+        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<S: Integrator<2> + Copy>(integrator: S) -> (Vec<f64>, Vec<Vector2<f64>>, Vec<f64>) {
+        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::<f64, 2>::repeat(1e-11)
+        );
+        assert_eq!(method.relative_tolerance(), 1e-9);
+
+        let per = SVector::<f64, 2>::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<Method<2>> = [
+            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: crate::scalar::Scalar>(_t: T, y: nalgebra::Vector1<T>) -> nalgebra::Vector1<T> {
+            -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::<f64, 2>::repeat(1e-9));
+    }
+}

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<T: Scalar, const D: usize> {
     /// [`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<Vec<SVector<T, D>>>,
 }
 
 /// 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 D: usize> {
-    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<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D>;
@@ -88,6 +110,46 @@ pub trait Integrator<const D: usize> {
     }
 }
 
+/// An adaptive integrator's error tolerances.
+///
+/// Split out from [`Integrator`] rather than folded into it, because a
+/// fixed-step method has no tolerances to set and `ADAPTIVE` exists to admit
+/// such methods. Every adaptive integrator in this crate implements it, so a
+/// function generic over `S: Tolerances<D>` can set and read the tolerances of
+/// whichever one it is handed.
+///
+/// # Example
+///
+/// ```rust
+/// use ordinary_diffeq::prelude::*;
+/// use nalgebra::SVector;
+///
+/// /// Tightens whatever integrator it is given.
+/// fn high_accuracy<const D: usize, S: Tolerances<D>>(integrator: S) -> S {
+///     integrator.a_tol(1e-12).r_tol(1e-12)
+/// }
+///
+/// let vern7 = high_accuracy(Vern7::<3>::new());
+/// assert_eq!(vern7.relative_tolerance(), 1e-12);
+/// assert_eq!(vern7.absolute_tolerance(), SVector::<f64, 3>::from_element(1e-12));
+/// ```
+pub trait Tolerances<const D: usize>: Integrator<D> {
+    /// The absolute tolerance, per component.
+    fn absolute_tolerance(&self) -> SVector<f64, D>;
+
+    /// The relative tolerance.
+    fn relative_tolerance(&self) -> f64;
+
+    /// Sets the absolute tolerance for every component.
+    fn a_tol(self, a_tol: f64) -> Self;
+
+    /// Sets the absolute tolerance per component.
+    fn a_tol_full(self, a_tol: SVector<f64, D>) -> Self;
+
+    /// Sets the relative tolerance.
+    fn r_tol(self, r_tol: f64) -> Self;
+}
+
 #[cfg(test)]
 mod tests {
     use super::dormand_prince::*;
@@ -119,3 +181,101 @@ mod tests {
         }
     }
 }
+
+#[cfg(test)]
+mod tolerance_tests {
+    use super::bs3::BS3;
+    use super::dormand_prince::DormandPrince45;
+    use super::rosenbrock::Rosenbrock23;
+    use super::vern7::Vern7;
+    use super::Tolerances;
+    use crate::controller::PIController;
+    use crate::problem::Problem;
+    use nalgebra::{SVector, Vector1};
+
+    /// Sets both tolerances and reads them back, for one integrator.
+    fn round_trip<S: Tolerances<3>>(integrator: S) {
+        let set = integrator.a_tol(1e-11).r_tol(1e-9);
+        assert_eq!(set.absolute_tolerance(), SVector::<f64, 3>::repeat(1e-11));
+        assert_eq!(set.relative_tolerance(), 1e-9);
+    }
+
+    /// `a_tol` sets every component; `a_tol_full` sets them independently.
+    fn per_component<S: Tolerances<3>>(integrator: S) {
+        let r_tol = integrator.relative_tolerance();
+        let uniform = integrator.a_tol(1e-7);
+        assert_eq!(
+            uniform.absolute_tolerance(),
+            SVector::<f64, 3>::repeat(1e-7)
+        );
+
+        let per = SVector::<f64, 3>::new(1e-6, 1e-9, 1e-12);
+        let full = uniform.a_tol_full(per);
+        assert_eq!(full.absolute_tolerance(), per);
+        // The relative tolerance is left alone by either setter.
+        assert_eq!(full.relative_tolerance(), r_tol);
+    }
+
+    #[test]
+    fn test_tolerances_round_trip() {
+        round_trip(BS3::<3>::new());
+        round_trip(DormandPrince45::<3>::new());
+        round_trip(Vern7::<3>::new());
+        round_trip(Rosenbrock23::<3>::new());
+    }
+
+    #[test]
+    fn test_tolerances_per_component() {
+        per_component(BS3::<3>::new());
+        per_component(DormandPrince45::<3>::new());
+        per_component(Vern7::<3>::new());
+        per_component(Rosenbrock23::<3>::new());
+    }
+
+    /// Solves `y' = -y` with whatever integrator it is handed, and reports how
+    /// many steps it took. Generic over `S`, which is the point of the trait.
+    fn steps_taken<S: Tolerances<1> + Copy>(integrator: S, a_tol: f64, r_tol: f64) -> usize {
+        fn derivative(_t: f64, y: Vector1<f64>, _p: &()) -> Vector1<f64> {
+            -y
+        }
+
+        let ode = crate::ode::ODE::new(&derivative, 0.0, 10.0, Vector1::new(1.0), ());
+        let mut problem = Problem::new(
+            ode,
+            integrator.a_tol(a_tol).r_tol(r_tol),
+            PIController::default(),
+        );
+        problem.solve().times.len()
+    }
+
+    #[test]
+    fn test_generic_tightening_costs_steps() {
+        for (name, loose, tight) in [
+            (
+                "BS3",
+                steps_taken(BS3::<1>::new(), 1e-4, 1e-4),
+                steps_taken(BS3::<1>::new(), 1e-9, 1e-9),
+            ),
+            (
+                "DormandPrince45",
+                steps_taken(DormandPrince45::<1>::new(), 1e-4, 1e-4),
+                steps_taken(DormandPrince45::<1>::new(), 1e-9, 1e-9),
+            ),
+            (
+                "Vern7",
+                steps_taken(Vern7::<1>::new(), 1e-4, 1e-4),
+                steps_taken(Vern7::<1>::new(), 1e-9, 1e-9),
+            ),
+            (
+                "Rosenbrock23",
+                steps_taken(Rosenbrock23::<1>::new(), 1e-4, 1e-4),
+                steps_taken(Rosenbrock23::<1>::new(), 1e-9, 1e-9),
+            ),
+        ] {
+            assert!(
+                tight > loose,
+                "{name}: tightening the tolerance should cost steps, got {loose} loose and {tight} tight"
+            );
+        }
+    }
+}

src/integrator/rosenbrock.rs

@@ -1,7 +1,7 @@
 use nalgebra::{SMatrix, SVector};
 
 use super::super::ode::{ProblemFunction, ODE};
-use super::{Integrator, StepResult};
+use super::{Integrator, StepResult, Tolerances};
 use crate::scalar::Scalar;
 
 /// Strategy for when to update the Jacobian matrix
@@ -143,8 +143,8 @@ pub struct Rosenbrock23<const D: usize> {
     c32: f64,
     /// Coefficient d = 1/(2+√2) ≈ 0.29289321881345254 (this is gamma!)
     d: f64,
-    /// Absolute tolerance for error estimation
-    a_tol: f64,
+    /// Absolute tolerance for error estimation, per component
+    a_tol: SVector<f64, D>,
     /// Relative tolerance for error estimation
     r_tol: f64,
     /// Strategy for updating the Jacobian
@@ -165,7 +165,7 @@ impl<const D: usize> Rosenbrock23<D> {
         Self {
             c32: 6.0 + 2.0_f64.sqrt(),
             d: 1.0 / (2.0 + 2.0_f64.sqrt()),
-            a_tol: 1e-6,
+            a_tol: SVector::<f64, D>::from_element(1e-6),
             r_tol: 1e-3,
             jacobian_strategy: JacobianUpdateStrategy::default(),
             cached_jacobian: None,
@@ -175,18 +175,6 @@ impl<const D: usize> Rosenbrock23<D> {
         }
     }
 
-    /// Set the absolute tolerance
-    pub fn a_tol(mut self, a_tol: f64) -> Self {
-        self.a_tol = a_tol;
-        self
-    }
-
-    /// Set the relative tolerance
-    pub fn r_tol(mut self, r_tol: f64) -> Self {
-        self.r_tol = r_tol;
-        self
-    }
-
     /// Set the Jacobian update strategy
     pub fn jacobian_strategy(mut self, strategy: JacobianUpdateStrategy) -> Self {
         self.jacobian_strategy = strategy;
@@ -209,6 +197,31 @@ impl<const D: usize> Rosenbrock23<D> {
     }
 }
 
+impl<const D: usize> Tolerances<D> for Rosenbrock23<D> {
+    fn absolute_tolerance(&self) -> SVector<f64, D> {
+        self.a_tol
+    }
+
+    fn relative_tolerance(&self) -> f64 {
+        self.r_tol
+    }
+
+    fn a_tol(mut self, a_tol: f64) -> Self {
+        self.a_tol = SVector::<f64, D>::from_element(a_tol);
+        self
+    }
+
+    fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
+        self.a_tol = a_tol;
+        self
+    }
+
+    fn r_tol(mut self, r_tol: f64) -> Self {
+        self.r_tol = r_tol;
+        self
+    }
+}
+
 impl<const D: usize> Default for Rosenbrock23<D> {
     fn default() -> Self {
         Self::new()
@@ -216,10 +229,21 @@ impl<const D: usize> Default for Rosenbrock23<D> {
 }
 
 impl<const D: usize> Integrator<D> for Rosenbrock23<D> {
-    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<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D> {
         let t = ode.t;
@@ -273,8 +297,7 @@ impl<const D: usize> Integrator<D> for Rosenbrock23<D> {
         // Compute scalar error estimate using weighted norm
         let mut error_sum = 0.0;
         for i in 0..D {
-            let scale =
-                self.a_tol + self.r_tol * uprev[i].re().abs().max(u_final[i].re().abs());
+            let scale = self.a_tol[i] + self.r_tol * uprev[i].re().abs().max(u_final[i].re().abs());
             let weighted_error = error_vec[i].re() / scale;
             error_sum += weighted_error * weighted_error;
         }

src/integrator/vern7.rs

@@ -1,9 +1,16 @@
 use nalgebra::SVector;
 
 use super::super::ode::ODE;
-use super::{Integrator, StepResult};
+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)
@@ -92,21 +99,31 @@ where
             r_tol: 1e-8,
         }
     }
+}
 
-    /// Set absolute tolerance (same value for all components)
-    pub fn a_tol(mut self, a_tol: f64) -> Self {
+impl<const D: usize> Tolerances<D> for Vern7<D>
+where
+    Vern7<D>: Integrator<D>,
+{
+    fn absolute_tolerance(&self) -> SVector<f64, D> {
+        self.a_tol
+    }
+
+    fn relative_tolerance(&self) -> f64 {
+        self.r_tol
+    }
+
+    fn a_tol(mut self, a_tol: f64) -> Self {
         self.a_tol = SVector::<f64, D>::from_element(a_tol);
         self
     }
 
-    /// Set absolute tolerance (different value per component)
-    pub fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
+    fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
         self.a_tol = a_tol;
         self
     }
 
-    /// Set relative tolerance
-    pub fn r_tol(mut self, r_tol: f64) -> Self {
+    fn r_tol(mut self, r_tol: f64) -> Self {
         self.r_tol = r_tol;
         self
     }
@@ -279,24 +296,40 @@ impl<'a, const D: usize> Integrator<D> for Vern7<D>
 where
     Vern7<D>: 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<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D> {
         // Allocate storage for the 10 stages
-        let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::zeros(); Self::STAGES];
+        let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::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
@@ -459,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];
@@ -553,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;

src/lib.rs

@@ -16,10 +16,11 @@ 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;
-    pub use super::integrator::{Integrator, StepResult};
+    pub use super::integrator::{Integrator, StepResult, Tolerances};
     pub use super::problem::{Problem, Solution};
     pub use super::scalar::Scalar;
 }

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(
@@ -219,6 +219,7 @@ mod tests {
     use crate::callback::{stop, FnCallback};
     use crate::controller::PIController;
     use crate::integrator::dormand_prince::DormandPrince45;
+    use crate::integrator::Tolerances;
     use approx::assert_relative_eq;
     use nalgebra::Vector3;