Many thanks for your swift answer!
Yes I also used the available PDMPs algorithms, with alterante results. Sometimes they work very well sometines they don’t. My problem is not having a code that works for a specific problem in itself.
As I understand from the documentation, my question comes down to the fact that the underlying logic of the algorithms used in DifferentialEquations.jl and in the other toolboxes appears to be different. Is this that I would like to understand better.
The integration strategy followed by QuTip and (more or less) by QuantumOprics.jl is the following:
-
Choose a random number r1 in [0,1] representing the probability that a quantum jump occurs.
-
Choose a random number r2 in [0,1] used to select which collapse operator was responsible for the jump.
-
Integrate the Schrödinger equation, using the effective Hamiltonian e.g.
G = i H + a C1 C1’ +b C2 C2’
until a time tau such that the norm of the wave function psi satisfies
sum(abs2, psi) = r1
at which time a jump occurs.
- The resultant jump projects the system at time tau into one of the renormalized states given by
psi_new = C1 \psi / sqrt(sum(abs2, C1 * psi) )
if
sum(abs2, C1 * psi) / (sum(abs2, C1 * psi) +sum(abs2, C2 * psi) ) >= r2
else
psi_new = C2 \psi / sqrt(sum(abs2, C2 * psi) )
5 Repeat until final time of the simulation
–
On a single realization I can implement this by introducing the condition 3
on the deterministic dissipative evolution. The condition triggers a continuous callback that implements 4.
Problem is how to pass r1 & r2 to the integrator. They must be independent random variables at every time step. Furthermore it must be possible to remake the simulation with different realizations of r1 & r2 for an arbitrary number of trajectories. A RODE algorithm with a custom noise seems perhaps the way to go but the problem is that the noise does not appear in the drift or the diffusion fields but must be only supplied to condition and callback.
Is it just me or does this mean that the exisisting tools cannot implement the QuTip algorithm?
Form this old discussion
[Code for simulating jump diffusion?]
it seems in fact that the DifferentialEquations algorithms have been developed having other applications in mind.
Many thanks!!!