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:

F0 = c1 r10 + c2 r20 = c1 + c2 = 1

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:

F0 = c1 r10 + c2 r20 + c = c1 + c2 - 1 = 1

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!