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
  • I submit a bug through MaplePrimes because I can't do it as usually (Hope some people understand me.). Let us consider

    M := Matrix(5, 5,  (i, j) -> (10*i+j)*sin((1/180)*Pi*(10*i+j))):
     #One sees a long and wrong output instead of the warning "Matrix M is singular"


    Digits := 500; evalf(Determinant(M), 495);
                               1.3 10 ^(-488)   

    I submit a bug through MaplePrimes because I can't do it as usually (Hope some people understand me.). Let us consider

    restart; pdsolve([diff(u(t, x), t, t) = diff(u(t, x), x, x), u(t, 0) = 0, u(t, Pi) = 0]);
    pdsolve([diff(u(t, x), t, t) = diff(u(t, x), x, x), u(t, 0) = 0, u(t, Pi) = 0], generalsolution);
    u(t, x) = Sum(sin(n*x)*(_C5(n)*cos(n*t)+_C1(n)*sin(n*t)), n = 1 .. infinity)
    u(t, x) = Sum(sin(n*x)*(_C5(n)*cos(n*t)+_C1(n)*sin(n*t)), n = 1 .. infinity)

    The question arises: what do these outputs mean? I don't see any explanation in ?pdsolve and ?examples,pdsolve_boundaryconditions. What are _C1(n) and _C5(n)? Under which conditions does the above series converge?


    pdetest(%, [diff(u(t, x), t, t) = diff(u(t, x), x, x), u(t, 0) = 0, u(t, Pi) = 0]);
                               [0, 0, 0]

    I think the above is simply a fake: it is possible to differentiate  a series only under certain conditions.

    Please, don't convert my post to a question. This is not correct and fair. Hope some people understand me.

    Inspired by this question. There is a simple iterated procedure that can generate the Puiseux expansion.

    Note the use of convert(..., rational, exact) to preserve the digits of the floating-point numbers.

    (proc () global COcrit, P; Digits := 30; COcrit := 1/convert(2*0.115e-12, rational, exact); P := (`@`(`@`(rcurry(unapply, [NO2, O3]), numer), rcurry(convert, rational, exact)))(1.8027350279544425625*O3^8/10^49+(2.982942852010948125*CO/10^49+2.27611508925183215625*NO2/10^47+3.754849807690565625/10^37)*O3^7+(1.2339511797887015625*CO^2/10^49+2.4836920553034140625*CO/10^37+(-1)*5.9759150186390484375*NO2/10^35+6.3862221928648528125*NO2^2/10^46+1.88311680046014890625*CO*NO2/10^47+(-1)*9.69704144499346875/10^24)*O3^6+((-1)*1.71193039685859375*CO^2/10^37+(-1)*5.7098636544065625*CO/10^24+(-1)*1.277325786277575*NO2/10^21+(-1)*1.0801570017944671875*NO2^2/10^32+(-1)*2.9081778565815421875*CO*NO2/10^34+(-1)*3.66453227489203125/10^11)*O3^5+(1.9152220505625*CO^2/10^25+8.1035862984375*CO/10^13+1.40846609345625*NO2/10^10+(-1)*5.19285353257125*NO2^2/10^21+(-1)*1.55036925507375*CO*NO2/10^21+(-1)*1.844759695120875*CO^2*NO2/10^34+(-1)*1.876842472427325*CO*NO2^2/10^32-7.201634275625)*O3^4+(1.793091274970625*NO2^2/10^7+8618.14231275*NO2+(-1)*2.298266460675*CO^2*NO2/10^22+(-1)*9.5902239009375*CO*NO2/10^10+(-1)*1.685705248305*CO*NO2^2/10^20+9.2666503797075*NO2^3/10^19)*O3^3+((-1)*2.5638555726*10^6*NO2^2+6.894799382025*NO2^2*CO^2/10^20+2.7563954788125*CO*NO2^2/10^7+3.5073544682475*NO2^3*CO/10^18+(-1)*0.604340578881e-4*NO2^3+(-1)*3.5683519372605*NO2^4/10^16)*O3^2+((-1)*8.75499651*10^6*NO2^3+0.482686765875e-5*NO2^3*CO+(-1)*0.98216263665e-4*NO2^4)*O3+5.4066609*10^7*NO2^4) end proc)()

    We're interested in the asymptotics of RootOf(P(NO2, _Z)) for small NO2.

    plot3d(RootOf(P(NO2, _Z)), CO = .5*COcrit .. 2*COcrit, NO2 = 10^(-3) .. 10^(-6))

    Start with P(NO2, Z) and look for the expansion for Z when NO2 is small.

    Expand P(NO2, Z) into monomials and convert each monomial a*Typesetting:-mi("NO2", italic = "true", mathvariant = "italic")^Typesetting:-mi("p", italic = "true", mathvariant = "italic")*Typesetting:-mi("Z", italic = "true", mathvariant = "italic")^Typesetting:-mi("q", italic = "true", mathvariant = "italic") into the point [p, q].

    ([op])(collect(P(NO2, Z), [NO2, Z], distributed, normal)); map(proc (e) options operator, arrow; `~`[degree](e, [NO2, Z]) end proc, %)

    [[4, 2], [4, 1], [4, 0], [3, 3], [3, 2], [3, 1], [2, 6], [2, 5], [2, 4], [2, 3], [2, 2], [1, 7], [1, 6], [1, 5], [1, 4], [1, 3], [0, 8], [0, 7], [0, 6], [0, 5], [0, 4]]


    The Newton polygon is the convex hull of the points.

    newtonPoly := proc (P, xy) local terms, pts; terms := ([op])(collect(P, xy, distributed, normal)); pts := map(proc (e) options operator, arrow; `~`[degree](e, xy) end proc, terms); plots:-display(plottools:-polygon(simplex:-convexhull(pts)), plottools:-point(pts), symbol = solidcircle, symbolsize = 15, style = line, thickness = 3, color = [magenta, black], axis = [gridlines = spacing(1)]) end proc

    newtonPoly(P(NO2, Z), [NO2, Z])


    The side closest to the origin will correspond to the asymptotics for small NO2.

    Sum the corresponding monomials and solve for Z.

    sideAsympt := proc (P, xy, pts) local lead; lead := add(coeff(coeff(P, xy[1], p[1]), xy[2], p[2])*xy[1]^p[1]*xy[2]^p[2], `in`(p, pts)); ([solve])(lead, xy[2]) end proc

    sideAsympt(P(NO2, Z), [NO2, Z], [[0, 4], [1, 3], [2, 2], [3, 1], [4, 0]])

    [600*NO2, 600*NO2, -57000000000000000*NO2/(4071*CO-17795000000000000), -228000000000000000*NO2/(4071*CO+35020000000000000)]


    We have obtained the first term in the expansion. If CO>COcrit, the smallest positive root is the one asymptotic to 600*NO2. That will be the value of RootOf(P(NO2, _Z)).

    The expansion can be continued in the same manner.

    newtonPoly(P(NO2, NO2*(600+Z)), [NO2, Z])


    In this case we don't have to choose between sides with different slopes. (Compare to "P(NO2,600*NO2+Z).)"

    sideAsympt(P(NO2, NO2*(600+Z)), [NO2, Z], [[4, 2], [5, 0]]); `assuming`([simplify(`assuming`([simplify(%)], [NO2 > 0]))], [CO > COcrit])

    [(531/430000)*354^(1/2)*NO2^(1/2)*(200000000000000+23*CO)^(1/2)/(23*CO-100000000000000)^(1/2), -(531/430000)*354^(1/2)*NO2^(1/2)*(200000000000000+23*CO)^(1/2)/(23*CO-100000000000000)^(1/2)]


    We get `~`[Z]*NO2^(1/2), so the next order in the expansion is NO2^(3/2).

    Again, if we want the solution that corresponds to RootOf(P(NO2, _Z)), we should take the negative term, as the principal root will be the smaller one.

    Find one more term. To avoid fractional powers, take NO2r = NO2^(1/2).

    approx := 600*NO2r^2+NO2r^3*(-(531/430000)*sqrt(354)*sqrt(200000000000000+23*CO)/sqrt(23*CO-100000000000000)+Z)

    newtonPoly(P(NO2r^2, approx), [NO2r, Z])


    sideAsympt(P(NO2r^2, approx), [NO2r, Z], [[10, 1], [11, 0]])



    We get `~`[Z]*NO2r, so we have obtained the coefficient at NO2r^3*NO2r = NO2^2.

    All of this can be done using the algcurves:-puiseux command. The only difficulty is the unwieldy expressions that algcurves:-puiseux will generate for P.

    Let's also find the NO2^2 term by the method of dominant balance. Suppose that we don't know p yet and are looking for the term alpha*NO2^p = alpha*NO2p.

    approx := 600*NO2r^2-(531/430000)*sqrt(354)*sqrt(200000000000000+23*CO)*NO2r^3/sqrt(23*CO-100000000000000)+alpha*NO2p

    terms := (`@`([op], collect))(P(NO2r^2, approx), [NO2r, NO2p], distributed, normal)

    Find all exponents of NO2, coming from NO2r and from NO2p.

    pows := map(proc (t) options operator, arrow; (1/2)*degree(t, NO2r)+p*degree(t, NO2p) end proc, terms)

    plot(pows, p = 1 .. 2, view = 4 .. 6, annotation = pows)


    At p=2, the two dominant terms NO2^(11/2) and NO2^(7/2+p) can balance each other.

    (`@`(normal, coeff))(subs(NO2p = NO2r^4, P(NO2r^2, approx)), NO2r, 11); ([solve])(%, alpha)






    Quantum Commutation Rules Basics


    Pascal Szriftgiser1 and Edgardo S. Cheb-Terrab2 

    (1) Laboratoire PhLAM, UMR CNRS 8523, Université Lille 1, F-59655, France

    (2) Maplesoft



    In Quantum Mechanics, in the coordinates representation, the component of the momentum operator along the x axis is given by the differential operator

     "`p__x`=-i `ℏ`(∂)/(∂x)  "


    The purpose of the exercises below is thus to derive the commutation rules, in the coordinates representation, between an arbitrary function of the coordinates and the related momentum, departing from the differential representation


    p[n] = -i*`ℏ`*`∂`[n]

    These two exercises illustrate how to have full control of the computational process by using different elements of the Maple language, including inert representations of abstract vectorial differential operators, Hermitian operators, algebra rules, etc.


    These exercises also illustrate a new feature of the Physics package, introduced in Maple 2017, that is getting refined (the computation below requires the Maplesoft updates of the Physics package) which is the ability to perform computations algebraically, using the product operator, but with differential operators, and transform the products into the application of the operators only when we want that, as we do with paper and pencil.


    %Commutator(g(x, y, z), p_) = I*`ℏ`*Nabla(F(X))


    restart; with(Physics); with(Physics[Vectors]); interface(imaginaryunit = i)


    Start setting the problem:


     all ofx, y, z, p__x, p__y, p__z are Hermitian operators


     all of x, y, z commute between each other


     tell the system only that the operators x, y, z are the differentiation variables of the corresponding (differential) operators p__x, p__y, p__z but do not tell what is the form of the operators


    Setup(mathematicalnotation = true, differentialoperators = {[p_, [x, y, z]]}, hermitianoperators = {p, x, y, z}, algebrarules = {%Commutator(x, y) = 0, %Commutator(x, z) = 0, %Commutator(y, z) = 0}, quiet)

    [algebrarules = {%Commutator(x, y) = 0, %Commutator(x, z) = 0, %Commutator(y, z) = 0}, differentialoperators = {[p_, [x, y, z]]}, hermitianoperators = {p, x, y, z}, mathematicalnotation = true]


    Assuming F(X) is a smooth function, the idea is to apply the commutator %Commutator(F(X), p_) to an arbitrary ket of the Hilbert space Ket(psi, x, y, z), perform the operation explicitly after setting a differential operator representation for `#mover(mi("p",mathcolor = "olive"),mo("→"))`, and from there get the commutation rule between F(X) and `#mover(mi("p",mathcolor = "olive"),mo("→"))`.


    Start introducing the commutator, to proceed with full control of the operations we use the inert form %Commutator

    alias(X = (x, y, z))


    ` F`(X)*`will now be displayed as`*F


    %Commutator(F(X), p_)*Ket(psi, X)

    Physics:-`*`(%Commutator(F(X), p_), Physics:-Ket(psi, x, y, z))


    For illustration purposes only (not necessary), expand this commutator

    Physics[`*`](%Commutator(F(X), p_), Physics[Ket](psi, x, y, z)) = expand(Physics[`*`](%Commutator(F(X), p_), Physics[Ket](psi, x, y, z)))

    Physics:-`*`(%Commutator(F(X), p_), Physics:-Ket(psi, x, y, z)) = Physics:-`*`(F(X), p_, Physics:-Ket(psi, x, y, z))-Physics:-`*`(p_, F(X), Physics:-Ket(psi, x, y, z))


    Note that  `#mover(mi("p",mathcolor = "olive"),mo("→"))`, F(X) and the ket Ket(psi, x, y, z) are operands in the products above and that they do not commute: we indicated that the coordinates x, y, z are the differentiation variables of `#mover(mi("p",mathcolor = "olive"),mo("→"))`. This emulates what we do when computing with these operators with paper and pencil, where we represent the application of a differential operator as a product operation.


    This representation can be transformed into the (traditional in computer algebra) application of the differential operator when desired, as follows:

    Physics[`*`](%Commutator(F(X), p_), Physics[Ket](psi, x, y, z)) = Library:-ApplyProductsOfDifferentialOperators(Physics[`*`](%Commutator(F(X), p_), Physics[Ket](psi, x, y, z)))

    Physics:-`*`(%Commutator(F(X), p_), Physics:-Ket(psi, x, y, z)) = Physics:-`*`(F(X), p_(Physics:-Ket(psi, x, y, z)))-p_(Physics:-`*`(F(X), Physics:-Ket(psi, x, y, z)))


    Note that, in `#mover(mi("p",mathcolor = "olive"),mo("→"))`(F(X)*Ket(psi, x, y, z)), the application of `#mover(mi("p",mathcolor = "olive"),mo("→"))` is not expanded: at this point nothing is known about  `#mover(mi("p",mathcolor = "olive"),mo("→"))` , it is not necessarily a linear operator. In the Quantum Mechanics problem at hands, however, it is. So give now the operator  `#mover(mi("p",mathcolor = "olive"),mo("→"))` an explicit representation as a linear vectorial differential operator (we use the inert form %Nabla, %Nabla, to be able to proceed with full control one step at a time)

    p_ := proc (f) options operator, arrow; -I*`ℏ`*%Nabla(f) end proc

    proc (f) options operator, arrow; -Physics:-`*`(Physics:-`*`(I, `ℏ`), %Nabla(f)) end proc


    The expression (1.5) becomes

    Physics[`*`](%Commutator(F(X), p_), Physics[Ket](psi, x, y, z)) = Physics[`*`](F(X), p_(Physics[Ket](psi, x, y, z)))-p_(Physics[`*`](F(X), Physics[Ket](psi, x, y, z)))

    Physics:-`*`(%Commutator(F(X), p_), Physics:-Ket(psi, x, y, z)) = -I*`ℏ`*Physics:-`*`(F(X), %Nabla(Physics:-Ket(psi, x, y, z)))+I*`ℏ`*%Nabla(Physics:-`*`(F(X), Physics:-Ket(psi, x, y, z)))


    Activate now the inert operator VectorCalculus[Nabla] and simplify taking into account the algebra rules for the coordinate operators {%Commutator(x, y) = 0, %Commutator(x, z) = 0, %Commutator(y, z) = 0}

    Simplify(value(Physics[`*`](%Commutator(F(X), p_), Physics[Ket](psi, x, y, z)) = -I*`ℏ`*Physics[`*`](F(X), %Nabla(Physics[Ket](psi, x, y, z)))+I*`ℏ`*%Nabla(Physics[`*`](F(X), Physics[Ket](psi, x, y, z)))))

    Physics:-`*`(Physics:-Commutator(F(X), p_), Physics:-Ket(psi, x, y, z)) = I*`ℏ`*_i*Physics:-`*`(diff(F(X), x), Physics:-Ket(psi, x, y, z))+I*`ℏ`*_j*Physics:-`*`(diff(F(X), y), Physics:-Ket(psi, x, y, z))+I*`ℏ`*_k*Physics:-`*`(diff(F(X), z), Physics:-Ket(psi, x, y, z))


    To make explicit the gradient in disguise on the right-hand side, factor out the arbitrary ket Ket(psi, x, y, z)

    Factor(Physics[`*`](Physics[Commutator](F(X), p_), Physics[Ket](psi, x, y, z)) = I*`ℏ`*_i*Physics[`*`](diff(F(X), x), Physics[Ket](psi, x, y, z))+I*`ℏ`*_j*Physics[`*`](diff(F(X), y), Physics[Ket](psi, x, y, z))+I*`ℏ`*_k*Physics[`*`](diff(F(X), z), Physics[Ket](psi, x, y, z)))

    Physics:-`*`(Physics:-Commutator(F(X), p_), Physics:-Ket(psi, x, y, z)) = I*`ℏ`*Physics:-`*`((diff(F(X), y))*_j+(diff(F(X), z))*_k+(diff(F(X), x))*_i, Physics:-Ket(psi, x, y, z))


    Combine now the expanded gradient into its inert (not-expanded) form

    subs((Gradient = %Gradient)(F(X)), Physics[`*`](Physics[Commutator](F(X), p_), Physics[Ket](psi, x, y, z)) = I*`ℏ`*Physics[`*`]((diff(F(X), y))*_j+(diff(F(X), z))*_k+(diff(F(X), x))*_i, Physics[Ket](psi, x, y, z)))

    Physics:-`*`(Physics:-Commutator(F(X), p_), Physics:-Ket(psi, x, y, z)) = I*`ℏ`*Physics:-`*`(%Gradient(F(X)), Physics:-Ket(psi, x, y, z))


    Since (1.10) is true for allKet(psi, x, y, z), this ket can be removed from both sides of the equation. One can do that either taking coefficients (see Coefficients ) or multiplying by the "formal inverse" of this ket, arriving at the (expected) form of the commutation rule between F(X) and `#mover(mi("p",mathcolor = "olive"),mo("→"))`

    (Physics[`*`](Physics[Commutator](F(X), p_), Ket(psi, x, y, z)) = I*`ℏ`*Physics[`*`](%Gradient(F(X)), Ket(psi, x, y, z)))*Inverse(Ket(psi, x, y, z))

    Physics:-Commutator(F(X), p_) = I*`ℏ`*%Gradient(F(X))


    Tensor notation, "[`X__m`,P[n]][-]=i `ℏ` g[m,n]"


    The computation rule for position and momentum, this time in tensor notation, is performed in the same way, just that, additionally, specify that the space indices to be used are lowercase latin letters, and set the relationship between the differential operators and the coordinates directly using tensor notation.

    You can also specify that the metric is Euclidean, but that is not necessary: the default metric of the Physics package, a Minkowski spacetime, includes a 3D subspace that is Euclidean, and the default signature, (- - - +), is not a problem regarding this computation.


    restart; with(Physics); interface(imaginaryunit = i)

    Setup(mathematicalnotation = true, coordinates = cartesian, spaceindices = lowercaselatin, algebrarules = {%Commutator(x, y) = 0, %Commutator(x, z) = 0, %Commutator(y, z) = 0}, hermitianoperators = {P, X, p}, differentialoperators = {[P[m], [x, y, z]]}, quiet)

    [algebrarules = {%Commutator(x, y) = 0, %Commutator(x, z) = 0, %Commutator(y, z) = 0}, coordinatesystems = {X}, differentialoperators = {[P[m], [x, y, z]]}, hermitianoperators = {P, p, t, x, y, z}, mathematicalnotation = true, spaceindices = lowercaselatin]


    Define now the tensor P[m]

    Define(P[m], quiet)

    {Physics:-Dgamma[mu], P[m], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-gamma_[a, b], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}


    Introduce now the Commutator, this time in active form, to show how to reobtain the non-expanded form at the end by resorting the operands in products

    Commutator(X[m], P[n])*Ket(psi, x, y, z)

    Physics:-`*`(Physics:-Commutator(Physics:-SpaceTimeVector[m](X), P[n]), Physics:-Ket(psi, x, y, z))


    Expand first (not necessary) to see how the operator P[n] is going to be applied

    Physics[`*`](Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]), Ket(psi, x, y, z)) = expand(Physics[`*`](Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]), Ket(psi, x, y, z)))

    Physics:-`*`(Physics:-Commutator(Physics:-SpaceTimeVector[m](X), P[n]), Physics:-Ket(psi, x, y, z)) = Physics:-`*`(Physics:-SpaceTimeVector[m](X), P[n], Physics:-Ket(psi, x, y, z))-Physics:-`*`(P[n], Physics:-SpaceTimeVector[m](X), Physics:-Ket(psi, x, y, z))


    Now expand and directly apply in one ago the differential operator P[n]

    Physics[`*`](Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]), Ket(psi, x, y, z)) = Library:-ApplyProductsOfDifferentialOperators(Physics[`*`](Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]), Ket(psi, x, y, z)))

    Physics:-`*`(Physics:-Commutator(Physics:-SpaceTimeVector[m](X), P[n]), Physics:-Ket(psi, x, y, z)) = Physics:-`*`(Physics:-SpaceTimeVector[m](X), P[n](Physics:-Ket(psi, x, y, z)))-P[n](Physics:-`*`(Physics:-SpaceTimeVector[m](X), Physics:-Ket(psi, x, y, z)))


    Introducing the explicit differential operator representation for P[n], here again using the inert %d_[n] to keep control of the computations step by step

    P[n] := proc (f) options operator, arrow; -I*`ℏ`*%d_[n](f) end proc

    proc (f) options operator, arrow; -Physics:-`*`(Physics:-`*`(I, `ℏ`), %d_[n](f)) end proc


    The expanded and applied commutator (2.5) becomes

    Physics[`*`](Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]), Ket(psi, x, y, z)) = Physics[`*`](Physics[SpaceTimeVector][m](X), P[n](Ket(psi, x, y, z)))-P[n](Physics[`*`](Physics[SpaceTimeVector][m](X), Ket(psi, x, y, z)))

    Physics:-`*`(Physics:-Commutator(Physics:-SpaceTimeVector[m](X), P[n]), Physics:-Ket(psi, x, y, z)) = -I*`ℏ`*Physics:-`*`(Physics:-SpaceTimeVector[m](X), %d_[n](Physics:-Ket(psi, x, y, z)))+I*`ℏ`*%d_[n](Physics:-`*`(Physics:-SpaceTimeVector[m](X), Physics:-Ket(psi, x, y, z)))


    Activate now the inert operators %d_[n] and simplify taking into account Einstein's rule for repeated indices

    Simplify(value(Physics[`*`](Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]), Ket(psi, x, y, z)) = -I*`ℏ`*Physics[`*`](Physics[SpaceTimeVector][m](X), %d_[n](Ket(psi, x, y, z)))+I*`ℏ`*%d_[n](Physics[`*`](Physics[SpaceTimeVector][m](X), Ket(psi, x, y, z)))))

    Physics:-`*`(Physics:-Commutator(Physics:-SpaceTimeVector[m](X), P[n]), Physics:-Ket(psi, x, y, z)) = I*`ℏ`*Physics:-g_[m, n]*Physics:-Ket(psi, x, y, z)


    Since the ket Ket(psi, x, y, z) is arbitrary, we can take coefficients (or multiply by the formal Inverse  of this ket as done in the previous section). For illustration purposes, we use   Coefficients  and note hwo it automatically expands the commutator

    Coefficients(Physics[`*`](Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]), Ket(psi, x, y, z)) = I*`ℏ`*Physics[g_][m, n]*Ket(psi, x, y, z), Ket(psi, x, y, z))

    Physics:-`*`(Physics:-SpaceTimeVector[m](X), P[n])-Physics:-`*`(P[n], Physics:-SpaceTimeVector[m](X)) = I*`ℏ`*Physics:-g_[m, n]


    One can undo this (frequently undesired) expansion of the commutator by sorting the products on the left-hand side using the commutator between X[m] and P[n]

    Library:-SortProducts(Physics[`*`](Physics[SpaceTimeVector][m](X), P[n])-Physics[`*`](P[n], Physics[SpaceTimeVector][m](X)) = I*`ℏ`*Physics[g_][m, n], [P[n], X[m]], usecommutator)

    Physics:-Commutator(Physics:-SpaceTimeVector[m](X), P[n]) = I*`ℏ`*Physics:-g_[m, n]


    And that is the result we wanted to compute.


    Additionally, to see this rule in matrix form,

    TensorArray(-(Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]) = I*`ℏ`*Physics[g_][m, n]))

    Matrix(%id = 18446744078261558678)


    In the above, we use equation (2.10) multiplied by -1 to avoid a minus sign in all the elements of (2.11), due to having worked with the default signature (- - - +); this minus sign is not necessary if in the Setup at the beginning one also sets  signature = `+ + + -`


    For display purposes, to see this matrix expressed in terms of the geometrical components of the momentum `#mover(mi("p",mathcolor = "olive"),mo("→"))` , redefine the tensor P[n] explicitly indicating its Cartesian components

    Define(P[m] = [p__x, p__y, p__z], quiet)

    {Physics:-Dgamma[mu], P[m], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-gamma_[a, b], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}


    TensorArray(-(Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]) = I*`ℏ`*Physics[g_][m, n]))

    Matrix(%id = 18446744078575996430)


    Finally, in a typical situation, these commutation rules are to be taken into account in further computations, and for that purpose they can be added to the setup via


    [algebrarules = {%Commutator(x, p__x) = I*`ℏ`, %Commutator(x, p__y) = 0, %Commutator(x, p__z) = 0, %Commutator(x, y) = 0, %Commutator(x, z) = 0, %Commutator(y, p__x) = 0, %Commutator(y, p__y) = I*`ℏ`, %Commutator(y, p__z) = 0, %Commutator(y, z) = 0, %Commutator(z, p__x) = 0, %Commutator(z, p__y) = 0, %Commutator(z, p__z) = I*`ℏ`}]


    For example, from herein computations are performed taking into account that

    (%Commutator = Commutator)(x, p__x)

    %Commutator(x, p__x) = I*`ℏ`







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

    So my latest project was a focus on the utilisation of the Readability command in StringTools, and was used to determine how "readable", or easy to understand, some of the past presidents were. 

    This program takes the Farewell Addresses of some of the past presidents (seen here), and using the Readability command included in the StringTools package, calculates the readability of the various Addresses using multiple methods, as will be outlined below. It then takes the different scores, and plots them on a heat map, to show how each president compares to one another.

    As we can see from the heat map, Harry Truman and Ronald Reagan seem to be the most readable overall, throughout all the methods. Whether this is a good thing or not, is really up to your interpretation. Is using bigger words a bad thing? 

    George Washington and Andrew Jackson have the worst readabiltiy scores, although this is probably in part due to the era of the addresses, and the change in language we've had since then.

    Download the application here: 

    The On-Line Encyclopedia of Integer Sequences (OEIS) is an online database of integer sequences. Originally founded by Neil Sloane in 1964, OEIS has evolved into a larger community based project that today contains information records on over 280,000 sequences. Each record contains information on the leading terms of the sequence, keywords, mathematical motivations, literature links, and more - there’s even Maple code to reproduce many of the sequences!

    You can read more about the OEIS project on the main site, or on Wikipedia. There have also been several articles written about the OEIS project; here’s a (somewhat) recent article on Building a Searching Engine for Mathematics that gives a little more history on the project.

    I wrote a simple package for Maple that can be used to search OEIS for a sequence of numbers or a string and retrieve information on various integer sequences. You can interact with the OEIS package using the commands: Get, Search or Print, or using the right-click context menu.

    Here are some examples:



    Search for a sequence of 6 or more numbers (note that you can also just right-click on an integer sequence to run the search):

    Search(0, 1, 1, 2, 4, 9, 20);
        6 results found
                 81, 255636, 35084, 58386, 5908, 14267


    Retrieve information on the integer sequence A000081. By default, this returns the data as a JSON string and stores the results in a table:

    results81 := Get(81):
        ["revision"], ["number"], ["formula"], ["offset"], ["example"], ["keyword"], ["id"], ["xref"], ["data"], ["reference"], ["mathematica"], ["maple"], ["created"], ["comment"], ["references"], ["time"], ["link"], ["author"], ["program"], ["name"]


    Entries in the table can be accessed as follows:

        "_N. J. A. Sloane_"


    I’ve also added some functionality for printing tables containing information on integer sequences. The Print command displays an embedded table containing all or selected output from records. For example, say we search for a sequence of numbers:

    Search(0, 1, 1, 2, 4, 9, 20);
        6 results found
                 81, 255636, 35084, 58386, 5908, 14267


    Let’s print out selected details on the first three of these sequences (including Maple code if available):

    Print([81, 255636, 35084], output = ["author", "name", "data", "maple"]);


    The OEIS API is somewhat limited, so I have had to put some restrictions in place for returning results. Any query that returns more than 100 entries will return an error and a message to provide a more precise specification for the query. Also, in order to reduce the number of search results, the Search command requires at least 6 integers.


    You can install the OEIS package directly from the MapleCloud or by running:



    Comments and feedback are welcome here, or on the project page:

    I write shell scripts that call Maple to automate frequent tasks. Because I prefer writing Maple code to shell code, I've created a Maple package, Bark, that generates a shell script from Maple source code. It provides a compact notation for defining both optional and positional command-line parameters, and a mechanism to print a help page, from the command-line, for the script. The optional parameters can be both traditional single letter Unix options, or the more expressive GNU-style long options.

    As an example, here is the Maple code, using Bark, for a hello-world script.

    hello := module()
        Parser := Bark:-ArgParser(NULL
                                  , 'prologue' = ( "Print `Hello, World!'" )
                                  , 'opts' = ['help' :: 'help' &c "Print this help page"]
        ModuleApply := proc(cmdline :: string := "")
            Bark:-ArgParser:-Parse(Parser, cmdline);
            Bark:-printf("Hello, World!\n");
        end proc;
    end module:

    The following command creates and installs the shell script in the user's bin directory.

    Bark:-CreateScript("hello", hello
                       , 'add_libname' = Bark:-SaveLib(hello, 'mla' = "hello.mla")

    The hello script is executed from the command-line as

    $ hello
    Hello,  World!

    Pass the -h (or --help) option to display the help.

    $ hello -h
    Usage: hello [-h|--help] [--]
    Print `Hello, World!'
    Optional parameters:
    -h, --help Print this help page

    CreateScript creates two files that are installed in the bin directory: the shell script and a Maple archive file that contains the Maple procedures. The shell script passes its argument in a call to the parser (a Maple procedure) saved in the archive file (.mla file). Here's the created shell script for the hello command:

    #!/usr/bin/env sh
    CMD_LINE=$(echo $0; for arg in "$@"; do printf '%s\n' "$arg"; done)
    echo "hello(\"$CMD_LINE\");" | "$MAPLE" -q -w2 -B --historyfile=none -b '/home/joe/bin/hello.mla'

    I've used Bark on Linux and Windows (with Cygwin tools). It should work on any unix-compatible OS with the Bash shell. If you use a different shell that does not work with it, let me know and I should be able to modify the CreateScript command to have options for specific shells.

    Bark is available on the MapleCloud. To install it, open the MapleCloud palette in Maple, select packages in the drop-down menu and go to the New tab (or possibly the Popular tab). You will also need the TextTools package which is also on the MapleCloud. The intro page for Bark has a command that automatically installs TextTools. Alternatively, executing the following commands in Maple 2017 should install both TextTools and Bark.


    The source for a few useful scripts are included in the examples directory of the installed Bark toolbox. Maple help pages are included with Bark, use "Bark" as the topic.

    I am pleased to announce that a new release of Maple T.A., our online testing and assessment system, is now available. Maple T.A. 2017 includes significant enhancements to learning management system integration, as well as security, performance, and other improvements. These same improvements are also available in a new version of the  Maple T.A. MAA Placement Test Suite.  For more information, see What’s New in Maple T.A. 


    As my first project as a Junior Applications Developer, I set out to learn to code in the best way I know how, by doing.  I ended up picking what was probably one of the hardest options I could pick, namely to replicate the sliding puzzle game, 2048. ( ) Of course I didn’t realize how hard it would be at the time, but after spending the first week alone working on the logic, I had already dug my hole.

    2048, the sliding puzzle game, basically starts with a 4x4 grid filled with zeros. As you swipe the grid, values move toward one of the up, down, left or right sides. With every subsequent swipe, a randomly placed value of 2 or 4 is added to the grid. Any neighbouring matching values in the direction of the swipe are added to one another. This was done by swiping 2 tiles of equal value into each other, creating a new tile with double the value. Two 2 tiles became a 4, two 4’s an 8, and so on.

    The goal? To create a 2048 tile.

    Overall the logic was probably the most challenging part of my task, once the framework was set. . The logic consisted of many if statements that made the numbers “slide” properly, ie not combining with another number. This was probably the hardest part. Troubleshooting and allowing for all the possible conditions also proved difficult. However, the user interface was probably the toughest part, figuring out the labels, making everything display correctly, and programming all the buttons to not break. That was fun.

    Anyway, it was a really fun project to work on, and I’m extremely happy for how it turned out, and I hope you enjoy playing it, as much as I did making it!

    You can try it out here:


    I often use in Maple the plot preview, but any time when i launch the plot preview, i must set the zoom factor to 100% to see all content of the document and further i have to maximize theplot-view window.

    It is very uncomfortable to use the plot preview.

    This should be fixed with the next upgrades of Maple.

    See here the effect, when launching the plot previewer:

    Thanks in advance




    With a member of the community I had some discussion about using Maple for limits of sequences.

    The specific task was about F:=4*sqrt(n)*sin(Pi*sqrt(4*n^2+sqrt(n))) and limit F for n --> infinity for integers, which is asserted to be Pi, see

    For moderate size of integers n it can be 'confirmed' by numerical evaluations. Formally it is not 'obvious' at all.

    However Maple answers by

    limit(F, n = infinity) assuming n::posint;


    This is 'explained' by the help about assuming/details:

    "The assuming command does not place assumptions on integration or summation dummy variables in definite integrals and sums, nor in limit or product dummy variables, because all these variables already have their domain restricted by the integration, summation or product range or by the method used to compute a limit. ..."
    The help continues with suggestions how to treat the situation.

    Which means that the limit is taken in the Reals, not in the discrete Naturals.

    A more simple example may be limit(sin(n*Pi), n = infinity) assuming n::posint which returns "-1 .. 1".

    But here we go:

    MultiSeries:-asympt(F, n):
    simplify(%) assuming n::posint: collect(%, n); #lprint(%);
    limit(%, n=infinity);

       O(1/2048*Pi^3/n^(5/2)) + Pi -1/96*Pi^3/n - 1/16*Pi/n^(3/2) + 1/30720*Pi^5/n^2


    As desired.

    As you know, the MapleCloud is a good way to share all sorts of interactive documents with others, in private groups or so they are accessible to everyone. We recently posted some new content that I thought people might be particularly interested in: a collection of Maple Assistants.


    Up until now, Maple Assistants were only available from within Maple, but now you can take advantage of these powerful tools wherever you are, using your web browser.


    Code Generation  - Translate Maple code to C, Java, Python, R, and more

    Scientific Constants – Explore over 20000 values of physical constants and properties of chemical elements, including units and uncertainty values

    Special Functions – Explore the properties of over 200 special functions, including the Hypergeometric, Bessel, Mathieu, Heun and Legendre families of functions.

    Units Converter – Convert between over 500 units of measurement. (In addition to the standard stuff, you can find out how many fortnights old you are, or how long your commute is in furlongs!)

    I have just posted an article with this title at Maplesoft Application Center here.
    It was motivated by a question posed by  Markiyan Hirnyk  here and a test problem proposed there by Kitonum.

    Now I just want to give the promissed complete solution to Kitonum's test:

    Compute the plane area of the region defined by the inequalities:

    R := [ (x-4)^2+y^2 <= 25, x^2+(y-3)^2 >= 9, (x+sqrt(7))^2+y^2 >= 16 ];

    plots:-inequal(R, x=-7..10, y=-6..6, scaling=constrained);

    The used procedures (for details see the mentioned article):

    ranges:=proc(simpledom::list(relation), X::list(name))
    local rez:=NULL, r,z,k,r1,r2;
    if nops(simpledom)<>2*nops(X) then error "Domain not simple!" fi;
    for k to nops(X) do    r1,r2:=simpledom[2*k-1..2*k][]; z:=X[k];
      if   rhs(r1)=z and lhs(r2)=z then rez:=z=lhs(r1)..rhs(r2),rez; #a<z,z<b
      elif lhs(r1)=z and rhs(r2)=z then rez:=z=lhs(r2)..rhs(r1),rez  #z<b,a<z
      else error "Strange order in a simple domain" fi
    end proc:
    MultiIntPoly:=proc(f, rels::list(relation(ratpoly)), X::list(name))
    local r,rez,sol,irr,wirr, rels1, w;
    rels1:=eval(rels, irr=~wirr);
    sol:=remove(hastype, eval(sol,wirr=~irr), `=`); 
    end proc:
    MeasApp:=proc(rel::{set,list}(relation), Q::list(name='range'(realcons)), N::posint)
    local r, n:=0, X, t, frel:=evalf(rel)[];
    if indets(rel,name) <> indets(Q,name)  then error "Non matching variables" fi;
    r:=[seq(rand(evalf(rhs(t))), t=Q)];
    to N do
      if evalb(eval(`and`(frel), X=~r())) then n:=n+1 fi;
    evalf( n/N*mul((rhs-lhs)(rhs(t)),t=Q) );
    end proc:

    Problem's solution:

    MultiIntPoly(1, R, [x,y]):  # Unfortunately it's slow; patience needed!

    evalf(%) = MeasApp(R, [x=-7..10,y=-6..6], 10000); # A rough numerical check
               61.16217534 = 59.91480000


    This presentation is about magnetic traps for neutral particles, first achieved for cold neutrons and nowadays widely used in cold-atom physics. The level is that of undergraduate electrodynamics and tensor calculus courses. Tackling this topic within a computer algebra worksheet as shown below illustrates well the kind of advanced computations that can be done today with the Physics package. A new feature minimizetensorcomponents and related functionality is used along the presentation, that requires the updated Physics library distributed at the Maplesoft R&D Physics webpage.


    Magnetic traps in cold-atom physics


    Pascal Szriftgiser1 and Edgardo S. Cheb-Terrab2 

    (1) Laboratoire PhLAM, UMR CNRS 8523, Université Lille 1, F-59655, France

    (2) Maplesoft


    We consider a device constructed with a set of electrical wires fed with constant electrical currents. Those wires can have an arbitrary complex shape. The device is operated in a regime such that, in some region of interest, the moving particles experience a magnetic field that varies slowly compared to the Larmor spin precession frequency. In this region, the effective potential is proportional to the modulus of the field: LinearAlgebra[Norm](`#mover(mi("B"),mo("&rarr;"))`(x, y, z)), this potential has a minimum and, close to this minimum, the device behaves as a magnetic trap.




    Figure 1: Schematic representation of a Ioffe-Pritchard magnetic trap. It is made of four infinite rods and two coils.



    Following [1], we show that:



    a) For a time-independent magnetic field  `#mover(mi("B"),mo("&rarr;"))`(x, y, z) in vacuum, up to order two in the relative coordinates X__i = [x, y, z] around some point of interest, the coefficients of orders 1 and 2 in this expansion, `v__i,j` and `c__i,j,k` , respectively the gradient and curvature, contain only 5 and 7 independent components.


    b) All stationary points of LinearAlgebra[Norm](`#mover(mi("B"),mo("&rarr;"))`(x, y, z))^2 (nonzero minima and saddle points) are confined to a curved surface defined by det(`&PartialD;`[j](B[i])) = 0.


    c) The effective potential, proportional to LinearAlgebra[Norm](`#mover(mi("B"),mo("&rarr;"))`(x, y, z)), has no maximum, only a minimum.


    Finally, we draw the stationary condition surface for the case of the widely used Ioffe-Pritchard magnetic trap.






    [1] R. Gerritsma and R. J. C. Spreeuw, Topological constraints on magnetostatic traps,  Phys. Rev. A 74, 043405 (2006)



    The independent components of `v__i,j` and `c__i,j,k` entering B[i] = u[i]+v[i, j]*X[j]+(1/2)*c[i, j, k]*X[j]*X[k]


    The stationary points are within the surface det(`&PartialD;`[j](B[i])) = 0


    U = LinearAlgebra[Norm](`#mover(mi("B",fontweight = "bold"),mo("&rarr;",fontweight = "bold"))`)^2 has only minima, no maxima


    Drawing the Ioffe-Pritchard Magnetic Trap


  or in pdf format with the sections open: MagneticTraps.pdf

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

    Much of this topic is developed using traditional techniques. Maple modernizes and optimizes solutions by displaying the necessary operators and simple commands to solve large problems. Using the conditions of equilibrium for both moment and force we find the forces and moments of reactions for any type of structure. In spanish.

    Lenin Araujo Castillo

    Ambassador of Maple

    First 29 30 31 32 33 34 35 Last Page 31 of 285