Implement constants module to store shared numerical values

Author: srenevey

src/constants.rs

@@ -0,0 +1,54 @@
+use crate::dop_shared::FloatNumber;
+
+/// Stiffness detection thresholds for different methods
+pub mod stiffness {
+    use super::FloatNumber;
+
+    /// Stiffness threshold for DOPRI5 method
+    /// Based on Hairer & Wanner, "Solving Ordinary Differential Equations II"
+    pub fn dopri5_threshold<T: FloatNumber>() -> T {
+        T::from(3.25).unwrap()
+    }
+
+    /// Stiffness threshold for DOP853 method
+    /// Based on Hairer & Wanner, "Solving Ordinary Differential Equations II"
+    pub fn dop853_threshold<T: FloatNumber>() -> T {
+        T::from(6.1).unwrap()
+    }
+
+    /// Maximum consecutive stiffness detections before error
+    pub const MAX_STIFF_ITERATIONS: u32 = 15;
+
+    /// Number of non-stiff steps needed to reset stiffness counter
+    pub const NON_STIFF_RESET_COUNT: u32 = 6;
+}
+
+/// Initial step size computation constants
+pub mod initial_step {
+    use super::FloatNumber;
+
+    /// Minimum tolerance for initial step computation
+    pub fn min_tolerance<T: FloatNumber>() -> T {
+        T::from(1.0e-10).unwrap()
+    }
+
+    /// Default initial step when tolerance conditions not met
+    pub fn default_initial_step<T: FloatNumber>() -> T {
+        T::from(1.0e-6).unwrap()
+    }
+
+    /// Safety factor for initial step estimation
+    pub fn safety_factor<T: FloatNumber>() -> T {
+        T::from(0.01).unwrap()
+    }
+}
+
+/// Dense output interpolation constants
+pub mod dense_output {
+    use super::FloatNumber;
+
+    /// Floating point tolerance for endpoint detection
+    pub fn endpoint_tolerance<T: FloatNumber>() -> T {
+        T::from(1.0e-9).unwrap()
+    }
+}

src/dop853.rs

@@ -1,6 +1,7 @@
 //! Explicit Runge-Kutta method with Dormand-Prince coefficients of order 8(5,3) and dense output of order 7.
 
 use crate::butcher_tableau::dopri853;
+use crate::constants::{dense_output, initial_step, stiffness};
 use crate::controller::Controller;
 use crate::dop_shared::*;
 
@@ -210,11 +211,11 @@ where
         }
 
         // Compute h0
-        let tol = T::from(1.0E-10).unwrap();
-        let mut h0 = if d0 < tol || d1 < tol {
-            T::from(1.0E-6).unwrap()
+        let tol = initial_step::min_tolerance();
+        let mut h0: T = if d0 < tol || d1 < tol {
+            initial_step::default_initial_step()
         } else {
-            T::from(0.01).unwrap() * (d0 / d1).sqrt()
+            initial_step::safety_factor::<T>() * (d0 / d1).sqrt()
         };
 
         h0 = h0.min(self.controller.h_max());
@@ -377,18 +378,18 @@ where
                         T::zero()
                     };
 
-                    if h_lamb > T::from(6.1).unwrap() {
+                    if h_lamb > stiffness::dop853_threshold() {
                         non_stiff = 0;
                         iasti += 1;
-                        if iasti == 15 {
+                        if iasti == stiffness::MAX_STIFF_ITERATIONS {
                             self.h_old = self.h;
                             return Err(IntegrationError::StiffnessDetected {
                                 x: f64::from(self.x),
                             });
                         }
                     } else {
                         non_stiff += 1;
-                        if non_stiff == 6 {
+                        if non_stiff == stiffness::NON_STIFF_RESET_COUNT {
                             iasti = 0;
                         }
                     }
@@ -528,7 +529,7 @@ where
                 }
 
                 // 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() {
+                if (self.xd - self.x_end).abs() < dense_output::endpoint_tolerance() {
                     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);

src/dopri5.rs

@@ -1,6 +1,7 @@
 //! Explicit Runge-Kutta method with Dormand-Prince coefficients of order 5(4) and dense output of order 4.
 
 use crate::butcher_tableau::dopri54;
+use crate::constants::{dense_output, initial_step, stiffness};
 use crate::continuous_output_model::ContinuousOutputModel;
 use crate::controller::Controller;
 use crate::dop_shared::*;
@@ -202,11 +203,11 @@ where
         }
 
         // Compute h0
-        let tol = T::from(1.0E-10).unwrap();
+        let tol = initial_step::min_tolerance();
         let mut h0 = if d0 < tol || d1 < tol {
-            T::from(1.0E-6).unwrap()
+            initial_step::default_initial_step()
         } else {
-            T::from(0.01).unwrap() * (d0 / d1).sqrt()
+            initial_step::safety_factor::<T>() * (d0 / d1).sqrt()
         };
 
         h0 = h0.min(self.controller.h_max());
@@ -383,18 +384,18 @@ where
                         T::zero()
                     };
 
-                    if h_lamb > T::from(3.25).unwrap() {
+                    if h_lamb > stiffness::dopri5_threshold() {
                         iasti += 1;
                         non_stiff = 0;
-                        if iasti == 15 {
+                        if iasti == stiffness::MAX_STIFF_ITERATIONS {
                             self.h_old = self.h;
                             return Err(IntegrationError::StiffnessDetected {
                                 x: f64::from(self.x),
                             });
                         }
                     } else {
                         non_stiff += 1;
-                        if non_stiff == 6 {
+                        if non_stiff == stiffness::NON_STIFF_RESET_COUNT {
                             iasti = 0;
                         }
                     }
@@ -468,7 +469,7 @@ where
             }
 
             // 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() {
+            if (self.xd - self.x_end).abs() < dense_output::endpoint_tolerance() {
                 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);

src/lib.rs

@@ -9,6 +9,7 @@ use nalgebra as na;
 
 // Declare modules
 pub mod butcher_tableau;
+mod constants;
 pub mod continuous_output_model;
 pub mod controller;
 pub mod dop853;