Instability and bifurcation of a cooperative system with periodic coefficients

1.   Introduction

For the linear autonomous system $\dot{x} = Ax, x\in\mathbb{R}^n$, it is well known that the zero solution of the system is asymptotically stable when the real parts of the eigenvalues of coefficient matrix $A$ are all negative. However, the nonautonomous linear system $\dot{x} = A(t)x,x\in\mathbb{R}^n$ does not possess this property in general. In fact, even if all the eigenvalues of $A(t)$ are negative and bounded away from zero, the zero solution can still be unstable; and moreover, the solution of system may grow to infinity (see, e.g. [4, 8, 14]). Very recently, Mierczyński [11] constructed a series of such examples in cooperative system by using an unusual and ingenious method. Consider a 2-D linear ODE

with $A(t)$ is denoted by

where

and $c$ is a positive number. For the sake of convenience, one can also rewrite the system (1) as:

where

and

Recall that a linear system is cooperative (resp., strongly cooperative) if, for each $t\in \mathbb{R}$, the matrix $A(t)$ has nonnegative (resp., positive) off-diagonal entries. Clearly, system (1) is a (strongly) cooperative periodic linear system with period $2$. In the biological science and population ecology, cooperative behavior means that an increase in the density of one species will enhance the growth of the other, which abounds across all domains of life, from animals to microbes. This type of interaction among species is widely used by theoretical ecologists to explain biological phenomena (see [2, 5, 10, 12]). For example, system (1) is usually described by some bacterial populations that they always switch between two states (dormant vs. active). If we let $x_{1}$ stand for the density of bacteria in the dormant stage and $x_{2}$ for the density of bacteria in the active stage, then $a_{12}(t)$ (resp., $a_{21}(t)$) describes the transition rate from the active state to the dormant state (resp., from the dormant state to the active state) at time $t$. It is now well-known that a (strongly) cooperative differential equation generates a classical (strongly) monotone dynamical system (see [7, 13]).

Clearly, for each $t$, the eigenvalues of the matrix $A(t)$ are both $-\frac{1}{2}$ and $-\frac{3}{2}$. Mierczyński [11] has shown that the zero solution of such a system can be unstable for any $c>0$ sufficiently large. More precisely, let $X(t;s)$ stand for the transition matrix of system (1), and let $P$ be the Poincaré map of system (1). Then

It follows that

Note that

So, the Poincaré map $P$ can be directly calculated as

As a consequence, the principal eigenvalue of $P$ is larger than 1 when $c>2.13834$, which implies that any nontrivial solution of system (1) will grow to infinity. Thus, the zero solution is unstable. In addition, Mierczyński [11] adopted another more easier way, Peano-Baker series, to estimate the principal eigenvalue of $P$ and obtain the similar results.

In the present paper, for the linear $2$-periodic system (1), we will focus on the optimization of a switching strategy for instability between the matrices $A^{(1)}$ and $A^{(2)}$ in a period of time. For this purpose, a switching strategy parameter $\lambda \in (0,2)$ is introduced into system (1) as:

where

The main purpose of this paper is to analyze the instability of the zero solution of the above system with respect to the switching parameter $\lambda\in(0,2)$. We obtain a complete distribution in $\lambda$-$c$ parameter plane corresponding to the instability of the zero solution in terms of simple inequality of the parameters $\lambda$ and $c$ (see Theorem 2.2). Furthermore, we present the bifurcation diagram (see Figure 2.1). The optimization value of the parameter $\lambda$ is also given in the unstable region (see Theorem 2.3).

The paper is organized as follows. Section 2 is devoted to study the relationship between parameters $\lambda$ and $c$ and the instability of the zero solution by using Poincaré map and transition matrix. Further, we gain the optimal strategy values of unstable branches of the zero solution and present the branch graph of instability. In Section 3, we provide some numerical simulations to illustrate the main theoretical results. This paper ends with a discussion in Section 4.

2.   Instability and bifurcation

Consider the following linear cooperative system with a parameter $\lambda$,

where

Here, the parameter $\lambda$ represents the strategy parameter switching between two systems in a period.

Clearly, system (1) is a special case of system (2) (with $\lambda$ = 1). The eigenvalues of $A_{\lambda}(t)$ are both $-\frac{1}{2}$ and $-\frac{3}{2}$ for each $t\in \mathbb{R}$. A direct calculation yields the following transition matrix associated with system (2) as

When restricted to the interval [0, 2), the transition matrix is reduced to

Let us first define the Poincaré map $Q$ of system (2), that is,

Since $X(2k+2;2k) = X(2;0)$ for any $k\in \mathbb{Z}$, we have

Let $\mathcal{M}$ stand for the family of real $2\times 2$ matrices with off-diagonal entries positive, and let $\mathcal{P}$ stand for the family of real $2\times2$ matrices with all entries positive. For any $A\in \mathcal{M}$, one has the following proposition.

Lemma 2.1. Let $A = [a_{ij}]_{i,j = 1}^{2}\in \mathcal{M}.$ Then the following three statements are true:

(i) $e^{tA}\in \mathcal{P}$ for all $t>0$.

(ii) $A$ has two real eigenvalues (denoted $\mu<\theta$). In particular, if $A \in \mathcal{P}$, then the principal eigenvalue $\theta$ is simple and $\theta>|\mu|$, where the principal eigenvalue is the largest of the eigenvalues of $A$.

(iii) An eigenvector $\xi$ corresponding to $\theta$ can be taken to have its coordinates positive.

Proof. For (i)-(iii), see Proposition 2.2 and Theorem 2.5 in [11]. If $A \in \mathcal{P}$, it easily follows from Theorem 1.2 in [3] that the principal eigenvalue $\theta$ is simple and $\theta>|\mu|$.

By virtue of Lemma 2.1, for any $\lambda \in (0,2)$, both $e^{\lambda A^{(1)}}$ and $e^{(2-\lambda)A^{(2)}}$ belong to $\mathcal{P}$; and hence, so is the Poincaré map $Q$. Furthermore, there exists a unique normalized principal eigenvector $\xi$ corresponding to the principal eigenvalue $\theta$ of $Q$ with $\xi \in \rm{Int}\mathbb{R}^{2}_{+}$ and $\theta >0$, where $\rm{Int}\mathbb{R}^{2}_{+}$ $: = \{(x_1,x_2)\in \mathbb{R}^2: x_1 > 0, x_2 > 0 \}$. Denote by $x(t)$ the solution of system (2) taking initial value $x$ at time $t = 0$. $x(t)$ is called to be Lyapunov stable (abbr. stable) if, given $\epsilon>0$, there exists some $\delta = \delta(\epsilon)>0$ such that, for any other solution $y(t)$ of system (2) satisfying $|x-y|<\delta$, then $|x(t)-y(t)|<\epsilon$ for any $t\geq0$. A solution which is not stable is said to be unstable.

For the special initial value $\xi$, one has

If $\theta>1$, then the solution $\xi(t)$ will satisfy:

which entails that the solution $\xi(t)$ is an unstable solution of system (2).

Therefore, one needs to find the relationship between parameters $\lambda$ and $c$ so that $\theta>1$. Our main result is given below, which provides a necessary and sufficient criterion for the instability of system (2).

Theorem 2.2. For $\lambda \in (0,2)$ and $c>0$, the zero solution of system $\rm(2)$ is unstable if and only if

Proof. Sufficiency. Recall that $Q = X(2;0) = e^{(2-\lambda)A^{(2)}}\bullet e^{\lambda A^{(1)}}$. A tedious calculation yields that

with

which are all positive entries. Since the matrix $Q$ belongs to $\mathcal{P}$, it follows from Lemma 2.1(ii)-(iii) that $\theta>0$ is a simple eigenvalue with a positive eigenvector $\xi$. Let $\mu$ be the other eigenvalue of $Q$. By Lemma 2.1(ii), we have $\theta>|\mu|$. More precisely, $\theta$ and $\mu$ satisfy the following quadratic equation related to the variable $m$,

Solving the above equation, we obtain

Note that

Then

Consequently, we obtain that $\theta>1$ if and only if

Simplifying the above inequality, we obtain that $\theta>1$ if and only if

Hence, when parameters $\lambda$ and $c$ satisfy the inequality (3), the solution $\xi(t)$ satisfies $\xi(2n) = Q^{n}\xi = \theta^{n}\xi\rightarrow\infty$ as $n\rightarrow\infty$. So, the zero solution of system (2) is unstable.

Necessity. If the zero solution of system (2) is unstable, we will show that parameters $\lambda$ and $c$ satisfy the inequality (3). Suppose on the contrary that $(2c-\frac{1}{2c})^{2}[\cosh1-\cosh(1-\lambda)]\leq 4(\cosh2-\cosh1)$. From the proof of the sufficiency, one can get $\theta\leq1$; and hence, $|\mu|<\theta\leq1$. Let $\zeta$ be the eigenvector corresponding to $\mu$, so $\xi$ and $\zeta$ are linearly independent. Thus for any initial value $x_{0}\in \mathbb{R}^{2}$, there exist $a,b\in\mathbb{R}$ such that $x_{0} = a\xi+b\zeta$. Then the solution $x(t;0,x_{0})$ of system (2) satisfies

As a consequence, it is easy to see that if $\theta<1$ then $x(2n;0,x_{0})\rightarrow 0$ as $n\rightarrow\infty$, which entails that the zero solution of system (2) is (asymptotically) stable. This is a contradiction to our assumption. On the other hand, if $\theta = 1$, then for any $\epsilon>0$, one can choose $a$ and $b$ such that $|x_{0}|<\delta = \epsilon$. Since $|\mu|<1$, there exists an integer $N>0$ such that $|b\mu^{n}\zeta|<\epsilon$ for any $n>N$. Let $P$ be the projection operator onto the $\xi$-direction. Noticing that the projection operator $P$ is bounded, there exists a positive constant $C$ such that $|x(2n;0,x_{0})| = |a\xi+b\mu^{n}\zeta|\leq|a\xi|+|b\mu^{n}\zeta| = |Px_{0}|+|b\mu^{n}\zeta|\leq C|x_{0}|+|b\mu^{n}\zeta|<(2C+1)\epsilon$ as $n>N$. Hence, it follows that the zero solution of system (2) is Lyapunov stable. This also leads to a contradiction. Thus, we have completed the proof.

According to Theorem 2.2, we write

as the unstable region. We have the following optimization of the strategy parameter $\lambda$ as:

Theorem 2.3. Let $U_1 = \{(\lambda,c)\in U: c<\frac{1}{2}\}$ and $U_2 = \{(\lambda,c)\in U: c>\frac{1}{2}\}$. For $\lambda\in(0,2)$ and $c>0$, we have the following statements:

(i) In unstable region $U_1$, $c$ attains the maximum at $\lambda = 1$.

(ii) In unstable region $U_2$, $c$ attains the minimum at $\lambda = 1$.

Proof. By the proof of Theorem 2.2, the critical relation between parameters $\lambda\in(0,2)$ and $c>0$ is

We rewrite as

Then,

For simplicity, we denote

(i). If $0<c<\frac{1}{2}$, then the equation (4) becomes the following form

which implies that $c = \frac{-m(\lambda)+\sqrt{m^{2}(\lambda)+1}}{2}$. Due to the continuity and monotonicity of $m(\lambda)$ with respect to the variable $\lambda$, it follows that if $\lambda\in(0,1)$ then $c$ is increasing as $\lambda$ increases. Meanwhile, if $\lambda\in(1,2)$ then $c$ is decreasing as $\lambda$ increases. So, when $\lambda = 1$, $c$ will attain the maximum in unstable region $U_{1}$. This implies that the statement (i) holds.

(ii). If $c>\frac{1}{2}$, then the equation (4) turns out to be

which implies that $c = \frac{m(\lambda)+\sqrt{m^{2}(\lambda)+1}}{2}$. Similarly, if $\lambda\in(0,1)$ then $c$ is decreasing as $\lambda$ increases. If $\lambda\in(1,2)$ then $c$ is increasing as $\lambda$ increases. So, when $\lambda = 1$, $c$ will attain the minimum in unstable region $U_{2}$, which implies the statement (ii).

Remark 1. By using Matlab tools, we obtain the $\lambda$-$c$ bifurcation diagram of instability for system (2) as demonstrated in Figure 2.1.

(i) When parameters $\lambda$ and $c$ are chosen in the red region $U$(= $U_{1}\cup U_{2}$), the zero solution of system (2) is unstable.

(ii) When parameters $\lambda$ and $c$ are chosen in the green region $S$, the zero solution of system (2) is (asymptotically) stable.

(iii) When parameters $\lambda$ and $c$ are chosen on the boundary between the red region $U$ and the green region $S$, the zero solution of system (2) is Lyapunov stable (but not asymptotically stable).

From such bifurcation diagram, one can see that the parameter $c$ has the maximum and minimum with respect to variable $\lambda \in (0,2)$ in the unstable region $U$.

3.   Numerical simulation

In this section, we provide several numerical simulations to illustrate our main results. First, we randomly select the initial values of parameters $\lambda$ and $c$ in unstable region $U = U_{1}\cup U_{2}$ to make the numerical fitting analysis.

When taking $\lambda = \frac{1}{2}<1, c = 0.1$ in the unstable region $U_{1}$ and initial value $\xi = (7.7,2.172)$, the solution $\xi(t)$ of system (2) is shown in Figure 3.1.

By taking $\lambda = \frac{4}{3}>1, c = 5$ in the unstable region $U_{2}$ and initial value $\xi = (2.85,9.585)$, the solution $\xi(t)$ of system (2) is shown in Figure 3.2.

From Figures 3.1-3.2, it is clear that the solutions will grow to infinity as $t$ increases gradually. We conclude the zero solution of system (2) is unstable.

Next, we will select parameters $\lambda$ and $c$ in stable region $S$. Taking $\lambda = \frac{1}{2}<1, c = 2$ and initial value $\xi = (3.209,9.471)$, the solution $\xi(t)$ of system (2) is shown in Figure 3.3.

Now, we take $\lambda = \frac{4}{3}>1, c = 2$ and initial value $\xi = (5.282,8.491)$, the solution $\xi(t)$ of system (2) is shown in Figure 3.4.

From Figure 3.3-3.4, one can see that all the solutions are asymptotic to zero as $t$ increases to infinity when the values of $\lambda$ and $c$ are chosen in region $S$.

4.   Discussion

In this paper, we analyze the instability of the zero solution of a cooperative differential system (2) with periodic coefficients. By introducing a switching strategy parameter $\lambda$, we investigate the bifurcation for instability of the zero solution. Theorem 2.2 gives a necessary and sufficient condition, which guarantees the instability of zero solution. By taking $\lambda \in (0,2)$ as a variable, we obtain that the parameter $c$ can attain the maximum at $\lambda = 1$ in unstable region $U_1$, while it attains the minimum at $\lambda = 1$ in unstable region $U_2$ in Theorem 2.3. In other words, $\lambda = 1$ will attain the optimal value either in $U_{1}$ or in $U_{2}$. In particular, when $\lambda = 1$, the system (2) is reduced to system (1), which has been investigated in Mierczyński [11].

The parameters $\lambda$-$c$ bifurcation graph in Figure 2.1 clearly reveals the unstable region and stable region of the zero solution. Moreover, we observe that the instability of the zero solution will disappear gradually when the parameter $\lambda$ tends to $0$ or $2$. In order to illustrate our results more vividly, we present numerical simulations to show the main conclusions. Based on the above analysis, we can also make suitable modification on the entries of $A^{(1)}$ and $A^{(2)}$ to obtain more analogous results.

In addition, we should bear in mind that the time-periodic differential equations are due to biological applications, such as the results of seasonal changes, availability of food. So it should be emphasized that analogous constructions could be made for quasi-periodic, almost periodic; and more general dependence on time. Meanwhile, in our system (2), the switching between the two matrices $A^{(1)}$ and $A^{(2)}$ in a period only occurs at one fixed time $\lambda$ in the interval (0, 2). If the switching time is allowed to be the random variable, it is also possible to construct more analogous systems (see, e.g. [1, 6, 9]). Another possible extension is to consider higher dimension dynamical systems. We will leave them for future research.

Acknowledgments

The authors are grateful to two anonymous referees for their valuable comments and suggestions which led to an improvement of our original manuscript.