4
$\begingroup$

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?

$\endgroup$
6
  • 1
    $\begingroup$ 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}$. $\endgroup$
    – Winther
    Feb 3, 2019 at 14:11
  • 2
    $\begingroup$ 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. $\endgroup$
    – Winther
    Feb 3, 2019 at 18:09
  • $\begingroup$ 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. $\endgroup$ Feb 4, 2019 at 10:58
  • 1
    $\begingroup$ 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. $\endgroup$ Feb 7, 2019 at 22:17
  • 1
    $\begingroup$ 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. $\endgroup$ Feb 8, 2019 at 9:11

1 Answer 1

1
+150
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ What do $k$ and $M$ refer to? $\endgroup$ Feb 9, 2019 at 15:18
  • $\begingroup$ $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. $\endgroup$ Feb 9, 2019 at 15:29

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .