In [1]:

    

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

from qutip import *
from numpy import log2, cos, sin
from scipy.integrate import odeint
from qutip.cy.spmatfuncs import cy_expect_rho_vec, spmv


    

In [3]:

    

th = 0.1 # Interaction parameter
alpha = cos(th)
beta = sin(th)
gamma = 1

# Exact steady state solution for Vc
Vc = (alpha*beta - gamma + sqrt((gamma-alpha*beta)**2 + 4*gamma*alpha**2))/(4*alpha**2)

#********* Model ************
NN = 200
tlist = linspace(0,50,NN)
Nsub = 200
N = 20
Id = qeye(N)
a = destroy(N)
s = 0.5*((alpha+beta)*a + (alpha-beta)*a.dag())
x = (a + a.dag())/sqrt(2)
H = Id
c_op = [sqrt(gamma)*a]
sc_op = [s]
e_op = [x, x*x]
rho0 = fock_dm(N,0) #initial vacuum state


    

In [4]:

    

# Solution of the differential equation for the variance Vc
y0 = 0.5
def func(y, t):
    return -(gamma - alpha*beta)*y - 2*alpha*alpha*y*y + 0.5*gamma
y = odeint(func, y0, tlist)


    

In [5]:

    

# list solver
help(stochastic_solvers)


    

    




Help on function stochastic_solvers in module qutip.stochastic:

stochastic_solvers()
    Available solvers for ssesolve and smesolve
        euler-maruyama:
            A simple generalization of the Euler method for ordinary
            differential equations to stochastic differential equations.
            Only solver which could take non-commuting sc_ops. *not tested*
            -Order 0.5
            -Code: 'euler-maruyama', 'euler', 0.5
    
        milstein, Order 1.0 strong Taylor scheme:
            Better approximate numerical solution to stochastic
            differential equations.
            -Order strong 1.0
            -Code: 'milstein', 1.0
            Numerical Solution of Stochastic Differential Equations
            Chapter 10.3 Eq. (3.1), By Peter E. Kloeden, Eckhard Platen
    
        milstein-imp, Order 1.0 implicit strong Taylor scheme:
            Implicit milstein scheme for the numerical simulation of stiff
            stochastic differential equations.
            -Order strong 1.0
            -Code: 'milstein-imp'
            Numerical Solution of Stochastic Differential Equations
            Chapter 12.2 Eq. (2.9), By Peter E. Kloeden, Eckhard Platen
    
        predictor-corrector:
            Generalization of the trapezoidal method to stochastic
            differential equations. More stable than explicit methods.
            -Order strong 0.5, weak 1.0
            Only the stochastic part is corrected.
                (alpha = 0, eta = 1/2)
                -Code: 'pred-corr', 'predictor-corrector', 'pc-euler'
            Both the deterministic and stochastic part corrected.
                (alpha = 1/2, eta = 1/2)
                -Code: 'pc-euler-imp', 'pc-euler-2', 'pred-corr-2'
            Numerical Solution of Stochastic Differential Equations
            Chapter 15.5 Eq. (5.4), By Peter E. Kloeden, Eckhard Platen
    
        platen:
            Explicit scheme, create the milstein using finite difference instead of
            derivatives. Also contain some higher order terms, thus converge better
            than milstein while staying strong order 1.0.
            Do not require derivatives, therefore usable for
            :func:`qutip.stochastic.general_stochastic`
            -Order strong 1.0, weak 2.0
            -Code: 'platen', 'platen1', 'explicit1'
            The Theory of Open Quantum Systems
            Chapter 7 Eq. (7.47), H.-P Breuer, F. Petruccione
    
        rouchon:
            Scheme keeping the positivity of the density matrix. (smesolve only)
            -Order strong 1.0?
            -Code: 'rouchon', 'Rouchon'
            Efficient Quantum Filtering for Quantum Feedback Control
            Pierre Rouchon, Jason F. Ralph
            arXiv:1410.5345 [quant-ph]
    
        taylor1.5, Order 1.5 strong Taylor scheme:
            Solver with more terms of the Ito-Taylor expansion.
            Default solver for smesolve and ssesolve.
            -Order strong 1.5
            -Code: 'taylor1.5', 'taylor15', 1.5, None
            Numerical Solution of Stochastic Differential Equations
            Chapter 10.4 Eq. (4.6), By Peter E. Kloeden, Eckhard Platen
    
        taylor1.5-imp, Order 1.5 implicit strong Taylor scheme:
            implicit Taylor 1.5 (alpha = 1/2, beta = doesn't matter)
            -Order strong 1.5
            -Code: 'taylor1.5-imp', 'taylor15-imp'
            Numerical Solution of Stochastic Differential Equations
            Chapter 12.2 Eq. (2.18), By Peter E. Kloeden, Eckhard Platen
    
        explicit1.5, Explicit Order 1.5 Strong Schemes:
            Reproduce the order 1.5 strong Taylor scheme using finite difference
            instead of derivatives. Slower than taylor15 but usable by
            :func:`qutip.stochastic.general_stochastic`
            -Order strong 1.5
            -Code: 'explicit1.5', 'explicit15', 'platen15'
            Numerical Solution of Stochastic Differential Equations
            Chapter 11.2 Eq. (2.13), By Peter E. Kloeden, Eckhard Platen
    
        taylor2.0, Order 2 strong Taylor scheme:
            Solver with more terms of the Stratonovich expansion.
            -Order strong 2.0
            -Code: 'taylor2.0', 'taylor20', 2.0
            Numerical Solution of Stochastic Differential Equations
            Chapter 10.5 Eq. (5.2), By Peter E. Kloeden, Eckhard Platen
    
        ---All solvers, except taylor2.0, are usable in both smesolve and ssesolve
        and for both heterodyne and homodyne. taylor2.0 only work for 1 stochastic
        operator not dependent of time with the homodyne method.
        The :func:`qutip.stochastic.general_stochastic` only accept derivatives
        free solvers: ['euler', 'platen', 'explicit1.5'].
    
    Available solver for photocurrent_sesolve and photocurrent_mesolve:
            Photocurrent use ordinary differential equations between
            stochastic "jump/collapse".
        euler:
            Euler method for ordinary differential equations.
            Default solver
            -Order 1.0
            -Code: 'euler'
    
        predictor–corrector:
            predictor–corrector method (PECE) for ordinary differential equations.
            -Order 2.0
            -Code: 'pred-corr'

In [6]:

    

# Euler-Maruyama
sol_euler = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
                   nsubsteps=Nsub, method='homodyne', solver='euler-maruyama',
                   options=Odeoptions(store_states=True, average_states=False))


    

    




Total run time:   1.09s

In [7]:

    

# Milstein
sol_mil = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
                   nsubsteps=Nsub, method='homodyne', solver='milstein',
                   options=Odeoptions(store_states=True, average_states=False))


    

    




Total run time:   1.55s

In [8]:

    

# predictor-corrector
sol_pc_euler = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
                   nsubsteps=Nsub, method='homodyne', solver='pc-euler',
                   options=Odeoptions(store_states=True, average_states=False))


    

    




Total run time:   2.42s

In [9]:

    

# predictor-corrector-2
sol_pc_euler2 = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
                   nsubsteps=Nsub, method='homodyne', solver='pred-corr-2',
                   options=Odeoptions(store_states=True, average_states=False))


    

    




Total run time:   3.77s

In [10]:

    

# Default: taylor1.5
sol_taylor15 = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
                   nsubsteps=Nsub, method='homodyne', solver='taylor1.5',
                   options=Odeoptions(store_states=True, average_states=False))


    

    




Total run time:   7.01s

In [11]:

    

# Taylor 2.0
sol_taylor20 = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
                   nsubsteps=Nsub, method='homodyne', solver='taylor2.0',
                   options=Odeoptions(store_states=True, average_states=False))


    

    




Total run time:   8.04s

In [12]:

    

plot(tlist,sol_euler.expect[1] - abs(sol_euler.expect[0])**2, label='Euler')
plot(tlist,sol_mil.expect[1] - abs(sol_mil.expect[0])**2, label='Milstein')
plot(tlist,sol_pc_euler.expect[1] - abs(sol_pc_euler.expect[0])**2, label='pc-euler')
plot(tlist,sol_pc_euler2.expect[1] - abs(sol_pc_euler2.expect[0])**2, label='pc-euler-2')
plot(tlist,sol_taylor15.expect[1] - abs(sol_taylor15.expect[0])**2, label='taylor1.5')
plot(tlist,sol_taylor20.expect[1] - abs(sol_taylor20.expect[0])**2, label='taylor2.0')
plot(tlist,Vc*ones(NN),"k:", label='exact steady state solution')
plot(tlist,y,"k", label="exact solution")
legend();
show()

#plot(tlist,sol_euler.expect[1] - abs(sol_euler.expect[0])**2, label='Euler')
plot(tlist,sol_mil.expect[1] - abs(sol_mil.expect[0])**2, label='Milstein')
plot(tlist,sol_pc_euler.expect[1] - abs(sol_pc_euler.expect[0])**2, label='pc-euler')
plot(tlist,sol_pc_euler2.expect[1] - abs(sol_pc_euler2.expect[0])**2, label='pc-euler-2')
plot(tlist,sol_taylor15.expect[1] - abs(sol_taylor15.expect[0])**2, label='taylor1.5')
plot(tlist,sol_taylor20.expect[1] - abs(sol_taylor20.expect[0])**2, label='taylor2.0')
plot(tlist,Vc*ones(NN),"k:", label='exact steady state solution')
plot(tlist,y,"k", label="exact solution")
legend();
xlim([0,25])
ylim([0.3238,0.325])
show()

#plot(tlist,sol_euler.expect[1] - abs(sol_euler.expect[0])**2, label='Euler')
#plot(tlist,sol_mil.expect[1] - abs(sol_mil.expect[0])**2, label='Milstein')
#plot(tlist,sol_pc_euler.expect[1] - abs(sol_pc_euler.expect[0])**2, label='pc-euler')
#plot(tlist,sol_pc_euler2.expect[1] - abs(sol_pc_euler2.expect[0])**2, label='pc-euler-2')
plot(tlist,sol_taylor15.expect[1] - abs(sol_taylor15.expect[0])**2, label='taylor1.5')
plot(tlist,sol_taylor20.expect[1] - abs(sol_taylor20.expect[0])**2, label='taylor2.0')
plot(tlist,Vc*ones(NN),"k:", label='exact steady state solution')
plot(tlist,y,"k", label="exact solution")
legend();
xlim([0,25])
ylim([0.3241,0.32418])
show()


    

In [13]:

    

th = 0.1
alpha = cos(th)
beta = sin(th)
gamma = 1
eps = 0.3

VcEps = ((1-2*eps)*alpha*beta - gamma + \
         sqrt((gamma-alpha*beta)**2 + 4*gamma*alpha*((1-eps)*alpha + eps*beta)))/(4*(1-eps)*alpha**2)
UcEps = (-(1-2*eps)*alpha*beta - gamma + \
         sqrt((gamma-alpha*beta)**2 + 4*eps*beta*gamma*(beta-alpha)))/(4*eps*beta**2)

NN = 200
tlist = linspace(0,3,NN)
Nsub = 200
N = 20
Id = qeye(N)
a = destroy(N)
s = 0.5*((alpha+beta)*a + (alpha-beta)*a.dag())
x = (a + a.dag())/sqrt(2)
H = Id
c_op = [sqrt(gamma)*a]
sc_op = [sqrt(1-eps)*s, sqrt(eps)*1j*s]
e_op = [x, x*x]
rho0 = fock_dm(N,0)
y0 = 0.5


    

In [14]:

    

def func(y, t):
    return -(gamma - (1-2*eps)*alpha*beta)*y - 2*(1-eps)*alpha*alpha*y*y + 0.5*(gamma + eps*beta*beta)
y = odeint(func, y0, tlist)

def funcZ(z, t):
    return -(gamma + (1-2*eps)*alpha*beta)*z - 2*eps*beta*beta*z*z + 0.5*(gamma + (1-eps)*alpha*alpha)
z = odeint(funcZ, y0, tlist)


    

In [15]:

    

# Euler-Maruyama
sol_euler = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
                   nsubsteps=Nsub, method='homodyne', solver='euler-maruyama',
                   options=Odeoptions(store_states=True, average_states=False))


    

    




Total run time:   1.49s

In [16]:

    

# Milstein
sol_mil = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
                   nsubsteps=Nsub, method='homodyne', solver='milstein',
                   options=Odeoptions(store_states=True, average_states=False))


    

    




Total run time:   3.05s

In [17]:

    

# predictor-corrector
sol_pc_euler = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
                   nsubsteps=Nsub, method='homodyne', solver='pc-euler',
                   options=Odeoptions(store_states=True, average_states=False))


    

    




Total run time:   4.21s

In [18]:

    

# predictor-corrector-2
sol_pc_euler2 = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
                   nsubsteps=Nsub, method='homodyne', solver='pred-corr-2',
                   options=Odeoptions(store_states=True, average_states=False))


    

    




Total run time:   7.26s

In [19]:

    

# Default: taylor1.5
sol_taylor15 = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
                   nsubsteps=Nsub, method='homodyne', solver='taylor1.5',
                   options=Odeoptions(store_states=True, average_states=False))


    

    




Total run time:  19.44s

In [20]:

    

plot(tlist,sol_euler.expect[1] - abs(sol_euler.expect[0])**2, label='Euler')
plot(tlist,sol_mil.expect[1] - abs(sol_mil.expect[0])**2, label='Milstein')
plot(tlist,sol_pc_euler.expect[1] - abs(sol_pc_euler.expect[0])**2, label='pc-euler')
plot(tlist,sol_pc_euler2.expect[1] - abs(sol_pc_euler2.expect[0])**2, label='pc-euler-2')
plot(tlist,sol_taylor15.expect[1] - abs(sol_taylor15.expect[0])**2, label='taylor1.5')
plot(tlist,Vc*ones(NN), label='exact steady state solution')
plot(tlist,y, label="exact solution")
legend();
show()

plot(tlist,sol_euler.expect[1] - abs(sol_euler.expect[0])**2, label='Euler')
plot(tlist,sol_mil.expect[1] - abs(sol_mil.expect[0])**2, label='Milstein')
plot(tlist,sol_pc_euler.expect[1] - abs(sol_pc_euler.expect[0])**2, label='pc-euler')
plot(tlist,sol_pc_euler2.expect[1] - abs(sol_pc_euler2.expect[0])**2, label='pc-euler-2')
plot(tlist,sol_taylor15.expect[1] - abs(sol_taylor15.expect[0])**2, label='taylor1.5')
plot(tlist,Vc*ones(NN), label='exact steady state solution')
plot(tlist,y, label="exact solution")
legend();
xlim([1.5,3.0])
ylim([0.347,0.356])


    

    













    
Out[20]:





(0.347, 0.356)






    

In [21]:

    

plot(tlist,sol_euler.expect[1] - abs(sol_euler.expect[0])**2-y.transpose()[0], label='Euler')
plot(tlist,sol_mil.expect[1] - abs(sol_mil.expect[0])**2-y.transpose()[0], label='Milstein')
plot(tlist,sol_pc_euler.expect[1] - abs(sol_pc_euler.expect[0])**2-y.transpose()[0], label='pc-euler')
plot(tlist,sol_pc_euler2.expect[1] - abs(sol_pc_euler2.expect[0])**2-y.transpose()[0], label='pc-euler-2')
plot(tlist,sol_taylor15.expect[1] - abs(sol_taylor15.expect[0])**2-y.transpose()[0], label='taylor1.5')
legend()
show()

plot(tlist,sol_mil.expect[1] - abs(sol_mil.expect[0])**2-y.transpose()[0], label='Milstein')
plot(tlist,sol_pc_euler.expect[1] - abs(sol_pc_euler.expect[0])**2-y.transpose()[0], label='pc-euler')
plot(tlist,sol_pc_euler2.expect[1] - abs(sol_pc_euler2.expect[0])**2-y.transpose()[0], label='pc-euler-2')
plot(tlist,sol_taylor15.expect[1] - abs(sol_taylor15.expect[0])**2-y.transpose()[0], label='taylor1.5')
legend()
show()


    

In [22]:

    

sc_ops_multiple = [sc_op[0]*0.5**0.5, sc_op[0]*0.5**0.5]
#Euler-Maruyama
sol_eul = smesolve(H, rho0, tlist, c_op, sc_ops_multiple, e_op,
                   nsubsteps=Nsub, method='homodyne', solver='euler',
                   options=Odeoptions(store_states=True, average_states=False))


    

    




Total run time:   1.74s

In [23]:

    

#Milstein
sol_milstein = smesolve(H, rho0, tlist, c_op, sc_ops_multiple, e_op,
                   nsubsteps=Nsub, method='homodyne', solver='milstein',
                   options=Odeoptions(store_states=True, average_states=False))


    

    




Total run time:   3.56s

In [24]:

    

#taylor1.5
sol_taylor = smesolve(H, rho0, tlist, c_op, sc_ops_multiple, e_op,
                   nsubsteps=Nsub, method='homodyne', solver='taylor1.5',
                   options=Odeoptions(store_states=True, average_states=False))


    

    




Total run time:  14.48s

In [25]:

    

plot(tlist,sol_eul.expect[1] - abs(sol_eul.expect[0])**2, label='Euler')
plot(tlist,sol_milstein.expect[1] - abs(sol_milstein.expect[0])**2, label='Milstein')
plot(tlist,sol_taylor.expect[1] - abs(sol_taylor.expect[0])**2, label='taylor1.5')
plot(tlist,Vc*ones(NN), label='exact steady state solution')
plot(tlist,y, label="exact solution")
legend()
show()


    

In [26]:

    

from qutip.ipynbtools import version_table

version_table()


    

    
Out[26]:




Software
Version
QuTiP
4.4.0.dev0+21a52e95
Numpy
1.16.0
SciPy
1.2.0
matplotlib
3.0.2
Cython
0.29.3
Number of CPUs
2
BLAS Info
OPENBLAS
IPython
7.2.0
Python
3.6.7 (default, Oct 22 2018, 11:32:17) 
[GCC 8.2.0]
OS
posix [linux]
Fri Feb 01 13:36:16 2019 JST

In [ ]: