Reference
Question
In the book Breuer, Heinz-Peter, and Francesco Petruccione. The Theory of Open Quantum Systems. Oxford: Clarendon Press, 2009., there is equation(10.17)
\[ \frac{\mathrm{d}}{\mathrm{d}t}c_1(t) = - \int_0^t \mathrm{d}t_1 f(t-t_1) c_1(t_1) \]where
\[ f(t) = \frac{1}{2} \gamma_0\lambda e^{-\lambda |t|} \]\(\gamma_0\) and \(\lambda\) are parameters, for example, we can set \(\gamma_0=1, \lambda=10\). For this \(f(t)\), we can solve it analytically, the solution is
\[ c_1(t) = c_1(0) e^{-\lambda t/2} \left[\cosh\left(\frac{dt}{2}\right) + \frac{\lambda}{d} \sinh\left(\frac{dt}{2}\right) \right] \]where \(d = \sqrt{\lambda^2 - 2\gamma_0\lambda}\).
My question is, how to numerically solve (for example, using python function scipy.integrate.solve_ivp, which is a function solving OEDs by Runge-Kutta method or other methods) this equation?
I want to numerically solve it, because
In [1]
import numpy as np
from scipy.integrate import quad, solve_ivp
import qutip as qutip
import matplotlib.pyplot as plt
In [2]
def ft(t, g0=1, lam=10):
"""Breuer Eq.(10.44)"""
return 1/2 * g0 * lam * np.exp(-lam*np.abs(t))
def dc_dt(t, c0=1, lam=10, g0=1):
"""Right hand side of Breuer Eq.(10.17),
chose f(t-t1) as Breuer Eq.(10.44)"""
res = quad(lambda t1: ft(t-t1)*ct(t1), 0, t)[0]
return -res
def ct(t, c0=1, lam=10, g0=1):
"""Breuer Eq.(10.45)"""
d = np.emath.sqrt(lam**2 - 2*g0*lam)
ct = c0 * np.exp(-lam*t/2)
ct *= np.cosh(d*t/2) + lam/d*np.sinh(d*t/2)
return ct
def ct0(t, c0=1, g0=1):
"""Markovian limit solution"""
return c0*np.exp(-g0*t/2)
In [3]
t = np.linspace(0, 10, 1000)
plt.plot(t, [np.abs(ct(ti, lam=1)) for ti in t], 'y')
plt.plot(t, np.abs(ct(t, lam=2.1)), 'r')
plt.plot(t, np.abs(ct(t, lam=10)), 'g')
plt.plot(t, np.abs(ct(t, lam=100)), 'b')
plt.plot(t[::10], ct0(t)[::10], 'kx', label='Markov')
plt.xlabel(r'$t$')
plt.legend()
<matplotlib.legend.Legend at 0x7bed65eda740>
<Figure size 640x480 with 1 Axes>
In [4]
# check Breuer Eq.(10.45)
dt = 1e-2
t = np.arange(0, 4, dt)
def ft(t, g0=1, lam=10):
"""Breuer Eq.(10.44)"""
return 1/2 * g0 * lam * np.exp(-lam*np.abs(t))
def dc_dt(t, c0=1, lam=10, g0=1):
"""Right hand side of Breuer Eq.(10.17),
chose f(t-t1) as Breuer Eq.(10.44)"""
res = quad(lambda t1: ft(t-t1)*ct(t1), 0, t)[0]
return -res
cts = ct(t)
plt.plot(t[1:], np.diff(cts)/dt)
plt.plot(t, np.array([dc_dt(ti) for ti in t]), 'x')
[<matplotlib.lines.Line2D at 0x7bed63534760>]
<Figure size 640x480 with 1 Axes>
In [5]
dt = 1e-2
ts = np.arange(0, 5, dt)
dcdt = [dc_dt(ti) for ti in ts]
plt.plot(ts, dcdt, 'r--')
plt.plot(ts[1:], np.diff(dcdt)/dt, 'g')
def myd(t):
return - ft(0)*ct(t) + ft(t)*ct(0)
plt.plot(ts, myd(ts), 'k')
[<matplotlib.lines.Line2D at 0x7bed63537e20>]
<Figure size 640x480 with 1 Axes>
In [2]
print('hello')