Speed test fo Milstein and Euler methods

Based on development-smesolve-milstein notebook

Denis V. Vasilyev

30 september 2013


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

Single stochastic operator


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))


Total run time:  14.85s

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))


Total run time:   3.47s

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))


Total run time:  14.88s

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))


Total run time:   4.46s

Variance $V_{\mathrm{c}}$ as a function of time


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();


Two stochastic operators


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))


Total run time:  28.98s

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))


Total run time:   3.99s

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))


Total run time:  43.02s

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))


Total run time:   7.52s

Variance $V_{\mathrm{c}}$ as a function of time


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();


Multiple stochastic operators


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))


Total run time:  56.55s

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))


Total run time:   4.70s

Software versions


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]:
SoftwareVersion
Python3.3.2+ (default, Oct 9 2013, 14:50:09) [GCC 4.8.1]
IPython2.0.0-dev
OSposix [linux]
Cython0.20.post0
numpy1.9.0.dev-d4c7c3a
matplotlib1.4.x
scipy0.15.0.dev-c04c506
qutip3.0.0.dev-8072d41
Fri Mar 07 16:26:14 2014 JST