numpy - Solve ode in python with complex matrix as initial value -

numpy - Solve ode in python with complex matrix as initial value -

i have von neumann equation looks like: dr/dt = - [h, r], r , h square matricies of complex numbers , need find r(t) using python script.

is there standart instruments integrate such equations?

when solving aquation vector initial value, schrodinger equation: dy/dt = - h y, used scipy.integrate.ode function ('zvode'), trying utilize same function von neumann equation gives me next 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**

here code:

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() , 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 homecoming 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) homecoming (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) homecoming r, t0 def main(): r, t0 = create_integrator() t1 = 10 ** -6 dt = 10 ** -11 e, g, u = integrate(r, t0, t1, dt) main()

i have created wrapper of scipy.integrate.odeint called odeintw can handle complex matrix equations such this. see how plot eigenvalues when solving matrix coupled differential equations in python? question involving matrix differential equation.

here's simplified version of code shows how utilize it. (for simplicity, got rid of of constants example).

import numpy np 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) homecoming (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 complex numpy array shape (101, 3, 3), holding solution rho(t). first index time index, , other 2 indices 3x3 matrix.

python numpy scipy ode quantum-computing