QuTiP mesolve: demonstate how to use the e_ops argument

In [1]:

    

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

from qutip import *


    

In [3]:

    

wc = 1.0  * 2 * pi  # cavity frequency
wa = 1.0  * 2 * pi  # atom frequency
g  = 0.05 * 2 * pi  # coupling strength
kappa = 0.005       # cavity dissipation rate
gamma = 0.05        # atom dissipation rate
N = 15              # number of cavity fock states
tlist = linspace(0,25,100)


    

In [4]:

    

a  = tensor(destroy(N), qeye(2))
sm = tensor(qeye(N), destroy(2))


    

In [5]:

    

psi0 = tensor(basis(N,0), basis(2,1))


    

In [6]:

    

H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())


    

In [7]:

    

c_ops = [sqrt(kappa) * a, sqrt(gamma) * sm]


    

In [8]:

    

e_ops = [a.dag() * a, sm.dag() * sm]


    

In [9]:

    

result = mesolve(H, psi0, tlist, c_ops, e_ops)


    

In [10]:

    

len(result.expect)


    

    
Out[10]:





2

In [11]:

    

result.expect[0].shape


    

    
Out[11]:





(100,)

In [12]:

    

plot_expectation_values(result);


    

In [13]:

    

e_ops = []


    

In [14]:

    

result = mesolve(H, psi0, tlist, c_ops, e_ops)


    

In [15]:

    

len(result.states)


    

    
Out[15]:





100

In [16]:

    

result.states[0]


    

    
Out[16]:





$\text{Quantum object: dims = [[15, 2], [15, 2]], shape = [30, 30], type = oper, isHerm = True}\\[1em]\begin{pmatrix}0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 1.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\\vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\\end{pmatrix}$

In [17]:

    

# loop over states and calculate the atom entropy
entropy_atom_vec = zeros(len(result.states))

for idx, state in enumerate(result.states):
    entropy_atom_vec[idx] = entropy_linear(ptrace(state, 0))


    

In [18]:

    

fig, ax = plt.subplots(1,1)

ax.plot(tlist, entropy_atom_vec)

ax.set_xlabel('time', fontsize=16)
ax.set_ylabel(r'$E(\rho_{\rm atom})$', fontsize=16);


    

In [19]:

    

entropy_cavity_vec = zeros(len(tlist))


    

In [20]:

    

# calculate the entropy of the cavity in a callback function

def e_ops_func(t, rho):
    # rho is only the data, so do Qobj(rho) to obtain a Qobj representation
    # of rho
    entropy_cavity_vec[e_ops_func.idx] = entropy_linear(ptrace(Qobj(rho), 1))
    e_ops_func.idx += 1
        
e_ops_func.idx = 0


    

In [21]:

    

result = mesolve(H, psi0, tlist, c_ops, e_ops_func)


    

In [22]:

    

fig, ax = plt.subplots(1,1)

ax.plot(tlist, entropy_atom_vec, label='atom')
ax.plot(tlist, entropy_cavity_vec, label='cavity')

ax.legend()
ax.set_xlabel('time', fontsize=16)
ax.set_ylabel(r'$E$', fontsize=16);


    

In [23]:

    

from qutip.ipynbtools import version_table
version_table()


    

    
Out[23]:




Software
Version
Python
3.3.1 (default, Apr 17 2013, 22:30:32) 
[GCC 4.7.3]
OS
posix [linux]
SciPy
0.14.0.dev-2a4ba40
QuTiP
2.3.0.dev-2283478
IPython
2.0.0-dev
Cython
0.19
matplotlib
1.4.x
Numpy
1.7.1
Sat Sep 28 11:41:47 2013 JST