Recursive Sequences as Linear Transformations

Fibonacci Spirals the Drain

February 02, 2024 · John Peach

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:

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,

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

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

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

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

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

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

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

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,

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

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

Expanding this gives the system of equations,

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

is

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

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

The Bretscher matrix also has this property, except that the eigenvectors aren’t orthogonal. Still, the eigenvectors found above 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.

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

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]

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.

Spiraling to the Solution: Where Mathematics Meets Beauty

This exploration of the Bretscher sequence demonstrates how linear algebra can be useful in solving seemingly complex recursive problems. We found that

• Recursive sequences can be reframed as linear transformations using matrices.
• The eigenvalues and eigenvectors of these matrices reveal the long-term behavior of the sequence.
• Complex eigenvalues, as in the Bretscher sequence, create spiral patterns in the sequence behavior.
• Matrix decomposition $(PDP^{-1})$ provides an elegant way to compute high-powered terms .

If you’re interested in exploring further, here are some suggested investigations:

• Experiment with different coefficients $(\alpha, \beta)$ in the general sequence $x_{n+2} = \alpha x_{n+1} + \beta x_n$ and observe how they affect the sequence behavior.
• Create visualizations of how different recursive sequences spiral in the complex plane.
• Investigate what happens when you extend the recursion to depend on three or more previous terms.
• Study how the initial conditions $(x_0,x_1)$ affect the long-term behavior of various recursive sequences.
• Explore connections between recursive sequences and natural phenomena like population growth or wave propagation.

The methods demonstrated here extend beyond just mathematical curiosity - they have applications in fields ranging from financial modeling to quantum mechanics. Understanding how to analyze recursive sequences through linear transformations provides a powerful tool for tackling complex problems in science and engineering.

Definitions

Linear combination: If vector $v_3 = \alpha v_1 + \beta v_2$ for some constants $\alpha$ and $\beta$, then $v_3$ is a linear combination of the two vectors $v_1$ and $v_2$.

Eigenvector, eigenvalue: If the vector $v$ can be written as $Av = \lambda v$ for some matrix $A$ and constant $\lambda$, then $v$ is an eigenvector of $A$, and $\lambda$ is the associated eigenvalue.

Orthogonal vectors: Vectors are orthogonal if the angle between them is $90 \degree$ in two or three dimensions. The better definition of orthogonality is that the dot product between the two vectors is zero, and this holds for all dimensions. To calculate the dot product, multiply the elements in one vector by the corresponding element in the second vector and add up the result. If it sums to zero, the vectors are orthogonal.

Orthogonal matrix: A square matrix with real valued elements and with the property that its inverse is equal to its transpose is orthogonal.

Diagonal matrix: All elements not on the main diagonal are zero.

Transition matrix: A matrix that transforms vectors from one basis to another. Also called a change of coordinates matrix.

Characteristic polynomial: The polynomial arising from $\det (A - \lambda I) = 0$ where $\det$ is the determinant and $I$ is the identity matrix.

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.
• JupyterLab - JupyterLab is the latest web-based interactive development environment for notebooks, code, and data.
• Wolfram Language - Wolfram Language is a symbolic language, deliberately designed with the breadth and unity needed to develop powerful programs quickly.
• 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.

References and 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

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