SDDE Problems do not guarantee maximum step sizes

[laserschwelle]

I tried to introduce noise to a laser system and this seems to break some behaviour regarding tstops, saveat and dt.

Given a SDDE dx(t) = f(x(t), x(t-τ))dt + g(x(t))dW and a corresponding DDE dx(t)/dt = f(x(t), x(t-τ)) I would assume that for g = 0 the solutions should be the same.

I also understand the documentation such that dt sets an upper bound to the step size irrespective of additional points set in tstops or the points set in saveat.

The following code reproduces the suspected bug:

using DifferentialEquations
using StochasticDelayDiffEq
using Plots

function f(du, u, h, p, t)    
    K, τ = p

    Eᵣ, Eᵢ, N, ρ = u
    S = Eᵣ^2 + Eᵢ^2

    ρᵗʰ = 14/23
    N₀ = 0.1 - ( ρᵗʰ / 400 )
    tmp = 1000*(ρᵗʰ + 0.2*(N-N₀)-ρ)

    du[1] = (2*ρ-1)*(230*Eᵣ + 115*Eᵢ) - 50*Eᵣ + K*h(p,t-τ)[1]*50    
    du[2] = (2*ρ-1)*(230*Eᵢ - 115*Eᵣ) - 50*Eᵢ + K*h(p,t-τ)[2]*50
    du[3] = 1.2 - N/0.17 - tmp
    du[4] = tmp - ρ/2 - 230*(2*ρ-1)*S
end

function g(du, u, h, p, t)
    du .= 0
end

h(p, t) = 1e-3 .* ones(4)

K = 0.05
τ = 1
p = [K, τ]
tspan = (0, 10)
dt = 1e-3

save_times = tspan[1]:0.2:tspan[2]
stop_times = tspan[1]:dt:tspan[2]

dde = DDEProblem(  f,  h(0, 0),h,tspan,p;constant_lags=[τ])
sdde = SDDEProblem(f,g,h(0, 0),h,tspan,p;constant_lags=[τ])

expected   = solve( dde,SROCK2(),dt=dt,saveat=save_times)
bug        = solve(sdde,SROCK2(),dt=dt,saveat=save_times)
no_fix     = solve(sdde,SROCK2(),dt=dt,saveat=save_times, tstops=stop_times)
workaround = solve(sdde,SROCK2(),dt=dt)

# helper function to plot laser intensities over time
function get_ts(sol)
    t = sol.t
    S = [u[1]^2+u[2]^2 for u in sol.u]
    return (t, S)
end

plot()
plot!(get_ts(expected),   label="expected")
plot!(get_ts(bug),        label="SDDE bug")
plot!(get_ts(no_fix),     label="SDDE no fix")
plot!(get_ts(workaround), label="SDDE workaround")
plot!(legend=:bottomright, xlabel="t", ylabel="S")

This produces the following plot. The orange and green curves lie on top of each other. The slight offset between the expected points and the workaround is also unexpected.

[devmotion]

I also understand the documentation such that dt sets an upper bound to the step size irrespective of additional points set in tstops or the points set in saveat.

No, you have to set dtmax to set an upper bound on the step size.

[devmotion]

expected = solve( dde,SROCK2(),dt=dt,saveat=save_times)

DDEProblems require a MethodOfSteps algorithm. I'm not sure why this does not throw an error.

[laserschwelle]

Using expected = solve(dde,MethodOfSteps(Tsit5()),dt=dt,saveat=save_times) fixes the offset between the expected result and the workaround, but the main problem remains.

[ChrisRackauckas]

DDEProblems require a MethodOfSteps algorithm. I'm not sure why this does not throw an error.

That's not so with SDDEProblems: it just dispatches directly on the problem + the original algorithm. But yeah, this is odd.

[ChrisRackauckas]

Forgot to mention that safety errors were added about a year back for this: https://github.com/SciML/StochasticDelayDiffEq.jl/pull/67