NOh

Dos and don'ts of Maple programming w.r....

Asked: NOh 55 Product: Maple 18

October 27 2014

1 4

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});

(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/((2j+1)!)'GetAdvCommutator0'(2j+1,n,H[0]),j=0..iquo(n-1,2)+1)+add(1/((2j+1)!)'GetAdvCommutator0'(2j+1,n,lambdaH[1]),j=0..iquo(n-2,2)+1)+add(1/((2j)!)'GetAdvCommutator0'(2j,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

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 ]]);

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 ;

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

>	HOnDiag:=n->eval(add(1/((2j)!)'GetAdvCommutator0'(2j,n),j=0..iquo(n,2)+1)+add(1/((2j)!)'GetAdvCommutator1'(2j,n),j=0..iquo(n-1,2)+1)+add(1/((2j+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:

>

>	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

Short-code way to obtain the 5 lowest or...

Asked: NOh 55 Product: Maple 18

October 23 2014

1 1

Hi all,

I have an expression of the form $Sum(a[l]*x^l,l=0..n)$ .

Is there a shorter way to obtain let's say the 5 lowest orders than $add(coeff(expression,x,l)*x^l,l=0..4)$ ?

Cheers

NOh

Insert matrices into expressions of quan...

Asked: NOh 55 Product: Maple 18

October 17 2014

1 1

Hi again,

I am lost... again. I worked with the Physics package and have obtained my result in form of a lengthy equation of commutators of non-commuting objects. So basically my problem is solved, except that I now want to compute some examples to see, if everything is correct. For that purpose I want to replace now the non-commuting objects with explicit matrices, compute the result and see if it agrees with numerical studies. However, I have no clue which command to use, because if I try it with subs, then the multiplication is of the wrong type and Maple won't multiply my matrices. I tried to use eval but wasn't successful.

Here is an exampe to see what I mean.

Test2.mw

Thank you for your time

NOh

trivial Commutator expression is not sim...

Asked: NOh 55 Product: Maple 18

October 17 2014

1 1

The following MWE shows what I mean:

with(Physics):Setup(mathematical=true):

Setup(noncommutativeprefix={MX,MY,MZ});

test:=proc()

local eq;

eq:=-Commutator(MX,MY)-Commutator(MZ,MY);

eq:=simplify(subs(MX=-MZ,eq));

return eq;

end proc:

test(); # yields -[-MZ,MY] - [MZ,MY]

% # yields 0

Any ideas how I can solve this? I would like to return the simplified version.

function call within sum, summation inde...

Asked: NOh 55 Product: Maple 18

October 16 2014

2 4

Hi guys

I have an if statement inside a function. This function then get's called in a sum. However, the result is not correct as the if statement checks wether the argument equals 0, but since the sum passes only the summation index name, it doesn't work as intended.

How can I tell Maple, that it has to evaluate the index before passing it on? MWE appended.

Cheers

NOh