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?