Schrödinger equation implementation using QuTiP.
QuTiP's sesolve() (schrödinger equation solver) integrates the equation directly in state-vector form. It takes the Hamiltonian H, an initial ket psi0, and a time array tlist, then returns a list of ket vectors $\lvert\psi(t)\rangle$ at each time step — the numerical solution to:
$$i\hbar\frac{\mathrm{d}}{\mathrm{d}t}\lvert\psi\rangle = H\lvert\psi\rangle$$
The e_ops argument is a list of operators whose expectation values $\langle O\rangle = \langle\psi(t)\lvert O\lvert\psi(t)\rangle$ are computed at every step and stored in result.expect.
# run: python3 schrodinger_qutip.py # desc: evolve a qubit in |+> under H = (omega/2)*Z and track <X> precessing on the Bloch sphere import numpy as np import qutip as qt # Initial state: |+> = (|0> + |1>) / sqrt(2) psi0 = (qt.basis(2, 0) + qt.basis(2, 1)).unit() # Hamiltonian: Larmor precession around z-axis at angular frequency omega omega = 2 * np.pi # one full Bloch-sphere rotation over t in [0, 1] H = 0.5 * omega * qt.sigmaz() tlist = np.linspace(0, 1.0, 300) # Integrate i * d/dt |psi> = H|psi> result = qt.sesolve(H, psi0, tlist, e_ops=[qt.sigmax(), qt.sigmaz()]) # <X> oscillates as cos(omega*t); <Z> stays 0 (precession in the xy-plane) print(result.expect[0][0], result.expect[0][-1]) # ~1.0 ~1.0 (full period) print(result.expect[1][0], result.expect[1][-1]) # ~0.0 ~0.0