### Author:  [@zekeriya.sari](https://discourse.julialang.org/u/zekeriya.sari)
#### Post date: [September 11, 2018, 11:58am UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/1 "2018-09-11T11:58:33Z")
I want to solve a stochastic differential equation in different time windows under Wiener process, e.g., instead of solving the system for a time interval of (0, 2) seconds, I want to solve the system for a time duration of (0, 1) first, and then for a time duration of (1, 2). However, although I feed the problem in the second time window with the properly modified Wiener process (I update the initial time and initial value of the wiener process in the second time window with the last time and value of the wiener process in the first time window), the solutions of the system solved in two time windows and the that of the system solved at once do not coinside.
Below are the script I wrote for this problem and the plot of the solutions I obtained. Any ideas what I am doing wrong?
```julia
using DifferentialEquations
using Random
using Plots; plotly()
# 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
Random.seed!(0)
x0 = [1 / 2]
tspan = (0., 2.)
noise = WienerProcess(0., 0.)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_full = solve(prob, EM(), dt=dt)
# Solve for 0-1 and 1-2 seconds
Random.seed!(0)
x0 = [1 / 2]
tspan = (0., 1.)
noise = WienerProcess(0., 0.)
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])
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part2 = solve(prob, EM(), dt=dt)
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))
gui()
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 11, 2018, 3:18pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/2 "2018-09-11T15:18:40Z")
> [@zekeriya.sari](#):
>
> Random.seed!(0)
You’re using global RNG seeding while the noise processes themselves have local RNGs.
>
The docs have been updated to reflect this. So what you’ll want to do is something like:
```julia
using RandomNumbers
rng = Xorshifts.Xoroshiro128Plus(0)
x0 = [1 / 2]
tspan = (0., 1.)
noise = WienerProcess(0., 0., rng = rng)
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)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part2 = solve(prob, EM(), dt=dt)
```
to have it continue with the same RNG.
---
### Author:  [@zekeriya.sari](https://discourse.julialang.org/u/zekeriya.sari)
#### Post date: [September 11, 2018, 3:41pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/3 "2018-09-11T15:41:21Z")
Thank you so much for your reply. But, I think that `SDEProblem` type does not accept `rng` as an argument and I think I should write
```julia
rng = Xorshifts.Xoroshiro128Plus(0)
noise = WienerProcess(0., 0., rng=rng)
```
since `NoiseProcess` has `rng` field. Athough I applied your suggestion, it did not change much, still the solutions do not coinside.
```julia
using DifferentialEquations
using RandomNumbers
using Plots; plotly()
# 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(0)
x0 = [1 / 2]
tspan = (0., 2.)
noise = WienerProcess(0., 0., rng=rng)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_full = solve(prob, EM(), dt=dt)
# Solve for 0-1 and 1-2 seconds
x0 = [1 / 2]
tspan = (0., 1.)
noise = WienerProcess(0., 0., rng=rng)
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)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part2 = solve(prob, EM(), dt=dt)
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))
gui()
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 11, 2018, 4:01pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/4 "2018-09-11T16:01:25Z")
> [@zekeriya.sari](#):
>
> But, I think that `SDEProblem` type does not accept `rng` as an argument and I think I should write
Yeah, that was a typo in my post. Edited to be fixed.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 11, 2018, 4:02pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/5 "2018-09-11T16:02:09Z")
You need to reset the RNG for the second run, like:
```julia
# Solve for 0-2 seconds
rng = Xorshifts.Xoroshiro128Plus(0)
x0 = [1 / 2]
tspan = (0., 2.)
noise = WienerProcess(0., 0., rng=rng)
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(0)
x0 = [1 / 2]
tspan = (0., 1.)
noise = WienerProcess(0., 0., rng=rng)
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)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part2 = solve(prob, EM(), dt=dt)
```
---
### Author:  [@zekeriya.sari](https://discourse.julialang.org/u/zekeriya.sari)
#### Post date: [September 11, 2018, 5:03pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/6 "2018-09-11T17:03:44Z")
Resetting the RNG in the second run did not solve the problem, either. Here is the case.
```julia
# Solve for 0-2 seconds
rng = Xorshifts.Xoroshiro128Plus(0)
x0 = [1 / 2]
tspan = (0., 2.)
noise = WienerProcess(0., 0., rng=rng)
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(0) # Reset the RNG!
x0 = [1 / 2]
tspan = (0., 1.)
noise = WienerProcess(0., 0., rng=rng)
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)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part2 = solve(prob, EM(), dt=dt)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 11, 2018, 6:12pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/7 "2018-09-11T18:12:04Z")
Oh you might need to set `reseed=false`.
```julia
# Solve for 0-2 seconds
rng = Xorshifts.Xoroshiro128Plus(0)
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(0) # Reset the RNG!
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)
```
That should be it. It automatically reseeds so that way giving the same noise process can give different solves which is normally wanted with a Monte Carlo, but here that’s not what you want.
---
### Author:  [@zekeriya.sari](https://discourse.julialang.org/u/zekeriya.sari)
#### Post date: [September 11, 2018, 6:29pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/8 "2018-09-11T18:29:40Z")
`reseed=false` made the solutions coincide, but this time I lost the WienerProcess. Here is the case,
```julia
# Solve for 0-2 seconds
rng = Xorshifts.Xoroshiro128Plus(0)
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(0)
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)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 11, 2018, 7:06pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/9 "2018-09-11T19:06:33Z")
Yeah you’re running into a bug here, but I’m not exactly sure what it is. I’ll need to sit down and really hammer through this. I’ll see if I find time later today.
---
### Author:  [@zekeriya.sari](https://discourse.julialang.org/u/zekeriya.sari)
#### Post date: [September 11, 2018, 7:19pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/10 "2018-09-11T19:19:57Z")
I am really grateful for your efforts. I am going through the documentation and source codes of the ecosystem to find a clue to cause of the problem. Thank you for all.
---
### Author:  [@mschauer](https://discourse.julialang.org/u/mschauer)
#### Post date: [September 11, 2018, 7:42pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/11 "2018-09-11T19:42:48Z")
If you have to handle the Wiener processes explicitly than you can look into using [GitHub - mschauer/Bridge.jl: A statistical toolbox for diffusion processes and stochastic differential equations. Named after the Brownian Bridge.](https://github.com/mschauer/Bridge.jl) :
```julia
using Bridge
using Random
using Plots
# Define the problem
f(t, x) = x
g(t, x) = x
dt = 2^(-6)
# Solve for 0-2 seconds
Random.seed!(0)
x0 = 1/2
t = 0.:dt:2.
W = sample(t, Wiener())
X = solve(Euler(), x0, W, (f, g))
Random.seed!(0)
x0 = 1/2
t1 = 0.:dt:1.
t2 = 1.:dt:2.
W1 = sample(t1, Wiener())
W2 = sample(t2, Wiener())
X1 = solve(EulerMaruyama(), x0, W1, (f, g))
X2 = solve(EulerMaruyama(), X1.yy[end]+0.1, W2, (f, g)) # shift by 0.1 to make path visible
plot(X)
plot!(X1)
plot!(X2)
```
---
### Author:  [@zekeriya.sari](https://discourse.julialang.org/u/zekeriya.sari)
#### Post date: [September 11, 2018, 8:37pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/12 "2018-09-11T20:37:51Z")
Bridge.jl did the job. This is the first time I tried Bridge.jl and it seems a very good tool to deal with stochastic differential equations.
The thing is that I need a syntax for drift and diffusion terms of the SDE as `f(dx, x, u, t), g(dx, x, u, t)` or `f(x, u, t) , g(x, u, t)`. Here `x` refers to the state, `dx` refers to the derivative of the state and `u` which is a function of `t` refers to the input and `t` is the current time of the state of a dynamical system modelled by stochastic differential equations. Currently, I do not know how to apply this syntax to the one required by Bridge.jl which, I think, requires drift and diffusion functions of the form
```julia
Bridge.b(t, x, P::OrnsteinUhlenbeck) = -P.β * x
Bridge.σ(t, x, P::OrnsteinUhlenbeck) = P.σ
Bridge.a(t, x, P::OrnsteinUhlenbeck) = P.σ^2
```
as its documentation says. Maybe I need more reading about Bridge.jl. Thank you for the advice.
---
### Author:  [@mschauer](https://discourse.julialang.org/u/mschauer)
#### Post date: [September 11, 2018, 9:41pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/13 "2018-09-11T21:41:57Z")
I am not sure what you need. Maybe the following is meaningful for you:
```julia
# Controlled example
struct ControlledOU{S} <: ContinuousTimeProcess{Float64}
t0::Float64
dt::Float64
u::Matrix{S}
end
ttoi(t, P) = 1 + round(Int, (t-P.t0)*P.dt)
function Bridge.b!(t, y, out, P::ControlledOU)
for i in 1:2
out[i] = -0.2 * y[i] + P.u[i, ttoi(t, P)]
end
out
end
Bridge.σ!(t, y, dw, out, P::ControlledOU) = out .= y.*dw
x0 = [1/2, 0.5]
dt = 2^(-6)
t = 0.:dt:10.
n = length(t)
u = rand(2,n)
P = ControlledOU(t[1], dt, u)
W = VSamplePath(t, zeros(2, n))
X = VSamplePath(t, zeros(2, n))
sample!(W, Wiener())
solve!(EulerMaruyama!(), X, x0, W, P)
plot(X)
```
---
### Author:  [@zekeriya.sari](https://discourse.julialang.org/u/zekeriya.sari)
#### Post date: [September 12, 2018, 8:28am UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/14 "2018-09-12T08:28:35Z")
@mschauer I read the docs and studied your code. I think the example you gave solves a differential equation of the form.
```julia
dx1 = (-0.2 *x1 + u1)dt + x1dW
dx2 = (-0.2 *x2 + u2)dt + x2dW
```
But, if I am not wrong the input function `u` is pre-generated as a vector in the with the line
```julia
u = rand(2, n)
```
and then the value of `u(t)` is then indexed using the function `ttoi`. This is not the case I wanted. But,I do not want the input function to be pre-generated but be a function that can be called at any time without being indexed by considering the initial time `t0` and time step `dt` of the Wiener process `P`. But, this is a very good example for a new starter of Bridge.jl. Thank you very much for the example.
By the way, I cannot find the docs for `VSamplePath` in the docs but it is in the source code. I wanted to let you know in case you missed it.
---
### Author:  [@mschauer](https://discourse.julialang.org/u/mschauer)
#### Post date: [September 13, 2018, 4:01pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/15 "2018-09-13T16:01:51Z")
Thank you for the feedback! I am not sure if I understand but you can of course play a bit with the design:
```julia
struct ControlledOU{S,T} <: ContinuousTimeProcess{Float64}
t0::Float64
dt::Float64
u::T
end
function Bridge.b!(t, y, out, P::ControlledOU)
u = P.u(t)
for i in 1:2
out[i] = -0.2 * y[i] + u[i]
end
out
end
P = ControlledOU(t[1], dt, sin)
```
…
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 30, 2018, 6:56am UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/16 "2018-09-30T06:56:20Z")
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`:
```julia
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](https://github.com/sunoru/RandomNumbers.jl/issues/48) . 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:
```julia
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))
```
Hope that helps. Sorry that took so long.
---
### Author:  [@zekeriya.sari](https://discourse.julialang.org/u/zekeriya.sari)
#### Post date: [September 30, 2018, 9:54am UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/17 "2018-09-30T09:54:08Z")
Thank you for your time and all your efforts @ChrisRackauckas . This solved the problem.
---
### Author:  [@ymardoukhi](https://discourse.julialang.org/u/ymardoukhi)
#### Post date: [November 24, 2023, 9:16pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/18 "2023-11-24T21:16:56Z")
@ChrisRackauckas I was trying to reproduce this plot by just rerunning your script. But it seems there is again a mismatch between the Wiener process in the second window compared to the full solution. This is what I get.
Notice that the Wiener processes in the first and second windows follow the same trend, which is odd. I tried different seeds, and the behaviour is still the same. Here is the plot when the seed is `2`.
I am using Julia 1.9.3; here is the version of the packages I use.
```julia
BenchmarkTools v1.3.2
DiffEqCallbacks v2.33.1
DiffEqNoiseProcess v5.19.0
DifferentialEquations v7.11.0
Distributions v0.25.103
ForwardDiff v0.10.36
Interpolations v0.14.7
RandomNumbers v1.5.3
StaticArrays v1.6.5
StochasticDiffEq v6.63.2
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 24, 2023, 10:14pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/19 "2023-11-24T22:14:28Z")
What is the script you’re running? The script has changed in the last 5 years, this was pre Julia v1.0. Check the docs on DiffEqNoiseProcess.
---
### Author:  [@ymardoukhi](https://discourse.julialang.org/u/ymardoukhi)
#### Post date: [November 26, 2023, 5:59pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/20 "2023-11-26T17:59:23Z")
The script is the same one you posted as the [solution](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807/16).