Fix ODESolver stopping one step short and composite Simpson tail-drop

CHANGELOG.md

@@ -4,6 +4,12 @@ All notable changes to this project will be documented in this file, starting fr
 The format is based on [Keep a Changelog](http://keepachangelog.com/)
 and this project adheres to [Semantic Versioning](http://semver.org/).
 
+## [Unreleased]
+
+### Fixed
+- `ODESolver::solve` stopped one step short, returning `y(x_end - h)` and degrading every method (including the default RK4) to first-order accuracy; it now integrates to `x_end` exactly, even when the step size does not divide the interval.
+- The composite Simpson's 1/3 and 3/8 rules in `NewtonCotes` dropped the trailing subinterval(s) when the subdivision count was not a multiple of 2 or 3 (including the default of 1000 for Simpson's 3/8); the count is now rounded up to a valid multiple.
+
 ## [0.4.1] - 2026-06-14
 
 ### Changed

src/integrators/newton_cotes.rs

@@ -24,6 +24,15 @@ pub enum Formula {
     SimpsonsThreeEighths,
 }
 
+/// Number of subintervals spanned by one panel of the given composite rule.
+const fn panel_width(formula: &Formula) -> usize {
+    match formula {
+        Formula::Trapezium => 1,
+        Formula::SimpsonsOneThird => 2,
+        Formula::SimpsonsThreeEighths => 3,
+    }
+}
+
 /// # Newton-Cotes
 ///
 /// A numerical integrator of functions `f: R -> R` based on the composite
@@ -48,6 +57,8 @@ pub struct NewtonCotes<T, F> {
     f: F,
     formula: fn(&Self, &[T], T) -> T,
     subdivisions: usize,
+    /// Number of subintervals spanned by one panel of the current composite rule.
+    panel: usize,
 }
 
 impl<T, F> NewtonCotes<T, F>
@@ -63,6 +74,7 @@ where
             f,
             subdivisions: DEFAULT_SUBDIVISIONS,
             formula: Self::simpsons_one_third,
+            panel: panel_width(&Formula::SimpsonsOneThird),
         }
     }
 
@@ -83,6 +95,7 @@ where
     /// assert!((integral - 1.631).abs() <= 1e-3);
     /// ```
     pub fn with_formula(&mut self, formula: Formula) -> &mut Self {
+        self.panel = panel_width(&formula);
         self.formula = match formula {
             Formula::Trapezium => Self::trapezium,
             Formula::SimpsonsOneThird => Self::simpsons_one_third,
@@ -93,6 +106,10 @@ where
 
     /// Specify the number of subintervals the integration interval is split in.
     ///
+    /// When the integrator is run, the number of subintervals will be adjusted to the next nearest
+    /// multiple that the given formula is expecting. In particular, since the default formula is
+    /// Simpson's 1/3 Method, the count in that case is adjusted to the next nearest even integer.
+    ///
     /// The default is 1000.
     ///
     /// ## Examples
@@ -128,8 +145,18 @@ where
             });
         }
 
-        let mut subdivision_values = vec![T::zero(); self.subdivisions + 1];
-        let Some(subdivisions_as_t) = T::from(self.subdivisions) else {
+        // The composite Simpson rules require the subdivision count to be a multiple of
+        // the panel width (2 for Simpson's 1/3, 3 for Simpson's 3/8). Round the requested
+        // count up to the next valid multiple so no trailing subinterval is dropped.
+        let remainder = self.subdivisions % self.panel;
+        let subdivisions = if remainder == 0 {
+            self.subdivisions
+        } else {
+            self.subdivisions + self.panel - remainder
+        };
+
+        let mut subdivision_values = vec![T::zero(); subdivisions + 1];
+        let Some(subdivisions_as_t) = T::from(subdivisions) else {
             return Err(SolverError::TypeConversionError);
         };
 
@@ -138,7 +165,7 @@ where
         for (i, item) in subdivision_values
             .iter_mut()
             .enumerate()
-            .take(self.subdivisions + 1)
+            .take(subdivisions + 1)
         {
             *item = (self.f)(from + T::from(i).unwrap() * delta);
         }

src/solvers/ode_solver.rs

@@ -67,8 +67,7 @@ pub struct ODESolver<T, V, F> {
     x0: T,
     y0: V,
     h: T,
-    half_h: T,
-    method: fn(&Self, T, V) -> V,
+    method: fn(&Self, T, V, T) -> V,
 }
 
 impl<T, V, F> ODESolver<T, V, F>
@@ -99,7 +98,6 @@ where
             x0,
             y0,
             h,
-            half_h: h / T::from(2_f64).unwrap(),
             method: Self::rk4_step,
         }
     }
@@ -144,18 +142,27 @@ where
     /// assert!((solution[0] - SOLUTION) <= 1e-3);
     /// ```
     pub fn solve(&self, x_end: T) -> SolverResult<V> {
-        let mut x = self.x0;
-        let mut y = self.y0;
-        let steps = T::to_usize(&((x_end - self.x0) / self.h)).unwrap_or(0);
+        let interval = x_end - self.x0;
+        let steps = T::to_usize(&(interval / self.h).round()).unwrap_or(0);
         if steps == 0 {
             return Err(SolverError::IncorrectInput {
                 details: "the number of steps should be positive",
             });
         }
 
-        for _ in 1..steps {
-            y = (self.method)(self, x, y);
-            x = x + self.h;
+        // Take a uniform step that lands exactly on `x_end`, even when the requested
+        // step size does not divide the interval evenly.
+        let h = if let Some(steps_casted) = T::from(steps) {
+            interval / steps_casted
+        } else {
+            return Err(SolverError::TypeConversionError);
+        };
+
+        let mut x = self.x0;
+        let mut y = self.y0;
+        for _ in 0..steps {
+            y = (self.method)(self, x, y, h);
+            x = x + h;
         }
 
         Ok(y)
@@ -174,16 +181,13 @@ where
     ///
     /// let mut solver = ODESolver::new(f, x0, y0, h);
     ///
-    /// # let solution = solver.solve(x_end).unwrap();
-    /// # assert!((solution - (-1_f64).exp()) > 1e-3); // Error too big!
     /// let solution = solver
     ///                 .with_step_size(0.001)
     ///                 .solve(x_end); // This changes solver's step size until changed again
     /// # assert!((solution.unwrap() - (-1_f64).exp()) <= 1e-3);
     /// ```
     pub fn with_step_size(&mut self, h: T) -> &mut Self {
         self.h = h;
-        self.half_h = h / T::from(2.).unwrap();
         self
     }
 
@@ -202,8 +206,6 @@ where
     /// let h = 0.1; // step size
     ///
     /// let mut solver = ODESolver::new(f, x0, y0, h);
-    /// # let solution = solver.solve(x_end).unwrap();
-    /// # assert!((solution - (-1_f64).exp()) > 1e-3); // Error too big!
     ///
     /// let solution = solver
     ///                 .with_steps(x_end, 1000)
@@ -212,7 +214,6 @@ where
     /// ```
     pub fn with_steps(&mut self, x_end: T, steps: usize) -> &mut Self {
         self.h = (x_end - self.x0) / T::from(steps).unwrap();
-        self.half_h = self.h / T::from(2.).unwrap();
         self
     }
 
@@ -256,21 +257,23 @@ where
 
     // === PRIVATE FUNCTIONS: A step in the different methods available ===
 
-    fn euler_step(&self, x: T, y: V) -> V {
-        y + (self.f)(x, y) * self.h
+    fn euler_step(&self, x: T, y: V, h: T) -> V {
+        y + (self.f)(x, y) * h
     }
 
-    fn heun_step(&self, x: T, y: V) -> V {
-        let y1 = y + (self.f)(x, y) * self.h;
-        y + (y1 + (self.f)(x + self.h, y1)) * self.half_h
+    fn heun_step(&self, x: T, y: V, h: T) -> V {
+        let half_h = h * T::from(0.5).unwrap();
+        let y1 = y + (self.f)(x, y) * h;
+        y + (y1 + (self.f)(x + h, y1)) * half_h
     }
 
-    fn rk4_step(&self, x: T, y: V) -> V {
+    fn rk4_step(&self, x: T, y: V, h: T) -> V {
+        let half_h = h * T::from(0.5).unwrap();
         let k1 = (self.f)(x, y);
-        let k2 = (self.f)(x + self.half_h, y + k1 * self.half_h);
-        let k3 = (self.f)(x + self.half_h, y + k2 * self.half_h);
-        let k4 = (self.f)(x + self.h, y + k3 * self.h);
+        let k2 = (self.f)(x + half_h, y + k1 * half_h);
+        let k3 = (self.f)(x + half_h, y + k2 * half_h);
+        let k4 = (self.f)(x + h, y + k3 * h);
 
-        y + (k1 + k2 + k2 + k3 + k3 + k4) * (self.h / T::from(6_f64).unwrap())
+        y + (k1 + k2 + k2 + k3 + k3 + k4) * (h / T::from(6_f64).unwrap())
     }
 }

tests/integrators.rs

@@ -50,7 +50,45 @@ fn simpsons_three_eighths() {
         .integrate(0., 1.)
         .unwrap();
 
-    assert!((result - APPROXIMATE_INTEGRAL_GAUSSIAN_0_TO_1).abs() <= 1e-3)
+    assert!((result - APPROXIMATE_INTEGRAL_GAUSSIAN_0_TO_1).abs() <= 1e-9)
+}
+
+#[test]
+fn simpsons_one_third_exact_for_cubics() {
+    // Simpson's 1/3 is exact for cubics; the composite rule must not drop the last
+    // subinterval when the subdivision count is odd (#15). ∫₀³ x³ dx = 20.25.
+    let f = |x: f64| x * x * x;
+    for n in [10usize, 11, 99, 100, 101, 1000, 1001] {
+        let result = NewtonCotes::new(f)
+            .with_formula(Formula::SimpsonsOneThird)
+            .with_subdivisions(n)
+            .integrate(0., 3.)
+            .unwrap();
+        assert!((result - 20.25).abs() < 1e-9, "n = {n}");
+    }
+}
+
+#[test]
+fn simpsons_three_eighths_exact_for_cubics() {
+    // Simpson's 3/8 is exact for cubics; the composite rule must not drop trailing
+    // subintervals when the count is not a multiple of three, including the default
+    // of 1000 (#15). ∫₀³ x³ dx = 20.25.
+    let f = |x: f64| x * x * x;
+    for n in [12usize, 13, 14, 100, 300, 301, 1000] {
+        let result = NewtonCotes::new(f)
+            .with_formula(Formula::SimpsonsThreeEighths)
+            .with_subdivisions(n)
+            .integrate(0., 3.)
+            .unwrap();
+        assert!((result - 20.25).abs() < 1e-9, "n = {n}");
+    }
+
+    // Default subdivisions (1000 is not a multiple of three) via the public default path.
+    let result = NewtonCotes::new(f)
+        .with_formula(Formula::SimpsonsThreeEighths)
+        .integrate(0., 3.)
+        .unwrap();
+    assert!((result - 20.25).abs() < 1e-9, "default subdivisions");
 }
 
 #[test]

tests/ode.rs

@@ -12,7 +12,7 @@ fn first_order_ode_solver() {
     let solver = ODESolver::new(f, x0, y0, step_size);
     let solution = solver.solve(x_end).unwrap();
 
-    assert!((solution - SOLUTION).abs() < 1e-2);
+    assert!((solution - SOLUTION).abs() < 1e-6);
     assert!(match solver.solve(-1.).unwrap_err() {
         SolverError::IncorrectInput { details: _ } => true,
         _ => false,
@@ -30,5 +30,43 @@ fn ode_system_solver() {
     let solver = ODESolver::new(f, x0, y0, step_size);
     let solution = solver.solve(x_end).unwrap();
 
-    assert!((solution[0] - SOLUTION).abs() < 1e-2);
+    assert!((solution[0] - SOLUTION).abs() < 1e-6);
+}
+
+#[test]
+fn rk4_order_of_accuracy() {
+    // RK4 is 4th order: halving the step size must cut the error by roughly 2^4 = 16.
+    // The off-by-one in solve dropped the final step, degrading this to 1st order (#15).
+    let f = |_t: f64, y: f64| -y;
+    let (x0, y0, x_end) = (0., 1., 1.);
+    let exact = (-1.0_f64).exp();
+
+    let err = |steps| {
+        (ODESolver::new(f, x0, y0, 1.)
+            .with_steps(x_end, steps)
+            .solve(x_end)
+            .unwrap()
+            - exact)
+            .abs()
+    };
+    let ratio = err(50) / err(100);
+
+    assert!(ratio > 8., "expected ~16 for 4th-order RK4, got {ratio}");
+}
+
+#[test]
+fn solve_returns_value_at_x_end() {
+    // solve must return y(x_end); the off-by-one returned y(x_end - h).
+    let f = |_t: f64, y: f64| -y;
+    let solution = ODESolver::new(f, 0., 1., 1e-3).solve(1.).unwrap();
+    assert!((solution - (-1.0_f64).exp()).abs() < 1e-8);
+}
+
+#[test]
+fn solve_reaches_x_end_for_indivisible_step() {
+    // A step size that does not divide the interval must still land on x_end,
+    // not stop a whole step short.
+    let f = |_t: f64, y: f64| -y;
+    let solution = ODESolver::new(f, 0., 1., 0.3).solve(1.).unwrap();
+    assert!((solution - (-1.0_f64).exp()).abs() < 1e-3);
 }