#hw5.txt. February 8, 2014. Katie McKeon Help5:=proc() print(`RPower(A,D1) Onsm(m,z) NumerOns(m,reso)`): end: read(`C5.txt`): #RPower(A,D1) estimates the dominant eigenvalue of matrix A with error #less than 10^-D1 using the power method RPower:=proc(A,D1) local k, x0,x: x0:=[seq([rand(1..100)()],k=1..nops(A))]: for k from 1 to 10000 do x:=RPow(A,k,x0): if abs(x-RPow(A,k+10,x0)) < 10^(-D1-2) then RETURN(x): fi: od: RETURN(FAIL): end: #Onsm(m,z) gives 2^m x 2^m Onsager matrix in variable z Onsm:=proc(m,z) local M, u,v,i,j,k: M:=[[0$2^m]$2^m]: for i from 1 to 2^m do u:=convert(i,base,2): u:=[0$(m-nops(u)),seq(u[k],k=1..nops(u))]: u:=subs(0=-1,u): for j from 1 to 2^m do v:=convert(j,base,2): v:=[0$(m-nops(v)),seq(v[k],k=1..nops(v))]: v:=subs(0=-1,v): M[i][j]:=z^(sum(u[k]*v[k], k=1..m)+sum(u[k]*u[k+1], k=1..(m-1))+u[m]*u[1]): od: od: M: end: #NumerOns(m,reso) gives the list fm(z) for z from .2 to 3 with increment reso #where fm(z) = log(dominant eisenvalue of Onsm(m,z))/m NumerOns:=proc(m,reso) local z, fmz: with(LinearAlgebra): z:=.2: fmz:=[]: while z<3 do fmz:=[op(fmz),log(Eigenvalues(convert(Onsm(m,z),matrix))[1])/m]: z:=z+reso: od: fmz: end: #If m is larger than 6, the size of the Onsager matrix is too large to be supported by the list structure #It appears that as m increases, there is a minimum around z=1 and a singularity at z=0