#hw5.txt, 2/9/14, Anthony Zaleski

Help:=proc():
Help1();
print('RPower(A,D1),Onsm(m,z),NumerOns(m,reso)');
end:

##########Problem 1: Rayleigh Power Method##########

#RPower(A,D1) inputs square matrix A, uses Rayleigh power method
#to get an estimate of the dominant eigenvalue with expected 
#error < 10^-D1.

RPower:=proc(A,D1) local x0,k,lambda0,lambda1,err:

x0:=RandMat(nops(A),1,5):

k:=1:
lambda0:=RPow(A,k,x0):
err:=infinity:

while err >= 10^(-(D1+2)) do
	
	k:=k+10:
	lambda1:=RPow(A,k,x0):
	err:=abs(lambda1-lambda0):
	lambda0:=lambda1

od:

lambda1:

end:

#Example application of RPower:
(*
A:=evalf(RandMat(3,3,4));
	A := [[-4., 1., -1.], [-4., 2., -3.], [4., 0., 4.]]
max(abs(Eigenvalues(A)));
	2.96736514139603
RPower(A,5);
	-2.967365142
*)


##########Problem 2: Onsager Matrix##########


#Onsm(m,z) outputs 2^mx2^m Onsager matrix in z.


Onsm:=proc(m,z) local n,M,i,j,k,u,v,pow:

n:=2^m:
M:=[[0$2^m]$2^m]:

for i from 0 to n-1 do
	
	u:=convert(i,base,2):
	u:=[op(u), 0$(m-nops(u))]: 
	u:=2*u-[1$(m)]:

	for j from 0 to n-1 do
		
		v:=convert(j,base,2):
		v:=[op(v), 0$(m-nops(v))]:
		v:=2*v-[1$(m)]:

		pow:=add(u[k]*v[k],k=1..m)+add(u[k]*u[k+1],k=1..m-1)+u[m]*u[1]:
		M[i+1][j+1]:=z^pow:
	od:
od:

M:

end: 

##########Problem 3: Phase Transition##########

#NumerOns(m,reso) outputs f_m(z)=log(largest eigenvalue of Onsm(m,z))/m
#for z between .2 and 3 in increments of reso.

NumerOns:=proc(m,reso) local M,z,Z,x,n,fm,i,MM:

M:=Onsm(m,z):

Z:=[seq(x,x=.2..3,reso)]:
n:=nops(Z):

fm:=[0$n]:

for i from 1 to n do 
	MM:=evalf(subs(z=Z[i],M)):
	fm[i]:=log(abs(RPower(MM,4)))/m:
od:

fm:

end:

###########RESULTS##########

#Also see f_2.png for a plot of z vs f_2(z).

#m  value of z where extremum occurs
#2  1
#3  1
#4  1




##########C5.txt####################

#C5.txt:
#Feb. 6, 2014: Raleigh's power method etc.
Help1:=proc(): print(` RPow(A,k,x0), Ons1(z) `): end:



#####stuff from Monday's "class"
#C4.txt: Feb. 3 2014. Distance-learning class (due to University closing due to snow storm)
Help4:=proc(): print(` MatAdd(A,B),MatMult(A,B), PadZeros(A), ,MatMultR(A,B) `) : 
print(`Tra(A), RandMat(m,n,R)`):
end:


#Tra(A): the transpose of the matrix A
Tra:=proc(A) local i,j:
[seq([seq(A[i][j],i=1..nops(A))],j=1..nops(A[1]))]:
end:


#RandMat(m,n,R): a random m by n matrix with entries between -R and R
#given as a list of lists
RandMat:=proc(m,n,R) local ra,i,j:
ra:=rand(-R..R):
[seq([seq(ra(),j=1..n)],i=1..m)]:
end:


#MatAdd(A,B): adds two matrices of the same size given as lists of lists
MatAdd:=proc(A,B) local i,j:

if nops(A)<>nops(B) then
 print(`Not the same number of rows`):
 RETURN(FAIL):
fi:

if nops({seq(nops(A[i]),i=1..nops(A))})<>1 then
 print(`Rows do not all the same lengths, hence`, A, `is not a matrix `):
 RETURN(FAIL):
fi:


if nops({seq(nops(B[i]),i=1..nops(B))})<>1 then
 print(`Rows do not all the same lengths, hence`, B, `is not a matrix `):
 RETURN(FAIL):
fi:

if nops(A)>0 and nops(A[1])<>nops(B[1]) then
 print(`Not the same number of columns`):
 RETURN(FAIL):
fi:


[seq([seq(A[i][j]+B[i][j],j=1..nops(A[i]))],i=1..nops(A))]

end:


#MatMult(A,B): multiplies two matrices A and B 
MatMult:=proc(A,B) local i,j,k:


if nops({seq(nops(A[i]),i=1..nops(A))})<>1 then
 print(`Rows do not all the same lengths, hence`, A, `is not a matrix `):
 RETURN(FAIL):
fi:


if nops({seq(nops(B[i]),i=1..nops(B))})<>1 then
 print(`Rows do not all the same lengths, hence`, B, `is not a matrix `):
 RETURN(FAIL):
fi:

if nops(A[1])<>nops(B) then
 print(`number of columns of`, A, `is not the same as the number of rows of`, B):
 print(`hence can't multiply them.`):
 RETURN(FAIL):
fi:


[seq([seq( add(A[i][k]*B[k][j],k=1..nops(A[i])), 
j=1..nops(B[1]))],i=1..nops(A))]:

end:

#PadZeros(A): given  a square matrix A pads it with zero to make its
#dimension a powr of 2
PadZeros:=proc(A) local i,n,n1:

if nops(A)<>nops(A[1]) then
 print(`not a square matrix`):
 RETURN(FAIL):
fi:

n:=nops(A):

if type(log[2](n),integer) then
A:
else
n1:=2^ceil(log[2](n)):
[seq([op(A[i]),0$(n1-n)],i=1..nops(A)),[0$n1]$(n1-n)]:
fi:

end:


#MatMultR(A,B): multiplies two square matrices A and B whose dimension
#is a power of 2
MatMultR:=proc(A,B) local n,i,A11,A12,A21,A22,B11,B12,B21,B22, 
M1,M2,M3,M4,M5,M6,M7,M8,C11,C12,C21,C22:
n:=nops(A):

if nops(B)<>n then
 RETURN(FAIL):
fi:

if not type(log[2](n),integer) then
     RETURN(FAIL):
fi:

if n=1 then
 RETURN([[A[1][1]*B[1][1]]]):
fi:

A11:=[seq([op(1..n/2,A[i])],i=1..n/2)]:
A12:=[seq([op(n/2+1..n,A[i])],i=1..n/2)]:
A21:=[seq([op(1..n/2,A[i])],i=n/2+1..n)]:
A22:=[seq([op(n/2+1..n,A[i])],i=n/2+1..n)]:



B11:=[seq([op(1..n/2,B[i])],i=1..n/2)]:
B12:=[seq([op(n/2+1..n,B[i])],i=1..n/2)]:
B21:=[seq([op(1..n/2,B[i])],i=n/2+1..n)]:
B22:=[seq([op(n/2+1..n,B[i])],i=n/2+1..n)]:



M1:=MatMultR(A11,B11):
M2:=MatMultR(A12,B21):
M3:=MatMultR(A11,B12):
M4:=MatMultR(A12,B22):


M5:=MatMultR(A21,B11):
M6:=MatMultR(A22,B21):
M7:=MatMultR(A21,B12):
M8:=MatMultR(A22,B22):
 
C11:=MatAdd(M1,M2):
C12:=MatAdd(M3,M4):
C21:=MatAdd(M5,M6):
C22:=MatAdd(M7,M8):

[seq([op(C11[i]),op(C12[i])],i=1..nops(C11)),
seq([op(C21[i]),op(C22[i])],i=1..nops(C21))]:

end:


#####end stuff from Monday's "class"


#RPow(A,k,x0): Raleigh quotient Power method for approximating
#the dominant eigenvalue of a matrix A (given as a list of lists)
RPow:=proc(A,k,x0) local x1,n,i,j,ma:

n:=nops(A):



 x1:=x0:

for i from 1 to k do
x1:=MatMult(A,x1):
ma:=max( seq( abs(x1[j][1]), j=1..n)):
x1:=[ seq([x1[j][1]/ma], j=1..n)]:
od:

MatMult(Tra(x1), MatMult(A,x1))[1][1]/MatMult(Tra(x1),x1)[1][1]:


end:



#Ons1(): the 2 by 2 transfer matrix for the 1 by infinity
 #Ising model columns are [1,-1], rows are [1,-1]
Ons1:=proc(z) 
[[z,1/z],[1/z,z]]:
end: