# OK to post homework
# Edna Jones, Math 640 (Spring 2018), Homework Assignment 25

# 1. Prove that 1 is the an eigenvalue of A(n) and the all 1 vector is an eigenvector

# Let v be the n-dimensional all 1's vector. Then A(n)*v = [1/2+1/2, 1/2+1/2,..., 1/2+1/2]^T = [1, 1,..., 1]^T = v = 1*v, where x^T is the transpose of x. Thus, 1 is an eigenvalue of A(n) and v is an eigenvector of A(n).


# 2. [Added April 20,2018] plot the eigenvectors corresponding to the second-and-third largest (i.e. the largest ones that complex) eigenvectors of A(n) for n=30. Recall that its components are complex numbers and the complex number x+y*I corresponds to the point [x,y]. What kind of shape emerges? can you conjecture an equation for it?

Digits := 20;

### From my hw24EdnaJones.txt
A := proc(n)
    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:
###

### From my C25.txt
with(linalg):

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]/sqrt(a),i=1..nops(v))]: # I replaced a with sqrt(a)
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, j:
    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[j], j in cha)];
        surv := {seq(j, j=1..nops(L))} minus cha;
        L := [seq(L[j], j in surv)];
    od:
    
    L1;
end:
###

### From my hw24EdnaJones.txt
#### From Dr. Z.'s C24.txt
#MtV(M,V): inputs a matrix M and a column vector V (expressed as a list), outputs the vector MV
#For example, try
#MtV([[1,2],[3,4]],[5,6]); MtV([[a11,a12],[a21,a22]],[b1,b2]);
MtV:=proc(M,V) local i,j:

if not (type(M,list) and {seq(type(M[i],list),i=1..nops(M))}={true} and nops({seq(nops(M[i]),i=1..nops(M))})=1 and type(V,list) and
         nops(M)>0 and nops(M[1])=nops(V)) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

[seq(add(M[i][j]*V[j],j=1..nops(M[i])),i=1..nops(M))]:

end:
####

Coeffs := proc(V,M)
    return MtV(convert(inverse(transpose(M)), list),V);
end:
###

# complexplot(EVC(evalf(A(30)))[2][2]) and complexplot(EVC(evalf(A(30)))[3][2]) look like the circle x^2 + y^2 = 1/30.



# 3. [Possibly a challenge] Explain why the limiting shape is an ellipse 
# [ Added April 21, 2018: previously I claimed that the inclination is 45 degrees, that was wrong. The inclination can be anything. The shape is an ellipse whose center is the origin, of course, but otherwise can be any shape and any orientation]

# I suspect this has something to do with the eigenvectors of the eigenvalues with the second largest magnitude having components close on a circle.
# Since the centroid of any good polygon is [0,0], the coefficient corresponding to the all 1's vector for a good polygon will be zero, so the limiting shape is determined for the eigenvectors of the eigenvalues with the second largest magnitude.


# 4. Write a program that inputs a good polygon and outputs the eccentricity [Added April 21, 2018: and orientation] of the limiting ellipse, by just using the initial points of the polygon (i.e. the vector of x-coordinates and the vector of y-coordinates, both of which add-up to 0). 
# [Correction April 21, previously I said "add up to 1", of course that was wrong, by construction of a good polygon the x and y vectors add up to 0]
# Hints: (1) use Coeffs (2) Convince yourself that if A and B are any complex numbers, the curve x=Real(A*exp(I*t)), y=Real(B*exp(I*t)) is an ellipse. By writing A=|A|*exp(I*arg(A)) and B=|B|*exp(I*arg(B)), can you find the eccentricity and orientation of that ellipse?

### From my hw22EdnaJones.txt
RandomGoodPoly := proc(n) local ra, xs, ys, aveX, aveY, j:
    ra := rand(-10.0..10.0);
    xs := [seq(ra(), j=1..n)];
    ys := [seq(ra(), j=1..n)];
    aveX := convert(xs, `+`)/n;
    aveY := convert(ys, `+`)/n;
    return [seq([xs[j]-aveX, ys[j]-aveY], j=1..n), [xs[1]-aveX, ys[1]-aveY]];
end:

AvePolygon := proc(P) local midP, j:
    midP := [];
    for j from 1 to nops(P)-1 do
        midP := [op(midP), [(P[j][1]+P[j+1][1])/2, (P[j][2]+P[j+1][2])/2]];
    od:
    return [op(midP), [(P[1][1]+P[2][1])/2, (P[1][2]+P[2][2])/2]];
end:

IterPolygon := proc(n,N) local P, j:
    P := RandomGoodPoly(n);
    for j from 1 to N do
        P := AvePolygon(P);
    od:
    plot(P);
end:
###

# EccOr(P): inputs a good polygon and (hopefully in the future) outputs the eccentricity and orientation of the limiting ellipse
EccOr := proc(P) local n, evalsevecs, evals, evecs, coefs:
    n := nops(P)-1;
    evalsevecs := EVC(evalf(A(n)));
    evals := [seq(evalsevecs[j][1], j=1..n)];
    evecs := [seq(evalsevecs[j][2], j=1..n)];
    coefs := Coeffs([seq(P[j][1]+I*P[j][2], j=1..n)], evecs);
    # I'm not exactly sure what to do from here on
end:


# 5. [Challenge, 10 dollars to be divided] Using experimentation, or otherwise, find an explicit expression, in terms of w:=exp(2*Pi*I/n), for the two largest eigenvalues of A(n) (they are complex-conjugates of each other) and if possible for the corresponding eigenvectors (they are also complex conjugates of each other).

# I claim that the two 2nd largest eigenvalues of A(n) are (1+exp(2*Pi*I/n))/2 = (1+w)/2 and (1+exp(-2*Pi*I/n))/2 = (1+w^(-1))/2 = (1+w^(n-1))/2.

# First of all, I claim that the set of all of the eigenvalues of A(n) is {(1+w^k)/2 : k=1..n} with a corresponding unit eigenvector of [seq(w^(j*k)/sqrt(n), j=1..n)]^T for (1+w^k)/2, where x^T is the transpose of the vector x. Note that each (1+w^k)/2, k = 1,...,n, is distinct since w is a primitive nth root of unity.
# Let k be an integer such that 1<=k<=n. Note that A(n)*[seq(w^(j*k)/sqrt(n), j=1..n)]^T = [seq((w^(j*k)+w^((j+1)*k))/sqrt(n)/2, j=1..n)]^T since w^n = 1 and so w^((n+1)*k) = w^(n*k+k) = w^(n*k)*w^k = w^k. Observe that [seq((w^(j*k)+w^((j+1)*k))/sqrt(n)/2, j=1..n)]^T = [seq((w^(j*k)+w^(j*k+k))/sqrt(n)/2, j=1..n)]^T = [seq(w^(j*k)*(1+w^k)/sqrt(n)/2, j=1..n)]^T = (1+w^k)/2*[seq(w^(j*k)/sqrt(n), j=1..n)]^T, so (1+w^k)/2 is an eigenvalue with a unit eigenvector of [seq(w^(j*k)/sqrt(n), j=1..n)]^T.

# Now that we have found all of the eigenvalues, we can determine the eigenvalues with the largest magnitude that is not 1. Note that the complex conjugate of w is w^(-1) = w^(n-1), so the complex conjugate of w^k is w^(-k) = w^(n-k).
# Let k be an integer such that 1<=k<=n. The magnitude of (1+w^k)/2 is sqrt((1+w^k)/2*conjugate((1+w^k)/2)) = sqrt((1+w^k)/2 * (1+w^(-k))/2) = sqrt(1 + w^k + w^(-k) + 1)/2 = sqrt(2 + 2*Re(w^k))/2. Now sqrt(2 + 2*Re(w^k))/2 = sqrt(2 + 2*cos(2*Pi*k/n))/2, which becomes larger as k approaches a multiple n (i.e., the local maxima all occur at multiples of n). Now sqrt(2 + 2*cos(2*Pi*n/n))/2 = sqrt(2 + 2*cos(2*Pi))/2 = sqrt(2 + 2)/2 = 1. The next closest k to multiples of n are k=1 and k=n-1, so the k that give the 2nd largest eigenvalues are k=1 and k=n-1. Thus, the two 2nd largest eigenvalues of A(n) are (1+w)/2 and (1+w^(n-1))/2.