1 Introduction

Here we use the following standard notation: \({\mathbb {N}}\), \({\mathbb {Z}}\), \({\mathbb {R}}\), and \({\mathbb {C}}\) are the sets of natural, integer, real, and complex numbers, respectively, and \({\mathbb {N}}_{0}={\mathbb {N}}\cup \{0\}\). If \(p,q\in {\mathbb {Z}}\), \(p\le q\), then we use the notation \(r=\overline {p,q}\) instead of writing \(p\le r\le q\), \(r\in {\mathbb {Z}}\).

In [10, p. 100], it is explained how the method of iteration can be used for approximating the values of a given function. Namely, if f is a continuous function, then it is supposed to compute

$$\begin{aligned} y=f(x) \end{aligned}$$
(1)

for a given value of variable x.

To compute the value of \(f(x)\), it is suggested that relation (1) is written in an implicit form \(F(x,y)=0\). An obvious way of writing relation (1) in an implicit form is by choosing \(F(x,y)=y-f(x)\) (or \(F(x,y)=f(x)-y\); for solving the equation, both choices are the same). However, depending on some properties of the function f, some other functions \(F(x,y)\) can be also chosen. For example, if f is invertible, then we can choose \(F(x,y)=x-f^{-1}(y)\). If the range of f does not contain the zero value, then, for example, we can also take \(F(x,y)=\frac{1}{f(x)}-\frac{1}{y}\), and so on. So, \(F(x,y)\) can be chosen in various ways, and the choice depends on its usefulness in solving the problem of computation of \(f(x)\). Two concrete examples of some suitable choices of function \(F(x,y)\) are given below.

Such a chosen function \(F(x,y)\) is a function of two variables, and the given x in relation (1) (a fixed number therein) belongs to the domain of definition of the function.

Assume that the functions \(F(x,y)\) and \(F_{y}'(x,y)\) are continuous and \(F_{y}'(x,y)\ne 0\). Let \(y_{n}\) be an approximation of y. Then by the Langrage mean-value theorem we have

$$ F(x,y_{n})=F(x,y_{n})-F(x,y)=(y_{n}-y)F_{y}'(x, \widehat{y}_{n}), $$

where \(\widehat{y}_{n}\) is a point between \(y_{n}\) and y, and, consequently,

$$\begin{aligned} y=y_{n}-\frac{F(x,y_{n})}{F_{y}'(x,\widehat{y}_{n})}. \end{aligned}$$
(2)

It is natural to assume that \(y_{n}\approx \widehat{y}_{n}\), so for computing \(f(x)\), from (2) we obtain the following recursive relation:

$$\begin{aligned} y_{n+1}=y_{n}-\frac{F(x,y_{n})}{F_{y}'(x,y_{n})},\quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(3)

This is the Newton method, applied to F as a function of y (here x is fixed).

Three nice examples for application of this method are given in [10]. Two of them are interesting for the present investigation. Hence we describe some relevant details.

The first example refers to computing reciprocals, that is, in this case, we have \(f(x)=\frac{1}{x}\). Hence, for a given x, relation (1) becomes \(y=\frac{1}{x}\).

To calculate it, choose the function \(F(x,y)=x-\frac{1}{y}\). It is supposed to find y for a given x such that \(F(x,y)=0\). Since \(F_{y}'(x,y)=\frac{1}{y^{2}}\), after some simple calculations, the recursive relation (3) becomes

$$\begin{aligned} y_{n+1}=y_{n}(2-xy_{n}), \quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(4)

The recursive relation (4) appears in many problem books. For example, it is Problem 639.1 in the well-known problem book [9]. The reader might have seen it in the literature, but probably many are not aware of the above procedure, which leads to getting the relation. The problem in [9] is the following:

Problem 1

Let \(x>0\), and let \((y_{n})_{n\in {\mathbb {N}}_{0}}\) satisfy (4). Show that if \(\min \{y_{0},y_{1}\}>0\), then \(y_{n}\) converges, and

$$ \lim_{n\to +\infty }y_{n}=\frac{1}{x}. $$

The problem is relatively simple and can be solved in several elementary ways. Some nice analyses related to the convergence of the sequence \(y_{n}\) are given in [10]. It is interesting that (4) is a difference equation solvable in closed form.

The hint given in [9] suggests consideration of the sequence \(\frac{1}{x}-y_{n}\). From (4) we have

$$ \frac{1}{x}-y_{n+1}=\frac{1}{x}-2y_{n}+xy_{n}^{2},\quad n\in {\mathbb {N}}_{0}, $$

and, consequently,

$$\begin{aligned} 1-xy_{n+1}=(1-xy_{n})^{2},\quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(5)

From (5) we easily obtain

$$ 1-xy_{n}=(1-xy_{0})^{2^{n}},\quad n\in {\mathbb {N}}_{0}, $$

from which it follows that

$$\begin{aligned} y_{n}=\frac{1-(1-xy_{0})^{2^{n}}}{x}, \quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(6)

Formula (6) is a closed-form formula for general solution to equation (4) in the case \(x\ne 0\). The case \(x=0\) is very simple, since in this case, (4) defines a geometric progression with quotient equal to two, and hence \(y_{n}=2^{n}y_{0}\), \(n\in {\mathbb {N}}_{0}\). However, from the practical point of view, this case is not interesting, since for \(x=0\), the function \(\frac{1}{x}\) is not defined. By using formula (6) the long-term behavior of solutions to equation (6) can be easily described.

The second example in [10] refers to computing the square roots of positive numbers. This means that \(f(x)=\sqrt{x}\) in (1). Further, the following function is chosen \(F(x,y)=y^{2}-x\). Hence, in this case, equation (3) becomes

$$\begin{aligned} y_{n+1}=y_{n}-\frac{y_{n}^{2}-x}{2y_{n}}=\frac{1}{2} \biggl(y_{n}+ \frac{x}{y_{n}} \biggr),\quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(7)

It is interesting that equation (7) is also solvable. Namely, we have

$$ y_{n+1}\pm \sqrt{x}=\frac{(y_{n}\pm \sqrt{x})^{2}}{2y_{n}},\quad n\in {\mathbb {N}}_{0}, $$

from which it follows that

$$ \frac{y_{n+1}+\sqrt{x}}{y_{n+1}-\sqrt{x}}= \biggl( \frac{y_{n}+\sqrt{x}}{y_{n}-\sqrt{x}} \biggr)^{2},\quad n \in {\mathbb {N}}_{0}, $$

and, consequently,

$$\begin{aligned} \frac{y_{n}+\sqrt{x}}{y_{n}-\sqrt{x}}= \biggl( \frac{y_{0}+\sqrt{x}}{y_{0}-\sqrt{x}} \biggr)^{2^{n}},\quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(8)

From (8) we obtain

$$ y_{n}=\sqrt{x} \frac{ (\frac{y_{0}+\sqrt{x}}{y_{0}-\sqrt{x}} )^{2^{n}}+1}{ (\frac{y_{0}+\sqrt{x}}{y_{0}-\sqrt{x}} )^{2^{n}}-1}, \quad n\in {\mathbb {N}}_{0}. $$

Now we know that equation (7) is a particular case of a wider class of solvable nonlinear difference equations (see, e.g., [39] and the related references therein).

Solvable difference equations have attracted attention of researchers in the eighteenth and nineteenth centuries. Classical books on difference equations and related topics such as [4, 5, 7, 8, 11, 13–17, 21–25] contain some information on the topic (see also the original sources [3, 6, 18–20]). In the last two decades, there has been some renewed interest in the topic (the reader can consult, e.g., [2, 26–42] and the related references therein).

Bearing in mind the above-mentioned facts related to equations (4) and (7), as well as some of our recent investigations on solvability of difference equations and systems, the following natural question arises.

Question 1

Is equation (4) also a particular case of a class of higher-order nonlinear difference equations solvable in closed form?

Here we give a positive answer to the question. We show that there is a class of theoretically solvable nonlinear difference equations that includes equation (4). We also show that some of them are practically solvable, presenting closed-form formulas for their general solutions.

Recall that we regard a difference equation as theoretically solvable if we know a form of its general solution, but some of the quantities appearing in the form of the general solution cannot be found in some suitable forms. For example, the homogeneous linear difference equations with constant coefficients

$$\begin{aligned} y_{n+k}-a_{k-1}y_{n+k-1}-\cdots -a_{1}y_{n+1}-a_{0}y_{n}=0,\quad n \in {\mathbb {N}}_{0}, \end{aligned}$$
(9)

where \(k\in {\mathbb {N}}\), \(a_{j}\in {\mathbb {C}}\), \(j=\overline {0,k-1}\), are theoretically solvable, but some of them are not practically solvable, since the roots of the associated characteristic polynomials cannot be always found by radicals [1]. The roots can be always found by radicals if the degree of the polynomials is less than or equal to four [12].

2 Main results

In this section, we present our main results. First, we quote the following auxiliary result, which is employed for several times in the proofs of some of the main results [38].

Lemma 1

Let \(m\in {\mathbb {N}}\), \(l\in {\mathbb {Z}}\), \(a_{j}\in {\mathbb {R}}\), \(j=\overline {1,m-1}\), \(a_{0}\in {\mathbb {R}}\setminus \{0\}\),

$$ R_{m}(s)=s^{m}-a_{m-1}s^{m-1}-\cdots -a_{1}s-a_{0}, $$

\(R_{m}(s_{k})=0\), \(k=\overline {1,m}\), \(s_{k}\ne s_{j}\), \(k\ne j\), and let \((z_{n})_{n\ge l-m}\) be the solution to the difference equation

$$ z_{n}=a_{m-1}z_{n-1}+\cdots +a_{1}z_{n-m+1}+a_{0}z_{n-m} $$

for \(n\ge l\) such that \(z_{j-m}=0\), \(j=\overline {l,l+m-2}\), \(z_{l-1}=1\). Then

$$ z_{n}=\sum_{k=1}^{m} \frac{s_{k}^{n+m-l}}{R_{m}'(s_{k})},\quad n\ge l-m. $$

Here we consider the higher-order nonlinear difference equation

$$\begin{aligned} y_{n+k}=y_{n+l}+y_{n}+ay_{n+l}y_{n},\quad n\in {\mathbb {N}}_{0}, \end{aligned}$$
(10)

where \(k\in {\mathbb {N}}\), \(l\in {\mathbb {N}}_{0}\), \(l< k\), \(a\in {\mathbb {C}}\), and \(y_{j}\in {\mathbb {C}}\), \(j=\overline {0,k-1}\).

First note that there are two cases to be considered. The case when \(a=0\), and the case when \(a\ne 0\).

Case \(a=0\). For \(a=0\), equation (10) becomes the following particular case of equation (9):

$$\begin{aligned} y_{n+k}=y_{n+l}+y_{n},\quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(11)

The equation is well known and appears from time to time in the literature. For example, it has appeared recently in [35].

The characteristic polynomial associated with equation (11) is

$$\begin{aligned} P_{k}(s)=s^{k}-s^{l}-1, \end{aligned}$$
(12)

and it can be certainly solved by radicals when \(k\le 4\). From this by a well-known theorem on the form of general solutions to homogeneous linear difference equations with constant coefficients [4, 13, 14, 21, 22, 24, 25] it follows that equation (11) is practically solvable in this case. By the Abel–Ruffini theorem [1] the roots of a polynomial of degree strictly greater than four cannot be always found by radicals. Hence some of the polynomials in (12) for \(k>5\) could be of this type.

Case \(a\ne 0\). Since \(a\ne 0\), we can multiply both sides of equation (10) by a. Then adding the unity to both sides of the obtained equation, we obtain the relation

$$ 1+ay_{n+k}=1+ay_{n+l}+ay_{n}+a^{2}y_{n+l}y_{n},\quad n\in {\mathbb {N}}_{0}. $$

Now note that the last relation can be written as

$$\begin{aligned} 1+ay_{n+k}=(1+ay_{n+l}) (1+ay_{n}), \quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(13)

Relation (13) is very important for solvability of equation (10). Namely, by the substitution

$$\begin{aligned} u_{n}=1+ay_{n}, \quad n\in {\mathbb {N}}_{0}, \end{aligned}$$
(14)

equation (13) becomes

$$\begin{aligned} u_{n+k}=u_{n+l}u_{n},\quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(15)

If the initial values \(u_{j}\), \(j=\overline {0,k-1}\), are positive numbers, then a simple inductive argument shows that all the members of such a solution are also positive. Hence, in this case, it is possible to take the logarithm of both sides of equation (15) and use the substitution \(v_{n}=\ln u_{n}\), \(n\in {\mathbb {N}}_{0}\), by which the equation is transformed into equation (11). Hence equation (15) is theoretically solvable in the case of positive initial values.

If some of the initial values are negative or complex numbers, then this method is not suitable for solving equation (15). However, there are some other methods for solving equations of this type, which are usually called product-type difference equations [33, 34].

If \(k\le 4\), then we can find the roots of the polynomial (12) by radicals, from which, together with the well-known theorem on the form of general solution to equation (9), we can find closed-form formulas for solutions to equation (11) and, consequently, closed-form formulas for solutions to equation (15). Such formulas can be used in the following consequence of relation (14):

$$\begin{aligned} y_{n}=\frac{u_{n}-1}{a},\quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(16)

Hence the following particular cases of equation (15) are certainly practically solvable:

$$\begin{aligned}& u_{n+2}=u_{n+1}u_{n},\quad n\in {\mathbb {N}}_{0}, \end{aligned}$$
(17)
$$\begin{aligned}& u_{n+3}=u_{n+1}u_{n},\quad n\in {\mathbb {N}}_{0}, \end{aligned}$$
(18)
$$\begin{aligned}& u_{n+3}=u_{n+2}u_{n},\quad n\in {\mathbb {N}}_{0}, \end{aligned}$$
(19)
$$\begin{aligned}& u_{n+4}=u_{n+1}u_{n},\quad n\in {\mathbb {N}}_{0}, \end{aligned}$$
(20)
$$\begin{aligned}& u_{n+4}=u_{n+2}u_{n},\quad n\in {\mathbb {N}}_{0}, \end{aligned}$$
(21)
$$\begin{aligned}& u_{n+4}=u_{n+3}u_{n},\quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(22)

From this it follows that to solve equation (10) in the cases \(k\le 4\), we should first solve equations (17)–(22) and then employ the obtained formulas for \(u_{n}\) in relation (16). To solve the equations, we use some methods for solving product-type difference equations (see, e.g., [30, 33, 34]).

2.1 Case \(k=2\), \(l=1\)

Let \(\alpha _{1}=1\) and \(\beta _{1}=1\). Then

$$ u_{n}=u_{n-1}^{\alpha _{1}}u_{n-2}^{\beta _{1}},\quad n\ge 2, $$

from which we obtain

$$ u_{n}=(u_{n-2}u_{n-3})^{\alpha _{1}}u_{n-2}^{\beta _{1}}=u_{n-2}^{ \alpha _{1}+\beta _{1}}u_{n-3}^{\alpha _{1}}=u_{n-2}^{\alpha _{2}}u_{n-3}^{ \beta _{2}},\quad n\ge 3, $$

where \(\alpha _{2}=\alpha _{1}+\beta _{1}\) and \(\beta _{2}=\alpha _{1}\).

By the same argument and induction we obtain that for each \(m\in {\mathbb {N}}\setminus \{1\}\),

$$\begin{aligned}& u_{n}=u_{n-m}^{\alpha _{m}}u_{n-m-1}^{\beta _{m}},\quad n\ge m+1, \end{aligned}$$
(23)
$$\begin{aligned}& \alpha _{m}=\alpha _{m-1}+\beta _{m-1},\quad \beta _{m}=\alpha _{m-1}. \end{aligned}$$
(24)

From (23) and (24) it easily follows that

$$\begin{aligned}& u_{n}=u_{1}^{\alpha _{n-1}}u_{0}^{\alpha _{n-2}}, \\& \alpha _{n}=\alpha _{n-1}+\alpha _{n-2}, \end{aligned}$$
(25)

and \(\alpha _{1}=1\) and \(\alpha _{2}=2\). Hence \(\alpha _{n}=f_{n+1}\) (\(f_{n}\) is the Fibonacci sequence [43]), and by (25) we obtain

$$\begin{aligned} u_{n}=u_{1}^{f_{n}}u_{0}^{f_{n-1}},\quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(26)

We obtain our first theorem by combining formulas (16) and (26).

Theorem 1

Let \(k=2\), \(l=1\), and \(a\in {\mathbb {C}}\setminus \{0\}\). Then the formula

$$ y_{n}=\frac{(1+ay_{1})^{f_{n}}(1+ay_{0})^{f_{n-1}}-1}{a},\quad n\in {\mathbb {N}}_{0}, $$

presents a general solution to equation (10) in this case.

2.2 Case \(k=3\), \(l=1\)

Let \(\alpha _{1}=1\), \(\beta _{1}=1\), and \(\gamma _{1}=0\). Then

$$\begin{aligned} u_{n}=u_{n-2}^{\alpha _{1}}u_{n-3}^{\beta _{1}}u_{n-4}^{\gamma _{1}},\quad n\ge 3, \end{aligned}$$
(27)

from which we obtain

$$ u_{n}=(u_{n-4}u_{n-5})^{\alpha _{1}}u_{n-3}^{\beta _{1}}u_{n-4}^{ \gamma _{1}}=u_{n-3}^{\beta _{1}}u_{n-4}^{\alpha _{1}+\gamma _{1}}u_{n-5}^{ \alpha _{1}} =u_{n-3}^{\alpha _{2}}u_{n-4}^{\beta _{2}}u_{n-5}^{ \gamma _{2}},\quad n\ge 5, $$

where \(\alpha _{2}=\beta _{1}\), \(\beta _{2}=\alpha _{1}+\gamma _{1}\), and \(\gamma _{2}=\alpha _{1}\).

By the same argument and induction we obtain that for each \(m\in {\mathbb {N}}\setminus \{1\}\),

$$\begin{aligned}& u_{n}=u_{n-m-1}^{\alpha _{m}}u_{n-m-2}^{\beta _{m}}u_{n-m-3}^{ \gamma _{m}},\quad n\ge m+3, \end{aligned}$$
(28)
$$\begin{aligned}& \alpha _{m}=\beta _{m-1},\qquad \beta _{m}=\alpha _{m-1}+\gamma _{m-1},\qquad \gamma _{m}=\alpha _{m-1}. \end{aligned}$$
(29)

From (29) it follows that

$$\begin{aligned} \alpha _{n}=\alpha _{n-2}+\alpha _{n-3}, \end{aligned}$$
(30)

and \(\alpha _{0}=0\), \(\alpha _{-1}=1\), \(\alpha _{-2}=\alpha _{-3}=0\), \(\alpha _{-4}=1\) ((29) can be also used for \(m\le 1\)).

Using (28)–(30), we get

$$\begin{aligned} u_{n}=u_{2}^{\alpha _{n-3}}u_{1}^{\beta _{n-3}}u_{0}^{\gamma _{n-3}}=u_{2}^{ \alpha _{n-3}}u_{1}^{\alpha _{n-2}}u_{0}^{\alpha _{n-4}},\quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(31)

We easily see that the roots \(s_{j}\), \(j=\overline {1,3}\), of the polynomial

$$\begin{aligned} Q_{3}(s)=s^{3}-s-1=0, \end{aligned}$$
(32)

are distinct. From this and Lemma 1 it follows that

$$\begin{aligned} \alpha _{n}=\sum_{j=1}^{3} \frac{s_{j}^{n+3}}{Q_{3}'(s_{j})}, \quad n \in {\mathbb {Z}}, \end{aligned}$$
(33)

is the solution to (30) with \(\alpha _{-3}=\alpha _{-2}=0\), and \(\alpha _{-1}=1\).

From (16) and (31) we obtain our second theorem.

Theorem 2

Let \(k=3\), \(l=1\), and \(a\in {\mathbb {C}}\setminus \{0\}\). Then the formula

$$\begin{aligned} y_{n}&= \frac{(1+ay_{2})^{\alpha _{n-3}}(1+ay_{1})^{\alpha _{n-2}}(1+ay_{0})^{\alpha _{n-4}}-1}{a},\quad n\in {\mathbb {N}}_{0}, \end{aligned}$$

where \(\alpha _{n}\) is defined in (33), presents a general solution to (10) in this case.

2.3 Case \(k=3\), \(l=2\)

Let \(\alpha _{1}=1\), \(\beta _{1}=0\), and \(\gamma _{1}=1\). Then

$$\begin{aligned} u_{n}=u_{n-1}^{\alpha _{1}}u_{n-2}^{\beta _{1}}u_{n-3}^{\gamma _{1}},\quad n\ge 3, \end{aligned}$$
(34)

from which we obtain

$$ u_{n}=(u_{n-2}u_{n-4})^{\alpha _{1}}u_{n-2}^{\beta _{1}}u_{n-3}^{ \gamma _{1}}=u_{n-2}^{\alpha _{1}+\beta _{1}}u_{n-3}^{\gamma _{1}}u_{n-4}^{ \alpha _{1}} =u_{n-2}^{\alpha _{2}}u_{n-3}^{\beta _{2}}u_{n-4}^{ \gamma _{2}},\quad n\ge 4, $$

where \(\alpha _{2}=\alpha _{1}+\beta _{1}\), \(\beta _{2}=\gamma _{1}\), and \(\gamma _{2}=\alpha _{1}\).

By the same argument and induction we obtain that for each \(m\in {\mathbb {N}}\setminus \{1\}\)

$$\begin{aligned}& u_{n}=u_{n-m}^{\alpha _{m}}u_{n-m-1}^{\beta _{m}}u_{n-m-2}^{\gamma _{m}},\quad n\ge m+2, \end{aligned}$$
(35)
$$\begin{aligned}& \alpha _{m}=\alpha _{m-1}+\beta _{m-1},\qquad \beta _{m}=\gamma _{m-1},\qquad \gamma _{m}=\alpha _{m-1}. \end{aligned}$$
(36)

From (36) we have

$$\begin{aligned} \alpha _{n}=\alpha _{n-1}+\alpha _{n-3}, \end{aligned}$$
(37)

and \(\alpha _{0}=1\), \(\alpha _{-1}=\alpha _{-2}=0\), \(\alpha _{-3}=1\), \(\alpha _{-4}=0\) ((36) can be also used for \(m\le 1\)).

Formulas (35) and (36) yield

$$\begin{aligned} u_{n}=u_{2}^{\alpha _{n-2}}u_{1}^{\beta _{n-2}}u_{0}^{\gamma _{n-2}}=u_{2}^{ \alpha _{n-2}}u_{1}^{\alpha _{n-4}}u_{0}^{\alpha _{n-3}},\quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(38)

It is not difficult to see that the roots \(s_{j}\), \(j=\overline {1,3}\), of the polynomial

$$\begin{aligned} R_{3}(s)=s^{3}-s^{2}-1=0 \end{aligned}$$
(39)

are different, from which, together with Lemma 1, it follows that

$$\begin{aligned} \alpha _{n}=\sum_{j=1}^{3} \frac{s_{j}^{n+2}}{R_{3}'(s_{j})},\quad n \in {\mathbb {Z}}, \end{aligned}$$
(40)

is the solution to (37) with \(\alpha _{-2}=\alpha _{-1}=0\), and \(\alpha _{0}=1\).

Using (38) in (16) we get the following theorem.

Theorem 3

Let \(k=3\), \(l=2\), and \(a\in {\mathbb {C}}\setminus \{0\}\). Then the formula

$$\begin{aligned} y_{n}&= \frac{(1+ay_{2})^{\alpha _{n-2}}(1+ay_{1})^{\alpha _{n-4}}(1+ay_{0})^{\alpha _{n-3}}-1}{a},\quad n\in {\mathbb {N}}, \end{aligned}$$

where \(\alpha _{n}\) is defined in (40), presents a general solution to (10) in this case.

2.4 Case \(k=4\), \(l=1\)

Let \(\alpha _{1}=\beta _{1}=1\) and \(\gamma _{1}=\delta _{1}=0\). Then

$$\begin{aligned} u_{n}=u_{n-3}^{\alpha _{1}}u_{n-4}^{\beta _{1}}u_{n-5}^{\gamma _{1}}u_{n-6}^{ \delta _{1}},\quad n\ge 4, \end{aligned}$$
(41)

from which we obtain

$$ u_{n}=(u_{n-6}u_{n-7})^{\alpha _{1}}u_{n-4}^{\beta _{1}}u_{n-5}^{ \gamma _{1}}u_{n-6}^{\delta _{1}}=u_{n-4}^{\beta _{1}}u_{n-5}^{ \gamma _{1}}u_{n-6}^{\alpha _{1}+\delta _{1}}u_{n-7}^{\alpha _{1}} =u_{n-4}^{ \alpha _{2}}u_{n-5}^{\beta _{2}}u_{n-6}^{\gamma _{2}}u_{n-7}^{ \delta _{2}} $$

for \(n\ge 7\), where \(\alpha _{2}=\beta _{1}\), \(\beta _{2}=\gamma _{1}\), \(\gamma _{2}=\alpha _{1}+\delta _{1}\), and \(\delta _{2}=\alpha _{1}\).

By the same argument and induction we obtain that for each \(m\in {\mathbb {N}}\setminus \{1\}\),

$$\begin{aligned}& u_{n}=u_{n-m-2}^{\alpha _{m}}u_{n-m-3}^{\beta _{m}}u_{n-m-4}^{ \gamma _{m}}u_{n-m-5}^{\delta _{m}},\quad n\ge m+5, \end{aligned}$$
(42)
$$\begin{aligned}& \alpha _{m}=\beta _{m-1},\qquad \beta _{m}=\gamma _{m-1}, \qquad \gamma _{m}=\alpha _{m-1}+\delta _{m-1},\qquad \delta _{m}=\alpha _{m-1}. \end{aligned}$$
(43)

From (43) we obtain

$$\begin{aligned} \alpha _{n}=\alpha _{n-3}+\alpha _{n-4}, \end{aligned}$$
(44)

and \(\alpha _{0}=\alpha _{-1}=0\), \(\alpha _{-2}=1\), \(\alpha _{-3}=\alpha _{-4}=\alpha _{-5}=0\), \(\alpha _{-6}=1\).

Using (42)–(44), we get

$$\begin{aligned} u_{n}=u_{3}^{\alpha _{n-5}}u_{2}^{\beta _{n-5}}u_{1}^{\gamma _{n-5}}u_{0}^{ \delta _{n-5}} =u_{3}^{\alpha _{n-5}}u_{2}^{\alpha _{n-4}}u_{1}^{ \alpha _{n-3}}u_{0}^{\alpha _{n-6}},\quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(45)

It is not difficult to see that the roots \(s_{j}\), \(j=\overline {1,4}\), of the polynomial

$$\begin{aligned} R_{4}(s)=s^{4}-s-1=0 \end{aligned}$$
(46)

are simple. From this by Lemma 1 it follows that

$$\begin{aligned} \alpha _{n}=\sum_{j=1}^{4} \frac{s_{j}^{n+5}}{R_{4}'(s_{j})},\quad n \in {\mathbb {Z}}. \end{aligned}$$
(47)

is the solution to (44) such that \(\alpha _{-5}=\alpha _{-4}=\alpha _{-3}=0\), \(\alpha _{-2}=1\).

The following result follows from (16) and (45).

Theorem 4

Let \(k=4\), \(l=1\), and \(a\in {\mathbb {C}}\setminus \{0\}\). Then the formula

$$\begin{aligned} y_{n}&= \frac{(1+ay_{3})^{\alpha _{n-5}}(1+ay_{2})^{\alpha _{n-4}}(1+ay_{1})^{\alpha _{n-3}}(1+ay_{0})^{\alpha _{n-6}}-1}{a},\quad n\in {\mathbb {N}}_{0}, \end{aligned}$$

where \(\alpha _{n}\) is defined in (47), presents a general solution to (10) in this case.

2.5 Case \(k=4\), \(l=2\)

In this case, we get (21), which is an equation with interlacing indices of second order [40, 42]. Thus the subsequences \(u_{2n}\) and \(u_{2n+1}\) are two solutions to equation (17). Employing Theorem 1, we obtain the following result.

Theorem 5

Let \(k=4\), \(l=2\), and \(a\in {\mathbb {C}}\setminus \{0\}\). Then the formulas

$$\begin{aligned} &y_{2n}=\frac{(1+ay_{2})^{f_{n}}(1+ay_{0})^{f_{n-1}}-1}{a}, \quad n\in {\mathbb {N}}_{0}, \\ &y_{2n+1}=\frac{(1+ay_{3})^{f_{n}}(1+ay_{1})^{f_{n-1}}-1}{a},\quad n \in {\mathbb {N}}_{0}, \end{aligned}$$

present a general solution to (10) in this case.

2.6 Case \(k=4\), \(l=3\)

Let \(\alpha _{1}=1\), \(\beta _{1}=0\), \(\gamma _{1}=0\), and \(\delta _{1}=1 \). Then

$$\begin{aligned} u_{n}=u_{n-1}^{\alpha _{1}}u_{n-2}^{\beta _{1}}u_{n-3}^{\gamma _{1}}u_{n-4}^{ \delta _{1}},\quad n\ge 4, \end{aligned}$$
(48)

from which we obtain

$$ u_{n}=(u_{n-2}u_{n-5})^{\alpha _{1}}u_{n-2}^{\beta _{1}}u_{n-3}^{ \gamma _{1}}u_{n-4}^{\delta _{1}}=u_{n-2}^{\alpha _{1}+\beta _{1}}u_{n-3}^{ \gamma _{1}}u_{n-4}^{\delta _{1}}u_{n-5}^{\alpha _{1}} =u_{n-2}^{ \alpha _{2}}u_{n-3}^{\beta _{2}}u_{n-4}^{\gamma _{2}}u_{n-5}^{ \delta _{2}} $$

for \(n\ge 5\), where \(\alpha _{2}=\alpha _{1}+\beta _{1}\), \(\beta _{2}=\gamma _{1}\), \(\gamma _{2}=\delta _{1}\), and \(\delta _{2}=\alpha _{1}\).

By the same argument and induction we obtain that for each \(m\in {\mathbb {N}}\setminus \{1\}\),

$$\begin{aligned}& u_{n}=u_{n-m}^{\alpha _{m}}u_{n-m-1}^{\beta _{m}}u_{n-m-2}^{\gamma _{m}}u_{n-m-3}^{ \delta _{m}},\quad n\ge m+3, \end{aligned}$$
(49)
$$\begin{aligned}& \alpha _{m}=\alpha _{m-1}+\beta _{m-1}, \qquad \beta _{m}=\gamma _{m-1},\qquad \gamma _{m}=\delta _{m-1},\qquad \delta _{m}=\alpha _{m-1}. \end{aligned}$$
(50)

From (50) we obtain

$$\begin{aligned} \alpha _{n}=\alpha _{n-1}+\alpha _{n-4}, \end{aligned}$$
(51)

and \(\alpha _{0}=1\), \(\alpha _{-1}=\alpha _{-2}=\alpha _{-3}=0\), \(\alpha _{-4}=1\), \(\alpha _{-5}=\alpha _{-6}=0\).

Relations (49) and (50) yield

$$\begin{aligned} u_{n}=u_{3}^{\alpha _{n-3}}u_{2}^{\beta _{n-3}}u_{2}^{\gamma _{n-3}}u_{0}^{ \delta _{n-3}}=u_{3}^{\alpha _{n-3}}u_{2}^{\alpha _{n-6}}u_{2}^{ \alpha _{n-5}}u_{0}^{\alpha _{n-4}},\quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(52)

It is not difficult to see that the roots \(s_{j}\), \(j=\overline {1,4}\), of the polynomial

$$\begin{aligned} S_{4}(s)=s^{4}-s^{3}-1=0 \end{aligned}$$
(53)

are simple. From this and Lemma 1 we have that

$$\begin{aligned} \alpha _{n}=\sum_{j=1}^{4} \frac{s_{j}^{n+3}}{S_{4}'(s_{j})},\quad n \in {\mathbb {Z}}, \end{aligned}$$
(54)

is the solution to (51) such that \(\alpha _{-3}=\alpha _{-2}=\alpha _{-1}=0\), and \(\alpha _{0}=1\).

From (16) and (52) we get our next result.

Theorem 6

Let \(k=4\), \(l=3\), and \(a\in {\mathbb {C}}\setminus \{0\}\). Then the formula

$$\begin{aligned} y_{n}&= \frac{(1+ay_{3})^{\alpha _{n-3}}(1+ay_{2})^{\alpha _{n-6}}(1+ay_{1}) ^{\alpha _{n-5}}(1+ay_{0})^{\alpha _{n-4}}-1}{a},\quad n\in {\mathbb {N}}_{0}, \end{aligned}$$

where \(\alpha _{n}\) is defined in (54), presents a general solution to (10) in this case.

Besides the above-described six cases where equation (10) is solvable in closed form, there are some other ones. Here we present one.

2.7 Case \(l=0\)

If \(l=0\), then (10) is an equation with interlacing indices [40, 42]. If

$$ y_{m}^{(i)}=y_{mk+i} $$

for \(m\in {\mathbb {N}}_{0}\), \(i=\overline {0,k-1}\), then

$$\begin{aligned} y_{m+1}^{(i)}=2y_{m}^{(i)}+a \bigl(y_{m}^{(i)}\bigr)^{2}, \quad m\in {\mathbb {N}}_{0}, \end{aligned}$$
(55)

for each \(i\in \{0,\ldots ,k-1\}\), implying that \((y_{m}^{(i)})_{m\in {\mathbb {N}}_{0}}\), \(i=\overline {0,k-1}\), are solutions to (10) with \(k=1\). Thus the solution to (10) consists of k solutions to the equation with \(k=1\). So, it is of interest to find a general solution in this case.

2.8 Case \(k=1\), \(l=0\)

Since \(k=1\), it follows that

$$\begin{aligned} 1+ay_{n+1}=(1+ay_{n})^{2},\quad n\in {\mathbb {N}}_{0}, \end{aligned}$$
(56)

that is, \(u_{n+1}=u_{n}^{2}\), \(n\in {\mathbb {N}}_{0}\), and, consequently,

$$\begin{aligned} u_{n}=u_{0}^{2^{n}}, \quad n\in {\mathbb {N}}_{0}. \end{aligned}$$
(57)

Combining (16) and (57), we obtain the following theorem.

Theorem 7

Let \(k=1\), \(l=0\), and \(a\in {\mathbb {C}}\setminus \{0\}\). Then the formula

$$\begin{aligned} y_{n}&=\frac{(1+ay_{0})^{2^{n}}-1}{a},\quad n\in {\mathbb {N}}_{0}, \end{aligned}$$

presents a general solution to (10) in this case.

2.9 Case \(k\in {\mathbb {N}}\), \(l=0\)

From Theorem 7 we obtain the following result in the general case.

Theorem 8

Let \(k\in {\mathbb {N}}\), \(l=0\), and \(a\in {\mathbb {C}}\setminus \{0\}\). Then the formulas

$$\begin{aligned} &y_{mk}=\frac{(1+ay_{0})^{2^{m}}-1}{a},\quad m\in {\mathbb {N}}_{0}, \\ &y_{mk+1}=\frac{(1+ay_{1})^{2^{m}}-1}{a},\quad m\in {\mathbb {N}}_{0}, \\ & \vdots \\ &y_{mk+k-1}=\frac{(1+ay_{k-1})^{2^{m}}-1}{a},\quad m\in {\mathbb {N}}_{0}, \end{aligned}$$

present a general solution to (10) in this case.

Remark 1

The roots of polynomials (32), (39), (46), and (53) can be found by radicals in some routine ways [12]. Hence we omit the calculations and leave them to the reader as simple exercises.

2.10 A brief overview of some classes of solvable difference equations and the place of equation (10) in the classes

The first nontrivial class of solvable difference equations appearing in the literature was (9) (the homogeneous linear with constant coefficients). The equations of small orders have been essentially solved for the first time by de Moivre [6–8] by using generating functions, whereas the equations of any order were solved by Bernoulli [3] by looking for solutions in the form \(y_{n}=\lambda ^{n}\). As an example, in [3], a closed-form formula for the Fibonacci sequence also appeared. Some formulas for solutions to the linear equations can be also found in Euler’s book [11]. Lagrange [19] proposed the method of decomposition of the linear equations in linear factors and solved therein a first-order nonhomogeneous equation. He also proposed the method of undetermined coefficients for solving nonhomogeneous equations [18]. Another method for solving first-order nonhomogeneous equations was given by Laplace [20]. Therein he also solved, among other equations, some nonlinear ones, a few cyclic systems, and several partial difference equations. The paper indirectly shows that he also knew how to solve bilinear/fractional difference equations. Reduction of product-type difference equations to linear ones by using the logarithm has been known to mathematicians of the eighteenth century (see, e.g., [5], [16, p. 204]). We should say that the method is justified only for positive solutions to product-type difference equations. Some methods dealing with complex-valued solutions to product-type equations are described, for example, in [30, 33, 34] (see also the related references therein). As we have shown here, equation (10) is transformed to a product-type difference equation, so it is a close relative of product-type equations, and, consequently, a bit farther relative of the homogeneous linear difference equations with constant coefficients.