In [1]:
import cutiepy
cutiepy.interactive.INTERACTIVE = False
import qutip
import numpy as np
import scipy
ts = np.linspace(0,10000*np.pi,10000)
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())
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))
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]))