Skip to content

Commit 2ad6f5a

Browse files
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 old `100eps(tstop)` snap in `fixed_t_for_floatingpoint_error!` that masked this. Fix: add a floating-point tolerance (`100 * eps(max(|t|, |tstop|))`) to the `dt < distance_to_tstop` comparison in `modify_dt_for_tstops!`. When `dt ≈ distance` within rounding the integrator now takes the tstop branch, and `fixed_t_for_tstop_error!` snaps `t` to the exact tstop value — eliminating the extra step. 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>
1 parent ba0dfa2 commit 2ad6f5a

2 files changed

Lines changed: 36 additions & 2 deletions

File tree

lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl

Lines changed: 13 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -178,11 +178,22 @@ function modify_dt_for_tstops!(integrator)
178178
tdir_t = integrator.tdir * integrator.t
179179
tdir_tstop = first_tstop(integrator)
180180
distance_to_tstop = abs(tdir_tstop - tdir_t)
181+
# Floating-point tolerance so that a dt whose nominal value matches
182+
# distance_to_tstop to within rounding still triggers the tstop
183+
# branch. Without this, accumulated `t + dt + dt + …` can drift
184+
# just past the last tstop and produce a spurious micro-step.
185+
tstop_tol = if integrator.t isa AbstractFloat && isfinite(tdir_tstop) &&
186+
isfinite(integrator.t)
187+
100 * eps(float(max(abs(integrator.t), abs(tdir_tstop)) /
188+
oneunit(integrator.t))) * oneunit(integrator.t)
189+
else
190+
zero(distance_to_tstop)
191+
end
181192

182193
if integrator.opts.adaptive
183194
original_dt = abs(integrator.dt)
184195
integrator.dtpropose = integrator.tdir * original_dt
185-
if original_dt < distance_to_tstop
196+
if original_dt + tstop_tol < distance_to_tstop
186197
_set_tstop_flag!(integrator, false)
187198
else
188199
_set_tstop_flag!(
@@ -198,7 +209,7 @@ function modify_dt_for_tstops!(integrator)
198209
elseif integrator.dtchangeable && !integrator.force_stepfail
199210
# always try to step! with dtcache, but lower if a tstop
200211
# however, if force_stepfail then don't set to dtcache, and no tstop worry
201-
if abs(integrator.dtcache) < distance_to_tstop
212+
if abs(integrator.dtcache) + tstop_tol < distance_to_tstop
202213
_set_tstop_flag!(integrator, false)
203214
else
204215
_set_tstop_flag!(

test/interface/ode_tstops_tests.jl

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -272,6 +272,29 @@ end
272272
@test sol.u[end] 0.0 atol = 1.0e-10
273273
end
274274

275+
@testset "Fixed-dt: no spurious tiny final step from accumulated drift" begin
276+
# Regression test: with fixed dt = 0.1 and tspan = (0.0, 1.0) the
277+
# accumulated `t + dt + dt + ...` drifts and used to produce a spurious
278+
# 12th step at 0.9999999999999999 followed by 1.0. A tolerance in
279+
# modify_dt_for_tstops! ensures the final tstop is hit cleanly.
280+
f!(du, u, p, t) = (du .= u; nothing)
281+
u0 = [1.0]
282+
prob = ODEProblem(f!, u0, (0.0, 1.0))
283+
sol = solve(prob, Euler(); dt = 0.1)
284+
@test length(sol.t) == 11
285+
@test sol.t[end] == 1.0
286+
287+
# Reverse direction
288+
prob_rev = ODEProblem(f!, u0, (1.0, 0.0))
289+
sol_rev = solve(prob_rev, Euler(); dt = -0.1)
290+
@test length(sol_rev.t) == 11
291+
@test sol_rev.t[end] == 0.0
292+
293+
# dt that does not evenly divide tspan still hits tspan[end] cleanly
294+
sol_uneven = solve(prob, Euler(); dt = 0.3)
295+
@test sol_uneven.t[end] == 1.0
296+
end
297+
275298
@testset "d_discontinuities: tprev advanced past discontinuity" begin
276299
# Directly verify the shift-past invariant using the integrator interface:
277300
# after stepping past the discontinuity, integrator.tprev should be

0 commit comments

Comments
 (0)