Help25a:=proc(): print(`RGx(n,N), Co(n,x), LE(P), LEeq(P,t), Com(n,N), RandomVeryGoodPoly(n), , IterPolygon1V(P,N), IterPolygonV(n,N) `): end: Help:=proc(): print(`RandM(m,n,N),IterPolygon(n,N), IterPolygon1(P,N), RandomGoodPoly(n), `): print(`A(n), Nor(v), EVC(M) , Champs(L) , CheckEVC1(M,lambda,v), CheckEVC(M,L) , Coeffs(V,M)`): print(`CheckCoeffs(V,M), CB(M,V) `): end: with(plots): with(linalg): #######From Yukun Yao's homework #OK to post homework #Yukun Yao, Apr 14, 2018, Assignment 22 ##Problem 2## For the random polygon project, get a feel for it, by #Writing a procedure #RandomGoodPoly(n) #that inputs a positive integer and outputs a random polygon with n vertices, whose centroid is the origin. #Hint: generate two random vectors, x and y, and substract from them their averages, also call these new vectors x and y to make their sum to zero and let the polygon be [[x1,y1], ..., [xn,yn],[x1,y1]] #Write a procedure #AvePolygon(P) #that inputs a polygon, (given by its list of vertices), and outputs the polygon whose vertices are the midpoints of the edges. #Write a procedure #IterPolygon(n,N) #that inputs a RandomGoodPoly(n) and a large pos. integer N, and performs AvePolygon(P) N times, and also plots it (use the command plot(P)). Perform IterPolygon(20,200) many times and look at the picture. #To generate a random polygon usually the x and y are rational. So by scaling let x,y be integers and give a range for them. RandomGoodPoly:=proc(n) local i,S,Ave,ra,L: S:=[]: ra:=rand(0..1000): for i from 1 to n do S:=[op(S), [ra(),ra()]]: od: Ave:=convert(S,`+`)/n: L:=[Ave$n]: S:=S-L: S:=[op(S), S[1]]: end: AvePolygon:=proc(P) local n,S,i: n:=nops(P)-1: S:=[]: for i from 1 to n do S:=[op(S),(P[i]+P[i+1])/2]: od: S:=[op(S), S[1]]: end: IterPolygon:=proc(n,N) local P,i: P:=RandomGoodPoly(n): for i from 1 to N do P:=AvePolygon(P): od: plot(P): end: IterPolygon1old:=proc(P,N) local i,P1: P1:=P: for i from 1 to N do P1:=AvePolygon(P1): od: plot(P1): end: IterPolygon1:=proc(P,N) local i,P1: P1:=P: for i from 1 to N do P1:=AvePolygon(P1): od: plot(P1,scaling=constrained): end: #It seems that the picture will converge to an oval. #######End From Yukun Yao's homework #From Edna Jones' hw24 #A(n): Inputs a positive integer n, and outputs the n by n matrix (given as a list of lists), corresponding to the "averaging process" [x1,x2, ..., xn] ->[(x1+x2)/2, (x2+x3)/2, ...., (xn+x1)/2] A:=proc(n) local j: if n=1 then return [[1]]; else return [seq([0$(j-1), 1/2, 1/2, 0$(n-j-1)], j=1..(n-1)), [1/2, 0$(n-2), 1/2]]; fi: end: #End from Edna Jones' hw24 #RandM(m,n,N): a random m by n matrix with entries from 1 to N RandM:=proc(m,n,N) local ra,i,j: ra:=rand(1..N): [seq([seq(ra(),i=1..n)],j=1..m)]: end: Nor:=proc(v) local i,a: a:=evalf(add(abs(v[i])^2,i=1..nops(v))): if a<10^(-Digits+2) then RETURN(FAIL): fi: [seq(v[i]/a,i=1..nops(v))]: end: #Champs(L): given a list of pairs outputs the set indices of pairs with the largest first component according to absolute value Champs:=proc(L) local cha,rec,i: if nops(L)=0 then RETURN(FAIL): fi: cha:={1}: rec:=evalf(abs(L[1][1])): for i from 2 to nops(L) do if evalf(abs(L[i][1]))>rec then cha:={i}: rec:= evalf(abs(L[i][1])): elif abs(L[i][1])=rec then cha:=cha union {i}: fi: od: cha: end: ## #EVC(M): the list of [eigenvalue,eigenvector] in weakly decreasing order of the absolute value of the eigenvalues EVC:=proc(M) local i,L,L1,cha,surv: L:=[eigenvectors(M)]: if {seq(L[i][2],i=1..nops(L))}<>{1} then print(`Some eigenspaces are not one-dimensional`): RETURN(FAIL): fi: L:=[seq([L[i][1],Nor( convert(L[i][3][1],list) )],i=1..nops(L))]: L1:=[]: while L<>[] do cha:=Champs(L): L1:=[op(L1),seq(L[i], i in cha)]: surv:={seq(i,i=1..nops(L))} minus cha: L:=[seq(L[i],i in surv)]: od: L1: end: #CheckEVC1(M,lambda,v): checks whether lambda is an eigenvalue of M with corr. eigenvector v CheckEVC1:=proc(M,lambda,v) local i,j: [seq(add(M[i][j]*v[j],j=1..nops(v))-lambda*v[i],i=1..nops(v))]: end: #CheckEVC(M,L): checks that EVC is (probably) correct CheckEVC:=proc(M,L) local i: {seq(CheckEVC1(M,op(L[i])),i=1..nops(L))}: end: #From Yonah B-A's hw24, corrected by Edna J. Coeffs:=proc(V,M) local n, vars, eqns, res,a,i,j: n:=nops(V): vars := [seq(a[i], i = 1 .. n)]: eqns := {seq(add(a[i]*M[i][j], i = 1 .. n) = V[j], j = 1 .. n)}: res := solve(eqns, {op(vars)}): subs(res, vars): end: CheckCoeffs:=proc(V,M) local Y,i: Y:=Coeffs(V,M): add(Y[i]*M[i], i=1..nops(Y))-V: end: #End Yonah B-A #CB(M,V): inputs a matrix M and a vector V and outputs the list of coefficients #expressing V in terms of the eigenvectors of M (arranged according tpo EVC) CB:=proc(M,V) local M1,i: M1:=EVC(M): M1:=[seq(M1[i][2],i=1..nops(M1))]: Coeffs(V,M1): end: #RGx(n,N): a random good vectors RGx:=proc(n,N) local ra,i,x,av: ra:=rand(1..N): x:=[seq(ra(),i=1..n)]: av:=convert(x,`+`)/n: [seq(x[i]-av,i=1..n)]: end: #Co(n,x): the coefficient of the second eigenvector in x Co:=proc(n,x) local M,gu,i: gu:=EVC(evalf(A(n))): M:=[seq(gu[i][2],i=1..nops(gu))]: Coeffs(x,M)[2]: end: #LE(P): inputs a good polygon and outputs the the eccenticity and the angle (in degrees) of the limiting ellipse in the everaging process LE:=proc(P) local n,x,y,lu,lu1,i: n:=nops(P)-1: print(`ERRONEOUS!`): x:=[seq(P[i][1],i=1..n)]: y:=[seq(P[i][2],i=1..n)]: lu:=evalf([Co(n,x),Co(n,y)]): if abs(lu[1])