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

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:

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

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))
```

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))
```

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))
```

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")
```

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]:

This function returns 2 things:

- A vector $d$
- 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]:

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

In [8]:

```
M @ D @ np.linalg.inv(M)
```

Out[8]:

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

In [9]:

```
M @ M.T
```

Out[9]:

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

In [10]:

```
M @ D @ M.T
```

Out[10]:

This week's homework is very simple: use `d`

and `M`

to compute $F_n$.

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!

.