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)
display(plt)
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