How to represent an Ornstein-like stochastic ODE in Julia?


1. The equation I want to solve

I have some vector stochastic ODE that has an almost Langevin form, but with an additional time dependence:

\frac{dy}{dt} = -\gamma (y - a - bt) + b + \xi(t)

where a,b,\gamma are constant parameters and \xi(t) is a white noise process characterised by variance \sigma (i.e. \langle \xi(t),\xi(t') \rangle = \sigma^2 \delta(t-t') )

This can also be written as
dy = \gamma (\mu(t) - y) dt + \xi dt = \gamma (\mu - y) dt + dW_t

for \mu(t) = a + b t + b and W_t a Wiener process

2. My attempt to implement in Julia

I want to integrate this equation in Julia.

A MWE might be something like

using DifferentialEquations, Noise, Plots 

#Vector of parameters
γ = [1e-13, 1e-14, 1e-15]
a = [300,400,500]
b = [-1e-15,2e-15,3e-15]

#Initial conditions on y 
y0 = [300,400,500]

#Define the ODEs 

f(du,u,p,t) = (du .= -γ.*u .+ γ.*(a .+ b*t) .+ b)
g(du,u,p,t) = (du .= 1.0)

#Specify the noise process to use 
noise = WienerProcess(0., 0.) # WienerProcess(t0,W0) where t0 is the initial value of time and W0 the initial value of the process

#Time over which to integrate 
t = range(0.0, 3e8, length=500)
tspan = (first(t),last(t))

#Setup SDE problem and solve
prob = SDEProblem(f,g,y0,tspan,tstops=t,noise=noise)
solution = solve(prob,EM())

#Plot it for interest
i = 1
solution_i = solution[i,:]
plt = plot(t,solution_i)

3. My question

How do I include \sigma as a parameter which can be set?

From the docs it suggests that the variance of the Wiener process is given by dt. I understand I could do something like noise.dt = σ or set the stepsize in the integrator accordingly (dropping tstops) but this then becomes a very long integration depending on choice of \sigma, and I only care about the value of the solution at tstops.

Another option is to include it in g(du,u,p,t), but then I am unsure how this “interacts” with the WienerProcess.

Thanks in advance for any help

When X ~ N(0,dt), sigma*X ~ N(0,sigma^2 * dt) (using N(mean, variance)). So yes, g(du,u,p,t) = (du .= p.σ) or however you want to write it would set the variance.

1 Like