I tried to write an example to explain my question, but another error shows up. Here is the code and error
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
wc = 1.0 * 2 * np.pi
wa = 1.0 * 2 * np.pi
g = 0.05 * 2 * np.pi
kappa = 0.005
gamma = 0.05
N = 15
step = 1001
vac = tensor(basis(N,0),basis(2,0))
a = tensor(destroy(N),qeye(2))
sm = tensor(qeye(N),destroy(2))
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
c_ops = [np.sqrt(kappa)*a,np.sqrt(gamma)*sm]
t_list = np.linspace(0,25.0,step)
psi0 = sm.dag() * vac
result1 = mesolve(H,psi0,t_list,c_ops,[a.dag()*a,sm.dag()*sm])
result0 = mcsolve(H,psi0,t_list,[],[a.dag()*a,sm.dag()*sm])
n_c_0 = result0.expect[0]
n_a_0 = result0.expect[1]
n_c_1 = result1.expect[0]
n_a_1 = result1.expect[1]
fig,axes = plt.subplots(1,1,figsize = (10,6))
axes.plot(t_list,n_c_0,label = 'Cavity(without dissipasion)')
axes.plot(t_list,n_a_0,label = 'Atom(without dissipasion)')
axes.plot(t_list,n_c_1,'-.',label = 'Cavity(with dissipasion)')
axes.plot(t_list,n_a_1,'-.',label = 'Atom(with dissipasion)')
axes.legend(loc = 'best')
axes.set_xlabel('Time')
axes.set_ylabel('Occupation probability')
axes.set_title('Vacuum Rabi oscillations')
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
This part works well. But when I tried to solve the time-evolution operator:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
U1 = propagator(H, t_list,c_ops)
U0 = propagator(H, t_list, [])
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-4-8134c024ad2f> in <module>
----> 1 U1 = propagator(H, t_list,c_ops)
2 U0 = propagator(H, t_list, [])
C:\anaconda3\lib\site-packages\qutip\propagator.py in propagator(H, t, c_op_list, args, options, unitary_mode, parallel, progress_bar, _safe_mode, **kwargs)
239 rho0 = Qobj(sp.csr_matrix(([1], ([row_idx], [col_idx])),
240 shape=(N, N), dtype=complex))
--> 241 output = mesolve(H, rho0, tlist, c_op_list, [], args, options,
242 _safe_mode=False)
243 for k, t in enumerate(tlist):
C:\anaconda3\lib\site-packages\qutip\mesolve.py in mesolve(H, rho0, tlist, c_ops, e_ops, args, options, progress_bar, _safe_mode)
263 raise Exception("Invalid H type")
264
--> 265 func, ode_args = ss.makefunc(ss, rho0, args, e_ops, options)
266
267 if _safe_mode:
C:\anaconda3\lib\site-packages\qutip\mesolve.py in _qobjevo_set(HS, rho0, args, e_ops, opt)
349 # Should be caught earlier in mesolve.
350 raise ValueError("rho0 must be a ket, density matrix or superoperator")
--> 351 _test_liouvillian_dimensions(H_td.cte.dims, rho0.dims)
352 return func, ()
353
C:\anaconda3\lib\site-packages\qutip\mesolve.py in _test_liouvillian_dimensions(L_dims, rho_dims)
330 raise ValueError("Liouvillian had nonsquare dims: " + str(L_dims))
331 if not ((L_dims[1] == rho_dims) or (L_dims[1] == rho_dims[0])):
--> 332 raise ValueError("".join([
333 "incompatible Liouvillian and state dimensions: ",
334 str(L_dims), " and ", str(rho_dims),
ValueError: incompatible Liouvillian and state dimensions: [[[15, 2], [15, 2]], [[15, 2], [15, 2]]] and [[30], [30]]