Say I have a Hamiltonian system of $N$ particles in 3D-3V phase space. I'm using some sort of update scheme taking the system from $t^{n-1}$ to $t^{n}$ to $t^{n+1}$. I want to know if the update scheme is symplectic, volume preserving, or neither, but I want to test this in an experimental way.* If it's symplectic, then the 2-form $\omega = dq \wedge dp = \sum_{i}^{N}{dq_i \wedge dp_i}$ is preserved over time (minor question here, are $dq$ and $dp$ vectors or scalars? I'm assuming vectors representing individual particles in phase space, but I could see an argument made for simply iterating over $3N$ $q,p$ pairs rather than $N$ vectors), which in turn implies phase space volume is preserved. My question is, how do I test this computationally? My current thought:
Compute $dq_i^n = q_i^{n}-q_i^{n-1}, dp_i^{n} = p_i^{n}-p_i^{n-1}$, which are vectors with three components. Taking the 2-form is, in the case of vectors in $\mathbb{R}^3$, just the cross product (from what I understand it's a bit more complicated than "just the cross product" in that the result is a two-vector or something to that effect, but I'm only concerned in having something to compare over the simulation, so as long as I'm on the right track here I'm content). Sum up all $dq_i^n \wedge dp_i^n$. Do the same for $dq_i^{n+1} \wedge dp_i^{n+1}$. If the scheme is symplectic, then the difference between these two values should be bound over the simulation (unsure if I'm comparing individual components or the norms of those values). For phase space volume preservation, simply compute $dV_i^n = dq_{1i}^ndq_{2i}^ndq_{3i}^ndp_{2i}^ndp_{1i}^ndp_{3i}^n$, take the product $\Pi_i^N{dV_i^n}$ up over all $i$, do the same for $dV_i^{n+1}$, compare over timesteps, the difference should be bound over the simulation.
Frustratingly, I cannot find any resources to confirm or deny this.** So I suppose my question is in two parts, first am I right/where am I wrong, second is what are some good resources for looking more into this? Incidentally, symplecticity $\implies$ phase space volume preservation, but whenever I try to prove that the above follows that rule I very quickly get lost in about a million cross terms, which either means I need to get back to it or I'm on the wrong track.
*I know if $\Omega$ is my update matrix, then $\Omega^TJ\Omega = J$ implies symplecticity and $det(\Omega)=1$ implies volume preservation. This is great, but I'd like to have a practical test for my code. I've been using the Hamiltonian as my litmus test, making sure it's bound over time, but from what I understand this doesn't necessarily track phase space volume preservation or symplecticity, just total energy in the system.
**Every paper, website, and youtube video I've seen has shown the same picture of a box or a smiley face or whatever in phase space being stretched with the area/volume being preserved over time, which is helpful in getting an intuitive idea, not helpful in actually solving the problem.