diff --git a/readme.md b/readme.md index ffb7480..554cfb1 100644 --- a/readme.md +++ b/readme.md @@ -23,7 +23,8 @@ See [benchmark report](VERN7_BENCHMARK_REPORT.md) for detailed performance analy ### Other Features - **Adaptive time stepping** with PI controller -- **Callback events** with zero-crossing detection +- **Callback events** with zero-crossing detection, whose effects can change the state, change the + parameters, or end the solution at the event - **Dense output interpolation** at any time point - **Parameters** in derivative and callback functions - **Lazy computation** of extra interpolation stages (Vern7) @@ -71,19 +72,68 @@ let ode = ODE::new(&derivative, 0.0, 6.3, y0, params); let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-6); let controller = PIController::default(); -let value_too_high = FnCallback { +let mut value_too_high = FnCallback { event: &|t: f64, _y: &Vector2, _p: &Params| { 5.0 - t }, effect: &stop, }; // Solve the problem -let mut problem = Problem::new(ode, dp45, controller).with_callback(&value_too_high); +let mut problem = Problem::new(ode, dp45, controller).with_callback(&mut value_too_high); let solution = problem.solve(); // Can interpolate solutions to whatever you want let _interpolated_answer = solution.interpolate(4.4); ``` +## Events that change the state + +A callback watches a function of the state for a zero crossing. When one is found the step is +shortened to land on the crossing, the state there is recorded, and the effect runs with the ODE at +that point — so the effect can move the state, and the integration continues from where it leaves it. + +That is what it takes to integrate through a discontinuity. A rotating body parameterized by +Modified Rodrigues Parameters, for instance, has a singularity at a full turn, avoided by swapping +the parameter vector for `-sigma / |sigma|^2` — the same attitude, different parameters — whenever +`|sigma|` passes 1: + +```rust +use nalgebra::{Matrix3, Vector3}; +use ordinary_diffeq::prelude::*; + +/// MRP kinematics for a constant angular velocity. +fn attitude_rate(_t: f64, s: Vector3, omega: &Vector3) -> Vector3 { + let skew = Matrix3::new(0.0, -s[2], s[1], s[2], 0.0, -s[0], -s[1], s[0], 0.0); + let b = Matrix3::identity() * (1.0 - s.dot(&s)) + skew * 2.0 + s * s.transpose() * 2.0; + b * omega * 0.25 +} + +let mut shadow_set = FnCallback { + event: &|_t: f64, s: &Vector3, _p: &Vector3| 1.0 - s.norm(), + effect: &|ode: &mut ODE<'_, f64, 3, Vector3>| { + ode.y = -ode.y / ode.y.norm_squared(); + }, +}; + +let ode = ODE::new(&attitude_rate, 0.0, 5.0, Vector3::zeros(), Vector3::new(0.0, 0.0, 1.0)); +let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12); +let mut problem = Problem::new(ode, vern7, PIController::default()) + .with_callback(&mut shadow_set); +let solution = problem.solve(); + +// The parameters never approach the singularity, having swapped at half a turn. +assert!(solution.states.iter().all(|s| s.norm() <= 1.0 + 1e-9)); +``` + +A jump like that one puts the same time in the solution twice, holding the state on each side of it, +and makes the step size controller forget a history that described the solution before the jump. +Interpolating at the jump reports the state the effect left, the one the trajectory goes on with. +`effect` takes `&mut self`, so a callback that implements the trait itself rather than going through +`FnCallback` can record that it has fired. `stop` ends the solution at the crossing, with the +crossing as its last point and the problem's own `t_end` left alone. + +If two events fall inside one step the earlier one fires and the step ends there; the later one is +found again on the step after it. + ## Choosing the integrator at runtime The integrator is normally a type parameter, fixed when the code is compiled. `Method` is an enum diff --git a/roadmap/README.md b/roadmap/README.md index 46c0a05..51fd345 100644 --- a/roadmap/README.md +++ b/roadmap/README.md @@ -8,7 +8,8 @@ The library currently has: - ✅ Dormand-Prince 4(5) adaptive explicit RK method with dense output - ✅ Tsit5 explicit RK method with dense output - ✅ PI step size controller -- ✅ Basic continuous callbacks with zero-crossing detection +- ✅ Continuous callbacks with zero-crossing detection, whose effects can change + the state, the parameters, or end the solution - ✅ Solution interpolation interface - ✅ Generic over array dimensions at compile-time - ✅ Support for ordinary and second-order ODE problems @@ -79,6 +80,18 @@ Each feature below links to a detailed implementation plan in the `features/` di ### Callbacks +- [x] **[Callback Effects That Can Change the State](features/42-callback-effects.md)** ✅ COMPLETED + - An effect's write to the state was overwritten by the step it was located in, + so a callback could stop a solution but not change one + - `effect(&mut self)`, termination at the crossing, both sides of a jump + recorded, the controller reset across one, and a defined order for two events + in one step + - Raised by `astrokit`, which needs a Modified Rodrigues Parameter shadow set + swap mid-propagation + - **Dependencies**: None + - **Effort**: Medium + - **Status**: All success criteria met; benches unchanged within run-to-run drift + - [ ] **[Discrete Callbacks](features/05-discrete-callbacks.md)** - Event detection based on conditions (not zero-crossings) - Useful for time-based events, iteration counts, etc. @@ -93,13 +106,13 @@ Each feature below links to a detailed implementation plan in the `features/` di ### Solution Interface -- [ ] **[Saveat Functionality](features/07-saveat.md)** +- [ ] **[Saveat Functionality](features/07-saveat-functionality.md)** - Save solution at specific timepoints - Dense vs sparse saving strategies - **Dependencies**: None - **Effort**: Medium -- [ ] **[Solution Derivatives](features/08-solution-derivatives.md)** +- [ ] **Solution Derivatives** (no plan written yet) - Access derivatives at any time point via interpolation - `solution.derivative(t)` interface - **Dependencies**: None @@ -111,7 +124,7 @@ Each feature below links to a detailed implementation plan in the `features/` di ### Algorithms -- [ ] **[DP8 (Dormand-Prince 8th order)](features/09-dp8-method.md)** +- [ ] **DP8 (Dormand-Prince 8th order)** (no plan written yet) - 8th order explicit RK for very high accuracy - Complements Vern7 for algorithm selection - **Dependencies**: None @@ -123,7 +136,7 @@ Each feature below links to a detailed implementation plan in the `features/` di - **Dependencies**: Linear solver infrastructure, Nordsieck representation - **Effort**: Large -- [ ] **[Rodas4/Rodas5P](features/11-rodas-methods.md)** +- [ ] **[Rodas4/Rodas5P](features/11-rodas4-rodas5p-methods.md)** - Higher-order Rosenbrock methods (4th/5th order) - Better accuracy for stiff problems - **Dependencies**: Rosenbrock23 @@ -138,7 +151,7 @@ Each feature below links to a detailed implementation plan in the `features/` di - **Dependencies**: At least one stiff solver (Rosenbrock23 or FBDF) - **Effort**: Large -- [ ] **[Default Algorithm Selection](features/13-default-algorithm.md)** +- [ ] **[Default Algorithm Selection](features/13-default-algorithm-selection.md)** - Smart defaults based on problem characteristics - `solve(problem)` without specifying algorithm - **Dependencies**: Auto-switching, multiple algorithms @@ -146,7 +159,7 @@ Each feature below links to a detailed implementation plan in the `features/` di ### Initialization -- [ ] **[Automatic Initial Step Size](features/14-initial-stepsize.md)** +- [ ] **[Automatic Initial Step Size](features/14-automatic-initial-step-size.md)** - Algorithm to determine good initial dt - Based on local Lipschitz estimate - **Dependencies**: None @@ -154,19 +167,19 @@ Each feature below links to a detailed implementation plan in the `features/` di ### Callbacks -- [ ] **[PresetTimeCallback](features/15-preset-time-callback.md)** +- [ ] **[PresetTimeCallback](features/15-presettimecallback.md)** - Trigger callbacks at specific predetermined times - Important for time-varying forcing functions - **Dependencies**: Discrete callbacks - **Effort**: Small -- [ ] **[TerminateSteadyState](features/16-terminate-steady-state.md)** +- [ ] **[TerminateSteadyState](features/16-terminatesteadystate.md)** - Auto-detect when solution reaches steady state - Stop integration early - **Dependencies**: Discrete callbacks - **Effort**: Small -- [ ] **[SavingCallback](features/17-saving-callback.md)** +- [ ] **[SavingCallback](features/17-savingcallback.md)** - Custom saving logic beyond default - For memory-efficient large-scale simulations - **Dependencies**: CallbackSet @@ -194,7 +207,7 @@ Each feature below links to a detailed implementation plan in the `features/` di ### Algorithms -- [ ] **[Low-Storage Runge-Kutta](features/20-low-storage-rk.md)** +- [ ] **[Low-Storage Runge-Kutta](features/20-low-storage-runge-kutta.md)** - 2N, 3N, 4N storage variants - Critical for large-scale problems - **Dependencies**: None @@ -212,7 +225,7 @@ Each feature below links to a detailed implementation plan in the `features/` di - **Dependencies**: Second-order ODE infrastructure (already exists) - **Effort**: Medium -- [ ] **[Verner Methods Suite](features/23-verner-methods.md)** +- [ ] **[Verner Methods Suite](features/23-verner-methods-suite.md)** - Vern6, Vern8, Vern9 - Complete Verner family for different accuracy needs - **Dependencies**: Vern7 @@ -258,19 +271,19 @@ Each feature below links to a detailed implementation plan in the `features/` di ### Advanced Callbacks -- [ ] **[VectorContinuousCallback](features/30-vector-continuous-callback.md)** +- [ ] **[VectorContinuousCallback](features/30-vectorcontinuouscallback.md)** - Multiple simultaneous event detection - More efficient than multiple callbacks - **Dependencies**: CallbackSet - **Effort**: Medium -- [ ] **[PositiveDomain](features/31-positive-domain.md)** +- [ ] **[PositiveDomain](features/31-positivedomain.md)** - Enforce positivity constraints - Important for physical systems - **Dependencies**: CallbackSet - **Effort**: Small -- [ ] **[ManifoldProjection](features/32-manifold-projection.md)** +- [ ] **[ManifoldProjection](features/32-manifoldprojection.md)** - Project solution onto constraint manifolds - For constrained mechanical systems - **Dependencies**: CallbackSet @@ -278,14 +291,14 @@ Each feature below links to a detailed implementation plan in the `features/` di ### Infrastructure -- [ ] **[Nonlinear Solver Infrastructure](features/33-nonlinear-solver.md)** +- [ ] **[Nonlinear Solver Infrastructure](features/33-nonlinear-solver-infrastructure.md)** - Newton's method - Quasi-Newton methods - Generic nonlinear solver interface - **Dependencies**: Linear solver infrastructure - **Effort**: Large -- [ ] **[Krylov Linear Solvers](features/34-krylov-solvers.md)** +- [ ] **[Krylov Linear Solvers](features/34-krylov-linear-solvers.md)** - GMRES, BiCGStab - For large sparse systems - **Dependencies**: Linear solver infrastructure @@ -355,12 +368,12 @@ Each algorithm implementation should include: ## Progress Tracking -Total Features: 41 -- Tier 1: 10 features (6/10 complete) ✅ +Total Features: 42 +- Tier 1: 11 features (7/11 complete) ✅ - Tier 2: 11 features (0/11 complete) - Tier 3: 20 features (0/20 complete) -**Overall Progress: 14.6% (6/41 features complete)** +**Overall Progress: 16.7% (7/42 features complete)** ### Completed Features 1. ✅ BS3 (Bogacki-Shampine 3/2) - Tier 1 (2025-10-23) @@ -369,5 +382,6 @@ Total Features: 41 4. ✅ Rosenbrock23 - Tier 1 (2025-10-24) 5. ✅ Tolerances Trait - Tier 1 (2026-09-30) 6. ✅ Runtime Integrator Selection - Tier 1 (2026-09-30) +7. ✅ Callback Effects That Can Change the State - Tier 1 (2026-10-01) -Last updated: 2026-09-30 +Last updated: 2026-10-01 diff --git a/roadmap/features/42-callback-effects.md b/roadmap/features/42-callback-effects.md new file mode 100644 index 0000000..6a39190 --- /dev/null +++ b/roadmap/features/42-callback-effects.md @@ -0,0 +1,231 @@ +# Feature: Callback Effects That Can Change the State + +**Status**: completed (2026-10-01) + +## Overview + +Continuous callbacks located their crossings correctly and then threw away what +the effect did about them. In `Problem::solve`: + +```rust +callback.effect(&mut self.ode); // writes to ode.y +} +} +self.ode.y = new_y; // and here it goes +self.ode.t += curr_step; +``` + +The effect ran on the pre-step state, and the two lines after the callback loop +overwrote `ode.y` with the state the integrator had produced. Anything the effect +wrote to the state was lost. `stop` worked only because it writes `t_end` +instead. + +Reproduced with `dy/dt = 1` from 0 to 10, `y0 = 0`, an event of `5 - y[0]` and an +effect of `ode.y = Vector1::new(100.0)`: the solution ended at `y = 10`, exactly +as it does with an empty effect, where it should end near 105. + +So the machinery was in place for the easy half of event handling — detect a +condition and stop — and could not do the other half: detect a condition, change +the state, and carry on. + +## Why This Feature Matters + +A rotating rigid body parameterized by Modified Rodrigues Parameters has a +singularity at a full turn. It is avoided by swapping the parameter vector for +the shadow set `-sigma / |sigma|^2` — the same attitude, different parameters — +whenever `|sigma|` exceeds 1. The event is `1 - |sigma|`, the effect is that +substitution, and the integration has to continue from the substituted state. +Any arc longer than a half turn crosses it, so this is the common case rather +than an edge case. + +The same shape of problem covers a body entering or leaving shadow, a vehicle +whose mass or area changes partway through a run, and a bouncing ball. All of +them are a state discontinuity located by an event, which is the one thing the +callbacks could not express. + +## What Was Done + +### The effect runs after the crossing is committed + +The event branch now commits the shortened step before calling the effect: + +1. the step is shortened to land on the crossing, and the state there is + recorded; +2. `ode.t` and `ode.y` are set to that crossing; +3. the effect is called, with the ODE at the crossing; +4. whatever the effect leaves in `ode.y` is what the integration continues from. + +Detection still reads the state at the start of the step — an event function is +asked about both ends — so nothing can be committed until every callback has been +asked. That is why this is a reordering of the branch rather than a swap of two +lines. + +### `effect(&mut self)` + +`Callback::effect` takes `&mut self`, so a callback can record that it has fired, +count its firings, or keep anything else it needs across them. `Problem` holds +`Vec<&'a mut dyn Callback>` and `with_callback` takes `&'a mut dyn +Callback`, which is the breaking part of the change: a registered callback has to +be declared `mut` and passed as `&mut`. + +`FnCallback` still holds its two closures by shared reference and so still has no +state of its own; a callback that needs state implements `Callback` on a struct, +which is what [Discrete Callbacks](05-discrete-callbacks.md) will do with its +`has_fired` field. + +`Problem` lost its `#[derive(Clone)]`, which a `&mut` field cannot support. It +was never used. + +### A terminal effect ends the solution at the crossing + +`ODE::finished` existed, was never set and was never read. It now means what its +name says: `ODE::terminate` sets it, the solve loop checks it, and `stop` is +`ode.terminate()`. + +Behaviorally this is close to what `stop` already did. Writing `t_end = ode.t` +also ended the loop, and because the step had already been shortened to the +crossing, the last recorded point was the crossing — the solution did not in fact +run a step past the event. What the old `stop` did do was mutate the problem's own +end time, so a solution could no longer report the `t_end` it was asked for. Now +`t_end` is left alone and `finished` says why the solution is short. + +### Both sides of a jump are recorded + +Where an effect changes the state, the trajectory genuinely has two values at +that time, and the solution holds both: the same `t` appears twice in `times`, +with the state at the crossing and the state the effect left. The `dense` entry +between them spans zero time and is never interpolated into — `interpolate` at +exactly that time hits the search's exact-match branch, and walks to the last +entry at that time, so it reports the state the effect left. That is the value the +trajectory carries forward, and it is the one a caller asking "where is the body +at this instant" wants; the value it arrived with is the entry before it in +`states`. + +This is what DifferentialEquations.jl does, where the two are `u_left` and +`u_right`. + +### The controller is reset across a jump + +`Controller` grew a required `reset` method, which discards the error history and +leaves the next step guess alone. `solve` calls it when an effect changes the +state. The history is a claim about how the error behaves along a smooth +solution, and across a discontinuity it is not about anything. + +Keyed on the state changing rather than on the effect running at all, so that a +callback which only observes does not throw away a good history every time it +fires. An effect that changes only the parameters is a discontinuity the solver +cannot see, and does not reset the controller. + +`PIDController::reset` already existed as an inherent method and became the trait +implementation; `PIController::reset` restores `factor_old` to its initial value. + +### Two events in one step + +Defined and tested: the earliest crossing in the step fires, the step ends there, +and any later event is still ahead and is found again on the following step. A +tie goes to registration order. + +The old behavior was to walk the callbacks in registration order, each shortening +what the last one left. That could apply a later event's effect and then truncate +the step to before it, so the effect was applied at a time the solution never +reached. One pass to find the earliest crossing costs an extra root find only +when more than one callback crosses in the same step. + +### Crossings at exactly zero, and not firing twice at one + +The test was `before * after < 0.0`, which misses a crossing that starts or ends +exactly on zero. Now: + +- `before == 0.0` is a crossing at the step's own start, located at `h = 0`; +- `after == 0.0` is a crossing at the step's end; +- opposite signs are located with Brent as before. + +The first of those needs a guard, because an effect is free to leave its event +function sitting on the zero it fired at. The shadow set swap does exactly that: +`|sigma|` is 1 on both sides of it, so `1 - |sigma|` is still zero after the +effect. Without a guard the callback is re-detected at the start of the next +step, located at `h = 0`, fired again, forever. + +The guard is that the step beginning at a firing searches from +`EVENT_PROBE_FRACTION` of itself — a millionth — rather than from its start, for +that callback alone. Both the `before` sign and Brent's bracket start there, so +the event just handled is out of reach while a second crossing later in the same +step is still found. `solve` tracks the time each callback last fired to know +when to do this. + +The first attempt at this guard was to skip a callback for the whole step that +begins at its firing, which is simpler and wrong: a bouncing ball bounced once +and then fell through the floor, because the step after the bounce contained the +next landing and the callback was not looking. That is in the test suite now. + +A fraction of a step is not a satisfying thing to have to choose, and this is the +least bad option found rather than a good one. Alternatives considered: an +absolute tolerance on the event value, which is unit-dependent and no better +founded; taking the sign from the step before the firing, which still leaves +Brent bracketing a root at `h = 0`; and refusing a located crossing below some +`h`, which is the same arbitrary fraction one step further down. If something +cleaner turns up, this is the place to change it. + +## Known Limitations + +- An even number of crossings inside one step is still missed: the two ends of + the step agree in sign and nothing looks in between. Documented on + `Problem::solve`. +- A crossing inside the first millionth of the step that begins at a firing is + missed, which is the cost of the no-refiring guard. It takes an event that + recurs almost immediately after one was handled. +- The crossing time is located in real arithmetic and carries no derivative. + Deliberate: differentiating through Brent's iterations gives the derivative of + the iterative procedure rather than the sensitivity of the event, which is + wrong without being loud about it. The state at the crossing comes from one + ordinary step of the integrator, so it does carry derivatives. +- An effect that changes the parameters rather than the state does not reset the + controller, and does not record two points. + +## Testing + +In `src/problem.rs`: + +- the reproduction above, ending near 105 rather than 10; +- both sides of the jump present in `times`/`states`, interpolation either side of + it following the branch it is on, and interpolation at the jump itself giving + the state the effect left; +- `stop` ending the solution at the crossing, with nothing past it, `t_end` + untouched and `finished` set; +- a callback counting its own firings and acting only on the first, which needs + the `&mut self`; +- two events in one step firing earliest-first, with the registration order + reversed to show it does not matter; +- an effect that lands back on its own event firing once there rather than + looping, and firing again when the event genuinely recurs; +- a bouncing ball, bouncing at the three analytic bounce times inside the + interval and interpolating correctly inside each arc; +- an event function already zero at the start of the integration firing there; +- the controller reset asked for once across a jump and not at all for an effect + that changes nothing, measured with a controller that counts the resets; +- a jump under dual numbers, where the derivative through the effect is exact + (in `src/lib.rs`, with the autodiff tests); +- the Modified Rodrigues Parameter swap itself, against the closed form for a + constant angular velocity: one swap, at exactly half a turn, `|sigma| <= 1` + throughout, and the final attitude equal to `-cot(theta / 4)` about the + rotation axis. + +## Success Criteria + +- An effect that writes to `ode.y` changes what is integrated from. ✅ +- An effect can record something about itself. ✅ +- A terminated solution ends at the crossing. ✅ +- A jump is recorded on both sides. ✅ +- The step size controller forgets its history across a jump. ✅ +- The order of two events in one step is defined and tested. ✅ +- No measurable change on the criterion benches. ✅ Measured before and after on + the same machine in one sitting: `bench_orbit` -1.02%, `bench_simple_1d` + +0.36%, `bench_interpolation_1d` -0.67%, all three within criterion's noise + threshold. Which is the right reading of them: re-running the *identical* + binary moved `bench_simple_1d` by -0.80% and `bench_interpolation_1d` by + +1.47%, so run-to-run drift here is about a percent either way and is larger + than the difference being looked for. No bench registers a callback, so what + they measure is the no-callback path, which gains a loop over an empty slice + and a `finished` check per step. + +**Effort**: Medium diff --git a/src/callback.rs b/src/callback.rs index 5643513..a740e2f 100644 --- a/src/callback.rs +++ b/src/callback.rs @@ -1,31 +1,96 @@ +//! Conditions that interrupt an integration, and what they do when they fire. + use nalgebra::SVector; use crate::ode::ODE; use crate::scalar::Scalar; -/// A callback function you can utilize to change your model mid-run without stopping. +/// A condition the solver watches for, and the change it makes when it is met. /// /// The solver evaluates [`event`](Callback::event) at each end of every accepted -/// step and looks for a sign change. When it finds one it locates the crossing -/// and applies [`effect`](Callback::effect) to the ODE there. +/// step and looks for a sign change. When it finds one it shortens the step to +/// land on the crossing, records the state there, and then calls +/// [`effect`](Callback::effect) on the ODE. Whatever the effect leaves in +/// `ode.y` is the state the integration carries on from. +/// +/// An effect may: +/// +/// - change the state, which records a jump: the solution holds both the value +/// at the crossing and the value the effect left, at the same time (see +/// [`Solution`](crate::problem::Solution)); +/// - change the parameters, which changes the right-hand side from the crossing +/// onward; +/// - end the solution at the crossing, by calling [`stop`] or +/// [`ODE::terminate`]; +/// - record something about itself, since it takes `&mut self`. +/// +/// # Example +/// +/// A fire-once callback, which needs the `&mut self` on the effect to remember +/// that it has run: +/// +/// ``` +/// use nalgebra::Vector1; +/// use ordinary_diffeq::prelude::*; +/// +/// struct Kick { +/// fired: usize, +/// } +/// +/// impl Callback for Kick { +/// fn event(&self, _t: f64, y: &Vector1, _p: &()) -> f64 { +/// 5.0 - y[0] +/// } +/// +/// fn effect(&mut self, ode: &mut ODE<'_, f64, 1, ()>) { +/// self.fired += 1; +/// ode.y = Vector1::new(100.0); +/// } +/// } +/// +/// fn derivative(_t: f64, _y: Vector1, _p: &()) -> Vector1 { +/// Vector1::new(1.0) +/// } +/// +/// let ode = ODE::new(&derivative, 0.0, 10.0, Vector1::new(0.0), ()); +/// let mut kick = Kick { fired: 0 }; +/// let mut problem = Problem::new(ode, DormandPrince45::new(), PIController::default()) +/// .with_callback(&mut kick); +/// let solution = problem.solve(); +/// +/// // The state jumped to 100 at t = 5 and kept integrating from there. +/// assert!((solution.states.last().unwrap()[0] - 105.0).abs() < 1e-9); +/// ``` +pub trait Callback { + /// The function whose zero crossings mark the event. + /// + /// Called repeatedly while the crossing is being located, so it must not + /// depend on how many times it has been called. + fn event(&self, t: T, y: &SVector, params: &P) -> T; + + /// What to do to the ODE once the crossing has been located. + /// + /// `ode.t` and `ode.y` are the time and state at the crossing when this is + /// called. + fn effect(&mut self, ode: &mut ODE<'_, T, D, P>); +} + +/// A [`Callback`] built from a pair of closures. +/// +/// The closures are held by shared reference, so neither can carry state of its +/// own. A callback that has to remember something — that it has already fired, +/// how many times, what the state was — implements [`Callback`] on a struct +/// instead; the example on that trait shows it. +/// /// ``` /// # use ordinary_diffeq::prelude::*; /// # use nalgebra::Vector2; /// # type Params = (); -/// let too_high = FnCallback { +/// let mut too_high = FnCallback { /// event: &|_t: f64, y: &Vector2, _p: &Params| 10.0 - y[0], /// effect: &stop, /// }; /// ``` -pub trait Callback { - /// The function whose zero crossings mark the event. - fn event(&self, t: T, y: &SVector, params: &P) -> T; - - /// What to do to the ODE once the crossing has been located. - fn effect(&self, ode: &mut ODE<'_, T, D, P>); -} - -/// A [`Callback`] built from a pair of closures. #[derive(Clone, Copy)] pub struct FnCallback<'a, T: Scalar, const D: usize, P> { /// The function to check for zero crossings @@ -42,14 +107,14 @@ impl Callback for FnCallback<'_, T, D, P> } #[inline] - fn effect(&self, ode: &mut ODE<'_, T, D, P>) { + fn effect(&mut self, ode: &mut ODE<'_, T, D, P>) { (self.effect)(ode) } } -/// A convenience function for stopping the integration +/// An effect that ends the solution at the crossing. pub fn stop(ode: &mut ODE<'_, T, D, P>) { - ode.t_end = ode.t; + ode.terminate(); } #[cfg(test)] diff --git a/src/controller.rs b/src/controller.rs index 6ac9ecc..b3b7323 100644 --- a/src/controller.rs +++ b/src/controller.rs @@ -34,6 +34,14 @@ pub trait Controller { fn next_step_guess(&self) -> TryStep; fn set_next_step_guess(&mut self, guess: TryStep); + + /// Discards the error history, leaving the next step guess alone. + /// + /// The history is an estimate of how the error behaves along a smooth + /// solution, so it says nothing across a discontinuity. + /// [`Problem::solve`](crate::problem::Problem::solve) calls this when a + /// callback effect changes the state. + fn reset(&mut self); } #[derive(Debug, Clone, Copy)] @@ -82,6 +90,9 @@ impl Controller for PIController { self.next_step_guess = guess; } + fn reset(&mut self) { + self.factor_old = 1.0e-4; + } } impl PIController { @@ -247,6 +258,11 @@ impl Controller for PIDController { self.next_step_guess = guess; } + fn reset(&mut self) { + self.err_old = 0.0; + self.err_older = 0.0; + self.h_old = 0.0; + } } impl PIDController { @@ -311,15 +327,6 @@ impl PIDController { 1e-4, // initial_h ) } - - /// Reset the controller's error history - /// - /// Useful when switching algorithms or restarting integration - pub fn reset(&mut self) { - self.err_old = 0.0; - self.err_older = 0.0; - self.h_old = 0.0; - } } impl Default for PIDController { diff --git a/src/lib.rs b/src/lib.rs index 36e7b3e..e6d2798 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -52,13 +52,13 @@ mod tests { let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-6); let controller = PIController::default(); - let value_too_high = FnCallback { + let mut value_too_high = FnCallback { event: &|t: f64, _y: &Vector2, _p: &Params| 5.0 - t, effect: &stop, }; // Solve the problem - let mut problem = Problem::new(ode, dp45, controller).with_callback(&value_too_high); + let mut problem = Problem::new(ode, dp45, controller).with_callback(&mut value_too_high); let solution = problem.solve(); // Can interpolate solutions to whatever you want @@ -322,6 +322,48 @@ mod autodiff_tests { } } + /// Derivatives survive a callback that changes the state. + /// + /// `y' = y` from `y0`, halved the moment it reaches 2, run to `t = 2`. It + /// reaches 2 twice in that interval — at `ln 2` and at `2 ln 2` — and `n` + /// halvings leave `y(2) = e^2 y0 / 2^n` whatever the times were, so the + /// solution and its sensitivity to `y0` are both `e^2 / 4` exactly. + /// + /// The crossing time is located in real arithmetic and carries no + /// derivative. This problem is one where that costs nothing: the state at + /// the crossing is fixed at 2 by the event itself, and the final state turns + /// out not to depend on when the crossing happened. A problem where the + /// event time does matter needs the sensitivity of the event condition, + /// which forward mode through Brent's iterations would not give. + #[test] + fn test_ad_through_a_callback_jump() { + fn derivative(_t: T, y: Vector1, _p: &()) -> Vector1 { + y + } + + let y0 = Dual64::from_re(1.0).derivative(); + let mut halve = FnCallback { + event: &|_t: Dual64, y: &Vector1, _p: &()| Dual64::lift(2.0) - y[0], + effect: &|ode: &mut ODE<'_, Dual64, 1, ()>| ode.y /= Dual64::lift(2.0), + }; + + let ode = ODE::new( + &derivative, + Dual64::lift(0.0), + Dual64::lift(2.0), + Vector1::new(y0), + (), + ); + let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12); + let mut problem = + Problem::new(ode, vern7, PIController::default()).with_callback(&mut halve); + let solution = problem.solve(); + let y = solution.states.last().unwrap()[0]; + + assert_relative_eq!(y.re, 2.0f64.exp() / 4.0, max_relative = 1e-10); + assert_relative_eq!(y.eps, 2.0f64.exp() / 4.0, max_relative = 1e-10); + } + /// The same, through Vern7 rather than Dormand-Prince, which exercises a /// path Dormand-Prince does not have: Vern7 keeps six of its sixteen stages /// back and computes them only when something interpolates inside a step. diff --git a/src/ode.rs b/src/ode.rs index 0fad207..5eeaf7c 100644 --- a/src/ode.rs +++ b/src/ode.rs @@ -23,6 +23,9 @@ pub struct ODE<'a, T: Scalar, const D: usize, P> { pub t0: T, pub t_end: T, pub h: T, + + /// Set when the solution has been ended short of `t_end`, by a callback + /// effect calling [`terminate`](ODE::terminate). pub finished: bool, } @@ -45,6 +48,15 @@ impl<'a, T: Scalar, const D: usize, P> ODE<'a, T, D, P> { finished: false, } } + + /// Ends the solution here, without disturbing `t_end`. + /// + /// Meant for a callback effect: `Problem::solve` records the state at the + /// crossing and then stops. [`stop`](crate::callback::stop) is this as a + /// ready-made effect. + pub fn terminate(&mut self) { + self.finished = true; + } } #[cfg(test)] diff --git a/src/problem.rs b/src/problem.rs index 6d903a1..e76ad90 100644 --- a/src/problem.rs +++ b/src/problem.rs @@ -1,4 +1,11 @@ //! Assembling an ODE, an integrator and a controller +//! +//! Event location is done in real arithmetic: the event condition is a +//! real-valued question, and Brent's method differentiated through would give +//! the derivative of its own iteration rather than the sensitivity of the +//! event, which is wrong in a way that does not announce itself. So the search +//! runs on `re()` projections and the state at the crossing is produced by one +//! ordinary step of the integrator, which does carry derivatives. use nalgebra::SVector; use roots::{find_root_brent, SimpleConvergency}; @@ -10,7 +17,22 @@ use super::integrator::Integrator; use super::ode::ODE; use super::scalar::Scalar; -#[derive(Clone)] +/// How far into a step the event function of a callback that just fired is read, +/// as a fraction of the step. +/// +/// An effect is free to leave its event function sitting on the zero it fired +/// at — the shadow set swap of a Modified Rodrigues Parameter vector does +/// exactly that — so reading it at the start of the next step would locate the +/// same event again, at no advance in time, forever. Reading a little way in +/// instead gives the sign the event function has after the effect, and the +/// search for a crossing starts from the same place, so a second crossing +/// further into the same step is still found. +/// +/// Small enough that a genuine crossing is very unlikely to hide inside it, and +/// large enough that the event function has moved clear of its own rounding at +/// the event. +const EVENT_PROBE_FRACTION: f64 = 1e-6; + pub struct Problem<'a, T: Scalar, const D: usize, S, P, C> where S: Integrator, @@ -19,7 +41,7 @@ where ode: ODE<'a, T, D, P>, integrator: S, controller: C, - callbacks: Vec<&'a dyn Callback>, + callbacks: Vec<&'a mut dyn Callback>, } impl<'a, T: Scalar, const D: usize, S, P, C> Problem<'a, T, D, S, P, C> @@ -36,6 +58,27 @@ where } } + /// Integrates from `t0` to `t_end`, or to wherever a callback ends it. + /// + /// Each accepted step is checked against every registered callback. When an + /// event function changes sign across a step, the step is shortened to land + /// on the crossing, the state there is recorded, and the callback's effect + /// runs with the ODE at that point — so an effect may change the state, the + /// parameters, or end the solution, and the integration carries on from + /// whatever it leaves behind. + /// + /// If more than one event falls inside one step, the earliest fires and the + /// step ends there; the others are still ahead and are found again on the + /// following step. Ties go to the order the callbacks were registered. + /// + /// A callback does not fire twice at the same time: an effect may leave its + /// event function sitting on the zero it fired at, and the step that begins + /// at a firing searches from a hair inside itself rather than from its start. + /// A second crossing later in that step is still found — a bouncing ball + /// keeps bouncing — but one inside that first millionth of the step is not. + /// + /// An even number of crossings inside one step is missed, in that the step's + /// two ends agree in sign and nothing looks in between. pub fn solve(&mut self) -> Solution<'_, T, S, D, P> { let mut convergency = SimpleConvergency { eps: 1e-12f64, @@ -44,7 +87,10 @@ where let mut times: Vec = vec![self.ode.t]; let mut states: Vec> = vec![self.ode.y]; let mut dense_coefficients: Vec>>> = Vec::new(); - while self.ode.t.re() < self.ode.t_end.re() { + // When each callback last fired, which is where its next search for a + // crossing starts from rather than at — see `EVENT_PROBE_FRACTION`. + let mut last_fired: Vec> = vec![None; self.callbacks.len()]; + while !self.ode.finished && self.ode.t.re() < self.ode.t_end.re() { if (self.ode.t + self.controller.next_step_guess().extract()).re() > self.ode.t_end.re() { @@ -52,7 +98,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 self.integrator.adaptive() { + let (new_y, curr_step, dense_option) = if self.integrator.adaptive() { // First, we try stepping with the "next step guess" to get the error let mut result = self .integrator @@ -85,32 +131,94 @@ where result.dense, ) }; - if !self.callbacks.is_empty() { - // Check for events occurring - for callback in &self.callbacks { - // The event condition is a real-valued question, so it is - // asked in the reals — see this module's header. - let before = callback.event(self.ode.t, &self.ode.y, &self.ode.params).re(); - let after = callback - .event(self.ode.t + curr_step, &new_y, &self.ode.params) - .re(); - if before * after < 0.0 { - // If the event crossed zero, then find the root - let f = |test_h: f64| { - let h = T::lift(test_h); - let test_y = self.integrator.step(&self.ode, h).y; - callback.event(self.ode.t + h, &test_y, &self.ode.params).re() - }; - let root = - find_root_brent(0.0f64, curr_step.re(), &f, &mut convergency).unwrap(); - curr_step = T::lift(root); - let result = self.integrator.step(&self.ode, curr_step); - new_y = result.y; - dense_option = result.dense; - callback.effect(&mut self.ode); + // Look for the earliest event in the step. Nothing is committed + // while looking, because an event function is asked about the state + // at the start of the step as well as at the end. + let mut earliest: Option<(usize, f64)> = None; + for (index, (callback, fired)) in + self.callbacks.iter().zip(last_fired.iter()).enumerate() + { + // A callback that fired at this step's start searches from a + // little way in instead, so that it cannot find the event it + // just handled again. + let search_start = if *fired == Some(self.ode.t.re()) { + curr_step * EVENT_PROBE_FRACTION + } else { + T::lift(0.0) + }; + // The event condition is a real-valued question, so it is + // asked in the reals — see this module's header. + let before = if search_start.re() == 0.0 { + callback.event(self.ode.t, &self.ode.y, &self.ode.params).re() + } else { + let probe = self.integrator.step(&self.ode, search_start).y; + callback + .event(self.ode.t + search_start, &probe, &self.ode.params) + .re() + }; + let after = callback + .event(self.ode.t + curr_step, &new_y, &self.ode.params) + .re(); + let crossing = if before == 0.0 { + // Sitting exactly on the event. At the start of a step that + // is a crossing, and it fires there; a little way into one + // it is an event function this step can do nothing with, + // because the effect has already run at its zero. + (search_start.re() == 0.0).then_some(0.0) + } else if after == 0.0 { + Some(curr_step.re()) + } else if before * after < 0.0 { + let f = |test_h: f64| { + let h = T::lift(test_h); + let test_y = self.integrator.step(&self.ode, h).y; + callback.event(self.ode.t + h, &test_y, &self.ode.params).re() + }; + Some( + find_root_brent(search_start.re(), curr_step.re(), &f, &mut convergency) + .unwrap(), + ) + } else { + None + }; + // Two events in one step: the earlier one fires and the step + // ends there, so the later one is still ahead and is found + // again on the next step. A tie goes to registration order. + if let Some(h_event) = crossing { + if earliest.is_none_or(|(_, h)| h_event < h) { + earliest = Some((index, h_event)); } } } + + if let Some((index, h_event)) = earliest { + // Land on the crossing and record it, before the effect runs: + // the trajectory does reach this state, and a terminal effect + // should leave it as the last point of the solution. + if h_event != 0.0 { + let result = self.integrator.step(&self.ode, T::lift(h_event)); + self.ode.t += T::lift(h_event); + self.ode.y = result.y; + times.push(self.ode.t); + states.push(self.ode.y); + dense_coefficients.push(RefCell::new(result.dense.unwrap())); + } + last_fired[index] = Some(self.ode.t.re()); + let before_effect = self.ode.y; + self.callbacks[index].effect(&mut self.ode); + if self.ode.y != before_effect { + // The state jumped. Both sides of it are part of the + // solution, so both are recorded at the same time, with a + // zero-length interval between them. The controller's error + // history describes the solution up to the jump and says + // nothing about the one after it. + times.push(self.ode.t); + states.push(self.ode.y); + dense_coefficients.push(RefCell::new(vec![before_effect, self.ode.y])); + self.controller.reset(); + } + continue; + } + self.ode.y = new_y; self.ode.t += curr_step; times.push(self.ode.t); @@ -127,12 +235,30 @@ where } } - pub fn with_callback(mut self, callback: &'a dyn Callback) -> Self { + /// Registers a callback. + /// + /// Callbacks are watched in the order they are registered, which decides + /// which one fires when two events land at the same time. + pub fn with_callback(mut self, callback: &'a mut dyn Callback) -> Self { self.callbacks.push(callback); self } } +/// What an integration produced: the steps it took, and the means to evaluate +/// between them. +/// +/// `times` and `states` run in step order and are the same length. `dense` has +/// one fewer entry, the coefficients for interpolating inside the step from +/// `times[i]` to `times[i + 1]`. +/// +/// A callback effect that changes the state puts the same time in `times` +/// twice, with the state before the jump and the state after it. The pair is +/// what the trajectory does: it reaches one value and leaves from the other. +/// [`interpolate`](Solution::interpolate) at exactly that time reports the +/// second of them, the state the effect left, since that is the one the +/// trajectory carries forward; the value it arrived with is the entry before it +/// in `states`. pub struct Solution<'a, T: Scalar, S, const D: usize, P> where S: Integrator, @@ -148,6 +274,10 @@ impl Solution<'_, T, S, D, P> where S: Integrator, { + /// The state at `t`, interpolated inside whichever step contains it. + /// + /// At the time of a jump this is the state the callback's effect left. Panics + /// if `t` is outside the solution. pub fn interpolate(&self, t: T) -> SVector { // First check that the t is within bounds let last = self.times.last().unwrap().re(); @@ -164,7 +294,16 @@ where // Then find the two t values closest to the desired t match times.binary_search_by(|x| x.total_cmp(&t.re())) { - Ok(index) => self.states[index], + Ok(index) => { + // A jump records its time twice and the search may land on + // either entry. Walk to the last of them, which is the state the + // effect left and the one the trajectory goes on with. + let mut index = index; + while index + 1 < times.len() && times[index + 1] == times[index] { + index += 1; + } + self.states[index] + } Err(end_index) => { // Then send that to the integrator let t_start = self.times[end_index - 1]; @@ -172,6 +311,14 @@ where let y_start = self.states[end_index - 1]; let h = t_end - t_start; + if h.re() == 0.0 { + // The two sides of a jump. Any `t` strictly inside this + // interval would have to be the jump's own time, which the + // search above answers exactly, so this is only reached if + // a step of zero length was recorded some other way. + return y_start; + } + // Check if we need to compute extra stages for lazy dense output let dense_cell = &self.dense[end_index - 1]; @@ -221,7 +368,16 @@ mod tests { use crate::integrator::dormand_prince::DormandPrince45; use crate::integrator::Tolerances; use approx::assert_relative_eq; - use nalgebra::Vector3; + use nalgebra::{Matrix3, Vector1, Vector2, Vector3}; + use std::cell::Cell; + use std::f64::consts::PI; + use std::rc::Rc; + + /// `dy/dt = 1`, so the state is the time elapsed and an event on the state + /// is an event on the clock. + fn unit_rate(_t: f64, _y: Vector1, _p: &()) -> Vector1 { + Vector1::new(1.0) + } #[test] fn test_problem() { @@ -259,12 +415,12 @@ mod tests { let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-5); let controller = PIController::default(); - let value_too_high = FnCallback { + let mut value_too_high = FnCallback { event: &|_: f64, y: &SVector, _: &Params| 10.0 - y[0], effect: &stop, }; - let mut problem = Problem::new(ode, dp45, controller).with_callback(&value_too_high); + let mut problem = Problem::new(ode, dp45, controller).with_callback(&mut value_too_high); let solution = problem.solve(); assert_relative_eq!( @@ -274,6 +430,410 @@ mod tests { ); } + /// An effect that writes to `ode.y` is what the integration carries on from. + #[test] + fn test_effect_changes_the_state() { + let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ()); + let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-10); + let mut jump = FnCallback { + event: &|_t: f64, y: &Vector1, _p: &()| 5.0 - y[0], + effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0), + }; + + let mut problem = + Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump); + let solution = problem.solve(); + + // Five units of integration before the jump, five after. + assert_relative_eq!(solution.times.last().unwrap(), &10.0, max_relative = 1e-12); + assert_relative_eq!( + solution.states.last().unwrap()[0], + 105.0, + max_relative = 1e-9 + ); + } + + /// A jump is recorded as two states at the same time, in the order the + /// trajectory visits them. + #[test] + fn test_jump_records_both_sides() { + let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ()); + let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-10); + let mut jump = FnCallback { + event: &|_t: f64, y: &Vector1, _p: &()| 5.0 - y[0], + effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0), + }; + + let mut problem = + Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump); + let solution = problem.solve(); + + let at_crossing: Vec = (0..solution.times.len()) + .filter(|&i| (solution.times[i] - 5.0).abs() < 1e-9) + .collect(); + assert_eq!(at_crossing.len(), 2, "the two sides of the jump"); + let (left, right) = (at_crossing[0], at_crossing[1]); + assert_eq!(right, left + 1); + assert_relative_eq!(solution.states[left][0], 5.0, max_relative = 1e-9); + assert_relative_eq!(solution.states[right][0], 100.0, max_relative = 1e-12); + + // Interpolation either side of the jump follows the branch it is on, + // and at the jump itself reports the state the effect left. + assert_relative_eq!(solution.interpolate(4.9)[0], 4.9, max_relative = 1e-9); + assert_relative_eq!(solution.interpolate(5.1)[0], 100.1, max_relative = 1e-9); + assert_relative_eq!( + solution.interpolate(solution.times[left])[0], + 100.0, + max_relative = 1e-12 + ); + } + + /// `stop` ends the solution at the crossing, and leaves `t_end` alone. + #[test] + fn test_stop_ends_at_the_crossing() { + let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ()); + let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-10); + let mut stopper = FnCallback { + event: &|_t: f64, y: &Vector1, _p: &()| 5.0 - y[0], + effect: &stop, + }; + + let mut problem = + Problem::new(ode, dp45, PIController::default()).with_callback(&mut stopper); + let solution = problem.solve(); + + assert_relative_eq!(solution.times.last().unwrap(), &5.0, max_relative = 1e-12); + assert_relative_eq!( + solution.states.last().unwrap()[0], + 5.0, + max_relative = 1e-12 + ); + // Nothing past the event, and the problem's own end time is unchanged. + assert!(solution.times.iter().all(|&t| t <= 5.0 + 1e-12)); + assert_eq!(solution.ode.t_end, 10.0); + assert!(solution.ode.finished); + } + + /// `effect` takes `&mut self`, so a callback can count its own firings and + /// act only on the first one. + #[test] + fn test_callback_can_remember_firing() { + struct FireOnce { + fired: usize, + } + + impl Callback for FireOnce { + fn event(&self, t: f64, _y: &Vector1, _p: &()) -> f64 { + // Crosses at t = 1, 3, 5 and 7. + (t * PI / 2.0).cos() + } + + fn effect(&mut self, ode: &mut ODE<'_, f64, 1, ()>) { + self.fired += 1; + if self.fired == 1 { + ode.y = Vector1::new(100.0); + } + } + } + + let ode = ODE::new(&unit_rate, 0.0, 9.0, Vector1::new(0.0), ()); + let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-10); + let mut once = FireOnce { fired: 0 }; + // `dy/dt = 1` has no local error at all, so the step size runs away to + // `h_max`; capped at half a period, since a step long enough to contain + // two crossings sees no sign change and finds neither. + let controller = PIController::new(0.17, 0.04, 10.0, 0.2, 1.0, 0.9, 1e-4); + + let mut problem = Problem::new(ode, dp45, controller).with_callback(&mut once); + let solution = problem.solve(); + let last = solution.states.last().unwrap()[0]; + + assert_eq!(once.fired, 4); + // Only the first firing moved the state: 1 unit before the jump, 8 after. + assert_relative_eq!(last, 108.0, max_relative = 1e-9); + } + + /// An event function that is already zero at the start of a step crosses + /// there, and fires there, rather than being passed over. + #[test] + fn test_event_already_zero_fires_at_the_step_start() { + // y starts at zero and the event is y, so the event is met at t = 0. + let ode = ODE::new(&unit_rate, 0.0, 1.0, Vector1::new(0.0), ()); + let mut at_zero = FnCallback { + event: &|_t: f64, y: &Vector1, _p: &()| y[0], + effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0), + }; + + let mut problem = Problem::new(ode, DormandPrince45::new(), PIController::default()) + .with_callback(&mut at_zero); + let solution = problem.solve(); + + assert_eq!(solution.times[0], 0.0); + assert_eq!(solution.states[0][0], 0.0); + assert_eq!(solution.times[1], 0.0); + assert_eq!(solution.states[1][0], 100.0); + assert_relative_eq!( + solution.states.last().unwrap()[0], + 101.0, + max_relative = 1e-9 + ); + } + + /// Two events inside one step: the earlier one fires, and the later one is + /// found again on the step after it. Which callback was registered first + /// does not change the order they fire in. + #[test] + fn test_earliest_event_in_a_step_fires_first() { + fn order_of_firings(three_first: bool) -> Vec { + let fired = Rc::new(Cell::new([0.0f64; 2])); + + let at_three = { + let fired = Rc::clone(&fired); + move |ode: &mut ODE<'_, f64, 1, ()>| { + let mut seen = fired.get(); + seen[0] = ode.t; + fired.set(seen); + } + }; + let at_five = { + let fired = Rc::clone(&fired); + move |ode: &mut ODE<'_, f64, 1, ()>| { + let mut seen = fired.get(); + seen[1] = ode.t; + fired.set(seen); + } + }; + + let mut three = FnCallback { + event: &|t: f64, _y: &Vector1, _p: &()| 3.0 - t, + effect: &at_three, + }; + let mut five = FnCallback { + event: &|t: f64, _y: &Vector1, _p: &()| 5.0 - t, + effect: &at_five, + }; + + let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ()); + // One step of 10 would swallow both events whole. + let controller = PIController::new(0.17, 0.04, 10.0, 0.2, 100000.0, 0.9, 10.0); + let mut problem = if three_first { + Problem::new(ode, DormandPrince45::new(), controller) + .with_callback(&mut three) + .with_callback(&mut five) + } else { + Problem::new(ode, DormandPrince45::new(), controller) + .with_callback(&mut five) + .with_callback(&mut three) + }; + problem.solve(); + fired.get().to_vec() + } + + for three_first in [true, false] { + let fired = order_of_firings(three_first); + assert_relative_eq!(fired[0], 3.0, max_relative = 1e-9); + assert_relative_eq!(fired[1], 5.0, max_relative = 1e-9); + } + } + + /// An effect that leaves the event function sitting on zero does not + /// re-trigger the callback at the same point. + #[test] + fn test_effect_landing_on_its_own_event_terminates() { + // The event is |y| - 1 and the effect negates y, so the event value is + // unchanged by the effect: the naive reading is that it has crossed + // again and should fire again, at no advance in time. + let ode = ODE::new(&unit_rate, 0.0, 3.0, Vector1::new(0.0), ()); + let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-10); + let mut flip = FnCallback { + event: &|_t: f64, y: &Vector1, _p: &()| y[0].abs() - 1.0, + effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = -ode.y, + }; + + let mut problem = + Problem::new(ode, dp45, PIController::default()).with_callback(&mut flip); + let solution = problem.solve(); + + // y rises to 1 at t = 1 and flips to -1, which the effect leaves on the + // event; it rises back to 1 at t = 3 and flips again. Two jumps, neither + // of them repeated at the time it happened. + let jumps: Vec = (1..solution.times.len()) + .filter(|&i| solution.times[i] == solution.times[i - 1]) + .map(|i| solution.times[i]) + .collect(); + assert_eq!(jumps.len(), 2); + assert_relative_eq!(jumps[0], 1.0, max_relative = 1e-9); + assert_relative_eq!(jumps[1], 3.0, max_relative = 1e-9); + assert_relative_eq!(solution.times.last().unwrap(), &3.0, max_relative = 1e-12); + assert_relative_eq!(solution.states.last().unwrap()[0], -1.0, max_relative = 1e-8); + } + + /// The controller's error history is discarded when an effect changes the + /// state, and left alone when it does not. + #[test] + fn test_controller_is_reset_across_a_jump() { + /// Counts the resets `solve` asks for; a `PIController` otherwise. + struct Spy { + inner: PIController, + resets: Rc>, + } + + impl Controller for Spy { + fn determine_step(&mut self, h: f64, err: f64) -> TryStep { + self.inner.determine_step(h, err) + } + + fn next_step_guess(&self) -> TryStep { + Controller::next_step_guess(&self.inner) + } + + fn set_next_step_guess(&mut self, guess: TryStep) { + self.inner.set_next_step_guess(guess) + } + + fn reset(&mut self) { + self.resets.set(self.resets.get() + 1); + Controller::reset(&mut self.inner); + } + } + + fn resets_during(effect: &dyn Fn(&mut ODE<'_, f64, 1, ()>)) -> usize { + let resets = Rc::new(Cell::new(0)); + let spy = Spy { + inner: PIController::default(), + resets: Rc::clone(&resets), + }; + let mut callback = FnCallback { + event: &|_t: f64, y: &Vector1, _p: &()| 5.0 - y[0], + effect, + }; + let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ()); + let mut problem = Problem::new(ode, DormandPrince45::new(), spy) + .with_callback(&mut callback); + problem.solve(); + resets.get() + } + + assert_eq!( + resets_during(&|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0)), + 1 + ); + // An effect that changes nothing leaves a smooth solution behind it. + assert_eq!(resets_during(&|_ode: &mut ODE<'_, f64, 1, ()>| {}), 0); + } + + /// A bouncing ball, which is the standard exercise of repeated jumps: an + /// event every time the ball reaches the floor, and an effect that reverses + /// and shrinks its velocity. + /// + /// The bounce times are known. A ball dropped from `h0` first lands at + /// `sqrt(2 h0 / g)` with speed `v = g t`, and each bounce after that is a + /// flight of `2 v / g` at the restitution-scaled speed. + #[test] + fn test_bouncing_ball() { + const G: f64 = 9.81; + const RESTITUTION: f64 = 0.8; + + fn free_fall(_t: f64, y: Vector2, g: &f64) -> Vector2 { + Vector2::new(y[1], -g) + } + + let ode = ODE::new(&free_fall, 0.0, 2.0, Vector2::new(1.0, 0.0), G); + let vern7 = crate::integrator::vern7::Vern7::new().a_tol(1e-12).r_tol(1e-12); + let mut bounce = FnCallback { + event: &|_t: f64, y: &Vector2, _p: &f64| y[0], + effect: &|ode: &mut ODE<'_, f64, 2, f64>| { + ode.y = Vector2::new(0.0, -RESTITUTION * ode.y[1]); + }, + }; + + let mut problem = + Problem::new(ode, vern7, PIController::default()).with_callback(&mut bounce); + let solution = problem.solve(); + + // Each bounce is a jump, so it is a repeated time in the solution. + let bounces: Vec<(f64, f64)> = (1..solution.times.len()) + .filter(|&i| solution.times[i] == solution.times[i - 1]) + .map(|i| (solution.times[i], solution.states[i][1])) + .collect(); + + // Three of them land before t = 2, and the fourth does not. + let mut expected_t = (2.0 / G).sqrt(); + let mut expected_v = RESTITUTION * G * expected_t; + assert_eq!(bounces.len(), 3); + for (t, v) in &bounces { + assert_relative_eq!(*t, expected_t, max_relative = 1e-9); + assert_relative_eq!(*v, expected_v, max_relative = 1e-9); + expected_t += 2.0 * expected_v / G; + expected_v *= RESTITUTION; + } + + // Interpolating inside an arc gives free fall from the bounce that + // started it, so the dense output survived the jumps. + for &(t, v) in &bounces[..2] { + for dt in [1e-3, 0.05, 0.1] { + assert_relative_eq!( + solution.interpolate(t + dt)[0], + v * dt - 0.5 * G * dt * dt, + epsilon = 1e-9 + ); + } + } + } + + /// The case the whole feature is for: a rotating rigid body parameterized by + /// Modified Rodrigues Parameters, whose parameters are swapped for the + /// shadow set before they run into the singularity at a full turn. + /// + /// For a constant angular velocity about one axis the answer is known. + /// Starting from zero, the MRP vector is `tan(theta / 4)` along the rotation + /// axis, with `theta = |omega| t`, so it reaches magnitude one at half a + /// turn. The shadow set `-sigma / |sigma|^2` describes the same attitude and + /// obeys the same kinematics, and past half a turn it is the branch that + /// stays bounded: `-cot(theta / 4)` along the axis. + #[test] + fn test_mrp_shadow_set_swap() { + /// MRP kinematics, `sigma' = B(sigma) omega / 4`. + fn mrp_rate(_t: f64, s: Vector3, omega: &Vector3) -> Vector3 { + let skew = Matrix3::new(0.0, -s[2], s[1], s[2], 0.0, -s[0], -s[1], s[0], 0.0); + let b = Matrix3::identity() * (1.0 - s.dot(&s)) + skew * 2.0 + s * s.transpose() * 2.0; + b * omega * 0.25 + } + + let omega = Vector3::new(0.0, 0.0, 1.0); + let t_end = 5.0; // past half a turn, short of a full one + let ode = ODE::new(&mrp_rate, 0.0, t_end, Vector3::zeros(), omega); + let vern7 = crate::integrator::vern7::Vern7::new().a_tol(1e-12).r_tol(1e-12); + + let mut shadow = FnCallback { + event: &|_t: f64, s: &Vector3, _p: &Vector3| 1.0 - s.norm(), + effect: &|ode: &mut ODE<'_, f64, 3, Vector3>| { + ode.y = -ode.y / ode.y.norm_squared(); + }, + }; + let mut problem = + Problem::new(ode, vern7, PIController::default()).with_callback(&mut shadow); + let solution = problem.solve(); + + // The swap happened at half a turn, and only once. A jump is the one + // thing that puts the same time in the solution twice. + let jumps: Vec = (1..solution.times.len()) + .filter(|&i| solution.times[i] == solution.times[i - 1]) + .map(|i| solution.times[i]) + .collect(); + assert_eq!(jumps.len(), 1); + assert_relative_eq!(jumps[0], PI, max_relative = 1e-9); + + // The parameters stayed away from the singularity throughout. + assert!(solution.states.iter().all(|s| s.norm() <= 1.0 + 1e-9)); + + // And the attitude at the end is the one the body is actually in. + let expected = -Vector3::new(0.0, 0.0, 1.0 / (t_end / 4.0).tan()); + for i in 0..3 { + assert_relative_eq!(solution.states.last().unwrap()[i], expected[i], epsilon = 1e-9); + } + } + #[test] fn test_with_interpolation() { type Params = ();