# 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]