Lagrange's planetary equations

Statement

Let $\mu > 0$ and let a body move by $$r'' = -\frac{\mu\, r}{\norm{r}^3} + \nabla_r \mathcal{R}(r, t),$$ where the disturbing function $\mathcal{R}$ is a smooth function and $\nabla_r$ is its gradient in $r$. Assume that the osculating orbit stays elliptic with $0 < e < 1$ and $0 < i < \pi$. Let $a, e, i, \Omega, \omega$ and the mean anomaly $M$ be its orbital elements, and let $n = \sqrt{\mu / a^3}$ be the mean motion. Write $\mathcal{R}$ as a function of these elements and $t$, with partial derivatives in $a$ taken at fixed $M$. Then the elements satisfy $$\frac{da}{dt} = \frac{2}{na} \frac{\partial \mathcal{R}}{\partial M},$$ $$\frac{de}{dt} = \frac{1 - e^2}{n a^2 e} \frac{\partial \mathcal{R}}{\partial M} - \frac{\sqrt{1 - e^2}}{n a^2 e} \frac{\partial \mathcal{R}}{\partial \omega},$$ $$\frac{di}{dt} = \frac{\cos i}{n a^2 \sqrt{1 - e^2} \sin i} \frac{\partial \mathcal{R}}{\partial \omega} - \frac{1}{n a^2 \sqrt{1 - e^2} \sin i} \frac{\partial \mathcal{R}}{\partial \Omega},$$ $$\frac{d\Omega}{dt} = \frac{1}{n a^2 \sqrt{1 - e^2} \sin i} \frac{\partial \mathcal{R}}{\partial i},$$ $$\frac{d\omega}{dt} = \frac{\sqrt{1 - e^2}}{n a^2 e} \frac{\partial \mathcal{R}}{\partial e} - \frac{\cos i}{n a^2 \sqrt{1 - e^2} \sin i} \frac{\partial \mathcal{R}}{\partial i},$$ $$\frac{dM}{dt} = n - \frac{1 - e^2}{n a^2 e} \frac{\partial \mathcal{R}}{\partial e} - \frac{2}{na} \frac{\partial \mathcal{R}}{\partial a}.$$