From a86dbfaf8264e689446a88258ffea7f2437dc368 Mon Sep 17 00:00:00 2001 From: gaoflow Date: Thu, 23 Jul 2026 09:29:55 +0200 Subject: [PATCH 1/2] Fix ODE solver stopping one step short and composite Simpson tail-drop ODESolver::solve looped `1..steps`, running one iteration too few, so it returned y(x_end - h) instead of y(x_end). The missing step is an O(h) error that dominates the result and reduces every method, including the default RK4, from its nominal order down to first order. solve now runs the full step count and uses a uniform step that lands exactly on x_end, so a step size that does not divide the interval no longer stops short. The composite Simpson's 1/3 and 3/8 rules in NewtonCotes iterate the panels with windows(3).step_by(2) / windows(4).step_by(3), which drop the trailing subinterval(s) whenever the subdivision count is not a multiple of 2 or 3. This silently truncated the integral, including for the default of 1000 subdivisions with Simpson's 3/8. The count is now rounded up to the next valid multiple for the selected rule. Adds order-of-accuracy and endpoint regressions for the solver and exactness regressions for both Simpson rules across even/odd and non-multiple-of-three subdivision counts, and tightens the existing tolerances that were loose enough to mask both bugs. --- CHANGELOG.md | 6 ++++ src/integrators/newton_cotes.rs | 25 +++++++++++++++-- src/solvers/ode_solver.rs | 49 ++++++++++++++++----------------- tests/integrators.rs | 40 ++++++++++++++++++++++++++- tests/ode.rs | 42 ++++++++++++++++++++++++++-- 5 files changed, 131 insertions(+), 31 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f5ece71..1c18f4f 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/src/integrators/newton_cotes.rs b/src/integrators/newton_cotes.rs index b8c1f9d..5cdd5b5 100644 --- a/src/integrators/newton_cotes.rs +++ b/src/integrators/newton_cotes.rs @@ -48,6 +48,9 @@ pub struct NewtonCotes { f: F, formula: fn(&Self, &[T], T) -> T, subdivisions: usize, + /// Number of subintervals spanned by one panel of the current composite rule: + /// 1 for the trapezium, 2 for Simpson's 1/3, 3 for Simpson's 3/8. + panel: usize, } impl NewtonCotes @@ -63,6 +66,7 @@ where f, subdivisions: DEFAULT_SUBDIVISIONS, formula: Self::simpsons_one_third, + panel: 2, } } @@ -83,6 +87,11 @@ where /// assert!((integral - 1.631).abs() <= 1e-3); /// ``` pub fn with_formula(&mut self, formula: Formula) -> &mut Self { + self.panel = match formula { + Formula::Trapezium => 1, + Formula::SimpsonsOneThird => 2, + Formula::SimpsonsThreeEighths => 3, + }; self.formula = match formula { Formula::Trapezium => Self::trapezium, Formula::SimpsonsOneThird => Self::simpsons_one_third, @@ -128,8 +137,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 +157,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); } diff --git a/src/solvers/ode_solver.rs b/src/solvers/ode_solver.rs index 735a467..32b2fdf 100644 --- a/src/solvers/ode_solver.rs +++ b/src/solvers/ode_solver.rs @@ -67,8 +67,7 @@ pub struct ODESolver { x0: T, y0: V, h: T, - half_h: T, - method: fn(&Self, T, V) -> V, + method: fn(&Self, T, V, T) -> V, } impl ODESolver @@ -99,7 +98,6 @@ where x0, y0, h, - half_h: h / T::from(2_f64).unwrap(), method: Self::rk4_step, } } @@ -144,18 +142,23 @@ where /// assert!((solution[0] - SOLUTION) <= 1e-3); /// ``` pub fn solve(&self, x_end: T) -> SolverResult { - 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 = interval / T::from(steps).unwrap(); + + 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,8 +177,6 @@ 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 @@ -183,7 +184,6 @@ where /// ``` 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 +202,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 +210,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 +253,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()) } } diff --git a/tests/integrators.rs b/tests/integrators.rs index 9561ac1..6c3992a 100644 --- a/tests/integrators.rs +++ b/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. ∫₀³ 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. ∫₀³ 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] diff --git a/tests/ode.rs b/tests/ode.rs index 0e6b62b..7bbca79 100644 --- a/tests/ode.rs +++ b/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. + 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); } From 64bf7260094deab35307ee97636dcf98872beb30 Mon Sep 17 00:00:00 2001 From: Vincent Gao Date: Sun, 26 Jul 2026 17:17:08 +0200 Subject: [PATCH 2/2] Address review: panel_width const fn, subdivision docs, typed cast error - extract the formula->panel-width match into a const fn and use it in new() and with_formula() - document that with_subdivisions rounds up to the formula's panel multiple - return SolverError::TypeConversionError instead of unwrapping the step cast - reference #15 from the regression tests --- src/integrators/newton_cotes.rs | 24 ++++++++++++++++-------- src/solvers/ode_solver.rs | 6 +++++- tests/integrators.rs | 4 ++-- tests/ode.rs | 2 +- 4 files changed, 24 insertions(+), 12 deletions(-) diff --git a/src/integrators/newton_cotes.rs b/src/integrators/newton_cotes.rs index 5cdd5b5..944bb45 100644 --- a/src/integrators/newton_cotes.rs +++ b/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,8 +57,7 @@ pub struct NewtonCotes { f: F, formula: fn(&Self, &[T], T) -> T, subdivisions: usize, - /// Number of subintervals spanned by one panel of the current composite rule: - /// 1 for the trapezium, 2 for Simpson's 1/3, 3 for Simpson's 3/8. + /// Number of subintervals spanned by one panel of the current composite rule. panel: usize, } @@ -66,7 +74,7 @@ where f, subdivisions: DEFAULT_SUBDIVISIONS, formula: Self::simpsons_one_third, - panel: 2, + panel: panel_width(&Formula::SimpsonsOneThird), } } @@ -87,11 +95,7 @@ where /// assert!((integral - 1.631).abs() <= 1e-3); /// ``` pub fn with_formula(&mut self, formula: Formula) -> &mut Self { - self.panel = match formula { - Formula::Trapezium => 1, - Formula::SimpsonsOneThird => 2, - Formula::SimpsonsThreeEighths => 3, - }; + self.panel = panel_width(&formula); self.formula = match formula { Formula::Trapezium => Self::trapezium, Formula::SimpsonsOneThird => Self::simpsons_one_third, @@ -102,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 diff --git a/src/solvers/ode_solver.rs b/src/solvers/ode_solver.rs index 32b2fdf..eec6e3e 100644 --- a/src/solvers/ode_solver.rs +++ b/src/solvers/ode_solver.rs @@ -152,7 +152,11 @@ where // Take a uniform step that lands exactly on `x_end`, even when the requested // step size does not divide the interval evenly. - let h = interval / T::from(steps).unwrap(); + 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; diff --git a/tests/integrators.rs b/tests/integrators.rs index 6c3992a..e6a32aa 100644 --- a/tests/integrators.rs +++ b/tests/integrators.rs @@ -56,7 +56,7 @@ fn simpsons_three_eighths() { #[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. ∫₀³ x³ dx = 20.25. + // 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) @@ -72,7 +72,7 @@ fn simpsons_one_third_exact_for_cubics() { 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. ∫₀³ x³ dx = 20.25. + // 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) diff --git a/tests/ode.rs b/tests/ode.rs index 7bbca79..c67cad3 100644 --- a/tests/ode.rs +++ b/tests/ode.rs @@ -36,7 +36,7 @@ fn ode_system_solver() { #[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. + // 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();