How to solve for Hamiltonian gradient?

[kulimer]

I am trying to understand how Hamiltonian gradient works.

H(q,p)=U(q)+K(p)​

U(q): potential energy
K(p): kinetic energy
q: position vector
p: momentum vector
both p and q are functions of time
H(q,p): total energy

\frac{d{{q}_{i}}}{dt}=\frac{\partial H}{\partial {{p}_{i}}}
\frac{d{{p}_{i}}}{dt}=-\frac{\partial H}{\partial {{q}_{i}}}
​

Now, I am trying to solve this (the technical name is called leapfrog method)

{p}_{i}(t+\varepsilon /2)={p}_{i}(t)-(\varepsilon /2)\frac{\partial U}{\partial {q}_{i}}(q(t))

{q}_{i}(t+\varepsilon )={q}_{i}(t)+\varepsilon \frac{{p}_{i}(t+\varepsilon /2)}{m}

{p}_{i}(t+\varepsilon )={p}_{i}(t+\varepsilon /2)-(\varepsilon /2)\frac{\partial U}{\partial {q}_{i}}(q(t+\varepsilon ))​

 

[t_evans]

Hi Kulimer,

It looks to me like your doing classical rather than quantum mechanics. But anyway, basically the equations you've got there are just expanding the derivatives in Hamilton's equations in a finite difference approximation. Anyway, the way to solve all such problems is iteratively. From the particular finite difference approximation taken, (and assuming you know U(q) at all q), you need a pair of boundary conditions, to specify vectors (p_i,q_i) at a particular point in time. Let's say (p_i(0),q_i(0)). Then start off at this time and use it with the finite difference equations to generate (p,q) at later times. So consider,

p_i(0+\epsilon/2)=p_i(0)-(\epsilon/2)*\frac{\partial U}{\partial q_i}(q_i(0))

Note that with the boundary condition you know the complete left hand side of this equation. Once you know this, you can solve the second equation at 0+\epsilon, then the third equation. Then you can use the solutions of these equations to solve the set of equations at the next step in time, and keep generating terms for as long as you like.

Note this isn't the sort of thing you'd want to do by hand. A few loops or a spreadsheet if you've not done programming before.

Out of interest, why is finite difference in momentum half that it is in position? If \epsilon it would be easier just to have do momentum and position on the same spacing, but then there may be a good reason.

Hope this helps.

 

[kulimer]

Thanks for your answer.
Should I move this thread to classical physics? Maybe you know how?

Out of interest, why is finite difference in momentum half that it is in position?

This is actually a discrete estimation of Hamiltonian using leapfrog.

As pointed out here on page 8
http://www.cs.toronto.edu/~radford/ftp/ham-mcmc.pdf

But it didn't explicitly say what are the partial derivatives.
I simulated 20 numbers, the pair of numbers (position, momentum), they keep on getting larger, which is weird.

Does this look right to you? I am no physics genius on this. Hope you could let me know if I am doing the right thing.

First column: position
Second column: momentum

[1,] 0.2194610 0.7644557
[2,] 0.4586734 0.8661759
[3,] 0.7391665 1.0458519
[4,] 1.0861846 1.3196545
[5,] 1.5309592 1.7122261
[6,] 2.1135202 2.2588980
[7,] 2.8862981 3.0088708
[8,] 3.9188427 4.0296419
[9,] 5.3040832 5.4130808
[10,] 7.1666912 7.2836969
[11,] 9.6743014 9.8098458
[12,] 13.0525987 13.2188808
[13,] 17.6056299 17.8176151
[14,] 23.7431677 24.0199348
[15,] 32.0175907 32.3840485
[16,] 43.1735969 43.6627267
[17,] 58.2152267 58.8710502
[18,] 78.4962270 79.3777683
[19,] 105.8418877 107.0284855
[20,] 142.7133183 144.3117664

 

[t_evans]

Yeah, I can see it could be more stable if the step size were smaller. It would especially make sense if the potential was a slower in q than in p. The thing is, for most introductory examples, computation time isn't really a problem so you can take /epsilon to be really small and do more steps then it is easier (and potentially stabler) to use the Euler method, because it's clearer precisely what's going on (it really is just a finite approximation of the derivatives in Hamilton's equations), it's also then easier to plot (p,q) behaviour.

I do concede though, as you get to more complicated examples (especially if you have more dimensions to worry about), computation speed may become an issue at some point and such optimisations may be useful (i.e. for precisely the sort of problems discussed in the rest of the article!)

As for making sense of the results, it's rather going to depend on the form of U! If you've used the Harmonic potential as in the article, so
U = q^2 /2
Then this is the incorrect behaviour, you would expect oscillations (like in the article, where (p,q) forms a circle, if you've got the other sign however, so
U = -q^2 /2
Then the exponential growth is the correct behaviour. So what is the potential in question?