# Lesson 31. Solving the Fibonacci Recurrence
# --------------------------------------------
#
> restart: read `doc/fib.def`;
Warning: new definition for norm
Warning: new definition for trace
[ 0 1 ]
M := [ ]
[ 1 1 ]
mpow := proc(n)
options remember;
if type(n,posint) then
if type(n,even) then evalm(mpow(1/2*n)^2)
else evalm(mpow(1/2*n-1/2)^2 &* M)
fi
elif type(n,negint) then evalm(1/mpow(-n))
elif type(n,`+`) then evalm(`&*`(op(map(mpow,n))))
elif type(n,`*`) then evalm(mpow(n/op(1,n))^op(1,n))
else matrix([[F[n-1],F[n]],[F[n],F[n]+F[n-1]]])
fi
end
fib := n -> mpow(n - 1)[2, 2]
# Another way to look at the Fibonacci numbers (and recurrence relations in general) is by using the eigenvalues
# of the matrix M.
> evs:= eigenvals(M);
1/2 1/2
evs := 1/2 + 1/2 5 , 1/2 - 1/2 5
# There are two distinct eigenvalues, so this 2 by 2 matrix is diagonalizable, i.e. there is a 2 by 2 matrix R so that
# R^(-1) M R is a diagonal matrix with the eigenvalues of M on the diagonal. The columns of matrix R are the
# corresponding eigenvectors. The "eigenvects" command in "linalg" tells you the eigenvectors as well as the
# eigenvalues. Specifying "radical" asks Maple to use explicit formulas involving square roots (if possible) rather
# than "RootOf" notation.
> ev:=eigenvects(M,radical);
1/2 1/2
ev := [1/2 + 1/2 5 , 1, {[ 1, 1/2 + 1/2 5 ]}],
1/2 1/2
[1/2 - 1/2 5 , 1, {[ 1, 1/2 - 1/2 5 ]}]
# This requires some untangling. The result is two lists. In each list, first comes the eigenvalue, then its
# multiplicity, then a set of vectors forming a basis for the eigenvectors for this eigenvalue. The columns of R are
# op(ev[1][3]) and op(ev[2][3]).
> op(ev[1][3]);
1/2
[ 1, 1/2 + 1/2 5 ]
# This rather strange format allows Maple to handle cases with eigenvalues of higher multiplicity. For example:
> matrix([[1,0,0],[1,1,0],[0,0,1]]);
[ 1 0 0 ]
[ ]
[ 1 1 0 ]
[ ]
[ 0 0 1 ]
> eigenvects(");
[1, 3, {[ 0, 1, 0 ], [ 0, 0, 1 ]}]
# This matrix has only one eigenvalue 1, of multiplicity 3. It has a two-dimensional space
# of eigenvectors for 1, with a basis consisting of [0,1,0] and [0,0,1].
#
# Now going back to M, we want to put the two eigenvectors together as columns to form a matrix. The
# "concat" command in "linalg" does this.
> R:= concat(op(ev[1][3]), op(ev[2][3]));
[ 1 1 ]
[ ]
R := [ 1/2 1/2 ]
[ 1/2 + 1/2 5 1/2 - 1/2 5 ]
> DM:= evalm(R^(-1) &* M &* R);
[ 1/2 ]
[ 1/2 + 1/2 5 0 ]
DM := [ ]
[ 1/2 ]
[ 0 1/2 - 1/2 5 ]
# Now the point of all this is that a power of a diagonal matrix is easy to compute.
> DM^n = map(t -> t^n, DM);
[ 1/2 n ]
n [ (1/2 + 1/2 5 ) 0 ]
DM = [ ]
[ 1/2 n ]
[ 0 (1/2 - 1/2 5 ) ]
# Since M = R DM R^(-1), we have M^n = R DM^n R^(-1).
> evalm(R &* rhs(") &* R^(-1));
1/2 1/2 1/2 1/2
[- 1/10 %2 5 + 1/2 %2 + 1/10 %1 5 + 1/2 %1, 1/5 %2 5 - 1/5 %1 5 ]
1/2 1/2 1/2 1/2
[1/5 %2 5 - 1/5 %1 5 , 1/10 %2 5 + 1/2 %2 - 1/10 %1 5 + 1/2 %1]
1/2 n
%1 := (1/2 - 1/2 5 )
1/2 n
%2 := (1/2 + 1/2 5 )
> soln:= F[n] = "[1,2];
1/2 n 1/2 1/2 n 1/2
soln := F[n] = 1/5 (1/2 + 1/2 5 ) 5 - 1/5 (1/2 - 1/2 5 ) 5
# This solves the recurrence relation, i.e. provides a formula for it. Let's check this formula for a few values of n.
> seq(expand(subs(n=kk, soln)), kk=0..5);
F[0] = 0, F[1] = 1, F[2] = 1, F[3] = 2, F[4] = 3, F[5] = 5
# The formula is not necessarily a good way to calculate particular values of F[n].
# If you want to do an exact calculation, you must expand out the n'th powers. Each
# will involve n+1 terms, of which the even numbered ones cancel and the odd ones are
# the same for both (...+...)^n and (...-...)^n, leaving about n/2 rational numbers to add.
# It's easier to use repeated squaring on the matrix, where the number of steps is essentially proportional to log(n).
# On the other hand, you could use "evalf", in which case the answer will only be approximate when n is large,
# unless you use a high value for Digits.
> evalf(subs(n=100,soln));
21
F[100] = .3542248538*10
> Digits:= 22: round(evalf(subs(n=100,rhs(soln))));
354224848179261915084
> fib(100);
354224848179261915075
> Digits:= 25: round(evalf(subs(n=100,rhs(soln))));
354224848179261915075
> Digits:= 10:
# The formula does, however, tell us about the asymptotics of the sequence.
> evalf(soln);
n n
F[n] = .4472135956 1.618033989 - .4472135956 (-.6180339890)
# Note that one eigenvalue is greater than 1, while the other has absolute value less than 1. As n -> infinity, the
# term for the first eigenvalue dominates, while the second tends to 0. So F[n] is approximately equal to the
# following:
> fapp:= unapply(op(1,rhs(soln)), n);
1/2 n 1/2
fapp := n -> 1/5 (1/2 + 1/2 5 ) 5
> fib(10), evalf(fapp(10));
55, 55.00363624
> fib(11), evalf(fapp(11));
89, 88.99775290
# For even n, F[n] is slightly less than fapp(n), for odd n slightly greater.
#
# We could use this to answer the following question:
# Given an integer x, find the greatest Fibonacci number F[n] <= x.
#
# We want F[n] <= x < F[n+1].
# We know that (for n >= 0) fapp(n) < F[n] + 1/2 < x + 1/2 while
# fapp(n+1) > F[n+1] - 1/2 > x - 1/2.
#
> napp:= solve(fapp(n)=x+1/2, n);
1/2
ln(1/2 (2 x + 1) 5 )
napp := ----------------------
1/2
ln(1/2 + 1/2 5 )
# In principle (at least for x > 2), the n we want should be the greatest integer <= napp.
# If x = F[n], fapp(n) is close to x while fapp(n+1) is greater than x+1, so napp should be between n and n+1.
# But roundoff error could result in napp being computed as slightly greater than an integer when it should really
# be slightly less. So we have to allow for the possibility that the answer is off by 1.
> invfib:= proc(x) local n;\
n:= floor(ln((2*x+1)*sqrt(5)/2)/ln((1+sqrt(5))/2));\
if fib(n) > x then n-1\
elif fib(n+1) <= x then n+1\
else n \
fi;\
end:
> invfib(200);
12
> invfib(0), invfib(1), invfib(2), invfib(3);
0, 2, 3, 4
# We can use this to represent a positive integer as a sum of Fibonacci numbers
# F[n]. This is called the Zeckendorf representation. The representation is unique if the
# Fibonacci numbers are non-consecutive (F[n]+F[n+1] can be replaced by F[n+2]), with 1 counted as F[2].
> zeck:= proc(x) local n;\
if x = 0 then 0 \
else \
n:= invfib(x);\
F[n] + zeck(x-fib(n));\
fi\
end:
> zeck(100);
F[11] + F[6] + F[4]
