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
  • Application to generate vector exercises using very easy to use patterns. These exercises are for sum of vectors in calculations of magnitude, direction, graph and projections. Basically to evaluate our students on the board and thus not lose the results. For university professors and engineering students.

    Generator_of_exercises_with_vectors.mw

    Lenin Araujo Castillo

    Ambassador of Maple

    We added a small new feature to MaplePrimes that will make it easier to tag other members in your questions, posts and comments.

    To use the feature, type "@" in a message. If you wait a moment, a list of names will be populated, starting with members who are participating in the current thread, followed by other members ordered by reputation. You can also begin to type after the "@" to search for a particular member. Similar to the above, members who match your search who are in the current thread will appear first, followed by others.

    If you are tagged in this way, an indication will appear in your notifications pane.

    We hope you find this new feature useful!

     

    Many of you enjoyed our profile on one of our developers, Paulina Chin, so we’re happy to bring you another one!

    Today, we’ll be talking with John May, Senior Developer of Maple. Let’s get started.

    1. What do you do at Maplesoft?
      Until recently I was consulting on-site at the NASA Jet Propulsion Laboratory helping people there more effectively solve their engineering problems using Maplesoft products.  But my main job that I am back to full time now is the development and maintenance of various parts of the Maple library.
       
    2. What did you study in school?
      I studied both Pure and Applied Mathematics at the University of Oregon,  focusing a lot on Abstract Algebra.  In graduate school, I specialized more in computation mathematics like computer algebra and numerical analysis.  My Ph.D. work focused on effective numerical algorithms for problems in polynomial algebra – with implementations in Maple!
       
    3. What area(s) of Maple are you currently focusing on in your development?
      Right now I am focused on addressing complaints I’ve gotten from engineers about the usability of units with other parts of the math library.
       
    4. What’s the coolest feature of Maple that you’ve had a hand in developing?
      A lot of the cool things I’ve built live pretty deep in the internals of Maple.  I’ve done a lot of meta-heuristic tuning to seamlessly integrate high-performance libraries into top-level Maple commands.

      I had a lot of fun developing a lot of the stuff for manipulation and visualization of colors in the ColorTools package.
       
    5. What do you like most about working at Maplesoft? How long have you worked here?
      I started working at Maple in 2007, but I’ve been a Maple user since 1997.  I love being part of the magic that brings powerful algorithmic mathematics to everyone.  The R&D team is also full of eccentric nerds who are great fun to work with.
       
    6. Favourite hobby?
      It varies by the season, but right now it is prime for mountain biking in southern California.  I ride my local trails a couple times a week, and when I get I chance, I love to get away on epic bikepacking adventures (like this one: https://www.bikemag.com/features/two-wheeled-escape-one-hour-from-l-a/  this is me: https://cdn.bikemag.com/uploads/2016/05/16File.jpg ).
       
    7. What do you like on your pizza?
      Anything and everything. Something different every time. My all-time favorite pie my from grad school days is the “Rio Rancho” from the dearly departed That’s Amore Pizza (which was next to the comic book store and across the street from North Carolina State University).  It was an olive oil and mozzarella pizza with chopped bacon that was covered in sliced fresh roma tomatoes and drizzled with ranch dressing when it came out of the oven. 
       
    8. What’s your favourite movie?
      It’s really hard to pick just one.  So, I’ll go with the safe answer and say the greatest movie of all time, and “Weird Al” Yankovic’s only foray into movies, UHF, is my favorite.
      http://www.imdb.com/title/tt0098546/
       
    9. What skill would you love to learn? (That you haven’t already) Why?
      Another hard one.  I feel like I’ve dabbled in lots of things that I would like to get better at.  At the top of the list is probably unicycling.  I’d love to get good enough to play Unicyle Football or do Muni (mountain unicyling).
      https://en.wikipedia.org/wiki/Mountain_unicycling
      http://www.unicyclefootball.com/
       
    10. Who’s your favourite mathematician?
      Batman. https://youtu.be/AcMEckOyoaM

     

    The Perimeter Institute for Theoretical Physics (PI) is a place where bold ideas flourish. It brings together great minds in a shared effort to achieve scientific breakthroughs that will transform our future. PI is an independent research center in foundational theoretical physics.

    One of the key mission objectives of Perimeter is to provide scientific training and educational outreach activities to the general public. Maplesoft is proud to be part of this endeavor, as PI’s Educational Outreach Champion. Maplesoft’s contributions support Perimeter’s Teacher Network training activities, core educational resources, development of online course material, support of events such as the EinsteinPlus workshop for high school teachers, the International Summer School for Young Physicists (ISSYP), and other initiatives.
     

    The annual International Summer School for Young Physicists (ISSYP) is a two-week camp that brings together 40 exceptional physics-minded students from high schools across the globe to PI. The program covers many different topics in physics such as quantum mechanics, special relativity, general relativity, and cosmology. Each year students receive a complimentary copy of Maple, and use the product to practice and strengthen their math skills. The program receives an average of three-hundred applications from students in grades eleven and twelve from around the world. Competition is intense, and students who are chosen for the program are extremely bright and advanced for their age; however there is some variation in their level of math and physics knowledge. . Students are asked to review a “math primer” document to prepare them with the background needed for the program.The ISSYP program now uses Möbius, Maplesoft’s online courseware platform to administer this primer. With Möbius, PI has moved from a pdf document primer to fully online material, which has motivated more students to complete the material and be more engaged in their courses. The interactivity and engagement that technology provides has made the summer program more productive and dynamic.

    EinsteinPlus is a one-week intensive workshop for Canadian and international high school teachers that focuses on modern physics, including quantum physics, special relativity, and cosmology. EinsteinPlus also provides unique opportunities to learn some of the latest developments in physics from expert researchers at the forefront of their fields. Maplesoft proudly supports this workshop by giving teachers access to and training in Maple.

     

    Perimeter Institute also organizes a lively program of seminars, regularly exposing researchers and students to current ideas in the wider theoretical physics community. The talks provide content outside of, but related to, core disciplines.  Recently Maplesoft’s own physics expert Dr. Edgardo Cheb-Terrab conducted a lecture and training session at PI on Computer Algebra for Theoretical Physics.

    Dr. Edgardo Cheb-Terrab is the force behind the algorithms and Maple libraries of the ODE and PDE symbolic solvers, the MathematicalFunctions package (an expert system in special functions) and the Physics package, among other things.  

    In his talk at PI, the Physics project at Maplesoft was presented and the resulting Physics package was illustrated through simple problems in classical field theory, quantum mechanics and general relativity, and through tackling the computations of some recent Physical Review papers in those areas.   In addition there was a hands-on workshop where attendees were offered four choices of activity:  follow the mini-course; explore items of the worksheet of the morning presentation; or bring their own problems so that Dr. Cheb-Terrab could guide them on how to tackle it using the Physics package and Maple in general.

    As a company that strives to continuously improve student learning, and empower instructors and researchers with the tools necessary to compete in an ever changing and demanding educational environment, Maplesoft’s partnership with the Perimeter Institute allows us to do just that. We take great pride and joy in bringing our technology to outreach programs for students and teachers, making these opportunities a more productive and dynamic experience for all.

     

    In a previous Mapleprimes question related to Dirac Matrices, I was asked how to represent the algebra of Dirac matrices with an identity matrix on the right-hand side of  %AntiCommutator(Physics:-Dgamma[j], Physics:-Dgamma[k]) = 2*g[j, k]. Since this is a hot-topic in general, in that, making it work, involves easy and useful functionality however somewhat hidden, not known in general, it passed through my mind that this may be of interest in general. (To reproduce the computations below you need to update your Physics library with the one distributed at the Maplesoft R&D Physics webpage.)

     

    restart

    with(Physics)

     

    First of all, this shows the default algebra rules loaded when you load the Physics package, for the Pauli  and Dirac  matrices

    Library:-DefaultAlgebraRules()

    %Commutator(Physics:-Psigma[j], Physics:-Psigma[k]) = (2*I)*(Physics:-Psigma[1]*Physics:-LeviCivita[4, j, k, `~1`]+Physics:-Psigma[2]*Physics:-LeviCivita[4, j, k, `~2`]+Physics:-Psigma[3]*Physics:-LeviCivita[4, j, k, `~3`]), %AntiCommutator(Physics:-Psigma[j], Physics:-Psigma[k]) = 2*Physics:-KroneckerDelta[j, k], %AntiCommutator(Physics:-Dgamma[j], Physics:-Dgamma[k]) = 2*Physics:-g_[j, k]

    (1)

    Now, you can always overwrite these algebra rules.

     

    For instance, to represent the algebra of Dirac matrices with an identity matrix on the right-hand side, one can proceed as follows.

    First create the identity matrix. To emulate what we do with paper and pencil, where we write I to represent an identity matrix without having to see the actual table 2x2 with the number 1 in the diagonal and a bunch of 0, I will use the old matrix command, not the new Matrix (see more comments on this at the end). One way of entering this identity matrix is

    `𝕀` := matrix(4, 4, proc (i, j) options operator, arrow; KroneckerDelta[i, j] end proc)

    array( 1 .. 4, 1 .. 4, [( 4, 1 ) = (0), ( 1, 2 ) = (0), ( 2, 3 ) = (0), ( 1, 3 ) = (0), ( 2, 2 ) = (1), ( 4, 2 ) = (0), ( 3, 4 ) = (0), ( 1, 4 ) = (0), ( 3, 1 ) = (0), ( 4, 4 ) = (1), ( 3, 2 ) = (0), ( 1, 1 ) = (1), ( 2, 1 ) = (0), ( 4, 3 ) = (0), ( 3, 3 ) = (1), ( 2, 4 ) = (0)  ] )

    (2)

    The most important advantage of the old matrix command is that I is of type algebraic and, consequently, this is the important thing, one can operate with it algebraically and its contents are not displayed:

    type(`𝕀`, algebraic)

    true

    (3)

    `𝕀`

    `𝕀`

    (4)

    And so, most commands of the Maple library, that only work with objects of type algebraic, will handle the task. The contents are displayed only on demand, for instance using eval

    eval(`𝕀`)

    array( 1 .. 4, 1 .. 4, [( 4, 1 ) = (0), ( 1, 2 ) = (0), ( 2, 3 ) = (0), ( 1, 3 ) = (0), ( 2, 2 ) = (1), ( 4, 2 ) = (0), ( 3, 4 ) = (0), ( 1, 4 ) = (0), ( 3, 1 ) = (0), ( 4, 4 ) = (1), ( 3, 2 ) = (0), ( 1, 1 ) = (1), ( 2, 1 ) = (0), ( 4, 3 ) = (0), ( 3, 3 ) = (1), ( 2, 4 ) = (0)  ] )

    (5)

    Returning to the topic at hands: set now the algebra the way you want, with an I matrix on the right-hand side, and without seeing a bunch of 0 and 1

    %AntiCommutator(Dgamma[mu], Dgamma[nu]) = 2*g_[mu, nu]*`𝕀`

    %AntiCommutator(Physics:-Dgamma[mu], Physics:-Dgamma[nu]) = 2*Physics:-g_[mu, nu]*`𝕀`

    (6)

    Setup(algebrarules = (%AntiCommutator(Physics[Dgamma][mu], Physics[Dgamma][nu]) = 2*Physics[g_][mu, nu]*`𝕀`))

    [algebrarules = {%AntiCommutator(Physics:-Dgamma[mu], Physics:-Dgamma[nu]) = 2*Physics:-g_[mu, nu]*`𝕀`}]

    (7)

    And that is all.

     

    Check it out

    (%AntiCommutator = AntiCommutator)(Dgamma[mu], Dgamma[nu])

    %AntiCommutator(Physics:-Dgamma[mu], Physics:-Dgamma[nu]) = 2*Physics:-g_[mu, nu]*`𝕀`

    (8)

    Set now a Dirac spinor; this is how you could do that, step-by-step.

     

    Again you can use {vector, matrix, array} or {Vector, Matrix, Array}, and again, if you use the Upper case commands, you always have the components visible, and cannot compute with the object normally since they are not of type algebraic. So I use matrix, not Matrix, and matrix instead of vector so that the Dirac spinor that is both algebraic and matrix, is also displayed in the usual display as a "column vector"

     

    _local(Psi)

    Setup(anticommutativeprefix = {Psi, psi})

    [anticommutativeprefix = {_lambda, psi, :-Psi}]

    (9)

    In addition, following your question, in this example I explicitly specify the components of the spinor, in any preferred way, for example here I use psi[j]

    Psi := matrix(4, 1, [psi[1], psi[2], psi[3], psi[4]])

    array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] )

    (10)

    Check it out:

    Psi

    Psi

    (11)

    type(Psi, algebraic)

    true

    (12)

    Let's see all this working together by multiplying the anticommutator equation by Psi

    (%AntiCommutator(Physics[Dgamma][mu], Physics[Dgamma][nu]) = 2*Physics[g_][mu, nu]*`𝕀`)*Psi

    Physics:-`*`(%AntiCommutator(Physics:-Dgamma[mu], Physics:-Dgamma[nu]), Psi) = 2*Physics:-g_[mu, nu]*Physics:-`*`(`𝕀`, Psi)

    (13)

    Suppose now that you want to see the matrix form of this equation

    Library:-RewriteInMatrixForm(Physics[`*`](%AntiCommutator(Physics[Dgamma][mu], Physics[Dgamma][nu]), Psi) = 2*Physics[g_][mu, nu]*Physics[`*`](`𝕀`, Psi))

    Physics:-`.`(%AntiCommutator(Physics:-Dgamma[mu], Physics:-Dgamma[nu]), array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] )) = 2*Physics:-g_[mu, nu]*Physics:-`.`(array( 1 .. 4, 1 .. 4, [( 4, 1 ) = (0), ( 1, 2 ) = (0), ( 2, 3 ) = (0), ( 1, 3 ) = (0), ( 2, 2 ) = (1), ( 4, 2 ) = (0), ( 3, 4 ) = (0), ( 1, 4 ) = (0), ( 3, 1 ) = (0), ( 4, 4 ) = (1), ( 3, 2 ) = (0), ( 1, 1 ) = (1), ( 2, 1 ) = (0), ( 4, 3 ) = (0), ( 3, 3 ) = (1), ( 2, 4 ) = (0)  ] ), array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] ))

    (14)

    The above has the matricial operations delayed; unleash them

    %

    Physics:-`.`(%AntiCommutator(Physics:-Dgamma[mu], Physics:-Dgamma[nu]), array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] )) = 2*Physics:-g_[mu, nu]*(array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] ))

    (15)

    Or directly perform in one go the matrix operations behind (13)

    Library:-PerformMatrixOperations(Physics[`*`](%AntiCommutator(Physics[Dgamma][mu], Physics[Dgamma][nu]), Psi) = 2*Physics[g_][mu, nu]*Physics[`*`](`𝕀`, Psi))

    Physics:-`.`(%AntiCommutator(Physics:-Dgamma[mu], Physics:-Dgamma[nu]), array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] )) = 2*Physics:-g_[mu, nu]*(array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] ))

    (16)

    REMARK: As shown above, in general, the representation using lowercase commands allows you to use `*` or `.` depending on whether you want to represent the operation or perform the operation. For example this represents the operation, as an exact mimicry of what we do with paper and pencil, both regarding input and output

    `𝕀`*Psi

    Physics:-`*`(`𝕀`, Psi)

    (17)

    And this performs the operation

    `𝕀`.Psi

    array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] )

    (18)

    Or to only displaying the operation

    Library:-RewriteInMatrixForm(Physics[`*`](`𝕀`, Psi))

    Physics:-`.`(array( 1 .. 4, 1 .. 4, [( 4, 1 ) = (0), ( 1, 2 ) = (0), ( 2, 3 ) = (0), ( 1, 3 ) = (0), ( 2, 2 ) = (1), ( 4, 2 ) = (0), ( 3, 4 ) = (0), ( 1, 4 ) = (0), ( 3, 1 ) = (0), ( 4, 4 ) = (1), ( 3, 2 ) = (0), ( 1, 1 ) = (1), ( 2, 1 ) = (0), ( 4, 3 ) = (0), ( 3, 3 ) = (1), ( 2, 4 ) = (0)  ] ), array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] ))

    (19)

    And to perform all these matricial operations within an algebraic expression,

    Library:-PerformMatrixOperations(Physics[`*`](`𝕀`, Psi))

    Matrix(%id = 18446744079185513758)

    (20)

    ``

     


     

    Download DiracAlgebraWithIdentityMatrix.mw

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

    Hey guys...

    One of my favorite new images rendered in Maple.....

    integrates mod congruency arguments like the following, into an 8 stage animated 3D cylindrical coordinates context::

     

    u*mod(x^(1/3)cos(x),23)

     

     

    As a powerlifter, I constantly find myself calculating my total between my competition lifts, bench, squat, and deadlift. Following that, I always end up calculating my wilks score. The wilks score is a score used to compare lifters between weight classes (https://en.wikipedia.org/wiki/Wilks_Coefficient ), as comparing someone who weighs 59kg in competition like me, to someone who weighs 120kg+, the other end of the spectrum; obviously the 120kg+ lifter is going to massively out-lift me.

    So I decided to program a wilks calculator for quick use, rather than needing to go search for one on the internet. For anyone curious about specific scoring, a score of 300+ is very strong for the average gym goer, and is about normal for a powerlifter. A score of 400+ makes you strong for a powerlifter, putting you in the top 250 powerlifters, while 500+ is the very top, as far as unequipped powerlifting goes, including the top 30. For anyone wondering, my best score at my best meet was 390, although given the progress I’ve made in the gym, should be above 400 by my next meet.

    Hope you all enjoy!

     Find it here: https://maple.cloud#doc=5687076260413440&key=301A440EFD2C4EDD8480D60B5E3147BF40CA460F842942449C939AB8D2E7D679

    The problem is to analyze the behavior of an involute of the cubical parabola near the singular points and at infinity.

    To do this efficiently, it is best to work with partially evaluated expressions.

    We want to investigate the singular points of an involute of the cubical parabola y = x^3.

    If the curve is given parametrically by "conjugate(r)(t)", an involute is given by

    "conjugate(E)(t)=conjugate(r)(t)-(conjugate(r)'(t))/(|conjugate(r)'(t)|)*(∫)[a]^(t)|conjugate(r)'(tau)| ⅆtau"

    What we want to do is to leave the integral unevaluated and compute only its derivatives. The simplest way is to define S := proc (a, t) options operator, arrow; Int(s(tau), tau = a .. t) end proc. Both diff and evalf will be handled automatically by Maple. We'll do it in a slightly more complicated way, leaving S(a, t) completely inert with the exception of the condition S(a, a)=0 and implementing a more efficient numerical evaluator.

    forget(diff); forget(series); r := proc (t) options operator, arrow; [t, t^3] end proc; s := unapply(sqrt((diff(r(t)[1], t))^2+(diff(r(t)[2], t))^2), t); S := proc (a, t) options operator, arrow; `if`(t = a, 0, 'S(a, t)') end proc; `diff/S` := proc (at, ft, t) options operator, arrow; s(ft)*(diff(ft, t))-s(at)*(diff(at, t)) end proc; E := unapply(r(t)-`~`[`*`](diff(r(t), t), S(a, t)/s(t)), [a, t])

    forget(evalf); `evalf/S` := (proc () local acur, Ssol; acur := undefined; Ssol := subs(dsolve({diff(f(t), t) = s(t), f(a) = 0}, f(t), numeric, parameters = [a], output = listprocedure), f(t)); proc (a, t) if [a, t]::(list(numeric)) then if a <> acur then Ssol(parameters = [a]); acur := a end if; Ssol(t) else 'S(a, t)' end if end proc end proc)()

    plot([[op(r(t)), t = -1 .. 3], [op(E(1, t)), t = -1.5 .. 4]], view = -1 .. 3, scaling = constrained)

     

    simplify(diff(E(a, t), t))

    [18*S(a, t)*t^3/(9*t^4+1)^(3/2), -6*t*S(a, t)/(9*t^4+1)^(3/2)]

    (1)

    The singular points are easily determined now: we need either S(a, t)=0, which implies t=a, or t=0. We assume a>0.

    Compute the Taylor series around t=a first.

    `~`[series](E(a, t+a), t = 0, 4); ser := collect(convert(%, polynom), t, normal)

    [-18*a^2*(6*a^4-1)*t^3/(9*a^4+1)^2+9*a^3*t^2/(9*a^4+1)+a, 2*(36*a^4-1)*t^3/(9*a^4+1)^2-3*a*t^2/(9*a^4+1)+a^3]

    (2)

    The center of the expansion E(a, a) is the point [a, a^3] on the cubical parabola.

    There are two quadratic terms, which give us the slope of the tangent line. The tangent coincides with the normal to the cubical parabola at this point.

    Rotate the axes to make the tangent vertical.

    rot := proc (e, a) options operator, arrow; convert(Student:-LinearAlgebra:-RotationMatrix(a).Vector(e), list) end proc

    rot(ser-E(a, a), (`@`(arctan, op))(`~`[coeff](ser, t, 2))); ser := `assuming`([collect(%, t, normal)], [a > 0])

    [-12*a^2*t^3/(9*a^4+1)^(3/2), -2*(162*a^8+9*a^4-1)*t^3/(9*a^4+1)^(5/2)+3*a*t^2/(9*a^4+1)^(1/2)]

    (3)

    In the new coordinates, the leading terms are t^3 and t^2, and the involute has the shape of a semicubical parabola. When t increases, the involute moves from the first to the second quadrant. In terms of x and y, `~`[x]*alpha*y^(3/2), and the coefficient alpha is found from the two leading terms. For t>0:

    `assuming`([simplify(coeff(ser[1], t, 3)/coeff(ser[2], t, 2)^(3/2))], [a > 0])

    -(4/3)*a^(1/2)*3^(1/2)/(9*a^4+1)^(3/4)

    (4)

    Now consider the second singular point t=0 in the original coordinate system.

    ser := convert(`~`[series](E(a, t), t = 0, 6), polynom)

    [-S(a, 0)+(9/2)*S(a, 0)*t^4+(18/5)*t^5, -3*S(a, 0)*t^2-2*t^3]

    (5)

    The center is the point [-S(a, 0), 0], which corresponds to the point [0, 0] on the cubical parabola. The leading terms t^4 and t^2 indicate that the tangent line is vertical and that this is a return point, with both branches lying in the second quadrant in the coordinate system shifted by [-S(a, 0), 0].

    In terms of x and y, `~`[x]*alpha*y^2, and to determine the next term, we need to analyze the two series expansions. t^4 gives a t^4 term and fixes alpha. t^2*(t^3) gives a t^5 term, but the coefficient doesn't match the coefficient at t^5 in x. So the next term has to be y^(5/2) in order to match t^5. For t>0, (b*t^3+a*t^2)^(5/2) = t^5*(b*t+a)^(5/2), which gives a^(5/2)*t^5 plus higher order terms. No other terms contribute, and we can equate the coefficients at t^4 and t^5 in x and y:

    collect(alpha*ser[2]^2+beta*coeff(ser[2], t, 2)^(5/2)*t^5, t); (`@`(simplify, solve))(`~`[coeff](%-ser[1], t, [4, 5]), {alpha, beta})

    4*alpha*t^6+(beta*(-3*S(a, 0))^(5/2)+12*alpha*S(a, 0))*t^5+9*alpha*S(a, 0)^2*t^4

     

    {alpha = (1/2)/S(a, 0), beta = -(4/45)*3^(1/2)/(-S(a, 0))^(5/2)}

    (6)

    We have obtained the coefficients in `~`[x]*alpha*y^2+beta*y^(5/2). alpha and beta are negative, while for t<0, beta has the opposite sign, and the branch corresponding to t<0 is the one which is closer to the vertical tangent.

    This could have been done by using Groebner:-Basis to eliminate t from the two series and then using algcurves:-puiseux. But it is still necessary to determine the truncation order.

    There is also a separate case a=0.

    ser := `~`[series](E(0, t), t = 0, 6)

    [series((18/5)*t^5+O(t^6),t,6), series(-2*t^3+O(t^7),t,7)]

    (7)

    Now the singularity is an inflection point, and the involute has a vertical tangent and moves from the second to the fourth quadrant. In the fourth quadrant, `~`[x]*alpha*(-y)^(5/3), and we find alpha in the same way as before:

    coeff(ser[1], t, 5)/(-coeff(ser[2], t, 3))^(5/3)

    (9/10)*2^(1/3)

    (8)

    Finally let's consider the asymptotic behavior of an involute at infinity. Now we do have to investigate S(a, t). One possible way is to apply the binomial formula to s(t) = 3*t^2*sqrt(1+1/(9*t^4))  and integrate term by term:

    `assuming`([int(3*tau^2*binomial(1/2, k)*(9*tau^4)^(-k), tau = a .. t)], [0 < a and a < t])

    3*binomial(1/2, k)*(9^(-k)*a^(-4*k+3)-9^(-k)*t^(-4*k+3))/(4*k-3)

    (9)

    This is the power series expansion of S(a, t) for large t, valid for 9*a^4 > 1. Substituting it into E(a, t), we obtain the required asymptotics.

    E(a, t)

    [-S(a, t)/(9*t^4+1)^(1/2)+t, -3*t^2*S(a, t)/(9*t^4+1)^(1/2)+t^3]

    (10)

    In the x component, we get -(1/3)*t+o(1)+t. In the y component, we get -t^3-C+o(1)+t^3. Thus the involute approaches a horizontal asymptote when t goes to +infinity, and the constant term -C gives the position of the asymptote:

    `C__+` := `assuming`([unapply(sum(3*binomial(1/2, k)*9^(-k)*a^(-4*k+3)/(4*k-3), k = 0 .. infinity), a)], [a > 1/sqrt(3)])

    proc (a) options operator, arrow; -a^3*hypergeom([-3/4, -1/2], [1/4], -(1/9)/a^4) end proc

    (11)

    That in fact extends to all positive a. For the asymptotics at -infinity, first we compute the integral over [a, -a], this time using the power series for small t. For the integral from -a to -t, S(-a, -t) = -S(a, t), and we need to subtract the value found above:

    `C__-` := `assuming`([unapply(sum(int(binomial(1/2, k)*(9*tau^4)^k, tau = a .. -a), k = 0 .. infinity)-`C__+`(a), a)], [0 < a and a < 1/sqrt(3)])

    proc (a) options operator, arrow; -2*a*hypergeom([-1/2, 1/4], [5/4], -9*a^4)+a^3*hypergeom([-3/4, -1/2], [1/4], -(1/9)/a^4) end proc

    (12)

    For negative a, `C__+` and `C__-` are swapped.

    This could have been done by converting S(a, t) into a difference of two integrals of hypergeometric type, for which Maple is able to find the asymptotics automatically.

    inv := proc (a, T) plots:-display(plot([op(r(t)), t = -1 .. 2], color = black), plot([op(E(a, t)), t = -1.5 .. T], color = red), plottools:-line(r(T), E(a, T), linestyle = dot), map(proc (y) options operator, arrow; plottools:-line([-1, y], [3, y], linestyle = spacedot) end proc, [-`C__+`(a), -`C__-`(a)]), scaling = constrained, view = [-1 .. 3, -1 .. 3]) end proc

    Explore(inv(a, t), a = -.3 .. 1.3, t = -1. .. 4., initialvalues = [a = 1., t = 3.])

    ``

    Download cubic.mw

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

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

    Indeed,

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

    Bug_in_MatrixInverse.mw

    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?

    Moreover,

    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.

    Bug_in_pdsolve.mw

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

    (1)

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

    (2)

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

    (3)

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

    [(4779/36980000000000)*(4341503*CO^2+29548445000000000000*CO-209981000000000000000000000000000)*NO2r/(529*CO^2-4600000000000000*CO+10000000000000000000000000000)]

    (4)

    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)

    [(4779/36980000000000)*(4341503*CO^2+29548445000000000000*CO-209981000000000000000000000000000)/(529*CO^2-4600000000000000*CO+10000000000000000000000000000)]

    (5)

    ``

    Download np.mw


     

    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

    NULL

    NULL

    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 `&hbar;`(&PartialD;)/(&PartialD;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*`&hbar;`*`&PartialD;`[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*`&hbar;`*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]

    (1.1)

    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("&rarr;"))`, and from there get the commutation rule between F(X) and `#mover(mi("p",mathcolor = "olive"),mo("&rarr;"))`.

     

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

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

    CompactDisplay(F(X))

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

    (1.2)

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

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

    (1.3)

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

    (1.4)

    Note that  `#mover(mi("p",mathcolor = "olive"),mo("&rarr;"))`, 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("&rarr;"))`. 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)))

    (1.5)

    Note that, in `#mover(mi("p",mathcolor = "olive"),mo("&rarr;"))`(F(X)*Ket(psi, x, y, z)), the application of `#mover(mi("p",mathcolor = "olive"),mo("&rarr;"))` is not expanded: at this point nothing is known about  `#mover(mi("p",mathcolor = "olive"),mo("&rarr;"))` , 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("&rarr;"))` 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*`&hbar;`*%Nabla(f) end proc

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

    (1.6)

    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*`&hbar;`*Physics:-`*`(F(X), %Nabla(Physics:-Ket(psi, x, y, z)))+I*`&hbar;`*%Nabla(Physics:-`*`(F(X), Physics:-Ket(psi, x, y, z)))

    (1.7)

    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*`&hbar;`*Physics[`*`](F(X), %Nabla(Physics[Ket](psi, x, y, z)))+I*`&hbar;`*%Nabla(Physics[`*`](F(X), Physics[Ket](psi, x, y, z)))))

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

    (1.8)

    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*`&hbar;`*_i*Physics[`*`](diff(F(X), x), Physics[Ket](psi, x, y, z))+I*`&hbar;`*_j*Physics[`*`](diff(F(X), y), Physics[Ket](psi, x, y, z))+I*`&hbar;`*_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*`&hbar;`*Physics:-`*`((diff(F(X), y))*_j+(diff(F(X), z))*_k+(diff(F(X), x))*_i, Physics:-Ket(psi, x, y, z))

    (1.9)

    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*`&hbar;`*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*`&hbar;`*Physics:-`*`(%Gradient(F(X)), Physics:-Ket(psi, x, y, z))

    (1.10)

    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("&rarr;"))`

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

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

    (1.11)

    Tensor notation, "[`X__m`,P[n]][-]=i `&hbar;` 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]

    (2.1)

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

    (2.2)

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

    (2.3)

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

    (2.4)

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

    (2.5)

    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*`&hbar;`*%d_[n](f) end proc

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

    (2.6)

    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*`&hbar;`*Physics:-`*`(Physics:-SpaceTimeVector[m](X), %d_[n](Physics:-Ket(psi, x, y, z)))+I*`&hbar;`*%d_[n](Physics:-`*`(Physics:-SpaceTimeVector[m](X), Physics:-Ket(psi, x, y, z)))

    (2.7)

    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*`&hbar;`*Physics[`*`](Physics[SpaceTimeVector][m](X), %d_[n](Ket(psi, x, y, z)))+I*`&hbar;`*%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*`&hbar;`*Physics:-g_[m, n]*Physics:-Ket(psi, x, y, z)

    (2.8)

    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*`&hbar;`*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*`&hbar;`*Physics:-g_[m, n]

    (2.9)

    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*`&hbar;`*Physics[g_][m, n], [P[n], X[m]], usecommutator)

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

    (2.10)

    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*`&hbar;`*Physics[g_][m, n]))

    Matrix(%id = 18446744078261558678)

    (2.11)

    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("&rarr;"))` , 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)}

    (2.12)

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

    Matrix(%id = 18446744078575996430)

    (2.13)

    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

    "Setup(?)"

    [algebrarules = {%Commutator(x, p__x) = I*`&hbar;`, %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*`&hbar;`, %Commutator(y, p__z) = 0, %Commutator(y, z) = 0, %Commutator(z, p__x) = 0, %Commutator(z, p__y) = 0, %Commutator(z, p__z) = I*`&hbar;`}]

    (2.14)

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

    (%Commutator = Commutator)(x, p__x)

    %Commutator(x, p__x) = I*`&hbar;`

    (2.15)

    NULL

    NULL


     

    Download DifferentialOperatorCommutatorRules.mw

     

    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: https://www.maplesoft.com/applications/view.aspx?SID=154353 

    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:

    with(OEIS):

     

    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):
    indices(results81);
        ["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:

    results81["author"];
        "_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:

    PackageTools:-Install(5693538177122304);

     

    Comments and feedback are welcome here, or on the project page: https://github.com/dskoog/Maple-OEIS

    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()
    export
        Parser := Bark:-ArgParser(NULL
                                  , 'prologue' = ( "Print `Hello, World!'" )
                                  , 'opts' = ['help' :: 'help' &c "Print this help page"]
                                 );
    export
        ModuleApply := proc(cmdline :: string := "")
            Bark:-ArgParser:-Parse(Parser, cmdline);
            Bark:-printf("Hello, World!\n");
            NULL;
        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
    MAPLE='/home/joe/maplesoft/sandbox/main/bin/maple'
    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.

    PackageTools:-Install~([5741316844552192,6273820789833728],'overwrite'):
    Bark:-InstallExamples('overwrite'):
    

    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.

    First 45 46 47 48 49 50 51 Last Page 47 of 301