Advanced Search

Question: Dos and don'ts of Maple programming w.r.t. the physics package

Question:Dos and don'ts of Maple programming w.r.t. the physics package

NOh 55 Maple 18
physics coeff commutator
October 27 2014
2

EDIT 2: Ok, here is what I am trying to do with Maple.

I have the following expression:

with

 

and

 (j times)

Now the above expression should be zero. I get a set of equations in orders of lambda. The H matrices are known. I only need to find the S matrices.

Anyways, the HOffDiag is the sum above and the GetAdvCommutator0 is this nested commutator. And I need this sum evaluated up to order n (n is a given input-number) in lambda. The rest is really not that performance hungry.

 

 

EDIT: Ok, it realy seems that this sum of commutators is something that takes Maple forever. I have written a new MWE. Is it possible to speed things up? Maybe using assumptions or so?

 

> 

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

[noncommutativeprefix = {H, M, S}]

 (1)
> 

GetAdvCommutator0:=proc(power::integer,order::integer,M)
    if power <= 0 then return M fi:
    if power = 1 then

        return mtaylor(Commutator(M,add(lambda^l*S[l],l=1..order)),lambda, order+1) fi:
    
    return mtaylor(Commutator(GetAdvCommutator0(power-1,order,M),add(lambda^l*S[l],l=1..order)),lambda,order+1):
end proc:
> 

HOffDiag:=n->eval(add(1/((2*j+1)!)*'GetAdvCommutator0'(2*j+1,n,H[0]),j=0..iquo(n-1,2)+1)+add(1/((2*j+1)!)*'GetAdvCommutator0'(2*j+1,n,lambda*H[1]),j=0..iquo(n-2,2)+1)+add(1/((2*j)!)*'GetAdvCommutator0'(2*j,n,lambda*H[2]),j=0..iquo(n-1,2)+1)):
> 

HOffDiag(3); # still fast

lambda^3*Physics:-Commutator(H[0], S[3])+lambda^2*Physics:-Commutator(H[0], S[2])+lambda*Physics:-Commutator(H[0], S[1])+(1/6)*lambda^3*Physics:-Commutator(Physics:-Commutator(Physics:-Commutator(H[0], S[1]), S[1]), S[1])+lambda^2*Physics:-Commutator(H[1], S[1])+lambda^3*Physics:-Commutator(H[1], S[2])+lambda*H[2]+(1/2)*lambda^3*Physics:-Commutator(Physics:-Commutator(H[2], S[1]), S[1])

 (2)
> 

HOffDiag(6); # takes a very long time, I didn't even try entering higher numbers

Warning,  computation interrupted

  
> 

 

 

SumOfCommutators.mw

 

Hi,

I know this question is really unspecific. However, since I am still new to Maple, and since you don't program with Maple as you would do with c++, I would appreciate some advise. I have written a program in Maple and I don't know if it works correctly. All I know is, it is pretty slow :(.

So if anyone of you could just have a short look at my attached worksheet and could tell me, if I at some point used an obvious "don't", I would greatly appreciate it.

 

Cheers

NOh

 

 

A Maple worksheet

NULL

Setup Maple
 

NULL

> 

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

[noncommutativeprefix = {H, MA, MB, S}]

 (1.1)

NULL

NULL

Input parameter
 

NULL

> 

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

H0 := Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = U, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = U})

 (2.1)
> 

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

H1 := Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0})

 (2.2)
> 

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

H2 := Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = -t, (1, 4) = -t, (2, 1) = 0, (2, 2) = 0, (2, 3) = t, (2, 4) = t, (3, 1) = -t, (3, 2) = t, (3, 3) = 0, (3, 4) = 0, (4, 1) = -t, (4, 2) = t, (4, 3) = 0, (4, 4) = 0})

 (2.3)
> 

orderNumber:= 4 ;

4

 (2.4)

NULL

 

 

 

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(power-1,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(power-1,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(power-1,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(n-1,2)+1)+add(1/((2*j+1)!)*'GetAdvCommutator1'(2*j+1,n),j=0..iquo(n-2,2)+1)+add(1/((2*j)!)*'GetAdvCommutator2'(2*j,n),j=0..iquo(n-1,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(n-1,2)+1)+add(1/((2*j+1)!)*'GetAdvCommutator2'(2*j+1,n),j=0..iquo(n-2,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 i-1 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 i-1 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:
> 

 

 

NULL

> 

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

Matrix(4, 4, {(1, 1) = -2*t^2/U, (1, 2) = 2*t^2/U, (1, 3) = 0, (1, 4) = 0, (2, 1) = 2*t^2/U, (2, 2) = -2*t^2/U, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = U+2*t^2/U, (3, 4) = 2*t^2/U, (4, 1) = 0, (4, 2) = 0, (4, 3) = 2*t^2/U, (4, 4) = U+2*t^2/U})

 (1)
> 

Eigenvalues(%);

Vector(4, {(1) = 0, (2) = U, (3) = -4*t^2/U, (4) = (U^2+4*t^2)/U})

 (2)
> 

 

Testing
 

``

#Fast

> 

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

Matrix(4, 4, {(1, 1) = -2*t^2/U, (1, 2) = 2*t^2/U, (1, 3) = 0, (1, 4) = 0, (2, 1) = 2*t^2/U, (2, 2) = -2*t^2/U, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = U+2*t^2/U, (3, 4) = 2*t^2/U, (4, 1) = 0, (4, 2) = 0, (4, 3) = 2*t^2/U, (4, 4) = U+2*t^2/U})

 (4.1)

#Not so fast

> 

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

Matrix(4, 4, {(1, 1) = -2*t^2/U-(8/3)*t^4/U^3, (1, 2) = 2*t^2/U+(8/3)*t^4/U^3, (1, 3) = 0, (1, 4) = 0, (2, 1) = 2*t^2/U+(8/3)*t^4/U^3, (2, 2) = -2*t^2/U-(8/3)*t^4/U^3, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = U+2*t^2/U+(8/3)*t^4/U^3, (3, 4) = 2*t^2/U+(8/3)*t^4/U^3, (4, 1) = 0, (4, 2) = 0, (4, 3) = 2*t^2/U+(8/3)*t^4/U^3, (4, 4) = U+2*t^2/U+(8/3)*t^4/U^3})

 (4.2)

#Already pretty slow :(

> 

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

Warning,  computation interrupted

  
> 

 

``

 

Download MapleWorksheet.mw
Please Wait...