## Linear 2nd order difference equation example: Fibonacci numbers

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 `Fn+2 = Fn+1 + Fn`
or, equivalently, one can also write this as `Fn = Fn-1 + Fn-2`
with initial conditions `F0=1, F1=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 `Fn` for any `n`. As you learned last week, such solution has the form:

#### Fn = c1 r1n + c2 r2n

where `r1` and `r2` are solutions of the quadratic equation obtained from the original difference equation:

#### r2 = r + 1

It is easy to see that `r1 = (1 + sqrt(5)) / 2` and `r2 = (1 - sqrt(5)) / 2`.
Now we only have to find `c1` and `c2`, using the initial conditions:

and

#### F1 = c1 r1 + c2 r2 = 1

The first equation yields `c2 = 1 - c1`. Plugging this value into the 2nd equation, we get `c2 = (1 - r2) / (r1 - r2) = r1 / sqrt(5)`

Now let us finally write a Matlab program that calculates our solution `Fn = c1 r1n + c2 r2n`, and compares it to the exact result, calculated directly from the definition of the Fibonacci sequence. In order to calculate `Fn` directly, let's introduce variables `Fnew=Fn`, `Fold=Fn-1` and `Folder=Fn-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 simpler way to directly calculate `Fn`, 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.

## Non-homogeneous linear 2nd order difference equation

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

#### Fn = Fn-1 + Fn-2 + 1

The solution will now have the form

#### Fn = c1 r1n + c2 r2n + c

where `r1` and `r2` 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 `Fn=Fn-1=Fn-2=c` into our equation, and obtain `c = c + c + 1`, from which we find `c=-1`.

Now, we have to once again determine `c1` and `c2` using the initial conditions:

and

#### F1 = c1 r1 + c2 r2 - 1 = 1

The first equation yields `c2 = 2 - c1`. Plugging this value into the 2nd equation, we get `c2 = 2 (1 - r2) / (r1 - r2) = 2 r1 / 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!