From da479a378e91182416e8ae5da68974b47304270f Mon Sep 17 00:00:00 2001 From: Connor Johnstone Date: Wed, 30 Sep 2026 13:48:20 -0400 Subject: [PATCH] Made the tolerance trait separable from the integrator itself --- roadmap/README.md | 30 ++- roadmap/features/40-tolerances-trait.md | 129 +++++++++++++ .../41-runtime-integrator-selection.md | 174 ++++++++++++++++++ src/integrator/bs3.rs | 24 ++- src/integrator/dormand_prince.rs | 24 ++- src/integrator/mod.rs | 138 ++++++++++++++ src/integrator/rosenbrock.rs | 48 +++-- src/integrator/vern7.rs | 24 ++- src/lib.rs | 2 +- src/problem.rs | 1 + 10 files changed, 548 insertions(+), 46 deletions(-) create mode 100644 roadmap/features/40-tolerances-trait.md create mode 100644 roadmap/features/41-runtime-integrator-selection.md diff --git a/roadmap/README.md b/roadmap/README.md index 34b7e8a..a551874 100644 --- a/roadmap/README.md +++ b/roadmap/README.md @@ -58,6 +58,23 @@ 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 + +- [ ] **[Runtime Integrator Selection](features/41-runtime-integrator-selection.md)** + - 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 + ### Callbacks - [ ] **[Discrete Callbacks](features/05-discrete-callbacks.md)** @@ -336,17 +353,18 @@ 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 (5/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: 12.2% (5/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) -Last updated: 2025-10-24 +Last updated: 2026-09-30 diff --git a/roadmap/features/40-tolerances-trait.md b/roadmap/features/40-tolerances-trait.md new file mode 100644 index 0000000..570e89e --- /dev/null +++ b/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 Vern7 where Vern7: Integrator { + pub fn a_tol(mut self, a_tol: f64) -> Self { .. } + pub fn a_tol_full(mut self, a_tol: SVector) -> 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` 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: Integrator { + fn with_a_tol(self, a_tol: f64) -> Self; + fn with_a_tol_per_component(self, a_tol: SVector) -> Self; + fn with_r_tol(self, r_tol: f64) -> Self; +} +impl Tolerances for Vern7 { .. } +``` + +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: Integrator { + /// The absolute tolerance, per component. + fn absolute_tolerance(&self) -> SVector; + + /// 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) -> 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` 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`, 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` 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` 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 diff --git a/roadmap/features/41-runtime-integrator-selection.md b/roadmap/features/41-runtime-integrator-selection.md new file mode 100644 index 0000000..aa19463 --- /dev/null +++ b/roadmap/features/41-runtime-integrator-selection.md @@ -0,0 +1,174 @@ +# Feature: Runtime Integrator Selection + +**Status**: proposed (2026-09-30) + +## Overview + +Which integrator runs a problem is fixed at compile time, because `Integrator` +keeps its structural facts in associated constants: + +```rust +pub trait Integrator { + 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`, because it would have to give one +`EXTRA_STAGES` for all of them, and Vern7 has six where Dormand-Prince has none. + +`Box>` is not available either, and will not become available: +`step` and `interpolate` 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` 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 { + 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 { + BS3(BS3), + DormandPrince45(DormandPrince45), + Vern7(Vern7), + Rosenbrock23(Rosenbrock23), +} + +impl Integrator for Method { + 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`. 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 diff --git a/src/integrator/bs3.rs b/src/integrator/bs3.rs index cecd0d1..1a1da76 100644 --- a/src/integrator/bs3.rs +++ b/src/integrator/bs3.rs @@ -1,7 +1,7 @@ use nalgebra::SVector; use super::super::ode::ODE; -use super::{Integrator, StepResult}; +use super::{Integrator, StepResult, Tolerances}; use crate::scalar::Scalar; /// Bogacki-Shampine 3/2 integrator trait for tableau coefficients @@ -80,21 +80,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 Tolerances for BS3 +where + BS3: Integrator, +{ + fn absolute_tolerance(&self) -> SVector { + self.a_tol + } + + fn relative_tolerance(&self) -> f64 { + self.r_tol + } + + fn a_tol(mut self, a_tol: f64) -> Self { self.a_tol = SVector::::from_element(a_tol); self } - /// Set absolute tolerance (different value per component) - pub fn a_tol_full(mut self, a_tol: SVector) -> Self { + fn a_tol_full(mut self, a_tol: SVector) -> 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 } diff --git a/src/integrator/dormand_prince.rs b/src/integrator/dormand_prince.rs index 5c39b4b..ec2e799 100644 --- a/src/integrator/dormand_prince.rs +++ b/src/integrator/dormand_prince.rs @@ -1,6 +1,6 @@ use nalgebra::SVector; -use super::{Integrator, StepResult}; +use super::{Integrator, StepResult, Tolerances}; use crate::ode::ODE; use crate::scalar::Scalar; @@ -73,25 +73,35 @@ where r_tol: 1e-8, } } +} - /// The same absolute tolerance on every component. - pub fn a_tol(self, a_tol: f64) -> Self { +impl Tolerances for DormandPrince45 +where + DormandPrince45: Integrator, +{ + fn absolute_tolerance(&self) -> SVector { + self.a_tol + } + + fn relative_tolerance(&self) -> f64 { + self.r_tol + } + + fn a_tol(self, a_tol: f64) -> Self { Self { a_tol: SVector::::from_element(a_tol), r_tol: self.r_tol, } } - /// A per-component absolute tolerance. - pub fn a_tol_full(self, a_tol: SVector) -> Self { + fn a_tol_full(self, a_tol: SVector) -> 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, diff --git a/src/integrator/mod.rs b/src/integrator/mod.rs index b5796c9..3cfba63 100644 --- a/src/integrator/mod.rs +++ b/src/integrator/mod.rs @@ -88,6 +88,46 @@ pub trait Integrator { } } +/// 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` 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>(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::::from_element(1e-12)); +/// ``` +pub trait Tolerances: Integrator { + /// The absolute tolerance, per component. + fn absolute_tolerance(&self) -> SVector; + + /// 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) -> Self; + + /// Sets the relative tolerance. + fn r_tol(self, r_tol: f64) -> Self; +} + #[cfg(test)] mod tests { use super::dormand_prince::*; @@ -119,3 +159,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>(integrator: S) { + let set = integrator.a_tol(1e-11).r_tol(1e-9); + assert_eq!(set.absolute_tolerance(), SVector::::repeat(1e-11)); + assert_eq!(set.relative_tolerance(), 1e-9); + } + + /// `a_tol` sets every component; `a_tol_full` sets them independently. + fn per_component>(integrator: S) { + let r_tol = integrator.relative_tolerance(); + let uniform = integrator.a_tol(1e-7); + assert_eq!( + uniform.absolute_tolerance(), + SVector::::repeat(1e-7) + ); + + let per = SVector::::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 + Copy>(integrator: S, a_tol: f64, r_tol: f64) -> usize { + fn derivative(_t: f64, y: Vector1, _p: &()) -> Vector1 { + -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" + ); + } + } +} diff --git a/src/integrator/rosenbrock.rs b/src/integrator/rosenbrock.rs index a620a56..c692fc7 100644 --- a/src/integrator/rosenbrock.rs +++ b/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 { 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, /// Relative tolerance for error estimation r_tol: f64, /// Strategy for updating the Jacobian @@ -165,7 +165,7 @@ impl Rosenbrock23 { Self { c32: 6.0 + 2.0_f64.sqrt(), d: 1.0 / (2.0 + 2.0_f64.sqrt()), - a_tol: 1e-6, + a_tol: SVector::::from_element(1e-6), r_tol: 1e-3, jacobian_strategy: JacobianUpdateStrategy::default(), cached_jacobian: None, @@ -175,18 +175,6 @@ impl Rosenbrock23 { } } - /// 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 Rosenbrock23 { } } +impl Tolerances for Rosenbrock23 { + fn absolute_tolerance(&self) -> SVector { + self.a_tol + } + + fn relative_tolerance(&self) -> f64 { + self.r_tol + } + + fn a_tol(mut self, a_tol: f64) -> Self { + self.a_tol = SVector::::from_element(a_tol); + self + } + + fn a_tol_full(mut self, a_tol: SVector) -> Self { + self.a_tol = a_tol; + self + } + + fn r_tol(mut self, r_tol: f64) -> Self { + self.r_tol = r_tol; + self + } +} + impl Default for Rosenbrock23 { fn default() -> Self { Self::new() @@ -273,8 +286,7 @@ impl Integrator for Rosenbrock23 { // 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; } diff --git a/src/integrator/vern7.rs b/src/integrator/vern7.rs index ac5ec91..bc377e6 100644 --- a/src/integrator/vern7.rs +++ b/src/integrator/vern7.rs @@ -1,7 +1,7 @@ use nalgebra::SVector; use super::super::ode::ODE; -use super::{Integrator, StepResult}; +use super::{Integrator, StepResult, Tolerances}; use crate::scalar::Scalar; /// Verner 7 integrator trait for tableau coefficients @@ -92,21 +92,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 Tolerances for Vern7 +where + Vern7: Integrator, +{ + fn absolute_tolerance(&self) -> SVector { + self.a_tol + } + + fn relative_tolerance(&self) -> f64 { + self.r_tol + } + + fn a_tol(mut self, a_tol: f64) -> Self { self.a_tol = SVector::::from_element(a_tol); self } - /// Set absolute tolerance (different value per component) - pub fn a_tol_full(mut self, a_tol: SVector) -> Self { + fn a_tol_full(mut self, a_tol: SVector) -> 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 } diff --git a/src/lib.rs b/src/lib.rs index 6e2e39e..b8b065a 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -19,7 +19,7 @@ pub mod prelude { 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; } diff --git a/src/problem.rs b/src/problem.rs index 885a65d..9d58f43 100644 --- a/src/problem.rs +++ b/src/problem.rs @@ -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;