(In)stability of a two-dimensional dynamical system

Question

Consider the following system of coupled differential equations \begin{eqnarray*} \dot{x}_1(t) & = & -x_1(t) - \cos(\omega t)x_1(t) + \cos(\omega t)x_2(t), \ x_1(0)\in\mathbb{R},\\ \dot{x}_2(t) & = & -\gamma x_2(t) - \cos(\omega t)x_2(t) + \cos(\omega t)x_1(t), \ x_2(0)\in\mathbb{R}, \end{eqnarray*} where $\omega$ and $\gamma$ are positive real constants. Observe that $\bar{x}=(\bar{x}_1,\ \bar{x}_2)=(0,0)$ is an equilibrium of the above system.

It is almost trivial to see that if $\gamma=1$ then $\bar x$ is attractive. Indeed, in this case, we have that $x(t)=[x_1(t), x_2(t)]^\top$ can be explicitly computed as $$ x(t) = \exp\left(\begin{bmatrix}-t &0\\ 0 & -t\end{bmatrix} + \frac{1}{\omega}\sin(\omega t)\begin{bmatrix}-1 &1\\ 1 & -1\end{bmatrix}\right)x(0), $$ so that $x(t)\to 0$ as $t\to \infty$.

However, in case $\gamma\ne 1$ proving the attractiveness of the origin is not obvious (and perhaps not even true!).

In particular, numerical simulations seem to suggest that for $\gamma$ and $\omega$ sufficiently small (e.g. $\gamma=0.001$ and $\omega=10$) the equilibrium $\bar{x}$ is not attractive.

I've struggled a lot to find a way of formally proving this, with no luck. So I decided to post the problem here hoping that some of you will provide some useful suggestions or tips. Thank you!

---

I post here the Mathematica code that I've used in my simulations:

(* nominal values for simulation *)
values = {gamma -> 0.001, w -> 10};

equations = {
   {x1'[t], x2'[t]} == {-x1[t] - Cos[w*t]*x1[t] + Cos[w*t]*x2[t], -gamma*x2[t] - Cos[w*t]*x2[t] + Cos[w*t]*x1[t]},
   {x1[0], x2[0]} == {0.1, 0.1}};

{x1t, x2t} = NDSolveValue[equations /. values, {x1[t], x2[t]}, {t, 0, 1000}];

Plot[x1t, {t, 0, 1000}, PlotRange -> {-0.2, 0.2}]
Plot[x2t, {t, 0, 1000}, PlotRange -> {-0.2, 0.2}]

---

Further remarks. Since the system is periodic, one could exploit Floquet theory to express the transition matrix of the system in the form $P(t)e^{Rt}$ where $P(t)$ is a periodic function and $R$ a constant matrix, whose eigenvalues determines the stability/instability of the system. Unfortunately, Floquet theory is not "constructive", so computing the latter decomposition is often a daunting task.

Answer 1

Okay, I think I can show that the origin is stable below. Sorry about the messy formatting.

First, let $$ y=\begin{bmatrix}0 &1\\ 1 & 0 \end{bmatrix}\exp\left( -\begin{bmatrix} -t & 0\\ 0 &-t \end{bmatrix} - \frac{1}{\omega} \sin (\omega t) \begin{bmatrix} -1 & 1 \\ 1 & -1 \end{bmatrix}\right)x$$ and cleaning this up we have $$y = e^{t + \frac{1}{\omega} \sin (\omega t)} \begin{bmatrix} -\sinh(\frac{1}{\omega} \sin(\omega t)) & \cosh(\frac{1}{\omega} \sin(\omega t)) \\ \cosh(\frac{1}{\omega} \sin(\omega t)) & -\sinh(\frac{1}{\omega} \sin(\omega t))\end{bmatrix} x. $$ So we are using the fundamental matrix $X(t) = e^{-t - \frac{1}{\omega} \sin (\omega t)} \begin{bmatrix} \cosh(\frac{1}{\omega} \sin(\omega t)) & -\sinh(\frac{1}{\omega} \sin(\omega t)) \\ -\sinh(\frac{1}{\omega} \sin(\omega t)) & \cosh(\frac{1}{\omega} \sin(\omega t))\end{bmatrix}$ of the solution to the original equation for $\gamma = 1$ to get a simpler equation equation for $y$, and then multiplying it by the matrix $\begin{bmatrix}0 &1\\ 1 & 0 \end{bmatrix}$ to simplify the anlysis below.

Taking the derivative of $y$ gives the messy expression: $$ \dot y = e^{t + \frac{1}{\omega} \sin (\omega t)} \begin{bmatrix} -\sinh(\frac{1}{\omega} \sin(\omega t)) & \cosh(\frac{1}{\omega} \sin(\omega t)) \\ \cosh(\frac{1}{\omega} \sin(\omega t)) & -\sinh(\frac{1}{\omega} \sin(\omega t))\end{bmatrix} \left( \begin{bmatrix} 1+ \cos(\omega t) & 0 \\ 0 & 1+\cos(\omega t)\end{bmatrix} + \cos(\omega t)\begin{bmatrix} \sinh(\frac{1}{\omega} \sin(\omega t)) & \cosh(\frac{1}{\omega} \sin(\omega t)) \\ \cosh(\frac{1}{\omega} \sin(\omega t)) & \sinh(\frac{1}{\omega} \sin(\omega t)) \end{bmatrix} \begin{bmatrix} -\cosh(\frac{1}{\omega} \sin(\omega t)) & \sinh(\frac{1}{\omega} \sin(\omega t)) \\ \sinh(\frac{1}{\omega} \sin(\omega t)) & -\cosh(\frac{1}{\omega} \sin(\omega t)) \end{bmatrix} + \begin{bmatrix} -1 & 0 \\ 0 & -\gamma \end{bmatrix} + \cos(\omega t) \begin{bmatrix} -1 & 1 \\ 1 & -1 \end{bmatrix}\right)x,$$ where the two matrices with hyperbolic sines and cosines multiply to give $\begin{bmatrix}0 & 1 \\ 1 & 0 \end{bmatrix}$, hence $$\dot{y} = \epsilon \begin{bmatrix} -\sinh(\frac{1}{\omega} \sin(\omega t)) & \cosh(\frac{1}{\omega} \sin(\omega t)) \\ \cosh(\frac{1}{\omega} \sin(\omega t)) & -\sinh(\frac{1}{\omega} \sin(\omega t)) \end{bmatrix} \begin{bmatrix} 0& 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} \sinh(\frac{1}{\omega} \sin(\omega t)) & \cosh(\frac{1}{\omega} \sin(\omega t)) \\ \cosh(\frac{1}{\omega} \sin(\omega t)) & \sinh(\frac{1}{\omega} \sin(\omega t))\end{bmatrix}y $$ where $\epsilon= 1-\gamma$. Multiplying this out, we finally get

$$ \dot{y} = \epsilon \begin{bmatrix} \cosh^2(\frac{1}{\omega} \sin(\omega t)) & \cosh(\frac{1}{\omega} \sin(\omega t)) \sinh(\frac{1}{\omega} \sin(\omega t)) \\ -\cosh(\frac{1}{\omega} \sin(\omega t)) \sinh(\frac{1}{\omega} \sin(\omega t)) & -\sinh^2(\frac{1}{\omega} \sin(\omega t))\end{bmatrix}y. $$

Let $$A(t) = \begin{bmatrix} \cosh^2(\frac{1}{\omega} \sin(\omega t)) & \cosh(\frac{1}{\omega} \sin(\omega t)) \sinh(\frac{1}{\omega} \sin(\omega t)) \\ -\cosh(\frac{1}{\omega} \sin(\omega t)) \sinh(\frac{1}{\omega} \sin(\omega t)) & -\sinh^2(\frac{1}{\omega} \sin(\omega t))\end{bmatrix}.$$ Note that $\epsilon = 0$ gives $y$ constant, which makes sense as it describes a perturbation of the original equation for $\gamma=1$, and further that $A(t)$ is periodic, so by Floquet's theorem we have a periodic (of period $T= \frac{2\pi}{\omega}$) matrix $P(t)$ and a constant matrix $B$ such that $y(t) = P(t)e^{tB}y(0)$ and $\text{tr}B = \epsilon$. Denote the fundamental matrix for this system by $Y(t) : = P(t)e^{tB}$.

Let $a := Y(t)a_0$ and $b := Y(-t)b_0$, for some constant vectors $a_0$, $b_0$. Since $A(-t) = A(t)^T$, we have $\dot{b}- = -\epsilon A(t)^T b $, and so the derivative of the inner product $\langle a,b\rangle$ is

$$ \frac{d}{dt} \langle a,b \rangle = \langle \epsilon A a,b\rangle + \langle a, -\epsilon A^T b\rangle = 0, $$

so $\langle a,b\rangle = \langle a_0, b_0 \rangle$. Hence, we have $P(t)e^{tB} (P(-t)e^{-tB})^T = I$, and therefore

$$e^{tB}e^{-tB^T} =P(t)^{-1} (P(t)^T)^{-1} .$$

Since both sides above must be periodic, both are equal to the identity $I$. Hence $B=B^T$ is symmetric and has real eigenvalues and eigenvectors. Since $\text{tr} B = \epsilon>0$ for $0<\gamma<1$, one eigenvalue is positive. Let us call $\lambda_+$ the largest eigenvalue of $B$, with eigenvector $v_+$. Similarly, call the other eigenvalue $\lambda_-$, with eigenvector $v_-$.

Now, note that the tangent of an integral curve (i.e. the graph of a solution) is given by $m(t):=\frac{\dot{y}_2}{\dot{y}_1} = - \tanh (\frac{1}{\omega}\sin(\omega t))$, and so $|m(t)| \leq |\tanh(\frac{1}{\omega})|<1$ for all $t \in \mathbb{R}$.

Edit to clarify this: the equations for $\dot{y}$ are

$$\dot{y}_1 = \epsilon \cosh(\frac{1}{\omega} \sin(\omega t))\left(\cosh(\frac{1}{\omega} \sin(\omega t))y_1 + \sinh(\frac{1}{\omega} \sin(\omega t))y_2 \right) $$

and

$$\dot{y}_2 = -\epsilon \sinh(\frac{1}{\omega} \sin(\omega t)) \left(\cosh(\frac{1}{\omega} \sin(\omega t))y_1 + \sinh(\frac{1}{\omega} \sin(\omega t))y_2\right).$$

Suppose $\lambda_-<0$. Then $P(t)e^{tB}v_- \to 0$ as $t \to \infty$. Now, since the slope of the integral curve is bounded, the origin must lie in the cone about a horizontal line through $v_-$ bounded by the lines of slope $\pm |\tanh(\frac{1}{\omega})|$ through $v_-$. Hence, $v_-$ lies inside a cone at the origin, which contains the $x$-axis and is bounded by lines of slope $\pm |\tanh(\frac{1}{\omega})|$ through the origin. Similarly, since $P(t)e^{tB}v_+ \to 0$ as $t\to -\infty$ and since its slope has the same bound, $v_+$ also lies in this cone at the origin. Hence, as the cone angle is less than $90^{\circ}$, the cosine of the angle between $v_+$ and $v_-$ is nonzero. So consider the solution $w(t):= P(t)e^{tB}(v_++ v_-)$. We now have $$\| v_+\|^2 +2 \langle v_+,v_-\rangle + \|v_-\|^2=\langle w(nT),w(-nT) \rangle = \langle e^{\lambda_+ nT}v_+ + e^{\lambda_- n T}v_-, e^{-\lambda_+ nT}v_+ + e^{-\lambda_- n T}v_- \rangle = \|v_+\|^2+2\cosh((\lambda_+ - \lambda_-)nT)\langle v_+,v_- \rangle+\|v_-\|^2$$ and so, since $\lambda_+ \neq \lambda_-$, we must have $\langle v_+,v_- \rangle=0$, a contradiction. Hence, $\lambda_- \geq 0$.

We are now done, as if both eigenvalues are positive and satisfy $\lambda_+ + \lambda_- = \epsilon$, both are smaller than $\epsilon$. So a solution $x$ to the perturbed system for $0<\gamma<1$ satisfies $\|x\| \leq e^{-t}e^{(1-\gamma)t}\|Q(t)x_0\|$, where $Q(t)$ is some periodic matrix, and hence $x(t) \to 0$ as $t \to \infty$.