I started reading chapter II.16 of Solving Ordinary Differential Equations I in order to understand nummerical methods more.
There, they state that symplectic methods don't always conserve the first integrals of the system. So for example, if a system is described by some Hamiltonian $H(p,q)$ then the nummerics don't need to conserve it.
What is conserved (Theorem 16.7 of said book or Sanz-Serna 1988) is the quadratic first integral that is also a first integral of the system.
Now, a quadratic first integral is a function of the shape $y^T C y$ for some symmetric $C$ and $y:=(p,q)$. Now, for this to be the first integral we would need to have $H(y) : = y^T C y$ right?
But then, the conservation of $H$ for most numerical cases wouldn't be possible! Take for example the Lennard-Jones potential. There we have: $$H(p, r) \approx \sum_ip_i^2 + \sum_{i,j} \frac{A}{r_{ij}^{12}}-\frac{B}{r_{ij}^6}$$ and we can't conserve this Hamiltonian because $y^TCy$ can only contain terms with $r$ to the power of $0,1,2$ and never negative. What is going on here?
Is it only possible to (nummericly) conserve Hamiltonians of the shape: $$ H(p,q) = \sum_{i,j} a_{ij} p_i p_j + b_{ij}p_{i}q_{j} + c_{ij} q_{i}q_j $$ with the additional constraint that the constants $a,b,c$ need to be symmetric (so $a_{ij} = a_{ji}, b_{ij} = b_{ji}, c_{ij} = c_{ji}$)?
My specific application is to show that: $$H(p,q) = \sum_i |p_i| - \sum_{ij}\frac{1}{r_{ij}}$$ doesn't need to be preserved by leapfrog algorithm and can become unstable (see related question)
What conservation laws can be derived here? What C can I find that will be conserved?
The claim that you cite is that if a Hamiltonian problem has a linear or quadratic first integral, then the numerical method will preserve it. An example is the angular momentum $\sum \vec p_i\times\vec q_i$ in rotationally invariant systems like a solar system simulation and also the Lennard-Jones system. Probably also your relativistic system, as it has no preferred direction or an-isotropic terms.
A non-quadratic first integral is preserved in the fashion that a perturbation of it is preserved exactly by the numerics. But that is only true where the perturbation series converges, which excludes the singularities of the system. For a truncation of the perturbation series, the corresponding high-order preservation is valid only under exclusion of a neighborhood of the singularities. For a higher degree truncation the radius gets smaller.
In Hairer/Lubich/Wanner (2003) "Geometric numerical integration as illustrated by the Störmer/Verlet method" it is put forward that $H^{(k)}\sim k!MR^{-k}$ at distance $R$ from a pole or singular set, then $h\ll\frac{\pi R}N$ is necessary for an expansion up to degree $N$ so that the perturbation is valid for solutions that come no closer than that distance $R$ to the singular set. They also write that getting concrete, quantitative results is very difficult and lengthy.