%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from qutip import *
from IPython.display import display
def solve(H, psi0, c_ops, a_ops, e_ops):
result_me = mesolve(H, psi0, times, c_ops, e_ops)
result_brme = brmesolve(H, psi0, times, a_ops, e_ops)
fig, ax = plot_expectation_values([result_me, result_brme])
display(fig)
plt.close(fig)
R, ekets = bloch_redfield_tensor(H, a_ops)
print("="* 20 + " Bloch-Redfield tensor: ")
display(R)
L = liouvillian(H, c_ops)
print("="* 20 + " Lindblad liouvilllian: ")
display(L)
print("="* 20 + " Bloch-Redfield steadystate dm")
R_rhoss_eb = steadystate(R)
R_rhoss = R_rhoss_eb.transform(ekets, True)
display(R_rhoss)
print("="* 20 + " Lindblad steadystate dm")
L_rhoss = steadystate(L)
display(L_rhoss)
print("="* 20 + " Steadystate expectation values")
print("R_ob: ", [expect(e, R_rhoss) for e in e_ops])
print("R_eb: ", [expect(e.transform(ekets), R_rhoss_eb) for e in e_ops])
print("L : ", [expect(e, L_rhoss) for e in e_ops])
print("="* 20 + " Dynamics final states")
print("R: ", [e[-1].real for e in result_brme.expect])
print("L: ", [e[-1] for e in result_me.expect])
delta = 0.0 * 2 * np.pi
epsilon = 0.5 * 2 * np.pi
gamma = 0.25
times = np.linspace(0, 50, 100)
H = delta/2 * sigmay() + epsilon/2 * sigmaz()
psi0 = (2 * basis(2, 0) + basis(2, 1)).unit()
c_ops = [np.sqrt(gamma) * sigmam()]
a_ops = [[sigmax(),lambda w : gamma * (w >= 0)]]
e_ops = [sigmax(), sigmay(), sigmaz()]
solve(H, psi0, c_ops, a_ops, e_ops)
==================== Bloch-Redfield tensor:
==================== Lindblad liouvilllian:
==================== Bloch-Redfield steadystate dm
==================== Lindblad steadystate dm
==================== Steadystate expectation values R_ob: [0.0, 0.0, -1.0] R_eb: [0.0, 0.0, -1.0] L : [0.0, 0.0, -1.0] ==================== Dynamics final states R: [0.0015427968283087687, 2.506068854864511e-07, -0.999994037354897] L: [0.001542796828343541, 2.5060686974779289e-07, -0.99999403735489634]
N = 10
w0 = 1.0 * 2 * np.pi
g = 0.05 * w0
kappa = 0.15
times = np.linspace(0, 50, 1000)
a = destroy(N)
H = w0 * a.dag() * a + g * (a + a.dag())
psi0 = ket2dm((basis(N, 4) + basis(N, 2) + basis(N,0)).unit())
a_ops = [[a + a.dag(),lambda w : kappa * (w >= 0)]]
e_ops = [a.dag() * a, a + a.dag()]
c_ops = [np.sqrt(kappa) * a]
solve(H, psi0, c_ops, a_ops, e_ops)
==================== Bloch-Redfield tensor:
==================== Lindblad liouvilllian:
==================== Bloch-Redfield steadystate dm
==================== Lindblad steadystate dm
==================== Steadystate expectation values R_ob: [0.0024999999999958857, -0.09999999999983651] R_eb: [0.0024999999999999467, -0.09999999999999895] L : [0.0024996438434595246, -0.09998575373839921] ==================== Dynamics final states R: [0.0034899627568894909, -0.097648226112852135] L: [0.0034896231351718215, -0.097634314487055257]
n_th = 1.5
c_ops = [np.sqrt(kappa * (n_th + 1)) * a, np.sqrt(kappa * n_th) * a.dag()]
w_th = w0/np.log(1 + 1/n_th)
def S_w_func(w):
if w >= 0:
return (n_th + 1) * kappa
else:
return (n_th + 1) * kappa * np.exp(w / w_th)
a_ops = [[a + a.dag(), S_w_func]]
solve(H, psi0, c_ops, a_ops, e_ops)
==================== Bloch-Redfield tensor:
==================== Lindblad liouvilllian:
==================== Bloch-Redfield steadystate dm
==================== Lindblad steadystate dm
==================== Steadystate expectation values R_ob: [1.441172521173384, -0.09593514655834637] R_eb: [1.4411725211733848, -0.0959351465583464] L : [1.4415837865925414, -0.09581420257265144] ==================== Dynamics final states R: [1.4412697964611352, -0.094946706436505962] L: [1.4416807381001973, -0.094807209086397559]
N = 10
a = tensor(destroy(N), identity(2))
sm = tensor(identity(N), destroy(2))
psi0 = ket2dm(tensor(basis(N, 1), basis(2, 0)))
a_ops = [[a + a.dag(),lambda w : kappa*(w > 0)]]
e_ops = [a.dag() * a, sm.dag() * sm]
w0 = 1.0 * 2 * np.pi
g = 0.05 * 2 * np.pi
kappa = 0.05
times = np.linspace(0, 150 * 2 * np.pi / g, 1000)
c_ops = [np.sqrt(kappa) * a]
H = w0 * a.dag() * a + w0 * sm.dag() * sm + g * (a + a.dag()) * (sm + sm.dag())
solve(H, psi0, c_ops, a_ops, e_ops)
==================== Bloch-Redfield tensor:
==================== Lindblad liouvilllian:
==================== Bloch-Redfield steadystate dm
==================== Lindblad steadystate dm
==================== Steadystate expectation values R_ob: [0.0006269556884685608, 0.0006253905020122077] R_eb: [0.0006269556893783309, 0.0006253905020103645] L : [0.0012531197519054338, 0.0012536697716992462] ==================== Dynamics final states R: [0.00062695568938809241, 0.00062539050201271512] L: [0.0012531197519057469, 0.0012536697716989808]
w0 = 1.0 * 2 * np.pi
g = 0.75 * 2 * np.pi
kappa = 0.05
times = np.linspace(0, 150 * 2 * np.pi / g, 1000)
c_ops = [np.sqrt(kappa) * a]
H = w0 * a.dag() * a + w0 * sm.dag() * sm + g * (a + a.dag()) * (sm + sm.dag())
solve(H, psi0, c_ops, a_ops, e_ops)
==================== Bloch-Redfield tensor:
==================== Lindblad liouvilllian:
==================== Bloch-Redfield steadystate dm
==================== Lindblad steadystate dm
==================== Steadystate expectation values R_ob: [0.267645307720865, 0.1560205514129042] R_eb: [0.267645307720865, 0.1560205514129042] L : [0.4465560981540205, 0.22225842348795147] ==================== Dynamics final states R: [0.26764721566971927, 0.15602106817515593] L: [0.44662603050034932, 0.22228388401493518]
from qutip.ipynbtools import version_table
version_table()
Software | Version |
---|---|
QuTiP | 4.2.0 |
Numpy | 1.13.1 |
SciPy | 0.19.1 |
matplotlib | 2.0.2 |
Cython | 0.25.2 |
Number of CPUs | 2 |
BLAS Info | INTEL MKL |
IPython | 6.1.0 |
Python | 3.6.1 |Anaconda custom (x86_64)| (default, May 11 2017, 13:04:09) [GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)] |
OS | posix [darwin] |
Wed Jul 19 22:06:57 2017 MDT |