在Python中用复数矩阵作为初始值解ODE

Question

我有一个冯·诺依曼方程，形式是：dr/dt = - i [H, r]，其中 r 和 H 是由复数组成的方阵，我需要用 Python 脚本找到 r(t)。

有没有什么标准工具可以用来求解这样的方程？

之前我在解决另一个方程时，那个方程的初始值是一个向量，比如薛定谔方程：dy/dt = - i H y，我使用了 scipy.integrate.ode 函数（'zvode'），但是尝试用同样的函数来处理冯·诺依曼方程时，出现了以下错误：

**scipy/integrate/_ode.py:869: UserWarning: zvode: Illegal input detected. (See printed message.)
ZVODE--  ZWORK length needed, LENZW (=I1), exceeds LZW (=I2)
self.messages.get(istate, 'Unexpected istate=%s' % istate))
  In above message,  I1 =        72   I2 =        24**

这是我的代码：

def integrate(r, t0, t1, dt):
  e = linspace(t0, t1, (t1 - t0) / dt + 10)
  g = linspace(t0, t1, (t1 - t0) / dt + 10)
  u = linspace(t0, t1, (t1 - t0) / dt + 10)
  while r.successful() and r.t < t1:
    r.integrate(r.t + dt)
    e[r.t / dt] = abs(r.y[0][0]) ** 2
    g[r.t / dt] = abs(r.y[1][1]) ** 2
    u[r.t / dt] = abs(r.y[2][2]) ** 2
  return e, g, u


# von Neumann equation's
def right_part(t, rho):
  hamiltonian = (h / 2) * array(
    [[delta, omega_s, omega_p / 2.0 * sin(t * w_p)],
    [omega_s, 0.0, 0.0],
    [omega_p / 2.0 * sin(t * w_p), 0.0, 0.0]],
    dtype=complex128)
  return (dot(hamiltonian, rho) - dot(rho, hamiltonian)) / (1j * h)


def create_integrator():
  r = ode(right_part).set_integrator('zvode', method='bdf', with_jacobian=False)
  psi_init = array([[1.0, 0.0, 0.0],
                   [0.0, 0.0, 0.0],
                   [0.0, 0.0, 0.0]], dtype=complex128)
  t0 = 0
  r.set_initial_value(psi_init, t0)
  return r, t0


def main():
  r, t0 = create_integrator()
  t1 = 10 ** -6
  dt = 10 ** -11
  e, g, u = integrate(r, t0, t1, dt)

main()

Answer 1

QuTiP 提供了一些很不错的工具，可以用来处理这类问题，主要是通过使用像主方程和林布拉德阻尼项这样的概念。QuTiP 本身其实只是对 scipy.odeint 的一个简单封装，但它让很多操作变得更加简单易懂，尤其是因为它有很好的文档说明。

Answer 2

我创建了一个叫做 odeintw 的工具，它是 scipy.integrate.odeint 的一个封装，可以处理复杂的矩阵方程，比如这个。你可以查看 如何在PYTHON中绘制解矩阵耦合微分方程时的特征值？，那也是一个关于矩阵微分方程的问题。

这里有一个简化版的代码，展示了你可以如何使用它。（为了简单起见，我去掉了你例子中的大部分常数）。

import numpy as np
from odeintw import odeintw


def right_part(rho, t, w_p):
    hamiltonian = (1. / 2) * np.array(
        [[0.1, 0.01, 1.0 / 2.0 * np.sin(t * w_p)],
        [0.01, 0.0, 0.0],
        [1.0 / 2.0 * np.sin(t * w_p), 0.0, 0.0]],
        dtype=np.complex128)
    return (np.dot(hamiltonian, rho) - np.dot(rho, hamiltonian)) / (1j)


psi_init = np.array([[1.0, 0.0, 0.0],
                     [0.0, 0.0, 0.0],
                     [0.0, 0.0, 0.0]], dtype=np.complex128)


t = np.linspace(0, 10, 101)
sol = odeintw(right_part, psi_init, t, args=(0.25,))

sol 将会是一个复杂的 numpy 数组，形状为 (101, 3, 3)，里面存储了解 rho(t) 的结果。第一个索引是时间索引，另外两个索引是 3x3 的矩阵。