6

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

是否有任何标准仪器来整合这些方程？

当我以向量作为初始值求解另一个 aquation 时，例如 Schrodinger 方程： 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()
4

2 回答 2

4

我创建了一个scipy.integrate.odeint被调用的包装器odeintw，可以处理诸如此类的复杂矩阵方程。请参阅在 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 矩阵。

于 2014-11-04T23:46:02.033 回答
2

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

于 2015-05-07T12:44:04.753 回答