Add continuous output type, fix unit test, add integration test

Author: srenevey

Cargo.lock

@@ -1,6 +1,6 @@
 [package]
 name = "ode_solvers"
-version = "0.6.0"
+version = "0.6.1"
 authors = ["Sylvain Renevey <[email protected]>"]
 description = "Numerical methods to solve ordinary differential equations (ODEs) in Rust."
 edition = "2021"

Cargo.toml

@@ -14,7 +14,7 @@ To start using the crate in a project, the following dependency must be added in
 
 ```rust
 [dependencies]
-ode_solvers = "0.6.0"
+ode_solvers = "0.6.1"
 ```
 
 Then, in the main file, add

README.md

@@ -75,12 +75,6 @@ where
         coefficients: Vec<OVector<T, D>>,
         step_size: T,
     ) {
-        debug_assert!(
-            self.breakpoints
-                .last()
-                .map_or(true, |&last| breakpoint > last),
-            "Breakpoints must be added in ascending order"
-        );
         self.breakpoints.push(breakpoint);
         self.intervals.push(Interval {
             coefficients,

src/continuous_output_model.rs

@@ -113,6 +113,7 @@ impl<T, V> Default for SolverResult<T, V> {
 pub enum OutputType {
     Dense,
     Sparse,
+    Continuous,
 }
 
 /// Enumeration of the errors that may arise during integration.

src/dop_shared.rs

@@ -246,11 +246,15 @@ where
         &mut self,
         continuous_output: &mut ContinuousOutputModel<T, OVector<T, D>>,
     ) -> Result<Stats, IntegrationError> {
+        self.out_type = OutputType::Continuous;
         self.integrate_core(Some(continuous_output))
     }
 
     /// Integrates the system.
     pub fn integrate(&mut self) -> Result<Stats, IntegrationError> {
+        if self.out_type == OutputType::Continuous {
+            panic!("Please use `integrate_with_continuous_output_model` to compute a continuous output model.");
+        }
         self.integrate_core(None)
     }
 
@@ -334,8 +338,8 @@ where
             k[1] = k[6].clone();
             self.stats.num_eval += 6;
 
-            // Prepare dense output
-            if self.out_type == OutputType::Dense {
+            // Prepare dense/continuous output
+            if self.out_type == OutputType::Dense || self.out_type == OutputType::Continuous {
                 self.rcont[4] = (&k[0] * dopri54::d::<T>(1)
                     + &k[2] * dopri54::d::<T>(3)
                     + &k[3] * dopri54::d::<T>(4)
@@ -396,8 +400,8 @@ where
                     }
                 }
 
-                // Prepare dense output
-                if self.out_type == OutputType::Dense {
+                // Prepare dense/continuous output
+                if self.out_type == OutputType::Dense || self.out_type == OutputType::Continuous {
                     let h = self.h;
 
                     let ydiff = &y_next - &self.y;
@@ -414,7 +418,7 @@ where
                 self.x += self.h;
                 self.h_old = self.h;
 
-                self.solution_output(y_next, &mut continuous_output_model);
+                self.solution_output(y_next, &mut continuous_output_model, &k[0]);
 
                 if self
                     .f
@@ -443,6 +447,7 @@ where
         &mut self,
         y_next: OVector<T, D>,
         continuous_output_model: &mut Option<&mut ContinuousOutputModel<T, OVector<T, D>>>,
+        dy: &OVector<T, D>,
     ) {
         if self.out_type == OutputType::Dense {
             while self.xd.abs() <= self.x.abs() {
@@ -453,6 +458,13 @@ where
                     self.results.push(self.xd, y_out);
                     self.xd += self.dx;
                 }
+
+                if self
+                    .f
+                    .solout(self.xd, self.results.get().1.last().unwrap(), dy)
+                {
+                    break;
+                }
             }
 
             // Ensure the last point is added if it's within floating point error of x_end.
@@ -514,7 +526,7 @@ mod tests {
     use crate::{OVector, System, Vector1};
     use nalgebra::{allocator::Allocator, DefaultAllocator, Dim};
 
-    // Same as Test3 from rk4.rs, but aborts after x is greater/equal than 0.5
+    // Same as Test3 from rk4.rs, but aborts after x is equal to or greater than 0.5
     struct Test1 {}
     impl<D: Dim> System<f64, OVector<f64, D>> for Test1
     where
@@ -536,9 +548,34 @@ mod tests {
         let _ = stepper.integrate();
 
         let x = stepper.x_out();
-        assert!((*x.last().unwrap() - 0.5).abs() < 1.0E-9); //
+        assert!((*x.last().unwrap() - 0.5).abs() < 1.0E-9);
 
         let out = stepper.y_out();
-        assert!((&out[5][0] - 0.913059243).abs() < 1.0E-9);
+        assert!((&out[5][0] - 0.9130611474392001).abs() < 1.0E-9);
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_integrate_when_continuous_output_type_panic() {
+        let system = Test1 {};
+        let mut stepper = Dopri5::from_param(
+            system,
+            0.,
+            1.,
+            0.1,
+            Vector1::new(1.),
+            1e-12,
+            1e-6,
+            0.1,
+            0.2,
+            0.3,
+            2.0,
+            5.0,
+            0.1,
+            10000,
+            1000,
+            OutputType::Continuous,
+        );
+        let _ = stepper.integrate();
     }
 }

src/dopri5.rs

@@ -3,6 +3,9 @@ use nalgebra::Vector1;
 use ode_solvers::continuous_output_model::ContinuousOutputModel;
 use ode_solvers::{Dopri5, System};
 
+const INTEGRATION_TOLERANCE: f64 = 1e-10;
+const STEP_SIZE: f64 = 0.1;
+
 type State = Vector1<f64>;
 type Time = f64;
 
@@ -15,13 +18,33 @@ impl System<f64, State> for Ode {
     }
 }
 
+struct OdeWithStoppingCondition;
+
+impl System<f64, State> for OdeWithStoppingCondition {
+    fn system(&self, t: Time, y: &State, dy: &mut State) {
+        dy[0] = -3.0 * t * t * y[0].exp();
+    }
+
+    fn solout(&mut self, x: f64, _y: &State, _dy: &State) -> bool {
+        x >= 4.1
+    }
+}
+
 #[test]
 fn test_dopri5() {
     let f = Ode {};
     let initial_value = State::new(0.0);
     let x_start = 0.0;
     let x_end = 6.0;
-    let mut stepper = Dopri5::new(f, x_start, x_end, 0.1, initial_value, 1e-10, 1e-10);
+    let mut stepper = Dopri5::new(
+        f,
+        x_start,
+        x_end,
+        STEP_SIZE,
+        initial_value,
+        INTEGRATION_TOLERANCE,
+        INTEGRATION_TOLERANCE,
+    );
     let _ = stepper.integrate();
 
     let (times_s, states) = stepper.results().get();
@@ -43,7 +66,15 @@ fn test_dopri5_with_continuous_output_model() {
     let initial_value = State::new(0.0);
     let x_start = 0.0;
     let x_end = 6.0;
-    let mut stepper = Dopri5::new(f, x_start, x_end, 0.1, initial_value, 1e-10, 1e-10);
+    let mut stepper = Dopri5::new(
+        f,
+        x_start,
+        x_end,
+        STEP_SIZE,
+        initial_value,
+        INTEGRATION_TOLERANCE,
+        INTEGRATION_TOLERANCE,
+    );
     let mut continuous_output_model = ContinuousOutputModel::default();
     let _ = stepper.integrate_with_continuous_output_model(&mut continuous_output_model);
 
@@ -63,6 +94,42 @@ fn test_dopri5_with_continuous_output_model() {
     }
 }
 
+#[test]
+fn test_dopri5_with_continuous_output_model_and_stopping_condition() {
+    let f = OdeWithStoppingCondition {};
+    let initial_value = State::new(0.0);
+    let x_start = 0.0;
+    let x_end = 6.0;
+    let mut stepper = Dopri5::new(
+        f,
+        x_start,
+        x_end,
+        STEP_SIZE,
+        initial_value,
+        INTEGRATION_TOLERANCE,
+        INTEGRATION_TOLERANCE,
+    );
+    let mut continuous_output_model = ContinuousOutputModel::default();
+    let _ = stepper.integrate_with_continuous_output_model(&mut continuous_output_model);
+
+    assert!(continuous_output_model.evaluate(-0.0001).is_none());
+    assert!(continuous_output_model.evaluate(0.0).is_some());
+    assert!(continuous_output_model.evaluate(4.1).is_some());
+    assert!(continuous_output_model.evaluate(4.1 + STEP_SIZE).is_none());
+    assert!(continuous_output_model.bounds().1 >= 4.1);
+    assert!(continuous_output_model.bounds().1 < 4.1 + STEP_SIZE);
+
+    let num_samples = 100;
+    let step_size = (x_end - 4.1) / num_samples as f64;
+    for i in 0..num_samples {
+        let t = i as f64 * step_size + x_start;
+        let continuous_output = continuous_output_model.evaluate(t);
+        let analytic = analytic_solution(t);
+        assert!(continuous_output.is_some());
+        assert_relative_eq!(continuous_output.unwrap(), analytic, epsilon = 1e-9);
+    }
+}
+
 /// Evaluates the analytic solution of the ODE at the given time.
 fn analytic_solution(t: Time) -> State {
     State::new(-(t.powi(3) + 1.0).ln())