Hello, I have recently started to play out with DifferentialEquations.jl; I’m gradually reading through the docs but there’s a lot of stuff in there which will keep me busy for a while, so I’m still at a beginner level.
I would like to understand what would be the preferred way of modeling the following problem:
I have a Langevin equation
\begin{cases}
\dot{x} = v \\
\dot{v} = -Av + B\xi_1(t)\xi_2(t)
\end{cases}
where \xi_1 is Gaussian white noise and \xi_2 is another stochastic process, switching between values 0 and 1 according to a Poisson process with distinct (but constant) forward and backward rates (i.e., \lambda_0 is the rate for the 0\rightarrow1 transition, while \lambda_1 is the rate for the 1\rightarrow0 transition).
Fundamentally, the system can switch randomly from a “diffusive” mode (when \xi_2=1) to a “non-diffusive” mode (when \xi_2=0).
If \xi_2 is not involved the representation is straightforward, I would use the following code (assume appropriate initial conditions u0, integration interval tspan and parameters p are given)
function f!(du, u, p, t)
A, B = p
du[1] = u[2]
du[2] = -A*u[2]
end # function
function g!(du, u, p, t)
A, B = p
du[1] = 0.0
du[2] = B
end # function
prob = SDEProblem(f!, g!, u0, tspan, p)
sol = solve(prob)
At this point, I had some trouble understanding how to “correctly” couple \xi_2 only with the stochastic term of this SDEProblem.
The easiest solution I conceived is to introduce an additional control parameter s (and two further parameters included in p that represent the forward and backward rates for \xi_2)which can be switched between values 0 and 1 via a ConstantRateJump (constant rate should be fine since the time dependence is embedded in s and there is no explicit dependence on u) then I would have
function f!(du, u, p, t)
A, B, s, λ0, λ1 = p
du[1] = u[2]
du[2] = -A*u[2]
end # function
function g!(du, u, p, t)
A, B, s, λ0, λ1 = p
du[1] = 0.0
du[2] = B*s
end # function
u0 = [0.0, 1.0]
tspan = (0.0, 10.0)
p = [1.0, 1.0, 1.0, 0.02, 0.1]
prob = SDEProblem(f!, g!, u0, tspan, p)
function rate(u,p,t)
A, B, s, λ0, λ1 = p
(1-s)*λ0 + s*λ1
end # function
affect!(integrator) = (integrator.p[3] = 1.0 - integrator.p[3])
jump = ConstantRateJump(rate, affect!)
jump_prob = JumpProblem(prob, Direct(), jump)
sol = solve(jump_prob)
This approach seems to work properly (both using ConstantRateJump and VariableRateJump, so I assume ConstantRate is good enough; does this mean the rate function is re-evaluated at every step?).
The question however is, is this the “best” approach for the problem?
Is there another way in which I could explicitly include the \xi_2(t) term (which in general could also be different from a Poisson process) in the equation and solve it like a SDEProblem rather than a JumpProblem?
Also, if this approach is adopted, then is it possible to obtain the timeseries for the parameter s throughout the integration?