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..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,6 +57,8 @@ 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. + panel: usize, } impl NewtonCotes @@ -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); } diff --git a/src/solvers/ode_solver.rs b/src/solvers/ode_solver.rs index 735a467..eec6e3e 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,27 @@ 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 = 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,8 +181,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 +188,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 +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()) } } diff --git a/tests/integrators.rs b/tests/integrators.rs index 9561ac1..e6a32aa 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 (#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] diff --git a/tests/ode.rs b/tests/ode.rs index 0e6b62b..c67cad3 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 (#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); }