Consider the following system of coupled differential equations $$ \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}, $$ 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.

  • $\begingroup$ Interesting! could you add some details on the numerical simulations you've used (some code for example) ? $\endgroup$ Jul 14 '18 at 19:07
  • 1
    $\begingroup$ @KonstantinosKanakoglou: I've edited the question adding the Mathematica code that I've used in my simulations. $\endgroup$
    – Ludwig
    Jul 14 '18 at 19:15
  • $\begingroup$ I think that in your edits the matrix $$\begin{bmatrix}e^t&0\\ 0&e^{\lambda t}\end{bmatrix}$$ to the left of $x$ is missing. Further, how do you have that a fundamental matrix is given by $\exp(\int_0^t A(\tau)\,d\tau)$? The matrices $$\begin{bmatrix}-1&e^{(1-\gamma)t}\\ e^{(\gamma-1)t}&-1\end{bmatrix}$$ don't commute for $t_1\ne t_2$ unless $\lambda=1$. Of course, commutation is not a necessary condition. $\endgroup$
    – user539887
    Jul 19 '18 at 12:54
  • 1
    $\begingroup$ A generic Floquet problem is not solvable analytically. $\endgroup$ Jul 21 '18 at 3:27
  • 1
    $\begingroup$ @PiyushGrover Is that a statement in differential algebra? Could you give some references? $\endgroup$
    – user539887
    Jul 21 '18 at 9:51

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


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

  • $\begingroup$ Thanks for your answer! Could you please elaborate a little more on the derivation of the closed-form expression of the matrix exponential in your second formula? (Because, applying the expression in Theorem 2.2 of this paper, I get $$e^{t+\frac{1}{\omega}\sin(\omega t)}\begin{bmatrix} \cosh(\sin(\omega t)/\omega) & -\sinh(\sin(\omega t)/\omega) \\ -\sinh(\sin(\omega t)/\omega) & \cosh(\sin(\omega t)/\omega) \end{bmatrix},$$ which looks different from yours..) $\endgroup$
    – Ludwig
    Jul 26 '18 at 3:31
  • $\begingroup$ Yes, I made a mistake. However, this can be fixed by multiplying $y$ by $\begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}$, I will add this later. $\endgroup$ Jul 26 '18 at 15:34
  • $\begingroup$ Thanks for fixing this. Also, it is not clear to me why $m(t):=\dot{y}_2/\dot{y}_1=-\tanh(\sinh(\omega t)/\omega)$. Could you please expand this a little bit? $\endgroup$
    – Ludwig
    Jul 26 '18 at 17:47
  • $\begingroup$ Sure, I have added the formula for $\dot{y}_1$ and $\dot{y}_2$. $\endgroup$ Jul 26 '18 at 21:44
  • $\begingroup$ Okay thanks (just forgot a square factor in the diagonal in my calculations...) I'll check the remaining part of the proof asap $\endgroup$
    – Ludwig
    Jul 26 '18 at 22:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.