Numerical solution of Hamiltonian systems

[fabsilfab]

The question is very general and could belong to another topic, but here it is.
Suppose one wants to solve the set of differential equations $$ \frac{\partial x}{\partial t}=\frac{\partial H(x,p)}{\partial p},$$ $$\frac{\partial p}{\partial t}=-\frac{\partial H(x,p)}{\partial x},$$ with some numerical method, e.g. a Newton-Rhapson iteration. Let us discretise the above system in the most straightforward way, $$ \frac{x^{n+1}-x^n}{\Delta t}=\frac{H(x,p^{n+1})-H(x,p^n)}{p^{n+1}-p^n}, $$ $$ \frac{p^{n+1}-p^n}{\Delta t}=\frac{H(x^{n+1},p)-H(x^n,p)}{x^{n+1}-x^n}.$$ From the numerical point of view, this system is clearly ill-defined when either $x$ or $p$ do not change over a finite time step $\Delta t$. Looking carefully at the equations, it is clear that the terms on the RHS represent an incremental ratio. In the limit of $\Delta x \rightarrow 0$ or $\Delta p\rightarrow 0$, such ratios reduce to the exact (analytic) derivative. This must hold, otherwise we would encounter non-physical situations. For example, suppose the velocity is constant, thus the position is changing: if the RHS in the position equation did not reduce to the analytic derivative, the position update would be undefined (the denominator being zero would make the equation blow up; or, the numerator being zero would mean no position change).

Now the question: from the point of view of implementing the numerical solution of such system in a code, there must be some sort of switch. When the difference between position or momentum becomes zero, the RHS of the corresponding equation must reduce to the analytic derivative of the Hamiltonian $H$. What I don't get is how to effectively due this, considering that in the most general case, these equations are implicit.

Any thoughts?

 

[DrClaude]

I don't see why you want to discretize the derivative of the Hamiltonian. What is usually done is that you calculate analytically the derivative of the Hamiltonian to get the equations of motion, and it is these equations of motion that you discretize.

 

[fabsilfab]

I don't see why you want to discretize the derivative of the Hamiltonian. What is usually done is that you calculate analytically the derivative of the Hamiltonian to get the equations of motion, and it is these equations of motion that you discretize.

Good point, but Hamiltonian-preserving numerical schemes are based on the direct discretisation of the Hamiltonian. Taking the analytic derivative and then discretising does not ensure that the Hamiltonian is preserved.

 

[DrClaude]

Good point, but Hamiltonian-preserving numerical schemes are based on the direct discretisation of the Hamiltonian.

Really? I'm surprised. Do you have examples of such methods? It is hard for me to see why doing the derivative numerically would give any sort of advantage.

Taking the analytic derivative and then discretising does not ensure that the Hamiltonian is preserved.

That's not correct. All the symplectic methods I know are based on analytical $H_p$ and $H_q$.

For more information, you can check out the book Geometric Numerical Integration.

 

[fabsilfab]

Really? I'm surprised. Do you have examples of such methods? It is hard for me to see why doing the derivative numerically would give any sort of advantage.That's not correct. All the symplectic methods I know are based on analytical $H_p$ and $H_q$.

For more information, you can check out the book Geometric Numerical Integration.

True. But symplectic methods do not necessarily preserve the Hamiltonian (e.g., energy). They do, however, preserve it to some high precision. But not exactly. Also, you are probably referring to the special case of separable Hamiltonians. For more general case, it is even more difficult to formulate symplectic methods.