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
  • The 196 algorithm goes like this.  Start with an integer.  Reverse the digits.  Add the reversed number to the integer.  For most numbers, this eventually leads to a palendrome.  That is to say the number is equal to the reversed number.  I wrote a little Maple procedure to explore 196, the smallest number that will probrably never become a palendrome when put into the algorithm.


    Let me know if you like my code.




    The material below was presented in the "Semantic Representation of Mathematical Knowledge Workshop", February 3-5, 2016 at the Fields Institute, University of Toronto. It shows the approach I used for “digitizing mathematical knowledge" regarding Differential Equations, Special Functions and Solutions to Einstein's equations. While for these areas using databases of information helps (for example textbooks frequently contain these sort of databases), these are areas that, at the same time, are very suitable for using algorithmic mathematical approaches, that result in much richer mathematics than what can be hard-coded into a database. The material also focuses on an interesting cherry-picked collection of Maple functionality, that I think is beautiful, not well know, and seldom focused inter-related as here.



    Digitizing of special functions,

    differential equations,

    and solutions to Einstein’s equations

    within a computer algebra system


    Edgardo S. Cheb-Terrab

    Physics, Differential Equations and Mathematical Functions, Maplesoft

    Editor, Computer Physics Communications



    Digitizing (old paradigm)



    Big amounts of knowledge available to everybody in local machines or through the internet


    Take advantage of basic computer functionality, like searching and editing



    Digitizing (new paradigm)


    By digitizing mathematical knowledge inside appropriate computational contexts that understand about the topics, one can use the digitized knowledge to automatically generate more and higher level knowledge




    1) how to identify, test and organize the key blocks of information,


    2) how to access it: the interface,


    3) how to mathematically process it to automatically obtain more information on demand





                                               Three examples

    Mathematical Functions


    "Mathematical functions, are defined by algebraic expressions. So consider algebraic expressions in general ..."

    The FunctionAdvisor (basic)


    "Supporting information on definitions, identities, possible simplifications, integral forms, different types of series expansions, and mathematical properties in general"



    General description





    Differential equation representation for generic nonlinear algebraic expressions - their use


    "Compute differential polynomial forms for arbitrary systems of non-polynomial equations ..."

    The Differential Equations representing arbitrary algebraic expresssions


    Deriving knowledge: ODE solving methods


    Extending the mathematical language to include the inverse functions


    Solving non-polynomial algebraic equations by solving polynomial differential equations





    Branch Cuts of algebraic expressions


    "Algebraically compute, and visualize, the branch cuts of arbitrary mathematical expressions"






    Algebraic expresssions in terms of specified functions


    "A conversion network for arbitrary mathematical expressions, to rewrite them in terms of different functions in flexible ways"



    General description





    Symbolic differentiation of algebraic expressions


    "Perform symbolic differentiation by combining different algebraic techniques, including functions of symbolic sequences and Faà di Bruno's formula"






    Ordinary Differential Equations


    "Beyond the concept of a database, classify an arbitrary ODE and suggest solution methods for it"

    General description







    Exact Solutions to Einstein's equations



    Lambda*g[mu, nu]+G[mu, nu] = 8*Pi*T[mu, nu]


    "The authors of "Exact solutions toEinstein's equations" reviewed more than 4,000 papers containing solutions to Einstein’s equations in the general relativity literature, organized the whole material into chapters according to the physical properties of these solutions. These solutions are key in the area of general relativity, are now all digitized and become alive in a worksheet"

    The ability to search the database according to the physical properties of the solutions, their classification, or just by parts of keywords (old paradigm) changes the game.

    More important, within a computer algebra system this knowledge becomes alive (new paradigm).


    The solutions are turned active by a simple call to one commend, called the g_  spacetime metric.


    Everything else gets automatically derived and set on the fly ( Christoffel symbols  , Ricci  and Riemann  tensors orthonormal and null tetrads , etc.)


    Almost all of the mathematical operations one can perform on these solutions are implemented as commands in the Physics  and DifferentialGeometry  packages.


    All the mathematics within the Maple library are instantly ready to work with these solutions and derived mathematical objects.


    Finally, in the Maple PDEtools package , we have all the mathematical tools to tackle the equivalence problem around these solutions.






    Download:,    Digitizing_Mathematical_Information.pdf

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

    A string is wound symmetrically around a circular rod. The string goes exactly
    4 times around the rod. The circumference of the rod is 4 cm and its length is 12 cm.
    Find the length of the string.
    Show all your work.

    (It was presented at a meeting of the European Mathematical Society in 2001,
    "Reference levels in mathematics in Europe at age16").

    Can you solve it? You may want to try before seing the solution.
    [I sometimes train olympiad students at my university, so I like such problems].

    eq:= 2/Pi*cos(t), 2/Pi*sin(t), 3/2/Pi*t; # The equations of the helix, t in 0 .. 8*Pi:
    p:=plots[spacecurve]([eq, t=0..8*Pi],scaling=constrained,color=red, thickness=5, axes=none):
    plots:-display(plottools:-cylinder([0,0,0], 2/Pi, 12, style=surface, color=yellow),
                             p, scaling=constrained,axes=none);

    VectorCalculus:-ArcLength(<eq>, t=0..8*Pi);



    Let's look at the first loop around the rod.
    If we develop the corresponding 1/4 of the cylinder, it results a rectangle  whose sides are 4 and 12/4 = 3.
    The diagonal is 5 (ask Pythagora why), so the length of the string is 4*5 = 20.


    This presentation is on an undergrad intermediate Quantum Mechanics topic. Tackling the problem within a computer algebra worksheet in the way shown below is actually the novelty, using the Physics package to formulate the problem with quantum operators and related algebra rules in tensor notation.


    Quantization of the Lorentz Force


    Pascal Szriftgiser1 and Edgardo S. Cheb-Terrab2 

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

    (2) Maplesoft


    We consider the case of a quantum, non-relativistic, particle with mass m and charge q evolving under the action of an arbitrary time-independent magnetic field "B=Curl(A(x,y,z)), "where `#mover(mi("A",mathcolor = "olive"),mo("&rarr;"))` is the vector potential. The Hamiltonian for this system is

    H = (`#mover(mi("p",mathcolor = "olive"),mo("&rarr;"))`-q*`#mover(mi("A",mathcolor = "olive"),mo("&rarr;"))`(X))^2/(2*m)

    where `#mover(mi("p",mathcolor = "olive"),mo("&rarr;"))` is the momentum of the particle, and the force acting in this particle, also called the Lorentz force, is given by


    `#mover(mi("F",mathcolor = "olive"),mo("&rarr;"))` = m*(diff(v(t), t))


    where `#mover(mi("v",mathcolor = "olive"),mo("&rarr;"))` is the quantized velocity of the particle, and all of  H, `#mover(mi("p",mathcolor = "olive"),mo("&rarr;"))`, `#mover(mi("v",mathcolor = "olive"),mo("&rarr;"))`, `#mover(mi("B",mathcolor = "olive"),mo("&rarr;"))`, `#mover(mi("A",mathcolor = "olive"),mo("&rarr;"))` and `#mover(mi("F",mathcolor = "olive"),mo("&rarr;"))` are Hermitian quantum operators representing observable quantities.


    In the classic (non-quantum) case, `#mover(mi("F"),mo("&rarr;"))` for such a particle in the absence of electrical field is given by


    `#mover(mi("F"),mo("&rarr;"))` = `&x`(q*`#mover(mi("v"),mo("&rarr;"))`, `#mover(mi("B"),mo("&rarr;"))`) ,


    Problem: Departing from the Hamiltonian, show that in the quantum case the Lorentz force is given by [1]


    `#mover(mi("F",mathcolor = "olive"),mo("&rarr;"))` = (1/2)*q*(`&x`(`#mover(mi("v",mathcolor = "olive"),mo("&rarr;"))`, `#mover(mi("B",mathcolor = "olive"),mo("&rarr;"))`)-`&x`(`#mover(mi("B",mathcolor = "olive"),mo("&rarr;"))`, `#mover(mi("v",mathcolor = "olive"),mo("&rarr;"))`))


    [1] Photons et atomes, Introduction à l'électrodynamique quantique, p. 179, Claude Cohen-Tannoudji, Jacques Dupont-Roc et Gilbert Grynberg - EDP Sciences janvier 1987.




    We choose to tackle the problem in Heisenberg's picture of quantum mechanices, where the state of a system is static and only the quantum operators evolve in time according to


    diff(O(t), t) = I*Physics:-Commutator(H, O(t))/`&hbar;`


    Also, the algebraic manipulations are simpler using tensor abstract notation instead of the standard 3D vector notation. We then start setting the framework for the problem, a system of coordinates X, indicating the dimension of the tensor space to be 3 and the metric Euclidean, and that we will use lowercaselatin letters to represent tensor indices. In addition, not necessary but for convenience, we set the lowercase latin i to represent the imaginary unit and we request automaticsimplification so that the output of everything comes automatically simplified in size.


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

    Setup(mathematicalnotation = true, automaticsimplification = true, coordinates = X, dimension = 3, metric = Euclidean, spacetimeindices = lowercaselatin, quiet)

    [automaticsimplification = true, coordinatesystems = {X}, dimension = 3, mathematicalnotation = true, metric = {(1, 1) = 1, (2, 2) = 1, (3, 3) = 1}, spacetimeindices = lowercaselatin]



    Next we indicate the letters we will use to represent the quantum operators with which we will work, and also the standard commutation rules between position and momentum, always the starting point when dealing with quantum mechanics problems


    Setup(quantumoperators = {F}, hermitianoperators = {A, B, H, p, r, v, x}, realobjects = {`&hbar;`, m, q}, algebrarules = {%Commutator(p[k], p[n]) = 0, %Commutator(x[k], p[l]) = I*`&hbar;`*KroneckerDelta[k, l], %Commutator(x[k], x[l]) = 0})

    [algebrarules = {%Commutator(p[k], p[n]) = 0, %Commutator(x[k], p[l]) = I*`&hbar;`*Physics:-KroneckerDelta[k, l], %Commutator(x[k], x[l]) = 0}, hermitianoperators = {A, B, H, p, r, v, x}, quantumoperators = {A, B, F, H, p, r, v, x}, realobjects = {`&hbar;`, m, q, x1, x2, x3, %dAlembertian, Physics:-dAlembertian}]



    Note that we start not indicating F as Hermitian, in order to arrive at that result. The quantum operators A, B, and F are explicit functions of X, so to avoid redundant display of this functionality on the screen we use


    CompactDisplay((A, B, F)(X))

    A(x1, x2, x3)*`will now be displayed as`*A


    B(x1, x2, x3)*`will now be displayed as`*B


    F(x1, x2, x3)*`will now be displayed as`*F


    Define now as tensors the quantum operators that we will use with tensorial notation (recalling: for these, Einstein's sum rule for repeated indices will be automatically applied when simplifying)


    Define(x, p, v, A, B, F, quiet)

    {A, B, F, p, v, x, Physics:-Dgamma[a], Physics:-Psigma[a], Physics:-d_[a], Physics:-g_[a, b], Physics:-KroneckerDelta[a, b], Physics:-LeviCivita[a, b, c], Physics:-SpaceTimeVector[a](X)}


    The Hamiltonian,

    H = (`#mover(mi("p",mathcolor = "olive"),mo("&rarr;"))`-q*`#mover(mi("A",mathcolor = "olive"),mo("&rarr;"))`(X))^2/(2*m)

    in tensorial notation, is given by

    H = (p[n]-q*A[n](X))^2/(2*m)

    H = (1/2)*Physics:-`^`(p[n]-q*A[n](X), 2)/m


    Generally speaking to arrive at  ```#mover(mi("F",mathcolor = "olive"),mo("&rarr;"))` = (1/2)*q*(`&x`(`#mover(mi("v",mathcolor = "olive"),mo("&rarr;"))`, `#mover(mi("B",mathcolor = "olive"),mo("&rarr;"))`)-`&x`(`#mover(mi("B",mathcolor = "olive"),mo("&rarr;"))`, `#mover(mi("v",mathcolor = "olive"),mo("&rarr;"))`)) what we now need to do is

    1) Express this Hamiltonian (5) in terms of the velocity


    And, recalling that, in Heisenberg's picture, quantum operators evolve in time according to

    diff(O(t), t) = I*Physics:-Commutator(H, O(t))/`&hbar;`


    2) Take the commutator of H with the velocity itself to obtain its time derivative and, from `#mover(mi("F",mathcolor = "olive"),mo("&rarr;"))` = m*(diff(v(t), t)) , that commutator is already the force up to some constant factors.


    To get in contact with the basic commutation rules between position and momentum behind quantum phenomena, the quantized velocity itself can be computed as the time derivative of the position operator, i.e as the commutator of x[k] with H

    I*Commutator(H = (1/2)*Physics[`^`](p[n]-q*A[n](X), 2)/m, x[k])/`&hbar;`

    I*Physics:-Commutator(H, x[k])/`&hbar;` = (1/2)*(I*q^2*Physics:-AntiCommutator(A[n](X), Physics:-Commutator(A[n](X), x[k]))-I*q*Physics:-AntiCommutator(p[n], Physics:-Commutator(A[n](X), x[k]))-2*(q*A[n](X)-p[n])*Physics:-KroneckerDelta[k, n]*`&hbar;`)/(`&hbar;`*m)


    This expression for the velocity, that involves commutators between the potential A[n](X), the position x[k] and the momentum p[n], can be simplified taking into account the basic quantum algebra rules between position and momentum. We assume that A[n](X)(X) can be decomposed into a formal power series (possibly infinite) of the x[k], hence all the A[n](X) commute between themselves as well as with all the x[k]


    {%Commutator(A[k](X), x[l]) = 0, %Commutator(A[k](X), A[l](X)) = 0}

    {%Commutator(A[k](X), x[l]) = 0, %Commutator(A[k](X), A[l](X)) = 0}


    (Note: in some cases, this is not true, but those cases are beyond the scope of this worksheet.)


    Add these rules to the algebra rules already set so that they are all taken into account when simplifying things


    Setup(algebrarules = {%Commutator(A[k](X), x[l]) = 0, %Commutator(A[k](X), A[l](X)) = 0})

    [algebrarules = {%Commutator(p[k], p[n]) = 0, %Commutator(x[k], p[l]) = I*`&hbar;`*Physics:-KroneckerDelta[k, l], %Commutator(x[k], x[l]) = 0, %Commutator(A[k](X), x[l]) = 0, %Commutator(A[k](X), A[l](X)) = 0}]


    Simplify(I*Physics[Commutator](H, x[k])/`&hbar;` = (1/2)*(I*q^2*Physics[AntiCommutator](A[n](X), Physics[Commutator](A[n](X), x[k]))-I*q*Physics[AntiCommutator](p[n], Physics[Commutator](A[n](X), x[k]))-2*(q*A[n](X)-p[n])*Physics[KroneckerDelta][k, n]*`&hbar;`)/(`&hbar;`*m))

    I*Physics:-Commutator(H, x[k])/`&hbar;` = (-A[k](X)*q+p[k])/m


    The right-hand side of (9) is then the kth component of the velocity tensor quantum operator, the relationship is the same as in the classical case

    v[k] = rhs(I*Physics[Commutator](H, x[k])/`&hbar;` = (-A[k](X)*q+p[k])/m)

    v[k] = (-A[k](X)*q+p[k])/m


    and with this the Hamiltonian (5) can now be rewritten in term of the velocity completing step 1)

    simplify(H = (1/2)*Physics[`^`](p[n]-q*A[n](X), 2)/m, {SubstituteTensorIndices(k = n, (rhs = lhs)(v[k] = (-A[k](X)*q+p[k])/m))})

    H = (1/2)*m*Physics:-`^`(v[n], 2)


    For step 2), to compute

     `#mover(mi("F",mathcolor = "olive"),mo("&rarr;"))` = m*(diff(v(t), t)) and m*(diff(v(t), t)) = I*m*Physics:-Commutator(H, v(t)[k])/`&hbar;` 


    we need the commutator between the different components of the quantized velocity which, contrary to what happens in the classical case, do not commute. For this purpose, take the commutator between (10) with itself after replacing the free index

    Commutator(v[k] = (-A[k](X)*q+p[k])/m, SubstituteTensorIndices(k = n, v[k] = (-A[k](X)*q+p[k])/m))

    Physics:-Commutator(v[k], v[n]) = -q*(Physics:-Commutator(A[k](X), p[n])+Physics:-Commutator(p[k], A[n](X)))/m^2


    To simplify (12), we use the fact that if f  is a commutative mapping that can be decomposed into a formal power series in all the complex plan (which is assumed to be the case for all A[n](X)(X)), then

    Physics:-Commutator(p[k], f(x, y, z)) = -I*`&hbar;`*`&PartialD;`[k](f(x, y, z))

    where p[k]"=-i `&hbar;` `&PartialD;`[k] " is the momentum operator along the x[k] axis. This relation reads in tensor notation:

    Commutator(p[k], A[n](X)) = -I*`&hbar;`*d_[k](A[n](X))

    Physics:-Commutator(p[k], A[n](X)) = -I*`&hbar;`*Physics:-d_[k](A[n](X), [X])


    Add this rule to the rules previously set in order to automatically take it into account in (12)

    Setup(Physics[Commutator](p[k], A[n](X)) = -I*`&hbar;`*Physics[d_][k](A[n](X), [X]))

    [algebrarules = {%Commutator(p[k], p[n]) = 0, %Commutator(p[k], A[n](X)) = -I*`&hbar;`*Physics:-d_[k](A[n](X), [X]), %Commutator(x[k], p[l]) = I*`&hbar;`*Physics:-KroneckerDelta[k, l], %Commutator(x[k], x[l]) = 0, %Commutator(A[k](X), x[l]) = 0, %Commutator(A[k](X), A[l](X)) = 0}]


    Physics[Commutator](v[k], v[n]) = -q*(Physics[Commutator](A[k](X), p[n])+Physics[Commutator](p[k], A[n](X)))/m^2

    Physics:-Commutator(v[k], v[n]) = -I*q*`&hbar;`*(Physics:-d_[n](A[k](X), [X])-Physics:-d_[k](A[n](X), [X]))/m^2


    Also add this other rule so that it is taken into account automatically

    Setup(Physics[Commutator](v[k], v[n]) = -I*q*`&hbar;`*(Physics[d_][n](A[k](X), [X])-Physics[d_][k](A[n](X), [X]))/m^2)

    [algebrarules = {%Commutator(p[k], p[n]) = 0, %Commutator(p[k], A[n](X)) = -I*`&hbar;`*Physics:-d_[k](A[n](X), [X]), %Commutator(v[k], v[n]) = -I*q*`&hbar;`*(Physics:-d_[n](A[k](X), [X])-Physics:-d_[k](A[n](X), [X]))/m^2, %Commutator(x[k], p[l]) = I*`&hbar;`*Physics:-KroneckerDelta[k, l], %Commutator(x[k], x[l]) = 0, %Commutator(A[k](X), x[l]) = 0, %Commutator(A[k](X), A[l](X)) = 0}]


    Recalling now the expression of the Hamiltonian (11) as a function of the velocity, one can compute the components of the force operator  "()Component(v*B,k)=m (v[k])=(i m [H,v[k]][-])/`&hbar;`"

    F[k](X) = I*m*%Commutator(rhs(H = (1/2)*m*Physics[`^`](v[n], 2)), v[k])/`&hbar;`

    F[k](X) = I*m*%Commutator((1/2)*m*Physics:-`^`(v[n], 2), v[k])/`&hbar;`


    Simplify this expression for the quantized force taking the quantum algebra rules (16) into account

    Simplify(F[k](X) = I*m*%Commutator((1/2)*m*Physics[`^`](v[n], 2), v[k])/`&hbar;`)

    F[k](X) = (1/2)*q*(-Physics:-`*`(Physics:-d_[n](A[k](X), [X]), v[n])+Physics:-`*`(Physics:-d_[k](A[n](X), [X]), v[n])-Physics:-`*`(v[n], Physics:-d_[n](A[k](X), [X]))+Physics:-`*`(v[n], Physics:-d_[k](A[n](X), [X])))


    It is not difficult to verify that this is the antisymmetrized vector product `&x`(`#mover(mi("v",mathcolor = "olive"),mo("&rarr;"))`, `#mover(mi("B",mathcolor = "olive"),mo("&rarr;"))`). Departing from `#mover(mi("B",mathcolor = "olive"),mo("&rarr;"))` = `&x`(VectorCalculus[Nabla], `#mover(mi("A",mathcolor = "olive"),mo("&rarr;"))`) expressed using tensor notation,

    B[c](X) = LeviCivita[c, n, m]*d_[n](A[m](X))

    B[c](X) = -Physics:-LeviCivita[c, m, n]*Physics:-d_[n](A[m](X), [X])


    and taking into acount that

     Component(`&x`(`#mover(mi("v",mathcolor = "olive"),mo("&rarr;"))`, `#mover(mi("B",mathcolor = "olive"),mo("&rarr;"))`), k) = `&epsilon;`[b, c, k]*v[b]*B[c](X) 

    multiply both sides of (19) by `&epsilon;`[b, c, k]*v[b], getting

    LeviCivita[k, b, c]*v[b]*(B[c](X) = -Physics[LeviCivita][c, m, n]*Physics[d_][n](A[m](X), [X]))

    Physics:-LeviCivita[b, c, k]*Physics:-`*`(v[b], B[c](X)) = -Physics:-LeviCivita[b, c, k]*Physics:-LeviCivita[c, m, n]*Physics:-`*`(v[b], Physics:-d_[n](A[m](X), [X]))


    Simplify(Physics[LeviCivita][b, c, k]*Physics[`*`](v[b], B[c](X)) = -Physics[LeviCivita][b, c, k]*Physics[LeviCivita][c, m, n]*Physics[`*`](v[b], Physics[d_][n](A[m](X), [X])))

    Physics:-LeviCivita[b, c, k]*Physics:-`*`(v[b], B[c](X)) = Physics:-`*`(v[m], Physics:-d_[k](A[m](X), [X]))-Physics:-`*`(v[n], Physics:-d_[n](A[k](X), [X]))


    Finally, replacing the repeated index m by n 

    SubstituteTensorIndices(m = n, Physics[LeviCivita][b, c, k]*Physics[`*`](v[b], B[c](X)) = Physics[`*`](v[m], Physics[d_][k](A[m](X), [X]))-Physics[`*`](v[n], Physics[d_][n](A[k](X), [X])))

    Physics:-LeviCivita[b, c, k]*Physics:-`*`(v[b], B[c](X)) = Physics:-`*`(v[n], Physics:-d_[k](A[n](X), [X]))-Physics:-`*`(v[n], Physics:-d_[n](A[k](X), [X]))


    Likewise, for

     Component(`&x`(`#mover(mi("v",mathcolor = "olive"),mo("&rarr;"))`, `#mover(mi("B",mathcolor = "olive"),mo("&rarr;"))`), k) = `&epsilon;`[b, c, k]*B[b]*B[c](X) 

    multiplying (19), this time from the right instead of from the left, we get

    Simplify(((B[c](X) = -Physics[LeviCivita][c, m, n]*Physics[d_][n](A[m](X), [X]))*LeviCivita[k, b, c])*v[b])

    Physics:-LeviCivita[b, c, k]*Physics:-`*`(B[c](X), v[b]) = Physics:-`*`(Physics:-d_[k](A[m](X), [X]), v[m])-Physics:-`*`(Physics:-d_[n](A[k](X), [X]), v[n])


    SubstituteTensorIndices(m = n, Physics[LeviCivita][b, c, k]*Physics[`*`](B[c](X), v[b]) = Physics[`*`](Physics[d_][k](A[m](X), [X]), v[m])-Physics[`*`](Physics[d_][n](A[k](X), [X]), v[n]))

    Physics:-LeviCivita[b, c, k]*Physics:-`*`(B[c](X), v[b]) = Physics:-`*`(Physics:-d_[k](A[n](X), [X]), v[n])-Physics:-`*`(Physics:-d_[n](A[k](X), [X]), v[n])


    Simplifying now the expression (18) for the quantized force taking into account (22) and (24) we get

    simplify(F[k](X) = (1/2)*q*(-Physics[`*`](Physics[d_][n](A[k](X), [X]), v[n])+Physics[`*`](Physics[d_][k](A[n](X), [X]), v[n])-Physics[`*`](v[n], Physics[d_][n](A[k](X), [X]))+Physics[`*`](v[n], Physics[d_][k](A[n](X), [X]))), {(rhs = lhs)(Physics[LeviCivita][b, c, k]*Physics[`*`](v[b], B[c](X)) = Physics[`*`](v[n], Physics[d_][k](A[n](X), [X]))-Physics[`*`](v[n], Physics[d_][n](A[k](X), [X]))), (rhs = lhs)(Physics[LeviCivita][b, c, k]*Physics[`*`](B[c](X), v[b]) = Physics[`*`](Physics[d_][k](A[n](X), [X]), v[n])-Physics[`*`](Physics[d_][n](A[k](X), [X]), v[n]))})

    F[k](X) = (1/2)*q*Physics:-LeviCivita[b, c, k]*(Physics:-`*`(v[b], B[c](X))+Physics:-`*`(B[c](X), v[b]))



    `#mover(mi("F",mathcolor = "olive"),mo("&rarr;"))` = (1/2)*q*(`&x`(`#mover(mi("v",mathcolor = "olive"),mo("&rarr;"))`, `#mover(mi("B",mathcolor = "olive"),mo("&rarr;"))`)-`&x`(`#mover(mi("B",mathcolor = "olive"),mo("&rarr;"))`, `#mover(mi("v",mathcolor = "olive"),mo("&rarr;"))`))

    in tensor notation. Finally, we note that this operator is Hermitian as expected

    (F[k](X) = (1/2)*q*Physics[LeviCivita][b, c, k]*(Physics[`*`](v[b], B[c](X))+Physics[`*`](B[c](X), v[b])))-Dagger(F[k](X) = (1/2)*q*Physics[LeviCivita][b, c, k]*(Physics[`*`](v[b], B[c](X))+Physics[`*`](B[c](X), v[b])))

    F[k](X)-Physics:-Dagger(F[k](X)) = 0


    Download:,   Quantization_of_the_Lorentz_force.pdf

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

    After days of fruitlessly searching the help files and the Internet for a means of converting a Dataseries of strings into numeric values in Maple I changed tack and determined how to do it myself.

    I was surprised there was no built in Maple function to do this. From searching the Internet I can see I am not alone.
    Since there are no other solutions on the net, here is mine:

    toNumeric := (y) -> map(x->parse(x),y);

    This creates the function toNumeric() that accepts a DataSeries of text values and returns it converted to a Dataseries of numeric values.

    thing := DataSeries(["10", "20", "30", "55.9"], 'labels' = ["a", "b", "c", "d"]);
    thing := toNumeric(thing); 

    dataframething := DataFrame([["cow", "sheep", "goat", "alpaca"], ["10", "20", "30", "55.9"]], rows = ["a", "b", "c", "d"]);
    dataframething[2] := toNumeric(dataframething[2]);

    If you want to see this in Maple, try:


    Students using Maple often have different needs than non-students. Students need more than just a final answer; they are looking to gain an understanding of the mathematical concepts behind the problems they are asked to solve and to learn how to solve problems. They need an environment that allows them to explore the concepts and break problems down into smaller steps.

    The Student packages in Maple offer focused learning environments in which students can explore and reinforce fundamental concepts for courses in Precalculus, Calculus, Linear Algebra, Statistics, and more. For example, Maple includes step-by-step tutors that allow students to practice integration, differentiation, the finding of limits, and more. The Integration Tutor, shown below, lets a student evaluate an integral by selecting an applicable rule at each step. Maple will also offer hints or show the next step, if asked.  The tutor doesn't only demonstrate how to obtain the result, but is designed for practicing and learning.

    For this blog post, I’d like to focus in on an area of great interest to students: showing step-by-step solutions for a variety of problems in Maple.

    Several commands in the Student packages can show solution steps as either output or inline in an interactive pop-up window. The first few examples return the solution steps as output.

    Precalculus problems:

    The Student:-Basics sub-package provides a collection of commands to help students and teachers explore fundamental mathematical concepts that are core to many disciplines. It features two commands, both of which return step-by-step solutions as output.

    The ExpandSteps command accepts a product of polynomials and displays the steps required to expand the expression:

    ExpandSteps( (a^2-1)/(a/3+1/3) );

    The LinearSolveSteps command accepts an equation in one variable and displays the steps required to solve for that variable.

    LinearSolveSteps( (x+1)/y = 4*y^2 + 3*x, x );

    This command also accepts some nonlinear equations that can be reduced down to linear equations.

    Calculus problems:

    The Student:-Calculus1 sub-package is designed to cover the basic material of a standard first course in single-variable calculus. Several commands in this package provide interactive tutors where you can step through computations and step-by-step solutions can be returned as standard worksheet output.

    Tools like the integration, differentiation, and limit method tutors are interactive interfaces that allow for exploration. For example, similar to the integration-methods tutor above, the differentiation-methods tutor lets a student obtain a derivative by selecting the appropriate rule that applies at each step or by requesting a complete solution all at once. When done, pressing “Close” prints out to the Maple worksheet an annotated solution containing all of the steps.

    For example, try entering the following into Maple:


    Next, right click on the Matrix and choose “Student Calculus1 -> Tutors -> Differentiation Methods…

    The Student:-Calculus1 sub-package is not alone in offering this kind of step-by-step solution finding. Other commands in other Student packages are also capable of returning solutions.

    Linear Algebra Problems:

    The Student:-LinearAlgebra sub-package is designed to cover the basic material of a standard first course in linear algebra. This sub-package features similar tutors to those found in the Calculus1 sub-package. Commands such as the Gaussian EliminationGauss-Jordan Elimination, Matrix Inverse, Eigenvalues or Eigenvectors tutors show step-by-step solutions for linear algebra problems in interactive pop-up tutor windows. Of these tutors, a personal favourite has to be the Gauss-Jordan Elimination tutor, which were I still a student, would have saved me a lot of time and effort searching for simple arithmetic errors while row-reducing matrices.

    For example, try entering the following into Maple:


    Next, right click on the Matrix and choose “Student Linear Algebra -> Tutors -> Gauss-Jordan Elimination Tutor

    This tutor makes it possible to step through row-reducing a matrix by using the controls on the right side of the pop-up window. If you are unsure where to go next, the “Next Step” button can be used to move forward one-step. Pressing “All Steps” returns all of the steps required to row reduce this matrix.

    When this tutor is closed, it does not return results to the Maple worksheet, however it is still possible to use the Maple interface to step through performing elementary row operations and to capture the output in the Maple worksheet. By loading the Student:-LinearAlgebra package, you can simply use the right-click context menu to apply elementary row operations to a Matrix in order to step through the operations, capturing all of your steps along the way!

    An interactive application for showing steps for some problems:

    While working on this blog post, it struck me that we did not have any online interactive applications that could show solution steps, so using the commands that I’ve discussed above, I authored an application that can expand, solve linear problems, integrate, differentiate, or find limits. You can interact with this application here, but note that this application is a work in progress, so feel free to email me (maplepm (at) any strange bugs that you may encounter with it.

    More detail on each of these commands can be found in Maple’s help pages.


    I don't understand why the spam issue has not been resolved yet.  A ton of spam this morning!

    Here's some more help. 

    Any title or question of a post containing http:// should be automatically deleted. 

    Clearly there is a name generator in place somewhere creating new accounts.  Is there not some kind of higher security name creation account done here at mapleprimes?  Ironically you would think this being a math type forum it would be hard to spam.

    While looking at the tags you will find a sizeable link to pages tagged garcinia ( There are no posts with this tag as the posts with this tag were indeed spam. Please remove all the related tags.

    We assume that the radius of the outer stationary circle is  1. If we set the radius  x  of the inner stationary circle, all the other circles are uniquely determined by solving the system Sys.  Should be  x<=1/3 . If  x=1/3  then all the inner circles have a radius  1/3 . The following picture explains the meaning of symbols in the procedure Circles:





    local OO, O1, O2, O3, O4, O2x, O2y, O3x, O3y, OT, T1, T2, T3, s, t, dist, Sys, Sol, sol, y, u, v, z, C0, R0, P;

    uses plottools, plots;

    OO:=[0,0]: O1:=[x+y,0]: O2:=[O2x,O2y]: O3:=[O3x,O3y]: O4:=[-x-z,0]: OT:=[x+2*y-1,0]:

    T1:=(O2*~y+O1*~u)/~(y+u): T2:=(O3*~u+O2*~v)/~(u+v): T3:=(O4*~v+O3*~z)/~(v+z):

    solve({(T2-T1)[1]*(s-((T1+T2)/2)[1])+(T2-T1)[2]*(t-((T1+T2)/2)[2])=0, (T3-T2)[1]*(s-((T2+T3)/2)[1])+(T3-T2)[2]*(t-((T3+T2)/2)[2])=0}, {s,t}):



    Sys:={dist(O1,O2)^2=(y+u)^2, dist(OO,O2)^2=(x+u)^2, dist(O2,O3)^2=(u+v)^2, dist(OO,O3)^2=(x+v)^2, dist(O3,O4)^2=(z+v)^2, x+y+z=1, dist(O2,OT)^2=(1-u)^2, dist(O3,OT)^2=(1-v)^2};


    sol:=select(i->is(eval(convert([y>0,u>0,v>0,z>0,O2y>0,x<=y,u<=y,v<=u,z<=v],`and`),i)), Sol)[];


    O1:=[x+y,0]: O2:=[O2x,O2y]: O3:=[O3x,O3y]: O4:=[-x-z,0]: OT:=[x+2*y-1,0]:




    local eq, r1, r, R, Ot, El, i, S, s, t, P1, P2;

    uses plots,plottools;







    for i from 2 to 6 do

    S:=[solve({1-dist(OT,[s,t])=dist(Ot[i-1],[s,t])-R[i-1], 1-dist(OT,[s,t])=dist(OO,[s,t])-x})];

    P1:=eval([s,t],S[1]); P2:=eval([s,t],S[2]);




    display(El,seq(disk(Ot[k],0.012),k=1..6),circle(C0,R0,color=gold,thickness=3),circle([x+2*y-1,0],1, color=blue,thickness=4), circle(OO,x, color=red,thickness=4), seq(circle(Ot[k],R[k], thickness=3),k=1..6), scaling=constrained, axes=none);

    end proc:

    animate(P,[phi], phi=0..Pi, frames=120);

    end proc:  


    Example of use (I got  x=0.22  just by measuring the ruler displayed original animation):





    The curve on the following animation is an astroid (a special case of hypocycloid). See wiki for details. Hypocycloid procedure creates animation for any hypocycloid.  Parameters of the procedure: R is the radius of the outer circle, r is the radius of the inner circle.


    local A, B, f, g, F;

    uses plots,plottools;






    animate(F,[t], t=0..2*Pi*denom(R/r), frames=90);

    end proc:


    Examples of use:










    The Saturday edition of our local newspaper (Waterloo Region Record) carries, as part of The PUZZLE Corner column, a weekly puzzle "STICKELERS" by Terry Stickels. Back on December 13, 2014, the puzzle was:

    What number comes next?

    1   4   18   96   600   4320   ?

    The solution given was the number 35280, obtained by setting k = 1 in the general term k⋅k!.

    On September 5, 2015, the column contained the puzzle:

    What number comes next?

    2  3  3  5  10  13  39  43  172  177  ?

    The proposed solution was the number 885, obtained as a10 from the recursion

    where a0 =2.

    As a youngster, one of my uncles delighted in teasing me with a similar question for the sequence 36, 9, 50, 55, 62, 71, 79, 18, 20. Ignoring the fact that there is a missing entry between 9 and 50, the next member of the sequence is "Bay Parkway," which is what 22nd Avenue is actually called in the Brooklyn neighborhood of my youth. These are subway stops on what was then called the West End line of the subway that went out to Stillwell Avenue in Coney Island.

    Armed with the skepticism inspired by this provincial chestnut, I looked at both of these puzzles with the attitude that the "next number" could be anything I chose it to be. After all, a sequence is a mapping from the (nonnegative) integers to the reals, and unless the mapping is completely specified, the function values are not well defined.

    Indeed, for the first puzzle, the polynomial f(x) interpolating the points

    (0, 1), (1, 4), (2, 18), (3, 93), (4, 600), (5, 4320)

    reproduces the first six members of the given sequence, and gives 18593 (not 35280) for f(7). In other words, the pattern k⋅k! is not a unique representation of the sequence, given just the first six members. The worksheet NextNumber derives the interpolating polynomial f and establishes that f(n) is an integer for every nonzero integer n.

    Likewise, for the second puzzle, the polynomial g(x) interpolating the points

    (1, 2) ,(2, 3) ,(3, 3) ,(4, 5) ,(5, 10) ,(6, 13) ,(7, 39) ,(8, 43), (9, 172) ,(10, 177)

    reproduces the first ten members of the given sequence, and gives -7331(not 885) for g(11). Once again, the pattern proposed as the "solution" is not unique, given that the worksheet NextNumber contains both g(x) and a proof that for integer n, all values of g(n) are integers.

    The upshot of these observations is that without some guarantee of uniqueness, questions like "what is the next number" are meaningless. It would be far better to pose such challenges with the words "Find a pattern for the given members of the following sequence" and warn that the function capturing that pattern might not be unique.

    I leave it to the interested reader to prove or disprove the following conjecture: Interpolate the first n terms of either sequence. The interpolating polynomial p will reproduce these n terms, but for k>n, p(k) will differ from the corresponding member of the sequence determined by the stated patterns. (Results of limited numerical experiments are consistent with the truth of this conjecture.)


    Bruce Jenkins is President of Ora Research, an engineering research and advisory service. Maplesoft commissioned him to examine how systems-driven engineering practices are being integrated into the early stages of product development, the results of which are available in a free whitepaper entitled System-Level Physical Modeling and Simulation. In this series of blog posts, Mr. Jenkins discusses the results of his research.

    This is the third entry in the series.

    My last post, System-level physical modeling and simulation: Adoption drivers vs. adoption constraints, described my firm’s research project to investigate the contemporary state of adoption and application of systems modeling software technologies, and their attendant methods and work processes, in the engineering design of off-highway equipment and mining machinery.

    In this project, I interviewed some half-dozen expert practitioners at leading manufacturers, including both engineering management and senior discipline leads, to identify key technological factors as well as business and competitive issues driving adoption and use of systems modeling at current levels.

    After identifying present-day adoption drivers as well as current constraints on adoption, finally I sought to learn practitioners’ visions, strategies and best practices for accelerating and institutionalizing the implementation and usage of systems modeling tools and practices in their organizations.

    I was strongly encouraged to find a wealth of avenues and opportunities for exploiting enterprise business drivers, current industry disruptions, and related internal realignments and change-management initiatives to help drive introduction—or proliferation—of these technologies and their associated new ways of working into engineering organizations:

    • Systems modeling essential to compete by creating differentiated products
    • Mechatronics revolution in off-highway equipment
    • Industry downturns and disruptions create opportunities for disruptive innovation
      • Opportunities to leverage change in underlying industry competitive dynamics
      • Mining industry down-cycle creates opportunity to innovate, find new ways of working
      • Some manufacturers are using current down-cycle in mining industry to change their product innovation strategy
    • Strategies of manufacturers pursuing disruptive innovation
      • Best odds are in companies with deep culture of continually inculcating new skills into their people, and rethinking methods and work processes
      • Some managements willing to take radical corporate measures to replace old-thinking engineering staff with “systems thinkers”
      • Downsizing in off-highway equipment manufacturers may push them to seek more systems-level value-add from their component suppliers
    • New technology opportunities inside manufacturers ready to move more deeply into systems modeling
      • Opportunities in new/emerging industries/companies without legacy investments in systems modeling tools and libraries
      • Best practice for introducing systems modeling: start with work process, then bring in software
      • Capitalizing on engineering’s leeway and autonomy in specifying systems modeling software compared with enterprise-standard CAD/PLM tools
    • Systems modeling technology advances anticipated by practitioner advocates
      • Improving software integration, interoperability, data interchange
      • Improving co-simulation across domain tools
      • Better, more complete FMI (Functional Mock-up Interface) implementation/compliance
      • Higher-fidelity versions of FMI or similar

    The white paper detailing the findings of this research is intended to offer guidance and advice for implementing change, as well as documentation to help convince colleagues, management and partners that new ways of working exist, and that the software technologies to support and enable them are available, accessible, and delivering payback and business advantage to forward-thinking engineering organizations today.

    My hope is that this research finds utility as a practical, actionable aid for engineers and engineering management in helping their organizations to adopt and implement—or to strengthen and deepen—a simulation-led, systems-driven approach to product development.

    You can download the full white paper reporting our findings here.

    Bruce Jenkins, Ora Research

    most effective built in operator code award goes to ppl that wrote the code for the union and intercect set operations for maple. Very important simple example below of  one of its applications.


    When i work with algorithms, probably one of my most primary ports of enquiry (figuratively jeez skynet)  is to set up and if statement triggered to terminate the loop once the operations performed for any further cycles is INDEMPOTENT. this doesnt always mean your output is convergent in every case but it allows you to minimize the amount of time the cpu needs to collect data( ie the point at which it would produce that same set as it did in the last most loop)



    Y := proc (X) local N, S1, `&Sopf;`; if X <> `union`(X, S1[N]) then N := (rand(1 .. NrANGE))(); S1[N] := {K[1](4+N), K[1](5+N), K[1](6+N), K[1](7+N), K[1](8+N), K[1](9+N), K[1](10+N), K[1](11+N), K[1](12+N), K[1](13+N), K[1](14+N), K[1](15+N)}; `&Sopf;` := `union`(X, S1[N]) else  end if end proc

    proc (X) local N, S1, `&Sopf;`; if X <> `union`(X, S1[N]) then N := (rand(1 .. NrANGE))(); S1[N] := {K[1](4+N), K[1](5+N), K[1](6+N), K[1](7+N), K[1](8+N), K[1](9+N), K[1](10+N), K[1](11+N), K[1](12+N), K[1](13+N), K[1](14+N), K[1](15+N)}; `&Sopf;` := `union`(X, S1[N]) else  end if end proc







    A prime producing polynomial.


    Observations on the trinomial n2 + n + 41.


    by Matt C. Anderson


    September 3, 2016


    The story so far


    We assume that n is an integer.  We focus our attention on the polynomial n^2 + n + 41.


    Furthur, we analyze the behavior of the factorization of integers of the form


    h(n) = n2 + n + 41                                          (expression 1)


    where n is a non-negative integer.  It was shown by Legendre, in 1798 that if 0 ≤ n < 40 then h(n) is a prime number.


    Certain patterns become evident when considering points (a,n) where


    h(n) ≡ 0 mod a.                                             (expression 2)


    The collection of all such point produces what we are calling a "graph of discrete divisors" due to certain self-similar features.  From experimental data we find that the integer points in this bifurcation graph lie on a collection of parabolic curves indexed by pairs of relatively prime integers.  The expression for the middle parabolas is –


    p(r,c) = (c*x – r*y)2 – r*(c*x – r*y) – x + 41*r2.           (expression 3)


    The restrictions are that 0<r<c and gcd(r,c) = 1 and all four of r,c,x, and y are integers.


    Each such pair (r,c) yields (again determined experimentally and by observation of calculations) an integer polynomial a*z2 + b*z + c, and the quartic h(a*z2 + b*z + c) then factors non-trivially over the integers into two quadratic expressions.  We call this our "parabola conjecture".  Certain symmetries in the bifurcation graph are due to elementary relationships between pairs of co-prime integers.  For instance if m<n are co-prime integers, then there is an observable relationship between the parabola it determines that that formed from (n-m, n).


    We conjecture that all composite values of h(n) arise by substituting integer values of z into h(a*z2 + b*z + c), where this quartic factors algebraically over Z for a*z2 + b*z + c a quadratic polynomial determined by a pair of relatively prime integers.  We name this our "no stray points conjecture" because all the points in the bifurcation graph appear to lie on a parabola.


    We further conjecture that the minimum x-values for parabolas corresponding to (r, c) with gcd(r, c) = 1 are equal for fixed n.  Further, these minimum x-values line up at 163*c^2/4 where c = 2, 3, 4, ...  The numerical evidence seems to support this.  This is called our "parabolas line up" conjecture.


    The notation gcd(r, c) used above is defined here.  The greatest common devisor of two integers is the smallest whole number that divides both of those integers.


    Theorem 1 - Consider h(n) with n a non negative integer. 

    h(n) never has a factor less than 41.


    We prove Theorem 1 with a modular construction.  We make a residue table with all the prime factors less than 41.  The fundamental theorem of arithmetic states that any integer greater than one is either a prime number, or can be written as a unique product of prime numbers (ignoring the order).  So if h(n) never has a prime factor less than 41, then by extension it never has an integer factor less than 41.


    For example, to determine that h(n) is never divisible by 2, note the first column of the residue table.  If n is even, then h(n) is odd.  Similarly, if n is odd then h(n) is also odd.  In either case, h(n) does not have factorization by 2.


    Also, for divisibility by 3, there are 3 cases to check.  They are n = 0, 1, and 2 mod 3. h(0) mod 3 is 2.  h(1) mod 3 is 1. and h(2) mod 3 is 2.  Due to these three cases, h(n) is never divisible by 3.  This is the second column of the residue table.


    The number 0 is first found in the residue table for the cases h(0) mod 41 and h(40) mod 41.  This means that if n is congruent to 0 mod 41 then h(n) will be divisible by 41.  Similarly, if n is congruent to 40 mod 41 then h(n) is also divisible by 41.

    After the residue table, we observe a bifurcation graph which has points when h(y) mod x is divisible by x.  The points (x,y) can be seen on the bifurcation graph.


    < insert residue table here >


    Thus we have shown that h(n) never has a factor less than 41.


    Theorem 2


    Since h(a) = a^2 + a + 41, we want to show that h(a) = h( -a -1).


    Proof of Theorem 2

    Because h(a) = a*(a+1) + 41,

    Now h(-a -1) = (-a -1)*(-a -1 +1) + 41.

    So h(-a -1) = (-a -1)*(-a) +41,

    And h(-a -1) = h(a).

    Which was what we wanted.

    End of proof of theorem 2.


    Corrolary 1

    Further, if h(b) mod c ≡ = then h(c –b -1) mod c ≡ 0.


    We can observe interesting patterns in the “graph of discrete divisors” on a following page.


    Walking into the big blue Maplesoft office on August 3rd was a bit nerve wracking. I had no idea who anyone was, what to expect, or even what I would be doing. As I sat in the front hall waiting for someone to receive me, I remember thinking, “What have I gotten myself into?”. Despite my worries on that first day, interning at Maplesoft has been a great experience! I never knew that I would be able to learn so much about programming and working in a company in such a short amount of time. Although Maple was a programming language that was foreign to me a couple weeks ago, I feel like I’m relatively well versed in it now. Trying to learn a new language in this short timespan hasn’t been easy, but I think that I picked it up quickly, even if I’ve had my fair share of frustrations.

    Chaos Game example on Rosetta Code

    At Maplesoft, I’ve been contributing to the Rosetta Code project by writing short programs using Maple. The Rosetta Code project is dedicated to creating programming examples for many different tasks in different programming languages. My summer project has been to create solutions using Maple for as many tasks as possible and to post these to Rosetta Code; the goal being to have the list of tasks without Maple implementation shrink with each passing day. It’s nice to feel like I’m leaving a mark in this world, even if it is in such a small corner of the internet.

    Flipping Bits example on Rosetta Code/MapleCloud

    This internship, of course, came with its share of challenges. During my work on the Rosetta Code project, I posted solutions for a total of 38 tasks. Some of them were easy, but some of them took days to complete. On some days, I felt like I was on top of the world. Everything I made turned out great and I knew exactly how to tackle each problem. Other days were slower. I’ve spent ages just staring at a computer monitor trying to figure out just how on earth I was going to make this machine do what I wanted it to do! The 24 Game task was particularly hard, but also very educational. Through this task, I learned about modules, a concept previously unknown to me. I’m fairly sure that the 24 Game also took me the longest, whereas the Increment a numerical string task took me no time at all. Despite it being easy, the Increment a numerical string task wasn’t particularly fun; a bit of a challenge is required for something to be entertaining, after all. My personal favourite was the Fibonacci n-step number sequences task. It was the first really challenging task I encountered, and for after which the feeling of finally completing a task that I spent so long on, of finally overcoming that mountain, was extremely satisfying. Not all challenges end in satisfaction, however. I often found myself accidentally doing something that made the window freeze. I would close the program, then cry a bit on the inside when I realized I just lost the past half an hour’s worth of unsaved work. Nevertheless, I’m glad I got to face all these obstacles because they have made me more resilient and a better programmer.

    The following is the code for the Fibonacci n-step number sequences task

    numSequence := proc(initValues :: Array)
    	local n, i, values;
    n := numelems(initValues);
    values := copy(initValues);
    for i from (n+1) to 15 do
    values(i) := add(values[i-n..i-1]);
    end do;
    return values;
    end proc:
    initValues := Array([1]):
    for i from 2 to 10 do
    initValues(i) := add(initValues):
    printf ("nacci(%d): %a\n", i, convert(numSequence(initValues), list));
    end do:
    printf ("lucas: %a\n", convert(numSequence(Array([2, 1])), list));

    Maple was a great software to program with and a fairly straightforward language to learn. Having previously programmed in Java, I found Maple similar enough that transitioning wasn’t too difficult. In fact, every once in a while when I didn`t know what to do for a task, I would take a look at the Java example in Rosetta Code and it would point me in a direction or give me some hints. While the two languages are similar, there are still many differences. For example, I liked the fact that in Maple, lists started at an index of 1 rather than 0 and arrays could an arbitrary starting index. Although it was different from what I was used to, I found that it made many things much less confusing. Another thing I liked was that the for loop syntax was very simple. I never once had to run through in my head how many times something would loop for. There were such a wide variety of commands in Maple. There was a command for practically anything, and if you knew that it existed and how to use it, then so much power could be at your fingertips. This is where the help system came in extremely handy. With a single search you might find that the solution to the exact problem you were trying to solve already existed as a Maple command. I always had a help window open when I was using Maple.

    Multiplication Tables example on Rosetta Code

    Spending my summer coding at Maplesoft has been fun, sometimes challenging, but an overall rewarding experience. Through contributing to the Rosetta Code project, I’ve learned so much about computer programming, and it certainly made the 45 minute drive out to Waterloo worth it!

    Yili Xu,
    Maplesoft SHAD Intern

    The presentation below is on undergrad Quantum Mechanics. Tackling this topic within a computer algebra worksheet in the way it's done below, however, is an exciting novelty and illustrates well the level of abstraction that is now possible using the Physics package.


    Quantum Mechanics: Schrödinger vs Heisenberg picture

    Pascal Szriftgiser1 and Edgardo S. Cheb-Terrab2 

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

    (2) Maplesoft


    Within the Schrödinger picture of Quantum Mechanics, the time evolution of the state of a system, represented by a Ket "| psi(t) >", is determined by Schrödinger's equation:

    I*`&hbar;`*(diff(Ket(psi, t), t)) = H*Ket(psi, t)

    where H, the Hamiltonian, as well as the quantum operators O__S representing observable quantities, are all time-independent.


    Within the Heisenberg picture, a Ket Ket(psi, 0) representing the state of the system does not evolve with time, but the operators O__H(t)representing observable quantities, and through them the Hamiltonian H, do.


    Problem: Departing from Schrödinger's equation,


    a) Show that the expected value of a physical observable in Schrödinger's and Heisenberg's representations is the same, i.e. that

    Bra(psi, t)*O__S*Ket(psi, t) = Bra(psi, 0)*O__H(t)*Ket(psi, 0)


    b) Show that the evolution equation of an observable O__H in Heisenberg's picture, equivalent to Schrödinger's equation,  is given by:

    diff(O__H(t), t) = (-I*Physics:-Commutator(O__H(t), H))*(1/`&hbar;`)

    where in the right-hand-side we see the commutator of O__H with the Hamiltonian of the system.

    Solution: Let O__S and O__H respectively be operators representing one and the same observable quantity in Schrödinger's and Heisenberg's pictures, and H be the operator representing the Hamiltonian of a physical system. All of these operators are Hermitian. So we start by setting up the framework for this problem accordingly, including that the time t and Planck's constant are real. To automatically combine powers of the same base (happening frequently in what follows) we also set combinepowersofsamebase = true. The following input/output was obtained using the latest Physics update (Aug/31/2016) distributed on the Maplesoft R&D Physics webpage.


    Physics:-Setup(hermitianoperators = {H, O__H, O__S}, realobjects = {`&hbar;`, t}, combinepowersofsamebase = true, mathematicalnotation = true)

    [combinepowersofsamebase = true, hermitianoperators = {H, O__H, O__S}, mathematicalnotation = true, realobjects = {`&hbar;`, t}]


    Let's consider Schrödinger's equation

    I*`&hbar;`*(diff(Ket(psi, t), t)) = H*Ket(psi, t)

    I*`&hbar;`*(diff(Physics:-Ket(psi, t), t)) = Physics:-`*`(H, Physics:-Ket(psi, t))


    Now, H is time-independent, so (2) can be formally solved: psi(t) is obtained from the solution psi(0) at time t = 0, as follows:

    T := exp(-I*H*t/`&hbar;`)



    Ket(psi, t) = T*Ket(psi, 0)

    Physics:-Ket(psi, t) = Physics:-`*`(exp(-I*t*H/`&hbar;`), Physics:-Ket(psi, 0))


    To check that (4) is a solution of (2), substitute it in (2):

    eval(I*`&hbar;`*(diff(Physics[Ket](psi, t), t)) = Physics[`*`](H, Physics[Ket](psi, t)), Physics[Ket](psi, t) = Physics[`*`](exp(-I*H*t/`&hbar;`), Physics[Ket](psi, 0)))

    Physics:-`*`(H, exp(-I*t*H/`&hbar;`), Physics:-Ket(psi, 0)) = Physics:-`*`(H, exp(-I*t*H/`&hbar;`), Physics:-Ket(psi, 0))


    Next, to relate the Schrödinger and Heisenberg representations of an Hermitian operator O representing an observable physical quantity, recall that the value expected for this quantity at time t during a measurement is given by the mean value of the corresponding operator (i.e., bracketing it with the state of the system Ket(psi, t)).

    So let O__S be an observable in the Schrödinger picture: its mean value is obtained by bracketing the operator with equation (4):

    Dagger(Ket(psi, t) = Physics[`*`](exp(-I*H*t/`&hbar;`), Ket(psi, 0)))*O__S*(Ket(psi, t) = Physics[`*`](exp(-I*H*t/`&hbar;`), Ket(psi, 0)))

    Physics:-`*`(Physics:-Bra(psi, t), O__S, Physics:-Ket(psi, t)) = Physics:-`*`(Physics:-Bra(psi, 0), exp(I*t*H/`&hbar;`), O__S, exp(-I*t*H/`&hbar;`), Physics:-Ket(psi, 0))


    The composed operator within the bracket on the right-hand-side is the operator O in Heisenberg's picture, O__H(t)

    Dagger(T)*O__S*T = O__H(t)

    Physics:-`*`(exp(I*t*H/`&hbar;`), O__S, exp(-I*t*H/`&hbar;`)) = O__H(t)


    Analogously, inverting this equation,

    (T*(Physics[`*`](exp(I*H*t/`&hbar;`), O__S, exp(-I*H*t/`&hbar;`)) = O__H(t)))*Dagger(T)

    O__S = Physics:-`*`(exp(-I*t*H/`&hbar;`), O__H(t), exp(I*t*H/`&hbar;`))


    As an aside to the problem, we note from these two equations, and since the operator T = exp((-I*H*t)*(1/`&hbar;`)) is unitary (because H is Hermitian), that the switch between Schrödinger's and Heisenberg's pictures is accomplished through a unitary transformation.


    Inserting now this value of O__S from (8) in the right-hand-side of (6), we get the answer to item a)

    lhs(Physics[`*`](Bra(psi, t), O__S, Ket(psi, t)) = Physics[`*`](Bra(psi, 0), exp(I*H*t/`&hbar;`), O__S, exp(-I*H*t/`&hbar;`), Ket(psi, 0))) = eval(rhs(Physics[`*`](Bra(psi, t), O__S, Ket(psi, t)) = Physics[`*`](Bra(psi, 0), exp(I*H*t/`&hbar;`), O__S, exp(-I*H*t/`&hbar;`), Ket(psi, 0))), O__S = Physics[`*`](exp(-I*H*t/`&hbar;`), O__H(t), exp(I*H*t/`&hbar;`)))

    Physics:-`*`(Physics:-Bra(psi, t), O__S, Physics:-Ket(psi, t)) = Physics:-`*`(Physics:-Bra(psi, 0), O__H(t), Physics:-Ket(psi, 0))


    where, on the left-hand-side, the Ket representing the state of the system is evolving with time (Schrödinger's picture), while on the the right-hand-side the Ket `&psi;__0`is constant and it is O__H(t), the operator representing an observable physical quantity, that evolves with time (Heisenberg picture). As expected, both pictures result in the same expected value for the physical quantity represented by O.


    To complete item b), the derivation of the evolution equation for O__H(t), we take the time derivative of the equation (7):

    diff((rhs = lhs)(Physics[`*`](exp(I*H*t/`&hbar;`), O__S, exp(-I*H*t/`&hbar;`)) = O__H(t)), t)

    diff(O__H(t), t) = I*Physics:-`*`(H, exp(I*t*H/`&hbar;`), O__S, exp(-I*t*H/`&hbar;`))/`&hbar;`-I*Physics:-`*`(exp(I*t*H/`&hbar;`), O__S, H, exp(-I*t*H/`&hbar;`))/`&hbar;`


    To rewrite this equation in terms of the commutator  Physics:-Commutator(O__S, H), it suffices to re-order the product  H  exp(I*H*t/`&hbar;`) placing the exponential first:

    Library:-SortProducts(diff(O__H(t), t) = I*Physics[`*`](H, exp(I*H*t/`&hbar;`), O__S, exp(-I*H*t/`&hbar;`))/`&hbar;`-I*Physics[`*`](exp(I*H*t/`&hbar;`), O__S, H, exp(-I*H*t/`&hbar;`))/`&hbar;`, [exp(I*H*t/`&hbar;`), H], usecommutator)

    diff(O__H(t), t) = I*Physics:-`*`(exp(I*t*H/`&hbar;`), H, O__S, exp(-I*t*H/`&hbar;`))/`&hbar;`-I*Physics:-`*`(exp(I*t*H/`&hbar;`), Physics:-`*`(H, O__S)+Physics:-Commutator(O__S, H), exp(-I*t*H/`&hbar;`))/`&hbar;`


    Normal(diff(O__H(t), t) = I*Physics[`*`](exp(I*H*t/`&hbar;`), H, O__S, exp(-I*H*t/`&hbar;`))/`&hbar;`-I*Physics[`*`](exp(I*H*t/`&hbar;`), Physics[`*`](H, O__S)+Physics[Commutator](O__S, H), exp(-I*H*t/`&hbar;`))/`&hbar;`)

    diff(O__H(t), t) = -I*Physics:-`*`(exp(I*t*H/`&hbar;`), Physics:-Commutator(O__S, H), exp(-I*t*H/`&hbar;`))/`&hbar;`


    Finally, to express the right-hand-side in terms of  Physics:-Commutator(O__H(t), H) instead of Physics:-Commutator(O__S, H), we take the commutator of the equation (8) with the Hamiltonian

    Commutator(O__S = Physics[`*`](exp(-I*H*t/`&hbar;`), O__H(t), exp(I*H*t/`&hbar;`)), H)

    Physics:-Commutator(O__S, H) = Physics:-`*`(exp(-I*t*H/`&hbar;`), Physics:-Commutator(O__H(t), H), exp(I*t*H/`&hbar;`))


    Combining these two expressions, we arrive at the expected result for b), the evolution equation of a given observable O__H in Heisenberg's picture

    eval(diff(O__H(t), t) = -I*Physics[`*`](exp(I*H*t/`&hbar;`), Physics[Commutator](O__S, H), exp(-I*H*t/`&hbar;`))/`&hbar;`, Physics[Commutator](O__S, H) = Physics[`*`](exp(-I*H*t/`&hbar;`), Physics[Commutator](O__H(t), H), exp(I*H*t/`&hbar;`)))

    diff(O__H(t), t) = -I*Physics:-Commutator(O__H(t), H)/`&hbar;`


    Download:     Schrodinger_vs_Heisenberg_picture.pdf

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

    First 8 9 10 11 12 13 14 Last Page 10 of 252