Fix spurious tiny final step for fixed-dt methods

When solving with a fixed-dt method (e.g. `solve(prob, Euler(); dt = 0.1)`),
the accumulated `t + dt + dt + ...` drifts past `tspan[end]` by one ulp,
producing a spurious trailing micro-step. PR #2869 removed the
`100eps(tstop)` snap hack that previously masked this.

Fix: for non-adaptive, dtchangeable algorithms with no user-supplied
tstops / d_discontinuities / callbacks, expand `tstops` to the range
`tspan[1]:dt:tspan[end]` whose TwicePrecision arithmetic gives exact
floating-point tstops, and inflate `dt` by 10 ulps so that
`modify_dt_for_tstops!` always takes the tstop branch and snaps `t`
exactly via `fixed_t_for_tstop_error!`.

The expansion is intentionally skipped for adaptive algorithms (even
when used with `adaptive=false`), `CompositeAlgorithm`, non-finite
tspans, and any solve with callbacks — preserving all existing
stepping semantics in those cases.

Adds a regression test covering forward, reverse, and non-evenly-
dividing `dt` on `tspan = (0.0, 1.0)`.

Fixes the PositiveIntegrators.jl CI failure noted in
NumericalMathematics/PositiveIntegrators.jl#192.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

Author: ChrisRackauckas

lib/OrdinaryDiffEqCore/src/solve.jl

@@ -401,6 +401,33 @@ function _ode_init(
         _tstops_cache = tstops
         tstops = ()
     end
+
+    # For fixed-timestep methods with a user-supplied dt and no user-supplied
+    # tstops/d_discontinuities/callbacks, expand tstops to the range
+    # tspan[1]:dt:tspan[end]. Julia's range semantics use TwicePrecision to
+    # avoid the floating-point drift of repeated `t + dt + dt + ...`, so the
+    # integrator hits each step exactly and lands on tspan[end] without a
+    # spurious tiny final step. Skipped whenever the user supplied tstops,
+    # d_discontinuities, or callbacks (which may mutate dt at runtime via
+    # `set_proposed_dt!`) to preserve existing "step at dtcache" semantics,
+    # and when tspan has non-finite endpoints.
+    if !isnothing(dt) && !iszero(dt) && !isadaptive(_alg) &&
+            isdtchangeable(_alg) && !(_alg isa CompositeAlgorithm) &&
+            sign(dt) == tdir &&
+            isempty(tstops) && isempty(d_discontinuities) &&
+            callback === nothing && all(isfinite, tspan)
+        tstops = ((tspan[1] + dt):dt:tspan[end]...,)
+        # Inflate dt so that dtcache is strictly larger than every range
+        # spacing.  Julia's TwicePrecision range arithmetic can produce
+        # spacings up to ~6 ulps above dt; inflating by 10 ulps ensures
+        # modify_dt_for_tstops! always takes the "hit tstop" branch,
+        # clipping dt to the exact range value and snapping t via
+        # fixed_t_for_tstop_error!.
+        if dt isa AbstractFloat
+            dt = tdir * nextfloat(abs(dt), 10)
+        end
+    end
+
     tstops_internal = initialize_tstops(tType, tstops, d_discontinuities, tspan)
     saveat_internal = initialize_saveat(tType, saveat, tspan)
     d_discontinuities_internal = initialize_d_discontinuities(

test/interface/ode_tstops_tests.jl

@@ -272,6 +272,32 @@ end
     @test sol.u[end] ≈ 0.0 atol = 1.0e-10
 end
 
+@testset "Fixed-dt: no spurious tiny final step from accumulated drift" begin
+    # Regression test: with fixed dt = 0.1 and tspan = (0.0, 1.0) the
+    # accumulated `t + dt + dt + ...` drifts and used to produce a spurious
+    # 12th step at 0.9999999999999999 followed by 1.0. The fix expands
+    # tstops to tspan[1]:dt:tspan[end], whose Julia range semantics land
+    # exactly on tspan[end].
+    f!(du, u, p, t) = (du .= u; nothing)
+    u0 = [1.0]
+    prob = ODEProblem(f!, u0, (0.0, 1.0))
+    sol = solve(prob, Euler(); dt = 0.1)
+    @test length(sol.t) == 11
+    @test sol.t[end] == 1.0
+    @test sol.t == collect(0.0:0.1:1.0)
+
+    # Reverse direction
+    prob_rev = ODEProblem(f!, u0, (1.0, 0.0))
+    sol_rev = solve(prob_rev, Euler(); dt = -0.1)
+    @test length(sol_rev.t) == 11
+    @test sol_rev.t[end] == 0.0
+    @test sol_rev.t == collect(1.0:-0.1:0.0)
+
+    # dt that does not evenly divide tspan still hits tspan[end] cleanly
+    sol_uneven = solve(prob, Euler(); dt = 0.3)
+    @test sol_uneven.t[end] == 1.0
+end
+
 @testset "d_discontinuities: tprev advanced past discontinuity" begin
     # Directly verify the shift-past invariant using the integrator interface:
     # after stepping past the discontinuity, integrator.tprev should be