# 1.31 Lesson 31

```# Lesson 31.  Solving the Fibonacci Recurrence
# --------------------------------------------
#
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.

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]

```