Liouville equation solution $\rho (q,p,t)$ for 1D-oscillator I'd like to consider the 1D-oscillator phase space probability density evolution problem with the ordinary 
Hamiltonian 
$$H = \frac{p^2}{2m} + \frac{kq^2}{2}.$$
Then the Liouville theorem is just 
$$\frac{{\partial \rho }}{{\partial t}} = { {\frac{{\partial \rho }}{{\partial {q}}}\frac{{\partial H}}{{\partial {p}}} - \frac{{\partial \rho }}{{\partial {p}}}\frac{{\partial H}}{{\partial {q}}}}} $$
for the case of one particle.
Assume that the initial distribution of $q$ and $p$ for a particle is normal
$\rho (q,p,t = 0) = \frac{1}{{2\pi }}{e^{ - {q^2} - {p^2}/2}}$
Now we can find the Hamiltonian partial derivatives
$\frac{\partial H}{\partial p} = p/m$ and $ \frac{\partial H}{\partial q} = -kq$.
After this step I'm confused how we can find the evolution of the probability?
How are $p$, $q$ related with distribution $\rho (q,p,t = 0)$? 
What is the $\frac{{\partial \rho }}{\partial {q}}$ and $\frac{{\partial \rho }}{\partial {p}}$ when
we even don't know $\rho (q,p,t)$ in arbitrary moment of time $t$?
Maybe my undestanding is wrong, but I just want to know $\rho (q,p,t)$.
 A: Consider a set of $N$ harmonic oscillators. $\rho(q,p,t=0)$ gives the probability to find one of these oscillators in the state $(q,p)$ at time $t$. Using the equations of motion (Hamilton equations in your case)
    $${dq\over dt}={p\over m},\hskip 1cm {dp\over dt}=-kq$$
you can tell that if there was an oscillator in the state $(q_0,p_0)$ then this
oscillator is in the state
    $$\left\{\eqalign{
q(t)&=q_0\cos\omega t+{p_0\over\omega}\sin\omega t\cr
p(t)&=-\omega x_0\sin\omega t+p_0\cos\omega t\cr
}\right.$$
at time $t$. Equivalently, you can tell that, if there is an oscillator in the state $(q,p)$ at time $t$, then it was in the state
    $$\left\{\eqalign{
    q_0&=q\cos\omega t-{p\over\omega}\sin\omega t\cr
    p_0&=\omega q\sin\omega t+p\cos\omega t\cr
    }\right.$$
at time 0. Therefore, the probability to find an oscillator in the state $(q,p)$ at time $t$ is
   $$\rho(q(t),p(t),t)=\rho(q(0),p(0),0)$$
which reads
   $$\rho(q,p,t)=\rho\Big(q_0=q\cos\omega t-{p\over\omega}\sin\omega t,
   p_0=\omega q\sin\omega t+p\cos\omega t,0\Big)$$
This is actually the solution of Liouville equation
   $${d\rho\over dt}={\partial\rho\over\partial t}
+\dot q{\partial\rho\over\partial q}
+\dot p{\partial\rho\over\partial p}=0$$
and this way of solving it is known as the method of characteristics.
A: According to Liouville's theorem,  
$$\frac{d\rho}{dt}=\frac{\partial \rho}{\partial t}+\{\rho,H\}=0$$ where {,} stand for the Poisson brackets. Now, in equilibrium statistical mechanics, we are concerned about the behaviour of systems in equilibrium when $\rho$ has no explicit time dependence, i.e., when$$\frac{\partial \rho}{\partial t}=0 .$$ So, to make these two equations valid simultaneously (or in other words, for the Liouville's theorem to hold when the system is in equilibrium), we need to have $\{\rho,H\}=0$ . This is possible if $\rho(q,p)$ is equal to some constant or more generally if $\rho=\rho[H(q,p)]$ , i.e., if $\rho$ has its dependence on q and p only through $H$.
To see this, you may read Statistical Mechanics by R. K. Pathria and Paule D. Beale (Second Chapter of Third Edition). 
Now, the form of $\rho$ depends on the type of the ensemble concerned; for example, $\rho$=constant for a microcanonical ensemble because all the microstates are equally probable here.  In a canonical ensemble, $\rho$ is proportional to $e^{-H(X)/k_BT}$ where $X=\{q_i,p_i\}$; $i=\{1,2,...,N\}$ , for N particles in one dimension.  
For your problem, I do not see how the initial distribution that you have mentioned is relevant. However, if we assume that there is an isolated system of energy U and volume V comprising of N non-interacting 1D linear harmonic oscillators, then it forms a microcanonical ensemble where the suitable form for $\rho(q,p)$ is as follows :
$$\rho(q,p)=\frac{\delta(H_N(X)-U)}{\Sigma}$$ where $\Sigma=\displaystyle{\int}_{\Gamma}\frac{d^{2N}X}{h^N}\delta(H_N(X)-U)$ is the microcanonical partition function. Here, $H_N(X)=\displaystyle{\sum_{i=1}^N}\frac{p_i^2}{2m}+\frac{kq_i^2}{2}$ . So,
$$\Sigma=\displaystyle{\int}_{\Gamma}\frac{d^{2N}X}{h^N}\delta(H_N(X)-U)= \displaystyle{\int_{U{\le}H_N(X){\le}U+\Delta}}\frac{d^{2N}X}{h^N}$$
$$=\displaystyle{\int_{U{\le}\sum_{i=1}^N\frac{p_i^2}{2m}+\frac{kq_i^2}{2}{\le}U+\Delta}}\frac{d^{2N}X}{h^N}  
={(2m)}^{N/2}{\left(\frac{2}{k}\right)}^{N/2}\displaystyle{\int_{U{\le}\sum_{i=1}^N(p_i^2+q_i^2){\le}U+\Delta}}\frac{{(dpdq)}^N}{h^N}$$
$$=\frac{{(2m)}^{N/2}}{h^N}{\left(\frac{2}{k}\right)}^{N/2}\times\,Volume\,of\, a \,shell\, between\, the\, 2N-dimensional\, hyperspheres\, of\, radii\, \sqrt{U}\, and\, \sqrt{U+\Delta}\,\,[where\, \Delta<<U]$$
$$=\frac{{(2m)}^{N/2}}{h^N}{\left(\frac{2}{k}\right)}^{N/2}\times\,surface\,of\,a\,2N-dimensional\,sphere\,of\,radius\, \sqrt{U}\,\times\,thickness\,of\,the\,shell$$
$$=\frac{{(2m)}^{N/2}}{h^N}{\left(\frac{2}{k}\right)}^{N/2}\frac{(2N){\pi}^{N}({\sqrt{U}}^{2N-1})}{\Gamma(N+1)}\times\frac{\Delta}{2\sqrt{U}}$$
which you can further simplify.
This is for distinguishable oscillators. If they are indistinguishable, $\Sigma$ should be divided by $N!$ . Putting the expression for $\Sigma$ in that of $\rho$, you will get the desired expression.
