# Stability, Neimark–Sacker Bifurcation, and Approximation of the Invariant Curve of Certain Homogeneous Second-Order Fractional Difference Equation

Authors: Mirela Garić-Demirović; Samra Moranjkić; Mehmed Nurkanović; Zehra Nurkanović
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

We investigate the local and global character of the unique equilibrium point and boundedness of the solutions of certain homogeneous fractional difference equation with quadratic terms. Also, we consider Neimark–Sacker bifurcations and give the asymptotic approximation of the invariant curve.


---

## Body

## 1. Introduction and Preliminaries

In this paper, the subject of our consideration is the following difference equation:(1)un+1=Aun2aun2+bunun−1+cun−12,n=0,1,2,…,with positive parameters A,a,b, and c where initial conditions u−1,u0 are positive numbers. By substituting xn=1/unn=0,1,2,…, equation (1) reduces to the following equation:(2)xn+1=axn−12+bxnxn−1+cxn2Axn−12,n=0,1,2,….Equations (1) and (2) are the special cases of the following homogeneous rational difference equation:(3)xn+1=Axn2+Bxnxn−1+Cxn−12axn2+bxnxn−1+cxn−12,n=0,1,2,….Equation (3) is a general homogeneous rational difference equation with quadratic terms and it is very complicated for investigation in many special cases. The function associated with the right side of equation (3) is of the form fu,v=Au2+Buv+Cv2/au2+buv+cv2, with the following monotonicity property: it is either monotonically decreasing in the first variable and monotonically increasing in the second one or monotonically increasing in the first variable and monotonically decreasing in the second variable. Using the theory of monotone maps, it was possible to investigate the global dynamics of some special cases of equation (3) [1–4] in the situations when corresponding function f is monotonically decreasing in the first variable and monotonically increasing in the second variable.In [4], the authors investigated the local and global character of the unique equilibrium point and the existence of Neimark–Sacker and period-doubling bifurcations of equation (3). They have also studied the local and global stability of the minimal period-two solutions for some special cases of the parameters. The special case when A=B=0 and C=1, i.e.,(4)xn+1=xn−12axn2+bxnxn−1+cxn−12,n=0,1,2,…,was considered in [2]. It was shown that equation (4) is characterized by three types of global behavior with respect to the existence of a unique positive equilibrium and existence of one or two minimal period-two solutions, one of which is locally asymptotically stable and the other is a saddle point. An important feature of this equation is the coexistence of an equilibrium and the minimal period-two solution which are both locally asymptotically stable. Also, the basins of attraction of these solutions are described in detail.In [3], a special case of equation (3) was considered when B=c=0, i.e.,(5)xn+1=Axn2+Cxn−12axn2+bxnxn−1,n=0,1,2,….It was shown that (5) exhibits period-two bifurcation and that stable manifolds of the minimal period-two solutions represent the boundaries of the basin of attraction of locally stable equilibrium point and the basins of attraction of points 0,∞ and ∞,0. Further, in the situations when equilibrium is a saddle point, corresponding stable manifold separates the basins of attraction of the points 0,∞ and ∞,0.The investigation of the special cases of equation (3) when corresponding function f is monotonically increasing in the first variable and monotonically decreasing in the second variable is significantly harder, since by now, there is no any general result for this type of monotonicity. Some initial steps, regarding this, were taken in [5], where the considered case was A=p, C=a=1, B=b=c=0, i.e.,(6)xn+1=p+xnxn−12,n=0,1,2,….In [5], the authors have successfully used the embedding method to demonstrate the boundedness of the solutions, and then they determined the invariant interval of equation (6). That was a crucial idea for proving that local asymptotic stability (which holds when p>1) implies global asymptotic stability of the unique positive equilibrium point when p>2. Also, the existence of Neimark–Sacker bifurcation is shown and asymptotic approximation of the invariant curve is computed.The investigation of the dynamics of equation (4) and its special cases has been the subject of many research studies for the last ten years. Some of these cases, as we have seen, have been successfully realized. However, only the ideas from [5] finally made the problem of investigating the behavior of equation (1) or equation (2) solvable. Namely, note that equation (2) has the form(7)xn+1=p+xnxn−1+qxnxn−12,n=0,1,2,…,where a/A=p, b/A=1, and c/A=q. Therefore, we will consider equation (7) instead of equation (1). However, the application of the embedding method on equation (7) is significantly harder compared with application on equation (6) (because the form is more complicated). We will investigate local and global stability of a unique equilibrium point and boundedness of the solutions of equation (7) and examine the existence of Neimark–Sacker bifurcation. Also, in the situation when Neimark–Sacker bifurcation appears, we will give the asymptotic approximation of the invariant curve.The special case of equation (3), when C=c=1 and B=b=0, was considered in [1]. In the region of parameters where a<A, corresponding function f is monotonically increasing by first variable and monotonically decreasing by second variable, and through very complicated calculations, using the so-called “M-m” theorems, it is shown that in some areas of parametric space of parameters A and a, unique positive equilibrium is globally stable. The existence of period-doubling bifurcation is proved in the case A<a, when corresponding function f is monotonically decreasing by first variable and monotonically increasing by second variable. Using the theory of monotone maps, the global stability of minimal period-two solution is shown for some special values of parameters.Notice that equation (3) is a special case (and probably the most complicated equation of the form (3, 3)) of the following general second-order rational difference equation with quadratic terms(8)xn+1=Axn2+Bxnxn−1+Cxn−12+Dxn+Exn−1+Faxn2+bxnxn−1+cxn−12+dxn+exn−1+f,n=0,1,2,…,that has caught the attention of mathematical researchers over the last ten years ([1–12]).The following lemma gives us the type of local stability of a unique positive equilibrium point of equation (7) depending on different values of parameters p and q.Lemma 1.
Equation (7) has a unique equilibrium point x¯=p+q+1, which is(a)
locally asymptotically stable ifp>q(b)
nonhyperbolic ifp=q(c)
a repeller ifp<qProof.
Denotegu,v=p+u/v+qu2/v2. The linearized equation associated with equation (7) about the equilibrium point has the form(9)zn+1=szn+tzn−1,where s=∂g/∂ux¯,x¯ and t=∂g/∂vx¯,x¯. Notice(10)s=−t=∂g∂ux¯,x¯=1+2qp+1+q>0.
Since(11)s<1−t<2⟺s<1+s<2⟺s<1⟺p>q,s=1−t∨t=−1∧s≤2⟺s=1+s∨s=1∧s≤2⟺s=1⟺p=q,s<1−t∧t>1⟺s<1+s∧s>1⟺p<q,the conclusion follows.This paper is organized as follows. In Section2, using the embedding method [5, 13] and the so-called “M-m” theorems [14–17], we prove global asymptotic stability of a unique positive equilibrium for p≥q+12q+1 and conduct the semicycle analysis as well. In Section 3, using Neimark–Sacker theorem [15, 18–21], we give reduction to the normal form and perform computation of the coefficients of the aforementioned bifurcation, based on the computational algorithm developed in [22]. Furthermore, we determine the asymptotic approximation of the invariant curve and give a visual evidence.


## 2. Global Asymptotic Stability

In this section, we will show that all solutions of equation (7) are bounded, and using the so-called “M-m” theorem, we will obtain sufficient conditions for unique positive equilibrium to be globally asymptotically stable. Similar to [5, 13], we apply the method of embedding. First, we substitute(12)xn=p+xn−1xn−2+qxn−1xn−22,in equation (7) and obtain(13)xn+1=p+pxn−1+1xn−2+qxn−1xn−22+qpxn−1+1xn−2+qxn−1xn−222.Then, by substituting(14)xn−1=p+xn−2xn−3+qxn−2xn−32,in equation (13), we have(15)xn+1=p+pxn−1+1xn−2+q2xn−32+pqxn−22+qxn−2xn−3+qpxn−1+1xn−2+q2xn−32+pqxn−22+qxn−2xn−32.Notice that the solutions of equation (15) are bounded:(16)p≤xn+1≤p+1+1p+q2p2+qp+qp2+q1+1p+q2p2+qp+qp22,n≥4,that is,(17)p≤xn+1≤p+q+1p+q+pq+p2+q2pq+p2+q2p4,n≥4.Since every solutionxn′n=−1∞ of equation (7) is also a solution xnn=−2∞ of equation (15) with initial values x−3=x−1′, x−2=x0′, x−1=p+x0′/x−1′+qx0′/x−1′2 and x0=p+x1′/x0′+qx1′/x0′2, we see that the solutions of equation (7) are also bounded.From (7) and (15), we get(18)xnxn−1=pxn−1+1xn−2+q2xn−32+pqxn−22+qxn−2xn−3,which implies(19)xn=p+xn−1xn−2+q2xn−1xn−32+pqxn−1xn−22+qxn−1xn−2xn−3,that is,(20)xn+1=p+xnxn−1+q2xnxn−22+pqxnxn−12+qxnxn−1xn−2.By replacingxn=p+xn−1/xn−2+qxn−1/xn−22 in (20), we obtain the following equation:(21)xn+1=p+xnxn−1+q2xnxn−22+pqp+xn−1/xn−2+qxn−1/xn−22xn−12+qxnxn−1xn−2,that is,(22)xn+1=p+xnxn−1+q2xnxn−22+p2qxn−12+pq2xn−22+pqxn−1xn−2+qxnxn−1xn−2.Remark 1.
Note that equation (22) has the unique equilibrium point x¯=p+q+1, which is the same equilibrium as in equation (7). It follows from(23)p+1+q2x+p2qx2+pq2x2+pqx2+px−x=p+q−x+1pq+qx+x2x2=0.
Furthermore, notice that every solutionxn′n=−1∞ of equation (7) is also a solution xnn=−2∞ of equation (22) with initial values x−2=x−1′, x−1=x0′ and x0=p+x0′/x−1′+qx0′/x−1′2 and that it is of the form xn+1=fxn,xn−1,xn−2, where(24)fu,v,w=p+uv+q2uw2+p2qv2+pq2w2+pqvw+quvw.Lemma 2.
Every intervalI of the form p,U, where(25)U=pq+pq+p2+q2p−q−1p+q,p>q+1,is an invariant interval for the function f.Proof.
As we know, forp<U, interval I=p,U is invariant for the function f if(26)∀u,v,w∈Ifu,v,w∈I.
Forp≤u,v,w≤U, we have that p≤fu,v,w≤p+U/p+q2U/p2+q+q2/p+q/p+qU/p2. If U satisfies(27)p+Up+q2Up2+q+q2p+qp+qUp2≤U,then we obtain fu,v,w≤U. It further implies that for every p>q+1, there exists such U, which means that I is invariant for the function f, where U≥pq+pq+p2+q2/p−q−1p+q. Since we can assume that p>q+1 and U=pq+pq+p2+q2/p−q−1p+q, I is the invariant for f.Lemma 3.
IntervalI=p,U is an attracting interval for equation (22).Proof.
It is clear that we need to show that every solution of equation (22) must enter interval I. Note that for arbitrary initial conditions x−2=x−1′, x−1=x0′, and x0=p+x0′/x−1′+qx0′/x−1′2, it holds(28)xn>pforn≥0.
Ifx1,x2,x3∈I, p>q+1, then xn∈I for n>3, by Lemma 2. Otherwise, if x3>pq+pq+p2+q2/p−q−1p+q, let us prove that there must be some k>3 such that xn∈I for all n≥k. Namely, suppose that x3>pq+pq+p2+q2/p−q−1p+q for arbitrary initial conditions x−2,x−1,x0. Then,(29)xn−3>p,xn−2>pforn≥3.
So, from equation (22), we obtain(30)xn=p+xn−1xn−2+q2xn−1xn−32+p2qxn−22+pq2xn−32+pqxn−2xn−3+qxn−1xn−2xn−3,which implies(31)xn<p+p2qp2+pq2p2+pqp2+xn−1p+q2xn−1p2+qxn−1p2,that is,(32)xn<q+pq+p2+q2p+p+q+q2p2xn−1.
By induction, we conclude(33)xn<q+pq+p2+q2p1−p+q+q2/p2n−31−p+q+q2/p2+x3p2/p+q+q2n−3forn>3.
Sincep2/p+q+q2>1 for p>q+1, the right side in (33) is decreasing sequence converging to(34)q+pq+p2+q2p11−p+q+q2/p2=pq+pq+p2+q2p−q−1p+q=U.
Then, from (33), we get(35)limn⟶∞xn≤U.
Since the caselimn⟶∞xn=U is not possible (otherwise there would exist another positive equilibrium different from x¯), it implies that there is some k>3 such that(36)p<xn<U,for all n≥k, i.e., every solution of equation (22) must enter the interval I.Theorem 1.
Ifp≥q+12q+1, then the equilibrium point x¯ of equation (7) is globally asymptotically stable.Proof.
It is enough to prove thatx¯ is an attractor of equation (22). Since there exists the invariant and attracting interval I=p,pq+pq+p2+q2/p−q−1p+q when p>q+1, we need to check the conditions of Theorem A.0.5 in [14]:(37)M=fM,m,mm=fm,M,M⟺M=p+p2qm2+Mm+qMm2+q2Mm2+pq2m2+pqm2,m=p+p2qM2+mM+qmM2+q2mM2+pq2M2+pqM2,i.e.,(38)M=fM,m,mm=fm,M,M⟺Mm2=pm2+p2q+Mm+qM+q2M+pq2+pq,mM2=pM2+p2q+mM+qm+q2m+pq2+pq,.
By subtracting the second equation in (38) from the first, we get(39)M−mq−Mp−mp+q2+Mm=0,from which m=M or for m≠M,(40)m−pM=mp−q−q2.
Ifm=p, then(41)0=p2−q−q2.
Sincep2−q−q2>0 for p>q+1, it implies that m≠p. Therefore,(42)M=mp−q−q2m−p,m≠p.
By substituting (42) into (38), we obtain the following quadratic equation:(43)m2p−q−1p+q−q2p−q+2pq+p2−1m+qq+p2q+p2+p3+2q2+q3=0,whose discriminant has the form(44)Dp,q=q3q−p2+2q2+1q+3p2q+4p3+3q2+2q3.
It is clear thatDp,q<0 if p>q+12q+1 and that there are no real solutions for equation (43). If p=q+12q+1, then m=M=p+q+1, which is a contradiction with assumption m≠M. When q+1<p<q+12q+1, equation (43) has two positive roots m±>p; it implies that M±>0 (for example, m=3.0032, M=952.12 for p=3, and q=1.99, where m,M is in invariant interval I). Therefore, the conditions of Theorem A.0.5 in [14] are satisfied for p≥q+12q+1 and every solution of equation (22) converges to x¯=p+q+1. By Remark 1, the unique equilibrium x¯=p+q+1 of equation (7) is an attractor. By using Lemma 1, we conclude that x¯ is globally asymptotically stable (see Figure 1).Figure 1
Local (a) and global (b) asymptotic stability inp,q plane.
(a)(b)Conjecture 1.
Ifp>q, then the equilibrium point x¯=p+q+1 is globally asymptotically stable.


## 3. Neimark–Sacker Bifurcations

The following results are obtained by applying the algorithm from Theorem1 and Corollary 1 in [5] (see also [22]). If we make a change of variable yn=xn−x¯, we will shift the equilibrium point to the origin. Then, the transformed equation is given by(45)yn+1=yn+p+1+qyn−1+p+1+q−1+qyn+p+1+q2yn−1+p+1+q2−1,n=0,1,….Set(46)un=yn−1,vn=yn,forn=0,1,…,and write equation (1) in the equivalent form:(47)un+1=vn,vn+1=vn+p+1+qun+p+1+q−1+qvn+p+1+q2un+p+1+q2−1.LetF denote the corresponding map defined by(48)Fuv=vv+p+1+qu+p+1+q−1+qv+p+1+q2u+p+1+q2−1.Then, the Jacobian matrix ofF is given by(49)JacFu,v=01−p+q+v+12qp+q+v+1+p+q+u+1p+q+u+132qp+q+v+1+p+q+u+1p+q+u+12.The eigenvalues ofJacF0,0 are μp and μp¯ where(50)μp=1+2q+i1+2q3+2q+4p2p+q+1,μp=1+2qp+q+1.Lemma 4.
Ifp=p0=q, then F has equilibrium point at 0,0 and eigenvalues of Jacobian matrix of F at 0,0 are μ and μ¯ where(51)μp0=12+i32.
Moreover,μ satisfies the following:(i)
μkp0≠1 for k=1,2,3,4(ii)
d/dpμp=dp0=−1/21+2q<0 at p=p0(iii)
Eigenvectors associated to theμ are(52)qp0=12−i32,1T,pp0=i3,163−i3,such that pA=μp, Aq=μq, and pq=1, where A=JacF0,0.Proof.
Letp=p0=q. Then, we obtain(53)A=JacF0,0=01−1+2qx¯1+2qx¯=01−11.
After straightforward calculation forp=p0=q, we obtain μp0=1 and(54)μp0=12+i32,μ2p0=−12+i32,μ3p0=−1,μ4p0=−12−i32,from which it follows that μkp0≠1 for k=1,2,3,4. Furthermore, we get(55)ddpμp=−1+2q2p+q+13/2,ddpμpp=p0=−121+2q<0.
It is easy to see thatpp0A=μpp0 and pp0qp0=1.Letp=p0+η, where η is a sufficiently small parameter. From Lemma 4, we can transform system (47) into the normal form(56)Fλ,x=ℱλ,x+Ox5,and there are smooth functions aλ,bλ and ωλ so that in polar coordinates, the function ℱλ,x is given by(57)rθ=μλr−aλr3θ+ωλ+bλr2.Now, we computeap0 following the procedure in [22]. Notice that p=p0 if and only if η=0. First, we compute K20,K11 and K02 defined in [22]. For p=p0=q, we have(58)Fuv=Auv+Guv,where(59)Guv≔0v+2q+1u+2q+1−1+qv+2q+1u+2q+12−1+u−v.Hence, forp=p0, system (47) is equivalent to(60)un+1vn+1=Aunvn+Gunvn.Define the basis ofℝ2 by Φ=q,q¯, where q=qp0=1/2−i3/2; then, we can represent u,v as(61)uv=Φzz¯=qz+q¯z¯=1+i32z¯+1−i32zz¯+z.Let(62)GΦzz¯=12g20z2+2g11zz¯+g02z¯2+Oz3.Now, we have(63)GΦzz¯=0z¯+z+2q+11+i3/2z¯+1−i3/2z+2q+1+qz¯+z+2q+11+i3/2z¯+1−i3/2z+2q+12−1−1−i32z¯−1+i32z−1.Since(64)g20=∂2∂z2GΦzz¯z=0=0−5q−2+i3q2q+12,g11=∂2∂z∂z¯GΦzz¯z=0=04q+12q+12,g02=∂2∂2z¯GΦzz¯z=0=0−5q−2−i3q2q+12,then(65)K20=μ2I−A−1g20=1+i3−5q−2+i3q42q+125q+2−i3q22q+12,K11=I−A−1g11=4q+12q+124q+12q+12,K02=μ¯2I−A−1g02=K¯20.By usingK20, K11, and K02, we have(66)g21=∂3∂z2∂z¯GΦzz¯+12K20z2+2K11zz¯+K02z¯2z=0=0−2i3q22q+14.Finally, we get(67)ap0=12Repg21μ¯=−q22q+14<0.Ifx¯ is fixed point of F, then invariant curve Γλ can be approximated by(68)x1x2≈x¯+2ρ0Reqeiθ+ρ02ReK20e2iθ+K11,where(69)d=ddημλλ=λ0,ρ0=−daη,θ∈ℝ.Thus, we prove the following result.Theorem 2.
Letx¯=p+q+1; then, there is a neighborhood U of the equilibrium point x¯ and δ>0 such that for p−q<δ and x0,x−1∈U, ω-limit set of solution of equation (7) is the equilibrium point if p>q and belongs to a closed invariant C1 curve Γp encircling the equilibrium point if p<q. Furthermore, Γq=0 and invariant curve Γp can be approximated by(70)x1x2≈1+p+q+q−p2q+132q2cost+3sint−p−q2q+18q−4q+1cos2t+32q+1sin2t+24q21+p+q+2q−p2q+13q2cost−p−q2q+13qsin2t+8q+5q+2cos2t+24q2.Sinceap0≠0, a nondegenerate Neimark–Sacker bifurcation occurs at the critical value p0=q. We proved ap0<0 and dp0<0, so it implies if p<p0, then an attracting closed curve exists, surrounding the unstable fixed point, when the parameter p crosses the bifurcation value p0 (supercritical Neimark–Sacker bifurcation). As p increases, the attracting closed curve decreases in size and merges with the fixed point at p=p0, leaving a stable fixed point (subcritical Neimark–Sacker bifurcation). All orbits starting outside or inside the closed invariant curve, except at the origin, tend to the attracting closed curve.The asymptotic approximation of the invariant curve is shown in Figure2, and some orbits and trajectories in the case where the unique equilibrium point is stable or nonhyperbolic are given in Figures 3 and 4. Using parameter value q=3 and decreasing the dynamical parameter p, the unique positive equilibrium point loses its stability via Neimark–Sacker bifurcation leading to chaos as depicted in Figure 5(a). Also, the Lyapunov exponents corresponding to Figure 5(a) are shown in Figure 5(b), which verifies the existence of chaos after the occurrence of Neimark–Sacker bifurcation.Figure 2
Trajectories (a) and invariant curve (b) forp=2.9 and q=3.
(a)(b)Figure 3
Orbits with initial conditions:x0,x−1=2.1,2.1 (a) and x0,x−1=6.8,7.3 (b), for p=2.9 and q=3.
(a)(b)Figure 4
Trajectories forp=3.2 and q=3 (a) and p=1.5 and q=1.5 (b).
(a)(b)Figure 5
Bifurcation diagrams inp,x plane (a) and corresponding Lyapunov exponent (b) for values of parameter q=3.
(a)(b)


## 4. Conclusion

The investigation of the dynamic stability of homogeneous difference equation (3) and all its special cases is very complicated. The corresponding function fu,v associated with the right side of equation (3) is monotonically increasing in the first variable and monotonically decreasing in the second variable or monotonically decreasing in the first variable and monotonically increasing in the second variable or it switches its type of monotonicity between the first and second case or vice versa, depending on the parameters which appear in the function.The theory of monotone maps (or more precisely, the theory of competitive maps) was used for determining the dynamics of equation (3) or some special case of equation (3), in the scenario when corresponding function fu,v is monotonically decreasing in the first variable and monotonically increasing in the second variable (see [2, 3]).The investigation of the special cases of equation (3) when corresponding function f is monotonically increasing in the first variable and monotonically decreasing in the second variable is significantly harder because there is no general result for this type of monotonicity.In every situation, when the corresponding equation does not possess minimal period-two solutions, the global stability of a unique equilibrium usually can be determined by applying the so-called “M-m” theorems and finding an invariant interval of the mapf before that, of course. However, it is very often impossible or extremely complicated to conduct, as we saw in [1]. That is the case with equation (1) which does not possess minimal period-two solution (since the corresponding function is monotonically increasing by first and monotonically decreasing by second variable). So, for that reason, we studied equation (7) instead of equation (1). By using the method of embedding, we were able to connect equation (7) with equation (15). Namely, we have shown that every solution xn′n=−1∞ of equation (7) is also a solution xnn=−2∞ of equation (15) with initial conditions x−3=x−1′, x−2=x0′, x−1=p+x0′/x−1′+qx0′/x−1′2, and x0=p+x1′/x0′+qx1′/x0′2. Furthermore, we have shown that every solution of equation (15) is bounded, which implies that every solution of equation (7) is also bounded. After that, we linked equation (15) with equation (22) and showed that every solution xn′n=−1∞ of equation (7) is also a solution of equation (22) with initial conditions x−2=x−1′, x−1=x0′, and x0=p+x0′/x−1′+qx0′/x−1′2. Additionally, we determined the invariant and attracting interval for the function that is associated with the right side of equation (22) and successfully applied “M-m” theorem to get conditions for parameters p and q under which the equilibrium of equation (22) and therefore of equation (7) is globally asymptotically stable. The area of the regions in p,q plane where unique equilibrium is locally asymptotically stable and is not globally asymptotically stable is small (see Figure 1). We expect to prove Conjecture 1 in some of our future studies.Finally, using Neimark–Sacker theorem [15, 18–21], we gave reduction to the normal form and performed computation of the coefficients of bifurcation, based on the computational algorithm developed in [22]. Furthermore, we determined the asymptotic approximation of the invariant curve and provided visual evidence (see Figures 2–5). Also, based on the computational algorithm in [23], we calculated the Lyapunov exponents corresponding to Figure 5(a) to confirm the existence of chaos after the occurrence of Neimark–Sacker bifurcation as in [24, 25] (see Figure 5(b)).



---

*Source: 6254013-2020-08-05.xml*