# Stability of Nonlinear Autonomous Quadratic Discrete Systems in the Critical Case

Authors: Josef Diblík; Denys Ya. Khusainov; Irina V. Grytsay; Zdenĕk Šmarda
Journal: Discrete Dynamics in Nature and Society
 (2010)

---

## Abstract

Many processes are mathematically simulated by systems of discrete equations with quadratic
right-hand sides. Their stability is thought of as a very important characterization of the process. In this paper, the method of Lyapunov functions is used to derive classes of stable quadratic discrete autonomous systems in a critical case in the presence of a simple eigenvalueλ=1 of the matrix of linear terms. In addition to the stability investigation, we
also estimate stability domains.


---

## Body

## 1. Introduction

The main results on the stability theory of difference equations are presented, for example, by Agarwal [1], Agarwal et al. [2], Chetaev [3], Elaydi [4], Halanay and Răsvan [5], Lakshmikantham and Trigiante [6], and Martynjuk [7]. Instability problems are considered, for example, in [8–10] by Slyusarchuk. Note that stability and instability results often have a local character and are usually obtained without any estimation of the stability domain, or without investigating the character of instability. Moreover, it should be emphasized that global instability questions have only been discussed for linear systems.Many processes and phenomena are described by differential or difference systems with quadratic nonlinearities. Among others, let us mention epidemic and populations models, models of chemical reactions, and models for describing convection currents in the atmosphere.The stability of a zero solution of difference systems(1.1)x(k+1)=f(x(k)),
where k=0,1,…, and x=(x1,…,xn)T with differentiable f=(f1,…,fn)T:ℝn→ℝn, is very often investigated by linearly approximating system (1.1) in question by using the matrix of linear terms(1.2)x(k+1)=Ax(k)+g(x(k)),
where A=f′(0,0,…,0) is the Jacobian matrix of f at (0,0,…,0), and g(x)=f(x)-Ax. This approach becomes unsuitable in what is called a critical case, that is, when the spectral radius of the matrix ρ(A)=1 because, among all systems (1.2), there are classes of stable systems as well as classes of unstable systems. Concerning this, we formulate the following known results (see, e.g., Corollary 4.34 [4, page 222] and Theorem 4.38 [4, page 226]).Theorem 1.1.
(1) If ρ(A)<1, then the zero solution of (1.2) is exponentially stable.
(2) If ρ(A)=1, then the zero solution of (1.2) may be stable or unstable.
(3) If ρ(A)>1 and g(x) is o(x) as ∥x∥→0, then the zero solution of (1.2) is unstable.In this paper, we consider a particular critical case when there exists a simple eigenvalueλ=1 of the matrix of linear terms and the remaining eigenvalues lie inside a unit circle centered at origin. The purpose of this paper is to obtain (using the method of Lyapunov functions) conditions for the stability of a zero solution of difference systems with quadratic nonlinearities in the above case and derive classes of stable systems. In addition to the stability investigation, we estimate the stability domains as well. The domains of stability obtained are also called guaranteed domains of stability. Preliminary results in this direction were published in [11].


### 1.1. Quadratic System and Preliminary Consideration

In the sequel, the norms used for vectors and matrices are defined as(1.3)∥x∥=(∑i=1nxi2)1/2
for a vector x=(x1,…,xn)T and (1.4)∥A∥=(λmax(ATA))1/2
for any m×n matrix A. Here and in the sequel, λmax(·) (or λmin(·)) is the maximal (or minimal) eigenvalue of the corresponding symmetric and positive- (semi-) definite matrix (see, e.g., [12]).Consider a nonlinear autonomous discrete system with a quadratic right-hand side(1.5)xi(k+1)=∑s=1naisxs(k)+∑s,q=1nbsqixs(k)xq(k),i=1…,n,
where k=0,1,… and the coefficients ais and bsqi (we assume that bsqi=bqsi) are constant. As emphasized, for example, in [3, 7, 12], system (1.5) can be written in a general vector-matrix form(1.6)x(k+1)=Ax(k)+XT(k)Bx(k),k=0,1,…,
where  (a)
A={ais}, i,s=1,2…,n, is an n×n constant square matrix,(b)
matrixXT={X1T,X2T,…,XnT} is n×n2 rectangular and all the elements of the n×n matrices XiT, i=1,…,n, are equal to zero except the ith row with entries xT=(x1,x2,…,xn), that is,
(1.7)XiT(k)=(00⋯0⋯⋯⋯⋯x1x2⋯xn⋯⋯⋯⋯00⋯0),(c)
matrixBT={B1,B2,…,Bn} is n2×n rectangular and the n×n constant matrices Bi={bsqi}, i,s,q=1,…,n, are symmetric.The stability of the zero solution of (1.6) depends on the stability of the matrix A. If ρ(A)<1, then the zero solution of (1.6) is exponentially stable for an arbitrary matrix B by Theorem 1.1. In this case, matrix B only impacts on the shape of the stability domain of the equilibrium state. If the zero solution of (1.6) is investigated on stability by the second Lyapunov method and an appropriate Lyapunov function is taken as the quadratic form V(x)=xTHx with a suitable n×n constant real symmetric positive-definite matrix H, which is defined below, then the first difference ΔV along the trajectories of (1.6) equals(1.8)ΔV(x(k))=V(x(k+1))-V(x(k))=xT(k+1)Hx(k+1)-xT(k)Hx(k)=[Ax(k)+XT(k)Bx(k)]TH[Ax(k)+XT(k)Bx(k)]-xT(k)Hx(k)=[xT(k)AT+xT(k)BTX(k)]H[Ax(k)+XT(k)Bx(k)]-xT(k)Hx(k)=xT(k)[(ATHA-H)+ATHXT(k)B+BTX(k)HA+BTX(k)HXT(k)B]x(k)=xT(k)[(ATHA-H)+2BTX(k)HA+BTX(k)HXT(k)B]x(k)
since (ATHXT(k)B)T=BTX(k)HA.Sinceρ(A)<1, for arbitrary positive-definite symmetric matrix C, the matrix Lyapunov equation (1.9)ATHA-H=-C
has a unique solution H—a positive-definite symmetric matrix (e.g., [4, Theorem 4.30, page 216]). We use such matrix H to estimate the stability domain. Then, as follows from (1.8),(1.10)ΔV(x(k))≤-[λmin(C)-2∥B∥·∥HA∥·∥x(k)∥-∥B∥2·∥H∥·∥x(k)∥2]·∥x(k)∥2.
Analysing (1.10), we deduce that the first difference ΔV(x(k)) will be negative definite if(1.11)∥B∥2·∥H∥·∥x(k)∥2+2∥B∥·∥HA∥·∥x(k)∥≤λmin(C),
that is, it will be negative definite in a neighborhood Uδ={x∈ℝn:∥x∥<δ} of the steady-state x(k)≡0, k=0,1,…, if δ is sufficiently small. In the case considered, the domain of stability can be described by means of two inequalities. The first inequality (1.11) defines a part of the space ℝn, where the first difference ΔV(x(k)) is negative definite. The second inequality(1.12)V(x)≤r,x∈ℝn,r>0,
describes points inside a level surface. The guaranteed domain of stability is given by inequality (1.12) if r is taken so small that the domain described by (1.12) is embedded in the domain described by inequality (1.11).Considering the investigated critical case, we will deal with a different structure of the right-hand side of the inequality from that seen in (1.10). Namely, we will show that, unlike the right-hand side of the inequality for ΔV(x(k)) that is multiplied by ∥x(k)∥2 with dimx(k)=n in (1.10), in the critical case considered, the right-hand side of the inequality (or equality) for ΔV(x(k)) will be multiplied only by a term ∥xn-1(k)∥2 with dimxn-1(k)=n-1<n (see (2.21) in the case n=2 and (2.69) in the general case below).


## 1.1. Quadratic System and Preliminary Consideration

In the sequel, the norms used for vectors and matrices are defined as(1.3)∥x∥=(∑i=1nxi2)1/2
for a vector x=(x1,…,xn)T and (1.4)∥A∥=(λmax(ATA))1/2
for any m×n matrix A. Here and in the sequel, λmax(·) (or λmin(·)) is the maximal (or minimal) eigenvalue of the corresponding symmetric and positive- (semi-) definite matrix (see, e.g., [12]).Consider a nonlinear autonomous discrete system with a quadratic right-hand side(1.5)xi(k+1)=∑s=1naisxs(k)+∑s,q=1nbsqixs(k)xq(k),i=1…,n,
where k=0,1,… and the coefficients ais and bsqi (we assume that bsqi=bqsi) are constant. As emphasized, for example, in [3, 7, 12], system (1.5) can be written in a general vector-matrix form(1.6)x(k+1)=Ax(k)+XT(k)Bx(k),k=0,1,…,
where  (a)
A={ais}, i,s=1,2…,n, is an n×n constant square matrix,(b)
matrixXT={X1T,X2T,…,XnT} is n×n2 rectangular and all the elements of the n×n matrices XiT, i=1,…,n, are equal to zero except the ith row with entries xT=(x1,x2,…,xn), that is,
(1.7)XiT(k)=(00⋯0⋯⋯⋯⋯x1x2⋯xn⋯⋯⋯⋯00⋯0),(c)
matrixBT={B1,B2,…,Bn} is n2×n rectangular and the n×n constant matrices Bi={bsqi}, i,s,q=1,…,n, are symmetric.The stability of the zero solution of (1.6) depends on the stability of the matrix A. If ρ(A)<1, then the zero solution of (1.6) is exponentially stable for an arbitrary matrix B by Theorem 1.1. In this case, matrix B only impacts on the shape of the stability domain of the equilibrium state. If the zero solution of (1.6) is investigated on stability by the second Lyapunov method and an appropriate Lyapunov function is taken as the quadratic form V(x)=xTHx with a suitable n×n constant real symmetric positive-definite matrix H, which is defined below, then the first difference ΔV along the trajectories of (1.6) equals(1.8)ΔV(x(k))=V(x(k+1))-V(x(k))=xT(k+1)Hx(k+1)-xT(k)Hx(k)=[Ax(k)+XT(k)Bx(k)]TH[Ax(k)+XT(k)Bx(k)]-xT(k)Hx(k)=[xT(k)AT+xT(k)BTX(k)]H[Ax(k)+XT(k)Bx(k)]-xT(k)Hx(k)=xT(k)[(ATHA-H)+ATHXT(k)B+BTX(k)HA+BTX(k)HXT(k)B]x(k)=xT(k)[(ATHA-H)+2BTX(k)HA+BTX(k)HXT(k)B]x(k)
since (ATHXT(k)B)T=BTX(k)HA.Sinceρ(A)<1, for arbitrary positive-definite symmetric matrix C, the matrix Lyapunov equation (1.9)ATHA-H=-C
has a unique solution H—a positive-definite symmetric matrix (e.g., [4, Theorem 4.30, page 216]). We use such matrix H to estimate the stability domain. Then, as follows from (1.8),(1.10)ΔV(x(k))≤-[λmin(C)-2∥B∥·∥HA∥·∥x(k)∥-∥B∥2·∥H∥·∥x(k)∥2]·∥x(k)∥2.
Analysing (1.10), we deduce that the first difference ΔV(x(k)) will be negative definite if(1.11)∥B∥2·∥H∥·∥x(k)∥2+2∥B∥·∥HA∥·∥x(k)∥≤λmin(C),
that is, it will be negative definite in a neighborhood Uδ={x∈ℝn:∥x∥<δ} of the steady-state x(k)≡0, k=0,1,…, if δ is sufficiently small. In the case considered, the domain of stability can be described by means of two inequalities. The first inequality (1.11) defines a part of the space ℝn, where the first difference ΔV(x(k)) is negative definite. The second inequality(1.12)V(x)≤r,x∈ℝn,r>0,
describes points inside a level surface. The guaranteed domain of stability is given by inequality (1.12) if r is taken so small that the domain described by (1.12) is embedded in the domain described by inequality (1.11).Considering the investigated critical case, we will deal with a different structure of the right-hand side of the inequality from that seen in (1.10). Namely, we will show that, unlike the right-hand side of the inequality for ΔV(x(k)) that is multiplied by ∥x(k)∥2 with dimx(k)=n in (1.10), in the critical case considered, the right-hand side of the inequality (or equality) for ΔV(x(k)) will be multiplied only by a term ∥xn-1(k)∥2 with dimxn-1(k)=n-1<n (see (2.21) in the case n=2 and (2.69) in the general case below).


## 2. Main Results

In this section we derive the classes of the stable systems (1.6) in a critical case when the matrix A has one simple eigenvalue λ=1.


### 2.1. Instability in One-Dimensional Case

We start by discussing a simple scalar equation with the eigenvalue of matrixA equaling one, that is, a11=1. Then (1.6) takes the form (2.1)x(k+1)=x(k)+bx2(k),k=0,…,
and it is easy to see that the trivial solution is unstable for an arbitrary b≠0 (to show this, we can apply, e.g., Theorem 1.15 [4, page 29]).This elementary example shows that stability in the case of system (1.6) has an extraordinary significance and the results on stability (for n≠1) lose their meaning for n=1 when we deal with instability. We show that, if n≠1 and B satisfies certain assumptions, the zero solution is stable. Moreover, the shape of the guaranteed domain of stability will be given.We divide our forthcoming analysis into two parts. In the first one we give an explicit coefficient criterion in the subcase ofn=2. Then we consider the general n-dimensional case.


### 2.2. Stability in the General Two-Dimensional Case

Letn=2. Then system (1.6) with the matrix A having a simple eigenvalue λ=1 reduces (after linearly transforming the dependent variables if necessary) to(2.2)x1(k+1)=ax1(k)+[b111x12(k)+2b121x1(k)x2(k)+b221x22(k)],x2(k+1)=x2(k)+[b112x12(k)+2b122x1(k)x2(k)+b222x22(k)].
We will assume that |a|<1. Define auxiliary numbers as follows: (2.3)α=h(1-a2),β1=hab111,β2=2hab121+b112,γ1=h(b111)2+(b112)2,γ2=4h(b121)2,δ1=2hb111b121,
where h is a positive number.Theorem 2.1.
Leth and r be positive numbers. Assume that |a|<1 and b122=b221=b222=0. Then the zero solution of system (2.2) is stable in the Lyapunov sense and a guaranteed domain of stability is given by the inequality
(2.4)hx12+x22≤r2
if r is taken so small that the domain described by (2.4) is embedded in the domain
(2.5)γ1x12+2δ1x1x2+γ2x22+2β1x1+2β2x2≤α.
If, moreover, b112·b121≠0, then a guaranteed domain of stability can be described using inequality
(2.6)hx12+x22≤(r*)2
with
(2.7)r*=min(x1,x2)hx12+x22,
where (x1,x2) runs over all real solutions of the nonlinear system with unknowns x1 and x2:
(2.8)γ1x12+2δ1x1x2+γ2x22+2β1x1+2β2x2=α,hx1(δ1x1+γ2x2+β2)-x2(γ1x1+δ1x2+β1)=0.Proof.
Define(2.9)B1=(b111b121b121b221),B2=(b112b122b122b222),x2=y,x=(x1y).
We rewrite system (2.2) as
(2.10)x1(k+1)=ax1(k)+xT(k)B1x(k),y(k+1)=y(k)+xT(k)B2x(k).
To investigate the stability of the zero solution, we use, in accordance with the direct Lyapunov method, an appropriate Lyapunov function V. Let a matrix H, defined as
(2.11)H=(hh12h12h22),
where instead of the entry h11 we put the number h, be positive definite. We set
(2.12)V(x(k))=V(x1(k),y(k)):=xT(k)Hx(k)=(x1(k),y(k))(hh12h12h22)(x1(k)y(k))=hx12(k)+2h12x1(k)y(k)+h22y2(k).
The first difference of the function V along the trajectories of system (2.10) equals
(2.13)ΔV(x(k))=V(x(k+1))-V(x(k))=hx12(k+1)+2h12x1(k+1)y(k+1)+h22y2(k+1)-hx12(k)-2h12x1(k)y(k)-h22y2(k)=h[ax1(k)+xT(k)B1x(k)]2+2h12[ax1(k)+xT(k)B1x(k)][y(k)+xT(k)B2x(k)]+h22[y(k)+xT(k)B2x(k)]2-hx12(k)-2h12x1(k)y(k)-h22y2(k)=h(a2-1)x12(k)+2h12(a-1)x1(k)y(k)+h[2ax1(k)[xT(k)B1x(k)]+[xT(k)B1x(k)]2]+2h12[[xT(k)B1x(k)][xT(k)B2x(k)]ax1(k)[xT(k)B2x(k)]+y(k)[xTB1x(k)]+[xT(k)B1x(k)][xT(k)B2x(k)]]+h22[2y(k)[xT(k)B2x(k)]+[xT(k)B2x(k)]2].
It is easy to see that ΔV does not preserve the sign if h12≠0. Therefore, we put h12=0 and ΔV reduces to
(2.14)ΔV(x(k))=h(a2-1)x12(k)+h[2ax1(k)[xT(k)B1x(k)]+[xT(k)B1x(k)]2]+h22[2y(k)[xT(k)B2x(k)]+[xT(k)B2x(k)]2].
In the polynomial ΔV, with respect to x1 and y, we will put together the third-degree terms (the expression F3(x1(k),y(k)) below) and the fourth-degree terms (the expression F4(x1(k),y(k)) below). In the computations we use the formulas
(2.15)xT(k)Bix(k)=b11ix12(k)+2b12ix1(k)y(k)+b22iy2(k),i=1,2,[xT(k)Bix(k)]2=(b11i)2x14(k)+4(b12i)2x12(k)y2(k)+(b22i)2y4(k)+4b11ib12ix13(k)y(k)+2b11ib22ix12(k)y2(k)+4b12ib22ix1(k)y3(k),i=1,2.
We get
(2.16)ΔV(x(k))=h(a2-1)x12(k)+F3(x1(k),y(k))+F4(x1(k),y(k)),
where
(2.17)F3(x1(k),y(k))=2hab111x13(k)+2[2hab121+h22b112]x12(k)y(k)+2[hab221+2h22b122]x1(k)y2(k)+2h22b222y3(k),F4(x1(k),y(k))=[h(b111)2+h22(b112)2]x14(k)+4[hb111b121+h22b112b122]x13(k)y(k)+2[2h(b121)2+2h22(b122)2+hb111b221+h22b112b222]x12(k)y2(k)+4[hb121b221+h22b122b222]x1(k)y3(k)+[h(b221)2+h22(b222)2]y4(k).
Analysing the increment of V, we see that, if |a|<1, ΔV will be nonpositive in a small neighborhood of the zero solution if the multipliers of the terms x1y2, y3 and x1y3 are equal to zero and the multiplier of the term y4 is nonpositive, that is, if
(2.18)hab221+2h22b122=0,h22b222=0,hb121b221+h22b122b222=0,h(b221)2+h22(b222)2≤0.
As long as the Lyapunov function is positive definite, h>0 and h22>0. Therefore, conditions (2.18) hold if and only if
(2.19)b221=0,b122=0,b222=0.
Then, system (2.2) turns into
(2.20)x1(k+1)=ax1(k)+[b111x12(k)+2b121x1(k)x2(k)],x2(k+1)=x2(k)+b112x12(k)
and ΔV (without loss of generality, we put h22=1, i.e., V(x1,y)=hx12+y2) into
(2.21)ΔV(x(k))=-[h(1-a2)-2hab111x1(k)-2[2hab121+b112]y(k)-[h(b111)2+(b112)2]x12(k)-4hb111b121x1(k)y(k)-4h(b121)2y2(k)]x12(k)=-[α-2β1x1(k)-2β2y(k)-γ1x12(k)-2δ1x1(k)y(k)-γ2y2(k)]x12(k).
The first difference of the Lyapunov function is nonpositive in a sufficiently small neighborhood of the origin (this is because h>0, |a|<1, and α=h(1-a2)>0). In other words, the zero solution is stable in the Lyapunov sense.
Now we will discuss the shape of the guaranteed domain of stability. It can be defined by the inequalities(2.22)γ1x12+2δ1x1y+γ2y2+2β1x1+2β2y≤α,hx12+y2≤r2,
where r>0. This means that inequalities (2.4) and (2.5) are correct. Both inequalities geometrically express closed ellipses if b112·b121≠0. For the second inequality, this is obvious. For the first one, this follows from the following inequalities: γ1>0, γ2>0 and
(2.23)γ1γ2-δ12=[h(b111)2+(b112)2]·[4h(b121)2]-4[hb111b121]2=4h(b112b121)2>0.
Moreover, for r→0, the ellipse (2.4)
(2.24)hx12+y2≤r2
is contained (because it shrinks to the origin) in the ellipse (2.5), that is, there exists such r=r* that, for r∈(0,r*), the ellipse (2.4) lies inside the ellipse (2.5) without any intersection points and, for r=r*, there exists at least one common boundary point of both ellipses. Let us find the value r*. It is characterized by the requirement that the slope coefficients k1 and k2 of both ellipses are the same at the point of contact. Therefore
(2.25)k1=-γ1x1+δ1y+β1δ1x1+γ2y+β2,k2=-hx1y,
where we assume (without loss of generality) that the denominators are nonzero. Thus, we get a quadratic system of two equations to find the contact points (x1,y):
(2.26)γ1x12+2δ1x1y+γ2y2+2β1x1+2β2y=α,hx1(δ1x1+γ2y+β2)-y(γ1x1+δ1y+β1)=0.
For the corresponding values of r, we have
(2.27)hx12+y2=r2.
In accordance with the geometrical meaning of the above quadratic system, we take such a solution (x1,y) as a defintion of the minimal positive value of r and set r*=r.Example 2.2.
Consider a system(2.28)x(k+1)=0.5x(k)+x2(k)-4x(k)y(k),y(k+1)=y(k)+x2(k).
In our case, n=2, a=0.5<1, and b221=b122=b222=0. Therefore, by Theorem 2.1, the zero solution of system (2.28) is stable in the Lyapunov sense. We will find the guaranteed domain of stability. We have
(2.29)b111=1,b121=-2,b112=1,
and b112·b121=-2≠0. Set h=2. Then
(2.30)α=h(1-a2)=2(1-0.25)=1.5,β1=hab111=2·0.5·1=1,β2=2hab121+b112=2·2·0.5·(-2)+1=-3,γ1=h(b111)2+(b112)2=2·12+12=3,γ2=4h(b121)2=4·2·(-2)2=32,δ1=2hb111b121=2·2·1·(-2)=-8.
That is, the guaranteed domain of stability is given by the inequalities
(2.31)3x2-16xy+32y2+2x-6y≤1.5,(2.32)2x2+y2≤r2
if r is so small that the domain described by inequality (2.32) is embedded in the domain described by inequality (2.31). We consider the case when the ellipse (2.32) is embedded in the ellipse (2.31) and the boundaries of both ellipses have only one intersection point. Solving the system (2.8), that is, the system
(2.33)3x2-16xy+32y2+2x-6y=1.5,(2.34)2x(-8x+32y-3)-y(3x-8y+1)=0,
with Mathematica software, we get the solutions (see Figure 1 where the x-axis is identified with the horizontal line and the y-axis is identified with the vertical line, the blue ellipse graphically depicts equation (2.33), and the red hyperbola graphically depicts equation (2.34)):
(2.35)(x,y)=(x1,y1)=(-1.60766,-0.31220),(x,y)=(x2,y2)=(-0.03568,-0.32187),(x,y)=(x3,y3)=(0.01952,-0.13664),(x,y)=(x4,y4)=(1.10728,0.37750).
Then, in accordance with (2.7),
(2.36)r*=mini=1,2,3,42xi2+yi2=2x32+y32≐0.1369,
and the guaranteed domain of stability
(2.37)2x2+y2≤(r*)2=(0.1369)2
obtained from (2.31), (2.32) is depicted in Figure 2 (as an ellipsoidal domain shaded in red and bounded by the thick red ellipse, with the identification of x-axis and y-axis being the same as before). Here, the domain (2.31) is bounded by the blue ellipse (2.33).Figure 1
Graphical solution of system (2.33) and (2.34).Figure 2
The guaranteed domain of stability.


### 2.3. Stability in the Generaln-Dimensional Case

Consider system (1.6) in ℝn. Assume that the matrix A has a simple eigenvalue that is equal to unity with the others lying inside the unit circle. After linearly transforming the dependent variables if necessary, we can assume, without loss of generality, that the matrix A of the linear terms in a block form, that is, (2.38)A=(A0θθT1),A0=(aij),i,j=1,2,…,n-1,
where θ=(0,0,…,0)T, is the (n-1)-dimensional zero vector and all the eigenvalues of the matrix A0 lie inside the unit circle. In order to formulate the next result and its proof, we have to introduce some new definitions (they copy the ones used in Section 1.1, but we use dimension or size n-1 instead of n and note this change as a subscript if necessary): (2.39)x(n-1)=(x1,x2,…,xn-1)T,y=xn,Bi0=(b11ib12i⋯b1,n-1ib21ib22i⋯b2,n-1i⋯⋯⋯⋯bn-1,1ibn-1,2i⋯bn-1,n-1i),i=1,2,…,n,B̃=(b1n1⋯bn-1,n1⋯⋯⋯b1nn-1⋯bn-1,nn-1),B¯T=(b111⋯b1,n-11b211⋯b2,n-11⋯bn-1,11⋯bn-1,n-11⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯b11n-1⋯b1,n-1n-1b21n-1⋯b2,n-1n-1⋯bn-1,1n-1⋯bn-1,n-1n-1).
Matrices Bi0, i=1,2,…,n, are symmetric since bsqi=bqsi,i,s,q=1,2,…,n (see Section 1.1). Moreover, we assume that there exists a symmetric positive definite (n-1)×(n-1) matrix H such that the symmetric matrix (2.40)C=H-A0THA0
is positive definite. Let h>0 be a positive number and (2.41)α=λmin(C),β1=12∥A0THB¯T∥,β2=12∥B̃HA0+3A0THB̃T+2h(Bn0)T∥,γ1=∥B¯HB¯T∥+h∥Bn0∥2,γ2=4∥B̃HB̃T∥,δ1=2∥B¯∥·∥H∥·∥B̃∥.Theorem 2.3.
Leth and r be positive numbers. Assume that
(2.42)bnn1=bnn2=⋯=bnnn=0,b1nn=b2nn=⋯=bn-1,nn=0.
Then the zero solution of system (1.6) is stable by Lyapunov and the guaranteed domain of stability is described by the inequalities
(2.43)γ1∥x∥2+2δ1∥x∥|y|+γ2y2+2β1∥x∥+2β2|y|≤α,(2.44)x(n-1)THx(n-1)+hy2≤r2
if r is so small that the domain described by inequality (2.44) is embedded into the domain described by inequality (2.43).Proof.
We will perform auxiliary matrix computations. With this in mind, we have defined an(n-1)2×(n-1) matrix X(n-1) as
(2.45)X(n-1)T=(X1(n-1)T,X2(n-1)T,…,Xn-1(n-1)T),
where all the elements of the (n-1)×(n-1) matrices Xi(n-1)T, i=1,2,…,n-1 are equal to zero except the row i, which equals x(n-1)T, that is,
(2.46)Xi(n-1)T=(00⋯0⋯⋯⋯⋯x1x2⋯xn-1⋯⋯⋯⋯00⋯0).
Moreover, we define  (a)
vectorsYi, i=1,2,…,n-1, as a row (n-1)-dimensional vector with coordinates equal to zero except the ith element, which equals xn, that is,
(2.47)Yi=(0,0,…,0,xn,0,…,0),(b)
(n-1)×(n-1) zero matrix Θ,(c)
vectorsbi=(b1ni,b2ni,…,bn-1,ni)T,i=1,2,…,n,(d)
vectorb̃=(bnn1,bnn2,…,bnnn-1)T.It is easy to see that
(2.48)XT(k)=(X1(n-1)T(k)Y1T(k)⋯Xn-1(n-1)T(k)Yn-1T(k)ΘθθT0⋯θT0x(n-1)T(k)y(k)),B=(B10b1b1Tbnn1⋯⋯Bn0bnbnTbnnn).
Now we are able to rewrite system (1.6) in an equivalent form
(2.49)(x(n-1)(k+1)y(k+1))=(A0θθT1)(x(n-1)(k)y(k))+(X1(n-1)T(k)Y1T(k)⋯Xn-1(n-1)T(k)Yn-1T(k)ΘθθT0⋯θT0x(n-1)T(k)y(k))×(B10b1b1Tbnn1⋯⋯Bn0bnbnTbnnn)(x(n-1)(k)y(k))=(A0+r11r12r211+r22)(x(n-1)(k)y(k)),
where
(2.50)r11=r11(x(n-1)(k),y(k))=∑j=1n-1[Xj(n-1)T(k)Bj0+YjT(k)bjT]=B¯TX(n-1)(k)+B̃y(k),r12=r12(x(n-1)(k),y(k))=∑j=1n-1[Xj(n-1)T(k)bj+YjT(k)bnnj]=B̃Tx(n-1)(k)+b̃y(k),r21=r21(x(n-1)(k),y(k))=x(n-1)T(k)Bn0+y(k)bnT,r22=r22(x(n-1)(k),y(k))=x(n-1)T(k)bn+y(k)bnnn.
Before the following computations, for the reader's convenience, we recall that for the (n-1)×(n-1) matrices 𝒜, 𝒜1, 1×(n-1) vectors ℓ, ℓ1, (n-1)×1 vectors 𝒞, 𝒞1 and 1×1 “matrices" m, m1, we have
(2.51)(𝒜𝒞ℓm)×(𝒜1𝒞1ℓ1m1)=(𝒜×𝒜1+𝒞×ℓ1A×𝒞1+𝒞×m1ℓ×𝒜1+m×ℓ1ℓ×𝒞1+m×m1).
To investigate the stability of system (1.6), we use the Lyapunov function
(2.52)V(x(n-1)(k),y(k))=x(n-1)T(k)Hx(n-1)(k)+hy2(k)=(x(n-1)T(k),y(k))(HθθTh)(x(n-1)(k)y(k)),
where H=H(n-1) is an (n-1)×(n-1) constant real symmetric and positive definite matrix. Let us find the first difference of the Lyapunov function (2.52) along the solutions of (2.49). We get
(2.53)ΔV(x(n-1)(k),y(k))=(x(n-1)T(k+1),y(k+1))(HθθTh)(x(n-1)(k+1)y(k+1))-(x(n-1)T(k),y(k))(HθθTh)(x(n-1)(k)y(k))=(x(n-1)T(k),y(k))(A0+r11r12r211+r22)T(HθθTh)(A0+r11r12r211+r22)(x(n-1)(k)y(k))-(x(n-1)T(k),y(k))(HθθTh)(x(n-1)(k)y(k))=(x(n-1)T(k),y(k)){(A0T+r11Tr21Tr12T1+r22)(HθθTh)×(A0+r11r12r211+r22)-(HθθTh)}×(x(n-1)(k)y(k))
or, using formula (2.51),
(2.54)ΔV(x(n-1)(k),y(k))=(x(n-1)T(k),y(k)){(A0TH+r11THhr21Tr12THh(1+r22))×(A0+r11r12r211+r22)-(HθθTh)}×(x(n-1)(k)y(k))=(x(n-1)T(k),y(k))×((A0T+r11T)H(A0+r11)+r21Thr21-H(A0T+r11T)Hr12+r21Th(1+r22)r12TH(A0+r11)+(1+r22)hr21r12THr12+(1+r22)h(1+r22)-h)×(x(n-1)(k)y(k)).
Using formulas (2.50), we have
(2.55)ΔV(x(n-1)(k),y(k))=(x(n-1)T(k),y(k))×(c11c12c21c22)(x(n-1)(k)y(k)),
where
(2.56)c11=c11(x(n-1)(k),y(k))=[A0+B¯TX(n-1)(k)+B̃y(k)]TH[A0+B¯TX(n-1)(k)+B̃y(k)]+h[(Bn0)Tx(n-1)(k)+bny(k),x(n-1)T(k)Bn0+y(k)bnT]-H,c12=c12(x(n-1)(k),y(k))=[A0+B¯TX(n-1)(k)+B̃y(k)]TH[B̃Tx(n-1)(k)+b̃y(k)]+h(1+x(n-1)T(k)bn+bnnny(k))[(Bn0)Tx(n-1)(k)+bny(k)],c21=c21(x(n-1)(k),y(k))=c12T,c22=c22(x(n-1)(k),y(k))=[B̃Tx(n-1)(k)+b̃y(k)]TH[B̃Tx(n-1)(k)+b̃y(k)]+h(1+x(n-1)T(k)bn+bnnny(k))2-h.
Then
(2.57)ΔV(x(n-1)(k),y(k))=x(n-1)T(k)c11x(n-1)(k)+2x(n-1)T(k)c12y(k)+y(k)c22y(k)=x(n-1)T[[A0+B¯TX(n-1)(k)+B̃y(k)]TH[A0+B¯TX(n-1)(k)+B̃y(k)]+h[(Bn0)Tx(n-1)(k)+bny(k)][x(n-1)T(k)Bn0+y(k)bnT]-H]x(n-1)(k)+2x(n-1)T(k)[[A0+B¯TX(n-1)(k)+B̃y(k)]TH[B̃Tx(n-1)(k)+b̃y(k)]+h(1+x(n-1)T(k)bn+bnnny(k))[(Bn0)Tx(n-1)(k)+bny(k)]]y(k)+y(k)[[B̃Tx(n-1)(k)+b̃y(k)]TH[B̃Tx(n-1)(k)+b̃y(k)]+h(1+x(n-1)T(k)bn+bnnny(k))2-h]y(k).
After further computation, we get
(2.58)ΔV(x(n-1)(k),y(k))=x(n-1)T(k)[A0THA0-H+A0THB¯TX(n-1)(k)+X(n-1)TB¯HA0+X(n-1)T(k)B¯HB¯TX(n-1)(k)+[X(n-1)T(k)B¯HB̃+B̃THB¯TX(n-1)(k)]y(k)+[B̃THA0+A0THB̃T]y(k)+B̃THB̃y2(k)+h(Bn0)Tx(n-1)(k)x(n-1)T(k)Bn0+2h(Bn0)Tx(n-1)(k)bnTy(k)+hbny2(k)bnT]x(n-1)(k)+2x(n-1)T(k)[A0THB̃Tx(n-1)(k)+A0THb̃y(k)+X(n-1)T(k)B¯HB̃Tx(n-1)(k)+X(n-1)T(k)B¯Hb̃y(k)+B̃HB̃Tx(n-1)(k)y(k)+B̃Hb̃y2(k)+h[(Bn0)Tx(n-1)(k)+bny(k)]+hx(n-1)T(k)bn[(Bn0)Tx(n-1)(k)+bny(k)]+hbnnny(k)[(Bn0)Tx(n-1)(k)+bny(k)]]y(k)+[x(n-1)T(k)B̃H[B̃Tx(n-1)(k)+b̃y(k)]+b̃TH[B̃Tx(n-1)(k)+b̃y(k)]y(k)+h(x(n-1)T(k)bn)2+h(bnnny(k))2+2hx(n-1)T(k)bn+2hbnnny(k)+2hx(n-1)T(k)bnbnnny(k)]y2(k).
Now we can represent ΔV as
(2.59)ΔV(x(n-1)(k),y(k))=F2(x(n-1)(k))+F3(x(n-1)(k),y(k))+F4(x(n-1)(k),y(k)),
where F2 contains only second-order polynomial terms, F3 third-order polynomial terms and F4 fourth-order polynomial terms with respect to the dependent variables. For F2, we have
(2.60)F2(x(n-1)(k))=-x(n-1)T(k)[H-A0THA0]x(n-1)(k)=-x(n-1)T(k)Cx(n-1)(k)≤-λmin(C)∥x(n-1)(k)∥2.
For F3, we get
(2.61)F3(x(n-1)(k),y(k))=F30(x(n-1)(k))+F21(x(n-1)(k),y(k))+F12(x(n-1)(k),y(k))+F03(y(k)),
where
(2.62)F30(x(n-1)(k))=x(n-1)T(k)[A0THB¯TX(n-1)(k)+X(n-1)T(k)B¯HA0]x(n-1)(k),F21(x(n-1)(k),y(k))=x(n-1)T(k)[B̃THA0+3A0THB̃T+2h(Bn0)T]x(n-1)(k)y(k),F12(x(n-1)(k),y(k))=2x(n-1)T(k)[A0THb̃+2hbn]y2(k),F03(y(k))=2hbnnny3(k).
Finally, F4 can be represented as
(2.63)F4(x(n-1)(k),y(k))=F40(x(n-1)(k))+F31(x(n-1)(k),y(k))+F22(x(n-1)(k),y(k))+F13(x(n-1)(k),y(k))+F04(y(k)),
where
(2.64)F40(x(n-1)(k))=x(n-1)T[X(n-1)T(k)B¯HB¯TX(n-1)(k)+h(Bn0)Tx(n-1)(k)x(n-1)T(k)Bn0]x(n-1)(k),F31(x(n-1)(k),y(k))=x(n-1)T[X(n-1)T(k)B¯HB̃+2X(n-1)T(k)B¯HB̃T+B̃THB¯TX(n-1)(k)+2h(Bn0)Tx(n-1)(k)bnT+hx(n-1)T(k)bn(Bn0)T]y(k)x(n-1)(k),F22(x(n-1)(k),y(k))=x(n-1)T[3B̃HB̃Tx(n-1)(k)+B̃THB̃x(n-1)(k)+4h(x(n-1)(k)bn)bn+2X(n-1)T(k)B¯Hb̃+2hbnnn(Bn0)Tx(n-1)(k)]y2(k),F13(x(n-1)(k),y(k))=x(n-1)T[3B̃Hb̃+4hbnnnbn]y3(k),F04(y(k))=[b̃THb̃+h(bnnn)2]y4(k).
Assume the parameters of the system and the coefficients of matrix H to be such that
(2.65)F12(x(n-1)(k),y(k))≡0,F03(y(k))≡0,F13(x(n-1)(k),y(k))≡0,F04(y(k))≡0.
Identities (2.65) hold if and only if
(2.66)A0THb̃+2hbn=0,hbnnn=0,3B̃Hb̃+4hbnnnbn=0,b̃THb̃+h(bnnn)2=0.
As the matrix H is positive definite, (2.66) hold if
(2.67)bnnn=0,b̃=0,bn=0
(these relations coincide with (2.42)). Hence
(2.68)F30(x(n-1)(k))≤∥A0THB¯T∥·∥x(n-1)(k)∥3,F21(x(n-1)(k),y(k))=∥B̃HA0+3A0THB̃T+2h(Bn0)T∥·∥x(n-1)(k)∥2|y(k)|,F40(x(n-1)(k))≤(∥B¯HB¯T∥+h∥Bn0∥2)·∥x(n-1)(k)∥4,F31(x(n-1)(k),y(k))≤4∥B¯∥·∥H∥∥B̃∥·∥x(n-1)(k)∥3·|y(k)|,F22(x(n-1)(k),y(k))≤4∥B̃HB̃T∥·∥x(k)∥2·|y(k)|2.
The first difference (2.59) of V(x(n-1)(k),y(k)) can be estimated as
(2.69)ΔV(x(n-1)(k),y(k))≤-{α-2β1∥x(n-1)(k)∥-2β2|y(k)|-γ1∥x(n-1)(k)∥2-2δ1∥x(n-1)(k)∥·|y(k)|-γ2|y(k)|2}∥x(n-1)(k)∥2
and is nonpositive if
(2.70)γ1∥x(n-1)(k)∥2+2δ1∥x(n-1)(k)∥·|y(k)|+γ2y2(k)+2β1∥x(n-1)(k)∥+2β2|y(k)|≤α.
The stability domain is defined by the inequality
(2.71)x(n-1)T(k)Hx(n-1)(k)+h|y(k)|≤r2
supposing that r is so small that the domain (2.71) is embedded in the domain (2.70).


## 2.1. Instability in One-Dimensional Case

We start by discussing a simple scalar equation with the eigenvalue of matrixA equaling one, that is, a11=1. Then (1.6) takes the form (2.1)x(k+1)=x(k)+bx2(k),k=0,…,
and it is easy to see that the trivial solution is unstable for an arbitrary b≠0 (to show this, we can apply, e.g., Theorem 1.15 [4, page 29]).This elementary example shows that stability in the case of system (1.6) has an extraordinary significance and the results on stability (for n≠1) lose their meaning for n=1 when we deal with instability. We show that, if n≠1 and B satisfies certain assumptions, the zero solution is stable. Moreover, the shape of the guaranteed domain of stability will be given.We divide our forthcoming analysis into two parts. In the first one we give an explicit coefficient criterion in the subcase ofn=2. Then we consider the general n-dimensional case.


## 2.2. Stability in the General Two-Dimensional Case

Letn=2. Then system (1.6) with the matrix A having a simple eigenvalue λ=1 reduces (after linearly transforming the dependent variables if necessary) to(2.2)x1(k+1)=ax1(k)+[b111x12(k)+2b121x1(k)x2(k)+b221x22(k)],x2(k+1)=x2(k)+[b112x12(k)+2b122x1(k)x2(k)+b222x22(k)].
We will assume that |a|<1. Define auxiliary numbers as follows: (2.3)α=h(1-a2),β1=hab111,β2=2hab121+b112,γ1=h(b111)2+(b112)2,γ2=4h(b121)2,δ1=2hb111b121,
where h is a positive number.Theorem 2.1.
Leth and r be positive numbers. Assume that |a|<1 and b122=b221=b222=0. Then the zero solution of system (2.2) is stable in the Lyapunov sense and a guaranteed domain of stability is given by the inequality
(2.4)hx12+x22≤r2
if r is taken so small that the domain described by (2.4) is embedded in the domain
(2.5)γ1x12+2δ1x1x2+γ2x22+2β1x1+2β2x2≤α.
If, moreover, b112·b121≠0, then a guaranteed domain of stability can be described using inequality
(2.6)hx12+x22≤(r*)2
with
(2.7)r*=min(x1,x2)hx12+x22,
where (x1,x2) runs over all real solutions of the nonlinear system with unknowns x1 and x2:
(2.8)γ1x12+2δ1x1x2+γ2x22+2β1x1+2β2x2=α,hx1(δ1x1+γ2x2+β2)-x2(γ1x1+δ1x2+β1)=0.Proof.
Define(2.9)B1=(b111b121b121b221),B2=(b112b122b122b222),x2=y,x=(x1y).
We rewrite system (2.2) as
(2.10)x1(k+1)=ax1(k)+xT(k)B1x(k),y(k+1)=y(k)+xT(k)B2x(k).
To investigate the stability of the zero solution, we use, in accordance with the direct Lyapunov method, an appropriate Lyapunov function V. Let a matrix H, defined as
(2.11)H=(hh12h12h22),
where instead of the entry h11 we put the number h, be positive definite. We set
(2.12)V(x(k))=V(x1(k),y(k)):=xT(k)Hx(k)=(x1(k),y(k))(hh12h12h22)(x1(k)y(k))=hx12(k)+2h12x1(k)y(k)+h22y2(k).
The first difference of the function V along the trajectories of system (2.10) equals
(2.13)ΔV(x(k))=V(x(k+1))-V(x(k))=hx12(k+1)+2h12x1(k+1)y(k+1)+h22y2(k+1)-hx12(k)-2h12x1(k)y(k)-h22y2(k)=h[ax1(k)+xT(k)B1x(k)]2+2h12[ax1(k)+xT(k)B1x(k)][y(k)+xT(k)B2x(k)]+h22[y(k)+xT(k)B2x(k)]2-hx12(k)-2h12x1(k)y(k)-h22y2(k)=h(a2-1)x12(k)+2h12(a-1)x1(k)y(k)+h[2ax1(k)[xT(k)B1x(k)]+[xT(k)B1x(k)]2]+2h12[[xT(k)B1x(k)][xT(k)B2x(k)]ax1(k)[xT(k)B2x(k)]+y(k)[xTB1x(k)]+[xT(k)B1x(k)][xT(k)B2x(k)]]+h22[2y(k)[xT(k)B2x(k)]+[xT(k)B2x(k)]2].
It is easy to see that ΔV does not preserve the sign if h12≠0. Therefore, we put h12=0 and ΔV reduces to
(2.14)ΔV(x(k))=h(a2-1)x12(k)+h[2ax1(k)[xT(k)B1x(k)]+[xT(k)B1x(k)]2]+h22[2y(k)[xT(k)B2x(k)]+[xT(k)B2x(k)]2].
In the polynomial ΔV, with respect to x1 and y, we will put together the third-degree terms (the expression F3(x1(k),y(k)) below) and the fourth-degree terms (the expression F4(x1(k),y(k)) below). In the computations we use the formulas
(2.15)xT(k)Bix(k)=b11ix12(k)+2b12ix1(k)y(k)+b22iy2(k),i=1,2,[xT(k)Bix(k)]2=(b11i)2x14(k)+4(b12i)2x12(k)y2(k)+(b22i)2y4(k)+4b11ib12ix13(k)y(k)+2b11ib22ix12(k)y2(k)+4b12ib22ix1(k)y3(k),i=1,2.
We get
(2.16)ΔV(x(k))=h(a2-1)x12(k)+F3(x1(k),y(k))+F4(x1(k),y(k)),
where
(2.17)F3(x1(k),y(k))=2hab111x13(k)+2[2hab121+h22b112]x12(k)y(k)+2[hab221+2h22b122]x1(k)y2(k)+2h22b222y3(k),F4(x1(k),y(k))=[h(b111)2+h22(b112)2]x14(k)+4[hb111b121+h22b112b122]x13(k)y(k)+2[2h(b121)2+2h22(b122)2+hb111b221+h22b112b222]x12(k)y2(k)+4[hb121b221+h22b122b222]x1(k)y3(k)+[h(b221)2+h22(b222)2]y4(k).
Analysing the increment of V, we see that, if |a|<1, ΔV will be nonpositive in a small neighborhood of the zero solution if the multipliers of the terms x1y2, y3 and x1y3 are equal to zero and the multiplier of the term y4 is nonpositive, that is, if
(2.18)hab221+2h22b122=0,h22b222=0,hb121b221+h22b122b222=0,h(b221)2+h22(b222)2≤0.
As long as the Lyapunov function is positive definite, h>0 and h22>0. Therefore, conditions (2.18) hold if and only if
(2.19)b221=0,b122=0,b222=0.
Then, system (2.2) turns into
(2.20)x1(k+1)=ax1(k)+[b111x12(k)+2b121x1(k)x2(k)],x2(k+1)=x2(k)+b112x12(k)
and ΔV (without loss of generality, we put h22=1, i.e., V(x1,y)=hx12+y2) into
(2.21)ΔV(x(k))=-[h(1-a2)-2hab111x1(k)-2[2hab121+b112]y(k)-[h(b111)2+(b112)2]x12(k)-4hb111b121x1(k)y(k)-4h(b121)2y2(k)]x12(k)=-[α-2β1x1(k)-2β2y(k)-γ1x12(k)-2δ1x1(k)y(k)-γ2y2(k)]x12(k).
The first difference of the Lyapunov function is nonpositive in a sufficiently small neighborhood of the origin (this is because h>0, |a|<1, and α=h(1-a2)>0). In other words, the zero solution is stable in the Lyapunov sense.
Now we will discuss the shape of the guaranteed domain of stability. It can be defined by the inequalities(2.22)γ1x12+2δ1x1y+γ2y2+2β1x1+2β2y≤α,hx12+y2≤r2,
where r>0. This means that inequalities (2.4) and (2.5) are correct. Both inequalities geometrically express closed ellipses if b112·b121≠0. For the second inequality, this is obvious. For the first one, this follows from the following inequalities: γ1>0, γ2>0 and
(2.23)γ1γ2-δ12=[h(b111)2+(b112)2]·[4h(b121)2]-4[hb111b121]2=4h(b112b121)2>0.
Moreover, for r→0, the ellipse (2.4)
(2.24)hx12+y2≤r2
is contained (because it shrinks to the origin) in the ellipse (2.5), that is, there exists such r=r* that, for r∈(0,r*), the ellipse (2.4) lies inside the ellipse (2.5) without any intersection points and, for r=r*, there exists at least one common boundary point of both ellipses. Let us find the value r*. It is characterized by the requirement that the slope coefficients k1 and k2 of both ellipses are the same at the point of contact. Therefore
(2.25)k1=-γ1x1+δ1y+β1δ1x1+γ2y+β2,k2=-hx1y,
where we assume (without loss of generality) that the denominators are nonzero. Thus, we get a quadratic system of two equations to find the contact points (x1,y):
(2.26)γ1x12+2δ1x1y+γ2y2+2β1x1+2β2y=α,hx1(δ1x1+γ2y+β2)-y(γ1x1+δ1y+β1)=0.
For the corresponding values of r, we have
(2.27)hx12+y2=r2.
In accordance with the geometrical meaning of the above quadratic system, we take such a solution (x1,y) as a defintion of the minimal positive value of r and set r*=r.Example 2.2.
Consider a system(2.28)x(k+1)=0.5x(k)+x2(k)-4x(k)y(k),y(k+1)=y(k)+x2(k).
In our case, n=2, a=0.5<1, and b221=b122=b222=0. Therefore, by Theorem 2.1, the zero solution of system (2.28) is stable in the Lyapunov sense. We will find the guaranteed domain of stability. We have
(2.29)b111=1,b121=-2,b112=1,
and b112·b121=-2≠0. Set h=2. Then
(2.30)α=h(1-a2)=2(1-0.25)=1.5,β1=hab111=2·0.5·1=1,β2=2hab121+b112=2·2·0.5·(-2)+1=-3,γ1=h(b111)2+(b112)2=2·12+12=3,γ2=4h(b121)2=4·2·(-2)2=32,δ1=2hb111b121=2·2·1·(-2)=-8.
That is, the guaranteed domain of stability is given by the inequalities
(2.31)3x2-16xy+32y2+2x-6y≤1.5,(2.32)2x2+y2≤r2
if r is so small that the domain described by inequality (2.32) is embedded in the domain described by inequality (2.31). We consider the case when the ellipse (2.32) is embedded in the ellipse (2.31) and the boundaries of both ellipses have only one intersection point. Solving the system (2.8), that is, the system
(2.33)3x2-16xy+32y2+2x-6y=1.5,(2.34)2x(-8x+32y-3)-y(3x-8y+1)=0,
with Mathematica software, we get the solutions (see Figure 1 where the x-axis is identified with the horizontal line and the y-axis is identified with the vertical line, the blue ellipse graphically depicts equation (2.33), and the red hyperbola graphically depicts equation (2.34)):
(2.35)(x,y)=(x1,y1)=(-1.60766,-0.31220),(x,y)=(x2,y2)=(-0.03568,-0.32187),(x,y)=(x3,y3)=(0.01952,-0.13664),(x,y)=(x4,y4)=(1.10728,0.37750).
Then, in accordance with (2.7),
(2.36)r*=mini=1,2,3,42xi2+yi2=2x32+y32≐0.1369,
and the guaranteed domain of stability
(2.37)2x2+y2≤(r*)2=(0.1369)2
obtained from (2.31), (2.32) is depicted in Figure 2 (as an ellipsoidal domain shaded in red and bounded by the thick red ellipse, with the identification of x-axis and y-axis being the same as before). Here, the domain (2.31) is bounded by the blue ellipse (2.33).Figure 1
Graphical solution of system (2.33) and (2.34).Figure 2
The guaranteed domain of stability.


## 2.3. Stability in the Generaln-Dimensional Case

Consider system (1.6) in ℝn. Assume that the matrix A has a simple eigenvalue that is equal to unity with the others lying inside the unit circle. After linearly transforming the dependent variables if necessary, we can assume, without loss of generality, that the matrix A of the linear terms in a block form, that is, (2.38)A=(A0θθT1),A0=(aij),i,j=1,2,…,n-1,
where θ=(0,0,…,0)T, is the (n-1)-dimensional zero vector and all the eigenvalues of the matrix A0 lie inside the unit circle. In order to formulate the next result and its proof, we have to introduce some new definitions (they copy the ones used in Section 1.1, but we use dimension or size n-1 instead of n and note this change as a subscript if necessary): (2.39)x(n-1)=(x1,x2,…,xn-1)T,y=xn,Bi0=(b11ib12i⋯b1,n-1ib21ib22i⋯b2,n-1i⋯⋯⋯⋯bn-1,1ibn-1,2i⋯bn-1,n-1i),i=1,2,…,n,B̃=(b1n1⋯bn-1,n1⋯⋯⋯b1nn-1⋯bn-1,nn-1),B¯T=(b111⋯b1,n-11b211⋯b2,n-11⋯bn-1,11⋯bn-1,n-11⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯b11n-1⋯b1,n-1n-1b21n-1⋯b2,n-1n-1⋯bn-1,1n-1⋯bn-1,n-1n-1).
Matrices Bi0, i=1,2,…,n, are symmetric since bsqi=bqsi,i,s,q=1,2,…,n (see Section 1.1). Moreover, we assume that there exists a symmetric positive definite (n-1)×(n-1) matrix H such that the symmetric matrix (2.40)C=H-A0THA0
is positive definite. Let h>0 be a positive number and (2.41)α=λmin(C),β1=12∥A0THB¯T∥,β2=12∥B̃HA0+3A0THB̃T+2h(Bn0)T∥,γ1=∥B¯HB¯T∥+h∥Bn0∥2,γ2=4∥B̃HB̃T∥,δ1=2∥B¯∥·∥H∥·∥B̃∥.Theorem 2.3.
Leth and r be positive numbers. Assume that
(2.42)bnn1=bnn2=⋯=bnnn=0,b1nn=b2nn=⋯=bn-1,nn=0.
Then the zero solution of system (1.6) is stable by Lyapunov and the guaranteed domain of stability is described by the inequalities
(2.43)γ1∥x∥2+2δ1∥x∥|y|+γ2y2+2β1∥x∥+2β2|y|≤α,(2.44)x(n-1)THx(n-1)+hy2≤r2
if r is so small that the domain described by inequality (2.44) is embedded into the domain described by inequality (2.43).Proof.
We will perform auxiliary matrix computations. With this in mind, we have defined an(n-1)2×(n-1) matrix X(n-1) as
(2.45)X(n-1)T=(X1(n-1)T,X2(n-1)T,…,Xn-1(n-1)T),
where all the elements of the (n-1)×(n-1) matrices Xi(n-1)T, i=1,2,…,n-1 are equal to zero except the row i, which equals x(n-1)T, that is,
(2.46)Xi(n-1)T=(00⋯0⋯⋯⋯⋯x1x2⋯xn-1⋯⋯⋯⋯00⋯0).
Moreover, we define  (a)
vectorsYi, i=1,2,…,n-1, as a row (n-1)-dimensional vector with coordinates equal to zero except the ith element, which equals xn, that is,
(2.47)Yi=(0,0,…,0,xn,0,…,0),(b)
(n-1)×(n-1) zero matrix Θ,(c)
vectorsbi=(b1ni,b2ni,…,bn-1,ni)T,i=1,2,…,n,(d)
vectorb̃=(bnn1,bnn2,…,bnnn-1)T.It is easy to see that
(2.48)XT(k)=(X1(n-1)T(k)Y1T(k)⋯Xn-1(n-1)T(k)Yn-1T(k)ΘθθT0⋯θT0x(n-1)T(k)y(k)),B=(B10b1b1Tbnn1⋯⋯Bn0bnbnTbnnn).
Now we are able to rewrite system (1.6) in an equivalent form
(2.49)(x(n-1)(k+1)y(k+1))=(A0θθT1)(x(n-1)(k)y(k))+(X1(n-1)T(k)Y1T(k)⋯Xn-1(n-1)T(k)Yn-1T(k)ΘθθT0⋯θT0x(n-1)T(k)y(k))×(B10b1b1Tbnn1⋯⋯Bn0bnbnTbnnn)(x(n-1)(k)y(k))=(A0+r11r12r211+r22)(x(n-1)(k)y(k)),
where
(2.50)r11=r11(x(n-1)(k),y(k))=∑j=1n-1[Xj(n-1)T(k)Bj0+YjT(k)bjT]=B¯TX(n-1)(k)+B̃y(k),r12=r12(x(n-1)(k),y(k))=∑j=1n-1[Xj(n-1)T(k)bj+YjT(k)bnnj]=B̃Tx(n-1)(k)+b̃y(k),r21=r21(x(n-1)(k),y(k))=x(n-1)T(k)Bn0+y(k)bnT,r22=r22(x(n-1)(k),y(k))=x(n-1)T(k)bn+y(k)bnnn.
Before the following computations, for the reader's convenience, we recall that for the (n-1)×(n-1) matrices 𝒜, 𝒜1, 1×(n-1) vectors ℓ, ℓ1, (n-1)×1 vectors 𝒞, 𝒞1 and 1×1 “matrices" m, m1, we have
(2.51)(𝒜𝒞ℓm)×(𝒜1𝒞1ℓ1m1)=(𝒜×𝒜1+𝒞×ℓ1A×𝒞1+𝒞×m1ℓ×𝒜1+m×ℓ1ℓ×𝒞1+m×m1).
To investigate the stability of system (1.6), we use the Lyapunov function
(2.52)V(x(n-1)(k),y(k))=x(n-1)T(k)Hx(n-1)(k)+hy2(k)=(x(n-1)T(k),y(k))(HθθTh)(x(n-1)(k)y(k)),
where H=H(n-1) is an (n-1)×(n-1) constant real symmetric and positive definite matrix. Let us find the first difference of the Lyapunov function (2.52) along the solutions of (2.49). We get
(2.53)ΔV(x(n-1)(k),y(k))=(x(n-1)T(k+1),y(k+1))(HθθTh)(x(n-1)(k+1)y(k+1))-(x(n-1)T(k),y(k))(HθθTh)(x(n-1)(k)y(k))=(x(n-1)T(k),y(k))(A0+r11r12r211+r22)T(HθθTh)(A0+r11r12r211+r22)(x(n-1)(k)y(k))-(x(n-1)T(k),y(k))(HθθTh)(x(n-1)(k)y(k))=(x(n-1)T(k),y(k)){(A0T+r11Tr21Tr12T1+r22)(HθθTh)×(A0+r11r12r211+r22)-(HθθTh)}×(x(n-1)(k)y(k))
or, using formula (2.51),
(2.54)ΔV(x(n-1)(k),y(k))=(x(n-1)T(k),y(k)){(A0TH+r11THhr21Tr12THh(1+r22))×(A0+r11r12r211+r22)-(HθθTh)}×(x(n-1)(k)y(k))=(x(n-1)T(k),y(k))×((A0T+r11T)H(A0+r11)+r21Thr21-H(A0T+r11T)Hr12+r21Th(1+r22)r12TH(A0+r11)+(1+r22)hr21r12THr12+(1+r22)h(1+r22)-h)×(x(n-1)(k)y(k)).
Using formulas (2.50), we have
(2.55)ΔV(x(n-1)(k),y(k))=(x(n-1)T(k),y(k))×(c11c12c21c22)(x(n-1)(k)y(k)),
where
(2.56)c11=c11(x(n-1)(k),y(k))=[A0+B¯TX(n-1)(k)+B̃y(k)]TH[A0+B¯TX(n-1)(k)+B̃y(k)]+h[(Bn0)Tx(n-1)(k)+bny(k),x(n-1)T(k)Bn0+y(k)bnT]-H,c12=c12(x(n-1)(k),y(k))=[A0+B¯TX(n-1)(k)+B̃y(k)]TH[B̃Tx(n-1)(k)+b̃y(k)]+h(1+x(n-1)T(k)bn+bnnny(k))[(Bn0)Tx(n-1)(k)+bny(k)],c21=c21(x(n-1)(k),y(k))=c12T,c22=c22(x(n-1)(k),y(k))=[B̃Tx(n-1)(k)+b̃y(k)]TH[B̃Tx(n-1)(k)+b̃y(k)]+h(1+x(n-1)T(k)bn+bnnny(k))2-h.
Then
(2.57)ΔV(x(n-1)(k),y(k))=x(n-1)T(k)c11x(n-1)(k)+2x(n-1)T(k)c12y(k)+y(k)c22y(k)=x(n-1)T[[A0+B¯TX(n-1)(k)+B̃y(k)]TH[A0+B¯TX(n-1)(k)+B̃y(k)]+h[(Bn0)Tx(n-1)(k)+bny(k)][x(n-1)T(k)Bn0+y(k)bnT]-H]x(n-1)(k)+2x(n-1)T(k)[[A0+B¯TX(n-1)(k)+B̃y(k)]TH[B̃Tx(n-1)(k)+b̃y(k)]+h(1+x(n-1)T(k)bn+bnnny(k))[(Bn0)Tx(n-1)(k)+bny(k)]]y(k)+y(k)[[B̃Tx(n-1)(k)+b̃y(k)]TH[B̃Tx(n-1)(k)+b̃y(k)]+h(1+x(n-1)T(k)bn+bnnny(k))2-h]y(k).
After further computation, we get
(2.58)ΔV(x(n-1)(k),y(k))=x(n-1)T(k)[A0THA0-H+A0THB¯TX(n-1)(k)+X(n-1)TB¯HA0+X(n-1)T(k)B¯HB¯TX(n-1)(k)+[X(n-1)T(k)B¯HB̃+B̃THB¯TX(n-1)(k)]y(k)+[B̃THA0+A0THB̃T]y(k)+B̃THB̃y2(k)+h(Bn0)Tx(n-1)(k)x(n-1)T(k)Bn0+2h(Bn0)Tx(n-1)(k)bnTy(k)+hbny2(k)bnT]x(n-1)(k)+2x(n-1)T(k)[A0THB̃Tx(n-1)(k)+A0THb̃y(k)+X(n-1)T(k)B¯HB̃Tx(n-1)(k)+X(n-1)T(k)B¯Hb̃y(k)+B̃HB̃Tx(n-1)(k)y(k)+B̃Hb̃y2(k)+h[(Bn0)Tx(n-1)(k)+bny(k)]+hx(n-1)T(k)bn[(Bn0)Tx(n-1)(k)+bny(k)]+hbnnny(k)[(Bn0)Tx(n-1)(k)+bny(k)]]y(k)+[x(n-1)T(k)B̃H[B̃Tx(n-1)(k)+b̃y(k)]+b̃TH[B̃Tx(n-1)(k)+b̃y(k)]y(k)+h(x(n-1)T(k)bn)2+h(bnnny(k))2+2hx(n-1)T(k)bn+2hbnnny(k)+2hx(n-1)T(k)bnbnnny(k)]y2(k).
Now we can represent ΔV as
(2.59)ΔV(x(n-1)(k),y(k))=F2(x(n-1)(k))+F3(x(n-1)(k),y(k))+F4(x(n-1)(k),y(k)),
where F2 contains only second-order polynomial terms, F3 third-order polynomial terms and F4 fourth-order polynomial terms with respect to the dependent variables. For F2, we have
(2.60)F2(x(n-1)(k))=-x(n-1)T(k)[H-A0THA0]x(n-1)(k)=-x(n-1)T(k)Cx(n-1)(k)≤-λmin(C)∥x(n-1)(k)∥2.
For F3, we get
(2.61)F3(x(n-1)(k),y(k))=F30(x(n-1)(k))+F21(x(n-1)(k),y(k))+F12(x(n-1)(k),y(k))+F03(y(k)),
where
(2.62)F30(x(n-1)(k))=x(n-1)T(k)[A0THB¯TX(n-1)(k)+X(n-1)T(k)B¯HA0]x(n-1)(k),F21(x(n-1)(k),y(k))=x(n-1)T(k)[B̃THA0+3A0THB̃T+2h(Bn0)T]x(n-1)(k)y(k),F12(x(n-1)(k),y(k))=2x(n-1)T(k)[A0THb̃+2hbn]y2(k),F03(y(k))=2hbnnny3(k).
Finally, F4 can be represented as
(2.63)F4(x(n-1)(k),y(k))=F40(x(n-1)(k))+F31(x(n-1)(k),y(k))+F22(x(n-1)(k),y(k))+F13(x(n-1)(k),y(k))+F04(y(k)),
where
(2.64)F40(x(n-1)(k))=x(n-1)T[X(n-1)T(k)B¯HB¯TX(n-1)(k)+h(Bn0)Tx(n-1)(k)x(n-1)T(k)Bn0]x(n-1)(k),F31(x(n-1)(k),y(k))=x(n-1)T[X(n-1)T(k)B¯HB̃+2X(n-1)T(k)B¯HB̃T+B̃THB¯TX(n-1)(k)+2h(Bn0)Tx(n-1)(k)bnT+hx(n-1)T(k)bn(Bn0)T]y(k)x(n-1)(k),F22(x(n-1)(k),y(k))=x(n-1)T[3B̃HB̃Tx(n-1)(k)+B̃THB̃x(n-1)(k)+4h(x(n-1)(k)bn)bn+2X(n-1)T(k)B¯Hb̃+2hbnnn(Bn0)Tx(n-1)(k)]y2(k),F13(x(n-1)(k),y(k))=x(n-1)T[3B̃Hb̃+4hbnnnbn]y3(k),F04(y(k))=[b̃THb̃+h(bnnn)2]y4(k).
Assume the parameters of the system and the coefficients of matrix H to be such that
(2.65)F12(x(n-1)(k),y(k))≡0,F03(y(k))≡0,F13(x(n-1)(k),y(k))≡0,F04(y(k))≡0.
Identities (2.65) hold if and only if
(2.66)A0THb̃+2hbn=0,hbnnn=0,3B̃Hb̃+4hbnnnbn=0,b̃THb̃+h(bnnn)2=0.
As the matrix H is positive definite, (2.66) hold if
(2.67)bnnn=0,b̃=0,bn=0
(these relations coincide with (2.42)). Hence
(2.68)F30(x(n-1)(k))≤∥A0THB¯T∥·∥x(n-1)(k)∥3,F21(x(n-1)(k),y(k))=∥B̃HA0+3A0THB̃T+2h(Bn0)T∥·∥x(n-1)(k)∥2|y(k)|,F40(x(n-1)(k))≤(∥B¯HB¯T∥+h∥Bn0∥2)·∥x(n-1)(k)∥4,F31(x(n-1)(k),y(k))≤4∥B¯∥·∥H∥∥B̃∥·∥x(n-1)(k)∥3·|y(k)|,F22(x(n-1)(k),y(k))≤4∥B̃HB̃T∥·∥x(k)∥2·|y(k)|2.
The first difference (2.59) of V(x(n-1)(k),y(k)) can be estimated as
(2.69)ΔV(x(n-1)(k),y(k))≤-{α-2β1∥x(n-1)(k)∥-2β2|y(k)|-γ1∥x(n-1)(k)∥2-2δ1∥x(n-1)(k)∥·|y(k)|-γ2|y(k)|2}∥x(n-1)(k)∥2
and is nonpositive if
(2.70)γ1∥x(n-1)(k)∥2+2δ1∥x(n-1)(k)∥·|y(k)|+γ2y2(k)+2β1∥x(n-1)(k)∥+2β2|y(k)|≤α.
The stability domain is defined by the inequality
(2.71)x(n-1)T(k)Hx(n-1)(k)+h|y(k)|≤r2
supposing that r is so small that the domain (2.71) is embedded in the domain (2.70).


## 3. Concluding Remarks

Since 1892, when the general problem of stability of differential equations by the linear approximation was considered by A.M. Lyapunov, investigation of stability by linear approximation has been attracting permanent interest. For example, Malkin [13, Chapter 3] considered a general case of stability by linear approximation and derived stability criteria. In our paper we deal with systems of difference equations of a special form (with quadratic right-hand sides). Comparing our results with possible extensions of Malkin's results to difference equations, we point out that the primary purpose of our results, unlike those of Malkin, is to select such terms of the quadratic right-hand sides as contributing as much as possible to the loss of stability. The domains of stability are described in terms of the coefficients.Since the matrixA has one simple eigenvalue λ=1, a question on asymptotic stability is not asked. In fact, we are looking for a domain V(x)<c, c>0 (embedded in the domain described by inequality ΔV(x)≤0), of admissible initial perturbations, that is, a set of invariance. All solutions generated by the initial data within this set remain in an ɛ neighborhood of the zero steady state provided that the initial perturbations are in a δ(ɛ) neighborhood of the zero steady state in the phase plane.The method presented can be extended to further classes of systems. It is, for example, possible to consider the case of the multiplicity of eigenvalueλ=1 being more than one or the case of |λ|=1 and get stability criteria. But the computations needed are too cumbersome. Other questions are related to the results of Theorem 2.1 and Theorem 2.3. Is the zero solution of system (2.2) or of system (1.6) unstable when the coefficients indicated are not necessarily zero coefficients? We formulate them as open problems—prove or disprove the following conjecture.Conjecture 3.1.
(a)  Let|a|<1 and (b122)2+(b221)2+(b222)2>0. Then the zero solution of system (2.2) is unstable.
(b) Let all the eigenvalues of the matrixA0 lie inside the unit circle and(3.1)∑j=1n(bnnj)2+∑j=1n-1(bjnn)2>0.
Then the zero solution of system (1.6) is unstable.Finally, we formulate an open problem related to the shape of the guaranteed domain of stability.Problem 1.
The guaranteed domain of stability of the zero solution of system (2.2) is described by inequality (2.6) with r* defined by (2.7) where (x1,x2) runs over all real solutions of the nonlinear system (2.8). Is it possible to derive an analogous shape of the guaranteed domain of stability for the zero solution of system (1.6) using inequalities (2.43), (2.44) and an analogy of method applied in the proof of Theorem 2.1?



---

*Source: 539087-2010-06-29.xml*