Averaging of the Van der Pol equation $\ddot V+V=\varepsilon (1-V^2)\dot V$. In this link we try to approximate the solution of van der Pol equation using averaging. It goes as follow :
Let $$\ddot V+V=\varepsilon (1-V^2)\dot V.$$
We transform it in $$\begin{cases}\dot V=I\\ \dot I=-V+\varepsilon (1-V^2)I.\end{cases}$$
We do the substitution $$\begin{pmatrix}X\\Y\end{pmatrix}=\begin{pmatrix}\cos t&-\sin t\\ \sin t&\cos t\end{pmatrix}\begin{pmatrix}V\\I\end{pmatrix}.$$
By setting $$f(X,Y,t,\varepsilon )=\begin{pmatrix}-\Big(1-(X\cos t+Y\sin t)^2\Big)(-X\sin t+Y\cos t)\sin t\\ \Big(1-(X\cos t+Y\sin t)^2\Big)(-X\sin t+Y\cos t)\cos t\end{pmatrix},$$
one can compute $$\bar f(X,Y):=\frac{1}{2\pi}\int_0^2\pi f(X,Y,t,0)dt=\frac{1}{8}(4-(X^2+Y^2))\begin{pmatrix}X\\Y\end{pmatrix}.$$
We have that $$\bar f(x,y)=0\iff V^2+I^2=X^2+Y^2=4,$$
which correspon to the limit cycl of the van der Pol oscillator. Set $\tau=X^2+Y^2$, this is an invariant of the flow of the linear equation.
Q1) What does it mean ?
Dropping all information on the phase, one can reduce the averaged equation to $$\dot \tau=\left(1-\frac{\tau}{4}\right)\tau.\tag{E}$$
Q2) How did they get this equation ? I don't really understand the trick since in one way we have a vectorial equation, whereas (E) is just a simple ODE. 
 A: @Q2) first-order system
This is the transformation of an higher order equation into a first order system.
General approach
The unperturbed problem, that is, for $ε=0$, is a simple harmonic oscillator $\ddot V+V=0$ with solution $V=R\cos(t+\phi)$. This leads to the idea to also set in the perturbed system $$V(t)=R(t)\cos(t+\phi(t))\tag1$$ where the functions $R(t),\phi(t)$ are nearly constant, or slowly changing relative to the period $2\pi$.
There are now two functions parametrizing one scalar solution, so one has one dependency free to select. Fix it so that $$\dot V(t)=-R(t)\sin(t+\phi(t))\tag2$$ so that $\dot V^2+V^2=R^2$. In the difference to the derivative formula this means that $$\dot R(t)\cos(t+\phi(t))-R(t)\sin(t+\phi(t))\dot\phi(t)=0\tag3$$ and insert into the equation
\begin{align}
\ddot V+V&=-\dot R(t)\sin(t+\phi(t))-R(t)\cos(t+\phi(t))\dot\phi(t)\\
=ε(1−V^2)\dot V&=-ε\Bigl(1-R(t)^2\cos^2(t+\phi(t))\Bigr)R(t)\sin(t+\phi(t)).\tag4
\end{align}
Now use the last two identities to isolate $\dot R$ and $\dot \phi$,
\begin{alignat}1
\dot R&=εR\Bigl(1-R^2\cos^2(t+\phi)\Bigr)\sin^2(t+\phi)
&=\fracε4R\Bigl(2-2\cos(2t+2\phi)-R^2\sin^2(2t+2\phi)\Bigr),\\
\dot\phi&=ε\Bigl(1-R^2\cos^2(t+\phi)\Bigr)\sin(t+\phi)\cos(t+\phi)
&=\fracε4\Bigl(2-R^2-R^2\cos(2t+2\phi)\Bigr)\sin(2t+2\phi).
\end{alignat}
@Q1) averaging the perturbation over a cycle
Now for the averaging process set the slowly varying functions on the right as constant and integrate over an interval of $2\pi$, followed by a division by $2\pi$. This then gives
\begin{alignat}1
\dot{\bar R}&=\fracε8R\Bigl(4-\bar R^2\Bigr),\\
\dot{\bar\phi}&=0.
\end{alignat}
so that with $\tau(εt)=\bar R(t)^2$ the equation $\tau'=\tau(1-\frac14\tau)$ follows.
