python - 以复数矩阵作为初始值求解 python 中的 ode

Question

我有一个冯诺依曼方程，它看起来像:dr/dt = - i [H, r]，其中 r 和 H 是复数的方阵，我需要使用 python 脚本找到 r(t)。

是否有任何标准仪器可以对这些方程进行积分？

当我求解另一个以向量为初始值的方程式时，例如薛定谔方程:dy/dt = - i H y，我使用了 scipy.integrate.ode 函数('zvode')，但尝试对 von Neumann 方程使用相同的函数时出现以下错误:

**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

我创建了一个 scipy.integrate.odeint 的包装器称为 odeintw可以处理像这样的复杂矩阵方程。参见 How to plot the Eigenvalues when solving matrix coupled differential equations in 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矩阵。