# Homework #2 > # Name: # # Prefered Email: # # Collaborators: # # Due February 16th at 23:45:00. Submit electronically via Wolfware @ # http://submit.ncsu.edu/. Neither email nor paper submissions will # not be accepted. > > #Hint: Make liberal use of evalm() > with(linalg): Warning, the protected names norm and trace have been redefined and unprotected # Problem the first: # # a) For which value(s) of 'a' does the following matrix have rank 2? # Explain. # [ 1 2 3 4] # [ 0 (a+1) a 1] # [ 0 a a 0] > > A1 := matrix([[1,2,3,4],[0,a+1,a,1],[0,a,a,0]]); [1 2 3 4] [ ] A1 := [0 a + 1 a 1] [ ] [0 a a 0] > gausselim(A1); [1 2 3 4] [ ] [0 a a 0] [ ] [0 0 -1 1] > # Two cases: Either a=-1 or not. If it is, the Gaussian elimination > above is not nessesarily valid, but its easy to see the matrix has > rank 3. Otherwise looking at the above REF its easy to see the rank is > 3 unless a=0. > # 3/3 # b) For which values of 'a' is the following matrix in row echelon # form? Explain. # [ 1 2 3 4] # [ 0 (a-1) a 1] # [ 0 a^2-1 a-1 0] > A2 := matrix([[1,2,3,4],[0,a-1,a,1],[0,a^2-1,a-1,0]]); [1 2 3 4] [ ] A2 := [0 a - 1 a 1] [ ] [ 2 ] [0 a - 1 a - 1 0] > # It is in Row echelon form when a^2-1 = 0, so a=1,-1 both work. > # 3/3 # Given the matrix # [ 1/2 6 1 ] # A = [ 2 0 3 ] # [ 4 1 0 ] # # c) Compute A', the transpose of A. > > A3 := matrix(3,3,[1/2, 6, 1, 2, 0, 3, 4, 1, 0]); [1/2 6 1] [ ] A3 := [ 2 0 3] [ ] [ 4 1 0] > transpose(A3); [1/2 2 4] [ ] [ 6 0 1] [ ] [ 1 3 0] > # 2/2 > # Problem the second: # # a) Is it true that for all square matrices A, B that (AB)^2 = A^2 B^2? # Explain your answer, and give examples/counter-examples. > > # It is not true > > # NOTE: The fact that AB != BA is not enough to say this isn't true. > It is possible that AB != BA but that > (AB)^2 = A^2 B^2. For example: > E := matrix([[0,0],[a,b]]);F:= matrix([[0,0],[c,0]]); [0 0] E := [ ] [a b] [0 0] F := [ ] [c 0] > evalm(E &* F - F &* E); # EF != FE [ 0 0] [ ] [b c 0] > evalm((E&*F)^2 - E^2 &* F^2); # (EF)^2 = E^2 F^2 [0 0] [ ] [0 0] > # Here is a counter example: > A4 := matrix(2,2,[0,1,1,0]); B4 := matrix(2,2,[1,1,0,1]); [0 1] A4 := [ ] [1 0] [1 1] B4 := [ ] [0 1] > evalm((A4 &*B4)^2);evalm(A4^2 &* B4^2); [1 1] [ ] [1 2] [1 2] [ ] [0 1] > # 3/3 # b) Given the binary operation [,]:R^(n,n) x R^(n,n) -> R^(n,n) defined # by [A,B] = A*B - B*A decide whether this operation is commutative. # (i.e. if it is, show it algebraically, otherwise give an example which # shows it is not). > brac := (A,B) -> evalm(A &* B - B &* A); brac := (A, B) -> evalm(`&*`(A, B) - `&*`(B, A)) > brac(A4,B4); > brac(B4,A4); [-1 0] [ ] [ 0 1] [1 0] [ ] [0 -1] > # No it is not. > # 3/3 > # [Bonus] Is the operation in part (b) associative? (if it is, show it # algebraically, otherwise give an example which shows it is not). > > # No it is not > C4 := matrix(2,2,[0,1,1,1]); [0 1] C4 := [ ] [1 1] > brac(brac(A4,B4),C4); [0 -2] [ ] [2 0] > brac(A4,brac(B4,C4)); [-1 -2] [ ] [ 2 1] > # 1 bonus > # Problem the third: # # Suppose the Fibonacci rabbits take 2 months to breed, and that each # pair gives birth to 2 new pairs of rabbits each time they breed. # # a) If we start with 1 pair of baby rabbits at month 0, how many # rabbits in months 1-5? > > 1 # (1 new) Month 0 > 1 # (1 1m old) Month 1 > 1 # (1 mature) Month 2 > 3 # (1 mature, 2 new) Month 3 > 5 # (1 mature, 2 1m old, 2 new) Month 4 > 7 # (3 mature, 2 1m old, 2 new) Month 5 > 13 # (5 mature, 2 1m old, 6 new) Month 6 > 23 # (7 mature, 6 1m old, 10 new) Month 7 > > # 2/2 > # b) Write a formula for the number of rabbits in month n in terms of # the number of rabbits in previous months (n >= 3). > > r[n] = r[n-1] + 2*r[n-3]; r[n] = r[n - 1] + 2 r[n - 3] > # 2/2 > # c) Write a matrix (it should be 3x3) which will compute the population # at month n given the population in the previous three months.Use # matrix powering to compute the population in month 20. > R := matrix(3,3,[1,0,2,1,0,0,0,1,0]); [1 0 2] [ ] R := [1 0 0] [ ] [0 1 0] > evalm(R^18 &* vector([1,1,1])); [21205, 12503, 7373] > # Population in month 20 is 21205 > # 4/4 # Problem the fourth: # # Given a polynomial f(x) = a_n * x^n + a_(n-1) * x^(n-1) + ... + a_2 # x^2 + a_1 x + a_0 (where a_i's are real numbers) we define the # polynomial evaluated on the m x m square matrix B as follows: f(B) = # a_n * B^n + ... + a_2 * B^2 + a_1 * B + a_0* I_m. # # e.g. If f(x) = 3*x^2 - 6, then f(A) = 3*A^2 - 6* I_m # # If B = matrix([[0, 0, -4], [1, 0, 1], [0, 1, 0]]): # # a) Compute g(B), where g(x) = x^2 - x + 1 > > B := matrix([[0, 0, -4], [1, 0, 1], [0, 1, 0]]); I3 := band([1],3); [0 0 -4] [ ] B := [1 0 1] [ ] [0 1 0] [1 0 0] [ ] I3 := [0 1 0] [ ] [0 0 1] > evalm(B^2 - B + I3); [ 1 -4 4] [ ] [-1 2 -5] [ ] [ 1 -1 2] > # b) Compute h(B), where h(x) = 3*x^3 - 3*x + 12 > > evalm(3*B^3 - 3*B + 12 * I3); [0 0 0] [ ] [0 0 0] [ ] [0 0 0] > > # If C = [ 0 -4 ] # [ 1 2 ] # # c) Find a polynomial f(x) = x^2 + b x + c with f(C) = 0^(2x2) (Hint: # write down f(C) as a matrix with b and c constants, Maple will do this # for you quite easily). Can you think of another polynomial which when # evaluated on C would give 0_(2x2)? > > C := matrix(2,2,[0,-4,1,2]);I2 := band([1],2); [0 -4] C := [ ] [1 2] [1 0] I2 := [ ] [0 1] > fC := evalm(C^2+b*C + c*I2); [-4 + c -8 - 4 b] fC := [ ] [2 + b 2 b + c ] > solve({fC[1,1] = 0,fC[1,2] = 0,fC[2,1] = 0,fC[2,2] = 0},{b,c}); {c = 4, b = -2} > # The polynomial x^2 - 2*x + 4 works. Any multiple of this polynomial > will work as well, e.g. x^3 - 2*x^2 + 4*x > evalm(C^2 - 2*C + 4*I2);#check [0 0] [ ] [0 0] > # Whole problem 3/3 if you tried all parts. # Problem the last: # # Given the matrix # [ 1/2 6 1 ] # A = [ 2 0 3 ] # [ 4 1 0 ] # # a) Is it true that B^(n+m) = B^n B^m for any matrix B, and all # integers n and m? If it is not, provide a counter example. If it is, # sketch a short argument why. > > # Yes it is > # B^n * B^m = B * B * ... * B (n times) * B * B * ... * B (m times) = > B* ... *B (n+m times) = B^(n+m) > # b) Compute A^2, A^4, A^6, and A^7 using exactly one matrix # multiplication for each (this means no using ^ either). > > A2 := evalm(A &* A); [65/4 4 37/2] [ ] A2 := [ 13 15 2 ] [ ] [ 4 24 7 ] > A4 := evalm(A2 &* A2); [ 6241 ] [ ---- 569 3505/8] [ 16 ] A4 := [ ] [1657/4 325 569/2 ] [ ] [ 405 544 171 ] > A6 := evalm(A2 &* A4); [991233 365473 ] [------ 82441/4 ------ ] [ 64 32 ] [ ] A6 := [193513 ] [------ 13360 82441/8] [ 16 ] [ ] [57349/4 13884 19555/2] > A7 := evalm(A &* A6); [12115025 4948401 ] [-------- 834793/8 ------- ] [ 128 64 ] [ ] A7 := [2367609 834793 ] [------- 165745/2 ------ ] [ 32 16 ] [ ] [592373/8 95801 223957/4] > > # These two parts briefly checked: 3/3 # d) [More Difficult, don't worry too much if you can't get it.] In # class we talked about how to compute A^(2^k) efficiently. Consider # the binary representation of an integer: # # m = b0 + b1 * 2 + b2 * 2^2 + ... + bn * 2^n where each bi is 0 # or 1. (i.e. (bn b(n-1) .... b1 b0) are the binary digits of m) # # Keeping in mind parts (b) and (c), how might you adapt the method from # class to compute A^m? > > # Compute A^(2^n) and save A^(2^i) for each i with b1 = 1, then > multiply them together. > # i.e. Compute: > A^b0 &* A^(b1*2) &* ... &* A^(2^n) = A^m > # [Bonus] In the case where m = 31, compare the number of operations it # takes to compute A^m your way versus computing A*A*A*...*A. > > A^31; # This is 30 matrix multiplications = 30 * 27 = 810 rational > number multiplications > # The fast way is floor(log[2](31)) + HammingLenth(31) - 1 > matrix multiplications = (4 + 5 - 1) * 27 = 216 rational number > multiplications. > # 1 bonus # [Bonus] Implement this method in Maple. The following is a Maple # function to get a list of the binary digits of m (not guaranteed to be # good code, but might be useful): > > binaryDigits := proc(m) > # Call as: binaryDigits(yournumber); > local binary_rep, binary_digits, bit; > binary_rep := convert(m,binary); # Maple converts the number to binary > # Note that it still treats > binary_rep as > # as a decimal number! > > binary_digits := []; # Initialize > while binary_rep >= 1 do > bit := irem(binary_rep,10); # Grab the last 1 or 0 > binary_digits := [op(binary_digits),bit]; # Store the bit in our > list > binary_rep := (binary_rep - bit) / 10; # Strip the last bit, > and shift > od; > RETURN(binary_digits); # Note that the digits are returned in reverse > order > end: > # 2 bonus > binaryDigits(31); [1, 1, 1, 1, 1] > > power := proc(A, m) > local digs,B,C,i,mults; > mults := 0; # Number of matrix mulitplications taken > digs := binaryDigits(m); > B := evalm(A); # This will hold A^(2^i) at step i > # C will hold A^m at the end > if digs[1] = 1 then C := evalm(A) else C := 0 fi; > for i from 2 to nops(digs) > do > B := evalm(B^2); mults := mults +1;# B was A^(2^(i-2)), now its, > A^(2^(i-1)) > if(digs[i] = 1) then if C =0 then C := evalm(B) else C := > multiply(B,C); mults := mults +1 fi; fi; > od; > print("Total number of matrix multiplications = ",mults); > RETURN(evalm(C)); > end: > power(evalm(A),31); "Total number of matrix multiplications = ", 8 [541338620673075562492504265020369 [--------------------------------- , [ 2147483648 39060583168916814167808311105113 -------------------------------- , 134217728 206263610728944399843073217015089] ---------------------------------] 1073741824 ] [102514414576907179015113649174009 [--------------------------------- , [ 536870912 7396988616877824850313184455105 ------------------------------- , 33554432 39060583168916814167808311105113] --------------------------------] 268435456 ] [27481349836945059734482756956757 [-------------------------------- , [ 134217728 1982932231499698901478291814653 ------------------------------- , 8388608 10471094060754098832794720250165] --------------------------------] 67108864 ] > >