Sampling a phase space uniformly for a given energy?

Question

The energy for a 1D harmonic oscillator can be written: $$ E = \frac{1}{2 m} p^2 + \frac{m \omega^2}{2} x^2 $$ where x is position and p is the momentum.

I would like to sample phase space $(x,p)$ uniformly (generating random numbers), in the case where: $$ \frac{1}{2 m} = \frac{m \omega^2}{2}, $$ trajectories with the same energy $E'$ in phase space lie in a circle centered at the origin with radius $\sqrt{2 m E'}$. Therefore, we can use the trick in this answer: uniform distribution on unit circle

In my particular problem, the energy is given by: $$ E = \frac{p_x^2 + p_y^2 + p_z^2}{2 m} + V(x,y,z) $$ where $V(x,y,z)$ is the potential energy of Cartesian co-ordinates. I would like to generate random numbers $(x,y,z,p_x,p_y,p_z)$ that are uniformly distributed in allowed phase space all with energy $E'$

Answer 1

Solving the one-dimensional problem is a first step toward solving the three-dimensional problem described in the question and comments. Without loss of generality, it can be written after renormalization as

1 == p^2 + Sqrt[x^2/4]

Noting that phase space is mirror-symmetric about both axes, we solve for the phase curve in the quadrant p > 0, x > 0.

xphase = x /. Assuming[x > 0, Solve[%, x]][[1, 1]];
Plot[xphase, {p, 0, 1}, AxesLabel -> {p, x}]

To "sample phase space uniformly", we compute the arc length of this curve and sample it uniformly.

dsdp = Sqrt[1 + (D[xphase, p])^2];
arc = Simplify[s[p] /. DSolve[{s'[p] == dsdp, s[0] == 0}, s[p], p][[1]]]

with solution

(4*p*Sqrt[1 + 16*p^2] + ArcSinh[4*p])/8

Knowing arc as a function of p, we compute the inverse function in the usual manner.

pinv[s_] := p /. FindRoot[arc - s, {p, .5}]

Finally, following 63805, we obtain

ListPlot[{pinv[#], xphase /. p -> pinv[#]} & /@ 
  RandomVariate[UniformDistribution[{0, arc /. p -> 1}], 100], AxesLabel -> {p, x}]

Presumably, the complete problem can be solved by generating a uniform distribution on the unit sphere in six dimensions and mapping it to the surface generated by the Hamiltonian described in the question and comments.

Answer 2

Alternative 2D Solution

A simpler and more readily generalizable alternative to the solution presented in my previous Answer can be obtained by noting that a UniformDistribution on the curve shown in the first figure is equivalent to a ProbabilityDistribution in p weighted by dsdp,

ProbabilityDistribution[dsdp/arcmax, {p, 0, 1}]

as can be can be demonstrated by

ListPlot[{#, 2 (1 - #^2)} & /@ RandomVariate[ProbabilityDistribution[dsdp/arcmax, {p, 0, 1}], 100],
  AxesLabel -> {p, x}]

which reproduces the second plot in the earlier Answer.

6D Solution (Almost)

Without loss of generality, the full problem can be written after renormalization as

1 == px^2 + py^2 + pz^2 + Sqrt[x^2/4 + y^2/4 + z^2]

Subsequent calculations can be simplified considerably by transforming momentum space to spherical coordinates and configuration space to cylindrical coordinates.

1 == p^2 + Sqrt[z^2 + r^2/4]

where p^2 = px^2 + py^2 + pz^2 and r^2 = x^2 + y^2. Noting that phase space is mirror-symmetric about all axes, we solve for the surface represented by this equation for {z > 0, r > 0}.

surface = z /. Assuming[{z > 0, r > 0}, Solve[%, z]][[2]];
Plot3D[surface, {p, 0, 1}, {r, 0, 2}, AxesLabel -> {p, r, s}]
(* Sqrt[4 - 8*p^2 + 4*p^4 - r^2]/2 *)

The differential area of this surface is given by

darea = Sqrt[1 + D[surface, p]^2 + D[surface, r]^2] // Simplify;

and the total area by

area = Integrate[p^2 r darea, {p, 0, 1}, {r, 0, 2 (1 - p^2)}] // N

(Its value is 0.401377.) The resulting probability function is simply

Piecewise[{{p^2 r darea/area, 0 <= r <= 2 (1 - p^2)}}, 0]

Finally, we transform back to the original coordinates to obtain

{px, py, pz, x, y, z} = 
  {p Sin[θ] Cos[ϕ], p Sin[θ] Sin[ϕ], p Cos[θ], r Cos[ξ], r Sin[ξ], 
   p Cos[θ], r Cos[ξ], r Sin[ξ], Sqrt[4 - 8*p^2 + 4*p^4 - r^2]/2}

Here, θ is drawn from

RandomVariate[ProbabilityDistribution[Sin[θ0], {θ0, 0, Pi/2}]]

ϕ and ξ are drawn from

RandomVariate[UniformDistribution[{0, Pi/2}]]

and we would like to draw {p, r} from

RandomVariate[ProbabilityDistribution[
  Piecewise[{{p^2 r darea/area, 0 <= r <= 2 (1 - p^2)}}, 0], {p, 0, 1}, {r, 0, 2}]]

Unfortunately, ProbabilityDistribution cannot handle this expression, as noted in 69150. Presumably, one of the approaches given in 2635 can.