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 and used to produce a spurious
trailing step at `1 - eps` followed by `1.0`. PR #2869 removed the
`100eps(tstop)` snap hack that previously masked this in
`fixed_t_for_floatingpoint_error!`.

Two coordinated changes restore the expected 11-step result without
reintroducing the old hack:

1. `_ode_init` expands `tstops` to `tspan[1]:dt:tspan[end]` when the user
   specifies `dt` for a fixed-time method and did not supply their own
   `tstops`/`d_discontinuities`. Julia's range semantics use TwicePrecision
   internally, so each range element is the exact floating-point
   representative of `k*dt` and lands on `tspan[end]` cleanly. The
   expansion is skipped whenever the user supplied any tstops, which
   preserves the existing "continue at dtcache between user tstops"
   behavior that tests like `sol.t == [0, 1/3, 1/2, 5/6, 1]` depend on.

2. `modify_dt_for_tstops!` compares `dt` against `distance_to_tstop` with
   a small floating-point tolerance (`100 * eps(max(t, tstop))`). Without
   this tolerance, `dtcache = 0.1` reads as strictly less than a
   `distance_to_tstop` of `0.10000000000000009` and the integrator takes
   a full step that overshoots the tstop by one ulp, then takes a
   matching tiny corrective step. The tolerance guard is gated on
   `isfinite(tdir_tstop) && isfinite(integrator.t)` so that semi-infinite
   `tspan = (0.0, Inf)` integrations still work.

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 where
`solve(prob, Euler(); dt = 0.1)` began returning 12 steps instead of 11.

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

Author: ChrisRackauckas

lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl

@@ -178,11 +178,22 @@ function modify_dt_for_tstops!(integrator)
         tdir_t = integrator.tdir * integrator.t
         tdir_tstop = first_tstop(integrator)
         distance_to_tstop = abs(tdir_tstop - tdir_t)
+        # Floating-point tolerance for the dt-vs-distance comparison so that a
+        # `dt` whose nominal value matches `distance_to_tstop` to within a few
+        # ulps still snaps to the tstop. Without this a fixed-dt run that should
+        # land exactly on a tstop drifts past it and takes a spurious tiny step.
+        tstop_tol = if integrator.t isa AbstractFloat && isfinite(tdir_tstop) &&
+                isfinite(integrator.t)
+            100 * eps(float(max(abs(integrator.t), abs(tdir_tstop)) /
+                            oneunit(integrator.t))) * oneunit(integrator.t)
+        else
+            zero(distance_to_tstop)
+        end
 
         if integrator.opts.adaptive
             original_dt = abs(integrator.dt)
             integrator.dtpropose = integrator.tdir * original_dt
-            if original_dt < distance_to_tstop
+            if original_dt + tstop_tol < distance_to_tstop
                 _set_tstop_flag!(integrator, false)
             else
                 _set_tstop_flag!(
@@ -198,7 +209,7 @@ function modify_dt_for_tstops!(integrator)
         elseif integrator.dtchangeable && !integrator.force_stepfail
             # always try to step! with dtcache, but lower if a tstop
             # however, if force_stepfail then don't set to dtcache, and no tstop worry
-            if abs(integrator.dtcache) < distance_to_tstop
+            if abs(integrator.dtcache) + tstop_tol < distance_to_tstop
                 _set_tstop_flag!(integrator, false)
             else
                 _set_tstop_flag!(

lib/OrdinaryDiffEqCore/src/solve.jl

@@ -401,6 +401,19 @@ function _ode_init(
         _tstops_cache = tstops
         tstops = ()
     end
+
+    # For fixed-timestep methods with a user-supplied dt and no user-supplied
+    # tstops, 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 when the user
+    # supplied tstops to preserve existing "step at dtcache between tstops"
+    # semantics.
+    if !isnothing(dt) && !iszero(dt) && (!isadaptive(_alg) || !adaptive) &&
+            sign(dt) == tdir && isempty(tstops) && isempty(d_discontinuities)
+        tstops = ((tspan[1] + dt):dt:tspan[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