Turn two ODEs first degree into a Hamiltonian with code

Question

This is more like a coding question, but I thought I'll post it here because of its mathematical foundation.

Two ODEs first degree can describe a Hamiltonian System. The connection between the ODEs and the Hamiltonian System is given by the canonical equations which are:

$ \dot{x}=\frac{\partial H}{\partial x} $

$ \dot{y}=-\frac{\partial H}{\partial y} $

For example:

Two ODEs:

$ \dot{x}= 2y $

$ \dot{y}= -2x $

Results in the following Hamiltonian:

$ H(x,y)=x^2+y^2 $

Is there a way to create a code/algorithm, that can calculate the Hamiltonian from these ODE?

Answer 1

Given the differential equations $$ \eqalign{\dot{x} &= f(x,y)\cr \dot{y} &= g(x,y)\cr} $$ you want $H= H(x,y)$ such that $$\eqalign{\frac{\partial H}{\partial x} &= -g(x,y)\cr \frac{\partial H}{\partial y} &= f(x,y)\cr} $$ That is, $H$ is a scalar potential for the vector field $\langle -g(x,y), f(x,y) \rangle$. Of course you need this vector field to be conservative. If this is the case, first take an antiderivative of the first component with respect to $x$: $$ V(x,y) = \int -g(x,y)\; dx$$ Then we should have $\dfrac{\partial}{\partial x} (H - V) = 0$, so $H - V$ should be a function of $y$ alone. And $$\dfrac{\partial}{\partial y} (H - V) = f(x,y) - \dfrac{\partial V}{\partial y}$$ so $$ H - V = \int \left(f(x,y) - \frac{\partial V}{\partial y} \right) \; dy + constant$$ i.e. a Hamiltonian is $$H = V + \int \left( f(x,y) - \frac{\partial V}{\partial y} \right)\; dy $$

Answer 2

You prob want a subroutine that can take a function as an argument. You might also want a variable step size.

If step size is $\Delta x>0, dy/dx=f(x,y)$ then $y_{i+1}\approx y_i + \Delta x f(x,y), x_{n+1}=x_n + \Delta x$. Then iterate.

That's not the most stable or accurate way of performing an integration. RK4 is much better in both regards.

As to taking derivatives, $df/dx \approx \frac{f(x+\Delta x)-f(x)}{\Delta x}$. You can just use the definition and use a small step size. Too big, obviously the calculation isn't accurately approximating, too small a step size, you could get a similar problem.

To reduce that issue, you want to use the symmetric derivative. $f'(x) \approx \frac{f(x+\Delta x/2)-f(x-\Delta x/2)}{\Delta x}$. This will give you a fairly accurate approximation to the derivative. One caveat among possibly many, it will give you a derivative where the derivative does not exist, cf. $|x|$ at $x=0$.

As to your specific problem. $H=x^2+y^2\implies \partial H/\partial x =2x $ yet it is given that $\dot{x}=2y$. Was that a mistake?