Any methods of solving this system of ODE's? I try to solve this system of ODE's: 
$$ \frac{dQ_1 (t)}{dt} =  - a \sin (\omega t) Q_2(t) + b \cos(\omega t) Q_3(t) $$
$$ \frac{dQ_2 (t)}{dt} = - a \sin (\omega t) Q_1 (t) - c Q_3(t) $$ 
$$ \frac{dQ_3 (t)}{dt} = c Q_2(t) - b \cos(\omega t) Q_1(t) $$
with constant $a, b, c$.
In terms of vectors RHS is just cross production of $P = \{ c, b \cos(\omega t), a \sin(\omega t)\}$ and $Q$ vectors.
Is there any analytical method to solve this system, at least approximately?.. For example, if $a, b, c \ll \omega$.
 A: This is a linear time-varying system of the form
$$\frac{dQ}{dt}=A(t)Q(t)$$
with solution
$$Q(t)=\Phi(t,0)Q(0)$$
where $\Phi(t,t_0)$ is the transition matrix that is calculated by the (convergent) series
$$\Phi(t,t_0)=\mathbb{I}+\int_{t_0}^t{A(s_1)ds_1}+\int_{t_0}^t{A(s_1)\int_{t_0}^{s_1}{A(s_2)ds_2}ds_1}+\cdots$$
Note that in this case 
$$A(t)=A_1+A_2\sin(\omega t)+A_3\cos(\omega t)$$ 
For times $t$ that are multiples of the period $2\pi/\omega$ all the integrals with $A_2, A_3$ terms are zero and the series above takes the form
$$\Phi\left(\frac{2k\pi}{\omega},0\right)=\mathbb{I}+A_1\frac{2k\pi}{\omega}+\frac{1}{2!}A_1^2\left(\frac{2k\pi}{\omega}\right)^2+\cdots=e^{\frac{2k\pi}{\omega}A_1}\:,\quad\forall k=0,1,\cdots $$
Thus,
$$Q\left(\frac{2k\pi}{\omega}\right)=e^{\frac{2k\pi}{\omega}A_1}Q(0)=\left[\matrix{Q_1(0) \\  m\cos(\frac{2ck\pi}{\omega}+\phi)\\ m\sin(\frac{2ck\pi}{\omega}+\phi)}\right]\:,\quad\forall k=0,1,\cdots $$
with $m,\phi$ depending on the initial values $Q_2(0),Q_3(0)$.
For the other points $t\neq 2k\pi/\omega$ we have
$$\Phi(t,0)=\Phi\left(t,\frac{2\pi}{\omega}\bigg\lfloor\frac{\omega t}{2\pi}\bigg\rfloor\right)\Phi\left(\frac{2\pi}{\omega}\bigg\lfloor\frac{\omega t}{2\pi}\bigg\rfloor,0\right)=\Phi\left(t,\frac{2\pi}{\omega}\bigg\lfloor\frac{\omega t}{2\pi}\bigg\rfloor\right)e^{\frac{2\pi}{\omega}\big\lfloor\frac{\omega t}{2\pi}\big\rfloor A_1}$$
So for the intermediate points we have to calculate the transition matrix within a time interval less than the period $2\pi/\omega$. A closed form solution may be tedious to derive but we can discuss the qualitative behavior of the solutions.
We can prove for example that $\|Q(t)\|$ is bounded. Using the fact that $\|A(t)\|<\delta$ for all $t\geq0$ for some $\delta>0$  we can prove that
$$\|\Phi(t_2,t_1)\|\leq e^{\delta(t_2-t_1)} \:,\quad \forall t_2\geq t_1\geq 0$$
Hence
$$\left\|\Phi\bigg(t,\frac{2\pi}{\omega}\bigg\lfloor\frac{\omega t}{2\pi}\bigg\rfloor\bigg)\right\|\leq \exp\left(\frac{2\pi\delta}{\omega}\bigg(\frac{\omega t}{2\pi}-\bigg\lfloor\frac{\omega t}{2\pi}\bigg\rfloor\bigg)\right)\leq \exp\left(\frac{2\pi\delta}{\omega}\right)$$
and therefore
$$\|Q(t)\|\leq \exp\left(\frac{2\pi\delta}{\omega}\right)\sqrt{Q_1^2(0)+m^2}$$
Also $\|dQ/dt\|$ is bounded by a constant that does not increase when $\omega$ increases. Actually
$$\bigg\|\frac{dQ}{dt}\bigg\|\leq \|A(t)\|\|Q(t)\|\leq \delta \exp\left(\frac{2\pi\delta}{\omega}\right)\sqrt{Q_1^2(0)+m^2}$$
If we assume now that $\omega$ increases then the points $2k\pi/\omega$ get closer while the allowed variations determined by $\|dQ/dt\|$ remain bounded. This explains why when $\omega\rightarrow\infty$ 
$$ Q(t)\rightarrow e^{A_1t}Q(0)=\left[\matrix{Q_1(0) \\  m\cos(ct+\phi)\\ m\sin(ct+\phi)}\right]$$
One can see in the Figures below  that as $\omega$ increases the described behavior occurs.



A: Since $\dot{Q}=P\wedge Q$, the norm of the vector $Q$ is conserved. Indeed, $$\frac{d}{dt}\left( Q\cdot Q\right)=Q\cdot(P\wedge Q)+(P\wedge Q)\cdot Q=0.$$
This allows to decrease the number of degrees of freedom by $1$. For example, we can parameterize $Q_{1,2,3}$ with spherical angles:
\begin{align*}
Q_1&=R\cos\theta,\\
Q_2&=R\cos\phi\sin\theta,\\
Q_3&=R\sin\phi \sin\theta,
\end{align*}
where $R$ is constant. Now the system reduces to two $1$st order equations for $\theta,\phi$:
\begin{align*}
\dot{\theta}&=a\sin\omega t\cos\phi-b\cos\omega t\sin \phi,\\
\dot{\phi}&=c-\left(a\sin\omega t\sin\phi+b\cos\omega t\cos \phi\right)\cot\theta.
\end{align*}
Further thoughts: One way to rewrite the last system that may turn out to be useful is as follows. Introduce two functions
\begin{align*}
A(t)&=\sqrt{a^2\sin^2\omega t+b^2\cos^2\omega t},\\
\phi_0(t)&=\arctan\left(\frac{a}{b}\tan\omega t\right),
\end{align*}
and use instead of $\phi$ a translated function $\varphi:=\phi-\phi_0(t)$. Then we can write
\begin{align*}
\dot{\theta}&=-A(t)\sin\varphi,\\
\dot{\varphi}&=-A(t)\cos\varphi\cot\theta+c-\dot{\phi}_0(t).
\end{align*}
Furthermore, if we change the independent variable
$$t\mapsto s(t)=-\int^tA(t)\,dt,$$
we end up with the system
\begin{align*}
\frac{d\theta}{ds}&=\sin\varphi,\\
\frac{d\varphi}{ds}&=\cos\varphi\cot\theta+C(s),
\end{align*}
where $C(s(t))=\displaystyle\frac{\dot{\phi}_0(t)-c}{A(t)}$. This is the simplest form I have found so far.
Special case $a=b$: Here the above expressions simplify to
$$A(t)=a,\qquad \phi_0(t)=\omega t, \qquad s(t)=-a t,\qquad C(s)=\frac{\omega-c}{a}.$$
The equations of motion for $\theta$ and $\varphi$ are
\begin{align*}
\dot{\theta}&=-a\sin\varphi,\\
\dot{\varphi}&=-a\cos\varphi\cot\theta+\left(c-\omega\right).
\end{align*}
Now introducing
\begin{align*}
M_1&=\cos\theta,\\
M_2&=\cos\varphi\sin\theta,\\
M_3&=\sin\varphi \sin\theta,
\end{align*}
these equations can be written as
$$\dot{M}=B\wedge M,\qquad B=\left(\begin{array}{c}  c-\omega \\  a \\ 0 \end{array}\right).$$
This problem can be easily solved. In physics language, it describes Larmor precession of a magnetic moment in the uniform magnetic field $B$. The idea is to rotate the coordinate system so that $B$ aligns along one of the axes, then the moment projection on this axis is conserved, and in the plane orthogonal to it we will have a uniform rotation. 
Here is the solution:
\begin{align*}
M_1&=a\left(\alpha\cos\omega_0t+\beta\sin\omega_0t\right)+\left(c-\omega\right)\gamma,\\
M_2&=-\left(c-\omega\right)\left(\alpha\cos\omega_0t+\beta\sin\omega_0t\right)+a\gamma,\\
M_3&=\omega_0\left(\beta\cos\omega_0t-\alpha\sin\omega_0t\right),
\end{align*}
where $\omega_0=\sqrt{a^2+\left(c-\omega\right)^2}$ and $\alpha,\beta,\gamma$ are arbitrary constants satisfying $$\alpha^2+\beta^2+\gamma^2=\omega_0^{-2}.$$
From this it is not difficult to get the solution for $Q_{1,2,3}$:
\begin{align*}
Q_1&=R\Bigl[ a\left(\alpha\cos\omega_0t+\beta\sin\omega_0t\right)+\left(c-\omega\right)\gamma\Bigr],\\
Q_2&=R\cos\omega t\Bigl[a\gamma-\left(c-\omega\right)\left(\alpha\cos\omega_0t+\beta\sin\omega_0t\right)\Bigr]-R\omega_0\sin\omega t\Bigl[\beta\cos\omega_0t-\alpha\sin\omega_0t\Bigr],\\
Q_3&=R\sin\omega t\Bigl[a\gamma-\left(c-\omega\right)\left(\alpha\cos\omega_0t+\beta\sin\omega_0t\right)\Bigr]+R\omega_0\cos\omega t\Bigl[\beta\cos\omega_0t-\alpha\sin\omega_0t\Bigr].
\end{align*}
