In [1]:

    

import cutiepy
cutiepy.interactive.INTERACTIVE = False
import qutip
import numpy as np
import scipy
N = 100
ts = np.linspace(0,1000,100)


    

In [2]:

    

init = qutip.coherent(N,2)
H = qutip.num(N)
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[:4], res.states[-1].norm())


    

    




CPU times: user 2 s, sys: 228 ms, total: 2.23 s
Wall time: 1.96 s
CPU times: user 1.94 s, sys: 0 ns, total: 1.94 s
Wall time: 1.94 s
[[ 0.13538037+0.j        ]
 [ 0.15226989-0.22388607j]
 [-0.14070393-0.35612185j]
 [-0.43138871-0.09691328j]] 1.0

In [3]:

    

init = cutiepy.coherent(N,2)
H = cutiepy.num(N)
%time res = cutiepy.sesolve(H, init, np.linspace(0,1000,100))
%time res = cutiepy.sesolve(H, init, np.linspace(0,1000,100))
print(res[-1].numerical[:4], np.linalg.norm(res[-1].numerical))


    

    




CPU times: user 2.78 s, sys: 28 ms, total: 2.81 s
Wall time: 6.69 s
CPU times: user 1.43 s, sys: 0 ns, total: 1.43 s
Wall time: 1.43 s
[[ 0.13533528+0.j        ]
 [ 0.20064000+0.18167593j]
 [ 0.03788156+0.38090694j]
 [-0.26279469+0.35539523j]] 0.999999596058

In [4]:

    

%%time
H = qutip.num(N).data.A
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=1000)                                      
state0_num = qutip.coherent(N,2).data.A
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][:4], np.linalg.norm(res[-1]))


    

    




[[ 0.13533528+0.j        ]
 [ 0.15221946-0.22381195j]
 [-0.14065830-0.3560059j ]
 [-0.43125246-0.09688203j]] 0.999703107528
CPU times: user 5.74 s, sys: 6 s, total: 11.7 s
Wall time: 3.09 s