In [1]:
import cutiepy
cutiepy.interactive.INTERACTIVE = False
import qutip
import numpy as np
import scipy
ts = np.linspace(0,10000*np.pi,10000)

qutip


In [2]:
init = qutip.basis(2,0)
H = qutip.sigmax()
opt = qutip.Options(rhs_reuse=True)
%time res = qutip.sesolve(H,init,ts,[], options=opt)
%time res = qutip.sesolve(H,init,ts,[], options=opt)
print(res.states[-1].data.A, res.states[-1].norm())


CPU times: user 9.03 s, sys: 20 ms, total: 9.05 s
Wall time: 9.03 s
CPU times: user 8.72 s, sys: 12 ms, total: 8.73 s
Wall time: 8.71 s
[[ 0.99986301+0.j        ]
 [ 0.00000000-0.01655151j]] 1.0

cutiepy


In [3]:
H = cutiepy.sigmax()
init = cutiepy.basis(2,0)
%time res = cutiepy.sesolve(H, init, ts)
%time res = cutiepy.sesolve(H, init, ts)
print(res[-1].numerical, np.linalg.norm(res[-1].numerical))


CPU times: user 3.39 s, sys: 136 ms, total: 3.52 s
Wall time: 7.2 s
CPU times: user 1.83 s, sys: 76 ms, total: 1.91 s
Wall time: 1.88 s
[[ 0.99411803+0.j        ]
 [ 0.00000000-0.02062065j]] 0.99433186891

scipy (dense)


In [4]:
%%time
H = np.array([[0,1],[1,0]],dtype='complex128')
op = -1j*H
f = lambda t, state: op.dot(state)
ode = scipy.integrate.ode(f)                                             
ode.set_integrator('zvode', method='adams', rtol=1e-6, atol=1e-8, nsteps=2000)     
state0_num = np.array([1,0],dtype='complex128')                                             
ode.set_initial_value(state0_num, ts[0])                                     
res = [state0_num]                                                              
for t in ts[1:]:                                                             
    ode.integrate(t)                                                            
    if not ode.successful():                                                    
        RuntimeError('The ODE solver failed.') 
    #y = ode.y                    
    #n = scipy.linalg.norm(y)     
    #y /= n                       
    #ode.set_initial_value(y, t)  
    res.append(ode.y)                                                           
print(res[-1], np.linalg.norm(res[-1]))


[ 0.99158045+0.j          0.00000000-0.01967444j] 0.991775620659
CPU times: user 1.29 s, sys: 0 ns, total: 1.29 s
Wall time: 1.28 s