mesolve issues

352 views
Skip to first unread message

Muir Kumph

unread,
Oct 25, 2017, 4:30:45 PM10/25/17
to QuTiP: Quantum Toolbox in Python
I am having a problem with mesolve (and sesolve) that might be caused by an underlying ODE integration bug.   I have tried to reproduce the issue here with the simplest code.  I took the time_dependent_interpolation notebook and it looks like I could reproduce the problem by making a few small tweaks.  It seems that if there is a Gaussian envelope on the time dependent part of a Hamiltonian, then at some critical value of the width of the Gaussian pulse, the simulation stops working.  If I use a tau of 1.02 - the attached notebook still integrates and gives a physical result.  If I use a value of tau=1.01, then the simulation goes south.  Note the right hand axes has a scale of 1e-102 when run with tau=1.01.  Its hard to believe that the pulse width should be so critical.  Is this a bug?  Is there a way to work around it?
time_dependent_interpolation-issues.ipynb

Paul Nation

unread,
Oct 25, 2017, 11:57:48 PM10/25/17
to QuTiP Group
When the pulse width is small, the solver has trouble dealing with the two different timescales in the problem.  Your system evolves much slower than the rate of change of the pulse, so it could be the case that the solver jumps over pulse when doing adaptive step sizing.   In your case, making the max step size small corrects the issue, even for the tau=0.1 case.  For example try:

options=Options(max_step=1e-4,nsteps=10000)

in the call to mesolve.

-P

On Oct 25, 2017, at 14:30, Muir Kumph <mku...@gmail.com> wrote:

I am having a problem with mesolve (and sesolve) that might be caused by an underlying ODE integration bug.   I have tried to reproduce the issue here with the simplest code.  I took the time_dependent_interpolation notebook and it looks like I could reproduce the problem by making a few small tweaks.  It seems that if there is a Gaussian envelope on the time dependent part of a Hamiltonian, then at some critical value of the width of the Gaussian pulse, the simulation stops working.  If I use a tau of 1.02 - the attached notebook still integrates and gives a physical result.  If I use a value of tau=1.01, then the simulation goes south.  Note the right hand axes has a scale of 1e-102 when run with tau=1.01.  Its hard to believe that the pulse width should be so critical.  Is this a bug?  Is there a way to work around it?

--
You received this message because you are subscribed to the Google Groups "QuTiP: Quantum Toolbox in Python" group.
To unsubscribe from this group and stop receiving emails from it, send an email to qutip+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
<time_dependent_interpolation-issues.ipynb>

Muir Kumph

unread,
Oct 26, 2017, 11:04:59 AM10/26/17
to QuTiP: Quantum Toolbox in Python
Thanks!  It seems strange that the solver could jump over a pulse that's so obvious to the eye.  A related problem I was working on had similar issues that were resolvable by doing the same trick.  It also seems that the solvers are very sensitive to the units being used.  If I try to work in nanoseconds (time units of 1e-9) with GHz frequencies there are a lot of un-real oscillations and it seems to miss the pulses.  If I uniformly scale everything up so that I am close to working with whole numbers, the simulation works well.  I imagine these are underlying problems with the ODE solver (ZVODE) in this case, yes?  It would be nice if we could use another solver to see if these problems go away.   I tried modifying sesolve.py to use the complex_ode class in scipy, but got an error in the integration step.  I guess ZVODE and complex_ode use slightly different dimensions on the result vectors.

Paul Nation

unread,
Oct 26, 2017, 12:32:09 PM10/26/17
to QuTiP Group
Well, if there is little change over a wide enough region, then it will adapt to large step sizes.  This is why the max_step option exists.

You should not be using any units in a simulation.  For example if you have a freq of 10* 2pi GHz, then just input 10 to the Hamiltonian, and the timescale of the evolution will naturally be nanosec.

None of these issues are problems with zvode, just how it is being used.  The solver itself has been in use for nearly thirty years, and is quite robust.



To unsubscribe from this group and stop receiving emails from it, send an email to qutip+unsubscribe@googlegroups.com.

Muir Kumph

unread,
Oct 26, 2017, 1:14:16 PM10/26/17
to qu...@googlegroups.com
You say to not use any units, but there are always some units, even if just implied.  In all the simulations I wrote before coming to Qutip, my data was in double floats (64 bit), so I could pretty much use whatever set of units made sense for the problem. SI is a popular choice if you want to compare to experiment.  Qutip seems to go south when I scale the whole Hamiltonian to work in frequencies of about 1e9.  The time scale goes to 1e-9, which means the timestep might go to 1e-13 or even smaller.  For double floats this should not matter.  They should work to about 1e-308.  Does ZVODE or Qutip use some smaller floats or fixed point arithmetic that I should be aware of?  Or maybe there is some hard coded rounding?  I am just trying to figure out what the limitations are so I can keep my simulations away from them.  One thing that is nice is this seems remarkably consistent across platforms (ie. Mac and Windows do the same thing).



--
You received this message because you are subscribed to a topic in the Google Groups "QuTiP: Quantum Toolbox in Python" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/qutip/slz4MCnkolw/unsubscribe.
To unsubscribe from this group and all its topics, send an email to qutip+unsubscribe@googlegroups.com.

Paul Nation

unread,
Oct 26, 2017, 3:31:48 PM10/26/17
to QuTiP Group
QuTiP uses double complex data, so two float64 numbers.  It maybe that you need to modify the atol and/or rtol values of the solver for you given parameters, but I think that that is besides the point.

Everything has units, but the computer does care, or know, about units.  Things like 1e9 are not needed, or particularly useful to have in a sim.  For example, an oscillator with angular frequency of 10ghz, or 10mhz, can be entered in the exact same way.

H = 10*a.dag()*a

The only difference comes from the interpretation of the units of time.  If, for example, you have

tlist = np.linspace(0,1,100)

Then the evolution corresponds to timescales of 1ns, or 1microsec, respectively, even though the simulation is the same.

One should also be cautious of the difference between the size of number that a float can hold, and the precision of a floating point number.
Reply all
Reply to author
Forward
0 new messages