If $f(x)=x^2-x-1$ and $f^n(x)=f(f(\cdots f(x)\cdots))$, find all $x$ for which $f^{3n}(x)$ converges.

Let $f:\mathbb{R}\to\mathbb{R}$ be the polynomial defined by $$f(x)=x^2-x-1$$ and let $$g_0(x)=f(x),\quad g_1(x)=f(f(x)),\quad\ldots\quad g_n(x)=f(f(f(\cdots f(x)\cdots)))$$ The positive root of $f(x)$ is the famous golden ratio $$\varphi=\frac{1+\sqrt{5}}{2}.$$

Problem: Find all $x\in\mathbb{R}$ for which the limit $$\lim_{n\to\infty}g_{3n}(x)$$ exists.

The answer appear to be not at all trivial. Based on numerical computations, I conjecture that the limit exists if and only if $$x\in[-\sqrt{2},1+\sqrt{2}]-\mathbb{N}.$$ But for some $x$, the convergence is very strange. For example, look at this plot of the sequence $\{g_{3n}(0.18)\}$ for $0\leq n\leq 200$:

EDIT: As pointed out by Zach Effman, there seem to be a fractal structure emerging from this problem. To better visualize it, I looked at the boundedness of the sequence $\{g_n(z)\}$ for $z\in\mathbb{C}$. More precisely, I made a grid of points $z\in\mathbb{C}$, colored black if $\{|g_n(z)|\}$ is bounded at $z$ and the more quickly $\{|g_n(z)|\}$ diverges to infinity the more white the pixel is. A dramatic fractal structure emerges: Does anybody recognize this fractal? Is this a known fractal?

The center of this picture reveals a very nice structure:

As noted by mercio, your function $f(x)=x^2-x-1$ has a three-cycle $$2 \cos \left(\frac{\pi }{7}\right) \mapsto 2 \sin\left(\frac{\pi }{14}\right) \mapsto -2 \sin\left(\frac{3 \pi }{14}\right) \mapsto 2 \cos \left(\frac{\pi }{7}\right).$$ Put another way, these points are fixed points of $F(x)=f(f(f(x)))$. Let us denote these points by $x_F^1$, $x_F^2$, and $x_F^3$. The cycle is neutral because $F'(x_F^i)=1$ for each $i$.

Now, suppose that $x_0>1+\sqrt{2}$. Then, \begin{align} f(x_0) &= x_0^2 - x_0 - 1 \\ &= \left(1+\sqrt{2}\right)+\left(x_0 - (\sqrt{2}+1)\right) \left(x_0+\sqrt{2}\right). \end{align} As a result, the orbit of $x_0$ diverges to $\infty$. Furthermore, if $x_0<-\sqrt{2}$, then $x_1=f(x_0)>1+\sqrt{2}$ so (again) the orbit diverges to $\infty$.

Thus, let $I=[-\sqrt{2},1+\sqrt{2}]$. We certainly know that every initial seed outside of $I$ diverges to $\infty$. As it turns out, almost every point inside $I$ converges to the fixed neutral orbit. The exceptional set of points forms a Cantor set $C$ in $I$ and the dynamics of $f$ are chaotic on $C$. It's not particularly easy to see this, however, without some fairly high powered (but famous) theorems.

To get a better grip on this, let's plot the function together with the line $y=x$.

We see quite clearly that the function maps $[-\sqrt{2},1+\sqrt{2}]$ into itself and there are a couple of fixed points which cross the line $y=x$ with slope larger than one in absolute value so that the will be repulsive. Since, we have an orbit of period 3, let's include a plot of $F(x)=f(f(f(x)))$.

Note that there are 3 points where the line $y=x$ is tangent to the graph of $F(x)$. These are exactly the points in the neutral orbit of period 3. Each of these is the left hand endpoint of an interval of points that converge to that fixed point under iteration of $F$. I've highlighted one such interval containing $1/2$ in red. Furthermore, the inverse image of this interval under $F$ consists of 6 more such intervals that map into the red interval. The inverses of these yield more points that converge to the middle fixed point under iteration of $F$, etc. The same can be done for the other fixed points. The complement of all these open intervals containing points converging to one of those fixed points is exactly the Cantor set of points that do not converge to any of those fixed points.

As mercio and Sheldon have both pointed out, the groovy fractal image is exactly the Julia set of $z^2-z-1$. To emphasize the relationship between this fact and my answer, let's look at the Julia set together with the red interval from my last image.

Thus, those intervals span the components of the Julia set with simply connected interior. The yellow point is your initial seed $x_0=0.18$. As it is not contained in any of the larger components of the Julia set, I would certainly expect it to jump around before it eventually settles down to an orbit.

It might also be worth pointing out that the true fractal Julia set is the boundary of the black region - thus, it looks something like so:

And here's a zoom onto the main component.

Finally, we have here a period $3$ orbit of the real map $f$ and there is a very famous paper entitled "Period three implies chaos". The theorems of that paper are exactly what one needs to definitely prove that $f$ is chaotic on a Cantor set $C$. I believe that $C$ has measure zero; it is certainly nowhere dense. That explains why it is hard to find point in $C$ by experiment. A randomly chosen point in $I$ generically and (I think) almost surely converges to the neutral period 3 orbit.

On the other hand, it is (in principle) feasible to find orbits of large prime period. Here's Mathemetica code that finds an orbit of period 41.

Clear[f];
f[x_] := x^2 - x - 1;
f[n_Integer, x_?NumericQ] := Nest[f, x, n];
per41 = x /. FindRoot[f[41, x] == x, {x, 0.1},
  WorkingPrecision -> 300]

(* Out: 
  0.09997111618903433899874853864417001759319794189148766236421781597929\
  2721939448495289732175832121442601846472892582352497849991420790365107\
  4378578596438323899498328663812648206777989645187097865648459105422784\
  2754855914820032740936983665825718876082741645285990326028932997157116\
  60320676920709780319319
*)

Note that high precision is necessary, since there will be points arbitraryily near by with other periods. Here's a sanity check:

Nest[f, per41, 41] - per41
(* Out: 0.*10^-283 *)

Here are some partial results. I show below that there are only three possible limits (when the sequence is nonconstant), and for each limit I find an interval such that if the initial value is in that interval then it converges to that limit.

Applying Zach Effman’s method to $g=f^3$ rather than to $f$, we have

$$ g(x)-x=(x^2-2x-1)(x^3-x^2-2x+1)^2 \tag{1} $$

The roots of $x^2-2x-1$ are

$$ \alpha_1=1-\sqrt{2},\ \ \alpha_2=1+\sqrt{2} \tag{2} $$

While the roots of $x^3-x^2-2x+1$ are

$$ \beta_1 \approx -1.246,\ \ \beta_2 \approx 0.445, \ \ \beta_3 \approx 1.801 \tag{3} $$

We have

$$ g'(\alpha_1)=25-22\sqrt{2}\approx -6.11, \ \ g'(\alpha_2)=25+22\sqrt{2}\approx 56.11 \tag{4} $$

So both those fixed points are repulsive. On the other hand,

$$ g'(\beta_1)=g'(\beta_2)=g'(\beta_3)=1 $$

(because of the identity $g'(x)-1=(x^3-x^2-2x+1)(8x^4 - 20x^3 - 4x^2 + 18x + 2)$)

There is a unique $\gamma_1<\beta_1$ such that $g(\gamma_1)=\beta_1$. It is easy to see then that the interval $I_1=[\gamma_1,\beta_1]$ is stable by $g$, and if $x\in I_1$ then the sequence converges to $\beta_1$.

Similarly : there is a unique $\gamma_2\in (\beta_2,\beta_3)$ such that $g(\gamma_2)=\beta_2$, and then the interval $I_2=[\beta_2,\gamma_2]$ is stable by $g$, and if $x\in I_2$ then the sequence converges to $\beta_2$, and there is a $\gamma_3>\beta_3$ such that the interval $I_3=[\beta_3,\gamma_3]$ is stable by $g$, and if $x\in I_3$ then the sequence converges to $\beta_3$.

This answer is inspired by comments with Mercio. First we map the problem to an equivalent $z \mapsto z^2+c\;$ Mandelbrot problem. And then we can use all of the known tools for analyzing the Mandelbrot set to analyze the Op's problem. For example, see Wolf Jung's mandelbrot link, which includes a Mandelbrot program with landing ray's.

$$z \mapsto z^2-\frac{7}{4}$$

This is a conjugation of the original problem with z=x-1/2, and $g(x)=f(x)-1/2$ For clarity, call the new function $g(z)=z^2-7/4$. We calculate the polynomial roots of $g^{o3}(z)-z=0$, which factors into a product of two identical cubic equations, and a quadratic. Here we show the parabolic 3-cyclic attracting fixed points, which are the roots of the cubic equation, followed by the roots of the 1-cyclic quadratic equation.

$$8z^3 + 4z^2 - 18z - 1=0\; \left(r_1= -1.7469796037; r_2=-0.054958132087; r_3= 1.3019377358\right)$$ $$4z^2 - 4z - 7\; z=\left(r_4= -0.914213562373, r_5=1.91421356237 \right)$$

If we iterate $z \mapsto g(z)$ an infinite number of times, every point that remains bounded is in the Julia set for $c=-7/4$. Points that are completely inside a bug will eventually you get to a 3-cyclic attracting period 3-cycling towards the three parabolic fixed points mentioned above. Otherwise, if there are point that escape to $\infty$ arbitrarily close to your point, then the point you picked is on the boundary of the Julia set. And then we need to figure out the rotation angle for the landing ray which lands at your point, perhaps using the Bottcher function.

But the simplest way to handle points exactly on the Julia set boundary is to look at values of $g^{on}(z)$ for finite values of n. If $g^{on}(z)=r_1,r_2,r_3, r_4, r_5$ for some finite value of n, then the 3-cycle converges, since once you land on one of these five points, iterating $g^{o3}(z)$ does not change the point. Other than these five cases, the function $g^{o3}(z)$ will not eventually repeat, nor will it attract towards a repeating point. The three cycle parabolic attracting points are only attracting for points in a connected region; in other words inside a Julia set bug, and not on the boundary.

One can also use the Bottcher function to put iterating $z \mapsto z^2$ into correspondence with iterating $z \mapsto z^2-7/4$, as both functions iterate to infinity. There will be an analytic Bottcher function mapping $\infty$ to zero. And then we generate the Bottcher function for $g(z)$

$$\beta\left(\frac{1}{g(1/z)}\right) = \beta(z)^2$$

The Bottcher function maps the unit circle to the boundary of the Julia set for $c=-7/4$. So if we have a point on the boundary of the Julia set, then the point corresponds to a particular number in the Boettcher set of the form $$\exp^{2\pi i r}$$ Where r is the rotation angle of the point z. We are particularly interested in these rotation angles: $$r = \left( \frac{\pm 1}{7}, \frac{\pm 2}{7}, \frac{\pm 4}{7}, \frac{\pm 1}{3}, 1 \right)$$

Since landing rays with these three rotation angles correspond to the three parabolic fixed points mentioned above. r=1 corresponds to the solution z=1.91421356237, and r=$\pm 1/3$ corresponds to -0.914213562373. Notice that iterating $z \mapsto z^2$ will double the rotation angle r. And if we start with r=1/7, then iterating $z^2-7/4$ will go through a 3-cycle iterating (1/7,2/7,4/7. There are two landing rays for each of the parabolic cusps; so (-1/7,-2/7,-4/7) is identical. In addition, we have the rotation angles for $r=(\pm 1/3, 1)$, which correspond to the two solutions for the primary fixed point.

If we look at the rational number for the rotation angle, $\frac{p}{q}$, for the point z, then we can tell how it will repeat. In particular, if $q=7\cdot 2^n$ or $q=3\cdot2^n$, or $q=2^n$ then the point is a pre-periodic point, and after the appropriate number of iterations, $g^{o3}(n)$ will repeat exactly. All other rotation angles will upon iterating $z \mapsto z^2 - 7/4$ will not go to one of these five points.

Not a complete answer, but this is what I can determine right away:

What you're describing is an iterated function. If the sequence $f^n(x)$ converges, then we must have $f(L) = L$ for the limit point. This gives $x^2-x-1=x$ with roots $1\pm\sqrt{2}$. At a fixed point, the absolute value of the derivative of our function, $|f'(x)|$ determines the local behavior of the iterated sequence. If the absolute value of the derivative is less than one, the fixed point is an attractor and starting points nearby converge to it. For the function here, $f'(x) = 2x-1$ giving $|f'(1-\sqrt{2})| = 2\sqrt{2}-1 \approx 1.8$ and $|f'(1+\sqrt{2}) = 1+2\sqrt{2} \approx 3.8$. Unfortunately, both are greater than one. This means the fixed points are repellors: starting points arbitrarily close to them are pushed away.

All I can say, then, is that the iterated function definitely converges for starting points $x_0 = 1\pm \sqrt{2}$ because in these cases the sequence is constant. Actually, that's not quite true: any sequence that hits those points will stay constant. So we could ask, what are the points $x$ such that $f(x) = 1 \pm \sqrt{2}$? They're $1\pm \sqrt{2}$ (of course) as well as $\pm\sqrt{2}$. Repeating this, these are the points I get:

I imagine there's some sort of fractal structure here.

I only attempted computational results as of now.

I first noticed one $x$ value out of your interval that gave a convergent sequence, which lead me to numerically test myself. Here is that value:

Let $\alpha$ be the root of $x^3 - x^2 - 2x + 1$ near $-1.24698$. Then for all $n$, $g_{3n}(\alpha) = \beta$, where $\beta$ is the root of $x^3 - x^2 - 2x + 1$ near $1.80194$.

From my own testing, I conjecture that the limit exists if and only if

$$ x \in [-\sqrt{2}, \; 1+\sqrt{2}] - S, $$

where $S$ is the set of points that form a cycle, i.e.

$$ S = \{x \in \mathbb{R} \mid g_{3k}(x) = g_0(x) \text{ for some } k, \; g_3(x) \neq g_0(x)\}. $$

The members of $S$ that correspond to $k = 2$ are approximately

$$ \{-1.118, -1.09847, -1, -0.962622, -0.677575, -0.601782, -0.496701, -0.367928, -0.30512, -0.256588, -0.0360765, 0, 0.110737, 0.136684, 0.863316, 0.889263, 1, 1.03608, 1.25659, 1.30512, 1.36793, 1.4967, 1.60178, 1.67758, 1.96262, 2, 2.09847, 2.118\}. $$

My methodology was for a fixed $x$, make a table of the first 2000 values of $g_{3n}$, find the power function of best fit, and assume the sequence converges if and only if the exponent of the fitted function is negative.

If you're curious, here is the Mathematica code I used:

f[x_] := x^2 - x - 1
g3n[n_, x_] := g3n[n, x] = f[f[f[g3n[n - 1, x]]]]
g3n[0, x_] := f[x]

SequenceTable[z_] := Table[
  If[# === Overflow[], Throw[False], #]&[g3n[n, z]], 
  {n, 0, 2000}
]

NegativeExponentQ[z_?NumericQ] := Catch@Quiet[
  Negative[b /. {
      FindFit[N@Differences@SequenceTable[z], a x^b, {a, b}, x]]
    }, 
  {General::ovfl, FindFit::sszero}
]

BinarySearch[l_, h_, {err_, dir_}] := Block[{lo = l, mid, hi = h, GoodResultQ},
  GoodResultQ = If[dir < 0, Not, Identity];
  mid = (lo + hi)/2;
  While[hi - lo > err,
    If[GoodResultQ[NegativeExponentQ[mid]],
      lo = mid,
      hi = mid
    ];
    mid = (lo + hi)/2;
  ];
  mid
]

And here are the commands I ran and their outputs:

Monitor[
  BinarySearch[-2.1`20, 0.1`20, {10^-10., -1}],
  {lo, mid, hi}
]

-1.4142134775087470189

Monitor[
  BinarySearch[0.1`20, 2.5`20, {10^-10., 1}],
  {lo, mid, hi}
]

2.4142135623493231833