In [1]:
%pylab inline
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)
# Righthand side for the Milstein method for a homodyne detection scheme
def rhs_milstein(L, rho_t, t, A_ops, dt, dW, d1, d2, args):
drho_t = spmv(L.data,
L.indices,
L.indptr, rho_t) * dt
A = A_ops[0]
M = A[0] + A[3]
e1 = cy_expect_rho_vec(M, rho_t)
d1_vec = spmv(A[7].data, A[7].indices, A[7].indptr, rho_t)
d2_vec = spmv(M.data, M.indices, M.indptr, rho_t)
d2_vec2 = spmv(M.data, M.indices, M.indptr, d2_vec)
e2 = cy_expect_rho_vec(M, d2_vec)
return rho_t + drho_t + d1_vec*dt + (d2_vec - e1*rho_t)*dW[0,0] + \
0.5 * (d2_vec2 - 2*e1*d2_vec + (-e2 + 2*e1*e1)*rho_t)*(dW[0,0]*dW[0,0] - dt)
In [5]:
#Default Euler-Maruyama
sol_eul = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
nsubsteps=Nsub, method='homodyne',
options=Odeoptions(store_states=True, average_states=False))
In [6]:
#Fast Euler-Maruyama
sol_eul_fast = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
nsubsteps=Nsub, method='homodyne', solver='euler-maruyama_fast',
options=Odeoptions(store_states=True, average_states=False))
In [7]:
#Default 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))
In [8]:
#Fast Milstein
sol_mil_fast = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
nsubsteps=Nsub, method='homodyne', solver='milstein_fast',
options=Odeoptions(store_states=True, average_states=False))
In [9]:
fig, ax = subplots()
ax.plot(tlist,sol_eul.expect[1] - abs(sol_eul.expect[0])**2, label='Euler')
ax.plot(tlist,sol_eul_fast.expect[1] - abs(sol_eul_fast.expect[0])**2, label='Euler-fast')
ax.plot(tlist,sol_mil_fast.expect[1] - abs(sol_mil_fast.expect[0])**2, label='Milstein-fast')
ax.plot(tlist,Vc*ones(NN), label='exact steady state solution')
ax.plot(tlist,y, label="exact solution")
ax.legend();
In [10]:
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 [11]:
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 [12]:
#Default Euler-Maruyama
sol_eul = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
nsubsteps=Nsub, method='homodyne',
options=Odeoptions(store_states=True, average_states=False))
In [13]:
#Fast Euler-Maruyama
sol_eul_fast = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
nsubsteps=Nsub, method='homodyne', solver='euler-maruyama_fast',
options=Odeoptions(store_states=True, average_states=False))
In [14]:
#Default 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))
In [15]:
#Fast Milstein
sol_mil_fast = smesolve(H, rho0, tlist, c_op, sc_op, e_op,
nsubsteps=Nsub, method='homodyne', solver='milstein_fast',
options=Odeoptions(store_states=True, average_states=False))
In [16]:
fig, ax = subplots()
ax.plot(tlist,sol_eul.expect[1] - abs(sol_eul.expect[0])**2, label='Euler')
ax.plot(tlist,sol_eul_fast.expect[1] - abs(sol_eul_fast.expect[0])**2, label='Euler-fast')
ax.plot(tlist,sol_mil_fast.expect[1] - abs(sol_mil_fast.expect[0])**2, label='Milstein-fast')
ax.plot(tlist,VcEps*ones(NN), label='Exact steady state solution')
ax.plot(tlist,y, label='Exact solution')
ax.legend();
In [17]:
#Default Euler-Maruyama
sol_eul = smesolve(H, rho0, tlist, c_op, sc_op + sc_op, e_op,
nsubsteps=Nsub, method='homodyne',
options=Odeoptions(store_states=True, average_states=False))
In [18]:
#Fast Euler-Maruyama
sol_eul_fast = smesolve(H, rho0, tlist, c_op, sc_op + sc_op, e_op,
nsubsteps=Nsub, method='homodyne', solver='euler-maruyama_fast',
options=Odeoptions(store_states=True, average_states=False))
In [19]:
#%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
In [20]:
%load_ext version_information
In [21]:
%version_information Cython, numpy, matplotlib, scipy, qutip
Out[21]: