odesolve raises NotImplemented error for time varying H as a lambda function

[GoogleCodeExporter]

The attached code snipped will throw an NotImplemented error from scipy's 
integrator. It is the basic example from the wiki on time varying Hamiltonians 
[1], except that the Hamiltonian is now a lambda function.

The lambda function returns the correct H(t) as an Qobj, so I would have 
expected this to work according to odesolve's docstring. 

The problem seems to be that the integrator from scipy got passed and Qobj 
instead of just its data.

I'm using

QuTip 1.0.0
Numpy 1.6.1
Scipy 0.9.0

on Mac OS X 10.5.8

[1] 
http://code.google.com/p/qutip/wiki/GuideDynamics#Time-dependent_Hamiltonians_%2
8unitary_and_non-unitary%29

Original issue reported on code.google.com by Markus.B...@gmail.com on 2 Oct 2011 at 3:37

Attachments:

• odesolve_lambda.py

[GoogleCodeExporter]

This is correct, the ODE solver gets passed a Qobj instead of a sparse matrix.  
At present, the Hamiltonian terms must be passed via the H_args argument.  They 
are then stripped down to sparse matrices and used in the callback function.  
Passing Qobjs would slow the ODE solvers down quite a bit as there is a large 
overhead when calling operations on the Qobj class.

Original comment by nonhermitian on 3 Oct 2011 at 1:05

• Changed state: Accepted

[GoogleCodeExporter]

Yes, the hamiltonian callback function must return the sparse matrix data for 
the Qobj since. You can access this data with the "data" method in the Qobj 
class, so that your lambda function could be defined like this:

hamiltonian_t = lambda t, args: -2*pi * 0.5  * sigmaz().data -2*pi * 0.05 * 
sigmax().data * sin(2*pi*t)

However, an additional complication is that the ODE solver expects the 
Hamiltonian in superoperator form when collapse operators are present, so in 
that case you would first have to convert each operator to a superoperator (you 
can use the function liouvillian for this):

hamiltonian_t = lambda t, args: -2*pi * 0.5  * liouvillian(sigmaz(), []).data 
-2*pi * 0.05 * liouvillian(sigmaz(), []).data * sin(2*pi*t)

It would be much more efficient to pass the sigmax and sigmay operators as 
arguments in the H_args list to the odesolve function, because then the 
conversion to superoperator form and sparse matrix extration is one 
automatically (and only once per operator). So you could do something like:

hamiltonian_t = lambda t, args: -2*pi * 0.5  * args[0] -2*pi * 0.05 * args[1] * 
sin(2*pi*t)
H_args = [sigmaz(), sigmax()]
expt_sz_td = odesolve(hamiltonian_t, psi0, tlist, c_op_list, [sigmaz()], H_args)

Original comment by jrjohans...@gmail.com on 3 Oct 2011 at 1:15

[GoogleCodeExporter]

Thanks for the explanation. I will review docstring on odesolve / wiki guide to 
time dependent Hamiltonians and clarify what the odesolve expects from the 
function.

Original comment by Markus.B...@gmail.com on 3 Oct 2011 at 9:32

• Changed state: Started
• Added labels: Component-Docs, Type-Enhancement
• Removed labels: Type-Defect

[GoogleCodeExporter]

Clarified docstring, probably will also clarify issue in wiki.

diff --git a/qutip/qutip/odesolve.py b/qutip/qutip/odesolve.py
index 6213196..b0ee298 100644
--- a/qutip/qutip/odesolve.py
+++ b/qutip/qutip/odesolve.py
@@ -53,7 +53,18 @@ def odesolve(H, rho0, tlist, c_op_list, expt_op_list, H_args=
     @param expect_ops *list/array* of expectation operators
     @param H_args *list/array* of arguments for time-dependent Hamiltonians
     @param options *Odeoptions* instance of ODE solver options
-    
+
+    Notes on using callback function:
+
+    odesolve transforms all Qobj objects to sparse matrices before
+    handing the problem to the integrator function. In order for
+    your callback function to work correctly, pass all Qobj objects
+    that are used in constructing the Hamiltonian via H_args. odesolve
+    will check for Qobj in H_args and handle the conversion to sparse
+    matrices. All other Qobj objects that are not passed via H_args
+    will be passed on to the integrator to scipy who will raise an
+    NotImplemented exception.
+
     """

Original comment by Markus.B...@gmail.com on 5 Oct 2011 at 8:40

[GoogleCodeExporter]

Original comment by nonhermitian on 22 Mar 2012 at 8:04

• Changed state: Done