von Neumann equation implementation using QuTiP.
QuTiP's mesolve() (master equation solver) integrates the Lindblad master equation. Passing an empty collapse-operator list c_ops=[] disables the dissipator, reducing it exactly to the von Neumann equation:
$$\frac{\mathrm{d}\rho}{\mathrm{d}t} = -\frac{i}{\hbar}[H,\,\rho]$$
Instead of a ket, mesolve() operates on a density matrix rho0 created with qt.ket2dm(). The expectation value of an observable $O$ is now $\langle O\rangle = \mathrm{tr}(\rho\, O)$, computed automatically when O is listed in e_ops. For a pure state with no dissipation, mesolve(c_ops=[]) produces identical trajectories to sesolve().
# run: python3 von_neumann_qutip.py # desc: evolve the density matrix rho=|+><+| under H=(omega/2)*Z with no dissipation import numpy as np import qutip as qt # Initial density matrix: rho = |+><+| psi0 = (qt.basis(2, 0) + qt.basis(2, 1)).unit() rho0 = qt.ket2dm(psi0) omega = 2 * np.pi H = 0.5 * omega * qt.sigmaz() tlist = np.linspace(0, 1.0, 300) # Integrate d/dt rho = -i[H, rho] (c_ops=[] removes the dissipator entirely) result = qt.mesolve(H, rho0, tlist, c_ops=[], e_ops=[qt.sigmax(), qt.sigmaz()]) # Identical to sesolve() output: <X> completes one full oscillation, <Z> stays 0 print(result.expect[0][0], result.expect[0][-1]) # ~1.0 ~1.0 print(result.expect[1][0], result.expect[1][-1]) # ~0.0 ~0.0