Any methods of solving this system of ODE's?

Question

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$.

Answer 1

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*}

Answer 2

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.