# Recursive Sequences as Linear Transformations

Subtitle: Fibonacci Spirals the Drain

This problem is from Otto K. Bretscher posted in the Classical Mathematics Facebook group.

We’re working on a new platform for scientific research called eurAIka that combines The Whiteboard for writing notes, equations, and pasting in graphics, The Librarian for AI-assisted literature search, and The Coder where you can write programs. This will be our first attempt at using eurAIka for a Wild Peaches article.

# # The Fibonacci Sequence

Let’s start with an easier problem, the Fibonacci sequence:

The Fibonacci Sequence

where the recursive relationship is $x_{n+2} = x_{n+1} + x_n$ and $x_0 = x_1 = 1$​. We could generalize this to include any sequence composed of multiples of the previous two terms,

$x_{n+2} = \alpha x_{n+1} + \beta x_{n}.$

The coefficients $\alpha = \beta = 1$ for the Fibonacci sequence, and for the problem we’re trying to solve $\alpha = 4$ and $\beta = -7$.

Fibonacci Spiral

It wouldn’t be difficult to write a function in a computer language that calculates all of the elements of the sequence up to the one we’re interested in. To see how the value for $x_1 = b$ changes the outcome, we’d have to run the code over many test cases to see what happened.

Another way to get to the value of $x_{2024}$ can be done in one step using linear algebra. Let’s call $\overrightarrow{v_{0}}$ the vector composed of the first two entries of the sequence, so

$\overrightarrow{v_{0}} = \begin{bmatrix} x_{1} \\ x_{0} \end{bmatrix}.$

The next step in the sequence puts $x_2$ in the first position and moves $x_1$ down giving

$\overrightarrow{v_{1}} = \begin{bmatrix} x_{2} \\ x_{1} \end{bmatrix} = \begin{bmatrix} \alpha x_{1} + \beta{x_{0}} \\ x_{1} \end{bmatrix}.$

The entries of this new vector are linear combinations of the previous two entries of the sequence, so we could define a matrix $A$ that multiplies $\overrightarrow{v_{0}}$ to get

$\overrightarrow{v_{1}} = A \overrightarrow{v_{0}} = \begin{bmatrix} \alpha & \beta \\ 1 & 0 \end{bmatrix} \begin{bmatrix} x_{1} \\ x_{0} \end{bmatrix}.$

The recursion step says that $\overrightarrow{v_2} = A \overrightarrow{v_1}$ and so on, or in terms of the original vector, $\overrightarrow{v_2} = A(A \overrightarrow{v_0}) = A^2 \overrightarrow{v_0}.$ For the Fibonacci sequence, $A = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}$. For example, we can calculate the 16th entry in the sequence by calculating $A^{14} = \begin{bmatrix} 610 & 377 \\ 377 & 233 \end{bmatrix}$, which gives $F_{15}$ in the $[1,1]$ location. A convenient trick to extract just the upper left entry using Julia is to write this as (A^14)[1,1].

If $Ax = \lambda x$ then $x$ is called an eigenvector and the constant $\lambda$ is an eigenvalue. For the Fibonacci matrix $A$ there are two eigenvalues, $\lambda = \frac{1 \pm \sqrt{5}}{2}$ (the positive one is the golden ratio, $\phi$) and if you take the ratio of $F_{n+1}$ to $F_n$ you’ll see that it converges to $\lambda$. The other eigenvalue is the negative inverse of the first one, so $\lambda_2 = -F_n / F_{n+1}$.

The matrix $A$ can be decomposed into $A = P D P^T$ where $P$ is the $2 \times 2$ array of eigenvectors and $D$ is a $2 \times 2$ diagonal array with $\lambda_1$ and $\lambda_2$ on the diagonal. The eigenvectors are unit vectors, so $P$ is orthogonal meaning that $P P^T = I$, the identity matrix. With this decomposition, powers of $A$ can be written as

$A^n = P D^n P^T.$

Since $D$ is diagonal, then $D^n$ is also diagonal with $\lambda_1^n$ and $\lambda_2^n$ on the diagonal, making matrix powers easy to calculate. Now the 15th iterate of $A$ can be found from

$A^{14} = P D^{14} P^T = \begin{bmatrix} 610.0 & 377.0 \\ 377.0 & 233.0 \end{bmatrix}.$

This is numerically unstable because $D$ will quickly overflow or underflow as powers of $\lambda_1$ and $\lambda_2$ approach either zero or $\pm \infty,$ but this idea will become theoretically useful for our problem.

As we saw above, the ratios of successive values of the Fibonacci sequence approach the golden ratio, $\phi$. If you plotted points $(F_n, F_{n+1})$ you would see that they all lie nearly on the line $y = \phi x$, and the fit gets better the farther out you go. The points $(F_n, F_{n+1})$ are the vectors which when multiplied by $A$ give the next point or vector in the sequence.

Geometrically, matrix multiplication transforms vectors through rotation and scaling. What we’re seeing here is that in the decomposition $P D^n P^T v_{n}$, the first multiplication $P^T v_n$ only performs a rotation since $P$ is unitary. The scaling occurs with $D^n (P^T v_n)$, and multiplication by $P$ restores the direction.

As an example, let $v_5 = [5, 3]^T$ so

\begin{aligned} P^T v_{5} &= \begin{bmatrix} 0.0767 \\ -5.830 \end{bmatrix} \\ DP^T v_{5} &= \begin{bmatrix} -0.0474 \\ -9.434 \end{bmatrix} \\ PDP^T v_{5} &= \begin{bmatrix} 8 \\ 5 \end{bmatrix} \\ \end{aligned}

Rotation and Scaling Effects

Notice that $v_5$ and $P^T v_5$ lie on the same circle as do the points $DP^Tv_5$ and $PDP^Tv_5$ showing that $P$ is a rotation, while the scaling between $F_5$ and $F_6$ is multiplication by $\lambda_1$​.

# # The Bretscher Sequence

Now let’s apply this same technique to the Bretscher sequence $x_{n+2} = 4 x_{n+1} - 7 x_n$. Define the transition matrix,

$B = \begin{bmatrix} 4 & -7 \\ 1 & 0 \end{bmatrix}$

which has the characteristic equation derived from the determinant of $|B - \lambda I|$

\begin{aligned} 0 &= (4 - \lambda)(-\lambda) + 7 \\ &= \lambda^2 - 4 \lambda + 7 \\ \Rightarrow \lambda &= 2 \pm \sqrt{ 3 } i. \end{aligned}

The eigenvectors are found by solving for $(B - \lambda I) u = 0$ for some $u$. This gives

$\begin{bmatrix} 4 - (2 + \sqrt{3}i) & -7 \\ 1 & -(2 + \sqrt{3}i) \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \end{bmatrix} = 0.$

Expanding this gives the system of equations,

\begin{aligned} (2 - \sqrt{3}i)u_1 - 7 u_2 &= 0 \\ u_1 - (2 + \sqrt{3}i) u_2 &= 0 \end{aligned}

giving eigenvectors $v = \begin{bmatrix} 2 \pm \sqrt{3}i & 1 \end{bmatrix}$. The inverse of

$P = \begin{bmatrix} 2+\sqrt{3}i & 2-\sqrt{3}i \\ 1 & 1 \end{bmatrix}$

is

$P^{-1} = \begin{bmatrix} -\frac{i}{2 \sqrt{3}} & \frac{1}{2} + \frac{i}{\sqrt{3}} \\ \frac{i}{2 \sqrt{3}} & \frac{1}{2} - \frac{i}{\sqrt{3}} \end{bmatrix}.$

Let $D$ be the diagonal matrix formed from the eigenvalues,

$D = \begin{bmatrix} 2+\sqrt{3}i & 0 \\ 0 & 2-\sqrt{3}i \end{bmatrix}$

and then verify that $B = PDP^{-1}$ and $P P^{-1} = I$, the identity matrix.

For the Fibonacci sequence, the eigenvectors formed an orthonormal basis, so we could quickly calculate powers of the matrix $A$ because

$A^n = PDP^{-1} PDP^{-1} \cdots PDP^{-1} = PDI \cdots I D P^{-1} = P D^nP^{-1}.$

The Bretscher matrix also has this property, except that the eigenvectors aren’t orthogonal. Still, the eigenvectors form a basis of the space spanned by $B$, so we can still calculate $B^n = PD^nP^{-1}$ for any power $n$. This wouldn’t be so bad for some small values of $n$, but we’ll need another method to answer the question of how to maximize $x_{2024}$.

Otto Bretscher says in response to his original question,

The roots of the characteristic polynomial are 2 ± i√3 = √7exp[±i*arctan(√3/2)]. We can write x(n) as a linear combination of (2 + i√3)^n and (2 - i√3)^n and bring this expression into real form. While these computations are straightforward, they are somewhat tedious, perhaps best left to a computing device; the answer provided by WolframAlpha is attached. Since the coefficient of b turns out to be negative for n = 2024, we maximize x(2024) by letting b = 0.

Characteristic Polynomial from Wolfram Alpha

Pari/GP can calculate numbers with very high precision, so let’s give it a try. The sequence begins $x_0 = 1, x_1 = b$, so the next element is found from

$B \begin{bmatrix} x_{1} \\ x_{0} \end{bmatrix} = \begin{bmatrix} 4 & -7 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} b \\ 1 \end{bmatrix} = \begin{bmatrix} 4b - 7 \\ b \end{bmatrix}.$

The fourth element is $B^2 \begin{bmatrix} b \\ 1 \end{bmatrix}$ so the $2024^{th}$ will be multiplied by $B^{2022}$. Enter the matrix $B$ in PARI/GP as

B = [4, -7; 1, 0]

PARI Solution

So, how do you choose $b \in [0,2 \pi]$ to maximize $x_{2024}$? Well you don’t need to know the value of the number that Pari just computed, all you need is the sign. Since the value of the entry in the $[1,1]$ position of this little $2 \times 2$ matrix is negative, the $b$ that maximizes $x_{2024}$​ is zero.

Recursive image of this problem in eurAIka

## # Code for this article

Problem2024.ipynb - JupyterLab notebook in Julia

Bretscher_sequence.ipynb - JupyterLab notebook using Wolfram Language 13.2

Fibonacci.html - Geogebra worksheet to demonstrate rotation and scaling effects of linear transformations (open in br)

## # Software

### # Julia

The Julia Project as a whole is about bringing usable, scalable technical computing to a greater audience: allowing scientists and researchers to use computation more rapidly and effectively; letting businesses do harder and more interesting analyses more easily and cheaply.

Posts using Julia

### # JupyterLab

JupyterLab is the latest web-based interactive development environment for notebooks, code, and data. Its flexible interface allows users to configure and arrange workflows in data science, scientific computing, computational journalism, and machine learning.

Posts using JupyterLab

### # Wolfram Language

Wolfram Language is a symbolic language, deliberately designed with the breadth and unity needed to develop powerful programs quickly. By integrating high-level forms—like Image, GeoPolygon or Molecule—along with advanced superfunctions—such as ImageIdentify or ApplyReaction—Wolfram Language makes it possible to quickly express complex ideas in computational form.

Posts using Wolfram Language

### # Pari/GP

PARI/GP is a widely used computer algebra system designed for fast computations in number theory (factorizations, algebraic number theory, elliptic curves, modular forms, L functions…), but also contains a large number of other useful functions to compute with mathematical entities such as matrices, polynomials, power series, algebraic numbers etc., and a lot of transcendental functions.

Posts using Pari/GP

See all software used on wildpeaches →