1) Pick a power law scaling for the density that I wish to study
2) Solve for the distribution of energies that is expected for an equilibrium solution with this scaling
3) Generate some random energies distributed according to this power law
4) Use these energies to generate positions, a) requiring that the orbits are equally destributed in time over a quarter orbit and b) using potential energies consistent with the expected power law scaling
5) generate velocities using v=sqrt(2E-2phi(x))
6) reflect the (x,v) pairs over the x and v axes to generate symmetric initial conditions
7) run the code to see if the scaling remains the same, or if it wanders to the attractor
Here is a plot of the initial phase space configuration for a $\rho \propto x^{-0.8}$.
This initial configuration, which has been constructed to be in equilibrium, actually explodes, the particles quickly go off to \pm inifinity essentially. Below is an intermediate phase space configuration, that to me seems clearly wrong:
And a little while later:
My next step, since we have an independent code, was to see if it is behaving in the same way. Here is the initial condition I fed to Jerry's code:
And here is what it looks like a little while later:
So they are both suffering from a similar problem, I think, although for Jerry, it is the velocities that are exploding (which must EVENTUALLY cause the system to explode, but I didn't run it that far). Some of the other power laws are not behaving quite so badly, but in terms of examining whether an equilibrium is stable or unstable (and running off to the attractor solution), I feel we need to understand what is going on here a little better.
I will email the initial condition file that causes the bad behavior, this time it is definitely reproducible.
No comments:
Post a Comment