我有一个冯诺依曼方程,它看起来像: 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()
最佳答案
我创建了一个 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矩阵。
关于python - 以复数矩阵作为初始值求解 python 中的 ode,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/26742818/