QuTiP example: compare time-dependence formats

J.R. Johansson and P.D. Nation

For more information about QuTiP see http://qutip.org

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 [ ]: