Fibonacci Numbers

The Fibonacci numbers are a famous sequence of numbers defined by the following recursion rule:

$$ \begin{matrix} F_0 &=& 0 &\\ F_1 &=& 1 &\\ F_n &=& F_{n-1} + F_{n-2} &\text{for } n>2 \end{matrix} $$

The first few Fibonacci numbers are: $0,1,1,2,3,5,8,13,21, \dots$

Fibonacci numbers are fascinating because they appear in many areas of mathematics. They are also closely related to the golden ratio:

$$ \varphi = \frac{1 + \sqrt{5}}{2} = 1.61803\dots $$

The ratio of successive Fibonacci numbers $F_{n+1}/F_n$ approximates $\varphi$ as $n$ gets larger.

Today, we will see three different ways of computing $F_n$ using Python:

  1. Directly, using the formula for $F_n$
  2. By taking powers of a matrix $F$
  3. By finding the diagonal decomposition of $F$, and then taking powers

Computing using recursion

The rule $F_n = F_{n-1} + F_{n-2}$ is sometimes called a recursive rule. This just means that future values depend on past values according to some fixed formula.

Let's implement this in Python, and use this to print the first 10 Fibonacci numbers. Can you understand what the code is doing?

In [1]:
def fib(n):
    if n == 0:
        return 0
    elif n == 1:
        return 1
    else:
        return fib(n-1) + fib(n-2)

for i in range(10):
    print(fib(i))
0
1
1
2
3
5
8
13
21
34

Computing using matrix powers

As we saw in our lecture, the Fibonacci numbers are also given by the matrix equation

$$ \begin{pmatrix} F_{n+1} \\ F_n \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} F_{n} \\ F_{n-1} \end{pmatrix} $$

We can also write this as:

$$ \begin{pmatrix} F_{n+1} \\ F_n \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n \begin{pmatrix} F_1 \\ F_0 \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n \begin{pmatrix} 1 \\ 0 \end{pmatrix} $$

Let $F$ be the matrix $\begin{pmatrix} 1& 1 \\ 1 & 0 \end{pmatrix}$.

Taking powers of $F$ gives us another way to computing $F_n$. We can do this in Python using the command np.linalg.matrix_power(F,n). In our code, the variable Fn_vector is the vector $(F_{n+1}, F_n)$, but since we only want $F_n$, we extract entry $1$ (recall that Python indexing starts at 0). Let's print the first 10 Fibonacci numbers to make sure we've got it right:

In [2]:
import numpy as np

F = np.array([[1,1], [1,0]])

def fib_mat(n):    
    Fn_vector = np.linalg.matrix_power(F,n) @ np.array([1,0]) 
    
    return Fn_vector[1]

for i in range(10):
    print(fib_mat(i))
0
1
1
2
3
5
8
13
21
34

This is nice. We can even compute negative Fibonacci numbers! (Why does this make sense?)

In [3]:
for i in range(-5, 5):
    print(fib_mat(i))
5.0
-3.0
2.0
-1.0
1.0
0
1
1
2
3

But calculating $F^n$ requires calculating the previous powers too, so fib_mat isn't really that different from fib. Here are the first ten powers of $F$:

In [4]:
for i in range(10):
    print(np.linalg.matrix_power(F,i), "\n")
[[1 0]
 [0 1]] 

[[1 1]
 [1 0]] 

[[2 1]
 [1 1]] 

[[3 2]
 [2 1]] 

[[5 3]
 [3 2]] 

[[8 5]
 [5 3]] 

[[13  8]
 [ 8  5]] 

[[21 13]
 [13  8]] 

[[34 21]
 [21 13]] 

[[55 34]
 [34 21]] 

Computing using the diagonal decomposition of $F$

It turns out that $F$ can be diagonally decomposed! We can compute the diagonal decomposition, also called the eigendecomposition, using the np.linalg.eig command:

In [5]:
np.linalg.eig(F)
Out[5]:
(array([ 1.61803399, -0.61803399]),
 array([[ 0.85065081, -0.52573111],
        [ 0.52573111,  0.85065081]]))

This function returns 2 things:

  1. A vector $d$
  2. A matrix $M$

Let's store these into variables d and M:

In [6]:
d, M = np.linalg.eig(F)

Let's also define a diagonal matrix D that given by d along the diagonal:

In [7]:
D = np.diag(d)
D
Out[7]:
array([[ 1.61803399,  0.        ],
       [ 0.        , -0.61803399]])

Then $F = M D M^{-1}$, as we may verify:

In [8]:
M @ D @ np.linalg.inv(M)
Out[8]:
array([[ 1.00000000e+00,  1.00000000e+00],
       [ 1.00000000e+00, -1.11022302e-16]])

Actually, it turns out that $M$ is orthogonal!

In [9]:
M @ M.T
Out[9]:
array([[1., 0.],
       [0., 1.]])

So $F$ is given by $M D M^T$:

In [10]:
M @ D @ M.T
Out[10]:
array([[ 1.00000000e+00,  1.00000000e+00],
       [ 1.00000000e+00, -5.55111512e-17]])

This week's homework is very simple: use d and M to compute $F_n$.

Homework

Copy the following script into Spyder, then modify it so that the function my_fib(n) gives $F_n$. Test your function by printing the first 10 Fibonacci numbers (don't worry about rounding errors).

In [70]:
import numpy as np

F = np.array([[1,1],[1,0]])
d,M = ????? # diagonal decomposition of F, stored into a vector v and a matrix M

def my_fib(n):
    dn = ????? # vector with entries of d to the power of n
    Dn = ????? # diagonal matrix with dn along the diagonal
    Fn = ????? # the matrix F to the power of n. Do not use np.linalg.matrix_power or np.linalg.inv!
    
    Fn_vector = Fn @ np.array([1,0])
    
    return Fn_vector[1]

for i in range(10):
    print(my_fib(i))

It may help to trying defining dn,Dn,Fn outside the function first, so you can see what they look like and whether they're correct.

Also remember that Python uses "v**n" to denote powers. This works even when v is a vector!

.