在python中以复数矩阵为初始值求解ode

Posted 2023-03-12

【中文标题】在python中以复数矩阵为初始值求解ode

【英文标题】：Solve ode in python with complex matrix as initial value

【发布时间】：2014-12-31 19:06:20

【问题描述】：

我有一个冯诺依曼方程，如下所示： 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()

【参考方案1】：

QuTiP 有一些很好的积分器可以做到这一点，使用诸如主方程和 Lindblad 阻尼项之类的东西。 QuTiP 本身只是 scipy.odeint 的一个薄包装，但它使许多机制变得更好，特别是因为它有很好的文档。

【讨论】：

貌似在张量运算方面可以替代numpy！很酷！

【参考方案2】：

我创建了一个名为odeintw 的scipy.integrate.odeint 包装器，它可以处理诸如此类的复杂矩阵方程。有关矩阵微分方程的另一个问题，请参阅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 矩阵。