# Demonstrate how to pass liouvillians to the mesolve¶

This notebook demonstrates how the mesolve solver can be used in different ways to solve the same problem, by passing a Hamiltonian and list of collapse operators, or by manually creating the Liouvillian and pass it to the solver in place of the Hamiltonian or collapse operators. These methods can be used to solve master equations that are not on standard Linblad form.

In [1]:
%pylab inline

In [2]:
from qutip import *


### Visualization of the mesolve results¶

In [3]:
def visualize(result):
fig, ax = subplots(1, 1, figsize=(8,5))

for e in result.expect:
ax.plot(result.times, e)

ax.set_xlabel('Time')
ax.set_ylabel('Occupation probability')
ax.set_ylim(0, 1)


### Parameters¶

In [4]:
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)


### Operators and initial states¶

In [5]:
psi0 = tensor(basis(N,0), basis(2,1))
a  = tensor(destroy(N), qeye(2))
sm = tensor(qeye(N), destroy(2))


### Standard method¶

In [6]:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())

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

result = mesolve(H, psi0, tlist, c_ops, [a.dag() * a, sm.dag() * sm])

visualize(result)


### Passing a liouvillian instead of an Hamiltonian¶

The Hamiltonian and the collapse operators are used to contruct a Liouvillian, which is passed to mesolve in place of the Hamiltonian (and with no collapse operators)

In [7]:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
c_ops = [sqrt(kappa) * a, sqrt(gamma) * sm]

L = liouvillian(H, c_ops)

result = mesolve(L, psi0, tlist, [], [a.dag() * a, sm.dag() * sm])

visualize(result)


### Pass liouvillian instead of collapse operators¶

Works with list-function and list-string time-dependencies too:

In [8]:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
c_ops = [sqrt(kappa) * a, sqrt(gamma) * sm]

L = liouvillian(None, c_ops)

result = mesolve(H, psi0, tlist, [[L, '1.0']], [a.dag() * a, sm.dag() * sm])

visualize(result)

In [9]:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
c_ops = [sqrt(kappa) * a, sqrt(gamma) * sm]

L = liouvillian(None, c_ops)

result = mesolve(H, psi0, tlist, [[L, lambda t, args: 1.0]], [a.dag() * a, sm.dag() * sm])

visualize(result)


### Liouvillian instead of a list of collapse operators¶

In [10]:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
c_ops = [sqrt(kappa) * a, sqrt(gamma) * sm]

L = liouvillian(None, c_ops)

result = mesolve(H, psi0, tlist, L, [a.dag() * a, sm.dag() * sm])

visualize(result)


### Mix normal collapse operators and Liouvillians, passed as the collapse-operator argument¶

In [11]:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
c_ops = [sqrt(kappa) * a, liouvillian(None, [sqrt(gamma) * sm])]

result = mesolve(H, psi0, tlist, c_ops, [a.dag() * a, sm.dag() * sm])

visualize(result)


### Make a custom non-standard Liouvillian¶

In [12]:
psi0 = tensor(basis(N,0), basis(2,1))    # start with an excited atom

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

H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
c_ops = [sqrt(kappa) * a, sqrt(gamma) * sm]

L = sum([spre(c)*spost(c.dag()) - 0.25 * spre(c.dag()*c) - 0.5 * spost(c.dag()*c) for c in c_ops])

result = mesolve(H, psi0, tlist, L, [a.dag() * a, sm.dag() * sm])

visualize(result)


### Software versions¶

In [13]:
from qutip.ipynbtools import version_table
version_table()

Out[13]:
SoftwareVersion
Cython0.16
SciPy0.10.1
QuTiP2.2.0.dev-793af74
Python2.7.3 (default, Sep 26 2012, 21:51:14) [GCC 4.7.2]
IPython0.13
OSposix [linux2]
Numpy1.7.0.dev-3f45eaa
matplotlib1.3.x
Mon Feb 18 16:06:16 2013