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?

• I haven't looked at the question to carefully, just one observation: from what I see it will be unstable no-matter what algorithm you use due to the $1/r_{ij}$ term (if the charge carriers can get arbitrarily close to each other). The standard way of going around this issue is to use a force softening $1/r_{ij} \to 1/\sqrt{r_{ij}^2+\epsilon^2}$. Feb 3, 2019 at 14:11
• As for your question: it is not true that "Now, for this to be the first integral we would need to have $H = y^TCy$ right?". For example if you have the system $\dot{y} = f(y)$ then $(y^TCy)' = 2y^TCf(y) = 0$ will be conserved in time if $y^TCf(y) = 0$ and for a non-linear $f(y)$ this will not correspond to the EOM for a quadratic Hamiltonian. Feb 3, 2019 at 18:09
• This is indeed an great point. Ofcourse, H is not the only conserved quantity when we transverse the EOM. In particular any C such that $y^T C f(y) = 0$ corresponds to a conserved quantity. Thank you it does help clear up the confusion a lot. Feb 4, 2019 at 10:58
• Not sure if this will help, but briefly reading through this made me think of how a quadratic programming problem is considered ill-formed when the symmetric matrix isn't positive semidefinite - the optimization problem no longer has a unique solution. Feb 7, 2019 at 22:17
• The claim is that if a problem has a linear or quadratic first integral, then the numerical method will preserve it. An example is angular momentum in rotationally invariant systems. A non-quadratic first integral is preserved in the fashion that a perturbation is preserved exactly by the numerics, but only 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. Feb 8, 2019 at 9:11

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.
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.
• What do $k$ and $M$ refer to? Feb 9, 2019 at 15:18
• $k$ is the order of the derivative, $M$ some scale factor of the function. Think of $f(x)=M/x$, then $f(R+s)=\sum [k!MR^{-1-k}]\frac{(-s)^k}{k!}$. Note that in the paper that reasoning is presented more as a hand-waving guideline than as an exact result. Feb 9, 2019 at 15:29