Complex Jump-Diffusion process (from heterodyne detection trajectory):

[WouterJRV]

Continuation of a discussion with Chris on stackexchange:

I am trying to find a good solver for the equation d v=(-1idt(H-0.5i(K+gamma+R)nop+0.5i(gamma+R)(nexpunit))+0.5Kdtconj(aexp)a-0.25Kdtabs(aexp).^2unit+sqrt(K)(a-0.5aexpunit)dZ-0.5sqrt(K)conj(aexp).conj(dZ)unit+(a/aexp-unit)dN+(adag/conj(aexp)-unit)dM)v;

where v is a complex-valued column-vector, a an annihilation operator(matrix representation) adag the hermitian conjugate of a (creation operator), nop=adaga the particle number operator, aexp=v^daggerav; nexp=v^daggernop*v.

As for the noise: dZ is a scalar such that dZ dZ=dt, dN, dM poissonian such that=nexpgamma and =(nexp+1)*R.

Regarding numerical values, gamma and R are of order 1 and K of order 0.001, we are interested in the regime with density 0<nexp<120, which should be contained if v has length~150.

The equation above should already be normconserving by itself, though there is also an alternative formulation that must be renormalized after each step (it won't hurt to do the same here, though) but that one is diagonal in the noise.

I've been trying plain Euler-Mayurama, but it forced me to take a timestep dt of 1e-5 or smaller in order to avoid numerical problems (and even then it was small enough only most of the time) when evolving it for hundreds of timeunits.

From here we found that for higher order methods the jump-parts can be separated from the drift-diffusion parts, at least for one jump (with both dN and dM present, perhaps their order should also be randomised?)

Regarding the drift-diffusion part, Chris found that the diffusion is commutative for at least the case of one equation (does that include matrix equations like here?). But even if it isn't exactly commutative, I wouldn't think that approximation would cause too much damage given the smallness of K.

Is JuliaDiffEq's RKMilCommutative method similar to Rössler's SRIC methods for this? or are there perhaps further suggestions?

Thanks for your advice.

[WouterJRV]

I can confirm that already separating the jumps from the drift-diffusion, caused large improvement to make a timestep 1e-4 work (haven't tested larger), even for EM-method. The parameters are as such that the probability of two jumps coinciding at any given time is low, but the probability a coincidence happens at some time is quite high.

[ChrisRackauckas]

A few things to unpack here. @onoderat should get in on this.

The equation above should already be normconserving by itself, though there is also an alternative formulation that must be renormalized after each step (it won't hurt to do the same here, though) but that one is diagonal in the noise.

I think this is what @onoderat was discussing. The noise looks diagonal, but since it's complex you're still having mixing between the real and imaginary parts, so it's not really diagonal. When you write it out as a 2n real system, you can see that for n=1 that holomorphic => it's commutative noise. I think that it holds for all n, and quick checking of cases makes it seem like that's the case, and @onoderat thinks he might be able to prove it.

So this puts it into the commutative noise category even when it looks diagonal. RKMilCommute is currently the best algorithm we have for that, and it's a lot like SRIC1 (it did come before, Rossler's paper makes an error here when he says there weren't any) and it's just from Kloeden and Platen (it's the same RKMil but with the commutative double integral handling). Going forward, for this to be efficient we want more commutative noise integrators since general noise is much much harder to do at order higher than 0.5 due to the Wiktorsson approximations. I don't think there's a silver bullet method right now. The low dt means we need more stability, but the higher stability and adaptive methods are currently for diagonal noise or are the S-ROCK methods. But S-ROCK (and Runge-Kutta Chebyshev in general) is notoriously unstable when there are complex values or oscillations and so this won't work for this problem. I am planning on working on a new method specifically for this now that I know there's an application. So yeah, that's a long way of saying that RKMilCommute (or if you know your derivative PCEuler is probably the best option right now but it's not good and I'm not happy about that!

From here we found that for higher order methods the jump-parts can be separated from the drift-diffusion parts, at least for one jump (with both dN and dM present, perhaps their order should also be randomised?)

When there are multiple jumps it's an SSA time stepping. The Gillespie method it's quite clear. The next jump is exponential with a rate of the sum of the rates, and which jump is a probability based on the rates. That's the gist of it.

The parameters are as such that the probability of two jumps coinciding at any given time is low, but the probability a coincidence happens at some time is quite high.

Technically the probability of two jumps happening at the same time is zero? They can be close but shouldn't ever be the same?

caused large improvement to make a timestep 1e-4 work (haven't tested larger), even for EM-method.

Is that 1e-4 limited by the jump timings or the SDE stability?

[WouterJRV]

Thanks, this Commutative Milstein case will be the first new algorithm if it turns out SDE stability is the limiting factor for the timestep. I've also found in Milstein-Tretyakov(Stochastic Numerics for mathematical physics) a chapter on 'small noise' problems, perhaps they could , alternatively, also be useful in this case.

Regarding the jumps, is the following what you have in mind for the Gillespie method at each timestep?

*First evolve the drift-diffusion part towards the 'almost sure left-hand limit t+dt-' of time t+dt *Pull a single Poisson number dA to determine the amount of jumps *if dA>0 for j=1:dA, -determine kind of jump j trough linear search or similar -let the jump j act v_old->v_new end end

I'd expect that would work indeed in the case above. Only in cases where dA> 1 would the for-loop make this much slower? Well, maybe not so much for Julia. Otherwise, the loop could perhaps be replaced by something like a random permutation of the jumps, from which their collective effect is constructed, acting at once on v_{t+dt-}.

[ChrisRackauckas]

Milstein methods give a lot less of a stability increase than you'd hope, so I wouldn't expect it to be that much better. But yes, I forgot about that book and the small noise methods can get around a lot of the Wiktorsson integral approximations. Those would be highly worthwhile to have.

As for Gillespie, I did a minimal SSA implementation over here:

https://github.com/JuliaDiffEq/DiffEqJump.jl/blob/master/src/simple_regular_solve.jl#L74-L88

That should explain it quite simply. What you do is calculate all of the rates, and then the time to the next jump is ttnj ~ rand(Exponential(lambda)) where lambda is that sum. When you hit t + ttnj, the jump occurs with the probability that is proportional to each rate, so do calculate that you just normalize the commutative sum of the rates and take a uniform random number find out in what interval it lies. That tells you when jumps occur and which jump it is. In every [t,t+ttnj] interval you solve the SDE (the one I linked to is a pure-Jump integrator, but it's clear how to stick an SDE in there). This is the form that mixes the time-adapted jump integrators of Platen with the Gillespie SSA. Note that this relies on the assumption that the jump rates are constant in the interval [t,t+ttnj]. If they are not constant then you need to solve a simultaneous ODE with event handling.

From this formulation, it should be pretty clear too that as the rate of jumps increases, your integrator dt becomes constrained by the choices of ttnj, and so this slows to a crawl if the jumps are frequent, but the jumps are exact so that's the tradeoff.

[ChrisRackauckas]

This discussion doesn't have anything actionable left so I'm closing it. More advanced jumping is going on in DiffEqJump.jl and more SDE integrators is StochasticDiffEq.jl.