SDDE Problems do not guarantee maximum step sizes #51
Description
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.
