# Multiple-scale analysis for 2D Hamiltonian?

Tags:
1. Sep 16, 2015

### Random137

I came across a technique called "multiple-scale analysis" https://en.wikipedia.org/wiki/Multiple-scale_analysis where the equation of motion involves a small parameter and it is possible to obtain an approximate solution in the time scale of $$\epsilon t$$.

I am wondering if it is possible to apply similar technique to a 2D Hamiltonian, in particular:
$$H = \frac{p_x^2 + p_y^2}{2 m} + c \sqrt{x^2 + \frac{y^2}{\epsilon}}$$
with the following equation of motion:
$$\ddot{x} = - \frac{c x}{m \sqrt{x^2 + \frac{y^2}{\epsilon}}}$$
$$\ddot{y} = - \frac{c y}{m \epsilon \sqrt{x^2 + \frac{y^2}{\epsilon}}}$$
where $$\epsilon$$ is the small parameter.

Last edited: Sep 16, 2015
2. Sep 16, 2015

### Dr. Courtney

Looks like a 3-D problem to me.

3. Sep 16, 2015

### Random137

Why? The Hamiltonian describes a particle moving in 2D?

4. Sep 16, 2015

### Dr. Courtney

The bottom equation of motion was z double dot before you edited it.

5. Sep 16, 2015

### Dr. Courtney

I'm not sure why you need to. If you can't simply beat it into submission with Runge-Kutta using a variable step size, then I think you could change coordinate systems and have something easier to integrate numerically.

6. Sep 16, 2015

### Random137

Ideally, I would like to estimate the following quantities:
$$\langle p_x^2 \rangle = \lim_{T \rightarrow \infty} \frac{1}{T} \int_0^T dt p_x^2 (t)$$
$$\langle p_y^2 \rangle = \lim_{T \rightarrow \infty} \frac{1}{T} \int_0^T dt p_y^2 (t).$$

So I need a way to see how those quantities depends (at least approximately) on the initial conditions.

7. Sep 16, 2015

### Dr. Courtney

OK, but you can approach that numerically. A numerical approach would also allow you to look at Poincare plots to see how the motion is confined to tori (or not) in phase space. I used to run Runge-Kutta routines to solve problems like this over a range of initial conditions to determine what portions of phase space were integrable for certain values of the parameters.

8. Sep 16, 2015

### Random137

Just to confirm, when you say Poincare plots, do you mean plotting x vs p_x when the y = 0 and p_y > 0 (or <)?

I can solve the trajectories numerically quite easily on Mathematica and create those Poincare plots but I don't really know how to do analysis with those plots and I don't seem to be able to find some straightforward methodology in the literatures.

9. Sep 17, 2015

### Dr. Courtney

Right.

If all the intersections of a given trajectory trace out a closed curve in the Poincare plot, then you are observing the result of that trajectory being confined to a torus, demonstrating the motion with that set of initial conditions is integrable. There must be some constant of motion that keeps the trajectory confined to the torus.

On the other hand, if of the intersections of a given trajectory seem randomly scattered throughout a region of the Poincare plot, then that trajectory is chaotic (not integrable). If the motion of all the trajectories for given values of the parameters go through all of the phase space that is energetically allowed, then the motion is completely chaotic. More likely, phase space is divided into some chaotic and some integrable regions by existing areas where the motion is still confined to tori.

Usually, an early approach to addressing Hamiltonians like the one you presented is to determine (for a given set of parameter values) which areas of phase space are integrable and which are chaotic. I expect there will be some parameter values for which most of the phase space is integrable and other values for which most of the phase space is chaotic.

10. Sep 17, 2015

### Random137

Here are some poincare plots I created when y = 0 and p_y > 0:

Can you conclude anything from these plots?

In this 2D Hamiltonian, the phase space is 4D, what's a good way to "sample" the phase space to see which areas are integrable and which are chaotic?

11. Sep 17, 2015

### Dr. Courtney

There do seem to be in tact tori. But I like to let the integrations run longer so the points fill out the tori more fully.

In the systems I've studied, I would start the trajectories at the origin and sweep through the possible launch angles systematically. I would then plot all the resulting Poincare plot data sets on the same graph. That's how the Poincare plot below was created. This plot shows that phase space is mixed in this case. There are areas where tori exist and areas where the orbits sample a broader region of phase space and are only confined between tori.

12. Sep 18, 2015

### Random137

I ran the integrations for 100 times longer than the previous plots and here are the results:

How can I find the other constant of the motion?

Just to clarify, do you mean picking x(0) = y(0) = 0 and for example set the magnitude of the momentum vector (p) as 1 and then try and generate angles from 0 to 2 pi uniformly?

$$p_x = p \cos \theta$$
$$p_y = p \sin \theta$$
This will make sure all the different trajectories have the same energy.

Last edited: Sep 18, 2015
13. Sep 18, 2015

### Dr. Courtney

So now you can see that the trajectories are staying confined to tori that do not trace out one big closed curve on the boundary of the energetically allowed phase space, but that the tori are really more complicated and interesting than a precessing orbit.

You approach is fine, but be aware that the nature of the dynamics likely changes with the energy.

I usually pick one energy to study at a time. Sweeping the initial momentum through the angles usually samples a representative portion of the phase space available at that given energy and generates informative Poincare plots.

14. Sep 18, 2015

### Random137

I did 20 different trajectories with x(0) = 0, y(0) = 1 with p_x^2 + p_y^2 = 1:

The potential that I am studying is a homogeneous function, so we can use mechanical similarity, so I think studying one energy is enough.

Last edited: Sep 18, 2015
15. Sep 18, 2015

### Dr. Courtney

That's pretty good work in such a short time with so little guidance.

One way to proceed from here would be to go ahead and compute (a numerical approximation for)
$$\langle p_x^2 \rangle = \lim_{T \rightarrow \infty} \frac{1}{T} \int_0^T dt p_x^2 (t)$$
$$\langle p_y^2 \rangle = \lim_{T \rightarrow \infty} \frac{1}{T} \int_0^T dt p_y^2 (t).$$

for each of the tori shown in the figure. You approximations will be most accurate if you end the integration from about the same place on the torus as you start it from. As T goes to infinity, the trajectory will just keep winding around the torus, so the more accurate estimates will be one that samples different parts of the torus evenly.

Another approach might be to first compute analagous Poincare plots for different values of epsilon. I expect the structure of the tori will be simpler as you increase epsilon from 0.25 to 1.0 and more complex (the tori will break up) as you reduce epison from 0.25 to 0.01 or 0.001.

16. Sep 18, 2015

### Dr. Courtney

I should explain why.

For epsilon = 1, the potential has a cylindrical symmetry (only depends on r in polar coordinates). Consequently, the angular momentum is conserved. Since there are two constants of motion (angular momentum and energy) and two degrees of freedom, all the motion is integrable (confined to tori). In fact, in this case it is trivial.

17. Sep 18, 2015

### Random137

For epsilon = 1, <p_x^2> = <p_y^2> (Bertrand's theorem) and using Virial Theorem, I can calculate <p_x^2> and <p_y^2> from energy conservation.

18. Sep 18, 2015

### Random137

Let me elaborate on this:

According to Bertrand's theorem, there are only 2 central potentials - k / r (Kepler) and k r^2 (harmonic) with all bound orbits which are also closed orbits. Therefore, one can assume/prove <p_x^2> = <p_y^2> for other central potentials.

According to Virial Theorem, for homogeneous potentials, there is a relationship between the average kinetic energy and average potential energy:
$$<T> = \alpha <V>$$
where \alpha is some constant which depends on the potential.

Therefore, the total energy (which is a constant of the motion) can be expressed purely in terms of the average kinetic energy:
$$E = <E> = <T> + <V> = \frac{1+\alpha}{\alpha}<T>$$

and the average kinetic energy is given by:
$$<T> = \frac{1}{2 m} (<p_x^2> + <p_y^2>) = \frac{1}{m} <p_x^2>$$
the last step using <p_x^2> = <p_y^2>. Therefore, <p_x^2> can be calculated from the total energy.

19. Sep 18, 2015

### Random137

I created the same plots for epsilon = 0.02, I have some numerical issues trying to reduce epsilon even more as it gets more difficult keeping the energy "constant" (within some reasonable amount of error).

It does seem to get simpler as you expected!

20. Sep 18, 2015

### Dr. Courtney

Great stuff. I think the outermost ring may be an ergotic region between two tori and not really be a tori itself. This may also be true for some of the inner rings.

21. Sep 18, 2015

### Dr. Courtney

//I created the same plots for epsilon = 0.02, I have some numerical issues trying to reduce epsilon even more as it gets more difficult keeping the energy "constant" (within some reasonable amount of error).//

Cases like this are where I've found Runge-Kutta with variable step size to be useful. Running the code in C was fairly fast, even on the 80386 with 4 MB RAM that I used back in the early 90s. I was living large when I bought by first Pentium. The ThinkPad on my lap these days blasts through all these.

22. Sep 18, 2015

### Random137

Currently, I am using something called Symplectic Partitioned Runge Kutta method in Mathematica (which I don't really understand how it works), it allows me to set invariant quantities and it will try and keep the invariant quantities almost constant throughout the integration over many time steps.

23. Sep 18, 2015

### Dr. Courtney

That method seems to be assuming that the given trajectory is either integrable or close to an integrable trajectory. In reality, this is not known a priori in your case.

I've always preferred a variable step size Runge-Kutta. You can write a bit of extra code to double check that energy is being conserved (to within your acceptable error) and make the step size (or error condition) smaller if it isn't.

Also, it is not clear (to me) if your code is integrating the four first order diff eqs that result from Hamiltonian mechanics or the two second order diff eqs that result from Newton's laws. When the numerical side gets squirrely, the four first order differential equations will be better behaved.

Any modern computer will be fast enough to implement an adaptive step size Runga-Kutta without making unwarranted assumptions. The ones in Numerical Recipes are straight forward and have worked well for us in a number of problems on both sides of the transitions to chaos and other squirrelly business.

24. Sep 18, 2015

### Random137

My understanding in this area of Hamiltonian dynamics (ergodicity, invariant torus, symplectic form etc) is lacking. I tried to read about this topic but most things I found are written in some rigorous mathematical language involving things like manifolds which I am not really familiar with. Do you have any references that are more suitable for an undergraduate in physics?

25. Sep 18, 2015

### Random137

I tried to implement the velocity-verlet algorithm by myself and I had the same problem with energy conservation when epsilon is small. Would Runge-Kutta fix such problem?

Do you happen to use Mathematica?