求解稳态解

介绍

对于具有衰减率大于相应激发率的与时间无关的开放量子系统，系统将趋向于一个稳态，当\(t\rightarrow\infty\)时，满足方程

\[\frac{d\hat{\rho}_{ss}}{dt}=\mathcal{L}\hat{\rho}_{ss}=0.\]

尽管对时间独立性的要求似乎相当严格，但人们通常可以采用相互作用绘景的变换来产生一个时间独立的哈密顿量。对于许多这些系统，求解渐近密度矩阵 \(\hat{\rho}_{ss}\) 可以通过直接或迭代解法比使用主方程或蒙特卡罗模拟更快地实现。尽管稳态方程具有简单的数学形式，但刘维尔算子的性质使得这个方程的解远非直截了当。

QuTiP中的稳态求解器

在QuTiP中，系统哈密顿量或Liouvillian的稳态解由steadystate给出。此函数实现了多种不同的方法来寻找稳态，每种方法都有其优缺点，可以通过method关键字参数选择使用的方法。

方法	关键词	描述
直接（默认）	‘direct’	直接求解 \(Ax=b\)。
特征值	‘eigen’	迭代地找到 \(\mathcal{L}\) 的零特征值。
逆幂法	‘power’	使用逆幂法求解。
SVD	‘svd’	通过Liouvillian的密集 SVD获得稳态解。

函数 steadystate 可以接受一个哈密顿量和一组崩溃算子作为输入，内部生成相应的Lindblad形式的Liouvillian超算子，或者用户传递的Liouvillian。

无论是 "direct" 还是 "power" 方法都需要解决一个线性方程组。为此，有多个可用的求解器：``

求解器	原始函数	描述
“solve”	`numpy.linalg.solve`	来自numpy的密集求解器。
“lstsq”	`numpy.linalg.lstsq`	密集最小二乘求解器。
“spsolve”	`scipy.sparse.linalg.spsolve`	来自scipy的稀疏求解器。
“gmres”	`scipy.sparse.linalg.gmres`	广义最小残差迭代求解器。
“lgmres”	`scipy.sparse.linalg.lgmres`	LGMRES 迭代求解器。
“bicgstab”	`scipy.sparse.linalg.bicgstab`	双共轭梯度稳定迭代求解器。
“mkl_spsolve”	`pardiso`	来自MKL的Intel Pardiso LU求解器

QuTiP 可以利用 Anacoda (2.5+) 和 Intel Python 发行版中附带的 Intel Math Kernel 库中的 Intel Pardiso LU 求解器。与 SciPy 使用的标准 SuperLU 方法相比，这显著提高了性能。要验证 QuTiP 是否可以找到必要的库，可以在 QuTiP 的关于框中检查 INTEL MKL Ext: True (about)。

使用稳态求解器

求解一般系统的Lindblad主方程的稳态解可以使用steadystate来完成：

>>> rho_ss = steadystate(H, c_ops)

其中 H 是表示系统哈密顿量的量子对象，c_ops 是系统崩溃算子的量子对象列表。输出标记为 rho_ss，是系统的稳态解。如果没有其他关键字传递给求解器，默认使用 'direct' 方法与 numpy.linalg.solve，生成一个精确到机器精度的解，但需要大量内存。然而，Liouvillians 通常非常稀疏，使用稀疏求解器可能更受欢迎：

rho_ss = steadystate(H, c_ops, method="power", solver="spsolve")

其中 method='power' 表示我们正在使用逆幂解法，而 solver="spsolve" 表示使用稀疏求解器。

稀疏求解器在分解矩阵时可能仍然会使用相当多的内存，因为Liouvillian通常具有较大的带宽。为了解决这个问题，steadystate 允许使用在 Additional Solver Arguments 中列出的带宽最小化算法。例如：

rho_ss = steadystate(H, c_ops, solver="spsolve", use_rcm=True)

其中 use_rcm=True 开启了一个带宽最小化程序。

虽然不明显，但'direct'、'eigen'和'power'方法在内部都使用了LU分解，因此可能会占用大量内存。相比之下，迭代求解器如'gmres'、'lgmres'和'bicgstab'不会对矩阵进行分解，因此比LU方法占用更少的内存，并且原则上允许处理非常大的系统规模。缺点是这些方法可能比直接方法花费更长的时间，因为Liouvillian矩阵的条件数较大，表明这些迭代方法需要大量的迭代才能收敛。为了克服这个问题，可以使用预条件子\(M\)来求解（修改后的）Liouvillian的近似逆，从而更好地调节问题，加快收敛速度。使用预条件子实际上可以使这些迭代方法比其他求解方法更快。预条件的问题在于它仅对Hermitian矩阵有明确的定义。由于Liouvillian是非Hermitian的，因此不能保证找到一个好的预条件子。此外，即使找到了预条件子，也不能保证它具有好的条件数。QuTiP在使用迭代求解器'gmres'、'lgmres'和'bicgstab'时，可以通过设置use_precond=True来使用不完全LU预条件子。预条件子可以选择性地使用对称和反对称矩阵排列的组合，以尝试改进预条件过程。这些特性在附加求解器参数部分中讨论。即使有了这些最先进的排列，由于在这一领域缺乏数学研究，为非对称矩阵生成成功的预条件子目前仍然是一个试错过程。建议在选择其他方法之前，首先使用没有附加参数的直接求解器。

寻找稳态解不仅限于主方程的Lindblad形式。任何由哈密顿量和坍缩算子构建的时间无关的Liouvillian都可以用作输入：

>>> rho_ss = steadystate(L)

其中 L 是 Louvillian。在这种情况下，所有额外的参数也可以使用。

额外的求解器参数

以下额外的求解器参数可用于稳态求解器：

关键词	默认值	描述
weight	None	设置在`'direct'`方法中使用的权重因子。
use_precond	False	在使用`'gmres'`和`'lgmres'`方法时生成预处理器。
use_rcm	False	使用反向Cuthill-Mckee重新排序，以最小化在LU分解中使用的修改后的Liouvillian的带宽。
use_wbm	False	使用加权二分匹配算法尝试使修改后的Liouvillian更加对角占优，从而更有利于预处理。
power_tol	1e-12	使用‘power’方法时的解决方案容差。
power_maxiter	10	幂方法的最大迭代次数。
power_eps	1e-15	在“power”方法中使用的小权重。
**kwargs	{}	传递给线性代数求解器的选项。有关完整列表，请参阅scipy的相应文档。

更多信息可以在 steadystate 文档字符串中找到。

示例：热浴中的谐振子

一个达到稳态的系统的简单例子是耦合到热环境的谐振子。下面我们考虑一个谐振子，最初处于\(\left|10\right>\)数态，并且弱耦合到一个以平均粒子期望值\(\left=2\)为特征的热环境。我们通过主方程和蒙特卡罗方法计算演化，并看到它们收敛到稳态解。这里我们选择只执行几条蒙特卡罗轨迹，以便我们可以将这种演化与主方程解区分开来。

import numpy as np
import matplotlib.pyplot as plt

import qutip

# Define paramters
N = 20  # number of basis states to consider
a = qutip.destroy(N)
H = a.dag() * a
psi0 = qutip.basis(N, 10)  # initial state
kappa = 0.1  # coupling to oscillator

# collapse operators
c_op_list = []
n_th_a = 2  # temperature with average of 2 excitations
rate = kappa * (1 + n_th_a)
if rate > 0.0:
    c_op_list.append(np.sqrt(rate) * a)  # decay operators
rate = kappa * n_th_a
if rate > 0.0:
    c_op_list.append(np.sqrt(rate) * a.dag())  # excitation operators

# find steady-state solution
final_state = qutip.steadystate(H, c_op_list)
# find expectation value for particle number in steady state
fexpt = qutip.expect(a.dag() * a, final_state)

tlist = np.linspace(0, 50, 100)
# monte-carlo
mcdata = qutip.mcsolve(H, psi0, tlist, c_op_list, e_ops=[a.dag() * a], ntraj=100)
# master eq.
medata = qutip.mesolve(H, psi0, tlist, c_op_list, e_ops=[a.dag() * a])

plt.plot(tlist, mcdata.expect[0], tlist, medata.expect[0], lw=2)
# plot steady-state expt. value as horizontal line (should be = 2)
plt.axhline(y=fexpt, color='r', lw=1.5)
plt.ylim([0, 10])
plt.xlabel('Time', fontsize=14)
plt.ylabel('Number of excitations', fontsize=14)
plt.legend(('Monte-Carlo', 'Master Equation', 'Steady State'))
plt.title(
    r'Decay of Fock state $\left|10\rangle\right.$'
    r' in a thermal environment with $\langle n\rangle=2$'
)
plt.show()