A Maple worksheet

Setup Maple


> 
restart;with(Physics):with(LinearAlgebra):Setup(mathematicalnotation=true):Setup(noncommutativeprefix={MA,MB,H,S});


(1.1) 


Input parameter


> 
H0:=Matrix([[ 0 , 0 , 0 , 0 ], [ 0 , 0 , 0 , 0 ], [ 0 , 0 , U , 0 ], [ 0 , 0 , 0 , U ]]);


(2.1) 
> 
H1:=Matrix([[ 0 , 0 , 0 , 0 ], [ 0 , 0 , 0 , 0 ], [ 0 , 0 , 0 , 0 ], [ 0 , 0 , 0 , 0 ]]);


(2.2) 
> 
H2:=Matrix([[ 0 , 0 , t , t ], [ 0 , 0 , t , t ], [ t , t , 0 , 0 ], [ t , t , 0 , 0 ]]);


(2.3) 

(2.4) 


Define all necessary functions


> 
GetAdvCommutator0:=proc(power::integer,order::integer) local Comm: if power <= 0 then return H[0] fi: if power = 1 then return add(coeff(Commutator(H[0],add(lambda^l*S[l],l=1..order)),lambda,n)*lambda^n,n=0..order) fi: return add(coeff(Commutator(GetAdvCommutator0(power1,order),add(lambda^l*S[l],l=1..order)),lambda,n)*lambda^n,n=0..order): end proc:

> 
GetAdvCommutator1:=proc(power::integer,order::integer) local Comm: if power <= 0 then return lambda*H[1] fi: if power = 1 then return add(coeff(Commutator(lambda*H[1],add(lambda^l*S[l],l=1..order)),lambda,n)*lambda^n,n=0..order) fi: return add(coeff(Commutator(GetAdvCommutator1(power1,order),add(lambda^l*S[l],l=1..order)),lambda,n)*lambda^n,n=0..order): end proc:

> 
GetAdvCommutator2:=proc(power::integer,order::integer) local Comm: if power <= 0 then return lambda*H[2] fi: if power = 1 then return add(coeff(Commutator(lambda*H[2],add(lambda^l*S[l],l=1..order)),lambda,n)*lambda^n,n=0..order) fi: return add(coeff(Commutator(GetAdvCommutator2(power1,order),add(lambda^l*S[l],l=1..order)),lambda,n)*lambda^n,n=0..order): end proc:

> 
HOffDiag:=n>eval(add(1/((2*j+1)!)*'GetAdvCommutator0'(2*j+1,n),j=0..iquo(n1,2)+1)+add(1/((2*j+1)!)*'GetAdvCommutator1'(2*j+1,n),j=0..iquo(n2,2)+1)+add(1/((2*j)!)*'GetAdvCommutator2'(2*j,n),j=0..iquo(n1,2)+1)):

> 
HOnDiag:=n>eval(add(1/((2*j)!)*'GetAdvCommutator0'(2*j,n),j=0..iquo(n,2)+1)+add(1/((2*j)!)*'GetAdvCommutator1'(2*j,n),j=0..iquo(n1,2)+1)+add(1/((2*j+1)!)*'GetAdvCommutator2'(2*j+1,n),j=0..iquo(n2,2)+1)):

> 
computeHBlockDiag:=proc( n::posint,firstBlockWidth::posint, H0::Matrix, H1::Matrix, H2::Matrix ) local expr, Orders, i, eq, j, eq2, orders,smatrices,rows::posint,evallist,m,l,cols,Hdiag,ff; expr:=HOffDiag(n); # !!!!! this one is one of the slowest parts I think Orders:=[$1..n]; for i from 1 to n do eq:=add(coeff(expr,lambda,l)*lambda^l,l=0..i)=0; for j from 1 to i1 do eq:=eval(eq,Commutator(H[0],S[j])=Orders[j]); od; eq:=subs(Commutator(H[0],S[i])=MX,eq); eq2:=solve(eq,MX); Orders[i]:=eq2; od; rows,cols:=Dimension(H0); smatrices:=[seq(H0,i=1..n)]; for i from 1 to n do evallist:=[`*`=`.`,H[1]=H1,H[2]=H2]; for j from 1 to i1 do evallist:=[op(evallist),S[j]=smatrices[j]]; od; smatrices[i]:=eval(Orders[i],evallist); for m from 1 to firstBlockWidth do for l from firstBlockWidth+1 to rows do smatrices[i][m,l]:=smatrices[i][m,l]*1/(H0[m][m]H0[l][l]); smatrices[i][l,m]:=smatrices[i][l,m]*1/(H0[m][m]H0[l][l]); od; od; od; expr:=HOnDiag(n); # !!!!! and this one of course (almost identical to the one above) Hdiag:=eval(add(coeff(expr,lambda,l)*lambda^l,l=0..n),lambda=1); evallist:= [op(evallist), H[0]=H0,S[n]=smatrices[n]]; Hdiag:=eval(Hdiag,evallist); return Hdiag; end proc:


> 
computeHBlockDiag(2,2,H0,H1,H2);


(1) 

(2) 

Testing


#Fast
> 
computeHBlockDiag(2,2,H0,H1,H2);


(4.1) 
#Not so fast
> 
computeHBlockDiag(4,2,H0,H1,H2);


(4.2) 
#Already pretty slow :(
> 
computeHBlockDiag(6,2,H0,H1,H2);


