Bridge SDE noise when a late tstop shortens dt

[ChrisRackauckas-Claude]

Note

Draft — please ignore until reviewed by @ChrisRackauckas.

Important

This does not fix #3175, and #3175 should stay open. The error that issue reports no longer reproduces on master — see the section at the bottom. This PR fixes a different defect found while reviewing #4041.

The bug

init draws the pending Brownian increment for the configured dt. If the step is then shortened before perform_step! runs — add_tstop! after init — the solver steps the new, shorter dt while the noise process still holds an increment drawn for the old one. accept_step! advances W.curt by that stale W.dt, so the noise grid slides off the solution grid and never returns:

f(u, p, t) = 1.0
g(u, p, t) = 0.1
prob = SDEProblem(f, g, [1.0], (0.0, 1.0))

i = init(prob, EM(); dt = 0.02, save_noise = true)
add_tstop!(i, 0.01)
solve!(i)

i.sol.t[1:4]   # [0.0, 0.01, 0.03, 0.05]
i.W.t[1:4]     # [0.0, 0.02, 0.04, 0.06]   <-- noise on the old grid
i.sol.t[end]   # 1.0
i.W.t[end]     # 1.0099999999999998

The first increment also carries variance 0.02 across a step of 0.01.

The fix

Shrinking a pending step from the same start time is exactly what a step rejection already describes, so loopheader! routes it through the existing reject_noise! hook. DiffEqNoiseProcess.reject_step! bridges the drawn increment down to the new dt and keeps the remainder on the RSWM stack; W.curt/W.curW are untouched.

function shrink_noise_to_integrator_dt!(integrator)
    W = _get_W(integrator)
    isnothing(W) && return nothing
    if abs(integrator.dt) < abs(W.dt)
        reject_noise!(W, integrator.dt, integrator.u, integrator.p)
        reject_noise!(_get_P(integrator), integrator.dt, integrator.u, integrator.p)
        integrator.sqdt = integrator.tdir * sqrt(abs(integrator.dt))
    end
    return nothing
end

Instrumenting the comparison shows it firing only for a late add_tstop!; plain fixed-step and adaptive solves, tstops passed up front, dtmax-limited solves, saveat, and reverse-time solves never reach it, so existing random streams are unchanged.

Distributional verification

du = dW, where Euler-Maruyama is exact and u(t) is the Brownian path. 20,000 seeded paths per configuration, tolerances at 4σ:

=== late add_tstop!(0.01), dt = 0.02, EM: per-interval increments ===
var   ΔW[0, 0.01]      0.01009   target 0.01000   ok
var   ΔW[0.01, 0.03]   0.02036   target 0.02000   ok
var   ΔW[0.03, 0.05]   0.02030   target 0.02000   ok
cor   ΔW1, ΔW2        -0.00081   target 0.00000   ok
cor   ΔW1, ΔW3        -0.00146   target 0.00000   ok
cor   ΔW2, ΔW3        -0.01148   target 0.00000   ok
kurtosis ΔW1           3.05895   target 3.00000   ok
var   u(1)             1.00661   target 1.00000   ok

The bridged increment carries the variance of the step it spans, the stashed remainder is reused without double-counting (uncorrelated successors), and u(T) == W(T) held on every one of the 20,000 paths for: no tstop, tstop at init, one late tstop, three late tstops, adaptive, and RSwM1/RSwM2/RSwM3.

Since du = dW produces no step rejections, the reject path was checked separately on multiplicative noise (du = 1.01u dt + 0.87u dW) where W(1) ~ N(0, 1) is grid-independent:

multiplicative LambaEM tol 1e-3            11491 rejections   var W(1) 1.00078  kurt 2.977  ok
multiplicative LambaEM tol 1e-4, 3 tstops  50733 rejections   var W(1) 0.98184  kurt 2.950  ok
multiplicative SOSRI tol 1e-4              48000 rejections   var W(1) 0.98104  kurt 2.930  ok

Known gap: the grow direction

fix_dt_at_bounds! can also raise dt (max(dt, dtmin)), and growth cannot be bridged, so this PR leaves that case exactly as master has it — which is wrong, by a factor of 10⁴ in the reproduction. Filed as #4097 rather than fixed here, to keep this PR to the direction it verifies.

Tests

Added to lib/StochasticDiffEq/test/tstops_tests.jl, covering both failure directions:

• the noise grid must track sol.t, and W.curt == integrator.t at the end — these fail on master ([0.0, 0.02, 1.01] vs [0.0, 0.01, 1.0]);
• the noise must not take a step the solver never took: W.iter == integrator.iter, and u(T) == W(T) for du = dW — these guard against "fixing" the grid by committing the stale increment instead of bridging it;
• all three RSWM variants, since only RSwM1 hands a stack chunk back through W.dt.

31/31 pass here, 27/31 on unmodified master. A 34-file local SDE sweep passes. Runic clean.

Why #3175 is not fixed by this

StochasticDiffEq.jl#413, which #3175 was migrated from, reports:

Something went wrong. Integrator stepped past tstops but the algorithm was dtchangeable. Please report this error.

That error does not reproduce on master. Running the MWE across 8 variants × EM/LambaEM — float and integer tspan, scalar and vector u0, rational dt, tstop on/off the step grid, tstop equal to dt, tstop at 1e-9 — all 16 land on the tstop and finish with ReturnCode.Success.

The only part of #413's MWE that still fails is u0 = [1], which throws InexactError: Int64(0.78…) — an integer state vector that cannot hold the value after a continuous noise increment. That is unrelated to tstops and reproduces identically on master.

#4041 targets the same area and asserts step! lands at t == 0.01; that assertion also passes on master and on #4041's own parent commit. Separately, its resync_noise! calls DiffEqNoiseProcess.accept_step!, which advances W.curt, increments W.iter, folds the stale increment into W.curW, and pushes a point into the saved path, despite documenting itself as resyncing "without advancing the noise clock". On that branch du = dW gives |u(T) - W(T)| = 0.128 — exactly the discarded increment — and W.iter == 52 against 51 solver steps.

🤖 Generated with Claude Code

https://claude.ai/code/session_01QSn6HPV8G8NWLpUeAEVA9b

[ChrisRackauckas-Claude]

Local regression sweep against this rebase: 34/34 SDE test files pass (tstops, callable tstops, saveat, events, noise types, scalar/nondiagonal/noncommutative noise, adaptive suites, Stratonovich and ODE convergence, IIF, mass matrix, static arrays, jump RNG, zeroed noise).

tstops_tests.jl: 31/31 on this branch, 27/31 on unmodified master — the four that fail on master are the new grid/clock assertions (W.t == sol.t, W.curt == integrator.t).

On CI, the red checks so far are all pre-existing and reproduce on #4093 / #4094 / master:

check	also fails on
Documentation	master (incl. base commit 86fb5301), #4093, #4094
PositiveIntegrators.jl/Downstream	#4041, #4093, #4094 — Unsatisfiable requirements ... LinearSolve
ProbNumDiffEq.jl/Downstream	#4041, #4093, #4094
downgrade-sublibraries (lib/StochasticDiffEq)	#4093, #4094 — precompile InitError: type Nothing has no field major
sublibrary-ci lib/DelayDiffEq [Interface]	#4093
sublibrary-ci lib/OrdinaryDiffEqLinear	#4093
sublibrary-ci lib/StochasticDiffEqHighOrder [QA]	#4093

The lib/StochasticDiffEq [Interface*] groups that contain tstops_tests.jl are still queued.