# QuTiP example: compare time-dependence formats¶

J.R. Johansson and P.D. Nation

In [3]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib

In [4]:
from qutip import *
import time


## System: driven qubit¶

Example comparing master equation with constant and time-dependent collapse Hamiltonian and collapse operators.

In [5]:
def gamma1_t(t, args):
""" a time-dependent rate """
return args['G1'] * exp(-0.1*t)

def gamma2_t(t, args):
""" a time-dependent rate """
return args['G2'] * exp(-0.2*t)

def H1_coeff_t(t, args):
return sin(args['w'] * t)

def hamiltonian_t(t, args):
""" evaluate the hamiltonian at time t. """
H0 = args[0]
H1 = args[1]
w  = args[2]
return H0 + H1 * sin(w * t)

In [6]:
def qubit_integrate(delta, eps0, A, w, gamma1, gamma2, psi0, tlist):

# Hamiltonian
sx = sigmax()
sz = sigmaz()
sm = destroy(2)

H0 = - delta/2.0 * sx - eps0/2.0 * sz
H1 = - A * sx

# --------------------------------------------------------------------------
# 1) time-dependent hamiltonian and collapse operators: using list-function
#    format
#

H = [H0, [H1, H1_coeff_t]]
args = {'w': w, 'G1': gamma1, 'G2': gamma2}

c_op_list = []
c_op_list.append([sm, gamma1_t]) # relaxation
c_op_list.append([sz, gamma2_t]) # dephasing
#c_op_list.append(sqrt(gamma1) * sm) # relaxation
#c_op_list.append(sqrt(gamma2) * sz) # dephasing

start_time = time.time()
output = mesolve(H, psi0, tlist, c_op_list, [sm.dag() * sm], args=args)
expt_list1 = output.expect
print('Method 1: time elapsed = ' + str(time.time() - start_time))

# --------------------------------------------------------------------------
# 2) time-dependent hamiltonian and collapse operators: using list-string
#    format
#

H = [H0, [H1, 'sin(w * t)']]
args = {'w': w, 'G1': gamma1, 'G2': gamma2}

c_op_list = []
c_op_list.append([sm, 'G1 * exp(-0.1*t)']) # relaxation
c_op_list.append([sz, 'G2 * exp(-0.2*t)']) # dephasing
#c_op_list.append(sqrt(gamma1) * sm) # relaxation
#c_op_list.append(sqrt(gamma2) * sz) # dephasing

start_time = time.time()
output = mesolve(H, psi0, tlist, c_op_list, [sm.dag() * sm], args=args)
expt_list2 = output.expect
print('Method 2: time elapsed = ' + str(time.time() - start_time))

# --------------------------------------------------------------------------
# 3) time-dependent hamiltonian and but time-independent collapse operators:
#    using function callback format
#

args = [H0, H1, w]

c_op_list = []
c_op_list.append(sqrt(gamma1) * sm) # relaxation
c_op_list.append(sqrt(gamma2) * sz) # dephasing

start_time = time.time()
output = mesolve(hamiltonian_t, psi0, tlist, c_op_list, [sm.dag() * sm], args=args)
expt_list3 = output.expect
print('Method 3: time elapsed = ' + str(time.time() - start_time))

# --------------------------------------------------------------------------
# 4) Constant hamiltonian and collapse operators
#

H_rwa = - delta/2.0 * sx - A * sx / 2

c_op_list = []
c_op_list.append(sqrt(gamma1) * sm) # relaxation
c_op_list.append(sqrt(gamma2) * sz) # dephasing

start_time = time.time()
output = mesolve(H_rwa, psi0, tlist, c_op_list, [sm.dag() * sm])
expt_list4 = output.expect
print('Method 4: time elapsed = ' + str(time.time() - start_time))

# --------------------------------------------------------------------------
# 5) Unitary evolution, constant hamiltonian
#

c_op_list = []

start_time = time.time()
output = mesolve(H_rwa, psi0, tlist, c_op_list, [sm.dag() * sm])
expt_list5 = output.expect
print('Method 5: time elapsed = ' + str(time.time() - start_time))

return expt_list1[0], expt_list2[0], expt_list3[0], expt_list4[0], expt_list5[0]

In [7]:
delta = 0.0 * 2 * pi   # qubit sigma_x coefficient
eps0  = 1.0 * 2 * pi   # qubit sigma_z coefficient
A     = 0.25 * 2 * pi  # driving amplitude (reduce to make the RWA more accurate)
w     = 1.0 * 2 * pi   # driving frequency
gamma1 = 0.2           # relaxation rate
gamma2 = 0.1           # dephasing  rate
psi0 = basis(2,1)      # initial state

tlist = linspace(0, 5.0 * 2 * pi / A, 250)

In [8]:
p_ex1, p_ex2, p_ex3, p_ex4, p_ex5 = qubit_integrate(delta, eps0, A, w, gamma1, gamma2, psi0, tlist)

Method 1: time elapsed = 0.824520111084
Method 2: time elapsed = 1.90881586075
Method 3: time elapsed = 0.380300045013
Method 4: time elapsed = 0.0217039585114
Method 5: time elapsed = 0.213099956512

In [12]:
figure(figsize=(12, 8))
plot(tlist, real(p_ex1), 'b',
tlist, real(p_ex2), 'g.',
tlist, real(p_ex3), 'r.-',
tlist, real(p_ex4), 'c--',
tlist, real(p_ex5), 'm')
xlabel('Time')
ylabel('Occupation probability')
title('Excitation probabilty of qubit')
legend(("Time-dependent function format",
"Time-dependent string format",
"Time-dependent func format",
"Const. collapse operators", "Unitary evolution"));


## System: Driven oscillator¶

In [13]:
def gamma_t(t, args):
""" a time-dependent rate """
return args['G'] * (1 + sin(args['w']*t)) * exp(-0.5*t)

def H1_coeff_t(t, args):
return sin(args['w'] * t)

def hamiltonian_t(t, args):
""" evaluate the hamiltonian at time t. """
H0 = args[0]
H1 = args[1]
w  = args[2]
return H0 + H1 * sin(w * t)

In [14]:
def oscillator_integrate(solver, N, w0, A, w, gamma, rho0, tlist):

# Hamiltonian: oscillator + driving on x
H0 = w0 * 2 * pi * num(N)
H1 = A * (destroy(N).dag() + destroy(N))

# --------------------------------------------------------------------------
# 1) time-dependent hamiltonian and collapse operators: using list-function
#    format
#

H = [H0, [H1, H1_coeff_t]]
args = {'w': w, 'G': gamma}

c_op_list = []
c_op_list.append([destroy(N), gamma_t]) # relaxation

start_time = time.time()
output = solver(H, rho0, tlist, c_op_list, [num(N)], args=args)
expt_list1 = output.expect
print('Method 1: time elapsed = ' + str(time.time() - start_time))

# --------------------------------------------------------------------------
# 2) time-dependent hamiltonian and collapse operators: using list-string
#    format
#

H = [H0, [H1, 'sin(w * t)']]
args = {'w': w, 'G': gamma}

c_op_list = []
c_op_list.append([destroy(N), 'G * (1 + sin(w*t))']) # relaxation

start_time = time.time()
output = solver(H, rho0, tlist, c_op_list, [num(N)], args=args)
expt_list2 = output.expect
print('Method 2: time elapsed = ' + str(time.time() - start_time))

# --------------------------------------------------------------------------
# 3) time-dependent hamiltonian and but time-independent collapse operators:
#    using function callback format
#
args = [H0, H1, w]

c_op_list = []
c_op_list.append(sqrt(gamma) * destroy(N)) # relaxation

start_time = time.time()
output = solver(hamiltonian_t, rho0, tlist, c_op_list, [num(N)], args=args)
expt_list3 = output.expect
print('Method 3: time elapsed = ' + str(time.time() - start_time))

# --------------------------------------------------------------------------
# 4) Constant hamiltonian and collapse operators
#
c_op_list = []
c_op_list.append(sqrt(gamma) * destroy(N)) # relaxation

start_time = time.time()
output = solver(H0, rho0, tlist, c_op_list, [num(N)])
expt_list4 = output.expect
print('Method 4: time elapsed = ' + str(time.time() - start_time))

# --------------------------------------------------------------------------
# 5) Unitary evolution, constant hamiltonian
#
start_time = time.time()
output = solver(H0, rho0, tlist, [], [num(N)])
expt_list5 = output.expect
print('Method 5: time elapsed = ' + str(time.time() - start_time))

return expt_list1[0], expt_list2[0], expt_list3[0], expt_list4[0], expt_list5[0]

In [15]:
N     = 10
w0    = 1.0 * 2 * pi   # oscillator energy
A     = 0.15 * 2 * pi  # driving amplitude (reduce to make the RWA more accurate)
w     = 1.0 * 2 * pi   # driving frequency
gamma = 0.2            # relaxation rate
rho0 = fock_dm(N,N-2)  # initial state

tlist = linspace(0, 10.0, 250)

In [17]:
solver = mesolve  #mesolve or mcsolve
p_ex1, p_ex2, p_ex3, p_ex4, p_ex5 = oscillator_integrate(solver, N, w0, A, w, gamma, rho0, tlist)

Method 1: time elapsed = 5.66927599907
Method 2: time elapsed = 1.95071101189
Method 3: time elapsed = 4.07780599594
Method 4: time elapsed = 0.020231962204
Method 5: time elapsed = 0.012766122818

In [21]:
figure(figsize=(12, 8))
plot(tlist, real(p_ex1), 'b',
tlist, real(p_ex2), 'g.',
tlist, real(p_ex3), 'r.-',
tlist, real(p_ex4), 'c--',
tlist, real(p_ex5), 'm')
xlabel('Time')
ylabel('Occupation probability')
title('Excitation probabilty of qubit')
legend(("Time-dependent function format",
"Time-dependent string format",
"Time-dependent func format. Const collapse ops.",
"Const. collapse operators", "Unitary evolution"));


## Versions¶

In [1]:
from qutip.ipynbtools import version_table

version_table()

Out[1]:
SoftwareVersion
Cython0.19
SciPy0.14.0.dev-2a4ba40
QuTiP2.3.0.dev-0fd8af4
Python2.7.4 (default, Apr 19 2013, 18:28:01) [GCC 4.7.3]
IPython2.0.0-dev
OSposix [linux2]
Numpy1.8.0.dev-928289b
matplotlib1.4.x
Sat Sep 28 23:14:28 2013 JST
In [ ]: