Solving Stochastic Differential Equations in Different Time Windows

Specific Domains Numerics
,
16

I traced the issue back to being an upstream issue with RandomNumbers.jl on v1.0. If you use a seed of zero you always get zeros out of randn:

using RandomNumbers
rng = Xorshifts.Xoroshiro128Plus(0)
randn(rng,Float64) # 0.0

That issue was upstreamed: zero seeds give zeros · Issue #48 · JuliaRandom/RandomNumbers.jl · GitHub . So if you don’t use a seed of zero and you set reseed=false, then the random process is not reseeded and you repeat the values:

using StochasticDiffEq, RandomNumbers, DiffEqNoiseProcess

# Define the problem
f(dx, x, u, t) = (dx .= x)
g(dx, x, u, t) = (dx .= x)
dt = 2^(-6)

# Solve for 0-2 seconds
rng = Xorshifts.Xoroshiro128Plus(1)
x0 = [1 / 2]
tspan = (0., 2.)
noise = WienerProcess(0., 0., rng=rng, reseed=false)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_full = solve(prob, EM(), dt=dt)

# Solve for 0-1 and 1-2 seconds
rng = Xorshifts.Xoroshiro128Plus(1)
x0 = [1 / 2]
tspan = (0., 1.)
noise = WienerProcess(0., 0., rng=rng, reseed=false)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part1 = solve(prob, EM(), dt=dt)

x0 = sol_part1.u[end]
tspan = (sol_part1.t[end], 2.)
noise = WienerProcess(sol_part1.W.t[end], sol_part1.W.u[end], rng=rng, reseed=false)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part2 = solve(prob, EM(), dt=dt)

using Plots
plt1 = plot(sol_full.t, vcat(sol_full.u...),
    linewidth=2, label="Full solution", linecolor=:red)
plot!(sol_part1.t, vcat(sol_part1.u...),
    linewidth=2, label="Solution part 1", linecolor=:green)
plot!(sol_part2.t, vcat(sol_part2.u...),
    linewidth=2, label="Solution part 2", linecolor=:blue)

plt2 = plot(sol_full.W.t, vcat(sol_full.W.u...),
    linewidth=2, label="Full solution Wiener", linecolor=:magenta)
plot!(sol_part1.W.t, vcat(sol_part1.W.u...),
    linewidth=2, label="Solution part 1 Wiener", linecolor=:cyan)
plot!(sol_part2.W.t, vcat(sol_part2.W.u...),
    linewidth=2, label="Solution part 2 Wiener", linecolor=:orange)

plot(plt1, plt2, layout=(2, 1))

res

Hope that helps. Sorry that took so long.

3 Likes