用复杂矩阵作为初始值在python中求解ode [英] Solve ode in python with complex matrix as initial value

问题描述

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

I have a von Neumann equation which looks like: dr/dt = - i [H, r], where r and H are square matricies of complex numbers and I need to find r(t) using python script.

是否有任何标准工具可以集成此类方程式?

Is there any standart instruments to integrate such equations?

当我用向量作为初始值求解另一个方程时，如薛定inger方程: dy/dt =-i H y ，我使用了scipy.integrate.ode函数('zvode')，但是尝试对冯·诺依曼方程使用相同的函数会出现以下错误:

When I was solving another aquation with a vector as initial value, like Schrodinger equation: dy/dt = - i H y, I used scipy.integrate.ode function ('zvode'), but trying to use the same function for von Neumann equation gives me the following error:

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

推荐答案

我创建了 odeintw 可以处理诸如此类的复杂矩阵方程.参见如何绘制在PYTHON中求解矩阵耦合微分方程时的特征值?另一个涉及矩阵微分方程的问题.

I have created a wrapper of scipy.integrate.odeint called odeintw that can handle complex matrix equations such as this. See How to plot the Eigenvalues when solving matrix coupled differential equations in PYTHON? for another question involving a matrix differential equation.

这是代码的简化版本，显示了如何使用它. (为简单起见，我从您的示例中删除了大多数常量.)

Here's a simplified version of your code that shows how you could use it. (For simplicity, I got rid of most of the constants from your example).

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将是形状为(101, 3, 3)的复杂numpy数组，其中包含解决方案rho(t).第一个索引是时间索引，其他两个索引是3x3矩阵.

sol will be a complex numpy array with shape (101, 3, 3), holding the solution rho(t). The first index is the time index, and the other two indices are the 3x3 matrix.

这篇关于用复杂矩阵作为初始值在python中求解ode的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！