Arab. J. Math. (2023) 12:541–551
https://doi.org/10.1007/s40065-023-00438-9 Arabian Journal of Mathematics

Mensah Folly-Gbetoula

On a family of higher order recurrence relations:
symmetries, formula solutions, periodicity and stability
analysis

Received: 21 March 2023 / Accepted: 24 July 2023 / Published online: 10 August 2023
© The Author(s) 2023

Abstract In this paper, we present formula solutions of a family of difference equations of higher order. We
discuss the periodic nature of the solutions and we investigate the stability character of the equilibrium points.
We utilize Lie symmetry analysis as part of our approach together with some number theoretic functions. Our
findings generalize certain results in the literature.

Mathematics Subject Classification 39A10 · 39A13 · 39A99

1 Introduction

There has been great attention and focus on difference equations. Just like differential equations, there are
techniques that one can use to solve difference equations. One method for solving them is to use Lie symmetry
analysis. In this approach, one finds an invariant which can be used to find a simpler form of the equation.
Amongst the first people to use Lie symmetry analysis to solve difference equations areMaeda [12] and Hydon
[9]. Examples of the use of difference equations in real-life include modeling the proliferation of disease, loan
payments, population studies, etc [4,11].

Recurrence equations of a general order have been investigated in the literature and from different
approaches by several authors [1–3,5,6,8]. In their work, Almatrafi and others [1] discussed the exact solutions,
stability, oscillation and periodic aspects of the the difference equations

ηn+1 = ηn−4k+1

±1 ±∏k
i=1 ηn−4i+1

. (1)

These equations are special cases of a more generalized setting

ηn+1 = ηn−4k+1

an + bn
∏k

i=1 ηn−4i+1
, (2)

where an and bn are arbitrary real sequences. One can readily see that with a suitable change of variables, the
above equations can be derived from the higher order Berverton–Holt difference [3] equation

rn+l = μnKnrn
Kn + (μn − 1)rn

, (3)

where μn > 1 denotes the growth rate, Kn is the carrying capacity and r0, r1, . . . , rl−1 are the positive initial
values.

M. Folly-Gbetoula (B)
School of Mathematics, University of the Witwatersrand (Wits), Johannesburg 2050, South Africa
E-mail: Mensah.Folly-Gbetoula@wits.ac.za

123

http://crossmark.crossref.org/dialog/?doi=10.1007/s40065-023-00438-9&domain=pdf
http://orcid.org/0000-0002-3046-0679


542 Arab. J. Math. (2023) 12:541–551

In this paper, we perform a symmetry analysis of the equivalent equation

un+4k = un

An + Bn
∏k

i=1 un+4(i−1)
, (4)

for some arbitrary real sequences An and Bn where u0, u1, . . . , u4k−1 are initial values. Using symmetries, we
obtain explicit formulas for the solution of (4) and we deduce the solutions of (2) from those of (4). We also
study the periodic nature of the solutions and a stability analysis of the difference equation is investigated.

The derivation of symmetries for higher order recurrences involves cumbersome calculations and, to the
best of our knowledge, there are no computer packages that generate these symmetries. For more on Lie
symmetry analysis of difference equations, the reader can refer to [9] and, among others, the articles [6,7,13].

1.1 Preliminaries

Consider the difference equation:

un+4k = G(n, un, un+4, un+8 . . . , un+4k−4), (5)

for some smooth function G satisfying ∂G/∂un �= 0. Symmetry groups are connected to the determination of
infinitesimal transformations. Let

ûn = un + εQ(n, un) + O(ε2), (6)

be the one parameter Lie group of transformations of (5) with the corresponding generator

V = Q(n, un)
∂

∂un
. (7)

Note that the knowledge of the characteristic Q = Q(n, un) requires the knowledge of the (4k − 4)-th
prolongation of V:

V [4k−4] = Q
∂

∂un
+ (S4Q)

∂

∂un+4
+ · · · + (S4k−4Q)

∂

∂un+4k−4
, (8)

where Si : n → n + i is the forward shift operator. It is known that (6) is a symmetry group if and only if the
condition

S4k Q(n, un) − V [4k−4](G) = 0
∣
∣
∣
un+5k=G(n,un ,un+4,...,un+4k−4)

(9)

holds. Suppose that the characteristic is obtained by solving the functional equation (9). One can use the
canonical coordinate [10]

cn =
∫

dun
Q(n, un)

(10)

to derive the invariantswhichmay be used to lower the order of the difference equations. In [9], the author attests
that with the choice of canonical coordinate (10), the recurrence equation can, without fail, be represented in
the form cn+1 − cn = dn whose solution takes the shape

cn =
n∑

k=n0

dk + w1 (11)

for some constant w1. From (11), it is not difficult to find the solution expressed in terms of the original
variables. In this paper, our solution is obtained via the use of the canonical coordinate through a different
methodology.

The following theorem and definitions [8] are useful for studying local and globally stability aspects of the
equilibrium point.

123



Arab. J. Math. (2023) 12:541–551 543

Definition 1.1 The equilibrium point ū of (5) is said to be locally stable if for any ε > 0 such that if {un}∞n=0
is a solution of (5) with

|u0 − ū| + |u1 − ū| + · · · + |u4k−2 − ū| + |u4k−1 − ū| < δ, (12)

then

|un − ū| < ε for all n ≥ 0. (13)

Definition 1.2 The equilibrium point ū of (5) is said to be a global attractor if for any solution {un}∞n=0 of (5),

lim
n←∞ un = ū. (14)

Definition 1.3 The equilibrium point ū of (5) is globally asymptotically stable if ū is locally stable and is a
global attractor of (5).

Let

pi = ∂ f

∂un+i
(ū, . . . , ū), i = 4r, r = 0, 1, . . . , k − 1. (15)

It follows that

λ4k − p4k−4λ
4k−4 − · · · − p4λ

4 − p0 = 0 (16)

is the corresponding characteristic equation of (5) about the equilibrium point ū.

Theorem 1.4 Suppose f is a smooth function defined on some open neighborhood of ū. Then the following
statements are true:

(i) The equilibrium point ū is locally asymptotically stable if all the roots of (16) have absolute value less
than one.

(ii) The equilibrium point ū is unstable if at least one root of (16) has absolute value greater than one.

Definition 1.5 The equilibrium point ū of (5) is called non-hyperbolic if there exists a root of (16) with
absolute value equal to one.

2 Lie analysis and solutions

To derive the characteristic function Q admitted by (4), that is,

un+4k = G = un

An + Bn
∏k

i=1 un+4(i−1)
, (17)

we apply the symmetry constraint equation (9) to (17) to get

Q(n + 4k, un+4k) −
k−1∑

i=0

G,un+4i Q(n + 4i, un+4i ) = 0, (18)

where f,x denotes the partial derivative of f with respect to x . To solve for Q, we first apply the differential
operator ∂/∂un+4 − (un+8/un+4)∂/∂un+8 on (18). This gives
(
un+8

un+4
G,un+8un+4 − G,un+4un+4

)

Q(n + 4, un+4) +
(
un+8

un+4
G,un+8un+8 − G,un+4un+8

)

Q(n + 8, un+8)

+
∑

i≥3

(
un+8

un+4
G,un+8un+4i − G,un+4un+4i

)

Q(n + 4i, un+4i ) − Bnu2n
∏k

i=3 un+4(i−1)
(
An + Bn

∏k
i=1 un+4(i−1)

)2
(
Q′(n + 8, un+8)

− Q′(n + 4, un+4)
)+

(
un+8

un+4
G,un+8un − G,un+4un

)

Q(n, un) = 0 (19)

123



544 Arab. J. Math. (2023) 12:541–551

which simplifies to

un+8Q(n + 4, un+4) − un+4Q(n + 8, un+8) + un+4un+8(Q
′(n + 8, un+8) − Q′(n + 4, un+4)) = 0. (20)

We then differentiate the above equation with respect to un+4 twice to get

(
un+4Q

′(n + 4, un+4) − Q(n + 4, un+4)
)′′ = 0. (21)

The general solution of (21) takes the form

Q(n, un) = αnun + βn(un ln un + un) + γn (22)

for some functions αn , βn and γn of n. Next, we substitute Q and its corresponding shifts in (18). Bearing in
mind that αn , βn and γn are independent of un and their shifts, we use the method of separation. It turns out
that βn is equal to zero and the system of overdetermined equations resulting from the separation is as follows:

1 : A2
nγn+4k − Anγn = 0

un : Anαn+4k − Anαn = 0

u2nun+4 . . . un+4k−4 : Bn(αn+4k + αn+4 + αn+8 + · · · + αn+4k−8 + αn+4k−4) = 0

unun+4 . . . un+4k−4 : 2AnBnγn+4k = 0 (23)

which reduces to

γn = 0, αn+4k − αn = 0, (24)

αn + αn+4 + αn+8 + · · · + αn+4k−8 + αn+4k−4 = 0. (25)

Equation (25) is a linear difference equation with constant coefficients and has the characteristic equation

1 + r4 + r8 + · · · + r4k−8 + r4k−4 = 0. (26)

It is well known that if r is a solution of (26), then αn = rn is a solution of (25). Multiplying (26) by 1 − r4

and solving the resulting equation gives

r1(s) = e
isπ
2k , r2(s) = −e

isπ
2k , r3(s) = ie

isπ
2k , r4(s) = −ie

isπ
2k , (27)

for 1 ≤ s ≤ k − 1. Thus, from (22), the finite dimensional Lie algebra is spanned by the vectors fields

X1(s) = e
insπ
2k un

∂

∂un
, X2(s) = (−1)ne

insπ
2k un

∂

∂un
, (28)

X3(s) = ine
insπ
2k un

∂

∂un
, X4(s) = (−i)ne

insπ
2k un

∂

∂un
(29)

for 1 ≤ s ≤ k − 1. To obtain the compatible variable, we use the canonical coordinate

C(n) =
∫

dun
αnun

= 1

αn
ln |un|, (30)

where αn satisfies (25). Replacing αn with C(n)αn in the left hand side of equation in (25) yields the group
invariant In = ln |unun+4 . . . un+4k−4| since (Xr (s))In = 0 for r = 1, 2, 3, 4. Hence, using (17), we have
e−In+4 = Ane−In + Bn . For the sake of simplicity, we instead use the invariant

rn = exp (−In) = 1

unun+4 . . . un+4k−4
. (31)

On one hand, shifting (31) four times and replacing un+4k in the resulting equation yields

rn+4 = Anrn + An (32)

123



Arab. J. Math. (2023) 12:541–551 545

whose iteration gives

r4n+ j = r j

⎛

⎝
n−1∏

k1=0

A4k1+ j

⎞

⎠+
n−1∑

l=0

⎛

⎝B4l+ j

n−1∏

k2=l+1

A4k2+ j

⎞

⎠ (33)

for j = 0, 1, 2, 3. On the other hand, using the same relation given in (31), we have

un+4k = rn
rn+4

un . (34)

We iterate (34) and its solution in closed form takes the form

u4kn+i = ui

(
n−1∏

s=0

r4ks+i

r4ks+4+i

)

, i = 0, . . . , 4k − 1. (35)

From the known fact that any integer r can be written as r = 4	r/4
 + τ(r), 0 ≤ τ(r) ≤ 3, where τ(r) is the
remainder when r is divided by 4, we can rewrite (35) as follows:

u4kn+i = ui

(
n−1∏

s=0

r4(ks+	 i
4 
)+τ(i)

r4(ks+	 i+4
4 
)+τ(i+4)

)

, (36)

where i = 0, . . . , 4k − 1 and 0 ≤ τ(i) ≤ 3. Using (35) in (36), we obtain

u4kn+i =ui

n−1∏

s=0

rτ(i)

⎛

⎝
∏

ks+
	 i
4 
−1

k1=0 A4k1+τ(i)

⎞

⎠ +∑
ks+

	 i
4 
−1

l=0

⎛

⎝B4l+τ(i)
∏

ks
+	 i

4 
−1
k2=l+1 A4k2+τ(i)

⎞

⎠

rτ(i+4)

⎛

⎝
∏

ks+
	 i+4

4 
−1
k1=0 A4k1+τ(i+4)

⎞

⎠ +∑
ks+

	 i+4
4 
−1

l=0

⎛

⎝B4l+τ(i+4)
∏

ks+
	 i+4

4 
−1
k2=l+1 A4k2+τ(i+4)

⎞

⎠

, (37)

i = 0, . . . , 4k − 1. Noting that τ(i + 4) = τ(i) and 1/ri = ∏k−1
j=0 ui+4 j , the above equation simplifies

u4kn+i =ui

n−1∏

s=0

⎛

⎝
∏

ks+
	 i
4 
−1

k1=0 A4k1+τ(i)

⎞

⎠+
(∏k−1

j=0 uτ(i)+4 j

)∑
ks+

	 i
4 
−1

l=0

⎛

⎝B4l+τ(i)
∏

ks
+	 i

4 
−1
k2=l+1 A4k2+τ(i)

⎞

⎠

(
∏ks+	 i

4 

k1=0 A4k1+τ(i)

)

+
(∏k−1

j=0 uτ(i)+4 j

)∑ks+	 i
4 


l=0

(

B4l+τ(i)
∏ks+	 i

4 

k2=l+1 A4k2+τ(i)

) , (38)

i = 0, . . . , 4k − 1. The solution of (2) is obtained by back shifting (38) 4k − 1 times. Hence, the closed form
solution of (2) is given by

η4kn−4k+1+i =ηi−4k+1

n−1∏

s=0

⎛

⎝
∏

ks+
	 i
4 
−1

k1=0 a4k1+τ(i)

⎞

⎠+
(∏k−1

j=0 ητ(i)−4k+1+4 j

)∑
ks+

	 i
4 
−1

l=0

⎛

⎝b4l+τ(i)
∏

ks+
	 i
4 
−1

k2 = l+1a4k2 + τ(i)

⎞

⎠

(
∏ks+	 i

4 

k1=0 a4k1+τ(i)

)

+
(∏k−1

j=0 ητ(i) −4k+1+4 j

)∑ks+	 i
4 


l=0

(

b4l+τ(i)
∏ks+	 i

4 

k2=l+1a4k2+τ(i)

) .

(39)

Observe that if {an} and {bn} are constant sequences, i.e. an = a for all n and bn = b for all n, then the
solutions of (2) and (4) are given by

η4kn−4k+1+i =ηi−4k+1

n−1∏

s=0

aks+	 i
4 
 + b

(∏k−1
j=0 ητ(i)−4k+1+4 j

)∑
ks+

	 i
4 
−1

l=0 al

aks+	 i
4 
+1 + b

(∏k−1
j=0 ητ(i)−4k+1+4 j

)∑ks+	 i
4 


l=0 al
(40)

123



546 Arab. J. Math. (2023) 12:541–551

and

u4kn+i =ui

n−1∏

s=0

Aks+	 i
4 
 + B

(∏k−1
j=0 uτ(i)+4 j

)∑
ks+

	 i
4 
−1

l=0 Al

Aks+	 i
4 
+1 + B

(∏k−1
j=0 uτ(i)+4 j

)∑
ks+
	 i
4 


l=0 Al

, (41)

respectively.
In the following section, we investigate some special cases. One of the aims is to realize some results in

[1].

3 Special cases

3.1 The case when a = 1 and b is a constant

We investigate the case when a = 1 and b constant. In this case, from (40), the solution is given by

η4kn−4k+1+i =ηi−4k+1

n−1∏

s=0

1 + b
(∏k−1

j=0 ητ(i)−4k+1+4 j

)
(ks + 	 i

4
)
1 + b

(∏k−1
j=0 ητ(i)−4k+1+4 j

)
(ks + 	 i

4
 + 1)
, (42)

i = 0, . . . , 4k − 1. Replacing b = ±1 in (42) yields the results in [1, see Theorems 1 and 6]. In fact, noting
that τ(4k − 1 − j) = 3 − τ( j) = 3 − j + 4	 j/4
, we have:

η4kn− j = η4kn−4k+1+(4k−1− j), j = 0, 1, . . . , 4k − 1 (43)

= η− j

n−1∏

s=0

1 + b
(∏k−1

r=0 ητ(4k−1− j)−4k+1+4r

)
(ks + 	 4k−1− j

4 
)
1 + b

(∏k−1
r=0 ητ(4k−1− j)−4k+1+4r

)
(ks + 	 4k−1− j

4 
 + 1)
(44)

= η− j

n−1∏

s=0

1 + b
(∏k−1

r=0 η− j+4	 j
4 
−4r

)
(ks + k − 1 − 	 j

4 
)
1 + b

(∏k−1
r=0 η− j+4	 j

4 
−4r

)
(ks + k − 	 j

4 
)
. (45)

3.2 The case when a �= 1 and b is a constant

Here, from (40), the solution is given by

η4kn−4k+1+i =ηi−4k+1

n−1∏

s=0

aks+	 i
4 
 + b

(∏k−1
j=0 ητ(i)−4k+1+4 j

)(
1−aks+	 i4 


1−a

)

aks+	 i
4 
+1 + b

(∏k−1
j=0 ητ(i)−4k+1+4 j

)(
1−aks+	 i4 
+1

1−a

) (46)

and, similarly, this can also be written in the form

η4kn− j =η− j

n−1∏

s=0

aks+k−1−	 i
4 
 + b

(
∏k−1

j=0 η− j+4	 j
4 
−4r

)(
1−aks+k−1+	 i4 


1−a

)

aks+k−	 i
4 
 + b

(
∏k−1

j=0 η− j+4	 j
4 
−4r

)(
1−aks+k−	 i4 


1−a

) . (47)

For a = −1, the above equation simplifies to

123



Arab. J. Math. (2023) 12:541–551 547

• For k even,

η4kn− j = η− j

n−1∏

s=0

−(−1)	 i
4 
 + b

(∏k−1
j=0 η− j+4	 j

4 
−4r

)(
1+(−1)	

i
4 


2

)

(−1)	 i
4 
 + b

(∏k−1
j=0 η− j+4	 j

4 
−4r

)(
1−(−1)	

i
4 


2

) (48)

= η− j

⎡

⎣−1 + b

⎛

⎝
k−1∏

j=0

η− j+4	 j
4 
−4r

⎞

⎠

⎤

⎦

(−1)	
j
4 
n

. (49)

This result was obtained in [1] for b = ±1 (see Theorems 11 and 18).
• For k odd,

η4kn− j = η− j

n−1∏

s=0

(−1)s−	 i
4 
 + b

(∏k−1
j=0 η− j+4	 j

4 
−4r

)(
1−(−1)s−	 i4 


2

)

−(−1)s−	 i
4 
 + b

(∏k−1
j=0 η− j+4	 j

4 
−4r

)(
1+(−1)s−	 i4 


2

) (50)

=

⎧
⎪⎪⎨

⎪⎪⎩

η− j , if n is even

η− j

[

−1 + b

(
k−1∏

r=0
η− j+4	 j

4 
−4r

)](−1)	
j
4 
+1

, if n is odd
(51)

for j = 0, 1, . . . , 4k − 1 and furthermore,

η8kn− j = η− j (52)

for all n. More explicitly, the 8k periodic solutions are as follows:

ηi = η4k(1)−(4k−i) = η−4k+i[
−1 + b

(∏k−1
r=0 η−4(r+1)+i

)] , i = 1, 2, 3, 4 (53)

η4+i = η4k(1)−(4k−i−4) = η−4k+i+4

⎡

⎣−1 + b

⎛

⎝
k−1∏

r=0

η−4(r+1)+i

⎞

⎠

⎤

⎦, i = 1, 2, 3, 4 (54)

...
...

...
...

η4k−8+i = η4k(1)−(8−i) = η8−i

⎡

⎣−1 + b

⎛

⎝
k−1∏

r=0

η−4(r+1)+i

⎞

⎠

⎤

⎦, i = 1, 2, 3, 4 (55)

η4k−4+i = η4k(1)−(4k−i−4) = η−4k+i+4[
−1 + b

(∏k−1
r=0 η−4(r+1)+i

)] , i = 1, 2, 3, 4 (56)

η4k+1+i = η4k(2)−(4k−1−i) = η−4k+1+i , i = 0, 1, . . . , 4k − 1. (57)

For this special case, the results were obtained in [1] for b = ±1 (see Theorems 9, 15 and 16).

4 Periodicity and behavior of the solutions

Theorem 4.1 Let un be a solution of

un+4k = un

A + B
∏k

i=1 un+4(i−1)
, (58)

for some non-zero constants A �= 1 and B. Suppose the initial conditions xi , i = 0, . . . , 4k − 1, are such that
∏k−1

j=0 u4 j+p = (1 − A)/B, p = 0, 1, 2, 3. Then the solution of (58) is periodic with period 4k.

123



548 Arab. J. Math. (2023) 12:541–551

Fig. 1 Graph of xn+8 = xn
(2 − xnxn+4)

, where x0 = −2, x1 = −3, x2 = −4, x3 = 1, x4 = −1/2, x5 = −1/3, x6 = −1/4, x7 =
1 and are such that x0x4 = x1x5 = x2x6 = x3x7 = (1 − A)/B

Proof Suppose
∏k−1

j=0 u4 j+p = (1 − A)/B, p = 0, 1, 2, 3. From (41), we get

u4kn+i = ui

n−1∏

s=0

Aks+	 i
4 
 + B

(∏k−1
j=0 uτ(i)+4 j

)∑
ks+

	 i
4 
−1

l=0 Al

Aks+	 i
4 
+1 + B

(∏k−1
j=0 uτ(i)+4 j

)∑ks+	 i
4 


l=0 Al
(59)

= ui
∏n−1

s=0

Aks+	 i
4 
 + B

( 1−A
B

)
(

1−Aks+	 i4 

1−A

)

Aks+	 i
4 
+1 + B

( 1−A
B

)
(

1−Aks+	 i4 
+1

1−A

) (60)

= ui , (61)

for all i = 0, 1, . . . , 4k − 1 since 0 ≤ τ(i) ≤ 3. ��
We plot Fig. 1 to illustrate Theorem 4.1. We note that for A = −1 and B = 1, we get the result in

Theorem 13 in [1] and the result’s restriction (
∏k−1

r=0 η− j+4	 j
4 
−4r = 2, j = 0, 1, . . . , 4k − 1 or simply

− j + 4	 j
4 
 = 0, 1, 2, 3) is a special case of the assumption in the above theorem (

∏k−1
j=0 u4 j+p = (1− A)/B,

that is,
∏k−1

r=0 ηp−4r	 j
4 
−4r = (1 − A)/B, p = 0, 1, 2, 3).

Observe that in Theorem 15 in [1], the authors ought to add the restriction v j �= −2. If this condition is
not satisfied, the period will be 4k and not 8k as they clearly stated in Theorem 20 in [1].

Theorem 4.2 Let un be a solution of

un+4k = un

1 + B
∏k

i=1 un+4(i−1)
, (62)

for some non-zero constant B. The zero equilibrium point is non hyperbolic. Furthermore, if the initial condi-
tions xi , i = 0, . . . , 4k − 1, and B are positive, then the solution converges to the zero equilibrium point.

Proof The equilibrium point of (62) is u = 0. Let

f (un, un+4, . . . , un+4(k−1)) = un+4k = un

1 + B
∏k

i=1 un+4(i−1)
. (63)

123



Arab. J. Math. (2023) 12:541–551 549

Fig. 2 Graph of xn+8 = xn/(2 + xnxn+4), where x0 = 2, x1 = 3, x2 = 4, x3 = 1, x4 = 1/2, x5 = 1/3, x6 = 1/4, x7 = 1

So,

f,un = 1
(
1 + B

∏k
i=1 un+4(i−1)

)2 , f,un+4 j = −Bu2n
(
1 + B

∏k
i=1 un+4(i−1)

)2

∏k−1
i=1
i �= j

un+4i , j = 1, 2, . . . , k − 1. (64)

We have f,un (0, 0, . . . , 0) = 1 and f,un+4 j (0, 0, . . . , 0) = 0, j = 1, 2, . . . , k − 1. Thus, the characteristic
equation associated with (58) is λ4k − 1 = 0 and therefore, |λi | = 1. Therefore, the zero equilibrium point is
non-hyperbolic.

Suppose the non-zero initial conditions are all positive. From (41), we get

u4kn+i = ui

n−1∏

s=0

1 + B
(∏k−1

j=0 uτ(i)+4 j

) (
ks +

⌊
i
4

⌋)

1 + B
(∏k−1

j=0 uτ(i)+4 j

) (
ks +

⌊
i
4

⌋
+ 1

) (65)

= ui

n−1∏

s=0

⎛

⎝1 −
B
(∏k−1

j=0 uτ(i)+4 j

)

1 + B
(∏k−1

j=0 uτ(i)+4 j

) (
ks +

⌊
i
4

⌋
+ 1

)

⎞

⎠ (66)

= ui

n−1∏

s=0


(s). (67)

If B is positive, 
(s) < 1, s = 0, 1, . . . , n − 1. Therefore, un tends to zero as n tends to infinity. ��
We plot Fig. 2 to illustrate Theorem 4.2.

Theorem 4.3 Assume that B is positive and ui ≥ 0, i = 0, 1, . . . , k−1.Then the zero equilibrium point u = 0
of (62) is globally asymptotically stable.

Proof The equilibrium point of (62) satisfies u(1 + Buk) = 0. Thus, u = 0. Let ε ≥ 0 and suppose the
u′
i s, i = 0, 1, . . . , k − 1 are such that

|ui | ≤ ε

4k(B + 1)
, i = 0, 1, . . . , k − 1.

We have that

|u0| + |u1| + · · · + |uk−1| ≤ ε

B + 1
(68)

and (see (67)),

u4kn+i ≤ ui , (69)

123



550 Arab. J. Math. (2023) 12:541–551

i = 0, 1, . . . , k − 1, for all n if B ≥ 0.
This implies that for |u4kn+i | ≤ |ui | ≤ ε

4k(B+1) ≤ ε, we have found δ = ε/(B + 1) such that |u0|+ |u1|+
. . . |uk−1| ≤ δ. Thus, the zero equilibrium point is locally stable. On the other hand (see Theorem 4.2), xn
tends to zero as n goes to infinity. The zero equilibrium being a global attractor and locally stable, it is globally
asymptotically stable. ��
Theorem 4.4 Assume A �= 1. The zero equilibrium point of (58) is asymptotically stable for |A| > 1 and
unstable for |A| < 1. Furthermore, all non zero equilibrium points of (58) are non-hyperbolic.

Proof The equilibrium points of (58) satisfy u(A + Buk − 1) = 0. Let

f (un, un+4, . . . , un+4(k−1)) = un+4k = un

A + B
∏k

i=1 un+4(i−1)
. (70)

We have

f,un = A
(
A + B

∏k
i=1 un+4(i−1)

)2 , f,un+4 j = −Bu2n
(
A + B

∏k
i=1 un+4(i−1)

)2

k−1∏

i=1
i �= j

un+4i , j = 1, 2, . . . , k − 1. (71)

• For the equilibrium point u = 0, we have f,un (0, 0, . . . , 0) = 1/A and f,un+4 j (0, 0, . . . , 0) = 0, j =
1, 2, . . . , k − 1. Thus, the characteristic equation associated with (58) is λ4k − 1

A = 0. Therefore, |λ| < 1
if |A| > 1 (that is, locally asymptotically stable) and |λ| > 1 if |A| < 1 (that is, unstable).

• The non-zero equilibrium points u satisfy A + Buk − 1 = 0. Then f,un (u, u, . . . , u) = A and
f,un+4 j (u, u, . . . , u) = A − 1, j = 1, 2, . . . , k − 1. Thus, the characteristic equation associated with
(58) is

λ4k − (A − 1)λ4k−4 − · · · − (A − 1)λ4 − A = 0. (72)

Multiplying the above equation by 1 − λ4, we get (after simplification)

(1 − λ4k)(λ4 − A) = 0. (73)

It follows that, for A < 0, the solutions of (73) are λr = −Aei
i(2r+1)π

4 , r = 0, 1, 2, 3 or λr = ei
2rπ
4k , p =

1, 2, . . . , k−1, k+1, . . . , 2k−1, 2k+1, . . . , 2k−1, 2k+1, . . . , 3k−1, 3k+1, . . . , 4k−1. For A > 0,
the solutions of (73) are λr = Aei

2rπ
4 , r = 0, 1, 2, 3 or λr = ei

2rπ
4k , p = 1, 2, . . . , k − 1, k + 1, . . . , 2k −

1, 2k + 1, . . . , 2k − 1, 2k + 1, . . . , 3k − 1, 3k + 1, . . . , 4k − 1. Therefore, for k > 1, there exists a root
of (72) with modulus equal to one.

��

5 Conclusion

We studied the difference equation ηn+1 = ηn−4k+1/(an + bn
∏k

i=1 ηn−4i+1) by performing its symmetry
analysis and we used the canonical coordinate to obtain its invariants. These invariants are utilized to derive the
solutions in closed form.We demonstrated that all the formula solutions in [1] are special cases of our findings.
Some conditions for existence of 4k and 8k periodic solutions were established. Finally, we investigated the
stability of the solution of the difference equation and proved the existence of non-hyperbolic and globally
asymptotically stable equilibrium points.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use,
sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original
author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other
third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit
line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted
by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To
view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

123

http://creativecommons.org/licenses/by/4.0/


Arab. J. Math. (2023) 12:541–551 551

Funding This research was funded by the National Research Foundation (NRF) of South Africa, Grant number: 132108.

Data availability Not applicable.

Declarations

Conflict of interest The author declares that there is no conflict of interest regarding the publication of this paper.

References

1. Aljoufi, L.S.; Almatrafi, M.B.; Seadawy, A.R.: Dynamical analysis of discrete time equations with a generalized order. Alex.
Eng. J. 64, 937–945 (2023)

2. Banasiak, J.: Mathematical Modelling in One Dimension: An Introduction via Difference and Differential Equation. Cam-
bridge University Press, Cambridge (2013)

3. Bohner, M.; Dannan, F.M.; Streipert, S.: A nonautonomous Beverton–Holt equation of higher order. J. Math. Anal. Appl.
457, 114–133 (2018)

4. Chen, S.; Jiang, X.: Modeling repayment behavior of consumer loan in portfolio across business cycle: a triplet Markov
model approach. Complexity 2020, 5458941, 11 (2020)

5. Elsayed, E.M.; Ibrahim,T.F.: Periodicity and solutions for some systems of non-linear rational difference equations.Hacettepe
J. Math. Stat. 44(6), 1361–1390 (2015)

6. Folly-Gbetoula, M.; Kgatliso Mkhwanazi, K.; Nyirenda, D.: On a study of a family of higher order recurrence relations.
Math. Probl. Eng. 2022, 6770105, 11 (2022)

7. Folly-Gbetoula, M.; Nyirenda, D.: A generalised two-dimensional system of higher order recursive sequences. J. Differ. Equ.
Appl. 26(2), 244–260 (2020)

8. Grove, E.A.; Ladas, G.: Periodicities in Nonlinear Difference Equations, vol. 4. Chapman & Hall/CRC, Boca Raton (2005)
9. Hydon, P.E.: Difference Equations by Differential Equation Methods. Cambridge University Press, Cambridge (2014)

10. Joshi, N.; Vassiliou, P.: The existence of Lie symmetries for first-order analytic discrete dynamical systems. J. Math. Anal.
Appl. 195, 872–887 (1995)

11. Li, Y.; Li, J.: Stage-structured discrete-time models for interacting wild and sterile mosquitoes with Beverton–Holt surviv-
ability. Math. Biosci. Eng. 16(2), 572–602 (2019)

12. Maeda, S.: The similarity method for difference equations. IMA J. Appl. Math. 38, 129–134 (1987)
13. Quispel, G.R.W.; Sahadevan, R.: Lie symmetries and the integration of difference equations. Phys. Lett. A 184, 64–70 (1993)

Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional
affiliations.

123


	On a family of higher order recurrence relations: symmetries, formula solutions, periodicity and stability analysis
	Abstract
	1 Introduction
	1.1 Preliminaries

	2 Lie analysis and solutions
	3 Special cases
	3.1 The case when a = 1 and b is a constant
	3.2 The case when a neq1 and b is a constant

	4 Periodicity and behavior of the solutions
	5 Conclusion
	References