The interplay of Hamiltonian and Lagrangian theory is based on the following general identities, where $L$ is the Lagrangian function of the system,
$$\dot{q}^k = \frac{\partial H}{\partial p_k}\:,\qquad(1)$$
$$\frac{\partial L}{\partial q^k} = -\frac{\partial H}{\partial q^k}\:.\qquad(2)$$
Above, the RH sides are functions of $t,q,p$ whereas the LH sides are functions of $t,q,\dot{q}$ and the two types of coordinates are related by means of the bijective smooth (with smooth inverse) relation,
$$t=t\:,\quad q^k=q^k\:,\quad p_k = \frac{\partial L(t,q,\dot{q})}{\partial \dot{q}^k}\:.\qquad(3)$$
Finally, the Hamiltonian function is defined as follows
$$H(t,q,p) = \sum_{k}\left.\frac{\partial L}{\partial \dot{q}^k}\right|_{(t,q,\dot{q}(t,q,p))}\dot{q}(t,q,p) - L(t,q,\dot{q}(t,q,p))\:.$$
Suppose that the following E-L hold,
$$\frac{d}{dq}\left(\frac{\partial L}{\partial \dot{q}^k}\right) - \frac{\partial L}{\partial q^k} = Q_k(t,q,\dot{q})\:, \quad \frac{d q^k}{dt}= \dot{q}^k \quad (4)\:.$$
The functions $Q_k$ take the (e.g. dissipative) forces into account, which cannot be included in the Lagrangian. For a system of $N$ points of matter with positions $\vec{x}_i$, if the degrees of freedom of the system are described by coordinates $q^1,\ldots,q^n$ such that $\vec{x}_i= \vec{x}_i(t,q^1,\ldots,q^n)$, one has:
$$Q_k = \sum_{i=1}^N \frac{\partial \vec{x}_i}{\partial q^k} \cdot \vec{f_i}$$
$\vec{f}_i$ being the total force, not described in the lagrangian, acting on the $i$th point.
One easily proves that, in view of the general identities (1) and (2), a curve $t \mapsto (t, q(t), \dot{q}(t))$ satisfies the EL equations (4), if and only if the corresponding curve $t \mapsto (t, q(t), p(t))$ (constructed out of the previous one via (3)), verifies the following equations:
$$\frac{dq^k}{dt} = \frac{\partial H}{\partial p_k}\:, \quad \frac{dp_k}{dt} = -\frac{\partial H}{\partial q^k} + Q_k\:.\quad(5)$$
In the absence of the terms $Q_k$, these are the standard Hamilton equation.
If $Q_k\equiv 0$, even if $H$ explicitely depend on time, the solutions of Hamilton equations preserve, in time, the canonical volume:
$$dq^1 \wedge \cdots \wedge dq^n \wedge dp_1 \wedge \cdots \wedge d p_n\:.$$
In the presence of dissipative forces which cannot be included in the Lagrangian, the term $Q_k$ show up and the volume above generally fails to be preserved. However this is not the whole story.
Let us consider the damped harmonic oscillator. In the absence of dissipative force, the Lagrangian reads
$$L(x, \dot{x}) = \frac{m}{2} \dot{x}^2 -\frac{k}{2} x^2\:.$$
The dissipative force $-\gamma \dot{x}$ takes place in the EL equations due to the presence of the term:
$$Q = - \gamma \dot{x}\:.$$
In this juncture, passing to the Hamiltonian formulation, the canonical volume is not preserved along the solutions of the equation of motion.
However, sticking to the damped oscillator, there is a way to include the dissipative force in the Lagrangian function. As a matter of fact, this new Lagrangian produces the correct equation of motion of a damped oscillator
$$L(t,q,\dot{q}) = e^{\gamma t/m}\left(\frac{m}{2} \dot{x}^2 -\frac{k}{2} x^2\right)\:.$$
In this case the Hamiltonian function turns out to be
$$H(t,q,p) = e^{-\gamma t/m}\frac{p^2}{2m} + e^{\gamma t/m}\frac{k}{2}x^2\:.$$
As the general theory proves, the canonical volume is preserved by the solutions of Hamilton equations referred to that Hamiltonian function, regardless the fact that the system is dissipative.
It is important to notice that, with the second Lagrangian, $p$ ceases to be the standard momentum $mv$ differently from the first case.
Here's something I believe is a simple proof. Unfortunately it uses a little bit of cohomology.
Consider the canonical 2-form in extended phase space $T^*M \times \mathbb{R}$
$$\omega = \sum_{i=1}^N dq_i \wedge dp_i - dH(\vec{q},\vec{p},t) \wedge dt ,$$
where $N = dim(M)$. A function $f: M \to M$ is said to be a canonical transformation iff $f^* \omega = \omega$. Thus, for any canonical transformation,
$$\sum_{i=1}^N dq_i \wedge dp_i - dH(\vec{q},\vec{p},t) \wedge dt = \sum_{i=1}^N dQ_i \wedge dP_i - dK(\vec{Q},\vec{P},T) \wedge dT ,$$
where we defined $(\vec{Q},\vec{P},T) = f^*(\vec{q},\vec{p},t)$. This means that
$$\sum_{i=1}^N q_i dp_i - H(\vec{q},\vec{p},t) dt - \left( \sum_{i=1}^N Q_i dP_i - dK(\vec{Q},\vec{P},T) dT \right) = dG ,$$
i.e, that the tautological forms associated to the canonical 2-forms are elements of the same de Rham cohomology class. By the homotopic invariance of the de Rham cohomology,
$$\oint_\gamma \left( \sum_{i=1}^N q_i dp_i - H(\vec{q},\vec{p},t) dt \right) = \oint_{\Gamma} \left( \sum_{i=1}^N Q_i dP_i - dK(\vec{Q},\vec{P},T) dT \right) ,$$
where $\Gamma$ is the image of the curve $\gamma$ by $f$.
As movement can be considered to be a canonical transformation, we just proved that the hamiltonian evolution in extended phase space needs to relate only homotopic curves. A corollary is that if you start with a compact simply-connected set in phase space and track its evolution in time, the curve that bounds this set can only evolve to others homotopic to it. As a simply connected set boundary is never homotopic to a non-simply connected one, we just showed that hamiltonian evolution takes simply-connected sets to simply-coonected sets.
I'm really sorry if this is too technical. I'm too ignorant to offer a proof with simpler arguments (which points out how limited is my knowledge on the subject).
EDIT: Small correction on dimension.
Best Answer
To obtain the result $\frac{\text d \rho }{\text d t}=0$ you need two facts: the first is that the hamiltonian flow preserves the volume of phase space. The second fact is the conservation of probability, that is, the probability that the system is found in a volume $U$ at time $t=0$ equals the probability of finding it within $\Phi _t U$ at time $t$, where $\Phi _t$ denotes the hamiltonian flow. This is a direct consequence of the deterministic nature of classical mechanics: the two propositions “$(p(0),q(0))\in U$” and “$(p(t),q(t))\in \Phi _t U$” are equivalent.
Using conservation of probability, for an arbitrary volume $U$ we can write an equation: $$\int _U \rho(p,q,0) \text d p \text d q=\int _{\Phi _t U} \rho(p,q,t)\text dp \text d q .$$ By Jacobi's theorem: $$\int _{\Phi_t U} \rho (p,q,t)\text d p \text d q=\int _U\rho (\Phi _t (p,q),t)\text J_{\Phi _t}d p \text d q.$$ The Jacobian $J_{\Phi _t}=1$, because the flow preserves volumes. It follows that: $$\int _U \rho (p,q,0)\text d p \text d q =\int _U \rho (\Phi _t (p,q),t)\text d p \text d q,$$ and, since the volume $U$ was arbitrary, $\rho (p,q,0)=\rho (\Phi _t (p,q),t)$, or $\text d\rho /\text d t=0$.