14

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?

Mohsen Shahriari
  • 6,718
templatetypedef
  • 10,206
  • You do realize that if your sequence has a limit $L$, then $1/L = 1/L + 1/L$? – JonathanZ Sep 12 '21 at 23:09
  • 1
    @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. – templatetypedef Sep 12 '21 at 23:15
  • 2
    Looks like a case for the Banach Fixed Point Theorem – Hyperplane Sep 12 '21 at 23:56
  • Have you checked cases where $a_0<0$ and $a_1>0?$ The sequence isn’t even defined for $a_0=-a_1,$ for example. – Thomas Andrews Sep 13 '21 at 00:11
  • @ThomasAndrews Yep - they converge to $-\sqrt{2}$, or at least the ones I checked do. (And yes, you’re right, the sequence isn’t defined for that case.) – templatetypedef Sep 13 '21 at 00:25
  • A limit exist if and only if the sequence is cauchy. You can try to see if if $m < n$ that $|a_n - a_m| \to 0$ as $m\to \infty$ so the sequence is cauchy so a limit must exist. And maybe if you consider $a_n = \sqrt 2 - d_n$ for some sequence of $d_n$ and we cann use $\frac 1{\sqrt 2 - d_n}> \frac 1{\sqrt 2} + \frac 1{d_n}$ if $d_n > 0$ etc.. to somehow bind the values of $|a_n - a_m|$. That's my idea for brainstorming. I have no idea if it will pan out. – fleablood Sep 13 '21 at 00:27
  • 2
    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. – fleablood Sep 13 '21 at 00:31
  • @fleablood Great call on the golden ratio - that gives some values that don’t stabilize at $\sqrt{2}$. Question updated accordingly. – templatetypedef Sep 13 '21 at 00:42
  • It was a wild guess. I imagine you can find some values will oscilate to values converging to $0$ and values diverging to infinity. – fleablood Sep 13 '21 at 00:48
  • I haven't thought this through for the case when the sign oscillates every term, but if there are 2 consecutive terms with the same sign, then the sequence is bounded. If fact, for any $a\ge 1$, if $a_n,a_{n+1}\in\left[\frac{a^{-1}}{\sqrt{2}},\frac{a}{\sqrt{2}}\right]$, then $a_{n+2}\in\left[\frac{a^{-1}}{\sqrt{2}},\frac{a}{\sqrt{2}}\right]$ as well. – Alexander Burstein Sep 13 '21 at 03:20
  • 2
    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}$. – JimmyK4542 Sep 13 '21 at 03:30
  • 1
    Related: Is the sequence defined by the recurrence $a_{n+2}=\frac{1}{a_{n+1}}+\frac{1}{a_n}$ convergent? and Proof of existence of a limit for the sequence recursively-defined with $a_1=1$, $a_2=1$ and $a_n=\frac{1}{a_{n-1}}+\frac{1}{a_{n-2}}$ for $n\ge2$. – Mohsen Shahriari Oct 02 '21 at 01:15

1 Answers1

16

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.

enter image description here

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.

enter image description here

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.

enter image description here

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])
Hyperplane
  • 11,659