本文整理汇总了Python中qutip.steadystate函数的典型用法代码示例。如果您正苦于以下问题:Python steadystate函数的具体用法?Python steadystate怎么用?Python steadystate使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了steadystate函数的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: _spectrum_es
def _spectrum_es(H, wlist, c_ops, a_op, b_op):
"""
Internal function for calculating the spectrum of the correlation function
:math:`\left<A(\\tau)B(0)\\right>`.
"""
if debug:
print(inspect.stack()[0][3])
# construct the Liouvillian
L = liouvillian(H, c_ops)
# find the steady state density matrix and a_op and b_op expecation values
rho0 = steadystate(L)
a_op_ss = expect(a_op, rho0)
b_op_ss = expect(b_op, rho0)
# eseries solution for (b * rho0)(t)
es = ode2es(L, b_op * rho0)
# correlation
corr_es = expect(a_op, es)
# covariance
cov_es = corr_es - np.real(np.conjugate(a_op_ss) * b_op_ss)
# spectrum
spectrum = esspec(cov_es, wlist)
return spectrum
示例2: test_ho_lgmres
def test_ho_lgmres():
"Steady state: Thermal HO - iterative-lgmres solver"
# thermal steadystate of an oscillator: compare numerics with analytical
# formula
a = destroy(40)
H = 0.5 * 2 * np.pi * a.dag() * a
gamma1 = 0.05
wth_vec = np.linspace(0.1, 3, 20)
p_ss = np.zeros(np.shape(wth_vec))
for idx, wth in enumerate(wth_vec):
n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature
c_op_list = []
rate = gamma1 * (1 + n_th)
c_op_list.append(np.sqrt(rate) * a)
rate = gamma1 * n_th
c_op_list.append(np.sqrt(rate) * a.dag())
rho_ss = steadystate(H, c_op_list, method='iterative-lgmres')
p_ss[idx] = np.real(expect(a.dag() * a, rho_ss))
p_ss_analytic = 1.0 / (np.exp(1.0 / wth_vec) - 1)
delta = sum(abs(p_ss_analytic - p_ss))
assert_equal(delta < 1e-3, True)
示例3: test_qubit_power
def test_qubit_power():
"Steady state: Thermal qubit - power solver"
# thermal steadystate of a qubit: compare numerics with analytical formula
sz = sigmaz()
sm = destroy(2)
H = 0.5 * 2 * np.pi * sz
gamma1 = 0.05
wth_vec = np.linspace(0.1, 3, 20)
p_ss = np.zeros(np.shape(wth_vec))
for idx, wth in enumerate(wth_vec):
n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature
c_op_list = []
rate = gamma1 * (1 + n_th)
c_op_list.append(np.sqrt(rate) * sm)
rate = gamma1 * n_th
c_op_list.append(np.sqrt(rate) * sm.dag())
rho_ss = steadystate(H, c_op_list, method='power')
p_ss[idx] = expect(sm.dag() * sm, rho_ss)
p_ss_analytic = np.exp(-1.0 / wth_vec) / (1 + np.exp(-1.0 / wth_vec))
delta = sum(abs(p_ss_analytic - p_ss))
assert_equal(delta < 1e-5, True)
示例4: steady
def steady(N = 20): # number of basis states to consider
n=num(N)
a = destroy(N)
H = a.dag() * a
print H.eigenstates()
#psi0 = basis(N, 10) # initial state
kappa = 0.1 # coupling to oscillator
c_op_list = []
n_th_a = 2 # temperature with average of 2 excitations
rate = kappa * (1 + n_th_a)
c_op_list.append(sqrt(rate) * a) # decay operators
rate = kappa * n_th_a
c_op_list.append(sqrt(rate) * a.dag()) # excitation operators
final_state = steadystate(H, c_op_list)
fexpt = expect(a.dag() * a, final_state)
#tlist = linspace(0, 100, 100)
#mcdata = mcsolve(H, psi0, tlist, c_op_list, [a.dag() * a], ntraj=100)
#medata = mesolve(H, psi0, tlist, c_op_list, [a.dag() * a])
#plot(tlist, mcdata.expect[0],
#plt.plot(tlist, medata.expect[0], lw=2)
plt.axhline(y=fexpt, color='r', lw=1.5) # ss expt. value as horiz line (= 2)
plt.ylim([0, 10])
plt.show()
示例5: find_expect
def find_expect(self, vg, pwr, fd):
if self.power_plot:
phi, pwr=vg
else:
phi, fd=vg
pwr_fridge=pwr-self.atten
lin_pwr=0.001*10**(pwr_fridge/10.0)
Omega=sqrt(lin_pwr/h*2.0)
gamma, Delta=self._get_GammaDelta(fd=fd, f0=self.f0, Np=self.Np, gamma=self.gamma)
g_el=self._get_Gamma_C(fd=fd)
wTvec=self._get_fTvec(phi=phi, gamma=gamma, Delta=Delta, fd=fd, Psaw=lin_pwr)
if self.acoustic_plot:
Om=Omega*sqrt(gamma/fd)
else:
Om=Omega*sqrt(g_el/fd)
wT = wTvec-fd*self.nvec #rotating frame of gate drive \omega_m-m*\omega_\gate
transmon_levels = Qobj(diag(wT[range(self.N_dim)]))
rate1 = (gamma+g_el)*(1.0 + self.N_gamma)
rate2 = (gamma+g_el)*self.N_gamma
c_ops=[sqrt(rate1)*self.a_op, sqrt(rate2)*self.a_dag]#, sqrt(rate3)*self.a_op, sqrt(rate4)*self.a_dag]
Omega_vec=-0.5j*(Om*self.a_dag - conj(Om)*self.a_op)
H=transmon_levels +Omega_vec
final_state = steadystate(H, c_ops) #solve master equation
fexpt=expect(self.a_op, final_state) #expectation value of relaxation operator
#return fexpt
if self.acoustic_plot:
return 1.0*gamma/Om*fexpt
else:
return 1.0*sqrt(g_el*gamma)/Om*fexpt
示例6: rhos
def rhos(self, nslice=None):
"""rho
return steadystate density matrices"""
self.precalc = True
if self.noisy:
print("N = {}, len() = {}".format(self.N_field_levels, len(self.long_range)))
def progress(*args):
if self.noisy:
print("|", sep="", end="", flush=True)
return args
if nslice is not None:
return np.asarray(
[progress(qt.steadystate(ham, self._c_ops())) for ham in list(self.hamiltonian())[nslice]]
).T
else:
return np.asarray([progress(qt.steadystate(ham, self._c_ops())) for ham in self.hamiltonian()]).T
示例7: test_dqd_current
def test_dqd_current():
"Counting statistics: current and current noise in a DQD model"
G = 0
L = 1
R = 2
sz = projection(3, L, L) - projection(3, R, R)
sx = projection(3, L, R) + projection(3, R, L)
sR = projection(3, G, R)
sL = projection(3, G, L)
w0 = 1
tc = 0.6 * w0
GammaR = 0.0075 * w0
GammaL = 0.0075 * w0
nth = 0.00
eps_vec = np.linspace(-1.5*w0, 1.5*w0, 20)
J_ops = [GammaR * sprepost(sR, sR.dag())]
c_ops = [np.sqrt(GammaR * (1 + nth)) * sR,
np.sqrt(GammaR * (nth)) * sR.dag(),
np.sqrt(GammaL * (nth)) * sL,
np.sqrt(GammaL * (1 + nth)) * sL.dag(),
]
I = np.zeros(len(eps_vec))
S = np.zeros(len(eps_vec))
for n, eps in enumerate(eps_vec):
H = (eps/2 * sz + tc * sx)
L = liouvillian(H, c_ops)
rhoss = steadystate(L)
I[n], S[n] = countstat_current_noise(L, [], rhoss=rhoss, J_ops=J_ops)
I2 = countstat_current(L, rhoss=rhoss, J_ops=J_ops)
assert_(abs(I[n] - I2) < 1e-8)
I2 = countstat_current(L, c_ops, J_ops=J_ops)
assert_(abs(I[n] - I2) < 1e-8)
Iref = tc**2 * GammaR / (tc**2 * (2 + GammaR/GammaL) +
GammaR**2/4 + eps_vec**2)
Sref = 1 * Iref * (
1 - 8 * GammaL * tc**2 *
(4 * eps_vec**2 * (GammaR - GammaL) +
GammaR * (3 * GammaL * GammaR + GammaR**2 + 8*tc**2)) /
(4 * tc**2 * (2 * GammaL + GammaR) + GammaL * GammaR**2 +
4 * eps_vec**2 * GammaL)**2
)
assert_allclose(I, Iref, 1e-4)
assert_allclose(S, Sref, 1e-4)
示例8: find_expect
def find_expect(vg): #phi=0.1, Omega_vec=3.0):
phi, Omega=vg#.shape
Omega_vec=- 0.5j*(Omega*adag - conj(Omega)*a)
Ej = Ejmax*absolute(cos(pi*phi)) #Josephson energy as function of Phi.
wTvec = -Ej + sqrt(8.0*Ej*Ec)*(nvec+0.5)+Ecvec #\omega_m
wT = wTvec-wdvec #rotating frame of gate drive \omega_m-m*\omega_\gate
transmon_levels = Qobj(diag(wT[range(N)]))
H=transmon_levels +Omega_vec #- 0.5j*(Omega_true*adag - conj(Omega_true)*a)
final_state = steadystate(H, c_op_list) #solve master equation
return expect( a, final_state) #expectation value of relaxation operator
示例9: find_expect
def find_expect(vg, self=a): #phi=0.1, Omega_vec=3.0):
phi, Omega=vg#.shape
Omega_vec=-0.5j*(Omega*self.a_dag - conj(Omega)*self.a_op)
Ej = self.Ejmax*absolute(cos(pi*phi)) #Josephson energy as function of Phi.
wTvec = (-Ej + sqrt(8.0*Ej*self.Ec)*(self.nvec+0.5)+self.Ecvec)/h #\omega_m
wT = wTvec-a.fdvec #rotating frame of gate drive \omega_m-m*\omega_\gate
transmon_levels = Qobj(diag(wT[range(self.N_dim)]))
H=transmon_levels +Omega_vec #- 0.5j*(Omega_true*adag - conj(Omega_true)*a)
final_state = steadystate(H, self.c_ops) #solve master equation
return expect( self.a_op, final_state) #expectation value of relaxation operator
示例10: test_driven_cavity_power_bicgstab
def test_driven_cavity_power_bicgstab():
"Steady state: Driven cavity - power-bicgstab solver"
N = 30
Omega = 0.01 * 2 * np.pi
Gamma = 0.05
a = destroy(N)
H = Omega * (a.dag() + a)
c_ops = [np.sqrt(Gamma) * a]
M = build_preconditioner(H, c_ops, method='power')
rho_ss = steadystate(H, c_ops, method='power-bicgstab', M=M, use_precond=1)
rho_ss_analytic = coherent_dm(N, -1.0j * (Omega)/(Gamma/2))
assert_((rho_ss - rho_ss_analytic).norm() < 1e-4)
示例11: test_driven_cavity_lgmres
def test_driven_cavity_lgmres():
"Steady state: Driven cavity - iterative-lgmres solver"
N = 30
Omega = 0.01 * 2 * np.pi
Gamma = 0.05
a = destroy(N)
H = Omega * (a.dag() + a)
c_ops = [np.sqrt(Gamma) * a]
rho_ss = steadystate(H, c_ops, method='iterative-lgmres')
rho_ss_analytic = coherent_dm(N, -1.0j * (Omega)/(Gamma/2))
assert_((rho_ss - rho_ss_analytic).norm() < 1e-4)
示例12: jc_steadystate
def jc_steadystate(self, N, wc, wa, g, kappa, gamma,
pump, psi0, use_rwa, tlist):
# Hamiltonian
a = tensor(destroy(N), identity(2))
sm = tensor(identity(N), destroy(2))
if use_rwa:
# use the rotating wave approxiation
H = wc * a.dag(
) * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
else:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (
a.dag() + a) * (sm + sm.dag())
# collapse operators
c_op_list = []
n_th_a = 0.0 # zero temperature
rate = kappa * (1 + n_th_a)
c_op_list.append(np.sqrt(rate) * a)
rate = kappa * n_th_a
if rate > 0.0:
c_op_list.append(np.sqrt(rate) * a.dag())
rate = gamma
if rate > 0.0:
c_op_list.append(np.sqrt(rate) * sm)
rate = pump
if rate > 0.0:
c_op_list.append(np.sqrt(rate) * sm.dag())
# find the steady state
rho_ss = steadystate(H, c_op_list)
return expect(a.dag() * a, rho_ss), expect(sm.dag() * sm, rho_ss)
示例13: _spectrum_pi
def _spectrum_pi(H, wlist, c_ops, a_op, b_op, use_pinv=False):
"""
Internal function for calculating the spectrum of the correlation function
:math:`\left<A(\\tau)B(0)\\right>`.
"""
#print issuper(H)
L = H if issuper(H) else liouvillian(H, c_ops)
tr_mat = tensor([qeye(n) for n in L.dims[0][0]])
N = prod(L.dims[0][0])
A = L.full()
b = spre(b_op).full()
a = spre(a_op).full()
tr_vec = transpose(mat2vec(tr_mat.full()))
rho_ss = steadystate(L)
rho = transpose(mat2vec(rho_ss.full()))
I = identity(N * N)
P = kron(transpose(rho), tr_vec)
Q = I - P
spectrum = zeros(len(wlist))
for idx, w in enumerate(wlist):
if use_pinv:
MMR = pinv(-1.0j * w * I + A)
else:
MMR = dot(Q, solve(-1.0j * w * I + A, Q))
s = dot(tr_vec,
dot(a, dot(MMR, dot(b, transpose(rho)))))
spectrum[idx] = -2 * real(s[0, 0])
return spectrum
示例14: find_expect
def find_expect(phi=0.1, Omega_vec=3.0, wd=wd):
wdvec=nvec*wd
Ej = Ejmax*absolute(cos(pi*phi)) #; % Josephson energy as function of Phi.
wq0=sqrt(8*Ej*Ec)
X=Np*pi*(wq0-f0)/f0
Ba=-1*gamma*(sin(2.0*X)-2.0*X)/(2.0*X**2.0)
Ecp=Ec*(1-2*Ba/wq0)
Ecvec=-Ecp*(6.0*nvec**2+6.0*nvec+3.0)/12.0
wTvec = -Ej + sqrt(8.0*Ej*Ecp)*(nvec+0.5)+Ecvec
#print diff(wTvec)
if phi==0.4:
print Ba/wq0
if 0:
zr=zeros(N)
X=Np*pi*(diff(wTvec)-f0)/f0
Ba=-0*gamma*(sin(2.0*X)-2.0*X)/(2.0*X**2.0)
ls=-Ba#/(2*Ct)/(2*pi)#/1e6
zr[1:]=ls
#print Ga0/(2*Ct)/(2*pi)
#print zr
wTvec=wTvec+zr
wTvec = wTvec-wdvec
#print wTvec+diff(wTvec)
transmon_levels = Qobj(diag(wTvec[range(N)]))
H=transmon_levels +Omega_vec #- 0.5j*(Omega_true*adag - conj(Omega_true)*a)
final_state = steadystate(H, c_op_list)
return expect( tm_l, final_state)
示例15: H
wVals = np.linspace(-2 * np.pi, 6 * np.pi, 100)
# Solving the dynamics of a spin 1/2 in B0 and Brf in the lab frame
def H(w):
""" Hamiltonian in the lab frame. """
return w0 / 2 * qt.sigmaz() + w1 * np.cos(w) * qt.sigmax()
c1 = np.sqrt(1.0/T1) * qt.destroy(2)
c2 = np.sqrt(1.0/T2) / 2 * qt.destroy(2)
c_op_list = [c1,c2]
wHVals = np.linspace(-2 * np.pi, 6 * np.pi, 51)
MxH = [qt.expect(qt.sigmax(), qt.steadystate(H(w), c_op_list, maxiter=100)) for w in wHVals]
MyH = [qt.expect(qt.sigmay(), qt.steadystate(H(w), c_op_list, maxiter=100)) for w in wHVals]
MzH = [qt.expect(qt.sigmaz(), qt.steadystate(H(w), c_op_list, maxiter=100)) for w in wHVals]
# Plot results of the algebraic formula and the solution given by qutip
plt.ylim([-1.0, 1.0])
plt.plot(wVals, Mx(wVals), label='Mx')
plt.plot(wVals, My(wVals), label='My')
plt.plot(wVals, Mz(wVals), label='Mz')
plt.plot(wHVals, MxH, 'bo', label='MxH')
plt.plot(wHVals, MyH, 'go', label='MyH')
plt.plot(wHVals, MzH, 'ro', label='MzH')