- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - MATH 2270, SEC. 3. 6TH LABORATORY ASSIGNMENT Treibergs Due Dec 10,1998. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > with(networks):with(linalg): Warning: new definition for charpoly Warning: new definition for rank Warning: new definition for charpoly Warning: new definition for rank In this lab assignment, we consider some fun applications of linear algebra, to networks and Markov chains. This will give us a chance to explore the networks package which deals with graph theory. This discussion follows Chapter 8.1 of Kolman closely, so it is a good idea to scan that part of the book before you proceed. 1. NETWORKS AND CLIQUES As our first application we consider graphs. Directed graphs may represent information flow from node to node in a network, or it may indicate influence relationships in any organization of individuals. In MAPLE, a graph is a new data structure. For example,we enter the directed graph pictured by figure 8.15 on p.429. First we name the graph G1. > new(G1): Then we enter the list of vertices of G1. > addvertex({1,2,3,4,5},G1); 1, 2, 3, 4, 5 Then we enter the directed edges, which we can visualize as an arrow connecting a tail vertex to a head vertex. The addedge command inserts the set of vertices, e.g., [1,2] means a directed edge from vertex 1 to vertex 2. Some edges indicate both arrows which means there are two directed edges, one from tail to head, the other from head to tail. If there is a two directional edge between vertices 1 and 2, then both directed edges [1,2] and [2,1] have to be entered in the list. > addedge({[1,3],[1,4],[1,5],[2,3],[2,5],[3,1],[3,5], [4,3],[4,5],[5,1],[5,2],[5,3]},G1); e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12 MAPLE will plot your graph. It will not draw the arrowheads, however, but you can add them by hand. > draw(G1); To apply linear algebra, we encode the digraph in a matrix. a[i,j] = 1 whenever there is a directed edge from vertex i to vertex j in the graph. Other entries are zero. > A:=adjacency(G1); [ 0 0 1 1 1 ] [ ] [ 0 0 1 0 1 ] [ ] A := [ 1 0 0 0 1 ] [ ] [ 0 0 1 0 1 ] [ ] [ 1 1 1 0 0 ] A CLIQUE (p. 425) in a digraph is a subset S1 of the vertices of G1 such that (a) S contains three or more vertices; (b) If Pi and Pj are in S then there is a directed edge from Pi to Pj and a directed edge from Pj to Pi; (c) there is no larger subset T of the vertices of G1 that satisfies property (b) and contains S. (So roughly speaking, cliques are the subsets of G1 such that every member of the clique can talk with every other member of the clique.) To find the ciques in the digraph, we follow the prescription of Theorem 8.2. First form a related matrix S such that for i ­ j, S[i,j] = 1 when both A[i,j]=1 and A[j,i]=1 but A[i,j]=0 otherwise. If you add A to its transpose, then S[i,j]=1 whenever the [i,j] element of the sum is 2. > add(A,transpose(A)); [ 0 0 2 1 2 ] [ ] [ 0 0 1 0 2 ] [ ] [ 2 1 0 1 2 ] [ ] [ 1 0 1 0 1 ] [ ] [ 2 2 2 1 0 ] If we divide each entry by 2 and take the greatest integer part we get S. The greatest integer is also called the floor function. floor(1)=1,floor(0.5)=0, floor(0)=0. > S:=map(floor,scalarmul(",1/2)); [ 0 0 1 0 1 ] [ ] [ 0 0 0 0 1 ] [ ] S := [ 1 0 0 0 1 ] [ ] [ 0 0 0 0 0 ] [ ] [ 1 1 1 0 0 ] The prescription is to cube the matrix S. Any nonzero diagonal element is in some clique. > evalm(S^3); [ 2 1 3 0 4 ] [ ] [ 1 0 1 0 3 ] [ ] [ 3 1 2 0 4 ] [ ] [ 0 0 0 0 0 ] [ ] [ 4 3 4 0 2 ] Thus vertices 1,3,5 are members of (the same) clique. 2. STRONGLY CONNECTED GRAPHS The property of strong connectedness is roughly that for any pair of individuals Pi and Pj in the organization, there is a sequence of distinct people (a path) [Pi,x1,x2,...,xk,Pj] such that Pi influences x1, x1 influences x2, and so on until xk influences Pj, and a possibly different set of distinct people [Pj,y1,y2,...,ym,Pi] such that Pj influences y1, y1 influences y2, and so on until ym influences Pi, More precisely, the digraph G is STRONGLY CONNECTED if for every two distinct vertices Pi, Pj of G, there is a directed path from Pi to Pj and a directed path from Pj to Pi. A DIRECTED PATH from Pi to Pj is a set of distinct vertices x1 ,x2, ...,xk such that there are directed edges [Pi,x1],[x1,x2],[x2,x3],...,[xk,Pj]. It turns out that a computation tells us if a digraph G is strongly connected by Theorem 8.3.One computes E = A + A^2 + ... + A^(n-1) where n is the number of vertices of G and A is the adjacency matrix of G. If E has no zero entries then G1 is strongly connected. In case of G1 above, n=5 so E is nonzero, thus G1 is strongly connected. > evalm(A+A^2+A^3+A^4); [ 16 8 21 7 21 ] [ ] [ 12 6 15 4 15 ] [ ] [ 14 7 16 5 16 ] [ ] [ 12 6 15 4 15 ] [ ] [ 16 8 21 7 21 ] Instead of having to type in the individual powers, we can run a little do-loop. j is a dummy variable which the do instruction causes to step for 2 to n-1. MS is a matrix which accumulates the sum. Let n denote the number of vertices. > n:=5; n := 5 > MS:=evalm(A); for j from 2 to n-1 do MS := evalm( MS + A^j); od; [ 0 0 1 1 1 ] [ ] [ 0 0 1 0 1 ] [ ] MS := [ 1 0 0 0 1 ] [ ] [ 0 0 1 0 1 ] [ ] [ 1 1 1 0 0 ] [ 2 1 3 1 3 ] [ ] [ 2 1 2 0 2 ] [ ] MS := [ 2 1 2 1 2 ] [ ] [ 2 1 2 0 2 ] [ ] [ 2 1 3 1 3 ] [ 6 3 8 3 8 ] [ ] [ 4 2 6 2 6 ] [ ] MS := [ 5 2 6 2 7 ] [ ] [ 4 2 6 2 6 ] [ ] [ 7 4 8 2 7 ] [ 16 8 21 7 21 ] [ ] [ 12 6 15 4 15 ] [ ] MS := [ 14 7 16 5 16 ] [ ] [ 12 6 15 4 15 ] [ ] [ 16 8 21 7 21 ] The next digraph G2 is figure 8.14 on p.428. > new(G2): > addvertex({1,2,3,4,5},G2); 1, 2, 3, 4, 5 > addedge({[1,2],[1,3],[1,5],[2,1],[2,5],[3,1],[3,4], [4,2],[4,3],[5,1]},G2); e1, e2, e3, e4, e5, e6, e7, e8, e9, e10 > A:=adjacency(G2); [ 0 1 1 0 1 ] [ ] [ 1 0 0 0 1 ] [ ] A := [ 1 0 0 1 0 ] [ ] [ 0 1 1 0 0 ] [ ] [ 1 0 0 0 0 ] > draw(G2); > n:=5;MS:=evalm(A): for j from 2 to n-1 do MS := evalm( MS + A^j): od:evalm(MS); n := 5 [ 15 6 6 5 10 ] [ ] [ 9 6 6 2 8 ] [ ] [ 8 9 9 3 8 ] [ ] [ 11 5 5 4 7 ] [ ] [ 5 5 5 1 5 ] The next is 8.17, p431. > new(G3): > addvertex({1,2,3,4,5},G3); 1, 2, 3, 4, 5 > addedge([[1,2],[1,4],[2,3],[2,4],[3,4],[3,5],[4,1],[4,3],[5,2],[5,4]],G3); e1, e2, e3, e4, e5, e6, e7, e8, e9, e10 > draw(G3); > A:=adjacency(G3); [ 0 1 0 1 0 ] [ ] [ 0 0 1 1 0 ] [ ] A := [ 0 0 0 1 1 ] [ ] [ 1 0 1 0 0 ] [ ] [ 0 1 0 1 0 ] > MS:=A:for j from 2 to 4 do MS:=evalm(MS+A^j): od:evalm(MS); [ 5 5 7 10 3 ] [ ] [ 5 4 8 10 3 ] [ ] [ 5 4 7 10 4 ] [ ] [ 5 4 7 10 4 ] [ ] [ 5 5 7 10 3 ] So G3 is strongly connected. > S:=map(floor, scalarmul(evalm(A+transpose(A)),1/2)); [ 0 0 0 1 0 ] [ ] [ 0 0 0 0 0 ] [ ] S := [ 0 0 0 1 0 ] [ ] [ 1 0 1 0 0 ] [ ] [ 0 0 0 0 0 ] > evalm(S^3); [ 0 0 0 2 0 ] [ ] [ 0 0 0 0 0 ] [ ] [ 0 0 0 2 0 ] [ ] [ 2 0 2 0 0 ] [ ] [ 0 0 0 0 0 ] But G3 has no cliques. PROBLEMS Description of Directed Graph The vertices form a square with labels P1,P3,P4,P6 and a triangle with labels P2,P5,P7. All six pairs of the square are connected by edges directed in both ways (P1=P3,P1=P4,P1=P6,P3=P4,P3=P6,P4=P6) and each edge of the triangle is connected with edges directed both ways (P2=P5,P2=P7,P5=P7). Finally there are two one way connecting edges from the square to the triangle (P1-->P2, P6-->P7). 1. Using MAPLE determine the cliques of this digraph. 2. Using MAPLE, determine if this digraph is strongly connected.> 3. MARKOV CHAINS It is probably a good idea to start the last part of the lab with a new file. > with(linalg): This question is based on chapter 8.3 of Kolman. Markov chains describe the evolution of a state vector. There are many applications as you'll see in the chapter but for concreteness, let's describe a model for jobs in a company. The state vector x = [x1,x2,...,xn] gives the number of people in job category 1 to job category n during this year. It is assumed that in one year there is a probability t[i,j] that a person from job j this year will be in job i next year. Since t[i,j] is a probability, 0 ² t[i,j] ² 1. Because a person is certain to be in one of the jobs next year we must have 1 = t[1,j]+...+t[n,j]. For example, if there are two job categories, technical and accounting, we can define the transition probabilities matrix > T:=matrix([[4/5,1/4],[1/5,3/4]]); [ 4/5 1/4 ] T := [ ] [ 1/5 3/4 ] This means that somebody in job 1 (technical) has .8 chance of being technical next year and 0.2 chance of switching to accounting and a accounting worker has 0.75 chance of remaning accounting and 0.25 chance of switching to technical. If x is the vector of workers in the technical and accounting jobs, then by Theorem 8.4, after k years, there will be T^k x technical workers in each department. For example, if there are initially > x:=vector([200.0,100.0]); x := [ 200.0, 100.0 ] technical and accounting workers then we can watch the numbers redistribute according to this model. > evalm(T&* x); [ 185.0000000, 115.0000000 ] > evalm(T^2 &* x); [ 176.7500000, 123.2500000 ] > seq(evalm(T^k&* x),k=1..8); [ 185.0000000, 115.0000000 ], [ 176.7500000, 123.2500000 ], [ 172.2125000, 127.7875000 ], [ 169.7168750, 130.2831250 ], [ 168.3442813, 131.6557188 ], [ 167.5893547, 132.4106453 ], [ 167.1741451, 132.8258549 ], [ 166.9457798, 133.0542202 ] It seems to converge. Accelerating the powers a bit. > seq(evalm(T^(2^k) &* x),k=3..7); [ 166.9457798, 133.0542202 ], [ 166.6690038, 133.3309962 ], [ 166.6666668, 133.3333332 ], [ 166.6666667, 133.3333333 ], [ 166.6666667, 133.3333333 ] By theorem 8.6, the sequence T^k x converges to a vector u which is a steady state vector provided that T is a regular matrix. To be a steady state vector means that Tu = u or (T-I)u=0. Thus we set M = T-I and we may solve for u in Mu=0: > M3:=evalm(T-band([1],2));linsolve(M3,vector([0,0])); [ -1/5 1/4 ] M3 := [ ] [ 1/5 -1/4 ] [ 5/4 _t[1], _t[1] ] This means that the company will converge to 55% technicians and 45% accountants. The matrices also converge: > seq(evalm(T^k),k=1..7); [ 69 31 ] [ 1259 741 ] [ 23849 16151 ] [ --- ---- ] [ ---- ---- ] [ ----- ----- ] [ 4/5 1/4 ] [ 100 80 ] [ 2000 1600 ] [ 40000 32000 ] [ ], [ ], [ ], [ ], [ 1/5 3/4 ] [ 31 49 ] [ 741 859 ] [ 16151 15849 ] [ --- ---- ] [ ---- ---- ] [ ----- ----- ] [ 100 80 ] [ 2000 1600 ] [ 40000 32000 ] [ 462339 337661 ] [ 9085729 6914271 ] [ 179943019 140056981 ] [ ------ ------ ] [ -------- -------- ] [ --------- --------- ] [ 800000 640000 ] [ 16000000 12800000 ] [ 320000000 256000000 ] [ ], [ ], [ ] [ 337661 302339 ] [ 6914271 5885729 ] [ 140056981 115943019 ] [ ------ ------ ] [ -------- -------- ] [ --------- --------- ] [ 800000 640000 ] [ 16000000 12800000 ] [ 320000000 256000000 ] > evalf("); [ .8000000000 .2500000000 ] [ .6900000000 .3875000000 ] [ ], [ ], [ .2000000000 .7500000000 ] [ .3100000000 .6125000000 ] [ .6295000000 .4631250000 ] [ .5962250000 .5047187500 ] [ ], [ ], [ .3705000000 .5368750000 ] [ .4037750000 .4952812500 ] [ .5779237500 .5275953125 ] [ .5678580625 .5401774219 ] [ ], [ ], [ .4220762500 .4724046875 ] [ .4321419375 .4598225781 ] [ .5623219344 .5470975820 ] [ ] [ .4376780656 .4529024180 ] Each of the columns converge to the solution of Tu=u which satisfies u[1]+...+u[n]=1.We can impose this equation as another row in the matrix equation. > MA:=stack(M3,matrix([[1,1]]));b:=vector([0,0,1]); [ -1/5 1/4 ] [ ] MA := [ 1/5 -1/4 ] [ ] [ 1 1 ] b := [ 0, 0, 1 ] > u:=linsolve(MA,b); u := [ 5/9, 4/9 ] The condition for regular transitions matrix is the following. Some power of T has to have all nonzero entries. This T is already regular since T^1 has nonzero entries. To end the discussion, lets consider another couple of transitions probability matrices. We check if they are regular and find the steady state probability distribution vector u. > T1:=matrix([[1/5,0,1/5,0],[0,3/4,0,1/6],[4/5,0,4/5,0],[0,1/4,0,5/6]]); [ 1/5 0 1/5 0 ] [ ] [ 0 3/4 0 1/6 ] T1 := [ ] [ 4/5 0 4/5 0 ] [ ] [ 0 1/4 0 5/6 ] First check if this is a transitions probability matrix by summing all the columns > ones:=matrix([[1,1,1,1]]);evalm(ones &* T1); ones := [ 1 1 1 1 ] [ 1 1 1 1 ] > seq(evalm(T1^k),k=2..6); [ 1/5 0 1/5 0 ] [ 1/5 0 1/5 0 ] [ ] [ ] [ 29 19 ] [ 299 277 ] [ 0 ---- 0 ---- ] [ 0 --- 0 --- ] [ 48 72 ] [ 576 864 ] [ ], [ ], [ 4/5 0 4/5 0 ] [ 4/5 0 4/5 0 ] [ ] [ ] [ 19 53 ] [ 277 587 ] [ 0 ---- 0 ---- ] [ 0 --- 0 --- ] [ 48 72 ] [ 576 864 ] [ 1/5 0 1/5 0 ] [ 1/5 0 1/5 0 ] [ ] [ ] [ 3245 3667 ] [ 36539 46405 ] [ 0 ---- 0 ----- ] [ 0 ----- 0 ------ ] [ 6912 10368 ] [ 82944 124416 ] [ ], [ ], [ 4/5 0 4/5 0 ] [ 4/5 0 4/5 0 ] [ ] [ ] [ 3667 6701 ] [ 46405 78011 ] [ 0 ---- 0 ----- ] [ 0 ----- 0 ------ ] [ 6912 10368 ] [ 82944 124416 ] [ 1/5 0 1/5 0 ] [ ] [ 421661 573667 ] [ 0 ------ 0 ------- ] [ 995328 1492992 ] [ ] [ 4/5 0 4/5 0 ] [ ] [ 573667 919325 ] [ 0 ------ 0 ------- ] [ 995328 1492992 ] So this powers of this matrix will never have all nozero entries so T1 is not regular. > T2:=matrix([[4/5,0,1/3,0],[0,3/4,0,1/6],[0,1/4,2/3,0],[1/5,0,0,5/6]]); [ 4/5 0 1/3 0 ] [ ] [ 0 3/4 0 1/6 ] T2 := [ ] [ 0 1/4 2/3 0 ] [ ] [ 1/5 0 0 5/6 ] > evalm(ones &* T2),evalm(T2^2),evalm(T2^3); [ 16 22 ] [ 64 133 364 ] [ ---- 1/12 ---- 0 ] [ --- --- --- 1/72 ] [ 25 45 ] [ 125 720 675 ] [ ] [ ] [ 19 ] [ 143 27 271 ] [ 1/30 9/16 0 ---- ] [ ---- ---- 1/90 --- ] [ 72 ] [ 1800 64 864 ] [ 1 1 1 1 ], [ ], [ ] [ 17 ] [ 217 ] [ 0 ---- 4/9 1/24 ] [ 1/120 --- 8/27 3/32 ] [ 48 ] [ 576 ] [ ] [ ] [ 49 25 ] [ 1801 23 125 ] [ --- 0 1/15 ---- ] [ ---- 1/60 --- --- ] [ 150 36 ] [ 4500 150 216 ] thus T2 is regular since T2^3 has all nonzero entries. We now find the probability distribution of the steady state u. > M2:=stack(evalm(T2-band([1],4)),matrix([[1,1,1,1]])); b2:=vector([0,0,0,0,1]); [ -1/5 0 1/3 0 ] [ ] [ 0 -1/4 0 1/6 ] [ ] M2 := [ 0 1/4 -1/3 0 ] [ ] [ 1/5 0 0 -1/6 ] [ ] [ 1 1 1 1 ] b2 := [ 0, 0, 0, 0, 1 ] > linsolve(M2,b2); [ 5/18, 2/9, 1/6, 1/3 ] PROBLEM 3. Let the transitions probability matrix be given by T. > T:=matrix([[9/10,0,0,1/7,0],[0,8/9,0,0,1/6],[1/10,0,7/8,0,0], [0,1/9,0,6/7,0],[0,0,1/8,0,5/6]]); [ 9/10 0 0 1/7 0 ] [ ] [ 0 8/9 0 0 1/6 ] [ ] T := [ 1/10 0 7/8 0 0 ] [ ] [ 0 1/9 0 6/7 0 ] [ ] [ 0 0 1/8 0 5/6 ] > > Check that T is a regular matrix and find the steady state vector u. Supposing that initially there are [100,200,300,400,500] people in the various jobs of this company, what is the limiting (steady state) distribution of people? (Hint: What happens to the total number of people going from year to year?)