# 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**

### # Image credits

- Hero: How I made money using The Golden Ratio Law of Design this Pandemic
- Problem statement: Otto K. Bretscher
- Fibonacci sequence: Fibonacci sequence - Wikipedia
- Fibonacci spiral: Wikimedia Commons, Fibonacci Spiral over tiled squares
- Rotation and Scaling Effects: Geogebra
- Characteristic Polynomial from Wolfram Alpha: Wolfram Alpha
- PARI Solution: PARI/GP
- Recursive image of this problem in eurAIka: 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)

## # Further reading

- Interactive Linear Algebra, Dan Margalit, Joseph Rabinoff. 5.5 Complex Eigenvalues
- Diagonalize a 2 by 2 Matrix $A$ and Calculate the Power $A^{100}$
- COMPLEX Eigenvalues, Eigenvectors & Diagonalization
**full example**- Dr. Trefor Bazett - Fibonacci sequence - Wikipedia
- Matrix multiplication - Wikipedia

## # 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.

### 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.

### 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 →