Method 1: Direct computation. The Hamiltonian vector field is $X_H = \left(\frac{\partial H}{\partial p_1}, \ldots, \frac{\partial H}{\partial p_n}, -\frac{\partial H}{\partial q_1}, \ldots, -\frac{\partial H}{\partial q_n}\right)$ in coordinates $(q_1,\ldots,q_n,p_1,\ldots,p_n)$.

The divergence is: $\nabla \cdot X_H = \sum_{i=1}^n \frac{\partial}{\partial q_i}\frac{\partial H}{\partial p_i} + \frac{\partial}{\partial p_i}\left(-\frac{\partial H}{\partial q_i}\right) = \sum_{i=1}^n\left(\frac{\partial^2 H}{\partial q_i\partial p_i} - \frac{\partial^2 H}{\partial p_i\partial q_i}\right) = 0$

by equality of mixed partials (assuming $H \in C^2$). By the divergence theorem for flows, $\nabla \cdot X_H = 0$ implies the flow preserves volume.

Method 2: Symplectic form. The Hamiltonian flow preserves the symplectic form $\omega = \sum dp_i \wedge dq_i$: $\mathcal{L}_{X_H}\omega = 0$ (where $\mathcal{L}$ is the Lie derivative), which follows from Cartan's formula $\mathcal{L}_X\omega = d(\iota_X\omega) + \iota_X(d\omega) = d(dH) + 0 = 0$. Since the volume form is $\Omega = \omega^n/n! = dq_1\wedge dp_1\wedge\cdots\wedge dq_n\wedge dp_n$: $\mathcal{L}_{X_H}\Omega = 0$, proving volume preservation.

Method 3: Jacobian. Let $\Phi_t$ be the flow map. The Jacobian matrix $J(t) = D\Phi_t$ satisfies $\dot{J} = A(t)J$ where $A = \begin{pmatrix} H_{qp} & H_{pp} \\ -H_{qq} & -H_{pq}\end{pmatrix}$. Therefore $\frac{d}{dt}\det(J) = \mathrm{tr}(A)\det(J) = 0$, since $\mathrm{tr}(A) = \mathrm{tr}(H_{qp}) - \mathrm{tr}(H_{pq}) = 0$. Hence $\det(J(t)) = \det(J(0)) = 1$ for all $t$. $\square$