Multiple-scale analysis for 2D Hamiltonian?

[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.

 

[Dr. Courtney]

Looks like a 3-D problem to me.

 

[Random137]

Looks like a 3-D problem to me.

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

 

[Dr. Courtney]

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

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

 

[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.

 

[Random137]

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.

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.

 

[Dr. Courtney]

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.

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.

 

[Random137]

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.

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.

 

[Dr. Courtney]

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 <)?

Right.

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.

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.

 

[Random137]

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.

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

Can you conclude anything from these plots?

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.

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?

Thanks in advance.

 

[Dr. Courtney]

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

Can you conclude anything from these plots?

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 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?

Thanks in advance.

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.

 

[Random137]

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.

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?

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.

View attachment 88801

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.

 

[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.

 

[Random137]

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.

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

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.

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

 

[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.

 

[Dr. Courtney]

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.

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.

 

[Random137]

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.

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.

 

[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.

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.

 

[Random137]

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.

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!

 

[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.

 

[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.

 

[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).//

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.

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.

 

[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.

 

[Random137]

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.

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?

 

[Random137]

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.

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?

 

[Dr. Courtney]

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?

Wow, you're an undergraduate. I congratulate you for jumping into the deep end of the pool. On the one hand, there isn't much stuff out there, and it's hard to recommend stuff without a better understanding of your background and goals. Have you asked the faculty member who you are doing this for? It's been a while since I've worked with undergrads on this stuff, and they are likely more aware of the resources than I am.

I made most of my progress in learning this material by doing it. Eventually the stuff I was reading began to make sense with my modelling efforts, but I was well into grad school by the time it did. I was figuring you for a grad student.

I'll poke around for materials on the theory side, but as a practical matter, let me encourage you to sample more of phase space (finer steps in launch angles) and zoom in and try to determine if a trajectory is really on a torus or merely squeezed in between two tori sampling all the phase space between them.

 

[Dr. Courtney]

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?

I don't have an MMa license now. I've used it in the past for certain things, but it is not the best tool for numerical integrations of dynamical systems. However, you have certainly made great progress with it and may not have time or inclination to reinvent your approach in C.

We've always been able to conserve energy with Runge-Kutta using a variable step size. If your faculty advisor recommended MMa, they can probably advise you better at this point. How small do you really need to go with epsilon? 0.01 can probably be made to work with your current method. 0.001 and smaller may require a different approach.

Were you assigned the problem to force you to apply classical perturbation theory? I guess it would help me to understand what you are really trying to accomplish and over what range of epsilon you need to understand the dynamics.

 

[Random137]

Wow, you're an undergraduate. I congratulate you for jumping into the deep end of the pool. On the one hand, there isn't much stuff out there, and it's hard to recommend stuff without a better understanding of your background and goals. Have you asked the faculty member who you are doing this for? It's been a while since I've worked with undergrads on this stuff, and they are likely more aware of the resources than I am.

I made most of my progress in learning this material by doing it. Eventually the stuff I was reading began to make sense with my modelling efforts, but I was well into grad school by the time it did. I was figuring you for a grad student.

I am a grad student but my knowledge in this area is a bit weak, I guess what I really need is something that bridge the gap between undergrad and grad level.

I'll poke around for materials on the theory side, but as a practical matter, let me encourage you to sample more of phase space (finer steps in launch angles) and zoom in and try to determine if a trajectory is really on a torus or merely squeezed in between two tori sampling all the phase space between them.

Can you explain a bit more about "if a trajectory is really on a torus or merely squeezed in between two tori"? Roughly how many trajectories do I need?

 

[Random137]

I don't have an MMa license now. I've used it in the past for certain things, but it is not the best tool for numerical integrations of dynamical systems. However, you have certainly made great progress with it and may not have time or inclination to reinvent your approach in C.

We've always been able to conserve energy with Runge-Kutta using a variable step size. If your faculty advisor recommended MMa, they can probably advise you better at this point. How small do you really need to go with epsilon? 0.01 can probably be made to work with your current method. 0.001 and smaller may require a different approach.

The reason I use Mathematica because it has built in function to solve differential equations numerically, so it doesn't require too much coding and understanding numerical algorithms.

Were you assigned the problem to force you to apply classical perturbation theory? I guess it would help me to understand what you are really trying to accomplish and over what range of epsilon you need to understand the dynamics.

What I am really interested is the case when epsilon = 0.25 (ultimate goal). The limit when epsilon = 1 is well understood (because of the symmetry), the idea is maybe in the limit as epsilon -> 0 (when epsilon is very small) we might be able to do some analytical analysis using some sort of perturbation theory.

The question that we wanted to answer is, for a given initial condition, the trajectory will have some values of <p_x^2> and <p_y^2> (which we still don't know how to predict when epsilon != 1), what happen if we change the initial momentum slightly such that:
$$
p_x(0) \rightarrow p_x(0) + q_x
$$
or
$$
p_y(0) \rightarrow p_y(0) + q_y
$$
by what amount will <p_x^2> and <p_y^2> change? How does the change depends on the parameter q_x and q_y?

From numerical integrations we notice the following when epsilon -> 0:

1. When q_x = 0, <p_x^2> almost do not change except when q_y becomes very large and <p_y^2> changes roughly quadratically with q_y:
   $$
   \Delta <p_y^2> \propto q_y^2
   $$
2. When q_y = 0, both <p_x^2> and <p_y^2> changes roughly quadratically with q_x but with different constant of proportionality.

We would like find a quantitative explanation and maybe a way to calculate those constant of proportionality.

When epsilon = 0.25, we do not have those "clean" dependence on q_x and q_y.

 

[Dr. Courtney]

I think all your questions regarding behavior close to epsilon = 0.25 can be answered with suitable numerical experiments using your current computational technique in MMa.

I would be confident that a variable step size Runge-Kutta would likely be able to investigate the dynamics for epsilon = 0.01 to 0.001 and likely even down to 0.0001.

 

[Random137]

I think all your questions regarding behavior close to epsilon = 0.25 can be answered with suitable numerical experiments using your current computational technique in MMa.

I would be confident that a variable step size Runge-Kutta would likely be able to investigate the dynamics for epsilon = 0.01 to 0.001 and likely even down to 0.0001.

The following contains 500 different trajectories with initial condition satisfying: x(0) = 0, y(0) = 1 with p_x^2(0) + p_y^2(0) = 1

 

[Dr. Courtney]

Very nice.

 

[Random137]

Doing the same thing with a smaller epsilon, it seems like the behaviour is quite different?

How can I extract information from these Poincare plots?

 

[Dr. Courtney]

I think you should compute your long time integrals for the squared momenta for different trajectories, overlay the values with the tori they represent in the graph and look for trends between the tori and resulting integrals.

 

[Random137]

I think you should compute your long time integrals for the squared momenta for different trajectories, overlay the values with the tori they represent in the graph and look for trends between the tori and resulting integrals.

Looking at the 500 trajectories, I calculated <p_x^2> and <p_y^2> for each of the trajectories by evolving the trajectories over a long period of time. Then I arrange the trajectories according to their <p_x^2> value (from smallest to largest).

There seems to be a pattern, similar values of <p_x^2> have similar poincare plots. As <p_x^2> increases, the poincare plots extend further.
http://gph.is/1OOtACG