# Global Dynamical Properties of Rational Higher-Order System of Difference Equations

Authors: A. Q. Khan; S. M. Qureshi
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

In this paper, global dynamical properties of rational higher-order system are explored in the interior ofℝ+3. It is explored that under certain parametric conditions, the discrete-time system has at most eight equilibria. By the method of linearization, local dynamics has been explored. It is explored that positive solution of the system is bounded, and moreover fixed point P000 is globally stable if α1/α2<1,α4/α5<1, α7/α8<1. It is also investigated that the positive solution of the system under consideration converges to P000. Lastly, theoretical results are confirmed by numerical simulation. The presented work is significantly extended and improves current results in the literature.


---

## Body

## 1. Introduction

It is a well-known fact that difference equations arise naturally as discrete analogues and as numerical solutions of differential as well as delay differential equations having applications in many fields like physics, biology, economy, and ecology. Recently, a lot of studies have been conducted concerning the global dynamics of difference equations and their systems [3–19]. It is really not easy to understand global dynamics of difference equations along their systems; particularly, investigating the global behavior of higher-order equations is a challenging job in recent years. Therefore, investigating the global dynamics of such difference equations along their systems is worth further consideration. For illustration, Gibbons et al. [20] explored global dynamics of the following equation:(1)xn+1=α1+α2xn−1α3+xn,where αℓℓ=1,2,3 and x−ℓℓ=1,0 are positive constants. Çinar [21] investigated the dynamics of the following equations:(2)xn+1=xn−11+∏i=01xn−i,where x−ℓℓ=1,0 are positive constants. Shojaei et al. [22] investigated global dynamics of the following difference equations:(3)xn+1=α1xn−2α2+α3∏i=02xn−i,where αℓℓ=1,2,3 and x−ℓℓ=2,1,0 are positive constants. Bajo and Liz [23] investigated the dynamics of the following difference equation:(4)xn+1=xn−1α1+α2∏i=01xn−i,where αℓℓ=1,2 and x−ℓℓ=1,0 are positive constants. Zhang et al. [24] have extended the work explored by numerous authors [21–23] to investigate the dynamics of the following rational system:(5)xn+1=xn−2α2+∏i=02yn−i,yn+1=yn−1α1+∏i=02xn−i,where αℓℓ=1,2 and xℓ,y−ℓℓ=2,1,0 are positive constants. Recently, Qureshi and Khan [25] investigated the global dynamics of the following rational system, which is extension of the work [21–24]:(6)xn+1=α1xn−1α2+α3∏i=01yn−i,yn+1=α4yn−1α5+α6∏i=01zn−i,zn+1=α7zn−1α8+α9∏i=01xn−i,where αℓℓ=1,2,…,9 and x−ℓ,y−ℓ,z−ℓℓ=1,0 are positive constants. Inspired from aforesaid studies, we will extend the work studied by numerous authors [21–25] to investigate global dynamics of the following 3k+3-dimensional system:(7)xn+1=α1xn−kα2+α3∏i=0kyn−i,yn+1=α4yn−kα5+α6∏i=0kzn−i,zn+1=α7zn−kα8+α9∏i=0kxn−i,where αℓℓ=1,2,…,9 and x−ℓ,ℓ−ℓ,z−ℓℓ=k,−k+1,…,1,0 are positive constants.The organization of this paper is as follows. In Section2, existence of equilibria in ℝ+3 and corresponding linearized form are investigated. Section 3 deals with the study of local dynamics about equilibrium points. Boundedness of positive solution for the discrete-time system is studied in Section 4. Further, global dynamics about P000 is explored in Section 5. In Section 6, we studied the rate of convergence which converges to P000 of the system. Theoretical results are numerically verified in Section 7, while concluding remarks are given in Section 8.


## 2. Equilibria and Linearized Form of System (7)

The existence of equilibrium solution in the interior ofℝ+3 and linearized form about Λ=x,y,z of system (7) are investigated in this section. So, existence of equilibrium solution can be summarized as the following lemma.Lemma 1.
In the interior ofℝ+3, discrete-time system (7) has at most eight equilibria. More precisely,(i)
∀αℓℓ=1,2,…,9, P0000,0,0 is the unique boundary point of discrete-time system (7).(ii)
Px0zk+1α7−α8/α9,0,k+1α4−α5/α6 is the boundary equilibrium point of system (7) if α1=α2,α4>α5andα7>α8.(iii)
Pxy0k+1α7−α8/α9,k+1α1−α2/α3,0 is the boundary equilibrium point of system (7) if α1>α2,α4=α5andα7>α8.(iv)
P0yz0,k+1α1−α2/α3,k+1α4−α5/α6 is the boundary equilibrium point of system (7) if α1>α2,α4>α5andα7=α8.(v)
Px00k+1α7−α8/α9,0,0 is the boundary equilibrium point of discrete-time system (7) if α1=α2,α4=α5andα7>α8.(vi)
P00z0,0,k+1α4−α5/α6 is the boundary equilibrium point of system (7) if α1=α2,α4>α5andα7=α8.(vii)
P0y00,k+1α1−α2/α3,0 is the boundary equilibrium point of discrete-time system (7) if α1>α2,α4=α5andα7=α8.(viii)
Pxyzk+1α7−α8/α9,k+1α1−α2/α3,k+1α4−α5/α6 is the unique positive equilibrium point of system (7) if α1>α2,α4>α5andα7>α8.Hereafter, we establish the corresponding linearized form of (7). For this, one has the following map in order to construct the corresponding linearized form:(8)β,βn,…,βn−k,γ,γn,…,γn−k,δ,δn,…,δn−k↦xn+1,xn,…,xn−k+1,yn+1,yn,…,yn−k+1,zn+1,zn,…,zn−k+1,where(9)β=α1xn−kα2+α3∏i=0kyn−i,βn=xn,…,βn−k=xn−k+1,γ=α4yn−kα5+α6∏i=0kzn−i,γn=yn,…,γn−k=yn−k+1,δ=α7zn−kα8+α9∏i=0kxn−i,δn=zn,…,δn−k=zn−k+1.Finally,J∧ about ∧ under map (8) becomes(10)J∧=00…0α1α2+α3yk+110…00⋮⋮⋱⋮⋮00…1000…0000…00⋮⋮⋱⋮⋮00…00−α7α9xkzα8+α9xk+12−α7α9xkzα8+α9xk+12…−α7α9xkzα8+α9xk+12−α7α9xkzα8+α9xk+1200…00⋮⋮⋱⋮⋮00…00−α1α3xykα2+α3yk+12−α1α3xykα2+α3yk+12…−α1α3xykα2+α3yk+12−α1α3xykα2+α3yk+1200…00⋮⋮⋱⋮⋮00…0000…0α4α5+α6zk+110…00⋮⋮⋱⋮⋮00…1000…0000…00⋮⋮⋱⋮⋮00…0000…0000…00⋮⋮⋱⋮⋮00…00−α4α6yzkα5+α6zk+12−α4α6yzkα5+α6zk+12…−α4α6yzkα5+α6zk+12−α4α6yzkα5+α6zk+1200…00⋮⋮⋱⋮⋮00…0000…0α7α8+α9xk+110…00⋮⋮⋱⋮⋮00…10.


## 3. Local Dynamics about Equilibria

By Theorem 1.5 of [1], the detailed local stability analysis about boundary equilibria P000, Px0z, Pxy0, P0yz, Px00, P00z, and P0y0 and the positive equilibrium point Pxyz will be investigated in this section.


### 3.1. Local Dynamics about Boundary Points

Theorem 1.
P000 of system (7) is a sink if(11)α1α2<1,α4α5<1,α7α8<1.Proof.
About the equilibrium pointP000, (10) takes the following form:(12)JP000=00…0α1α200…0000…0010…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…1000…0000…0000…0000…0α4α500…0000…0010…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…1000…0000…0000…0000…0α7α800…0000…0010…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…10.
Now, ifν1,ν2,…,νk+1,…,ν2k+1,…,ν3k+3 denotes characteristic roots of P000, the diagonal matrix Ω=diagρ1,ρ2,…,ρk+1,…,ρ2k+1,…,ρ3k+3, where(13)ρ1=ρk+2=ρ2k+3=1,ρ1+ς=ρk+2+ς=ρ2k+3+ς=1−ςε,1≤ς≤k,0<ε<1,(14)0<ε<min1k1−α1α2,1k1−α4α5,1k1−α7α8.
Now,(15)ΩJP000Ω−1=00…0ρ1ρk+1−1α1α200…00ρ2ρ1−10…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…ρk+1ρk−1000…0000…0000…0ρk+2ρ2k+2−1α4α500…00ρk+3ρk+2−10…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…ρ2k+2ρ2k+1−1000…0000…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…0000…00⋮⋮⋱⋮⋮00…0000…0000…00⋮⋮⋱⋮⋮00…0000…0ρ2k+3ρ3k+3−1α7α8ρ3k+2ρ3k+1−10…00⋮⋮⋱⋮⋮00…ρ3k+3ρ3k+2−10.
So,(16)0<ρk+1<⋯<ρ2<ρ1,0<ρ2k+2<⋯<ρk+2,0<ρ3k+3<⋯<ρ2k+3.
From (16), one obtains(17)ρ2ρ1−1<1,…,ρk+1ρk−1<1,ρk+3ρk+2−1<1,…,ρ2k+2ρ2k+1−1<1,ρ2k+4ρ2k+3−1<1,…,ρ3k+3ρ3k+2−1<1.
From (13) and (14), one obtains(18)ρ1ρk+1−1α1α2=α1α211−kε<1,ρk+2ρ2k+2−1α4α5=α4α511−kε<1,ρ2k+3ρ3k+3−1α7α8=α7α811−kε<1.
In view of (17) and (18), one obtains(19)max1≤m≤3k+3ρm=ΩJP000Ω−1=maxρ1ρk+1−1α1α2,ρ2ρ1−1,…,ρk+1ρk−1,ρk+3ρk+2−1,…,ρ2k+2ρ2k+1−1,ρ3k+2ρ3k+1−1,…,ρ3k+3ρ3k+2−1,ρk+2ρ2k+2−1α4α5,ρ2k+3ρ3k+3−1α7α8<1.
From (19), one gets the required statement.Theorem 2.
Px0z of (7) is unstable.Proof.
About the equilibrium pointPx0z, (10) takes the following form:(20)J|Px0z=00…0α1α200…0000…0010…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…1000…0000…0000…0000…0100…0000…0010…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…1000…00ΓΓ…ΓΓ00…0000…0100…0000…0010…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…10,where(21)Γ=−α9α7k+1α4−α5α6kk+1α7−α8α9.
Now, ifν1,ν2,…,νk+1,…,ν2k+1,…,ν3k+3 denotes characteristic roots of Px0z, the diagonal matrix Ω=diagρ1,ρ2,…,ρk+1,…,ρ2k+1,…,ρ3k+3, where (13) holds. Moreover,(22)Ψ=00…0ρ1ρk+1−1α1α200…00ρ2ρ1−10…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…ρk+1ρk−1000…0000…0000…0ρk+2ρ2k+2−100…00ρk+3ρk+2−10…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…ρ2k+2ρ2k+1−10ρ2k+3ρ1−1Γρ2k+3ρ2−1Γ…ρ2k+3ρk−1Γρ2k+3ρk+1−1Γ00…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…0000…00⋮⋮⋱⋮⋮00…0000…0000…00⋮⋮⋱⋮⋮00…0000…0ρ2k+3ρ3k+3−1ρ3k+2ρ3k+1−10…00⋮⋮⋱⋮⋮00…ρ3k+3ρ3k+2−10,where Φ=ΩJ|Px0zΩ−1. From (16) and (17), we obtain(23)ρ1ρk+1−1α1α2=α1α211−kε<1.
From (23), one gets ε<1−α1/α2. But ρk+2ρ2k+2−1=ρ2k+2−1=1/1−kε>1, and hence Px0z of (7) is unstable.
Similarly, local dynamics aboutPxy0 and P0yz of system under consideration can be summarized as follows.Theorem 3.
(i)
Pxy0 of system (7) is unstable.(ii)
P0yz of system (7) is unstable.Proof.
Same as proof of Theorems1 and 2.Theorem 4.
Px00 of system (7) is unstable.Proof.
About the equilibrium pointPx00, (10) takes the following form:(24)JPx00=00…0α1α200…0000…0010…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…1000…0000…0000…0000…0α4α500…0000…0010…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…1000…0000…0000…0000…0100…0000…0010…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…10.
Now, ifν1,ν2,…,νk+1,…,ν2k+1,…,ν3k+3 denotes characteristic roots of Px00, the diagonal matrix Ω=diagρ1,ρ2,…,ρk+1,…,ρ2k+1,…,ρ3k+3, where (13) holds. Moreover,(25)ΩJPx00Ω−1=00…0ρ1ρk+1−1α1α200…00ρ2ρ1−10…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…ρk+1ρk−1000…0000…0000…0ρk+2ρ2k+2−1α4α500…00ρk+3ρk+2−10…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…ρ2k+2ρ2k+1−1000…0000…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…0000…00⋮⋮⋱⋮⋮00…0000…0000…00⋮⋮⋱⋮⋮00…0000…0ρ2k+3ρ3k+3−1ρ3k+2ρ3k+1−10…00⋮⋮⋱⋮⋮00…ρ3k+3ρ3k+2−10.
From (16) and (17), we obtain(26)ρ2k+3ρ3k+3−1=ρ3k+3−1=11−kε>1,which gives Px00 of (7) is unstable.
Similarly, local dynamics aboutP00z,P0y0 of system under consideration can be summarized as follows.Theorem 5.
(i)
P00z of system (7) is unstable.(ii)
P0y0 of system (7) is unstable.Proof.
Same as proof of Theorems1, 2, and 4.


### 3.2. Local Dynamics about Positive Point

Theorem 6.
Pxyzof (7) is unstable.Proof.
About the equilibrium pointPxyz, (10) takes the following form:(27)JPxyz=00…01Λ1Λ1…Λ1Λ100…0010…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…1000…0000…0000…0000…01Λ2Λ2…Λ2Λ200…0010…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…1000…00ΓΓ…ΓΓ00…0000…0100…0000…0010…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…10.where(28)Λ1=−α3α1kk+1α1−α2α3k+1α7−α8α9,Λ2=−α6α4k+1α1−α2α3kk+1α4−α5α6,and Γ is depicted in (21). Moreover,(29)χ=00…0ρ1ρk+1−1ρ1ρk+2−1Λ1ρ1ρk+3−1Λ1…ρ1ρ2k+1−1Λ1ρ1ρ2k+2−1Λ1ρ2ρ1−10…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…ρk+1ρk−1000…0000…0000…0ρk+2ρ2k+2−100…00ρk+3ρk+2−10…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…ρ2k+2ρ2k+1−10ρ2k+3ρ1−1Γρ2k+3ρ2−1Γ…ρ2k+3ρk−1Γρ2k+3ρk+1−1Γ00…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…0000…00⋮⋮⋱⋮⋮00…00ρk+2ρ2k+3−1Λ2ρk+2ρ2k+4−1Λ2…ρk+2ρ3k+2−1Λ2ρk+2ρ3k+3−1Λ200…00⋮⋮⋱⋮⋮00…0000…0ρ2k+3ρ3k+3−1ρ3k+2ρ3k+1−10…00⋮⋮⋱⋮⋮00…ρ3k+3ρ3k+2−10,where χ=ΩJPxyzΩ−1. From (16) and (17), we obtain(30)ρ1ρk+1−1+ρ1ρk+2−1+ρk+3−1+⋯+ρ2k+1−1+ρ2k+2−1Λ1=11−kε+1+11−ε+⋯+11−k−1ε+11−kεΛ1>1,ρk+2ρ2k+2−1+ρk+2ρ2k+3−1+ρ2k+4−1+⋯ρ3k+2−1+ρ3k+3−1Λ2=11−kε+1+11−ε+⋯+11−k−1ε+11−kεΛ2>1,ρ2k+3ρ3k+3−1+ρ2k+3ρ1−1+ρ2−1+⋯+ρk−1+ρk+1−1Γ=11−kε+1+11−ε+⋯+11−k−1ε+11−kεΛ2>1.From (30), one can conclude that Pxyz of (7) is unstable.


## 3.1. Local Dynamics about Boundary Points

Theorem 1.
P000 of system (7) is a sink if(11)α1α2<1,α4α5<1,α7α8<1.Proof.
About the equilibrium pointP000, (10) takes the following form:(12)JP000=00…0α1α200…0000…0010…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…1000…0000…0000…0000…0α4α500…0000…0010…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…1000…0000…0000…0000…0α7α800…0000…0010…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…10.
Now, ifν1,ν2,…,νk+1,…,ν2k+1,…,ν3k+3 denotes characteristic roots of P000, the diagonal matrix Ω=diagρ1,ρ2,…,ρk+1,…,ρ2k+1,…,ρ3k+3, where(13)ρ1=ρk+2=ρ2k+3=1,ρ1+ς=ρk+2+ς=ρ2k+3+ς=1−ςε,1≤ς≤k,0<ε<1,(14)0<ε<min1k1−α1α2,1k1−α4α5,1k1−α7α8.
Now,(15)ΩJP000Ω−1=00…0ρ1ρk+1−1α1α200…00ρ2ρ1−10…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…ρk+1ρk−1000…0000…0000…0ρk+2ρ2k+2−1α4α500…00ρk+3ρk+2−10…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…ρ2k+2ρ2k+1−1000…0000…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…0000…00⋮⋮⋱⋮⋮00…0000…0000…00⋮⋮⋱⋮⋮00…0000…0ρ2k+3ρ3k+3−1α7α8ρ3k+2ρ3k+1−10…00⋮⋮⋱⋮⋮00…ρ3k+3ρ3k+2−10.
So,(16)0<ρk+1<⋯<ρ2<ρ1,0<ρ2k+2<⋯<ρk+2,0<ρ3k+3<⋯<ρ2k+3.
From (16), one obtains(17)ρ2ρ1−1<1,…,ρk+1ρk−1<1,ρk+3ρk+2−1<1,…,ρ2k+2ρ2k+1−1<1,ρ2k+4ρ2k+3−1<1,…,ρ3k+3ρ3k+2−1<1.
From (13) and (14), one obtains(18)ρ1ρk+1−1α1α2=α1α211−kε<1,ρk+2ρ2k+2−1α4α5=α4α511−kε<1,ρ2k+3ρ3k+3−1α7α8=α7α811−kε<1.
In view of (17) and (18), one obtains(19)max1≤m≤3k+3ρm=ΩJP000Ω−1=maxρ1ρk+1−1α1α2,ρ2ρ1−1,…,ρk+1ρk−1,ρk+3ρk+2−1,…,ρ2k+2ρ2k+1−1,ρ3k+2ρ3k+1−1,…,ρ3k+3ρ3k+2−1,ρk+2ρ2k+2−1α4α5,ρ2k+3ρ3k+3−1α7α8<1.
From (19), one gets the required statement.Theorem 2.
Px0z of (7) is unstable.Proof.
About the equilibrium pointPx0z, (10) takes the following form:(20)J|Px0z=00…0α1α200…0000…0010…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…1000…0000…0000…0000…0100…0000…0010…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…1000…00ΓΓ…ΓΓ00…0000…0100…0000…0010…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…10,where(21)Γ=−α9α7k+1α4−α5α6kk+1α7−α8α9.
Now, ifν1,ν2,…,νk+1,…,ν2k+1,…,ν3k+3 denotes characteristic roots of Px0z, the diagonal matrix Ω=diagρ1,ρ2,…,ρk+1,…,ρ2k+1,…,ρ3k+3, where (13) holds. Moreover,(22)Ψ=00…0ρ1ρk+1−1α1α200…00ρ2ρ1−10…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…ρk+1ρk−1000…0000…0000…0ρk+2ρ2k+2−100…00ρk+3ρk+2−10…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…ρ2k+2ρ2k+1−10ρ2k+3ρ1−1Γρ2k+3ρ2−1Γ…ρ2k+3ρk−1Γρ2k+3ρk+1−1Γ00…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…0000…00⋮⋮⋱⋮⋮00…0000…0000…00⋮⋮⋱⋮⋮00…0000…0ρ2k+3ρ3k+3−1ρ3k+2ρ3k+1−10…00⋮⋮⋱⋮⋮00…ρ3k+3ρ3k+2−10,where Φ=ΩJ|Px0zΩ−1. From (16) and (17), we obtain(23)ρ1ρk+1−1α1α2=α1α211−kε<1.
From (23), one gets ε<1−α1/α2. But ρk+2ρ2k+2−1=ρ2k+2−1=1/1−kε>1, and hence Px0z of (7) is unstable.
Similarly, local dynamics aboutPxy0 and P0yz of system under consideration can be summarized as follows.Theorem 3.
(i)
Pxy0 of system (7) is unstable.(ii)
P0yz of system (7) is unstable.Proof.
Same as proof of Theorems1 and 2.Theorem 4.
Px00 of system (7) is unstable.Proof.
About the equilibrium pointPx00, (10) takes the following form:(24)JPx00=00…0α1α200…0000…0010…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…1000…0000…0000…0000…0α4α500…0000…0010…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…1000…0000…0000…0000…0100…0000…0010…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…10.
Now, ifν1,ν2,…,νk+1,…,ν2k+1,…,ν3k+3 denotes characteristic roots of Px00, the diagonal matrix Ω=diagρ1,ρ2,…,ρk+1,…,ρ2k+1,…,ρ3k+3, where (13) holds. Moreover,(25)ΩJPx00Ω−1=00…0ρ1ρk+1−1α1α200…00ρ2ρ1−10…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…ρk+1ρk−1000…0000…0000…0ρk+2ρ2k+2−1α4α500…00ρk+3ρk+2−10…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…ρ2k+2ρ2k+1−1000…0000…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…0000…00⋮⋮⋱⋮⋮00…0000…0000…00⋮⋮⋱⋮⋮00…0000…0ρ2k+3ρ3k+3−1ρ3k+2ρ3k+1−10…00⋮⋮⋱⋮⋮00…ρ3k+3ρ3k+2−10.
From (16) and (17), we obtain(26)ρ2k+3ρ3k+3−1=ρ3k+3−1=11−kε>1,which gives Px00 of (7) is unstable.
Similarly, local dynamics aboutP00z,P0y0 of system under consideration can be summarized as follows.Theorem 5.
(i)
P00z of system (7) is unstable.(ii)
P0y0 of system (7) is unstable.Proof.
Same as proof of Theorems1, 2, and 4.


## 3.2. Local Dynamics about Positive Point

Theorem 6.
Pxyzof (7) is unstable.Proof.
About the equilibrium pointPxyz, (10) takes the following form:(27)JPxyz=00…01Λ1Λ1…Λ1Λ100…0010…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…1000…0000…0000…0000…01Λ2Λ2…Λ2Λ200…0010…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…1000…00ΓΓ…ΓΓ00…0000…0100…0000…0010…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…10.where(28)Λ1=−α3α1kk+1α1−α2α3k+1α7−α8α9,Λ2=−α6α4k+1α1−α2α3kk+1α4−α5α6,and Γ is depicted in (21). Moreover,(29)χ=00…0ρ1ρk+1−1ρ1ρk+2−1Λ1ρ1ρk+3−1Λ1…ρ1ρ2k+1−1Λ1ρ1ρ2k+2−1Λ1ρ2ρ1−10…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…ρk+1ρk−1000…0000…0000…0ρk+2ρ2k+2−100…00ρk+3ρk+2−10…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…ρ2k+2ρ2k+1−10ρ2k+3ρ1−1Γρ2k+3ρ2−1Γ…ρ2k+3ρk−1Γρ2k+3ρk+1−1Γ00…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…0000…00⋮⋮⋱⋮⋮00…00ρk+2ρ2k+3−1Λ2ρk+2ρ2k+4−1Λ2…ρk+2ρ3k+2−1Λ2ρk+2ρ3k+3−1Λ200…00⋮⋮⋱⋮⋮00…0000…0ρ2k+3ρ3k+3−1ρ3k+2ρ3k+1−10…00⋮⋮⋱⋮⋮00…ρ3k+3ρ3k+2−10,where χ=ΩJPxyzΩ−1. From (16) and (17), we obtain(30)ρ1ρk+1−1+ρ1ρk+2−1+ρk+3−1+⋯+ρ2k+1−1+ρ2k+2−1Λ1=11−kε+1+11−ε+⋯+11−k−1ε+11−kεΛ1>1,ρk+2ρ2k+2−1+ρk+2ρ2k+3−1+ρ2k+4−1+⋯ρ3k+2−1+ρ3k+3−1Λ2=11−kε+1+11−ε+⋯+11−k−1ε+11−kεΛ2>1,ρ2k+3ρ3k+3−1+ρ2k+3ρ1−1+ρ2−1+⋯+ρk−1+ρk+1−1Γ=11−kε+1+11−ε+⋯+11−k−1ε+11−kεΛ2>1.From (30), one can conclude that Pxyz of (7) is unstable.


## 4. Boundedness

The boundedness of positive solution of (7) is investigated in this section, as follows.Theorem 7.
IfΩnn=−k∞=xn,yn,znn=−k∞ is a positive solution of (7), then following holds for υ≥0:(31)0<xn≤α1α2υ+1x−k,ifn=k+1υ+1,α1α2υ+1x−k+1,ifn=k+1υ+2,⋮α1α2υ+1x−1,ifn=k+1υ+k,α1α2υ+1x0,ifn=k+1υ+k+1,(32)0<yn≤α4α5υ+1y−k,ifn=k+1υ+1,α4α5υ+1y−k+1,ifn=k+1υ+2,⋮α4α5υ+1y−1,ifn=k+1υ+k,α4α5υ+1y0,ifn=k+1υ+k+1,(33)0<zn≤α7α8υ+1z−k,ifn=k+1υ+1,α7α8υ+1z−k+1,ifn=k+1υ+2,⋮α7α8υ+1z−1,ifn=k+1υ+k,α7α8υ+1z0,ifn=k+1υ+k+1.Proof.
Forυ=0, (31)–(33) are true trivially. Suppose that (31)–(33) are true for υ=σ≥1, that is,(34)0<xn≤α1α2σ+1x−k,ifn=k+1σ+1,α1α2σ+1x−k+1,ifn=k+1σ+2,⋮α1α2σ+1x−1,ifn=k+1σ+k,α1α2σ+1x0,ifn=k+1σ+k+1,0<yn≤α4α5σ+1y−k,ifn=k+1σ+1,α4α5σ+1y−k+1,ifn=k+1σ+2,⋮α4α5σ+1y−1,ifn=k+1σ+k,α4α5σ+1y0,ifn=k+1σ+k+1,0<zn≤α7α8σ+1z−k,ifn=k+1σ+1,α7α8σ+1z−k+1,ifn=k+1σ+1,⋮α7α8σ+1z−1,ifn=k+1σ+k,α7α8σ+1z0,ifn=k+1σ+k+1.
Finally forυ=σ+1 and using (7), one obtains(35)xk+1σ+1+1=α1xk+1σ+1α2+α3∏i=0kyk+1σ+k+1−i≤α1α2xk+1σ+1≤α1α2σ+2x−k,xk+1σ+1+2=α1xk+1σ+2α2+α3∏i=0kyk+1σ+k+2−i≤α1α2xk+1σ+2≤α1α2σ+2x−k+1,⋮xk+1σ+1+k=α1xk+1σ+kα2+α3∏i=0kyk+1σ+2k−i≤α1α2xk+1σ+k≤α1α2σ+2x−1,xk+1σ+1+k+1=α1xk+1σ+k+1α2+α3∏i=0kyk+1σ+2k+1−i≤α1α2xk+1σ+k+1≤α1α2σ+2x0,(36)yk+1σ+1+1=α4yk+1σ+1α5+α6∏i=0kzk+1σ+k+1−i≤α4α5yk+1σ+1≤α4α5σ+2y−k,yk+1σ+1+2=α4yk+1σ+2α5+α6∏i=0kzk+1σ+k+2−i≤α4α5yk+1σ+2≤α4α5σ+2y−k+1,⋮yk+1σ+1+k=α4yk+1σ+kα5+α6∏i=0kzk+1σ+2k−i≤α4α5yk+1σ+k≤α4α5σ+2y−1,yk+1σ+1+k+1=α4yk+1σ+k+1α5+α6∏i=0kzk+1σ+2k+1−i≤α4α5yk+1σ+k+1≤α4α5σ+2y0,(37)zk+1σ+1+1=α7zk+1σ+1α8+α9∏i=0kxk+1σ+k+1−i≤α7α8zk+1σ+1≤α7α8σ+2z−k,zk+1σ+1+2=α7zk+1σ+2α8+α9∏i=0kxk+1σ+k+2−i≤α7α8zk+1σ+2≤α7α8σ+2z−k+1,⋮zk+1σ+1+k=α7zk+1σ+kα8+α9∏i=0kxk+1σ+2k−i≤α7α8zk+1σ+k≤α7α8σ+2z−1,zk+1σ+1+k+1=α7zk+1σ+k+1α8+α9∏i=0kxk+1σ+2k+1−i≤α7α8zk+1σ+k+1≤α7α8σ+2z0.
From (37)–(39), one has required result.Lemma 2.
If (11) holds, then Ωnn=−k∞ of (7) is bounded.Proof.
It is consequence of Theorem7.Remark 1.
From above theoretical studies, we conclude that no equilibria are locally stable exceptP000, and hence in the following section we will investigate that it is globally asymptotically stable under certain parametric conditions.


## 5. Global Dynamics aboutP000

Theorem 8.
P000 of system (7) is globally asymptotically stable if conditions, which are depicted in (11), hold.Proof.
By Theorem1, P000 of (7) is a sink and also Ωnn=−k∞ is bounded by Lemma 2. In order to get the required statement, it is adequate to prove xn,yn,zn is nonincreasing. In view of (11), from (7), one gets xn+1<xn−k, which again gives xk+1n+1<xk+1n−k and xk+1n+k+2<xk+1n+1. Thus, xn is decreasing because xk+1n+1,xk+1n+2,…,xk+1n+k+1 are nonincreasing. Similarly, it is easy to investigate that yn and zn are also decreasing. This implies that limn⟶∞xn,yn,zn=P000.Lemma 3.
If(38)α1α2>1,α4α5>1,α7α8>1,then for Ωnn=−k∞ of (7), the following holds:(i)
xn⟶∞ if yn⟶0.(ii)
yn⟶∞ if zn⟶0.(iii)
zn⟶∞ if xn⟶0.Lemma 4.
If (38) holds, then for ℓ=k,−k+1,…,1,0, one gets the following invariant intervals:(i)
x−ℓ,y−ℓ,z−ℓ∈0,k+1α7−α8/α9×0,k+1α1−α2/α3×k+1α4−α5/α6,∞⇒xn,yn,zn∈0,k+1α7−α8/α9×0,k+1α1−α2/α3×k+1α4−α5/α6,∞,∀n≥1.(ii)
x−ℓ,y−ℓ,z−ℓ∈k+1α7−α8/α9,∞×k+1α1−α2/α3,∞×0,k+1α4−α5/α6⇒xn,yn,zn∈k+1α7−α8/α9,∞×k+1α1−α2/α3,∞×0,k+1α4−α5/α6,∀n≥1.Proof.
This follows by induction.


## 6. Convergence Rate

Theorem 9.
If (11) is true, then error vector(39)ϱn=ϱn1⋮ϱn−k1ϱn2⋮ϱn−k2ϱn3⋮ϱn−k3,of positive solution Ωnn=−k∞ of (7) satisfies the following relations:(40)limn⟶∞nϱn=νJP000,limn⟶∞ϱn+1ϱn=νJP000,where νJP000 is equivalent to modulus of one of the characteristic roots of JP000 calculated at P000.Proof.
IfΩnn=−k∞ is an arbitrary solution of (7) such that limn⟶∞xn=x,limn⟶∞yn=y and limn⟶∞zn=z, then(41)xn+1−x=α1xn−kα2+α3∏i=0kyn−i−α1xα2+α3yk+1=α1α2+α3∏i=0kyn−ixn−k−x−α1α3x∏i=1kyn−iα2+α3∏i=0kyn−iα2+α3yk+1yn−y−α1α3xy∏i=2kyn−iα2+α3∏i=0kyn−iα2+α3yk+1yn−1−y−⋯−α1α3xyk−1yn−kα2+α3∏i=0kyn−iα2+α3yk+1yn−k+1−y−α1α3xykα2+α3∏i=0kyn−iα2+α3yk+1yn−k−y.
Similarly from (7), we obtain(42)yn+1−y=α4α5+α6∏i=0kzn−iyn−k−y−α4α6y∏i=1kzn−iα5+α6∏i=0kzn−iα5+α6zk+1zn−z−α4α6yz∏i=2kzn−iα5+α6∏i=0kzn−iα5+α6zk+1zn−1−z−⋯−α4α6yzk−1zn−kα5+α6∏i=0kzn−iα5+α6zk+1zn−k+1−z−α4α6yzkα5+α6∏i=0kzn−iα5+α6zk+1zn−k−z,zn+1−z=−α7α9z∏i=1kxn−iα8+α9∏i=0kxn−iα8+α9zk+1xn−x−α7α9zx∏i=2kxn−iα8+α9∏i=0kxn−iα8+α9xk+1xn−1−x−⋯−α7α9zxk−1zn−kα8+α9∏i=0kxn−iα8+α9xk+1xn−k+1−x−α7α9zxkα8+α9∏i=0kxn−iα8+α9xk+1xn−k−x+α7α8+α9∏i=0kxn−izn−k−z.
Now, set(43)ϱn1=xn−x,ϱn2=yn−y,ϱn3=zn−z.
Utilizing (43) in (41) and (42), we get(44)ϱn+11=Ann−kϱn−k1+Bnnϱn2+Bnn−1ϱn−12+⋯+Bnn−k+1ϱn−k+12+Bnn−kϱn−k2,ϱn+12=Cnn−kϱn−k2+Dnnϱn3+Dnn−1ϱn−13+⋯+Dnn−k+1ϱn−k+13+Dnn−kϱn−k3,ϱn+13=Ennϱn1+Enn−1ϱn−11+⋯+Enn−k+1ϱn−k+11+Enn−kϱn−k1+Fnn−kϱn−k3,where(45)Ann−k=α1α2+α3∏i=0kyn−i,Bnn=−α1α3x∏i=1kyn−iα2+α3∏i=0kyn−iα2+α3yk+1,Bnn−1=−α1α3xy∏i=2kyn−iα2+α3∏i=0kyn−iα2+α3yk+1,⋮Bnn−k+1=−α1α3xyk−1yn−kα2+α3∏i=0kyn−iα2+α3yk+1,Bnn−k=−α1α3xykα2+α3∏i=0kyn−iα2+α3yk+1,Cnn−k=α4α5+α6∏i=0kzn−i,Dnn=−α4α6y∏i=1kzn−iα5+α6∏i=0kzn−iα5+α6zk+1,Dnn−1=−α4α6yz∏i=2kzn−iα5+α6∏i=0kzn−iα5+α6zk+1,(46)⋮Dnn−k+1=−α4α6yzk−1zn−kα5+α6∏i=0kzn−iα5+α6zk+1,Dnn−k=−α4α6yzkα5+α6∏i=0kzn−iα5+α6zk+1,Enn=−α7α9z∏i=1kxn−iα8+α9∏i=0kxn−iα8+α9zk+1,Enn−1=−α7α9zx∏i=2kxn−iα8+α9∏i=0kxn−iα8+α9xk+1,⋮Enn−k+1=−α7α9zxk−1zn−kα8+α9∏i=0kxn−iα8+α9xk+1,Enn−k=−α7α9zxkα8+α9∏i=0kxn−iα8+α9xk+1,Fnn−k=α7α8+α9∏i=0kxn−i.
From (45) and (46), we get(47)limn⟶∞Ann−k=α1α2+α3yk+1,limn⟶∞Bnn=limn⟶∞Bnn−1=⋯=limn⟶∞Bnn−k+1=limn⟶∞Bnn−k=−α1α3xykα2+α3yk+12,limn⟶∞Cnn−k=α4α5+α6zk+1,limn⟶∞Dnn=limn⟶∞Dnn−1=⋯=limn⟶∞Dnn−k+1=limn⟶∞Dnn−k=−α4α6yzkα5+α6zk+12,limn⟶∞Enn=limn⟶∞Enn−1=⋯=limn⟶∞Enn−k+1=limn⟶∞Enn−k=−α7α9xkzα8+α9xk+12,limn⟶∞Fnn−k=α7α8+α9xk+1.
So, we have the limiting system [26]:(48)ϱn+1=Kϱn,where ϱn is depicted in (39) and K is same as J|∧ about ∧. In particular, about P000 it becomes(49)ϱn+11⋮ϱn−k+11ϱn+12⋮ϱn−k+12ϱn+13⋮ϱn−k+13=00…0α1α200…0000…0010…0000…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…1000…0000…0000…0000…0α4α500…0000…0010…0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…1000…0000…0000…0000…0α7α800…0000…0010…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00…0000…0000…10ϱn1⋮ϱn−k1ϱn2⋮ϱn−k2ϱn3⋮ϱn−k3,which is same as linearized system of (7) about P000.


## 7. Numerical Simulations

We will present numerous numerical examples, which represent different dynamical behaviors of system (7), in this section.Example 1.
Figure1 indicates the dynamics of (7) about P000 if k=2, αℓℓ=1,2,…,9, respectively, are 14,15,3,6,5,1,11,13,and3 and initial values x−ℓ, ℓ−ℓ, z−ℓℓ=2,1,0, respectively, are 0.07, 0.9, 0.02, 0.9, 0.9, 1.4, 4.9, 0.9, and 0.8. More precisely, Figures 1(a)–1(c) show that P000 of (7) is a sink while its corresponding global attractor is shown is Figure 1(d). Figure 1(d) clearly indicates that all solutions go towards P000; this means that P000 is globally asymptotically stable. Hence, this simulation agrees with the obtained results.Figure 1
Dynamics of (7).
(a)(b)(c)(d)Example 2.
Figure2 indicates the same dynamics of (7) about P000 if k=4, αℓℓ=1,2,…,9, respectively, are 14, 15, 33, 6, 5, 2, 11, 13, and 3 and initial values x−ℓ, ℓ−ℓ, z−ℓℓ=4,3,2,1,0, respectively, are 0.07, 0.9, 0.07, 0.9, 0.02, 0.9, 0.9, 0.9, 0.9, 1.4, 4.9, 1.4, 4.9, 0.9, and 0.8.Figure 2
Dynamics of (7).
(a)(b)(c)(d)Example 3.
Figure3 indicates the dynamics of (7) about P000 if k=6, αℓℓ=1,2,…,9, respectively, are 14, 15, 33, 6, 5, 2, 11, 13, and 3 and initial values x−ℓ, ℓ−ℓ, z−ℓℓ=6,5,4,3,2,1,0, respectively, are 0.7, 1.9, 0.07, 0.9, 0.07, 0.9, 0.02, 1.2, 2.9, 0.9, 1.4, 2.9, 0.29, 4.9, 1.4, 4.9, 0.9, and 0.8.Figure 3
Dynamics of (7).
(a)(b)(c)(d)


## 8. Concluding Remarks

The presented work is a natural extension of [21–25]. In this work, we have investigated the global dynamics of higher-order rational system (7) which is our main contribution to dynamical system. We explored that system (7) has eight equilibria P000,Px0z,Pxy0,P0yz,Px00,P00z,P0y0, and Pxyz under certain conditions on parameters αℓℓ=1,2,…,9. We investigated that only P000 is locally stable under the condition, which is depicted in (11), while the rest of the equilibria Px0z,Pxy0,P0yz,Px00,P00z,P0y0,andPxyz are locally unstable. Further, we have studied the global dynamics, boundedness, and rate of convergence of positive solution of (7). Finally, we have verified theoretical results numerically [2].



---

*Source: 3696874-2020-07-15.xml*