Experimenting with Discrete Dynamical Systems

George SPAHN and Doron ZEILBERGER

Dedicated to Saber Elaydi on his 80^{th} birthday, and to Gerry Ladas on his 85^{th} birthday

Abstract: We demonstrate the power of Experimental Mathematics and Symbolic Computation to study intriguing problems on rational difference equations, studied extensively by Difference Equations giants, Saber Elaydi and Gerry Ladas (and their students and collaborators). In particular we rigorously prove some fascinating conjectures made by Amal Amleh and Gerry Ladas back in 2000. For other conjectures we are content with semi-rigorous proofs. We also extend the work of Emilie Purvine (formerly Hogan) and Doron Zeilberger for rigorously and semi-rigorously proving global asymptotic stability of arbitrary rational difference equations (with positive coefficients), and more generally rational transformations of the positive orthant of R^{k} into itself.

Keywords: Non-Linear Difference Equations, Discrete Dynamical Systems, Semi-Rigorous Proofs

2020 MSC: 39A23, 39A30

Preface: Difference Equations, linear, but especially non-linear, are not only so useful in mathematical biology and elsewhere, but are fascinating to study for their own sake. Recall that already the first order logistic difference equation

x_{n+1} = λx_{n}(1 − x_{n}) , 0 ≤ λ ≤ 4 ,

introduced by Sir Robert May, lead to chaos, one of the central paradigms of our time, as well as to period-doubling, and the Feigenbaum constants. Here the action takes place in the finite interval 0 ≤ x ≤ 1.

Often in population dynamics one encounters higher-order difference equation where x_{n+1} is a rational function of the previous k values.

x_{n+1} = \frac{a_{0} + ∑_{i=1}^{k} a_{i}x_{n+1−i}}{b_{0} + ∑_{i=1}^{k} b_{i}x_{n+1−i}} ,

where the coefficients a_{0}, a_{1}, …, a_{k} and b_{0}, b_{1}, …, b_{k} are assumed non-negative (and of course at least one of them is strictly positive at the bottom), and with positive initial conditions x_{1}, …, x_{k}. Here the action takes place in the infinite interval 0 < x < ∞.

These difference equations have been studied extensively, with a deep theory, by the two ‘birthday boys’ and their many disciples, see [AL],[CL], [E1], [E2], [KoL], [KuL], and references thereof.

Twelve years ago one of us (DZ), in collaboration with his then PhD student, Emilie Hogan (now Purvine), initiated the use of symbolic computation, and computer-generated (rigorous!) proofs to study such difference equations, but only those whose solutions always converge to a unique equilibrium point, i.e. for which there exists a unique x̄ > 0 such that for any (positive) initial conditions, x_{1}, …, x_{k}, one has lim_{n→∞} x_{n} = x̄, and the challenge was to prove it rigorously. Note that to prove that a candidate fixed point x̄ is locally stable can be easily done using standard techniques (see below).

In [AL] several fascinating conjectures were made (with monetary prizes offered by Gerry Ladas for some of them, alas with the very ungenerous deadline of Jan. 1, 2002, long passed). Here is the first one:

Conjecture ([AL]): Prove that every positive solution of the difference equation

x_{n+1} = \frac{x_{n−1}}{x_{n−1} + x_{n−2}} ,

converges to a period two solution of that equation of the form

…, ϕ, 1 − ϕ, …,

with 0 ≤ ϕ ≤ 1.

In this paper we extend, and improve, the method of [HZ], and then extend this methodology to rigorously prove (with the aid of our beloved silicon servants) some of the fascinating conjectures in [AL]. For other ones we will be content with semi-rigorous proofs (see [Z] for the concept). We will argue that often such proofs suffice, since we know that there exists a rigorous proof (or disproof), but it would be a waste of the computer’s time and memory to find it, since the probability that the semi-rigorous proof was a false positive is negligible.

Accompanying Maple packages

This article is accompanied by two Maple packages.

• DRDS.txt, to experiment, numerically, and symbolically, with solutions of rational difference equations of any order, and for proving, if possible, rigorously, but more often (if the order is three and up) semi-rigorously, global asymptotic stability. This latter part is a continuation, and improvement, of the pioneering work in [HZ], that focused on second-order difference equations.

• AmalGerry.txt, to prove rigorously (and in some cases, semi-rigorously), some of the intriguing and tantalizing conjectures made by Amal Amleh and Gerry Ladas in 2001 [AL].

Both packages, and numerous input and outputs files, are viewable (and downloadable!) from the front of this article

https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/dds.html .

Experimenting with Some Random Rational Difference Equations One of the purposes of this article is to serve as a tutorial on our Maple package DRDS.txt, that in addition to proving can be also used as a calculator for exploring and experimenting, numerically, with rational difference equations. To get a feel of how a typical rational difference equation looks like, use procedure RRDE(x,n,k,d,A), where the input parameters are:

• x and n are symbols that correspond to the notation x_{n} in ‘humanese’.

• k is the order of the difference equation

• d is the degree of the numerator and denominator (in this paper we will focus on the degree one case, but the Maple package can handle any degree).

• A is a positive integer. The output is such a random rational difference equation whose numerator and denominators have integer coefficients drawn from {1, …, A}. For example, typing

F:=RRDE(x,n,3,2,30); ,

may give you something like this (of course, every time you get a different difference equation, since this is random)

x_{n+1} = \frac{17x_{n}^{2} + 2x_{n}x_{n−2} + 20x_{n}x_{n−1} + 21x_{n−2}^{2} + 25x_{n−2}x_{n−1} + 17x_{n−1}^{2} + 8x_{n} + 23x_{n−2} + 6x_{n−1} + 13}{4x_{n}^{2} + 4x_{n}x_{n−2} + 29x_{n}x_{n−1} + 10x_{n−2}^{2} + 4x_{n−2}x_{n−1} + 4x_{n−1}^{2} + 9x_{n} + 6x_{n−2} + 19x_{n−1} + 9} .

(1)

To get the first N + k terms of such a difference equation with given initial conditions, type:

L:=OrbD(F,x,n,INI,N);

where F is the difference equation in the above format x_{n+1} = F, x and n are the symbols used to express it, INI is a list of of length k, and N is the number of extra terms (so the output list is of length N + k). For example, typing:

OrbD(x[n]+x[n-1],x,n,[1,1],10);

would yield

[1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144] .

Going back to the above complicated (random) third-order difference equation denoted by F above, typing

L:=OrbD(F,x,n,evalf([21,27,39]),1000):

would give you in floating-point, the first 1003 terms of the sequence with initial conditions x_{1} = 21, x_{2} = 27, x_{3} = 39. Note the colon :, as opposed to the semi-colon, ;, since we really don’t want to see all of them, we are only interested in the long-term behavior, i.e. whether ultimately it converges to one number (or in the long-run to a periodic solution (see below)).

Typing L[-1]; and L[-2]; would give respectively, the 1002^{th} and 1003^{th} terms:

1.6358881124186260402 …, 1.6358881124186260402 …,

indicating that they are very close, hence, at least with the above initial conditions, the sequence seems to have a limit. Experimenting with randomly chosen (positive) initial conditions, we see again and again, that it seems to converge to the same number. Hence, completely based on numerics, we are safe to make the following conjecture.

Random Conjecture: For any positive initial conditions x_{1}, x_{2}, x_{3}, the terms of the sequence satisfying the difference equation (1) converge to a certain algebraic number (that can be easily found), whose floating-point value is 1.6358881124186260402 .

So far we only used numerical computations. Later on we will use symbolic computation to actually prove it, either rigorously or semi-rigorously.

Experimenting with many other random cases, and several initial conditions, you get again and again this phenomenon that the sequences seem to converge to the same number, let’s call it x̄, i.e. that it is a global equilibrium. To find out its value, one replaces x_{n+1}, x_{n}, x_{n−1}, …, x_{n−k+1} by x̄ getting a one-variable polynomial equation in x̄, asking Maple to solve it, and if all goes well only getting one real and positive solution (of course, the resulting equation would have d + 1 complex roots).

Let’s pick another random example, this time a third-order difference equation, and to make it simpler, let’s have the numerator and denominator of degree 1. Typing

T:=RRDE(x,n,3,1,30);

yielded (this time)

x_{n+1} = \frac{17 + 5x_{n} + 24x_{n−1} + 16x_{n−2}}{23 + 4x_{n} + 19x_{n−1} + 2x_{n−2}} .  (2)

Taking random initial conditions (in this case x_{1} = 11, x_{2} = 27, x_{3} = 37, and typing

L:=OrbD(T,x,n,evalf([11,27,37]),1000):  L[-1],L[-2];

gives:

1.37466571564383380885297479 …, 1.37466571564383380885297479 …,

and similarly for many other randomly chosen initial conditions. To find the exact value of this (so far conjectured) equilibrium point, solve

x̄ = \frac{17 + 5x̄ + 24x̄ + 16x̄}{23 + 4x̄ + 19x̄ + 2x̄} ,

getting

x̄ = \frac{17 + 45x̄}{23 + 25x̄} ,

that simplifies to the quadratic equation:

25x̄^{2} − 22x̄ − 17 = 0 ,

whose roots are:

[\frac{11}{25} + \frac{\sqrt{546}}{25}, \frac{11}{25} − \frac{\sqrt{546}}{25}] ,

that in decimals are:

[1.374665716, −0.4946657156] ,

discarding the negative root, we got the exact value, namely \frac{11}{25} + \frac{\sqrt{546}}{25}. The Amleh-Ladas Fascinating conjectures Procedure RRDE artificially made all coefficients strictly positive, and hence it turns out that, generically, one gets rather boring limiting behavior, i.e. convergence to a unique positive equilibrium point, or phrased otherwise, a limiting period-one solution. In [AL], Amal Amleh and Gerry Ladas made the following intriguing conjectures, that exhibited far more interesting long-term behavior. Conjecture 1: For any positive initial conditions x_{1}, x_{2}, x_{3},

x_{n+1} = \frac{x_{n−1}}{x_{n−1} + x_{n−2}} , n ≥ 3 ,

the sequence {x_{n}} converges to a period-two solution of the form

…, ϕ , 1 − ϕ , …,

with 0 ≤ ϕ ≤ 1. Conjecture 2: For any positive initial conditions x_{1}, x_{2}, x_{3},

x_{n+1} = \frac{x_{n} + x_{n−2}}{x_{n−1}} , n ≥ 3 ,

the sequence {x_{n}} converges to a period-four solution of the form

…, ϕ , ψ , \frac{ϕ + ψ^{2}}{ϕψ − 1} , \frac{ϕ^{2} + ψ}{ϕψ − 1} …,

with ϕ, ψ ∈ (0, ∞), and ϕψ > 1. Conjecture 3: For any positive initial conditions x_{1}, x_{2}, x_{3},

x_{n+1} = \frac{1 + x_{n−2}}{x_{n}} , n ≥ 3 ,

the sequence {x_{n}} converges to a period-five solution of the form

…, ϕ , ψ , \frac{1 + ϕ}{ϕψ − 1} , ϕψ − 1 , \frac{1 + ψ}{ϕψ − 1} …

with ϕ, ψ ∈ (0, ∞), and ϕψ > 1. Conjecture 4: For any positive initial conditions x_{1}, x_{2}, x_{3},

x_{n+1} = \frac{1 + x_{n}}{x_{n−1} + x_{n−2}} , n ≥ 3

The sequence {x_{n}} converges to a period-six solution of the form

…, ϕ , ψ , \frac{ϕ}{ψ} , \frac{1}{ϕ} , \frac{1}{ψ} , \frac{ϕ}{ψ} , …,

with ϕ, ψ ∈ (0, ∞). We will later show how to prove these using symbolic computation, but for now, let’s confirm them numerically, using procedure OrbD. These difference equations are hard-coded in procedure LadadDB(x,n), that contains 15 interesting difference equations, the first four, namely

LadasDB(x,n)[1]   ,   LadasDB(x,n)[2]   ,   LadasDB(x,n)[3]   ,   LadasDB(x,n)[4],

correspond to the difference equations featured in the above four conjectures. For example, typing

T:=LadasDB(x,n)[1]:  evalf(OrbD(T,x,n,evalf([1,2,3]),1000)[-3..-1],10);

gives

[0.7012220196, 0.2987779809, 0.7012220196] ,

while

T:=LadasDB(x,n)[1]:  evalf(OrbD(T,x,n,evalf([11,25,34]),1000)[-3..-1],10);

yields

[0.9348089961, 0.06519100396, 0.9348089961] .

Readers are welcome to experiment with many other initial conditions, and similarly for the other difference equations featured in Conjectures 2,3, and 4, to numerically (empirically) confirm these intriguing conjectures. Of course this ‘only’ gives numerical confirmation (that the authors of [AL] must have already done, but probably using numerical software rather than Maple). We will later see how to prove them either rigorously or semi-rigorously.

Recalling some basics and More Numerical Explorations

While it is highly non-trivial, in general, to prove that every choice of initial conditions will make the solution sequence converge to an equilibrium, it is purely routine, today, to decide whether it is true when the initial conditions are not too far from that equilibrium.

The first step (already recalled in [HZ]) is to convert a k-th order difference equation in (0, ∞) to a first-order difference equation in (0, ∞)^{k}.

The k-th order difference equation

x_{n+1} = F (x_{n}, x_{n−1}, …, x_{n−k+1})

becomes the transformation

[ \begin{matrix}x_{1} \\ x_{2} \\ … \\ x_{k} \end{matrix} ] → [ \begin{matrix}x_{2} \\ x_{3} \\ … \\ x_{k} \\ F (x_{k}, x_{k−1}, x_{k−2}, …, x_{1}) \end{matrix} ] .

From this point of view, the more general problem is to investigate whether, given a general rational transformation from (0, ∞)^{k} into (0, ∞)^{k} of the form

[ \begin{matrix}x_{1} \\ x_{2} \\ … \\ x_{k} \end{matrix} ] → [ \begin{matrix}R_{1}(x_{1}, …, x_{k}) \\ R_{2}(x_{1}, …, x_{k}) \\ … \\ R_{k}(x_{1}, …, x_{k}) \end{matrix} ] ,

where R_{1}(x_{1}, …, x_{k}), …, R_{k}(x_{1}, …, x_{k}) are rational functions of their arguments. In order to prove that x̄ is the limit of every solution of the original difference equation, one has to prove that the orbit, starting at any point in (0, ∞)^{k} of the resulting transformation (as constructed above) converges to (x̄, x̄, …, x̄).

Procedure Targem in the Maple package DRDS.txt converts a difference equation to a transformation.

Recall that for any transformation x → F(x) in R^{k}, a point x_{0} is called an equilibrium point if it is a fixed point.

F(x_{0}) = x_{0} .

If the transformation is, like in our case, rational, this gives a system of k polynomial equations with k unknowns, that at least, in principle, but often also in practice (for small k) is fully solvable (in the realm of algebraic numbers, if all the coefficients are integers). In order to generate random examples to experiment with, use procedure

RRT(x,k,d,A),

where x is a symbol, k is the dimension of the space, d is the degree of the numerator and the denominator, and A is a positive integer, such that the coefficients are randomly chosen from {1, 2, …, A}. For example,

T:=RRT(x,2,1,30); ,

might give

[\frac{18 + 20x_{1} + 24x_{2}}{11 + 19x_{1} + 25x_{2}} , \frac{26 + 29x_{1} + 28x_{2}}{29 + 14x_{1} + 18x_{2}} ] .

To get all the fixed points (including complex and with negative coordinates), type

FP(T,x);

In this example, you would get a big mess, so let’s take the floating-point version, and type

evalf( FP(T,x),10);

getting

[[0.5983411214, -1.834086341], [1.100408318, 1.394961226], [-0.9483476940, 0.159782893],
[-100.8951386, 76.30710943]] .

We are only interested in points in (0, ∞)^{2}, so the only point that we are interested in is the second one. To get such points with all positive coordinates right away, type

evalf( FPp(T,x),10);

getting, indeed,

{[1.100408318 …, 1.394961226 …]} .

We next ask whether it is locally stable. As is well-known, and fairly easy to see (e.g. [KL]), one computes the Jacobian matrix, then for each of the candidate points (those with positive coordinates) plugs-it-in then computes the eigenvalues of this numerical matrix, and if all of them have absolute value less than one, then we know for sure that the examined equilibrium point is locally stable. This is an important first step before we can prove global stability, since every globally stable fixed point must be, first of all, a local one. Of course if there are more than one locally stable fixed points in (0, ∞)^{k} there is no hope that any of them is global, but it is still nice to know all of them. This is all done, thanks to Maple, automatically, with procedure LSFP(T,x). For example, with the above T, we would get that this point is indeed locally stable. It also gives you the (floating-point) versions of the eigenvalues (just for the record). So typing

evalf(LSFP(T,x),10);

gives

{[[1.100408318, 1.394961226], [0.01399544579+0.07999899783 i, 0.01399544579−0.07999899783 i]]} ,

indicating that indeed there is only one locally stable fixed point of our transformation, and also giving the eigenvalues. Before trying to prove rigorously (that takes lots of effort!), or even semi-rigorously (that also takes some effort, see below) it is very easy to conjecture whether there (most probably) is a globally stable fixed point, and to actually find it. In fact, we don’t need to find the locally stable fixed points. The completely numeric and empirical procedure

CoGSFP(T,x,K1,K2) ,

takes as input a transformation T and picks K2 random points in (0, ∞)^{k}, and for each of them computes the orbit of length K1. If for each of these orbits the difference between the last two terms is very small, and they all give the same point, we know non-rigorously, but (almost) certainly that there is a globally stable fixed point, and we know its floating point version. To get its exact value, as an algebraic number, you need FPp mentioned above. The Proof Strategy This is essentially what was done in [HZ], but with a new implementation, and with a semi-rigorous option. Suppose that we have a (rational) transformation, T, from (0, ∞)^{k} to (0, ∞)^{k} , and that we already have ample empirical/numerical evidence that a certain candidate point, (gotten either from LSFP or CoGSFP), x̄, is a globally stable fixed point of T in (0, ∞)^{k}. Our goal in life is to prove that for every initial point, x_{0} ∈ (0, ∞)^{k},

lim_{n→∞} T^{n}(x_{0}) = x̄ .

Let |x| be the usual Euclidean norm. The above statement is equivalent to

lim_{n→∞} |T^{n}(x_{0}) − x̄|^{2} = 0 .

Suppose that we can come-up with a positive integer r, and some real α > 1, such that, for any point x ∈ (0, ∞)^{r}, we have the inequality

α |T^{r}(x) − x̄|^{2} ≤ |x − x̄|^{2} ,

then we know that this sequence of ‘distance-squared from the fixed point’ shrinks (by at least the factor 1/α), every r-th iteration compared to what it was. This would automatically entail what we want. Note that this is only a sufficient condition, and there is no a priori reason (that we know of, at least), that such a real α > 1 and an integer r ≥ 1 exist, but if we are lucky enough to find a candidate, and then, succeed in proving it, we are be done!

Again, we first investigate things numerically. We start with r = 1, and then see whether for many initial points, the resulting orbit has the property that the distance-squared from x̄ shrinks every r-th iteration. Once we get a successful candidate r ≥ 1, we have to prove it. For the sake of definiteness, and not to clutter the Maple code with another parameter, we decided to take α = \frac{101}{100}.

So we have to prove (either rigorously or semi-rigorously) that for any x ∈ (0, ∞)^{k}, we have the inequality

α|T^{r}(x) − x̄|^{2} ≤ |x − x̄|^{2} .

in other words

|x − x̄|^{2} − α|T^{r}(x) − x̄|^{2} ≥ 0 .

Now the left-side is a (usually rather complicated) rational function of x_{1}, …, x_{k}, (recall that x = (x_{1}, …, x_{k})). Simplifying, we get a denominator that is a perfect square, and hence automatically positive. The numerator is a certain (usually complicated) polynomial, and everything boils down to proving that this polynomial, let’s call it P (x_{1}, …, x_{k}), is non-negative. In other words we need to minimize P in the region (0, ∞)^{k} and prove that it is ≥ 0 (in fact if it is ≥ 0 it must be 0, since the value of the above left side is 0 when x = x̄).

But this is a routine multivariable calculus exercise!, that Maple (and Mathematica, and Sage), know how to do. Alas, in applications to our problems, this polynomial turns out to be too complicated for k > 2, so we can rigorously prove global stability for k = 2, but it would take too long (on our modest laptops) to do it for k ≥ 3. So we opt to do it numerically, checking it for many random points, and since we know that there exists a way, at least in principle, to prove it rigorously, why bother? This semi-rigorous approach to mathematical deduction was proposed by one of us thirty years ago [Z]. To get a fully rigorous proof, use procedure GSFPv(T,x,K);, where K is the maximum r we are willing to take. That works well with two dimensions, but for higher dimensions, it takes way too long. One should use instead GSFPvSR(T,x,K);, in order to get a semi-rigorous proof. If none is found, it returns FAIL. Sample output for proving (Rigorously and Semi-Rigorously Global Stability)

• If you want to see 20 theorems that state that certain second-order difference equations always converge to the unique stable equilibrium, with fully rigorous proofs, look here:

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oDRDS1.txt .

• If you want to see 10 theorems that state that certain third-order difference equations always converge to the unique stable equilibrium, with semi-rigorous proofs, look here:

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oDRDS2.txt .

• If you want to see 5 theorems that state that certain fourth-order difference equations always converge to the unique stable equilibrium, with semi-rigorous proofs, look here:

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oDRDS3.txt .

(Warning: the file is large!) The Amleh-Ladas Conjectures We now apply the proof strategy described above to the fascinating Amleh-Ladas conjectures. Recall conjecture 1: Conjecture 1: For any positive initial conditions x_{1}, x_{2}, x_{3},

x_{n+1} = \frac{x_{n−1}}{x_{n−1} + x_{n−2}} , n ≥ 3 ,

the sequence {x_{n}} converges to a period-two solution of the form

…, ϕ , 1 − ϕ , …,

with 0 ≤ ϕ ≤ 1. To tackle this problem, we first view the sequence as a map from R^{3} to R^{3}.

T : (x, y, z) → (y, z, \frac{y}{x + y})

We hope to show that for any start point (x_{0}, y_{0}, z_{0}) with x_{0}, y_{0}, z_{0} > 0, this dynamical system, T converges to a point somewhere on the line segment parametrized by (t, 1 − t, t), t ∈ [0, 1].

To measure how close points are to this line, we use the norm

v(x, y, z) = 1 + x^{2} + y^{2} + z^{2} − x − 2y − z + xy − xz + yz

. This is the square of the Euclidean distance from (x, y, z) to the line. To show that T converges to a point with v-norm equal to 0, we consider the objective function

F (x, y, z) = v(x, y, z) − v(T (x, y, z))

If we can show that this objective function is always ≥ 0 for all positive (x, y, z), then we conclude that applying T cannot increase the v-norm. Unfortunately this objective function is sometimes negative for this particular T and v. To fix it, we replace T with T^{3}, three consecutive applications of T. We also apply T to both points for smoothing. This gives

F (x, y, z) = v(T (x, y, z)) − v(T^{4}(x, y, z))

.

It is now time to set maple to work and show that the objective function is nonnegative. We use the function NLPSolve from the Optimization package, to minimize the objective function. This function requires an initial point, and uses iterative methods to search for improvements on the initial point using floating point arithmetic. The maple documentation for this function is given here:

https://www.maplesoft.com/support/help/maple/view.aspx?path=Optimization%2FNLPSolve

We ran the solver with 36 different initial points and each time it returned that the minimum was within an acceptable range (±10^{−6}) of 0. This completes a semi-rigorous proof of Conjecture 1. To run our code yourself, download the maple package AmalGerry.txt, and type:

run nlp(T4,v4,1,4,3); .

The first argument asks for the transformation on R^{3}. T4 is the transformation T defined above; the 4 is a reference to the fact that this transformation was equation 4 on the Amleh-Ladas paper. The second argument is the norm to be used; v4 corresponds to T4. The third and fourth arguments specify how many iterations of T should be applied to the initial point when creating the objective function. The last argument determines the amount of different initial points that are tested.

We used this method to semi-rigorously prove conjectures 1 through 4. Take a look at the maple package AmalGerry.txt for more details!

Using Maple’s symbolic minimizer

The NLPSolve function is not a mathematical proof of minimization. Perhaps it could be made into one by taking a very close look at the algorithm in the code, but this would be very tedious. In this section we attempt to use maple’s built-in symbolic minimizer to produce a fully rigorous proof of the minimization. The main limitation here is computational resources, so we experiment with different norms and parameters.

Still looking at the T and v from the previous section, maple is not able to symbolically compute the minimum of

F (x, y, z) = v(T (x, y, z)) − v(T^{4}(x, y, z))

in a reasonable amount of time. Instead we define simpler norms

v_{xy} = (x + y − 1)^{2}
v_{yz} = (y + z − 1)^{2}

Let

F_{xy} = v_{xy}(x, y, z) − v_{xy}(T^{4}(x, y, z))

The denominator of F_{xy} turns out to be

(xz + yz + y)^{2}(y + z)^{2}

which is never negative! Thus we just instruct maple to minimize the numerator, which is a degree 8 polynomial in the variables x, y, z. Maple symbolically computes that the minimum of this polynomial is 0, so we have a rigorous proof that T converges to a point with v_{xy} norm equal to 0. The same process works for v_{yz}. The only points where v_{xy} and v_{yz} are both 0 is exactly the line (t, 1 − t, t). This completes a rigorous mathematical proof of Conjecture 1. The authors attempted to use a similar approach for the other conjectures however we lacked the computational resources to execute the code.

Periodic Difference Equations

As mentioned in [AL] (p. 71), and [KL] (p. 628), the following (second-order) Lynnes difference Equation

x_{n+1} = \frac{1 + x_{n}}{x_{n−1}}

always has period five, regardless of the initial conditions. This is easily confirmed with our Maple package:

Entering:

OrbD((1+x[n])/x[n-1],x,n,[a,b],5);

immediately gives

[a, b, \frac{1 + b}{a}, \frac{a + 1 + b}{ab}, \frac{a + 1}{b}, a, b] ,

meaning that for symbolic, i.e. general, initial conditions, things get repeated every five iterations. Of course this is easy enough to do by hand. Also mentioned in [KL] (p. 628, Eq. (25) there), that the following (third-order) difference equation

x_{n+1} = \frac{1 + x_{n} + x_{n−1}}{x_{n−2}} ,

is always of period eight, regardless of the initial conditions. Indeed, typing:

OrbD((1+x[n]+x[n-1])/x[n-2],x,n,[a,b,c],8);

yields

[a, b, c, \frac{1 + c + b}{a}, \frac{ca + a + b + c + 1}{ab}, \frac{ab + ca + b^{2} + bc + a + 2b + c + 1}{abc}, \frac{ca + a + b + c + 1}{bc}, \frac{a + 1 + b}{c}, a, b, c] .

This gave us the hope that the fourth-order difference equation

x_{n+1} = \frac{1 + x_{n} + x_{n−1} + x_{n−2}}{x_{n−3}} ,

is perhaps periodic? Alas, entering

L:=OrbD((1+x[n]+x[n-1]+x[n-2])/x[n-3],x,n,[1,1,1,1],1000):  member(1,op(5..nops(L),L));

gives false, meaning, that if there is a period, it would be larger than 1000, so this difference equation is unlikely to be periodic. It would be very interesting to discover such rational difference equations with higher periods, that do not trivially follow from the known ones by ‘merging’. Automated Discovery of Invariants that Imply that Every Solution is Bounded In the Kulenovic-Ladas fascinating book [KL], it is mentioned that the generalized Lynnes Equation

x_{n+1} = \frac{p + x_{n}}{x_{n−1}} ,

for any positive p has the following invariant:

I_{n} = (p + x_{n−1} + x_{n}) (1 + \frac{1}{x_{n−1}} ) (1 + \frac{1}{x_{n}} ) .

We could have found it, ab initio, using procedure FindInv:

FindInv((p+x[n])/x[n-1],x,n,2,3);

We were able to find invariants for the higher-order difference equations

x_{n+1} = \frac{p + x_{n} + x_{n−1} + … + x_{n−k+2}}{x_{n−k+1}} ,

for k ≤ 8. They do get more and more complicated, see the output file

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oDRDS6.txt .

They all turn out to have positive coefficients. As noted in [KL] for the generalized Lynnes equation, but is also true in general, the existence of an invariant of the form

\frac{P (x_{n}, x_{n−1}, …, x_{n−k+1})}{x_{n}x_{n−1} … x_{n−k+1}} = Constant ,

with the coefficients of P all positive (all the ones we found were of that form) immediately implies that for any positive initial conditions, the solution sequence is always bounded. In fact, for the generalized Lynnes equation x_{n+1} = \frac{p+x_{n}}{x_{n−1}} , one can use discriminants to predict, a priori, these lower and upper bounds, for any positive p and any positive initial conditions x_{1} = a_{1}, x_{2} = a_{2}. See the output file

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oDRDS4.txt .

For the analogous third-order difference equation x_{n+1} = \frac{p+x_{n}+x_{n−1}}{x_{n−2}} , see:

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oDRDS5.txt .

Conclusion Using the great human-generated theory developed by Saber Elaydi, Gerry Ladas, and their many collaborators and disciples, and bringing into the game both numeric and symbolic computation, we hope that we demonstrated the power of computer-kind to extend the human efforts. Alas, even computers have their limits, and we advocate that often a semi-rigorous proof suffices, as first preached in 1993 in [Z]. References

[AL] A.M. Amleh and G. Ladas, Convergence to Periodic Solutions, Journal of Difference Equations and Applications 7(2001), 621-631.

[CL] E. Camouzis and G. Ladas, “Dynamics of Third Order Rational Difference Equations”, Chapman and Hall/CRC press (2008).

[E1] Saber Elaydi, “Introduction to Difference Equations”, Springer (2000).

[E2] Saber Elaydi, Global dynamics of discrete dynamical systems and difference equations, Difference equations, discrete dynamical systems and applications, 51-81, Springer Proc. Math. Stat. 287, Springer, 2019.

[HZ] Emilie Hogan (now Purvine) and Doron Zeilberger, A New Algorithm for Proving Global Asymptotic Stability of Rational Difference Equations, Journal of Difference Equations and Applications 18(2012), 1853-1873.

https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/gas.html .

[KoL] V. Kocic and G. Ladas, “Global Behavior of Nonlinear Difference Equations of Higher Order with Applications”, Kluwer Academic Publishers (1993).

[KuL] M. R. S. Kulenovic and G. Ladas, “Dynamics of Second Order Rational Difference Equations, Chapman and Hall/CRC press (2001).

[Z] Doron Zeilberger, Theorems for a price: tomorrow’s semi-rigorous mathematical culture, Notices of the Amer. Math. Soc. 40, 978-981. Reprinted in: Math. Intell. 16 (1994) 11-14.

https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/priced.html .

George Spahn and Doron Zeilberger, Department of Mathematics, Rutgers University (New Brunswick), Hill Center-Busch Campus, 110 Frelinghuysen Rd., Piscataway, NJ 08854-8019, USA.

Email: gs828 at math dot rutgers dot edu , DoronZeil] at gmail dot com .

Written: June 21, 2023.