*Article* **Global Stability of Delayed Ecosystem via Impulsive Differential Inequality and Minimax Principle**

**Ruofeng Rao 1,2**


**Abstract:** This paper reports applying Minimax principle and impulsive differential inequality to derive the existence of multiple stationary solutions and the global stability of a positive stationary solution for a delayed feedback Gilpin–Ayala competition model with impulsive disturbance. The conclusion obtained in this paper reduces the conservatism of the algorithm compared with the known literature, for the impulsive disturbance is not limited to impulsive control.

**Keywords:** Minimax principle; linear approximation theory; ecosystem; steady state solution

#### **1. Introduction**

It is well known that Gilpin–Ayala competition model (GACM) has been hotly discussed (see in [1–7]) due to its importance in simulating two or more competing biological populations in nature. As diffusion is an essential characteristic of most biological populations, Ling Bai and Ke Wang began to investigate the global stability of reaction-diffusion Gilpin–Ayala ecosystem under Neumann zero boundary value in 2005 (see in [8]), and obtained good results. Actually, Neumann zero boundary value means that the populations do not migrate beyond the biosphere boundary. However, many animal populations are at the edge of the biosphere, where the population density is usually zero, which is not reflected by Neumann zero boundary value. Thus, the Dirichlet zero boundary value was considered in recent literature [6,7]. In recent years, linear impulsive control and nonlinear impulsive control technology are widely used in ordinary differential dynamical systems and infinite dimensional dynamical systems [9–16]. For example, in [14], event-triggered nonlinear impulsive control to stabilize damped wave equations was designed, and rapid exponential stabilization was achieved. However, in this paper, linear impulse is relatively simple and feasible in practical ecological management because ecological management is a natural system of impulse artificial intervention. Therefore, the linear impulsive control is considered in this paper. Note that impulse control is employed to make the GACM stable globally in [6,7], but this paper involves the impulsive disturbance, which is not limited to impulsive control. Minimax principle will be employed to derive the existence of multiple stationary solutions, which improve the method of Mountain Pass Lemma in [6]. On the other hand, impulsive disturbance is considered in this paper, not just impulse control. In fact, some impulse management measures other than impulsive control sometimes occur in ecological management due to accidents, such as releasing animals, hunting animals harmful to the population, and so on. These pulse measures mean that the pulse intensity is not necessarily less than 1 based on system stability.

#### **2. Preparatory Knowledge**

Consider the following reaction-diffusion Gilpin–Ayala competition model (RDGACM) with delayed feedback under Dirichlet boundary value

**Citation:** Rao, R. Global Stability of Delayed Ecosystem via Impulsive Differential Inequality and Minimax Principle. *Mathematics* **2021**, *9*, 1943. https://doi.org/10.3390/math9161943

Academic Editor: Mustafa R.S. Kulenovic

Received: 11 July 2021 Accepted: 12 August 2021 Published: 14 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

$$\begin{cases} \begin{aligned} \frac{\partial u\_1}{\partial t} &= d\_1 \Delta u\_1 + u\_1 (b\_1 - a\_{11} u\_1^{\theta\_1} - a\_{12} u\_2) - k\_1 (u\_1 - u\_1 (t - \tau\_1, x)), \quad t \geqslant 0, \, x \in \Omega, \\\ \frac{\partial u\_2}{\partial t} &= d\_2 \Delta u\_2 + u\_2 (b\_2 - a\_{21} u\_1 - a\_{22} u\_2^{\theta\_2}) - k\_2 (u\_2 - u\_2 (t - \tau\_2, x)), \quad t \geqslant 0, \, x \in \Omega, \\\ u\_1(t, x) &= u\_2(t, x) = 0, \quad t \geqslant 0, \, x \in \partial\Omega, \end{aligned} \end{cases} \tag{1}$$

equipped with the initial value

$$
\mu\_1(\mathbf{s}, \mathbf{x}) = \xi\_1(\mathbf{x}), \ \mu\_2(\mathbf{s}, \mathbf{x}) = \xi\_2(\mathbf{x}), \ \mathbf{s} \in [0, \pi\_0], \ \mathbf{x} \in \Omega,\tag{2}
$$

where Ω is a bounded domain in R*N*(1 *N* 3) with smooth boundary *∂*Ω. Time delays *τ*1, *τ*<sup>2</sup> ∈ [0, *τ*0], *ui*(*t*, *x*) represents the population density of the *i*th population at time *t* and the spatial location *x*, *bi* > 0 represents the birth rate of the population of the *i*th species, and *aij* > 0 represents the competition parameter between the species *i* and the species *j*. *di* > 0 represents the diffusion coefficient for the species *i*. Initial value function *ξ*(*x*)=(*ξ*1(*x*), *ξ*2(*x*))*<sup>T</sup>* is bounded and continuous.

Assume that

(A1) For *<sup>i</sup>* <sup>=</sup> 1, 2, setting <sup>−</sup><sup>1</sup> <sup>&</sup>lt; *<sup>θ</sup><sup>i</sup>* <sup>&</sup>lt; 4, and *<sup>s</sup>*2+*θ<sup>i</sup>* -0 for all *<sup>s</sup>* ∈ *<sup>R</sup>*1.

(A2) For *i* = 1, 2, there exist positive constants *Mi* > 0 such that

$$0 \ll u\_i \ll M\_i.$$

(A3) For *i* = 1, 2, |∇*ui*(*t*, *x*)| is bounded for all *x* ∈ Ω.

Due to the limited natural resources, it is reasonable to assume in (A2) that each population density is limited. Besides, the limited natural resources imply that the boundedness assumption of (A3) is suitable to the real state of nature.

**Remark 1.** *(A1) expands greatly the allowable range of parameters θi, compared with the previous related literature (see, e.g., in [6,7]). For example, the harsh condition "*0 < *θ<sup>i</sup>* < 1 *with* ˆ *θ<sup>i</sup> being an even number, and* ˇ *θ<sup>i</sup> being an odd number" is deleted.*

**Lemma 1** (see, e.g., in [17])**.** *Let <sup>J</sup>* ∈ *<sup>C</sup>*1(*H*<sup>1</sup> <sup>0</sup> (Ω), <sup>R</sup>1). *If there is an upper boundness of <sup>J</sup> in H*<sup>1</sup> <sup>0</sup> (Ω)*, and J satisfies the (PS) condition, then c*<sup>∗</sup> = sup *<sup>v</sup>*∈*H*<sup>1</sup> <sup>0</sup> (Ω) *J*(*v*) *is a critical value of J.*

Here, the (*PS*)*<sup>c</sup>* condition may be found in [18] (Definition 2). Actually, the (PS) condition is equivalent to the (*PS*)*<sup>c</sup>* condition. For convenience, the author describes the (PS) condition as follows:

**Definition 1** ([17])**.** *Let ψ be a real C*<sup>1</sup> *functional defined on a Banach space X. If any sequence* {*un*} *in X with ψ* (*un*) *<sup>X</sup>*<sup>∗</sup> <sup>→</sup> <sup>0</sup> *and the bounded sequence* {*ψ*(*un*)}<sup>∞</sup> *<sup>n</sup>*=<sup>1</sup> *has a convergent subsequence in X, then ψ is called satisfying the (PS) condition.*

**Lemma 2** ([7], Theorem 3.1)**.** *Set u*∗(*x*)=(*u*∗ <sup>1</sup> (*x*), *u*<sup>∗</sup> <sup>2</sup> (*x*))*T. Suppose that the condition (A2) holds, and* 0 < *θ<sup>i</sup>* < 1 *for i* = 1, 2. *Moreover, if there exists a positive constant c*<sup>∗</sup> > 0 *such that*

$$0 \ll h(\mu^\*(x)) \lesssim c\_\* D \mathfrak{H}\_\prime \tag{3}$$

*then there are at least a positive bounded equilibrium solution u*∗(*x*) *for the RDAGCM (1), where* H = (1, 1)*T, h*(*u*)=(*h*1(*u*1, *u*2), *h*2(*u*1, *u*2))*<sup>T</sup> with u* = (*u*1, *u*2)*<sup>T</sup> and*

$$h\_1(u\_1, u\_2) = u\_1(b\_1 - a\_{11}u\_1^{\theta\_1} - a\_{12}u\_2), \quad h\_2(u\_1, u\_2) = u\_2(b\_2 - a\_{21}u\_1 - a\_{22}u\_2^{\theta\_2}),$$

$$D = \left(\begin{array}{cc} d\_1 & 0\\ 0 & d\_2 \end{array}\right) > 0.$$

The conditions of Lemma 2 guarantee the existence of a positive stationary solution (*u*∗ <sup>1</sup> (*x*), *u*<sup>∗</sup> <sup>2</sup> (*x*)) for the delayed feedback system (1). Set

$$\begin{cases} \mathcal{U}\_1 = u\_1 - u\_1^\*(\mathfrak{x}) \\ \mathcal{U}\_2 = u\_2 - u\_2^\*(\mathfrak{x}) \end{cases}$$

and the stationary solution (*u*∗ <sup>1</sup> (*x*), *u*<sup>∗</sup> <sup>2</sup> (*x*)) of the system (1) corresponds to the zero solution (0, 0)*<sup>T</sup>* of the following system:

$$\begin{cases} \frac{\partial \mathcal{U}\_{1}}{\partial t} = d\_{1} \Delta \mathcal{U}\_{1} + b\_{1} \mathcal{U}\_{1} - \Phi\_{1}(\mathcal{U}\_{1}, \mathcal{U}\_{2}) - k\_{1} [\mathcal{U}\_{1} - \mathcal{U}\_{1}(t - \tau\_{1}, x)], \quad t \geqslant 0, \, x \in \Omega, \\\frac{\partial \mathcal{U}\_{2}}{\partial t} = d\_{2} \Delta \mathcal{U}\_{2} + b\_{2} \mathcal{U}\_{2} - \Phi\_{2}(\mathcal{U}\_{1}, \mathcal{U}\_{2}) - k\_{2} [\mathcal{U}\_{2} - \mathcal{U}\_{2}(t - \tau\_{2}, x)], \quad t \geqslant 0, \, x \in \Omega, \\\mathcal{U}\_{1}(t, \boldsymbol{x}) = \mathcal{U}\_{2}(t, \boldsymbol{x}) = 0, \quad t \geqslant 0, \, \boldsymbol{x} \in \partial \Omega, \end{cases}$$

or

$$\begin{cases} \frac{\partial Ll\_1}{\partial t} = d\_1 \Delta l l\_1 + (b\_1 - k\_1)lL\_1 - \Phi\_1(lL\_1, lL\_2) + k\_1 lL\_1(t - \tau\_1, \mathbf{x}), \quad t \gg 0, \, \mathbf{x} \in \Omega, \\\frac{\partial Ll\_2}{\partial t} = d\_2 \Delta l l\_2 + (b\_2 - k\_2)lL\_2 - \Phi\_2(lL\_1, lL\_2) + k\_2 lL\_2(t - \tau\_2, \mathbf{x}), \quad t \gg 0, \, \mathbf{x} \in \Omega, \\\ l l\_1(t, \mathbf{x}) = ll\_2(t, \mathbf{x}) = 0, \quad t \gg 0, \, \mathbf{x} \in \partial\Omega, \end{cases} (4)$$

where we denote *U* = (*U*1, *U*2)*T*, and

$$\begin{aligned} \Phi\_1(\mathcal{U}) &= (\mathcal{U}\_1 + u\_1^\*(\mathbf{x})) [a\_{11}(\mathcal{U}\_1 + u\_1^\*(\mathbf{x}))^{\theta\_1} + a\_{12}(\mathcal{U}\_2 + u\_2^\*(\mathbf{x}))] - u\_1^\*(\mathbf{x}) (a\_{11} u\_1^\*(\mathbf{x})^{\theta\_1} + a\_{12} u\_2^\*(\mathbf{x})), \\ \Phi\_2(\mathcal{U}) &= (\mathcal{U}\_2 + u\_2^\*(\mathbf{x})) [a\_{21}(\mathcal{U}\_1 + u\_1^\*(\mathbf{x})) + a\_{22}(\mathcal{U}\_2 + u\_2^\*(\mathbf{x}))^{\theta\_2}] - u\_2^\*(\mathbf{x}) (a\_{21} u\_1^\*(\mathbf{x}) + a\_{22} u\_2^\*(\mathbf{x})^{\theta\_2}). \end{aligned} \tag{5}$$

The following system is the system (4) in form of vector-matrix:

$$\begin{cases} \frac{\partial \mathcal{U}}{\partial t} = D \Delta \mathcal{U} + (B - K)\mathcal{U} - \Phi(\mathcal{U}) + KL(t - \tau, \mathbf{x}), \quad t \gg 0, \ \mathbf{x} \in \Omega, \\\mathcal{U}(t, \mathbf{x}) = 0, \quad t \gg 0, \mathbf{x} \in \partial \Omega, \end{cases} \tag{6}$$

where*<sup>U</sup>* = (*U*1, *<sup>U</sup>*2)*T*, *<sup>U</sup>*(*<sup>t</sup>* − *<sup>τ</sup>*, *<sup>x</sup>*)=(*U*(*<sup>t</sup>* − *<sup>τ</sup>*1, *<sup>x</sup>*), *<sup>U</sup>*(*<sup>t</sup>* − *<sup>τ</sup>*2, *<sup>x</sup>*))*T*, <sup>Φ</sup>(*U*)=(Φ1(*U*), <sup>Φ</sup>2(*U*))*<sup>T</sup>* and

$$D = \begin{pmatrix} d\_1 & 0 \\ 0 & d\_2 \end{pmatrix}, \ B = \begin{pmatrix} b\_1 & 0 \\ 0 & b\_2 \end{pmatrix}, \ K = \begin{pmatrix} k\_1 & 0 \\ 0 & k\_2 \end{pmatrix}. \tag{7}$$

Considering the impulse disturbance on (6), one can get the following system:

$$\begin{cases} \frac{\partial \mathcal{U}}{\partial t} = D \Delta \mathcal{U} + (B - K)\mathcal{U} - \Phi(\mathcal{U}) + K\mathcal{U}(t - \tau, \mathbf{x}), \quad t \geqslant 0, \ t \neq t\_k, \mathbf{x} \in \Omega, \\\mathcal{U}(t\_k^+, \mathbf{x}) = A\_k \mathcal{U}(t\_k^-, \mathbf{x}), \quad k = 1, 2 \cdot \cdots \\\mathcal{U}(t, \mathbf{x}) = 0, \quad t \geqslant 0, \ \mathbf{x} \in \partial \Omega, \end{cases} \tag{8}$$

where {*tk*}<sup>∞</sup> *<sup>k</sup>*=<sup>1</sup> is a sequence of fixed impulsive instants, satisfying 0 < *t*<sup>1</sup> < *t*<sup>2</sup> < ··· < *tk* < *tk*<sup>+</sup><sup>1</sup> < ··· and lim *<sup>k</sup>*→<sup>∞</sup> *tk* = +∞. Besides, *Ui*(*<sup>t</sup>* + *<sup>k</sup>* , *x*) = *Ui*(*tk*, *x*), *Ui*(*t* − *<sup>k</sup>* , *x*) = lim *t*→*t* − *k Ui*(*t*, *x*) for all *i* = 1, 2, *k* = 1, 2, ··· .

**Definition 2.** (*u*∗ <sup>1</sup> (*x*), *u*<sup>∗</sup> <sup>2</sup> (*x*))*<sup>T</sup> is said to be globally exponentially stable under impulsive disturbances if the zero solution of the system (8) is globally exponentially stable.*

**Lemma 3** (see [19])**.** *Consider the following differential inequality:*

$$\begin{cases} \begin{array}{ll} D^+v(t) & \leq -av(t) + b[v(t)]\_{\mathbb{T}'} \; t \neq t\_k \\ v(t\_k) & \leq a\_k v(t\_k^-) + b\_k[v(t\_k^-)]\_{\mathbb{T}'} \end{array} \end{cases}$$

*where v*(*t*) ≥ 0, [*v*(*tk*)]*<sup>τ</sup>* = sup *t*−*τ*≤*s*≤*t v*(*s*), [*v*(*t* − *<sup>k</sup>* )]*<sup>τ</sup>* = sup *t*−*τ*≤*s*<*t v*(*s*) *and v*(*t*) *is continuous except tk*, *k* = 1, 2, ··· , *where it has jump discontinuities. The sequence tk satisfies* 0 = *t*<sup>0</sup> < *t*<sup>1</sup> < *t*<sup>2</sup> < ··· < *tk* < *tk*<sup>+</sup><sup>1</sup> < ··· *, and* lim *<sup>k</sup>*→<sup>∞</sup> *tk* <sup>=</sup> <sup>∞</sup>. *Suppose that (1) a* > *b* ≥ 0;

$$\text{In (2) } t\_k - t\_{k-1} > \delta \text{\textquotedblleft, where } \delta > 1 \text{, and there exist constants } \gamma > 0 \text{, M > 0 such that}$$

$$
\rho\_1 \rho\_2 \cdots \rho\_{k+1} e^{k\lambda \tau} \precsim M e^{\gamma t\_k},\tag{9}
$$

*where <sup>ρ</sup><sup>i</sup>* <sup>=</sup> *max*{1, *ai* <sup>+</sup> *bieλτ*}*, <sup>λ</sup>* <sup>&</sup>gt; <sup>0</sup> *is the unique solution of equation <sup>λ</sup>* <sup>=</sup> *<sup>a</sup>* <sup>−</sup> *beλτ*; *then*

$$v(t) \ll M[v(0)]\_\tau e^{-(\lambda - \gamma)t}.$$

*In addition, if θ* = sup *k*∈*Z* {1, *ak* <sup>+</sup> *bkeλτ*}, *then*

$$w(t) \ll \theta[v(0)]\_\tau e^{-(\lambda - \frac{\ln(\theta e^{\lambda r})}{\delta \tau})t}, \quad t \ge 0.$$

**Notation 1.** *Denote by λ*<sup>1</sup> *the first positive eigenvalue of the operator* −Δ *in the Sobolev space H*<sup>1</sup> <sup>0</sup> (Ω) *equipped with the norm v* = <sup>Ω</sup> |∇*v*|2*dx for any <sup>v</sup>*(*x*) ∈ *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (Ω). *Denote by E*(*λ*1) *the eigenfunction space of λ*1*. Denote by ϕ*1(*x*) > 0 *the positive eigenfunction corresponding to E*(*λ*1) *with ϕ*1(*x*) = 1. *Besides, I represents the identity matrix. Denote by λ*max(*A*) *the maximum eigenvalue of symmetric matrix A, and by λ*min(*A*) *the minimum eigenvalue of symmetric matrix A.*

#### **3. Main Results**

**Theorem 1.** *Suppose the conditions (A1)–(A3) and (3) hold, and if the following conditions are satisfied:*

$$b\_1 < d\_1 \lambda\_1 \tag{10}$$

$$b\_2 < d\_2 \lambda\_1 \tag{11}$$

*then the system (1) owns multiple stationary solutions, including the positive solution* (*u*∗ <sup>1</sup> (*x*), *u*<sup>∗</sup> <sup>2</sup> (*x*))*T.*

**Proof.** To complete the proof of Theorem 1, the author needs to do it step by step.

*Step 1.* Under the condition (10), there is at least a stationary solution (*α*∗(*x*), 0) for the system (1).

Let (*α*(*x*), 0)*<sup>T</sup>* be a stationary solution of the system (1), satisfying

$$a\_1 \Delta a(\mathbf{x}) + a(\mathbf{x})(b\_1 - a\_{11}a(\mathbf{x})^{\theta\_1} - a\_{12} \cdot 0) = 0, \quad \mathbf{x} \in \Omega; \quad a(\mathbf{x})|\_{\partial \Omega} = 0,\tag{12}$$

whose functional is

$$\Psi(\mathbf{a}) = \frac{1}{2} \int\_{\Omega} |\nabla a(\mathbf{x})|^2 d\mathbf{x} - \frac{b\_1}{2d\_1} \int\_{\Omega} |a(\mathbf{x})|^2 d\mathbf{x} + \frac{a\_{11}}{(2 + \theta\_1)d\_1} \int\_{\Omega} a(\mathbf{x})^{2 + \theta\_1} d\mathbf{x},\tag{13}$$

It is obvious that <sup>Ψ</sup>(0) = 0 and <sup>Ψ</sup> ∈ *<sup>C</sup>*1(*H*<sup>1</sup> <sup>0</sup> (Ω), <sup>R</sup>1).

In fact, for example, when *N* = 3, the assumption (A1) yields that there exists a real number *c* > 0 big enough that

$$|a(\mathbf{x})(b\_1 - a\_{11}\mathbf{a}(\mathbf{x})^{\theta\_1} - a\_{12}\mathbf{\cdot}\mathbf{0})| \leqslant c(1 + |a(\mathbf{x})|^{1+\theta\_1}) \leqslant c(1 + |a(\mathbf{x})|^{2^\*-1}),$$

where 2<sup>∗</sup> = <sup>2</sup>*<sup>N</sup> <sup>N</sup>*−<sup>2</sup> is the Sobolev critical exponent in the case of <sup>Ω</sup> <sup>⊂</sup> *<sup>R</sup>N*. This means <sup>Ψ</sup> ∈ *<sup>C</sup>*1(*H*<sup>1</sup> <sup>0</sup> (Ω), <sup>R</sup>1), and then a critical point of the functional <sup>Ψ</sup> is corresponding to the solution of the Equation (12).

Next, the author claim that Ψ satisfies the (PS) condition.

Indeed, if there exists a real number *<sup>a</sup>* and {*αn*} ⊂ *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (Ω), satisfying |Ψ(*αn*)| *<sup>a</sup>* ∈ *<sup>R</sup>*<sup>1</sup> and Ψ (*αn*) (*H*<sup>1</sup> <sup>0</sup> (Ω))<sup>∗</sup> <sup>→</sup> 0, it means that when *<sup>n</sup>* is big enough,

$$\Psi(a\_n) = \frac{1}{2} \int\_{\Omega} |\nabla a\_n(\mathbf{x})|^2 d\mathbf{x} - \frac{b\_1}{2d\_1} \int\_{\Omega} |a\_n(\mathbf{x})|^2 d\mathbf{x} + \frac{a\_{11}}{(2 + \theta\_1)d\_1} \int\_{\Omega} a\_n(\mathbf{x})^{2 + \theta\_1} d\mathbf{x} \lessgtr a,\tag{14}$$

which together with (A1) and the poincare inequality means

$$\frac{1}{2}(1 - \frac{b\_1}{d\_1 \lambda\_1}) \int\_{\Omega} |\nabla a\_n(\mathbf{x})|^2 d\mathbf{x} \lessapprox a \tag{15}$$

Equation (15) implies the boundedness of {*αn*} in the Sobolev space *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (Ω). Moreover, Hilbert space *H*<sup>1</sup> <sup>0</sup> (Ω) yields that there exists a subsequence, say {*αn*}, such that *αn*(*x*) - *α*(*x*) in *H*<sup>1</sup> <sup>0</sup> (Ω). Rellich Theorem means that *<sup>α</sup>n*(*x*) → *<sup>α</sup>*(*x*) in *<sup>L</sup>p*(Ω) with 1 *<sup>p</sup>* < <sup>2</sup>∗. Therefore,

$$\int\_{\Omega} |\alpha\_{\mathfrak{n}} - \mathfrak{a}|^2 d\mathfrak{x} \to 0, \quad \int\_{\Omega} |\alpha\_{\mathfrak{n}} - \mathfrak{a}|^{2 + \theta\_1} d\mathfrak{x} \to 0, \quad \mathfrak{n} \to \infty.$$

Thus, as *n* → ∞,

$$\|\|a\_{\boldsymbol{n}} - \boldsymbol{a}\|\|^2 = \langle \Psi'(\boldsymbol{a}\_{\boldsymbol{n}}) - \Psi'(\boldsymbol{a}), \boldsymbol{a}\_{\boldsymbol{n}} - \boldsymbol{a} \rangle + \frac{b\_1}{d\_1} \int\_{\Omega} |\boldsymbol{a}\_{\boldsymbol{n}} - \boldsymbol{a}|^2 d\mathbf{x} - \frac{a\_{11}}{d\_1} \int\_{\Omega} |\boldsymbol{a}\_{\boldsymbol{n}} - \boldsymbol{a}|^{2 + \theta\_1} d\mathbf{x} \to 0,$$

which verifies that the (PS) condition is satisfied.

Next, the author claims that there is an upper boundedness for Ψ. In fact, (A1) and (A2) yields

$$\begin{aligned} \Psi(a) &= \frac{1}{2} \int\_{\Omega} |\nabla a(\mathbf{x})|^2 d\mathbf{x} - \frac{b\_1}{2d\_1} \int\_{\Omega} |a(\mathbf{x})|^2 d\mathbf{x} + \frac{a\_{11}}{(2 + \theta\_1)d\_1} \int\_{\Omega} a(\mathbf{x})^{2 + \theta\_1} d\mathbf{x} \\ &\le \frac{1}{2} \int\_{\Omega} |\nabla a(\mathbf{x})|^2 d\mathbf{x} + \frac{a\_{11}}{(2 + \theta\_1)d\_1} M\_1^{2 + \theta\_1} m \text{es}(\Omega), \end{aligned}$$

which together with (A3) means that there exists an upper boundedness for Ψ.

According to Lemma 1, there exists *α*∗(*x*) such that

$$J(\mathfrak{a}\_\*(x)) = \sup\_{\upsilon \in H\_0^1(\Omega)} J(\upsilon)$$

and (*α*∗(*x*), 0)*<sup>T</sup>* is a stationary solution of the system (1).

*Step 2.* The author claims that the system (1) owns multiple stationary solutions, including the positive solution.

First, the condition (3) and Lemma 2 guarantee the existence of a positive stationary solution for the system (1). Second, zero solution (0, 0)*<sup>T</sup>* is obviously another stationary solution for the system (1). Next, (*α*∗(*x*), 0)*<sup>T</sup>* is the third stationary solution thanks to Step 1. In fact, the continuity of *ϕ*1(*x*) yields

$$J(a\_\*(x)) = \sup\_{v \in H\_0^1(\Omega)} J(v) \gg \quad J(\varphi\_1) \gg \frac{a\_{11}}{(2 + \theta\_1)d\_1} \int\_{\Omega} \varphi\_1(x)^{2 + \theta\_1} dx > 0,$$

which means that (*α*∗(*x*), 0)*<sup>T</sup>* is a nontrivial stationary solution for the system (1). Finally, one can similarly prove that there exists a nontrivial stationary solution (0, *<sup>β</sup>*∗(*x*))*<sup>T</sup>* for the system (1).

**Theorem 2.** *Suppose that all the conditions of Theorem 1 are satisfied. Assume, in addition,*

*(B1) there exist three positive constants pm*, *pM, ε, and a positive definite diagonal matrix P* = *diag*(*p*1, *p*1) > 0 *such that the following LMI conditions hold:*

$$2\lambda\_1 PD - 2P(B - K) - p\_M \Theta - \varepsilon PK > 0\tag{16}$$

$$P < p\_M I \tag{17}$$

$$
\bar{p}\_M I < P \tag{18}
$$

*where*

$$\Theta = \begin{pmatrix} 2\left(a\_{11}(1+\theta\_1)(2M\_1)^{\theta\_1} + a\_{12}M\_2\right) & a\_{12}M\_1 + a\_{21}M\_2\\ \ast & 2\left(a\_{22}(1+\theta\_2)(2M\_2)^{\theta\_2} + a\_{21}M\_1\right) \end{pmatrix}.$$

$$\text{(B2) } a > b \gtrless 0,\\
\text{where } a = \frac{\lambda\_{\text{min}} \left(2\lambda\_1 PD - 2P(B-K) - p\_M \Theta - \varepsilon PK\right)}{p\_M},\\
b = \frac{\lambda\_{\text{max}}(K)}{\varepsilon}$$

*(B3) there exists a constant <sup>δ</sup>* <sup>&</sup>gt; <sup>1</sup> *such that* inf*k*∈*Z*(*tk* <sup>−</sup> *tk*−1) <sup>&</sup>gt; *δτ and <sup>λ</sup>* <sup>&</sup>gt; ln(*ρeλτ* ) *δτ , where ρ* = sup *<sup>j</sup>*∈<sup>Z</sup> {1, *aj* <sup>+</sup> *bjeλτ*} *with aj* <sup>=</sup> *<sup>λ</sup>*max(*A<sup>T</sup> <sup>j</sup> PAj*) *pm and bj* ≡ 0, *and λ* > 0 *is the unique solution of the equation <sup>λ</sup>* <sup>=</sup> *<sup>a</sup>* <sup>−</sup> *beλτ.*

*then the zero solution of the system (8) is globally exponentially stable with convergence rate* 1 <sup>2</sup> (*<sup>λ</sup>* <sup>−</sup> ln(*ρeλτ* ) *δτ* )*, and* (*u*<sup>∗</sup> <sup>1</sup> (*x*), *u*<sup>∗</sup> <sup>2</sup> (*x*))*<sup>T</sup> is said to be globally exponentially stable under impulsive disturbances with convergence rate* <sup>1</sup> <sup>2</sup> (*<sup>λ</sup>* <sup>−</sup> ln(*ρeλτ* ) *δτ* ).

**Proof.** Consider the following Lyapunov function:

$$V(t) = \int\_{\Omega} \mathcal{U}^T(t, \mathbf{x}) P \mathcal{U}(t, \mathbf{x}) d\mathbf{x} = \int\_{\Omega} |\mathcal{U}(t, \mathbf{x})|^T P |\mathcal{U}(t, \mathbf{x})| d\mathbf{x}$$

then for *t* -0, *t* = *tk*, the Poincare inequality yields

$$\begin{split} D^{+}V &= 2 \int\_{\Omega} \mathcal{U}^{T} P \Big( D\Delta \mathcal{U} + (B - K)\mathcal{U} - \Phi(\mathcal{U}) + K\mathcal{U}(t - \tau, \mathbf{x}) \Big) d\mathbf{x} \\ &\leqslant \int\_{\Omega} \mathcal{U}^{T} P \Big( -2\lambda\_{1}D + 2(B - K) \Big) \mathcal{U} d\mathbf{x} + \int\_{\Omega} \Big( -2\mathcal{U}^{T} P \Phi(\mathcal{U}) + 2\mathcal{U}^{T} P K \mathcal{U}(t - \tau, \mathbf{x}) \Big) d\mathbf{x} \\ &\leqslant \int\_{\Omega} |\mathcal{U}|^{T} P \Big( -2\lambda\_{1}D + 2(B - K) \Big) |\mathcal{U}| d\mathbf{x} + \int\_{\Omega} \Big( 2|\mathcal{U}|^{T} P |\Phi(\mathcal{U})| + 2|\mathcal{U}|^{T} P K |\mathcal{U}(t - \tau, \mathbf{x})| \Big) d\mathbf{x} \end{split} \tag{19}$$

On the other hand, it follows from (5) that Φ1(0, 0) = 0 = Φ2(0, 0), and

$$\Phi\_1(0, \mathsf{U}2) = a\_{12}u\_1^\*(\mathbf{x})(\mathsf{U}2 + u\_2^\*(\mathbf{x})) - a\_{12}u\_1^\*(\mathbf{x})u\_2^\*(\mathbf{x}) = a\_{12}u\_1^\*(\mathbf{x})\mathsf{U}2 \tag{20}$$

and thus differential mean value theorem and (A2) yield

$$\begin{split} |\Phi\_{1}(\mathcal{U})| &= |\Phi\_{1}(\mathcal{U}) - \Phi\_{1}(0)| \leqslant |\Phi\_{1}(\mathcal{U}\_{1}, \mathcal{U}\_{2}) - \Phi\_{1}(0, \mathcal{U}\_{2})| + |\Phi\_{1}(0, \mathcal{U}\_{2}) - \Phi\_{1}(0, 0)| \\ &\leqslant \Big( a\_{11}(1 + \theta\_{1})(2M\_{1})^{\theta\_{1}} + a\_{12}M\_{2} \Big) |\mathcal{U}\_{1}| + a\_{12}M\_{1} |\mathcal{U}\_{2}|. \end{split} \tag{21}$$

Similarly,

$$|\Phi\_2(\mathcal{U})| \lesssim a\_{21} M\_2 |\mathcal{U}\_1| + \left( a\_{22} (1 + \theta\_2)(2M\_2)^{\theta\_2} + a\_{21} M\_1 \right) |\mathcal{U}\_2| \tag{22}$$

Thus,

$$\begin{split} 2|\mathcal{U}|^{T}P|\Phi(\mathcal{U})| &\leqslant p\_{M}(2|\mathcal{U}\_{1}|\cdot|\Phi\_{1}(\mathcal{U})| + 2|\mathcal{U}\_{2}|\cdot|\Phi\_{2}(\mathcal{U})|) \\ &\leqslant p\_{M}|\mathcal{U}|^{T}\Theta|\mathcal{U}| \end{split} \tag{23}$$

$$2|\mathcal{U}|^T PK|\mathcal{U}(t-\tau,\mathbf{x})| \lessapprox \varepsilon |\mathcal{U}|^T (PK)|\mathcal{U}| + \frac{1}{\varepsilon} \lambda\_{\text{max}}(K) |\mathcal{U}(t-\tau,\mathbf{x})|^T P |\mathcal{U}(t-\tau,\mathbf{x})| \tag{24}$$

Combining (19)–(24) results in

$$\begin{split} \label{eq:1} D^{+}V(t) &\leqslant \int\_{\Omega} |\mathcal{U}|^{T} P \Big( -2\lambda\_{1}D + 2(B-K) \Big) |\mathcal{U}| d\mathbf{x} + \int\_{\Omega} \Big( 2|\mathcal{U}|^{T} P |\Phi(\mathcal{U})| + 2|\mathcal{U}|^{T} P K |\mathcal{U}(t-\tau, \mathbf{x})| \Big) d\mathbf{x} \\ &\leqslant -\frac{\lambda\_{\text{min}} \Big( 2\lambda\_{1}PD - 2P(B-K) - p\_{M}\Theta - \varepsilon PK \Big)}{p\_{M}} \int\_{\Omega} |\mathcal{U}|^{T} P |\mathcal{U}| d\mathbf{x} + \frac{\lambda\_{\text{max}}(K)}{\varepsilon} V(t-\tau) \\ &\leqslant -av(t) + b[v(t)]\_{\tau}, t \neq t\_{k}. \end{split} \tag{25}$$

On the other hand, letting *<sup>γ</sup>* <sup>=</sup> ln(*ρeλτ* ) *δτ* , one can conclude from Lemma 3 that

$$V(t) \ll (\rho^2 e^{\lambda \tau}) [V(0)] \, \_1e^{-(\lambda - \gamma)t}, \quad t \gg t\_{0\prime} \tag{26}$$

or equivalently,

$$V(t) \ll (\rho^2 e^{\lambda \tau}) [V(0)]\_{\tau} e^{-(\lambda - \frac{\ln(\rho e^{\lambda \tau})}{\delta \tau})t}, \quad t \gg t\_{0\prime} \tag{27}$$

Indeed,

$$\begin{split} V(t\_k) &= \int\_{\Omega} \mathcal{U}^T(t\_k, \mathbf{x}) P \mathcal{U}(t\_k, \mathbf{x}) d\mathbf{x} \\ &\leqslant \frac{\lambda\_{\text{max}}(A\_k^T P A\_k)}{p\_{\text{m}}} \int\_{\Omega} \mathcal{U}^T(t\_k^-, \mathbf{x}) P \mathcal{U}(t\_k^-, \mathbf{x}) d\mathbf{x} \\ &= a\_k V(t\_k^-). \end{split}$$

According the conditions (B1)–(B3), one can see it from Lemma 3 that (26) and (27) holds if the condition (9) is verified. In fact, in Lemma 3, let *M* = *ρ*2*eλτ*, then

$$\begin{aligned} Me^{\gamma t\_k} &= (\rho^2 e^{\lambda \tau}) e^{\gamma (t\_k - t\_0)} \\ &\geqslant (\rho^2 e^{\lambda \tau}) (\rho e^{\lambda \tau})^{k-1} \\ &= (\rho^{k+1} e^{k\lambda \tau})\_{\prime} \end{aligned}$$

which means that the condition (9) is satisfied, and then Lemma 3 makes (26) and (27) hold. Moreover, (27) yields

$$\begin{split} \|p\_{\mathfrak{m}}\| \|\mathcal{U}\|\_{L^{2}(\Omega)} \leqslant & V(t) \leqslant (\rho^{2} e^{\lambda \tau}) [V(0)]\_{\mathfrak{r}} e^{-\left(\lambda - \frac{\ln(\rho e^{\lambda \tau})}{\delta \overline{\tau}}\right) t} \\ \leqslant & (\rho^{2} e^{\lambda \tau}) p\_{M} \|\zeta(\mathsf{s}, \mathsf{x}) - \mathsf{u}^{\*}(\mathsf{x})\|\_{\mathfrak{r}}^{2} e^{-\left(\lambda - \frac{\ln(\rho e^{\lambda \tau})}{\delta \overline{\tau}}\right) t}, \quad t \geqslant t\_{0}. \end{split} \tag{28}$$

where *<sup>ξ</sup>*(*s*, *<sup>x</sup>*) − *<sup>u</sup>*∗(*x*) <sup>2</sup> *<sup>τ</sup>* = sup *s*∈[−*τ*,0] <sup>Ω</sup>[*ξ*(*s*, *<sup>x</sup>*) − *<sup>u</sup>*∗(*x*)]*T*[*ξ*(*s*, *<sup>x</sup>*) − *<sup>u</sup>*∗(*x*)]*dx* with *<sup>ξ</sup>* = (*ξ*1, *<sup>ξ</sup>*2)*<sup>T</sup>* and *u*∗ = (*u*∗ <sup>1</sup>, *u*<sup>∗</sup> <sup>2</sup> )*T*. Obviously, (28) completes the proof.

**Remark 2.** *Theorem 2 offers a better stabilization criterion than the previous literature ([6,7]), which reduces the conservatism of the algorithm. In fact, in Theorem 2, the impulse condition λ*min*Ak may not be smaller than 1, which implies that this paper deletes the harsh restrictions on small impulse of the related literature [6,7].*

**Remark 3.** *Compared with the previous related literature [6,7], Theorem 2 expands the range of the parameters θ*1, *θ*<sup>2</sup> *from* (0, 1) *to* (−1, 4)*.*

#### **4. Numerical Examples**

First, the following example shows the effectiveness of Theorem 1.

**Example 1.** *Let θ*<sup>1</sup> = <sup>2</sup> <sup>3</sup> , *<sup>θ</sup>*<sup>2</sup> <sup>=</sup> <sup>4</sup> <sup>5</sup> *, bi* = 0.13 + 0.0001*i*, *di* = 0.1 + 0.0001*i*, *i* = 1, 2, *and* Ω = (0, 1) × (0, 1)*. Direct calculation yields that λ*<sup>1</sup> = 19.7392 *([14] (Remark 14)), and b*<sup>1</sup> = 0.1301 < 0.1001 × 19.7392 = *d*1*λ*<sup>1</sup> *and b*<sup>2</sup> = 0.1302 < 0.1002 × 19.7392 = *d*2*λ*1*. Furthermore, set Mi* = 2 + 0.1*i*, *i* = 1, 2, *and a*<sup>11</sup> = 0.03, *a*<sup>12</sup> = 0.02, *a*<sup>21</sup> = 0.025, *a*<sup>22</sup> = 0.03*, k*<sup>1</sup> = 0.15, *k*<sup>2</sup> = 0.12*. An accurate calculation can verify that the condition (3) is satisfied if letting c*<sup>∗</sup> = 100000*. Now, one can conclude from Theorem 1 that there is a positive stationary solution* (*u*∗ <sup>1</sup> (*x*), *u*<sup>∗</sup> <sup>2</sup> (*x*))*<sup>T</sup> and other three stationary solutions for the ecosystem (1).*

Below, the feasibility of Theorem 2 need be verified, too.

**Example 2.** *All the data of Example 1 are employed in this example, then an accurate calculation yields that*

$$
\Theta = \begin{pmatrix}
0.3483 & 0.0970 \\
0.0970 & 0.4583
\end{pmatrix}
$$

*(B2) a* > *b* - 0*. Furthermore, using computer Matlab LMI toolbox to solve LMI condition (16)–(18) yields the following feasible data:*

$$P = \begin{pmatrix} 0.9998 & 0\\ 0 & 1.0013 \end{pmatrix}, \varepsilon = 0.9996, \ p\_M = 1.0015, \ p\_{\mathfrak{m}} = 0.9973.$$

*Then, a direct calculation obtains a* = 3.3046, *b* = 0.1501, *and thus a* > *b* - 0. *Let τ* = 0.5, *solving the equality <sup>λ</sup>* <sup>=</sup> *<sup>a</sup>* <sup>−</sup> *beλτ reaches <sup>λ</sup>* <sup>=</sup> 2.7199*. Set*

$$A\_j \equiv \begin{pmatrix} 1.0603 & 0 \\ 0 & 1.0783 \end{pmatrix}, \forall j \in \mathbb{Z}, \tag{29}$$

*which together the above data derives that aj* ≡ 1.1674*, and thus ρ* = 1.1674. *Set δ* = 2, *then an immediate calculation yields <sup>λ</sup>* <sup>−</sup> ln(*ρeλτ* ) *δτ* = 1.2052 > <sup>0</sup>*, and* <sup>1</sup> <sup>2</sup> [*<sup>λ</sup>* <sup>−</sup> ln(*ρeλτ* ) *δτ* ] = 0.6026. *According to Theorem 2, the zero solution of the system (8) is globally exponentially stable with convergence rate* 60.26%*.*

**Remark 4.** *Example 2 verifies the advantages described in Remarks 2–3.*

#### **5. Conclusions and Further Considerations**

Compared with the known literature, this paper has double advantages in method and conclusion. On one hand, employing the Minimax principle and impulsive differential inequality improves the methods in [6,7]. For example, in deriving the existence of multiple stationary solutions of RDGACM, the methods involved in Minimax principle is more simpler than those in applying Mountain Pass Lemma of [6]. Besides, in stabilizing globally the ecosystem, utilizing the impulsive differential inequality makes the impulse range wider. Especially, an impulse range means that people can adjust and manage the ecosystem more flexibly.

For *<sup>v</sup>* ∈ *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (Ω), the norm

$$\|v\| = \sqrt{\int\_{\Omega} |\nabla v|^2 dx} \tag{30}$$

in this paper is simpler than the norm

$$\|v\| = \sqrt{\int\_{\Omega} (|\nabla v|^2 + \frac{\mathbb{C}}{D} v^2) dx} \tag{31}$$

in [18] (Statement 2). Now, with the help of such a simple norm (30), some further considerations are presented below.

In fact, in [18] (Statement 2), the following ordinary differential equation and its corresponding partial differential equation were considered:

$$\frac{d\mathbf{x}(t)}{dt} = -\mathbb{C}\mathbf{x}(t) + A f(\mathbf{x}(t)) + B f(\mathbf{x}(t - \mathbf{r}(t))) + f,\qquad \text{and } \mathbf{x} \in \mathbb{R}^1,\tag{32}$$

and its corresponding reaction-diffusion cellular neural networks

$$\begin{cases} \frac{\partial u(t, \mathbf{x})}{\partial t} = \mathcal{D}\Lambda u(t, \mathbf{x}) - \mathbb{C}u(t, \mathbf{x}) + A f(u(t, \mathbf{x})) + \mathcal{B}f(u(t - \tau(t), \mathbf{x})) + \mathbb{I}, \quad \text{and } t \gg 0, \mathbf{x} \in \Omega, \\\ u(t, \mathbf{x}) = 0, \qquad \mathbf{x} \in \partial\Omega, \end{cases} \tag{33}$$

where <sup>Ω</sup> is an open bounded domain in <sup>R</sup><sup>3</sup> with smooth boundary *<sup>∂</sup>*Ω, *<sup>D</sup>* <sup>∈</sup> <sup>R</sup><sup>1</sup> is the diffusion coefficient with *D* > 0, and *C*, *A* both are positive real numbers, *J* = 0, *B* = 0, the function *f* is defined as follows:

$$f(u) = \begin{cases} \frac{3D}{A} \mu\_1 u^{\frac{1}{3}} + \frac{2D}{A} \mu\_{1\prime} & u \leqslant -1; \\ \frac{D}{A} \mu\_1 u\_{\prime} & u \in [-1, 1]; \\ \frac{3D}{A} \mu\_1 u^{\frac{1}{3}} - \frac{2D}{A} \mu\_{1\prime} & u \geqslant 1. \end{cases} \tag{34}$$

Here, we denote by *μ<sup>i</sup>* the *i*th positive eigenvalue of the following eigenvalue problem:

$$\begin{cases} -\Delta u(\mathbf{x}) + \frac{\mathcal{C}}{D} u(\mathbf{x}) = \mu u(\mathbf{x}), \mathbf{x} \in \Omega, \\ \qquad u(\mathbf{x}) = 0, \qquad u \in \partial \Omega, \end{cases} \tag{35}$$

Particularly in the case of *C* = 0, the norm (31) is just that of (30), and then *μ*<sup>1</sup> = *λ*<sup>1</sup> > 0 is the first positive eigenvalue of the operator −<sup>Δ</sup> in *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (Ω) with the norm (30). Thus, in the case of *C* = 0, the following theorem holds:

**Theorem 3.** *If zero solution is the global stable unique equilibrium point of the following ordinary differential system*

$$\frac{d\mathbf{x}(t)}{dt} = A f(\mathbf{x}(t)) + B f(\mathbf{x}(t - \tau(t))), \qquad \text{and } \mathbf{x} \in \mathbb{R}^1,\tag{36}$$

*where*

$$f(u) = \begin{cases} \frac{3D}{A} \lambda\_1 u^{\frac{1}{3}} + \frac{2D}{A} \lambda\_1 & u \leqslant -1; \\ \frac{D}{A} \lambda\_1 u, & u \in [-1, 1]; \\ \frac{3D}{A} \lambda\_1 u^{\frac{1}{3}} - \frac{2D}{A} \lambda\_1, & u \geqslant 1, \end{cases} \tag{37}$$

*then its corresponding reaction-diffusion system:*

$$\begin{cases} \frac{\partial u(t, \mathbf{x})}{\partial t} = \mathcal{D} \Delta u(t, \mathbf{x}) + A f(u(t, \mathbf{x})) + B f(u(t - \tau(t), \mathbf{x})), \quad \text{and} \ t \geqslant 0, \mathbf{x} \in \Omega, \\\ u(t, \mathbf{x}) = 0, \qquad \mathbf{x} \in \partial \Omega, \end{cases} \tag{38}$$

*owns zero solution and other stationary solutions which are at least two non-zero functions or infinitely many positive functions and negative functions.*

**Remark 5.** *It follows from [18](Remark 11) that zero solution is actually the global stable unique equilibrium point of the ordinary differential system (36). In fact, according to the Introduction in [20], the function f defined by (37) satisfies the conditions [20] (Equation (7)) and [20] (Equation (8)), and thus the zero solution is actually the unique equilibrium point of the ordinary differential system (36). That is, there exists such an example that under the influence of diffusion, the unique equilibrium point of the ordinary differential system (36) with the Lipschitz activation function f can become at least three equilibrium points of its corresponding reaction-diffusion system (38).*

Now, in view of Theorem 3 and Remark 5, the author wants to know whether an example can be designed such that the global stable unique equilibrium point *x*∗ of the ordinary differential system can become multiple equilibrium points *u*∗ *<sup>i</sup>* (*x*)(*i* ∈ Λ) of its corresponding reaction-diffusion system under the influence of diffusion? Here, Λ is a finite index set or infinite index set. Furthermore, is the diffusion coefficient related to the number of the index set Λ? Is the smaller the diffusion coefficient, the fewer the number of the index set Λ? Moreover, if the diffusion coefficient is small enough, is the norm *u*<sup>∗</sup> *<sup>i</sup>* (*x*) − *x*∗ ∗ is also small? Here, · ∗ may be *u*<sup>∗</sup> *<sup>i</sup>* (*x*) − *x*∗ ∗ = sup *x*∈Ω |*u*∗ *<sup>i</sup>* (*x*) − *x*∗|. All

these problems are interesting.

**Funding:** The research was supported by the Application Basic Research Project of Science and Technology Department of Sichuan Province in China (No. 2020YJ0434).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**

