Resuelve oda en python con matriz compleja como valor inicial

numpy scipy

Question

Tengo una ecuación de von Neumann que se ve como: dr / dt = - i [H, r] , donde r y H son matrices cuadradas de números complejos y necesito encontrar r (t) usando la secuencia de comandos python.

¿Hay algún instrumento estándar para integrar tales ecuaciones?

Cuando estaba resolviendo otra aquation con un vector como valor inicial, como la ecuación de Schrodinger: dy / dt = - i H y , utilicé la función scipy.integrate.ode (''zvode''), pero tratando de usar la misma función para von Neumann La ecuación me da el siguiente 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**

Aquí está el código:

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

QuTiP tiene algunos buenos integradores para hacer justamente esto, usando cosas como ecuaciones maestras y términos de amortiguación de Lindblad. QuTiP es solo una capa delgada alrededor de scipy.odeint, pero hace que muchos de los mecanismos sean mucho más agradables, sobre todo porque tiene una gran documentación.

Answer 2

He creado un contenedor de scipy.integrate.odeint llamado odeintw que puede manejar ecuaciones matriciales complejas como esta. Vea Cómo graficar los autovalores al resolver ecuaciones diferenciales matriciales en PYTHON? para otra pregunta que implica una ecuación diferencial matricial.

Aquí hay una versión simplificada de su código que muestra cómo podría usarlo. (Por simplicidad, me deshice de la mayoría de las constantes de su ejemplo).

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 será una matriz numpy compleja con forma (101, 3, 3) , sosteniendo la solución rho(t) . El primer índice es el índice de tiempo, y los otros dos índices son la matriz de 3x3.