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

[gaoflow]

Two related correctness bugs, both cases where an iteration loop stops short of covering the full integration domain.

ODESolver::solve returns y(x_end - h)

solve loops over 1..steps, so it runs one iteration too few and returns y(x_end - h) instead of y(x_end). The dropped step is an O(h) error that dominates the result and pulls every method — including the default RK4 — down to first-order accuracy.

let f = |_x: f64, y: f64| -y;                    // y' = -y, exact y(1) = e^-1
let sol = ODESolver::new(f, 0., 1., 1e-3).solve(1.).unwrap();
// sol = 0.36824750…  == e^-0.999 == y(1 - h)

Halving h then only halves the error (first order) where RK4 should cut it by ~16x. 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 a whole step short either.

Composite Simpson rules drop the trailing subinterval(s)

simpsons_one_third (windows(3).step_by(2)) and simpsons_three_eighths (windows(4).step_by(3)) skip the last subinterval(s) when the subdivision count is not a multiple of 2 or 3. Since the default is 1000, Simpson's 3/8 is wrong out of the box:

let f = |x: f64| x * x * x;                       // ∫₀³ x³ = 20.25, exact for Simpson
NewtonCotes::new(f)
    .with_formula(Formula::SimpsonsThreeEighths)
    .integrate(0., 3.).unwrap();                  // 20.169…, off by ~0.08

The count is now rounded up to the next valid multiple for the selected rule. (AdaptiveNewtonCotes evaluates each rule over a single panel, so it is unaffected.)

Notes

• Behavior: solve integrates to x_end using round((x_end - x0) / h) uniform steps; for a step size that already divides the interval — all documented usage and anything from with_steps — this is unchanged. Simpson may use up to two extra subintervals to reach a valid count.
• Tests: adds order-of-accuracy + endpoint regressions for the solver and exactness regressions for both Simpson rules across even/odd and non-multiple-of-three counts, and tightens the existing ODE/Simpson tolerances that were loose enough to pass the wrong values. cargo test, cargo fmt -- --check, and cargo clippy all pass.

[AzeezDa]

Hi @gaoflow,

Thank you very much for catching these errors; I completely missed them!
The changes look good, but I only have a few comments on the documentation and readability. Feel free to fix them yourself, but if you want I can fix them as well.

Otherwise, when my comments are amended, the pull request is approved to be merged 😃!

Thanks again!

PS: Feel free to mention/link to this pull request in the tests you added as you seem to hint at there being an error that was fixed.

[AzeezDa]

Instead of using the raw panel width immediately in e.g. line 69, I suggest moving the match statement of line 90-94 here as a const function as below, then using that function in the rest of the file (see comments for at line 69 and 90-94).

(Feel free to change the function name and document it.)

Suggested change

	const fn panel_width(formula: &Formula) -> usize {
	match formula {
	Formula::Trapezium => 1,
	Formula::SimpsonsOneThird => 2,
	Formula::SimpsonsThreeEighths => 3,
	}
	}

Another alternative is to keep the match statement as it is in your original pull request, but instead of writing 1, 2, 3 (in lines 69, 90-94), constants (e.g. const TRAPEZIUM_PANEL_WIDTH: usize = 1;) are instead added to the file and used later on. I have no preference with either of the alternatives.

[AzeezDa]

If the previous comment is implemented, then this line would become:

Suggested change

	panel: 2,
	panel: panel_width(&Formula::SimpsonsOneThird)

[AzeezDa]

If the first comment is implemented this line would simply be:

Suggested change

	self.panel = match formula {
	Formula::Trapezium => 1,
	Formula::SimpsonsOneThird => 2,
	Formula::SimpsonsThreeEighths => 3,
	};
	self.panel = panel_width(&formula);

[AzeezDa]

Since the subdivisions internally will be adjusted according to the formula, I suggest writing this explicitly here to increase transparency for the users.

Suggested change

	/// 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.

[AzeezDa]

From what I can see, this unwrap should be not be problematic because it would've been caught in the previous unwrap_or and the if-statement under it. Nonetheless, if there is a case where it would fail, I suggest adding the below. (What do you think?)

Suggested change

	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);
	};

[gaoflow]

Thanks — all five are in (64bf726).

• panel_width(&Formula) as a const fn next to the enum, used by both new() and with_formula(); dropped the now-redundant comment on the panel field.
• with_subdivisions docs use your wording about rounding up to the formula's panel multiple.
• solve() returns SolverError::TypeConversionError instead of unwrapping the T::from(steps) cast. You're right that it's unreachable today given the unwrap_or(0) + zero check above, but a T with a narrower range than usize would make it reachable, so the explicit branch is worth having.
• Added (#15) to the three regression-test comments that reference the fixed behaviour.

cargo clippy, cargo fmt -- --check and cargo test (135 pass) are clean.

[AzeezDa]

Thank you!