Closed form of $x_{n+1}=x_n-1/x_n$ with $x \neq 0,1$

Question

I have the sequence

$x_{n+1}=x_n-1/x_n$ with $x \neq 0,1$

For certain values of $x_0$ this sequence will be periodic. To find these values one can look at the root of: $x_{p} = x_0$ where $p$ is the period. For $p=2$ this can be solved easily.

For $p>10$ it's not viable to define the function recursively and search for roots numerically. To look for the roots at higher values of $p$ I'm searching for a closed form of this formula.

I already found an answer to a similar problem, but failed to transfer the solution to my different sequence. I tried using the same approach.

By $x=\cot(\theta)$ you get: $x_1=2\cot(2\theta)=\cot(\theta)-\frac{1}{\cot(\theta)}$

I calculated $x_2$ by hand which is: $x_2=\cot(4\theta)+\frac{3}{2}\cot(2\theta)$
But I wasn't able to find a simple formula for $x_3$

Any help with finding a closed form (or an idea how to find the roots more easily) would be appreciated!

EDIT:
I moved my follow up questions to a seperate topic

Answer 1

The equations $x_{n+1} = x_n - \frac 1{x_n}$ and $x_{n+1} = \frac 12 \left( x_n - \frac 1{x_n} \right)$ are maybe similar in appearance but very different qualitatively. For instance, the first one has no fixed points while the second one does (in the complex numbers): $$ x = \frac 12 \left( x - \frac 1x \right) = \frac{x^2-1}{2x} \quad \Longrightarrow \quad 2x^2 = x^2 - 1 \quad \Longrightarrow \quad x^2 + 1 = 0 $$ and so it has the solutions $\pm i$. If the first one had fixed points, it would imply $x = x - \frac 1x$ has a solution, but this clearly has no solution.

Working on our original problem again but considering $x_{n+2}$ as a function of $x_n$ instead, we get $$ x_{n+2} = x_{n+1} - \frac 1{x_{n+1}} = \frac{x_{n+1}^2 - 1}{x_{n+1}} = \frac{\left(\frac{x_n^2-1}{x_n}\right)^2 - 1}{\frac{x_n^2 -1}{x_n}} = \frac{x_n^4 - 3x_n^2 + 1}{x_n(x_n^2-1)} $$ and if the recurrence equation $y_n = x_{2n}$ has a fixed point, then it satisfies $y_{n+1} = x_{2n+2} = x_{2n} = y_n$, so we can use the above equation and deduce $$ -\frac 14 = y^4 - y^2 = y^2(y^2-1) = y^4 - 3y^2 + 1 = \frac 14 - \frac 32 + 1 \quad \Longrightarrow \quad 2y^2 - 1 = 0 $$ which means $y = \pm \frac 1{\sqrt 2}$ are the possible values of $x_0$ for which the sequence $x_{2n}$ will be periodic.

You might wonder, if starting with $x_0 = \frac 1{\sqrt 2}$ means $x_0 = x_2 = \cdots$, what happens to $x_1,x_3,\cdots$? Well, since $x_{2n+1} = f(x_{2n})$, and $x_{2n} = x_0$ for all $n \ge 0$, it follows that $x_{2n+1} = x_1$ for all $n \ge 0$. This means that if $x$ is a solution to $f(f(x)) = x$, then $f(x)$ is also a solution to $f(f(x)) = x$, and because we know there are no solution to $f(x) = x$, it has to be a distinct one! We deduce the following result:

Let $p$ be a positive integer. If $x$ is a solution to the equation $f^{(p)}(x) = x$ (where $f^{(p)}(x)$ is obtained by applying $f$ on $x$ precisely $p$ times), then $x, f(x), \cdots, f^{(k)}(x), \cdots, f^{(p-1)}(x)$ are all solutions to the equation $f^{(p)}(x)$. Furthermore,

• if $q > 1$ is a positive integer that divides $p$, then $f^{(p/q)}(x)$ is a solution to the equation $f^{(q)}(x) = x$, so that if $x$ has period $p$, $f^{(p/q)}(x)$ has period $q$
• if $p$ is prime, then the set $\{x, f(x), \cdots, f^{(k)}(x), \cdots, f^{(p-1)}(x)\}$ are distinct solutions to the equation $f^{(p)}(x)$, and therefore are periodic solutions of period $p$

To be able to say a little bit more about the periods of order $p$, consider the following general pattern. Write $x_{n+p} = \frac{f_p(x_n)}{g_p(x_n)}$ where $f_p$ and $g_p$ are the polynomials in $x_n$ we obtain by iterating the above idea $p$ times instead of $2$ times. We saw that $$ \deg f_1 = 2, \quad \deg g_1 = 1, \quad \deg f_2 = 4, \quad \deg g_2 = 3. $$ Let us prove by induction on $n$ that $\deg f_p = 2^p$ and $\deg g_p = 2^p-1$. The base case is proved already, and $$ \frac{f_{p+1}(x_n)}{g_{p+1}(x_n)} = \frac{f_p(x_n)}{g_p(x_n)} - \frac{g_p(x_n)}{f_p(x_n)} = \frac{f_p(x_n)^2 - g_p(x_n)^2}{f_p(x_n)g_p(x_n)}. $$ Since $\deg g_p < \deg f_p$ by induction, no cancellation occurs in the numerator, so $\deg f_{p+1} = 2 \deg f_p = 2^{p+1}$. Then $\deg g_{p+1} = \deg f_p + \deg g_p = 2^p + 2^p - 1 = 2^{p+1}-1$.

This proves our induction hypothesis, so we can deduce the following result:

Let $p > 1$ be a positive integer. There are at most $2^p - 2$ solutions to the equation $f^{(p)}(x) = x$.

The proof is simple; the solutions to that equation are solution to the polynomial equation $f^p(x) - xg_p(x) = 0$, and this polynomial has degree at most $2^p$.

We can slightly improve our result by noticing that there is some cancellation involved in that subtraction. By induction on $n$, one easily sees that the leading coefficients of $f_p$ and $g_p$ are both $1$, so the leading term of $f_p(x)$ and $xg_p(x)$ are the same.

We can also see by induction on $n$ that $f_p(x)$ is an even polynomial function and $g_p(x)$ is an odd polynomial function, so $x g_p(x)$ is even. Even polynomial functions can be written as polynomials in $x^2$, meaning that the terms of odd degree of $f_p(x)$ and $xg_p(x)$ are all zero. This brings our upper bound down on the number of solutions to $2^p-2$.

I believe this upper bound is sharp (but this is conjectural) since using sympy, I computed
$$ f_3(x) = x^8 - 7x^6 + 13x^4 - 7x^2 + 1, \quad xg_3(x) = x^8 - 4x^6 + 4x^4 - x^2 $$ and so we see that there is no cancellation. I repeated the reasoning computationally and it seems to hold for $p \le 15$, where the polynomial has degree $2^p-2$.

Note that nothing says the polynomial doesn't have double roots, so we can't claim that there are always exactly $2^p-2$ solutions. Some of them might be repeated roots of that polynomial.

If I find more interesting results, I'll add to this answer.

EDIT: Adding here my sympy code to generate the polynomials $f_p$ and $g_p$.

from sympy.abc import x
def recurrence_relation():
    f = x
    g = 1
    while True:
            yield (f,g)
            f,g = f**2 - g**2, f*g

def first_n_terms(n):
    for (f,g), i in zip(recurrence_relation(), range(n)):
            if i <= 1:
                    continue
            print(f"Polynomial f_p: {f}")
            print(f"Polynomial xg_p: {x*g}")
            print(f"Degree of f_p - xg_p: {poly(f-x*g, x).degree()}")

I hope this tells you roughly how to operate with sympy for your own experiments.

Hope that helps,