Fibonacci number sequence is a sequence where each number is a sum of the previous two numbers:

`1, 1, 2, 3, 5, 8, 13, 21, ....`

Therefore, it corresponds to a difference equation

`F`_{n+2} = F_{n+1} + F_{n}

or, equivalently, one can also write this as

`F`_{n} = F_{n-1} + F_{n-2}

with initial conditions

`F`_{0}=1, F_{1}=1.

Fibonacci is an Italian mathematician who came up with this sequence in year 1202, as a simple model of rabbit reproduction. Fibonacci sequence has a lot of interesting properties. If you are curious, there is a lot of information on the web, for instance here and also here

Now that we know how to solve linear 2-nd order difference equations, we can find

`F`_{n}

for
any `n`

. As you learned last week, such solution has the form:`r`_{1}

and `r`_{2}

are solutions of the quadratic equation obtained from the
original difference equation:
`r`_{1} = (1 + sqrt(5)) / 2

and `r`_{2} = (1 - sqrt(5)) / 2

.Now we only have to find

`c`_{1}

and `c`_{2}

, using the initial conditions:
`c`_{2} = 1 - c_{1}

. Plugging this value into the 2nd equation, we get
`c`_{2} = (1 - r_{2}) / (r_{1} - r_{2}) = r_{1} / sqrt(5)

Now let us finally write a Matlab program that calculates our solution

`F`_{n} = c_{1} r_{1}^{n} + c_{2} r_{2}^{n}

, and compares
it to the exact result, calculated directly from the definition of the Fibonacci sequence. In order to calculate `F`_{n}

directly, let's introduce
variables `Fnew=F`_{n}

, `Fold=F`_{n-1}

and `Folder=F`_{n-2}

, and swap
("rotate") the values of these variables after each iteration:
```
Folder = 1; Fold = 1;
```

r1 = (1 + sqrt(5)) / 2; r2 = (1 - sqrt(5) ) / 2;

c1 = r1 / sqrt(5); c2 = 1 - c1;

fprintf('\n n Direct_F(n) Solution_F(n) error \n\n');

for n=2:30

Fnew = Fold + Folder;

Fn = c1 * r1^n + c2 * r2^n;

fprintf('%2d %12d %12.0f %9.1e\n', n, Fnew, Fn, abs(Fnew-Fn));

Folder = Fold;

Fold = Fnew;

end;

Note that this is the most computationally efficient method, but not the simplest one. Here is a

`F`_{n}

, which avoids additional variables `Fold`

and `Folder`

:
```
F(1) = 1; F(2) = 2;
```

% Note that F2 = F1 + F0 = 2

r1 = (1 + sqrt(5)) / 2; r2 = (1 - sqrt(5) ) / 2;

c1 = r1 / sqrt(5); c2 = 1 - c1;

fprintf('\n n Direct_F(n) Solution_F(n) error \n\n');

for n=3:30

F(n) = F(n-1) + F(n-2);

Fn = c1 * r1^n + c2 * r2^n;

fprintf('%2d %12d %12.0f %9.1e\n', n, F(n), Fn, abs(F(n)-Fn));

end;

The only disadvantage of this second program is that it creates an array

`F(n)`

, where the Fibonacci numbers
will be stored; this uses up memory resources when `n`

is very large.
Note also that in this case we have to start with term `n=3`

instead of `n=2`

, since array indices
must be greater than zero (so an error would result in `F(n-2)`

when `n=2`

).
When running these programs, note how quickly the round-off error accumulates.

To practice solving non-homogeneous linear difference equations, we will change slightly the Fibonacci equation by adding a constant

`1`

to this equation:
`r`_{1}

and `r`_{2}

are solutions of the same quadratic equation that we already
solved, while the additional constant `c`

is the solution of the non-homogeneous equation. Therefore, to find `c`

,
we plug in `F`_{n}=F_{n-1}=F_{n-2}=c

into our equation, and obtain
`c = c + c + 1`

, from which we find `c=-1`

.Now, we have to once again determine

`c`_{1}

and `c`_{2}

using the initial conditions:
`c`_{2} = 2 - c_{1}

. Plugging this value into the 2nd equation, we get
`c`_{2} = 2 (1 - r_{2}) / (r_{1} - r_{2}) = 2 r_{1} / sqrt(5)

.
So our program will now read:```
Folder = 1; Fold = 1;
```

r1 = (1 + sqrt(5)) / 2; r2 = (1 - sqrt(5) ) / 2;

c1 = r1 / sqrt(5); c2 = 1 - c1;

fprintf('\n n Direct_F(n) Solution_F(n) error \n\n');

for n=2:30

Fnew = Fold + Folder;

Fn = c1 * r1^n + c2 * r2^n;

fprintf('%2d %12d %12.0f %9.1e\n', n, Fnew, Fn, abs(Fnew-Fn));

Folder = Fold;

Fold = Fnew;

end;

or, using the simpler but more memory-consuming method:

```
F(1) = 1; F(2) = 3; % Note that F2 = F1 + F0 + 1 = 3
```

r1 = (1 + sqrt(5)) / 2; r2 = (1 - sqrt(5) ) / 2;

c1 = 2 * r1 / sqrt(5); c2 = 2 - c1;

fprintf('\n n Direct_F(n) Solution_F(n) error \n\n');

for n=3:30

F(n) = F(n-1) + F(n-2) + 1;

Fn = c1 * r1^n + c2 * r2^n - 1;

fprintf('%2d %12d %12.0f %9.1e\n', n, F(n), Fn, abs(F(n)-Fn));

end;

If you understand these steps, homework #6 will be easy!