Double pendulum probability distribution Double Pendulum has a very beautiful stochastic trajectory. Is there any way to calculate the distribution of probability of finding the end of pendulum at each point? 
Link to formulations.

 A: What you could try is to first write the system in Hamiltonian form.  Then quote https://en.wikipedia.org/wiki/Liouville%27s_theorem_(Hamiltonian) which says that the distribution on phase space is uniform.  You could assume ergodicity, that is, the solution path is dense in the energy submanifold $H(q,p) = E$.  Then project this onto your 2-D space.
I think it would be a lot of work, but definitely doable.
So making a start, the potential energy is
$$ V = - m_1 g L_1 \cos(\theta_1) - m_2 g L_2 \cos(\theta_2) ,$$
and the Kinetic Energy is
$$ T = \tfrac12 [\dot \theta_1, \dot \theta_2] A [\dot \theta_1, \dot \theta_2]^T = \tfrac12 [\eta_1, \eta_2] A^{-1} [\eta_1,  \eta_2]^T $$
where
$$ A = \begin{bmatrix} m_1 L_1^2 + m_2 L_1^2 & m_2 L_1 L_2 \cos(\theta_1-\theta_2) \\ m_2 L_1 L_2 \cos(\theta_1-\theta_2) & m_2 L_2^2 \end{bmatrix} $$
and $\eta_1$ and $\eta_2$ are the generalized momentums.
See https://physics.stackexchange.com/questions/142238/non-integrability-of-the-2d-double-pendulum
The Hamiltonian $H = V + T$ is conserved, so calculate its value using initial conditions.  For a given $[\theta_1,\theta_2]$, the probability density will be proportional to the volume of the ellipsoid
$$ \{[\eta_1,\eta_2] \in \mathbb R^2 : \tfrac12 [\eta_1, \eta_2] A^{-1} [\eta_1,  \eta_2]^T \le H - V\}$$
And this volume will be proportional to
$$ \begin{cases}\det(A) (H-V)^2 & \text{if $H > V$}\\ 0 &\text{if $H \le V$}\end{cases} $$
Remember both $V$ and $A$ are functions of $\theta_1$ and $\theta_2$.  Note that for any place where the pendulum can go, there will generally be two values of $\theta_1$ and $\theta_2$ for that point, where for one point the value of $\theta_1 - \theta_2$ will be exactly the negative of the value for the other position.
Remember, this will give you the density in $[\theta_1,\theta_2]$ space (the green picture on the left side of the web page you provided), so you will have to divide this by the absolute value of determinant of the Jacobian of the map that takes $[\theta_1,\theta_2]$ to $[x,y]$ (that is, $L_2L_2|\sin(\theta_1-\theta_2)|$).
Now my calculations might be all wrong because maybe there is another constant of the motion (which means I lose the ergodicity property).  The web page https://physics.stackexchange.com/questions/142238/non-integrability-of-the-2d-double-pendulum suggests there isn't another constant of motion.  But in the green picture in http://www.myphysicslab.com/pendulum/double-pendulum/double-pendulum-en.html, I can plainly see that for any choice of $[\theta_1, \theta_2]$ that the velocity has a particular direction (or rather a choice of two directions).  Conservation of the Hamiltonian (energy) would only restrict the magnitude of the velocity going through each point!
Added later: I asked someone more knowledgeable than I about ergodicity.  There are even situations where the set of points reached in phase space has positive measure on a level set of the energy function, but the double pendulum is still not ergodic.  Apparently many books have been written on the subject.  So my final answer is "I don't know."  I'll not delete this answer, because I think the discussion is worthwhile.  But I don't expect to receive any upvotes.
