MaplePrimes Posts

MaplePrimes Posts are for sharing your experiences, techniques and opinions about Maple, MapleSim and related products, as well as general interests in math and computing.

Latest Post
  • Latest Posts Feed
  • Version 2 do not enable to expand multiple product like A*A*B*E
    Version 3 will now do that.
    I just forgot to add this feature.

    LL.
     

    ######################################################################
    # NOTICE                                                             #
    # Author: Louis Lamarche                                             #
    #         Institute of Research of Hydro-Quebec (IREQ)               #
    #         Science des données et haute performance                   #
    #         2018, March 7                                              #
    #                                                                    #
    # Function name: ExpandQop (x)                                       #
    #       Purpose: Compute the tensor product of two quantum           #
    #                operators in Dirac notations                        #
    #      Argument: x: a quantum operator                               #
    #  Improvements: Manage all +, -, *, /, ^, mod  operations           #
    #                in the argument. Manages multiple tensor products   #
    #                like A*B*C*F                                        #
    #       Version: 3                                                   #
    #                                                                    #
    #  Copyrigth(c) Hydro-Quebec.                                        #
    #        Note 1: Permission to use this softwate is granted if you   #
    #                acknowledge its author and copyright                #
    #        Note 2: Permission to copy this softwate is granted if you  #
    #                leave this 21 lines notice intact. Thank you.       #
    ######################################################################
    restart;

    with(Physics):
    interface(imaginaryunit=i):
    Setup(mathematicalnotation=true);

    [mathematicalnotation = true]

    (1)

    Setup(quantumoperators={A,B,C,Cn});
    Setup(noncommutativeprefix={a,b,q});

    [quantumoperators = {A, B, C, Cn}]

     

    [noncommutativeprefix = {a, b, q}]

    (2)

    opexp:= op(0,10^x):            # exponentiation id
    opnp := op(0,10*x):            # normal product id
    optp := op(0,Ket(q0)*Ket(q1)): # tensor product id
    opdiv:= `Fraction`:            # fraction       id          
    opsym:= op(0,x):               # symbol         id
    opint:= op(0,10):              # integer        id
    opflt:= op(0,10.0):            # float          id
    opcpx:= op(0,i):               # complex        id
    opbra:= op(0,Bra(q)):          # bra            id
    opket:= op(0,Ket(q)):          # ket            id
    opmod:= op(0, a mod b):        # mod            id
    ExpandQop:=proc(x)
        local nx,ret,j,lkb,cbk,rkb,no,lop,success;
        lop:=op(0,x);
        no:=nops(x);
        if lop = opsym or lop = opint or lop = opflt or
           lop = opbra or lop = opket or lop = opcpx then
             ret:=x;
        else
        if lop = opexp then
            ret:=x;
        else       
        if lop = opnp then
            ret:=1;
            for j from 1 to no do
                ret:=ret*ExpandQop(op(j,x));
            end do;        
        else
        if lop = `+` then
            ret:=0;
            for j from 1 to no do
                ret:=ret+ExpandQop(op(j,x));
            end do;
        else
        if lop = `-` then
            ret:=0;
            for j from 1 to no do
                ret:=ret-ExpandQop(op(j,x));
            end do;
        else
        if lop = opdiv then
           ret:=1;
           for j from 1 to no do
               ret:=ret/ExpandQop(op(j,x));
           end do;
        else
        if lop = opmod then
           ret:=x;
        else
        if lop = optp then
           if (no > 3 ) then
               success:=false;
               nx:=x;
               while not success do
                 lkb:=0; cbk:=0; rkb:=0;ret:=1;
                 for j from 1 to no do
                     if (j>1) then
                          if(lkb=0) then
                              if( type(op(j-1,nx),Ket) and
                                  type(op(j,nx),Bra) ) then lkb:=j-1; fi;
                          else
                              if( type(op(j-1,nx),Ket) and
                                  type(op(j,nx),Bra) ) then rkb:=j;   fi;
                          fi;
                          if( type(op(j-1,nx),Bra) and type(op(j,nx),Ket) )
                                                       then cbk:=j;   fi;
                     fi;
                 end do;
                 if ( (lkb < cbk) and (cbk<rkb) ) then
                     for j from 1     to lkb   do ret := ret*op(j,nx); end do;
                     for j from cbk   to no    do ret := ret*op(j,nx); end do;
                     for j from lkb+1 to cbk-1 do ret := ret*op(j,nx); end do;
                 else
                   ret:=nx;
                 fi;
                 
                 if nx = ret then
                    success := true;
                 else
                    nx := ret;
                 fi
               end do;
           else
               ret:=x;
           fi;
        else ret:=x;
        fi; # optp
        fi; # opmod
        fi; # opdiv
        fi; # `-`
        fi; # `+`
        fi; # `opnp
        fi; # `opexp`
        fi; # opsym, opint, opflt, opbra, opket, opcpx

        return ret;
    end proc:

    # For saving
    # save opexp,opnp,optp,opdiv,opint,opflt,opcpx,opbra,opket,opmod, ExpandQop,"ExpandQop.m"

    # Let A be an operator in a first Hilbert space of dimension n
    #  using the associated orthonormal ket and bra vectors
    #
    #
    kets1:=Ket(a1)*Ket(a2)*Ket(a3)*Ket(a4)*Ket(a5):
    A:=kets1*Dagger(kets1);


    # Let B be an operator in a second Hilbert (Ketspace of dimension m
    #  using the associated orthonormal ket and bra vectors
    #
    #
    kets2:=Ket(b1)*Ket(b2)*Ket(b3):
    B:=kets2*Dagger(kets2);


    # The tensor product of the two operators acts on a n+m third
    # Hilbert space   unsing the appropriately ordered ket
    # and bra  vectors of the two preceding spaces. The rule for
    # building this operator in Dirac notation is as follows,
    #
    #


    print("Maple do not compute the tensor product of operators,");
    print("C=A*B gives:");
    C:=A*B;

    print("ExpandQop(C) gives the expected result:");
    Cn:=ExpandQop(C);

    Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

     

    Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

     

    "Maple do not compute the tensor product of operators,"

     

    "C=A*B gives:"

     

    Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

     

    "ExpandQop(C) gives the expected result:"

     

    Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

    (3)

    kets3:=kets1*kets2;
    bras3:=Dagger(kets3);
    print("Matrix elements computed with C appears curious");
    'bras3.C. kets3'="...";
    bras3.C.kets3;
    print("Matrix elements computed with Cn as expected");
    'bras3.Cn.kets3'=bras3.Cn.kets3;

    Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3))

     

    Physics:-`*`(Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

     

    "Matrix elements computed with C appears curious"

     

    Physics:-`.`(bras3, C, kets3) = "..."

     

    Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(b3))*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))

     

    "Matrix elements computed with Cn as expected"

     

    Physics:-`.`(bras3, Cn, kets3) = Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))^2*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))^2*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))^2*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))^2*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))^2*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))^2

    (4)

    print("Example");
    En:=ExpandQop(10*(1-x+y+z)*i*(1/sqrt(2))*A*B);

    "Example"

     

    -(5*I)*2^(1/2)*(-1+x-y-z)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

    (5)

    print("Another example");
    'F'='A*B/sqrt(2)+B*A/sqrt(2)';
    F:=A*B/sqrt(2)+B*A/sqrt(2):
    'op(1,F)'=op(1,F);
    'op(2,F)'=op(2,F);

    'Fn'='ExpandQop(F)';
    Fn:=ExpandQop(F):
    'op(1,Fn)'=op(1,Fn);
    'op(2,Fn)'=op(2,Fn);

    "Another example"

     

    F = Physics:-`*`(Physics:-`*`(A, B), Physics:-`^`(sqrt(2), -1))+Physics:-`*`(Physics:-`*`(B, A), Physics:-`^`(sqrt(2), -1))

     

    op(1, F) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

     

    op(2, F) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

     

    Fn = ExpandQop(F)

     

    op(1, Fn) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

     

    op(2, Fn) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

    (6)

    print("Final example, multiple products");
    G:=B*B*B;
    'G'=ExpandQop(G);

    "Final example, multiple products"

     

    Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

     

    G = Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

    (7)

     


     

    Download ExpandQopV3.mw

     

    Automatic handling of collision of tensor indices in products

     

     

    The design of products of tensorial expressions that have contracted indices got enhanced. The idea: repeated indices in certain subexpressions are actually dummies. So suppose T[a, b] and B[b] are tensors, then in T[trace] = T[a, `~a`], a is just a dummy, therefore T[a, `~a`]*B[a] = T[b, `~b`]*B[a] is a well defined object. The new design automatically maps input like T[a, `~a`]*B[a] into T[b, `~b`]*B[a]. This significantly simplifies the manipulation of indices when working with products of tensorial expressions: each tensorial expression being multiplied has its repeated indices automatically transformed into different ones when they would collide with the free or repeated indices of the other expressions being multiplied.

     

    This functionality is available within the Physics update distributed at the Maplesoft R&D Physics webpage (but for what you see under Preview that makes use of a new kernel feature of the Maple version under development).

     

    restart

    with(Physics); Setup(spacetimeindices = lowercaselatin, quiet)

    [spacetimeindices = lowercaselatin]

    (1)

    Define(T[a, b], B[b])

    `Defined objects with tensor properties`

     

    {B[b], Physics:-Dgamma[a], Physics:-Psigma[a], T[a, b], Physics:-d_[a], Physics:-g_[a, b], Physics:-KroneckerDelta[a, b], Physics:-LeviCivita[a, b, c, d]}

    (2)

    This shows the automatic handling of collision of indices

    T[a, a]*B[a]

    T[b, `~b`]*B[a]

    (3)

    T[a, a]^2

    T[a, `~a`]*T[b, `~b`]

    (4)

    ``

    Preview only in the upcomming version under development

     

    Consider now the case of three tensors

    Define(A[a], C[a])

    `Defined objects with tensor properties`

     

    {A[a], B[b], C[a], Physics:-Dgamma[a], Physics:-Psigma[a], T[a, b], Physics:-d_[a], Physics:-g_[a, b], Physics:-KroneckerDelta[a, b], Physics:-LeviCivita[a, b, c, d]}

    (5)

    A[a]*B[a]*C[a]

    A[a]*B[a]*C[a]

    (6)

    The product above has indeed the index a repeated more than once, therefore none of its occurrences got automatically transformed into contravariant in the output, and Check  detects the problem interrupting with an error  message

    Check(A[a]*B[a]*C[a])

    Error, (in Physics:-Check) wrong use of the summation rule for repeated indices: `a repeated 3 times`, in A[a]*B[a]*C[a]

     

     

    However, it is now also possible to indicate, using parenthesis, that the product of two of these tensors actually form a subexpression, so that the following two tensorial expressions are well defined, and the colliding dummy is automatically replaced making that explicit

    A[a]*B[a]*C[a]

    A[b]*B[`~b`]*C[a]

    (7)

    A[a]*B[a]*C[a]

    A[a]*B[b]*C[`~b`]

    (8)

     

     

    This change in design makes concretely simpler the use of indices in that it eliminates the need for manually replacing dummies. For example, consider the tensorial expression for the angular momentum in terms of the coordinates and momentum vectors, in 3 dimensions

     

    Setup(coordinates = cartesian, dimension = 3, metric = euclidean, quiet)

    [coordinatesystems = {X}, dimension = 3, metric = {(1, 1) = 1, (2, 2) = 1, (3, 3) = 1}]

    (9)

    Define L[j], p[k] respectively representing angular and linear momentum

    Define(L[j], p[k])

    `Defined objects with tensor properties`

     

    {Physics:-Dgamma[a], L[j], Physics:-Psigma[a], Physics:-d_[a], Physics:-g_[a, b], p[k], Physics:-KroneckerDelta[a, b], Physics:-LeviCivita[a, b, c], Physics:-SpaceTimeVector[a](X)}

    (10)

    Introduce the tensorial expression for L[a]

    L[a] = LeviCivita[a, b, c]*X[b]*p[c]

    L[a] = Physics:-LeviCivita[a, b, c]*Physics:-SpaceTimeVector[b](X)*p[c]

    (11)

    The left-hand side has one free index, a, while the right-hand side has two dummy indices b and c

    Check(L[a] = Physics[LeviCivita][a, b, c]*Physics[SpaceTimeVector][b](X)*p[c], all)

    `The repeated indices per term are: `[{`...`}, {`...`}, `...`]*`; the free indices are: `*{`...`}

     

    ([{}], {a}) = ([{b, c}], {a})

    (12)

    If we want to compute`#mrow(msup(mfenced(mover(mi("L"),mo("&rarr;")),open = "&Vert;",close = "&Vert;"),mn("2")),mo("&equals;"),msubsup(mi("L"),mi("a"),mn("2")))`we can now take the square of (11) directly, and the dummy indices on the right-hand side are automatically handled, there is now no need to manually substitute the repeated indices to avoid their collision

    (L[a] = Physics[LeviCivita][a, b, c]*Physics[SpaceTimeVector][b](X)*p[c])^2

    L[a]^2 = Physics:-LeviCivita[a, b, c]*Physics:-SpaceTimeVector[b](X)*p[c]*Physics:-LeviCivita[a, d, e]*Physics:-SpaceTimeVector[d](X)*p[e]

    (13)

    The repeated indices on the right-hand side are now a, b, c, d, e

    Check(L[a]^2 = Physics[LeviCivita][a, b, c]*Physics[SpaceTimeVector][b](X)*p[c]*Physics[LeviCivita][a, d, e]*Physics[SpaceTimeVector][d](X)*p[e], all)

    `The repeated indices per term are: `[{`...`}, {`...`}, `...`]*`; the free indices are: `*{`...`}

     

    ([{a}], {}) = ([{a, b, c, d, e}], {})

    (14)

    NULL


     

    Download AutomaticHandlingCollisionOfTensorIndices.mw

    Edgardo S. Cheb-Terrab
    Physics, Differential Equations and Mathematical Functions, Maplesoft

    Here is a major upgrade of the procedure I submitted in february.

    https://www.mapleprimes.com/posts/209030-Procedure-For-Computing-The-Tensor-Product

    There is a line after the procedure to save it in the file "ExpandQop.m"
    In future post I will use it in order to minimize the size of the examples.

    Louis Lamarche
     

    ######################################################################
    # NOTICE                                                             #
    # Author: Louis Lamarche                                             #
    #         Institute of Research of Hydro-Quebec (IREQ)               #
    #         Science des données et haute performance                   #
    #         2018, March 7                                              #
    #                                                                    #
    # Function name: ExpandQop (x)                                       #
    #       Purpose: Compute the tensor product of two quantum           #
    #                operators in Dirac notations                        #
    #      Argument: x: a quantum operator                               #
    #  Improvements: Manage all +, -, *, /, ^, mod  operations           #
    #                in the argument                                     #
    #       Version: 2                                                   #
    #                                                                    #
    #  Copyrigth(c) Hydro-Quebec.                                        #
    #        Note 1: Permission to use this softwate is granted if you   #
    #                acknowledge its author and copyright                #
    #        Note 2: Permission to copy this softwate is granted if you  #
    #                leave this 21 lines notice intact. Thank you.       #
    ######################################################################
    restart;

    with(Physics):
    interface(imaginaryunit=i):
    Setup(mathematicalnotation=true);

    [mathematicalnotation = true]

    (1)

    Setup(quantumoperators={A,B,C,Cn});
    Setup(noncommutativeprefix={a,b,q});

    [quantumoperators = {A, B, C, Cn}]

     

    [noncommutativeprefix = {a, b, q}]

    (2)

    opexp:= op(0,10^x):            # exponentiation id
    opnp := op(0,10*x):            # normal product id
    optp := op(0,Ket(q0)*Ket(q1)): # tensor product id
    opdiv:= `Fraction`:            # fraction       id          
    opsym:= op(0,x):               # symbol         id
    opint:= op(0,10):              # integer        id
    opflt:= op(0,10.0):            # float          id
    opcpx:= op(0,i):               # complex        id
    opbra:= op(0,Bra(q)):          # bra            id
    opket:= op(0,Ket(q)):          # ket            id
    opmod:= op(0, a mod b):        # mod            id
    ExpandQop:=proc(x)
        local ret,j,lkb,cbk,rkb,no,lop;
        lkb:=0; cbk:=0; rkb:=0;
        lop:=op(0,x);
        no:=nops(x);
        if lop = opsym or lop = opint or lop = opflt or
           lop = opbra or lop = opket or lop = opcpx then
             ret:=x;
        else
        if lop = opexp then
            ret:=x;
        else       
        if lop = opnp then
            ret:=1;
            for j from 1 to no do
                ret:=ret*ExpandQop(op(j,x));
            end do;        
        else
        if lop = `+` then
            ret:=0;
            for j from 1 to no do
                ret:=ret+ExpandQop(op(j,x));
            end do;
        else
        if lop = `-` then
            ret:=0;
            for j from 1 to no do
                ret:=ret-ExpandQop(op(j,x));
            end do;
        else
        if lop = opdiv then
           ret:=1;
           for j from 1 to no do
               ret:=ret/ExpandQop(op(j,x));
           end do;
        else
        if lop = opmod then
           ret:=x;
        else
        if lop = optp then
            ret:=1;
           if (no > 3 ) then
               for j from 1 to no do
                   if (j>1) then
                        if(lkb=0) then
                            if( type(op(j-1,x),Ket) and
                                type(op(j,x),Bra) ) then lkb:=j-1; fi;
                        else
                            if( type(op(j-1,x),Ket) and
                                type(op(j,x),Bra) ) then rkb:=j;   fi;
                        fi;
                        if( type(op(j-1,x),Bra) and type(op(j,x),Ket) )
                                                    then cbk:=j;   fi;
                   fi;
               end do;
               if ( (lkb < cbk) and (cbk<rkb) ) then
                   for j from 1     to lkb   do ret := ret*op(j,x); end do;
                   for j from cbk   to no    do ret := ret*op(j,x); end do;
                   for j from lkb+1 to cbk-1 do ret := ret*op(j,x); end do;
               else
                   ret:=x;
               fi;
           else
               ret:=x;
           fi;
        else ret:=x;
        fi; # optp
        fi; # opmod
        fi; # opdiv
        fi; # `-`
        fi; # `+`
        fi; # `opnp
        fi; # `opexp`
        fi; # opsym, opint, opflt, opbra, opket, opcpx

        return ret;
    end proc:

    # For saving
    # save opexp,opnp,optp,opdiv,opint,opflt,opcpx,opbra,opket,opmod, ExpandQop,"ExpandQop.m"

    # Let A be an operator in a first Hilbert space of dimension n
    #  using the associated orthonormal ket and bra vectors
    #
    #
    kets1:=Ket(a1)*Ket(a2)*Ket(a3)*Ket(a4)*Ket(a5):
    A:=kets1*Dagger(kets1);


    # Let B be an operator in a second Hilbert (Ketspace of dimension m
    #  using the associated orthonormal ket and bra vectors
    #
    #
    kets2:=Ket(b1)*Ket(b2)*Ket(b3):
    B:=kets2*Dagger(kets2);


    # The tensor product of the two operators acts on a n+m third
    # Hilbert space   unsing the appropriately ordered ket
    # and bra  vectors of the two preceding spaces. The rule for
    # building this operator in Dirac notation is as follows,
    #
    #


    print("Maple do not compute the tensor product of operators,");
    print("C=A*B gives:");
    C:=A*B;

    print("ExpandQop(C) gives the expected result:");
    Cn:=ExpandQop(C);

    Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

     

    Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

     

    "Maple do not compute the tensor product of operators,"

     

    "C=A*B gives:"

     

    Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

     

    "ExpandQop(C) gives the expected result:"

     

    Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

    (3)

    kets3:=kets1*kets2;
    bras3:=Dagger(kets3);
    print("Matrix elements computed with C appears curious");
    'bras3.C. kets3'="...";
    bras3.C.kets3;
    print("Matrix elements computed with Cn as expected");
    'bras3.Cn.kets3'=bras3.Cn.kets3;

    Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3))

     

    Physics:-`*`(Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

     

    "Matrix elements computed with C appears curious"

     

    Physics:-`.`(bras3, C, kets3) = "..."

     

    Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(b3))*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))

     

    "Matrix elements computed with Cn as expected"

     

    Physics:-`.`(bras3, Cn, kets3) = Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))^2*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))^2*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))^2*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))^2*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))^2*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))^2

    (4)

    print("Example");
    En:=ExpandQop(10*(1-x+y+z)*i*(1/sqrt(2))*A*B);

    "Example"

     

    -(5*I)*2^(1/2)*(-1+x-y-z)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

    (5)

    print("Another example");
    'F'='A*B/sqrt(2)+B*A/sqrt(2)';
    F:=A*B/sqrt(2)+B*A/sqrt(2):
    'op(1,F)'=op(1,F);
    'op(2,F)'=op(2,F);

    'Fn'='ExpandQop(F)';
    Fn:=ExpandQop(F):
    'op(1,Fn)'=op(1,Fn);
    'op(2,Fn)'=op(2,Fn);

    "Another example"

     

    F = Physics:-`*`(Physics:-`*`(A, B), Physics:-`^`(sqrt(2), -1))+Physics:-`*`(Physics:-`*`(B, A), Physics:-`^`(sqrt(2), -1))

     

    op(1, F) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

     

    op(2, F) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

     

    Fn = ExpandQop(F)

     

    op(1, Fn) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

     

    op(2, Fn) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

    (6)

     


     

    Download ExpandQopV2.mw

    --- Prolog.ue ---

    The best things in life come free of charge.

    Happiness, love, and wisdom of expertise are first few that hit my mind.

    As a business economist, I keep my eyes keenly open to opportunities for growth; such as Maple 2017 training session.

    It was a Saturday afternoon in Waterloo, ON, this chilly Feburary which was blessed by snowstorm warning.

     

    --- Encountering with Maple ---

    I was aware of Maple for many years back when my academic career began.

    In fact, Maple was available in the lab computers at university. 

    But I did not know what to do with it.

    Nor did I use any mathematics softwares until recently, but I had this thought : one day I could learn.

    The motivation for this `learn how to use it' did not occur to me for a long time (14 years!!).

    Things changed this year when I enrolled to an Electrical Engineering program at Lassonde.

    Mind you, I have already been using various types of languages and tools such as: Python, C, Java, OpenOfficeSuites, Stata, SAS, Latex just to mention a few.

    These stuffs also run on multiple platforms which I am sure you have heard of if you're reading this post; Windows, OSX and Linux. And Maple supports all these major operating systems.

     

    --- Why do I like Maple ---

    During the first week of school, Dr. Smith would ask us to purchase and practice using MATLAB because it had a relatively easy learning curve for beginners like python and we were going to use it for labs.

    Furthermore, students get a huge discount (i.e. $1500 to $50) for these softwares.

    Then, the professor also added; "Maple is also a great tool to use, but we won't use it for this class".

    ME: ' Why not ? '

    The curiosity inside me gave in and I decided to try both!

    After all, my laziness in taking derivatives by hand or the possibility of making mistake would disappear if I can verify results using software.

    That's it...!

    Being able to check correct answer was already worth more than $50.

    I can not emphasize this point enough; 

    For people in the industry being paid for their time, or students like me who got a busy schedule can not afford to waste any time. (i.e. need to minimize homework effort & frustration, while maximizing the educational attainment & final grades)

    Right? Time is money.

    Don't we all just want more spare time for things we care?

    Googling through many ambiguous Yahoo Answers or online forums like Stackoverflow replies are often misleading and time consuming. 

    I have spent years (estimated 3000+ hours) going through those wildly inaccurate webpages hoping for some clearly written information with sub-optimal outcome.

    Diverting many hours of study time is not something a first year S.T.E.M. students can afford.

     

    --- Maple Training ---

    Now you know about my relationships with Maple; Let me describe how the training session went.

    I will begin with the sad news first, =(

    First of all, there was no coffee available when I arrived. It arrived only after lunch.

    Although it was a free event aside other best things in life, this was only a material factor I didn't enjoy at the site. 

    Still a large portion of Canadians start their work with a zolt of caffeine in my defence.

    Secondly, there was a kind of assumption which expected attendee were familiar with software behavior.

    A handful of people were having trouble opening example file, perhaps because of their browser setting or link to preferred software by OS.

    Not being able to follow the tutorials as the presenter demonstrated various facets of software substantially diminished the  efficacy of training session for those who could not be on the same page.

    These minor annoyances were the only drawbacks I experinced from the event.

     

    Here comes the happy side, =)

    1. The staffs were considerate enough to provide vegetarion options for inclusive lunch as well as answering all my curious, at times orthogonal questions regarding Maplesoft company.

    2. Highly respectable professionals were presenting themselves; 

    That is, Prof. Illias Kotsireas, Dr. Erik Postma and Dr. Jürgen Gerhard.

    I can not appreciate enough of their contribution for the training in an eloquent and humble manners.

    To put it other way, leading of the presentation was well structured and planned out.

    In the beginning, Prof. Kotsireas presented `Introduction to Maple' which included terminology and basic behaviors of Maple (i.e. commands and features) with simple examples you can quickly digest. Furthermore, Maple has internal function to interface with Latex! No more typing hours of $$s and many frac{}{}, \delta_{} to publish. In order for me to study all this would have been two-weeks kind of commitment in which he summarized in a couple of hours time. Short-cut keys that are often used by his project was pretty interesting, which will improve work efficiency.

    After a brief lunch, which was supplied more than enough for all, Dr. Erik Postma delivered a critical component of simluation. That is, `Random Number Generation'. Again, he showed us some software-related tricks such as `Text mode' vs. `Math mode'.  The default RNG embedded in the software allows reproducible results unless we set seed and randomize further. Main part of the presentation was regarding `Optimization of solution through simulation'. He iteratively improved efficiency of test model, which I will not go in depth here. However, visually and quantitatively showing the output was engaging the attendees and Maple has modularized this process (method available for all the users!!).

    Finally, we got some coffee break that allowed to me to push through all the way to the end. I believe if we had some coffee earlier less attendees would have left.

    The last part of the training was presented by Dr. Jürgen Gerhard. In this part, we were using various applications of Maple in solving different types of problems. We tackled combinatorics of Fibonacci sequence by formula manipulation. In this particular example he showed us how to optimize logic of a function that made a huge impact in processing time and memory usage. Followed by graph theory example, damped harmonic oscillator, 2 DOF chaotic system, optimization and lastly proof of orthocentre by coding. I will save the examples for you to enjoy in future sessions. 

    The way they went through examples were super easy to follow. This can only be done with profound understanding of the subject and a lot of prior effort in preparing the presentation.
     

    I appreciate much efforts put together by whom organized this event, allocating their own precious weekend time and allowing many to gain opportunity to learn directly from the person in the house.

     

    --- Epilogue ---

    My hope for Maple usage lies in enhancing education outcome for first year students, especially in the field of Science and Economics. This is a free opportunity for economic empowerment which is uncaptured.

    Engineering students are already pretty good at problem solving, and will figure things out as I witnessed my colleagues have.

    However, students of natural sciences and B.A. programs tend to skimp on utilizing tools due to lack of exposure.

    Furthermore, I am supporting their development of SaaS, software as service, which delivers modules like gRPC does.

    Also, I hope the optimization package from prior version written by Dr. Postma will become available to public sometime.

    Here's a BIG thank you to staffs once again, and forgive me for any grammatical errors from rushed writing. I tried to incorporate as much observation as possible gathered from the event.

    To contact me, my email is hyonwoo.kee (at) gmail.com;

     

    Minimize the number of tensor components according to its symmetries
    (and relabel, redefine or count the number of independent tensor components)

     

     

    The nice development described below is work in collaboration with Pascal Szriftgiser from Laboratoire PhLAM, Université Lille 1, France, used in the Mapleprimes post Magnetic traps in cold-atom physics

     

    A new keyword in Define  and Setup : minimizetensorcomponents, allows for automatically minimizing the number of tensor components taking into account the tensor symmetries. For example, if a tensor with two indices in a 4D spacetime is defined as antisymmetric using Define with this new keyword, the number of different tensor components will be exactly 6, and the elements of the diagonal are automatically set equal to 0. After setting this keyword to true with Setup , all subsequent definitions of tensors automatically minimize the number of components while using this keyword with Define  makes this minimization only happen with the tensors being defined in the call to Define .

     

    Related to this new functionality, 4 new Library routines were added: MinimizeTensorComponents, NumberOfIndependentTensorComponents, RelabelTensorComponents and RedefineTensorComponents

     

    Example:

    restart; with(Physics)

     

    Define an antisymmetric tensor with two indices

    Define(F[mu, nu], antisymmetric)

    `Defined objects with tensor properties`

     

    {Physics:-Dgamma[mu], F[mu, nu], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu]}

    (1.1)

    Although the system knows that F[mu, nu] is antisymmetric, you need to use Simplify to apply the (anti)symmetry

    F[mu, nu]+F[nu, mu]

    F[mu, nu]+F[nu, mu]

    (1.2)

     

    Simplify(F[mu, nu]+F[nu, mu])

    0

    (1.3)

    so by default the components of F[mu, nu] do not automatically reflect the (anti)symmetry; likewise

    F[1, 2]+F[2, 1]

    F[1, 2]+F[2, 1]

    (1.4)

    Simplify(F[1, 2]+F[2, 1])

    0

    (1.5)

    and computing the array form of F[mu, nu]we do not see the elements of the diagonal equal to zero nor the lower-left triangle equal to the upper-right triangle but for a different sign:

    TensorArray(F[mu, nu])

    Matrix(%id = 18446744078270093062)

    (1.6)

     

    On the other hand, this new functionality, here called minimizetensorcomponents, makes the symmetries of the tensor be explicitly reflected in its components.

     

    There are three ways to use it. First, one can minimize the number of tensor components of a tensor previously defined. For example

     

    Library:-MinimizeTensorComponents(F)

    Matrix(%id = 18446744078270064630)

    (1.7)

    After this, both (1.2) and (1.3) are automatically equal to 0 without having to use Simplify

    F[mu, nu]+F[nu, mu]

    0

    (1.8)

    0

    0

    (1.9)

    And the output of TensorArray  in (1.6) becomes equal to (1.7).

     

    NOTE: in addition, after using minimizetensorcomponents in the definition of a tensor, say F, all the keywords implemented for Physics tensors are available for F:

     

    F[]

    F[mu, nu] = Matrix(%id = 18446744078247910206)

    (1.10)

    F[trace]

    0

    (1.11)

    F[nonzero]

    F[mu, nu] = {(1, 2) = F[1, 2], (1, 3) = F[1, 3], (1, 4) = F[1, 4], (2, 1) = -F[1, 2], (2, 3) = F[2, 3], (2, 4) = F[2, 4], (3, 1) = -F[1, 3], (3, 2) = -F[2, 3], (3, 4) = F[3, 4], (4, 1) = -F[1, 4], (4, 2) = -F[2, 4], (4, 3) = -F[3, 4]}

    (1.12)

    "F[~1,mu,matrix]"

    F[`~1`, mu] = Vector[row](%id = 18446744078247885990)

    (1.13)

    Alternatively, one can define a tensor, specifying that the symmetries should be taken into account to minimize the number of its components passing the keyword minimizetensorcomponents to Define .

     

    Example:

     

    Define a tensor with the symmetries of the Riemann  tensor, that is, a tensor of 4 indices that is symmetric with respect to interchanging the positions of the 1st and 2nd pair of indices and antisymmetric with respect to interchanging the position of its 1st and 2nd indices, or 3rd and 4th indices, and define it minimizing the number of tensor components

     

    Define(R[alpha, beta, mu, nu], symmetric = {[[1, 2], [3, 4]]}, antisymmetric = {[1, 2], [3, 4]}, minimizetensorcomponents)

    `Defined objects with tensor properties`

     

    {Physics:-Dgamma[mu], F[mu, nu], Physics:-Psigma[mu], R[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu]}

    (1.14)

    We now have

    R[1, 2, 3, 4]+R[2, 1, 3, 4]

    0

    (1.15)

    R[alpha, beta, mu, nu]-R[mu, nu, alpha, beta]

    0

    (1.16)
    • 

    One can always retrieve the symmetry properties in the abstract notation used by the Define command using the new Library:-GetTensorSymmetryProperties, its output is ordered, first the symmetric then the antisymmetric properties

     

    Library:-GetTensorSymmetryProperties(R)

    {[[1, 2], [3, 4]]}, {[1, 2], [3, 4]}

    (1.17)
    • 

    After making the symmetries explicit (and also before that), it is frequently useful to know the number of independent components of a given tensor. For this purpose you can use the new Library:-NumberOfIndependentTensorComponents

     

    Library:-NumberOfIndependentTensorComponents(R)

    21

    (1.18)

    and besides taking into account the symmetries, in the case of the Riemann  tensor, after taking into account the first Bianchi identity this number of components is further reduced to 20.

     

    A third way of using the new minimizetensorcomponents functionality is using Setup , so that, automatically, every subsequent definition of tensors with symmetries is performed minimizing the number of its components using the indicated symmetries

     

    Example:

    Setup(minimizetensorcomponents = true)

    [minimizetensorcomponents = true]

    (1.19)

    So from hereafter you can define tensors taking into account their symmetries explicitly and without having to include the keyword minimizetensorcomponents at each definition

     

    Define(C[alpha, beta], antisymmetric)

    `Defined objects with tensor properties`

     

    {C[mu, nu], Physics:-Dgamma[mu], F[mu, nu], Physics:-Psigma[mu], R[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu]}

    (1.20)

     

    C[]

    C[mu, nu] = Matrix(%id = 18446744078408747598)

    (1.21)
    • 

    Two new related functionalities are provided via Library:-RelabelTensorComponents and Library:-RedefineTensorComponent, the first one to have the number of tensor components directly reflected in the names of the components, the second one to redefine only one of these components

    Library:-RelabelTensorComponents(C)

    Matrix(%id = 18446744078408729774)

    (1.22)

     

    Suppose now we want to make one of these components equal to 1, say C__2

    Library:-RedefineTensorComponent(C[1, 2] = 1)

    C[mu, nu] = Matrix(%id = 18446744078270104390)

    (1.23)

    This nice development is work in collaboration with Pascal Szriftgiser from Laboratoire PhLAM, UMR CNRS 8523, Université Lille 1, F-59655, France.

    ``

     

    Download MinimizeTensorComponents.mw

     

    Edgardo S. Cheb-Terrab
    Physics, Differential Equations and Mathematical Functions, Maplesoft

     

    Hello, 

    I study mainly subjects that fall under umbrella of number theory, but i have specified a little further in the worksheet. This is really a request for assistance, because in as much as i have met so many brilliant people online via social media etc,  I would always love to meet more, and especially ones who are more experienced in this field. 

     

    Basically i am too cheap and old to think about going to a good university, so I am trying to get free advice from the people who have probably completed doctorates in the relevant field. Got to be honest I say.

     

    Anyway my contact email is at the top of the attached worksheet.

     

    First thing that stood out to me about the distributions produced in this worksheet is how sparse the number of points is for N=17 relative to all the other values of N.

    EXAMPLE_FOR_MAPLE3.mw

     

    Edit: Another example worksheet added.

    MAPLE_EXAMPLE_13.mw

    Dear all , I' would like to join a group to produce quantum information tools in Maple

     

    I wanted to use MAPLE to preform symbolic quantum computations. The role
    of quantum operators and their tensor product is very important in simplying
    the understating of such new calculus at least for the beginners. For instance,
    (using "o" for the tensor product and "." for the scalar product, H being the Hadamard
    operator on a qubit, I the identity operator, and CNOT the 2 qubit controled not
    operator)
    1) generating the Bells states |Bxy> two stages of operators are needed
         (CNOT) .  (H o  I)  . |x> o |y>

    2) performing quantum teleportation of |psi>
         (H o I o I) . (CNOT o I ) . |psi>o |B00>
        followed by a measurements on the first two qubits for driving the application of
        quantum gates to the third qubit.

    All these tensor products of operators can be easily written in MAPLE.

    Here is a first version of the ExpandQop procedure that will be usefull the purpose of
    expanding correctly the tensor product of two quantum operator expressed in Dirac notation.

    I hope this is usefull.

    LL 

     

    ######################################################################
    # Author: Louis Lamarche                                             #
    #         Institute of Research of Hydro-Quebec (IREQ)               #
    #         Science des données et haute performance                   #
    #         2018/02/20                                                 #
    #                                                                    #
    #         Function name: ExpandQop (x)                               #
    #               Purpose: Compute the tensor product of two quantum   #
    #                        operators in Dirac notations                #
    #              Argument: x - a simple quantum operator               #
    #   Future improvements: Manage all +, -, *, /, ^, mod  operations   #
    #                        in the argument                             #
    #               Version: 1.0                                         #
    ######################################################################
    restart;

    with(Physics):
    interface(imaginaryunit=i):
    Setup(mathematicalnotation=true);

    [mathematicalnotation = true]

    (1)

    Setup(quantumoperators={A,B,C,Cn});
    Setup(noncommutativeprefix={a,b});

    [quantumoperators = {A, B, C, Cn}]

     

    [noncommutativeprefix = {a, b}]

    (2)

    ExpandQop:=proc(x)
        local ret,j,lkb,cbk,rkb,no;
        ret:=1; lkb:=0; cbk:=0; rkb:=0; no:=nops(x);
        if (no > 3 ) then
            for j from 1 to no do
                if (j>1) then
                     if(lkb=0) then
                         if( type(op(j-1,x),Ket) and
                             type(op(j,x),Bra) ) then lkb:=j-1; fi;
                     else
                         if( type(op(j-1,x),Ket) and
                             type(op(j,x),Bra) ) then rkb:=j;   fi;
                     fi;
                     if( type(op(j-1,x),Bra) and type(op(j,x),Ket) )
                                                 then cbk:=j;   fi;
                fi;
            end do;
            if ( (lkb < cbk) and (cbk<rkb) ) then
                for j from 1     to lkb   do ret := ret*op(j,x); end do;
                for j from cbk   to no    do ret := ret*op(j,x); end do;
                for j from lkb+1 to cbk-1 do ret := ret*op(j,x); end do;
            else
                ret:=x;
            fi;
        else
            ret:=x;
        fi;
        return ret;
    end proc:

    # Let A be an operator in a first Hilbert space of dimension n
    #  using the associated orthonormal ket and bra vectors
    #
    #
    kets1:=Ket(a1)*Ket(a2)*Ket(a3)*Ket(a4)*Ket(a5):
    A:=kets1*Dagger(kets1);


    # Let B be an operator in a second Hilbert (Ketspace of dimension m
    #  using the associated orthonormal ket and bra vectors
    #
    #
    kets2:=Ket(b1)*Ket(b2)*Ket(b3):
    B:=kets2*Dagger(kets2);


    # The tensor product of the two operators acts on a n+m third
    # Hilbert space   unsing the appropriately ordered ket
    # and bra  vectors of the two preceding spaces. The rule for
    # building this operator in Dirac notation is as follows,
    #
    #


    print("Maple do not compute the tensor product of operators,");
    print("C=A*B gives:");
    C:=A*B;

    print("ExpandQop(C) gives the expected result:");
    Cn:=ExpandQop(C);

    Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

     

    Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

     

    "Maple do not compute the tensor product of operators,"

     

    "C=A*B gives:"

     

    Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

     

    "ExpandQop(C) gives the expected result:"

     

    Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

    (3)

    kets3:=kets1*kets2;
    bras3:=Dagger(kets3);
    print("Matrix element computed with C appears curious");
    'bras3.C. kets3'="...";
    bras3.C.kets3;
    print("Matrix element computed with Cn as expected");
    'bras3.Cn.kets3'=bras3.Cn.kets3;

    Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3))

     

    Physics:-`*`(Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

     

    "Matrix element computed with C"

     

    Physics:-`.`(bras3, C, kets3) = "..."

     

    Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(b3))*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))

     

    "Matrix element computed with Cn as expected"

     

    Physics:-`.`(bras3, Cn, kets3) = Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))^2*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))^2*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))^2*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))^2*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))^2*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))^2

    (4)

     


     

    Download ExpandQop.mw

     

     

    Lesson_on_functions.mws

    As the title says, a lesson on functions:  eg the -> operator, f(2), eval, evalf etc 

    Using the learning sequence as an alternative to learn problems related to "balance of a body" is shown in this video; thanks to the kindness that Maple offers us in its fundamental programming syntax.

    Balance_of_a_body_with_learning_sequence.mw

    Lenin Araujo

    Ambassador Of Maple

     

     

    One case where an "expansion beyond all orders" may be needed is investigating the asymptotic behavior of the difference of two functions with coinciding dominant series.

    We are interested in the asymptotic behavior of F(z) for large positive z:

    h1 := proc (z) options operator, arrow; hypergeom([1/2], [5/4, 3/2, 7/4], z) end proc; h2 := proc (z) options operator, arrow; (3/4)*sqrt(2*Pi)*hypergeom([1/4], [3/4, 5/4, 3/2], z)/(GAMMA(3/4)*z^(1/4)) end proc; F := proc (z) options operator, arrow; h1(z)-h2(z)+(3/8)*sqrt(Pi)/sqrt(z) end proc

    series does not succeed:

    series(F(z), z = infinity, 20)

    O((1/z)^(23/3))*exp(3/(1/z)^(1/3))

    (1)

    The reason is that the dominant terms containing the factor exp(3*z^(1/3)) in the two hypergeometric functions cancel exactly, and we have to look for the subdominant terms.

    The order of the leading terms can be found from DETools:-formal_sol:

    deq1 := FunctionAdvisor(DE, h1(z), f(z))[2, 1]

    diff(diff(diff(diff(f(z), z), z), z), z) = -(15/2)*(diff(diff(diff(f(z), z), z), z))/z-(195/16)*(diff(diff(f(z), z), z))/z^2+(1/32)*(32*z-105)*(diff(f(z), z))/z^3+(1/2)*f(z)/z^3

    (2)

    DETools:-formal_sol(deq1, f(z), z = infinity, order = 0)

    [(1/z)^(1/2), -exp(-3/(-1/z)^(1/3))/z, -exp(3*(-1)^(1/3)/(-1/z)^(1/3))/z, -exp(-3*(-1)^(2/3)/(-1/z)^(1/3))/z]

    (3)

    As expected, one of the solutions (the third one for positive z) contains the exp(3*z^(1/3)) factor, the leading term being of the order exp(3*z^(1/3))/z.

    Another, subdominant, solution is algebraic and, in fact, is a series containing only one term, as 1/z^(1/2) is an exact solution. It will turn out that the algebraic part in F(z) also cancels out.

    Thus we have to look for the subsubdominant terms, which contain decaying exponentials. We will accomplish this by applying the steepest descent method to the integral representations of h1(z) and h2(z).

    ms := convert([h1(z), h2(z)], MeijerG)

    [(3/32)*Pi*2^(1/2)*MeijerG([[1/2], []], [[0], [-1/4, -1/2, -3/4]], -z), (3/32)*2^(1/2)*Pi*MeijerG([[3/4], []], [[0], [1/4, -1/4, -1/2]], -z)/z^(1/4)]

    (4)

    m2g := proc (m, y) local a, b, c, d; a, b := op(op(1, m)); c, d := op(op(2, m)); -((1/2)*I)*mul(`~`[GAMMA](`~`[`-`](1+y, a)))*mul(`~`[GAMMA](`~`[`-`](c, y)))*op(3, m)^y/(Pi*mul(`~`[GAMMA](`~`[`-`](b, y)))*mul(`~`[GAMMA](`~`[`-`](1+y, d)))) end proc

    gs := applyrule(conditional(e::anything, _op(0, e) = MeijerG) = 'm2g(e, y)', ms)

    [-((3/64)*I)*2^(1/2)*GAMMA(1/2+y)*GAMMA(-y)*(-z)^y/(GAMMA(5/4+y)*GAMMA(3/2+y)*GAMMA(7/4+y)), -((3/64)*I)*2^(1/2)*GAMMA(1/4+y)*GAMMA(-y)*(-z)^y/(z^(1/4)*GAMMA(3/4+y)*GAMMA(5/4+y)*GAMMA(3/2+y))]

    (5)

    gs[2] := combine(eval(gs[2], [1/z^(1/4) = exp(I*Pi*(1/4))/(-z)^(1/4), y = y+1/4]), power)

    -((3/64)*I)*GAMMA(1/2+y)*((1/2)*2^(1/2)+((1/2)*I)*2^(1/2))*(-z)^y*GAMMA(-1/4-y)*2^(1/2)/(GAMMA(1+y)*GAMMA(3/2+y)*GAMMA(7/4+y))

    (6)

    h1(z) and h2(z)are the integrals of gs[1] and of gs[2] over the same path, which is a loop encircling the poles ofGAMMA(-y) and of GAMMA(-1/4-y). Now a standard technique is to extend the integration contour far to the left, while still keeping both endpoints at "+infinity". Then the arguments of the gamma functions can be made large everywhere on the integration path, and the gamma functions can be replaced by their asymptotic approximations.

    When moving the contour, we have to take into account the pole of the integrand at y = -1/2. The other poles of GAMMA(1/2+y) will be cancelled by the zeros of 1/GAMMA(3/2+y), which is why the algebraic part of the expansion will contain the single term of the order 1/z^(1/2).

    This is the negative of the third term in F(z):

    `assuming`([simplify((2*Pi*I)*residue(gs[1]-gs[2], y = -1/2))], [z > 0])

    -(3/8)*Pi^(1/2)/z^(1/2)

    (7)

    Expanding the gamma functions produces terms containing exp(-I*Pi*y) and exp(I*Pi*y)

    `assuming`([simplify(convert(MultiSeries:-series((gs[1]-gs[2])/z^y, y = -infinity, 1), polynom))], [z > 0]); collect(convert(%, exp), exp)

    (-3/64-(3/64)*I)*(-1/y)^(1/2)*exp(-3*(ln(-y)-1)*y)*exp(-I*Pi*y)/(y^3*Pi^(1/2))+(3/64-(3/64)*I)*(-1/y)^(1/2)*exp(-3*(ln(-y)-1)*y)*exp(I*Pi*y)/(y^3*Pi^(1/2))

    (8)

    As we shall see, those terms have saddle points y0(z) = exp(`&+-`((1/3)*(2*Pi*I)))*z^(1/3) located in the left half-plane and contribute exponentially small factors exp(3*y0(z)). The terms for which the saddle point would be located at y = z^(1/3) have cancelled out, thus cancelling the exponentially large contributions. Another possible way to achieve the same result was to write h1(z)-h2(z) as a single Meijer G-function -(3/32)*MeijerG([[1/2], []], [[-1/4, 0], [-3/4, -1/2]], z).

    We write the first term above in the form g(y)*exp(f(y)):

    f := proc (z, y) options operator, arrow; -3*y*(ln(-y)-1)-I*Pi*y+y*ln(z) end proc

    g := proc (y) options operator, arrow; (-3/64-(3/64)*I)*sqrt(-1/y)/(sqrt(Pi)*y^3) end proc

    diff(f(z, y), y)

    -3*ln(-y)-I*Pi+ln(z)

    (9)

    For this to become zero, we need argument(-y) = -(1/3)*Pi, and thus y = exp((1/3)*(2*I)*Pi)*z^(1/3). We can visualize the paths where the imaginary part of f(z, y) stays constant. The path of the steepest descent is the one that goes through the saddle point in the direction exp(I*Pi*(1/3)); the blue color indicates smaller values of the real part of f(z, y):

    y0 := proc (z) options operator, arrow; exp(((2/3)*I)*Pi)*z^(1/3) end proc

    (proc () local z; z := 2; plots:-display(plots:-contourplot(Re(f(z, u+I*v)), u = -5 .. 5, v = -5 .. 5, contours = ([seq])(Re(f(z, y0(z)))+i, i = -30 .. 6, 6), filledregions, coloring = [blue, red], grid = [100, 100]), plots:-implicitplot(Im(f(z, u+I*v)-f(z, y0(z))), u = -5 .. 5, v = -5 .. 5, gridrefine = 5, color = green), plot([cos((1/3)*Pi)*xi+Re(y0(z)), sin((1/3)*Pi)*xi+Im(y0(z)), xi = -3 .. 3], linestyle = dot, color = white), axes = boxed) end proc)()

     

    The real part of f(z, y) has a maximum along this path at y0(z).

    `assuming`([(`@`(`@`(simplify, evalc), series))(f(z, y0(z)+exp(I*Pi*(1/3))*xi), xi = 0, 3)], [z > 0]); quad := convert(%, polynom)

    series((3/2)*(-1+I*3^(1/2))*z^(1/3)-((3/2)/z^(1/3))*xi^2+O(xi^3),xi,3)

    (10)

    Now we can compute the lead asymptotic term contributed by the saddle point y0(z):

    lt1 := `assuming`([(`@`(simplify, evalc))(g(y0(z))*exp(I*Pi*(1/3))*(int(exp(quad), xi = -infinity .. infinity)))], [z > 0])

    -(1/64)*exp(-(3/2)*z^(1/3))*3^(1/2)*((-1+I)*cos((3/2)*z^(1/3)*3^(1/2))+(-1-I)*sin((3/2)*z^(1/3)*3^(1/2)))*2^(1/2)/z

    (11)

    We repeat the same procedure for the second term of the integrand.

    f := proc (z, y) options operator, arrow; -3*y*(ln(-y)-1)+I*Pi*y+y*ln(z) end proc

    g := proc (y) options operator, arrow; (3/64-(3/64)*I)*sqrt(-1/y)/(sqrt(Pi)*y^3) end proc

    diff(f(z, y), y)

    -3*ln(-y)+I*Pi+ln(z)

    (12)

    y0 := proc (z) options operator, arrow; exp(-((2/3)*I)*Pi)*z^(1/3) end proc

    The direction should be chosen as exp((1/3)*(2*I)*Pi) to be consistent with the direction of the integration contour, which goes from the lower to the upper half-plane.

    lterm := proc (gy, fy, eq, dir) options operator, arrow; (eval(gy*exp(fy), eq))*dir*sqrt(-2*Pi/((eval(diff(fy, `$`(y, 2)), eq))*dir^2)) end proc

    lt2 := `assuming`([(`@`(simplify, evalc))(lterm(g(y), f(z, y), y = y0(z), exp((1/3)*(2*I)*Pi)))], [z > 0])

    (1/64)*exp(-(3/2)*z^(1/3))*((1+I)*cos((3/2)*z^(1/3)*3^(1/2))+(1-I)*sin((3/2)*z^(1/3)*3^(1/2)))*3^(1/2)*2^(1/2)/z

    (13)

    Combining the two results yields the leading term of F(z). The next terms can be obtained by expanding gs[1] and gs[2] to higher orders.

    Fasympt := unapply(simplify(lt1+lt2), z)

    proc (z) options operator, arrow; (1/32)*exp(-(3/2)*z^(1/3))*3^(1/2)*2^(1/2)*(cos((3/2)*z^(1/3)*3^(1/2))+sin((3/2)*z^(1/3)*3^(1/2)))/z end proc

    (14)

    (proc () Digits := 50; plot(`~`[`*`](exp((3/2)*z^(1/3)), [F(z), Fasympt(z)]), z = 1000 .. 10000, linestyle = [solid, dot], thickness = [1, 5], axes = frame) end proc)()

     

    Download steep.mw

    Just a simple little worksheet to see if I have enough propane to heat my house for the rest of the winter.


     

    Do I have enough propane for the winter?

    NULL

    I've taken some measurements from my propane tank throughout the winter.  Now we can use Maple to see if we have enough to last the rest of the winter.

    ``

    a := [["nov 27, 2017", 73.5], ["dec 9, 2017", 72], ["dec 16, 2017", 69], ["dec 31, 2017", 62], ["jan 12, 2018", 60], ["jan 19, 2018", 56], ["jan 26, 2018", 54], ["feb 4,2018", 51]]

    [["nov 27, 2017", 73.5], ["dec 9, 2017", 72], ["dec 16, 2017", 69], ["dec 31, 2017", 62], ["jan 12, 2018", 60], ["jan 19, 2018", 56], ["jan 26, 2018", 54], ["feb 4,2018", 51]]

    (1)

    with(Finance)  ``

    pts := [seq([DayCount(a[1, 1], a[i, 1]), a[i, 2]], i = 1 .. nops(a))]

    [[0, 73.5], [12, 72], [19, 69], [34, 62], [46, 60], [53, 56], [60, 54], [69, 51]]

    (2)

    with(plots)

    listplot(pts)

     

    Adding a 30% and 20% level to the graph.  We probably shouldn't be too worried about the cold in June so DayCount("Nov 27, 2017", "Jun 1, 2018") = 186 we'll extend these reference lines out to 186.

    plot({pts, [[0, 20], [186, 20]], [[0, 30], [186, 30]]}, view = [default, 0 .. 80])

     

     

    30% is the recommended level your propane company wants you to fill up at.  The technician who installed the tank said 20% is all right.  It's up to you if you want to go to 10% but if you run out of propane the company has to come in and do a leak test on your system which is an added cost you don't want.  So let's predict at what point we need to start worrying about filling up our propane tank.  To do that, of course, all we need is a forecast line.  For that we'll just calculate a best fit.

     

    a1 := [seq(DayCount(a[1, 1], a[i, 1]), i = 1 .. nops(a))]

    [0, 12, 19, 34, 46, 53, 60, 69]

    (3)

    a2 := a[() .. (), 2]

    [73.5, 72, 69, 62, 60, 56, 54, 51]

    (4)

    X := convert(a1, Vector)

    Y := convert(a2, Vector)

    with(Statistics)

    L1 := LinearFit([1, x], X, Y, x)

    HFloat(74.79237702730747)-HFloat(0.34416046490941915)*x

    (5)

    Plotting it all together

    plot({L1, pts, [[0, 20], [186, 20]], [[0, 30], [186, 30]]}, x = 0 .. 200, y = 0 .. 80, labels = ["Days", ""], tickmarks = [default, [seq(10*i = cat(10*i, "%"), i = 1 .. 8)]])

     

    Projecting the line to 30% we get

    solve(L1 = 30)

    130.1496877

    (6)

    AdvanceDate(a[1, 1], trunc(solve(L1 = 30)))

    Record(monthDay = 6, month = 4, year = 2018, format = "%B %e, %Y", ModulePrint = proc (m) Finance:-FormatDate(m) end proc)

    (7)

    April is still a bit chilly so maybe if we wait until 20%, of course it's getting warmer all this time so our usage should go down.  

    AdvanceDate(a[1, 1], trunc(solve(L1 = 20)))

    Record(monthDay = 5, month = (), year = 2018, format = "%B %e, %Y", ModulePrint = proc (m) Finance:-FormatDate(m) end proc)

    (8)

    It isn't warm enough to turn off the furnace yet but it looks like we'll have enough to get us into the warm months

    AdvanceDate(a[1, 1], trunc(solve(L1 = 10)))

    Record(monthDay = 3, month = 6, year = 2018, format = "%B %e, %Y", ModulePrint = proc (m) Finance:-FormatDate(m) end proc)

    (9)

    We'll hit 10% well into late spring and almost right into summer of course it's a rough estimate however it looks like we won't have to fill up during the high price winter season.  I can tell my wife to relax, we should have enough propane for the winter.

     

     

    NULL


     

    Download Propane_usage-.mw

    On the example of a manipulator with three degrees of freedom.
    A mathematical model is created that takes into account degrees of freedom of the manipulator and the trajectory of the movement from the initial point to the final one (in the figure, the ends of the red curve). In the text of the program, these are the equations fi, i = 1..5.
    Obviously, the straight line could be the simplest trajectory, but we will consider a slightly different variant. The solution of the system of equations is the coordinates of the points of the manipulator (x1, x2, x3) and (x4, x5, x6) in all trajectory. After that, knowing the lengths of the links and the coordinates of the points at each moment of time, any angles of the manipulator are calculated. The same selected trajectory is reproduced from these angles. The possible angles are displayed by black color.
    All the work on creating a mathematical model and calculating the angles can be done without the manipulator itself, is sufficient to have only the instruction with technical characteristics.
    To display some angles, the procedure created by vv is used.
    MAN_2.mw

    We’re kicking off 2018 right, with another Meet Your Developers interview! This edition comes from Erik Postma, Manager of the Mathematical Software Group.

    To catch up on previous interviews, search the “meet-your-developers” tag.

    Without further ado…

     

    1. What do you do at Maplesoft?
      I’m the manager of the mathematical software group, a team of 7 mathematicians and computer scientists working on the mathematical algorithms in Maple (including myself). So my work comes in two flavours: I do the typical managerial things, involving meetings to plan new features and solve my team’s day to day problems, and in the remaining time I do my own development work.
       
    2. What did you study in school?
      I studied at Eindhoven University of Technology in the Netherlands. The first year, I took a combined program of mathematics and computer science; then for the rest of my undergrad, I studied mathematics. The program was called Applied Mathematics, but with the specialization I took it really wasn’t all that applied at all. Afterwards I continued in the PhD program at the same university, where my thesis was on a subject in abstract algebra (Lie algebras over finite fields).
       
    3. What area(s) of Maple are you currently focusing on in your development?
      I’ve spent quite a bit of time over the past two years making the facilities for working with units of measurement in Maple easier to use. There is a very powerful package for doing this that has been part of Maple for many years, but we keep hearing from our users it’s difficult to use. So I’ve worked on keeping the power of the package but making it easier to use.
       
    4. What’s the coolest feature of Maple that you’ve had a hand in developing?
      This was actually working on a problem in a part of the code that existed long before I started with Maplesoft. We have a very clever algorithm for drawing random numbers according to a custom, user-specified probability distribution. I wrote about it on MaplePrimes in a series of four blog posts, here. I’ve talked at various workshops and the like about this algorithm and how it is implemented in Maple.
       
    5. What do you like most about working at Maplesoft? How long have you worked here?
      I love working at the crossroads of mathematics and computer science; there aren’t many places in the world where you can do that as much as at Maplesoft. But the best thing is the people I work with: us mathematicians are all crazy in slightly different ways, and that makes for a very interesting working environment.
       
    6. Favourite hobby?
      Ultimate frisbee. I captain a mixed (i.e., coed) team called The Clockwork. (We play in orange jerseys – it references the book/movie A Clockwork Orange.) We play in a couple of local leagues, and some of the other members also work here. We don’t win much – but we work hard and have fun!
       
    7. What do you like on your pizza?
      Mushrooms. Mushrooms on everything!
       
    8. What’s your favourite movie?
      Probably Black Book, a dark movie about the Dutch resistance in the second world war from 2006, directed by Paul Verhoeven. I think what I like best about it is that it highlights the moral shades of grey in even so morally elevated a group as the resistance.
       
    9. What skill would you love to learn? Why?
      I’d love to learn to speak Russian! I’m trying, but I have a very hard time with it. It would allow me to communicate with my in-laws more easily; they speak Russian.
       
    10. Who’s your favourite mathematician?
      Oh, so many to choose from! I’m torn between:
    • Ada Lovelace (1815-1852), known as the first programmer.
    • Felix Klein (1849-1925), driving force behind a lot of research into geometries and their underlying symmetry groups.
    • Wilhelm Killing (1847-1923), a secondary school teacher who made big contributions to the theory of Lie algebras.

    Or wait, can I choose my wife?

    Please take a look at the attached document, a partial design for a power supply I'm working on.  I find I am spending a lot of time reformatting results with units to look as nice as what you see here.  For every result, I need to do Units Formatting, change to a sensible unit like uH instead of 10^-6 H, and then do Numeric Formatting to change the number to show just three significant digits.  That requires from 0 to 2 decimals, in fixed point.

    This is the way engineering documents should look.  You want to see a fixed point number from 1.00 to 999, with a certain number of significant digits (not decimal points), and have the unit scaled accordingly.  You want to see 12.3 uA, not 1.23402 x 10^-5 A.

    I would like to see Maple add "N significant digits" to its Numeric Formatting options and auto-scale results with units to the appropriate multiplier.  If I could set that as my default result formatting it would save a huge amount of work.  Often as a design progresses the multiplier will change, also.  A result may initially come out in mA but later change to uA.  Not only do I have to do them all manually now, but I have to go back and change them.  Automating all that would be a great help.

    (You may also notice that my vector results with units are not scaled like I describe here.  If anyone can tell me how to do that I would appreciate it.  Otherwise, it looks like a bug to me.)

    Example_Document.zip

    First 46 47 48 49 50 51 52 Last Page 48 of 304