Python QuTiP 中的集成未成功

Question

我一直在尝试使用 QuTiP 来求解量子力学矩阵微分方程（Lindblad 方程）。这是代码：

from qutip import *
from matplotlib import *
import numpy as np

hamiltonian = np.array([[215, -104.1, 5.1, -4.3  ,4.7,-15.1 ,-7.8 ],
[-104.1,  220.0, 32.6 ,7.1, 5.4, 8.3, 0.8],
      [ 5.1, 32.6, 0.0, -46.8, 1.0 , -8.1, 5.1],
     [-4.3, 7.1, -46.8, 125.0, -70.7, -14.7, -61.5],
       [ 4.7, 5.4, 1.0, -70.7, 450.0, 89.7, -2.5],
    [-15.1, 8.3, -8.1, -14.7, 89.7, 330.0, 32.7],
     [-7.8, 0.8, 5.1, -61.5,  -2.5, 32.7,  280.0]])
H=Qobj(hamiltonian)
ground=Qobj(np.array([[ 0.0863685 ],
  [ 0.17141713],
  [-0.91780802],
  [-0.33999268],
  [-0.04835763],
  [-0.01859027],
  [-0.05006013]]))

rho0 = ground*ground.dag()
from scipy.constants import *
ktuple=physical_constants['Boltzmann constant in eV/K']
k = ktuple[0]* 8065.6
htuple = physical_constants['Planck constant in eV s'] 
hbar = (htuple[0]* 8065.6)/(2*pi)
gamma=(2*pi)*((k*300)/hbar)*(35/(150*hbar))

L1 = Qobj(np.array([[1,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L2 = Qobj(np.array([[0,0,0,0,0,0,0],[0,1,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L3 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,1,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L4 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,1,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L5 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,1,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L6 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,1,0],[0,0,0,0,0,0,0]]))
L7 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,1]]))
#Since our gamma variable cannot be directly applied onto the Lindblad operator, we must multiply it with the collapse operators:

L1=gamma*L1
L2=gamma*L2
L3=gamma*L3
L4=gamma*L4
L5=gamma*L5
L6=gamma*L6
L7=gamma*L7

options = Options(nsteps=100000)

results = mesolve(H, rho0, np.linspace(0.0, 1000, 20),[L1,L2,L3,L4,L5,L6,L7],[],options=options)
print results

此代码应该求解以下等式：

其中L_i是矩阵（列表中：[L1,L2,L3,L4,L5,L6,L7]），H是哈密顿矩阵，另一个矩阵，是密度矩阵，是等于的常数，其中 T 是温度，k 是玻尔兹曼常数，，其中 h 是普朗克常数。

每次我运行代码，它都会给我以下错误：

ZVODE--  At T (=R1) and step size H (=R2), the
       corrector convergence failed repeatedly
       or with abs(H) = HMIN
      In above,  R1 =  0.0000000000000D+00   R2 =  0.1202322246215D-36
/usr/local/lib/python2.7/dist-packages/scipy/integrate/_ode.py:853: UserWarning: zvode: Repeated convergence failures. (Perhaps bad Jacobian supplied or wrong choice of MF o
r tolerances.)
  'Unexpected istate=%s' % istate))
Traceback (most recent call last):
  File "lindbladqutip.py", line 48, in <module>
    results = mesolve(H, rho0, np.linspace(0.0, 1000, 20),[L1,L2,L3,L4,L5,L6,L7],[],options=options)
  File "/projects/d6138712-e5f4-4d85-9d4d-77ce0a7b4a61/.local/lib/python2.7/site-packages/qutip/mesolve.py", line 264, in mesolve
    progress_bar)
  File "/projects/d6138712-e5f4-4d85-9d4d-77ce0a7b4a61/.local/lib/python2.7/site-packages/qutip/mesolve.py", line 692, in _mesolve_const
    return _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar)
  File "/projects/d6138712-e5f4-4d85-9d4d-77ce0a7b4a61/.local/lib/python2.7/site-packages/qutip/mesolve.py", line 866, in _generic_ode_solve
    raise Exception("ODE integration error: Try to increase "
Exception: ODE integration error: Try to increase the allowed number of substeps by increasing the nsteps parameter in the Options class.

调试分析了一下，好像第一次或第二次集成失败了。该错误告诉我增加我试过的 nsteps 参数。即使那样它也失败了。更改时间列表（np.linspace 函数生成时间列表）也没有效果。

我非常想知道我能做些什么来修复这个错误。如果大家需要更多详细信息，请发表评论。

感谢您的帮助！

Answer 1

从编程的角度来看，问题似乎出在您的 gamma 值以及折叠运算符的大小上。打印出 gamma - 它的顺序是 10**25 - 这似乎是阻止求解器收敛的原因。

只是为了测试（我是工程师，不是量子物理学家...），我输入了较小的 gamma 值（例如 0.1），求解器似乎有效并且给出了明显合理的值在 results.states

中输出

我不太明白你的gamma - 好像的单位是cm^-1s^-2 就像你设置的那样。我想知道您是否只想除以 hbar 一次，也许吧。正如我所说，我不是量子物理学家，所以我只是根据使编程挂在一起的原因以及一些维度分析来猜测。

编辑

OP 在评论中指出 gamma 的错误数量级/单位似乎是编程问题（即防止数值微积分收敛），但并不完全清楚如何计算 gamma .在此阶段，可能值得在 http://physics.stackexchange.com or http://math.stackexchange.com 上发布一个关于此的问题 - 如有必要，请参考此问题作为上下文。

编辑 2

我注意到你问了这个 related question on the physics site. This makes it clear where the expression for gamma comes from，从而澄清了在这个问题中简单地表示为 30 和 150 的常数项实际上有单位（分别是能量和频率）。这改变了量纲分析——gamma 的单位是 s^-1 或者，通过适当的转换，cm^-1.

它还显示了您在评论中提到的值 - 300 cm^-1.

Python QuTiP 中的集成未成功

Integration not successful in Python QuTiP

python

numpy

scipy

integrate

qutip