# Why does the process defined with $a_{n+2} = \frac{1}{a_n} + \frac{1}{a_{n+1}}$ converge to $\pm\sqrt{2}$ for most choices of the starting values?

I was tinkering around with Fibonacci-style recurrences and, on a whim, tried this one:

$$a_0 = 1 \qquad a_1 = 1 \qquad a_{n+2} = \frac{1}{a_n} + \frac{1}{a_{n+1}} \text.$$

I had no idea what this would look like, so I whipped up a short Python script:

a = 1
b = 1
for i in range(100):
print(a)
a, b = b, 1/a + 1/b


I was surprised to see that this sequence quickly converged on $$\sqrt{2}$$. Wondering if that was an artifact of the starting conditions, I tried changing $$a_0$$ and $$a_1$$ to other values to see what happened. Assuming $$a_0 \ne -a_1$$ and neither value was zero, I found that the numbers always seemed to converge on either $$\sqrt{2}$$ or $$-\sqrt{2}$$, depending on the signs of the initial values.

It makes sense to me why once the two most recent values are close to the square root of two, the next term will also be close to the square root of two:

$$\frac{1}{\sqrt{2}} + \frac{1}{\sqrt{2}} = \frac{2}{\sqrt{2}} = \sqrt{2}.$$

However, I’m not sure why values keep getting pulled in toward $$\pm\sqrt{2}$$. I was surprised that the values didn’t, say, oscillate out toward infinity.

I have a suspicion that this recurrence might be a hidden application of something like the Babylonian method for finding square roots, but my background on these sorts of things is somewhat limited and I’m not seeing a quick way to map this recurrence onto it.

Is there a straightforward explanation as to why this recurrence tends towards $$\pm\sqrt{2}$$ for just about all initial values?

• @JonathanZ Yes, that makes sense to me. It would explain why, if we get values close to $\sqrt{2}$, then future values will stay there as well. But I’m not sure why we get values close to $\sqrt{2}$ in the first place, especially for arbitrary starting values. Sep 12, 2021 at 23:15
• Looks like a case for the Banach Fixed Point Theorem Sep 12, 2021 at 23:56
• Alternatively, maybe there are values where the sequence will oscillate but chaotically they are very fair and few between and the system is "unstable" that these oscillating initial values will be very complicated to calculate and find (It's plausible that they must be weird irrational values) but do exist. For gits and shiggles have you tried anything using the golden ratio.... (doubt it will work)... I think it is safe to say this may be a hard problem. Sep 13, 2021 at 0:31
• If $a_0 = 1$ and $a_1 = \tfrac{1+\sqrt{5}}{2}$, then $a_2 = \tfrac{1+\sqrt{5}}{2}$, and $a_3 = \sqrt{5}-1$, so $a_n$ does not toggle back and forth between $1$ and the golden ratio. In particular, the sequence $\{a_n\}$ cannot alternate between two distinct values because if $a_{n+3} = a_{n+1}$ and $a_{n+2} = a_n$, then $a_{n+3} = \tfrac{1}{a_{n+1}}+\tfrac{1}{a_{n+2}} = \tfrac{1}{a_{n+1}}+\tfrac{1}{a_n} = a_{n+2}$. Sep 13, 2021 at 3:30

EDIT: I uploaded a higher-resolution, correctly flipped image for convenience. Upon further inspection, it turns out there is a hidden small-scale fractal structure. An uncompressed 10k zoom on the interesting part can be found here.

EDIT: This is what the basins of attraction look like. It seems that indeed there is always converges against $$x=y=±\sqrt{2}$$, except for the lines $$x=-y$$ and $$y=0$$ and a complicated bow-tie figure. Python source code, using JAX at the bottom.

The line $$x=0$$ is a well-behaved singularity, because

$$\begin{pmatrix}±0\\y\end{pmatrix} ⟼ \begin{pmatrix}y \\ ±∞ \end{pmatrix} ⟼ \begin{pmatrix}±∞\\ 1/y\end{pmatrix} ⟼ \begin{pmatrix}1/y\\ y\end{pmatrix}$$

The iteration can be encoded in the function $$f(x, y) = (y, \tfrac{1}{x} + \tfrac{1}{y})$$. Note that

$$\|f(x,y)\|_2^2 ≤ \|(x,y)\|_2^2 ⟺ y^2 + (\tfrac{1}{x} + \tfrac{1}{y})^2 ≤ x^2 + y^2 ⟺ (\tfrac{1}{x} + \tfrac{1}{y})^2 ≤ x^2$$

Within regions where this equality is strict, the function is contractive. However, both a plot and analysis reveal that the critical points $$x=y=±\sqrt{2}$$ lie exactly on the boundary.

However, if we move to the second iteration we have:

$$\|f^{(2)}(x,y)\|_2^2 ≤ \|(x,y)\|_2^2 ⟺ \Big(\tfrac{1}{\frac{1}{x} + \frac{1}{y}}\Big)^2 + \Big(\tfrac{1}{y} + \tfrac{1}{\frac{1}{x} + \frac{1}{y}}\Big)^2 ≤ x^2 + y^2$$

Now, we see the critical point lie well inside the region and thus we have convergence by Banach's Fixed Point Theorem.

There are still some technical holes to fill (*), but I think this should give a rough idea of why there is convergence.

Hypothesis/Conjecture: Let $$C_n = \{(x, y) ∣ \|f^{∘n}(x,y)\|_2^2 < \|(x,y)\|_2^2\}$$ be the region for which $$f^{∘n}=f∘f\circ\ldots∘f$$ is contractive. Then, in the set-theoretic limit

$$\lim\limits_{n\to ∞} C_n = ℝ^2\setminus \{(x,y) \mid x=-y ∨ x=0 ∨ y=0\}$$

(*) Basically, we can formulate an argument like this: We already know $$x=y=\sqrt{2}$$ is a fixed point. By (Lipschitz-) continuity, $$f$$ maps a ball around the fixed point to a nbhd. of the fixed point. Choose the original ball small enough such that it is contained in the contractive region. Then its image must be a subset of itself since $$\|f(x)-x^*\| =\|f(x)-f(x^*)\|≤\|x-x^*\|$$, hence $$f$$ restricted to this ball satisfies all technical necessary conditions for applying Banach's Fixed Point Theorem.

# Source Code

import jax
import jax.numpy as np
from tqdm import trange
import matplotlib.pyplot as plt
jax.config.update("jax_enable_x64", True)

@jax.jit
def f(a):
x, y = a[0], a[1]
return np.array([y, 1/x + 1/y])

N = 2**12  # 4k resolution
xmin, xmax = -.5, +.5
ymin, ymax = -.5, +.5
X = np.linspace(xmin, xmax, num=N)
Y = np.linspace(ymin, ymax, num=N)
Z = np.stack(np.meshgrid(X, Y))

for k in trange(1000):
Z = f(Z)

fig, axes = plt.subplots(ncols=2, figsize=(12, 6), tight_layout=True,
sharex=True, sharey=True, subplot_kw=dict(box_aspect=1))

for ax, im, title in zip(axes, Z, ("x", "y")):
img = ax.imshow(im, interpolation="none", vmin=-1.5, vmax=1.5)
converged = np.mean(np.isclose(im, np.sqrt(2)) | np.isclose(im, -np.sqrt(2)))
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.set_xticks(np.linspace(0, N, 9))
ax.set_yticks(np.linspace(0, N, 9))
ax.set_xticklabels(np.linspace(xmin, xmax, 9))
ax.set_yticklabels(np.linspace(ymin, ymax, 9))
ax.set_title(title+f" {converged=:.2%}")

fig.savefig(f"4248805_{N}px.png", dpi=N/fig.get_size_inches()[1])