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'$

  • $\begingroup$ ...This sounds like a hard research problem. Even in your 1D example, if you don't have $\frac1{2m}=\frac{m\omega^2}2$ you're trying to sample uniformly from the boundary of an ellipse, and that's not easy. $\endgroup$ – user484 Dec 14 '14 at 17:50
  • $\begingroup$ Do you have a particular function in mind for V(x,y,z)? If not, can you state whether E = E' is a continuously differentiable closed surface? $\endgroup$ – bbgodfrey Dec 14 '14 at 19:13
  • $\begingroup$ The potential I have in mind is $V(x,y,z) = c \sqrt{\frac{x^2+y^2}{4 a^2} + \frac{z^2}{a^2}}$ where $a$ and $c$ are constants. When I want to generate the random numbers, $a, c, m$ and $E$ will be some numerical value. $\endgroup$ – user29165 Dec 14 '14 at 19:48
  • $\begingroup$ @user29165, do you really mean to include the Sqrt, which causes V not to be analytical at the origin? I would have expected something like V [x, y, z] = c ((x^2 + y^2)/(4 a^2) + z^2/a^2), analogous to your 1D example. $\endgroup$ – bbgodfrey Dec 14 '14 at 20:18
  • $\begingroup$ In my case, I do have the sqrt in the potential, I was just listing the 1d harmonic oscillator as an example as it would be familiar to most people. $\endgroup$ – user29165 Dec 14 '14 at 20:24

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}]

phase space

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.

|improve this answer|||||

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

Phase surface

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

  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.

|improve this answer|||||

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.