Add final point of the integration interval if within rounding error

- Updated the `solution_output` function to ensure that the last point is added if it's within floating point error of  x_end.
- Moved the `sign` function to dop_shared.
- Fixed typos.

Author: srenevey

src/dop853.rs

@@ -4,7 +4,7 @@ use crate::butcher_tableau::dopri853;
 use crate::controller::Controller;
 use crate::dop_shared::*;
 
-use nalgebra::{allocator::Allocator, DefaultAllocator, Dim, OVector};
+use nalgebra::{allocator::Allocator, DefaultAllocator, Dim, MatrixSum, OVector, U1};
 
 trait DefaultController<T: FloatNumber> {
     fn default(x: T, x_end: T) -> Self;
@@ -123,9 +123,9 @@ where
     /// * `fac_min` - Minimum factor between two successive steps. Default is 0.333
     /// * `fac_max` - Maximum factor between two successive steps. Default is 6.0
     /// * `h_max`   - Maximum step size. Default is `x_end-x
-    /// * `h`       - Initial value of the step size. If h = 0.0, the intial value of h is computed automatically
+    /// * `h`       - Initial value of the step size. If h = 0.0, the initial value of h is computed automatically
     /// * `n_max`   - Maximum number of iterations. Default is 100000
-    /// * `n_stiff` - Stifness is tested when the number of iterations is a multiple of n_stiff. Default is 1000
+    /// * `n_stiff` - Stiffness is tested when the number of iterations is a multiple of n_stiff. Default is 1000
     /// * `out_type`    - Type of the output. Must be a variant of the OutputType enum. Default is Dense
     ///
     #[allow(clippy::too_many_arguments)]
@@ -190,7 +190,7 @@ where
         }
     }
 
-    /// Compute the initial stepsize
+    /// Computes the initial step size.
     fn hinit(&self) -> T {
         let (rows, cols) = self.y.shape_generic();
         let mut f0 = OVector::zeros_generic(rows, cols);
@@ -367,7 +367,7 @@ where
                 self.f.system(self.x + self.h, &y_tmp, &mut k[3]);
                 self.stats.num_eval += 1;
 
-                // Stifness detection
+                // Stiffness detection
                 if self.stats.accepted_steps % self.n_stiff == 0 || iasti > 0 {
                     let num = T::from((&k[3] - &k[2]).dot(&(&k[3] - &k[2]))).unwrap();
                     let den = T::from((&k[4] - &y_next).dot(&(&k[4] - &y_next))).unwrap();
@@ -521,30 +521,42 @@ where
                     if self.x_old.abs() <= self.xd.abs() && self.x.abs() >= self.xd.abs() {
                         let theta = (self.xd - self.x_old) / self.h_old;
                         let theta1 = T::one() - theta;
-
-                        let y_out = &self.rcont[0]
-                            + (&self.rcont[1]
-                                + (&self.rcont[2]
-                                    + (&self.rcont[3]
-                                        + (&self.rcont[4]
-                                            + (&self.rcont[5]
-                                                + (&self.rcont[6] + &self.rcont[7] * theta)
-                                                    * theta1)
-                                                * theta)
-                                            * theta1)
-                                        * theta)
-                                    * theta1)
-                                * theta;
+                        let y_out = self.compute_y_out(theta, theta1);
                         self.results.push(self.xd, y_out);
                         self.xd += self.dx;
                     }
                 }
+
+                // Ensure the last point is added if it's within floating point error of x_end.
+                if (self.xd - self.x_end).abs() < T::from(1e-9).unwrap() {
+                    let theta = (self.x_end - self.x_old) / self.h_old;
+                    let theta1 = T::one() - theta;
+                    let y_out = self.compute_y_out(theta, theta1);
+                    self.results.push(self.x_end, y_out);
+                    self.xd += self.dx;
+                }
             }
         } else {
-            self.results.push(self.x0, y_next);
+            self.results.push(self.x0, y_next.clone());
         }
     }
 
+    /// Computes the value of y for given theta and theta1 values.
+    fn compute_y_out(&mut self, theta: T, theta1: T) -> MatrixSum<T, D, U1, D, U1> {
+        &self.rcont[0]
+            + (&self.rcont[1]
+                + (&self.rcont[2]
+                    + (&self.rcont[3]
+                        + (&self.rcont[4]
+                            + (&self.rcont[5]
+                                + (&self.rcont[6] + &self.rcont[7] * theta) * theta1)
+                                * theta)
+                            * theta1)
+                        * theta)
+                    * theta1)
+                * theta
+    }
+
     /// Getter for the independent variable's output.
     pub fn x_out(&self) -> &Vec<T> {
         self.results.get().0
@@ -572,14 +584,6 @@ where
     }
 }
 
-fn sign<T: FloatNumber>(a: T, b: T) -> T {
-    if b > T::zero() {
-        a.abs()
-    } else {
-        -a.abs()
-    }
-}
-
 #[cfg(test)]
 mod tests {
     use super::*;

src/dop_shared.rs

@@ -14,7 +14,7 @@ use thiserror::Error;
 /// The type parameter T should be either `f32` or `f64`, the trait [FloatNumber] is used
 /// internally to allow generic code.
 ///
-/// The type parameter V is a state vector. To have an easy start it is recommend to use [nalgebra] vectors.
+/// The type parameter V is a state vector. To have an easy start it is recommended to use [nalgebra] vectors.
 ///
 /// ```
 /// use ode_solvers::{System, SVector, Vector3};
@@ -34,6 +34,7 @@ where
 {
     /// System of ordinary differential equations.
     fn system(&self, x: T, y: &V, dy: &mut V);
+
     /// Stop function called at every successful integration step. The integration is stopped when this function returns true.
     fn solout(&mut self, _x: T, _y: &V, _dy: &V) -> bool {
         false
@@ -141,12 +142,6 @@ impl Stats {
             rejected_steps: 0,
         }
     }
-
-    /// Prints some statistics related to the integration process.
-    #[deprecated(since = "0.2.0", note = "Use std::fmt::Display instead")]
-    pub fn print(&self) {
-        println!("{}", self);
-    }
 }
 
 impl fmt::Display for Stats {
@@ -156,3 +151,11 @@ impl fmt::Display for Stats {
         write!(f, "Number of rejected steps: {}", self.rejected_steps)
     }
 }
+
+pub(crate) fn sign<T: FloatNumber>(a: T, b: T) -> T {
+    if b > T::zero() {
+        a.abs()
+    } else {
+        -a.abs()
+    }
+}