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.
%pylab inline
from qutip import *
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)
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)
psi0 = tensor(basis(N,0), basis(2,1))
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]
result = mesolve(H, psi0, tlist, c_ops, [a.dag() * a, sm.dag() * sm])
visualize(result)
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)
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)
Works with list-function and list-string time-dependencies too:
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)
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)
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)
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)
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)
from qutip.ipynbtools import version_table
version_table()
Software | Version |
---|---|
Cython | 0.16 |
SciPy | 0.10.1 |
QuTiP | 2.2.0.dev-793af74 |
Python | 2.7.3 (default, Sep 26 2012, 21:51:14) [GCC 4.7.2] |
IPython | 0.13 |
OS | posix [linux2] |
Numpy | 1.7.0.dev-3f45eaa |
matplotlib | 1.3.x |
Mon Feb 18 16:06:16 2013 |