Applications, Examples and Libraries

Share your work here

December 2017: This is, perhaps, one of the most complicated computations done in this area using the Physics package. To the best of my knowledge never before performed on a computer algebra worksheet. It is exciting to present a computation like this one. At the end the corresponding worksheet is linked so that it can be downloaded and the sections be opened, the computation be reproduced. There is also a link to a pdf with everything open.  Special thanks to Pascal Szriftgiser for bringing this problem. To reproduce the computations below, please update the Physics library with the one distributed from the Maplesoft R&D Physics webpage.

June 17, 2020: updated taking advantage of new features of Maple 2020.1 and Physics Updates v.705. Submitted to arxiv.org

January 25, 2021: updated in the arXiv and submitted for publication in Computer Physics Communications

 

 

Quantum Runge-Lenz Vector and the Hydrogen Atom,

the hidden SO(4) symmetry using Computer Algebra

 

Pascal Szriftgiser1 and Edgardo S. Cheb-Terrab2 

(1) University of Lille, CNRS, UMR 8523 - PhLAM - Physique des Lasers, Atomes et Molécules, F-59000 Lille, France

(2) Maplesoft

NULL

 

Abstract

 

Pauli first noticed the hidden SO(4) symmetry for the Hydrogen atom in the early stages of quantum mechanics [1]. Departing from that symmetry, one can recover the spectrum of a spinless hydrogen atom and the degeneracy of its states without explicitly solving Schrödinger's equation [2]. In this paper, we derive that SO(4) symmetry and spectrum using a computer algebra system (CAS). While this problem is well known [3, 4], its solution involves several steps of manipulating expressions with tensorial quantum operators, simplifying them by taking into account a combination of commutator rules and Einstein's sum rule for repeated indices. Therefore, it is an excellent model to test the current status of CAS concerning this kind of quantum-and-tensor-algebra computations. Generally speaking, when capable, CAS can significantly help with manipulations that, like non-commutative tensor calculus subject to algebra rules, are tedious, time-consuming and error-prone. The presentation also shows a pattern of computer algebra operations that can be useful for systematically tackling more complicated symbolic problems of this kind.

 

Introduction

 

The primary purpose of this work is to derive, step-by-step, the SO(4) symmetry of the Hydrogen atom and its spectrum using a computer algebra system (CAS). To the best of our knowledge, such a derivation using symbolic computation has not been shown before. Part of the goal was also to see whether this computation can be performed entering only the main definition formulas, followed by only simplification commands, and without using previous knowledge of the result. The intricacy of this problem is in the symbolic manipulation and simplification of expressions involving noncommutative quantum tensor operators. The simplifications need to take into account commutator rules, symmetries under permutation of indices of tensorial subexpressions, and use Einstein's sum rule for repeated indices.

We performed the derivation using the Maple 2020 system with the Maplesoft Physics Updates v.705. Generally speaking, the default computational domain of CAS doesn't include tensors, noncommutative operators nor related simplifications. On the other hand, the Maple system is distributed with a Physics package that extends that default domain to include those objects and related operations. Physics includes a Simplify command that takes into account custom algebra rules and the sum rule for repeated indices, and uses tensor-simplification algorithms [5] extended to the noncommutative domain.

 

A note about notation: when working with a CAS, besides the expectation of achieving a correct result for a complicated symbolic calculation, readability is also an issue. It is desired that one be able to enter the definition formulas and computational steps to be performed (the input, preceded by a prompt >, displayed in black) in a way that resembles as closely as possible their paper and pencil representation, and that the results (the output, computed by Maple, displayed in blue) use textbook mathematical-physics notation. The Physics package implements such dedicated typesetting. In what follows, within text and in the output, noncommutative objects are displayed using a different color, e.g. H, vectors and tensor indices are displayed the standard way, as in `#mover(mi("L",mathcolor = "olive"),mo("→"))`, and L[q], and commutators are displayed with a minus subscript, e.g. "[H,L[q]][-]". Although the Maple system allows for providing dedicated typesetting also for the input, we preferred to keep visible the Maple input syntax, allowing for comparison with paper and pencil notation. We collected the names of the commands used and a one line description for them in an Appendix at the end. Maple also implements the concept of inert representations of computations, which are activated only when desired. We use this feature in several places. Inert computations are entered by preceding the command with % and are displayed in grey. Finally, as is usual in CAS, every output has an equation label, which we use throughout the presentation to refer to previous intermediate results.

 

In Sec.1, we recall the standard formulation of the problem and present the computational goal, which is the derivation of the formulas representing the SO(4) symmetry and related spectrum.

 

In Sec.2, we set tensorial non-commutative operators representing position and linear and angular momentum, respectively X[a], p[a] and L[a], their commutation rules used as departure point, and the form of the quantum Hamiltonian H. We also derive a few related identities used in the sections that follow.

 

In Sec.3, we derive the conservation of both angular momentum and the Runge-Lenz quantum operator, respectively "[H,L[q]][-]=0" and "[H,Z[k]][-]=0". Taking advantage of the differentialoperators functionality in the Physics package, we perform the derivation exploring two equivalent approaches; first using only a symbolic tensor representation p[j] of the momentum operator, then using an explicit differential operator representation for it in configuration space, p[j] = -i*`ℏ`*`∂`[j].  With the first approach, expressions are simplified only using the departing commutation rules and Einstein's sum rule for repeated indices. Using the second approach, the problem is additionally transformed into one where the differentiation operators are applied explicitly to a test function G(X). Presenting both approaches is of potential interest as it offers two partly independent methods for performing the same computation, which is helpful to provide confidence on in the results when unknown, a relevant issue when using computer algebra.

 

In Sec. 4, we derive %Commutator(L[m], Z[n]) = I*`ℏ`*`ε`[m, n, u]*Z[u] and show that the classical relation between angular momentum and the Runge-Lenz vectors,  "L *"`#mover(mi("Z"),mo("→"))` = 0, due to the orbital momentum being perpendicular to the elliptic plane of motion while the Runge-Lenz vector lies in that plane, still holds in quantum mechanics, where the components of these quantum vector operators do not commute but "L *"`#mover(mi("Z",mathcolor = "olive"),mo("→"))` = "(Z) *"`#mover(mi("L",mathcolor = "olive"),mo("→"))` = 0.

 

In Sec. 5, we derive "[Z[a],Z[b]][-]=-(2 i `ℏ` `ε`[a,b,c] (H L[c]))/`m__e`" using the two alternative approaches described for Sec.3.

In Sec. 6, we derive the well-known formula for the square of the Runge-Lenz vector, Z[k]^2 = 2*H*(`ℏ`^2+L[a]^2)/m__e+kappa^2.

 

Finally, in Sec. 7, we use the SO(4) algebra derived in the previous sections to obtain the spectrum of the Hydrogen atom. Following the literature, this approach is limited to the bound states for which the energy is negative.

 

Some concluding remarks are presented at the end, and input syntax details are summarized in an Appendix.

 

1. The hidden SO(4) symmetry of the Hydrogen atom

 

 

Let's consider the Hydrogen atom and its Hamiltonian

H = LinearAlgebra[Norm](`#mover(mi("p"),mo("→"))`)^2/(2*m__e)-kappa/r,

 

where `#mover(mi("p"),mo("→"))`is the electron momentum, m__e its mass, κ a real positive constant, r = `≡`(LinearAlgebra[Norm](`#mover(mi("r"),mo("→"))`), sqrt(X[a]^2)) the distance of the electron from the proton located at the origin, and X[a] is its tensorial representation with components ["x, y,z]". We assume that the proton's mass is infinite. The electron and nucleus spin are not taken into account. Classically, from the potential -kappa/r, one can derive a central force `#mover(mi("F"),mo("→"))` = -kappa*`#mover(mi("r"),mo("∧"))`/r^2 that drives the electron's motion. Introducing the angular momentum

 

`#mover(mi("L"),mo("→"))` = `&x`(`#mover(mi("r"),mo("→"))`, `#mover(mi("p"),mo("→"))`),

 

one can further define the Runge-Lenz vector `#mover(mi("Z"),mo("→"))`

 

"Z=1/(`m__e`) (L)*(p)+kappa ( r)/r."

 

It is well known that `#mover(mi("Z"),mo("→"))` is a constant of the motion, i.e. diff(`#mover(mi("Z"),mo("→"))`(t), t) = 0. Switching to Quantum Mechanics, this condition reads

 

%Commutator(H, Z_) = 0.

 

where, for hermiticity purpose, the expression of `#mover(mi("Z",mathcolor = "olive"),mo("→"))` must be symmetrized

 

`#mover(mi("Z",mathcolor = "olive"),mo("→"))` = (`&x`(`#mover(mi("L",mathcolor = "olive"),mo("→"))`, `#mover(mi("p",mathcolor = "olive"),mo("→"))`)-`&x`(`#mover(mi("p",mathcolor = "olive"),mo("→"))`, `#mover(mi("L",mathcolor = "olive"),mo("→"))`))/(2*m__e)+kappa*`#mover(mi("r",mathcolor = "olive"),mo("→"))`/r.

 

In what follows, departing from the Hamiltonian H, the basic commutation rules between position`#mover(mi("r",mathcolor = "olive"),mo("→"))`, momentum `#mover(mi("p",mathcolor = "olive"),mo("→"))` and angular momentum `#mover(mi("L",mathcolor = "olive"),mo("→"))` in tensor notation, we derive the following commutation rules between the quantum Hamiltonian, angular momentum and Runge-Lenz vector `#mover(mi("Z",mathcolor = "olive"),mo("→"))`

 

 

"[H,L[n]][-]"

=

0

"[H,Z[n]][-]"

=

0

" [L[m],Z[n]][-]"

=

I*`ℏ`*`ε`[m, n, o]*Z[o]

" [Z[m],Z[n]][-]"

=

-(2*(I*`ℏ`/m__e))*H*`ε`[m, n, o]*L[o]

 

 

Since H commutes with both `#mover(mi("L",mathcolor = "olive"),mo("→"))`NULL and `#mover(mi("Z",mathcolor = "olive"),mo("→"))`, defining

 

"`M__n`=sqrt(-(`m__e`)/(2 H)) `Z__n`,"

these commutation rules can be rewritten as

 

"[L[m],L[n]][-]"

=

I*`ℏ`*`ε`[m, n, o]*L[o]

" [L[m],M[n]][-]"

=

I*`ℏ`*`ε`[m, n, o]*M[o]

"[M[m],M[n]][-]"

=

I*`ℏ`*`ε`[m, n, o]*L[o]

 

 

  

This set constitutes the Lie algebra of the SO(4) group.

  

 

2. Setting the problem, commutation rules and useful identities

   

3. Commutation rules between the Hamiltonian and each of the angular momentum and Runge-Lenz tensors

   

4. Commutation rules between the angular momentum L[q]and the Runge-Lenz Z[k]tensors

   

5.  Commutation rules between the components of the Runge-Lenz tensor

   

6. The square of the norm of the Runge-Lenz vector

   

7. The atomic hydrogen spectrum

   

Conclusions

 

 

In this presentation, we derived, step-by-step, the SO(4) symmetry of the Hydrogen atom and its spectrum using the symbolic computer algebra Maple system. The derivation was performed without departing from the results, entering only the main definition formulas in eqs. (1), (2) and (5), followed by using a few simplification commands - mainly Simplify, SortProducts and SubstituteTensor - and a handful of Maple basic commands, subs, lhs, rhs and isolate. The computational path that was used to get the results of sections 2 to 7 is not unique. Instead of searching for the shortest path, we prioritized clarity and illustration of the techniques that can be used to crack problems like this one.

This problem is mainly about simplifying expressions using two different techniques. First, expressions with noncommutative operands in products need reduction with respect to the commutator algebra rules that have been set. Second, products of tensorial operators require simplification using the sum rule for repeated indices and the symmetries of tensorial subexpressions. Those techniques, which are part of the Maple Physics simplifier, together with the SortProducts and SubstituteTensor commands for sorting the operands in products to apply tensorial identities, sufficed. The derivations were performed in a reasonably small number of steps.

Two different computational strategies - with and without differential operators - were used in sections 3 and 5, showing an approach for verifying results, a relevant issue in general when performing complicated algebraic manipulations. The Maple Physics ability to handle differential operators as noncommutative operands in products (as frequently done in paper and pencil computations) facilitates readability and ease in entering the computations. The complexity of those operations is then handled by one Physics:-Library command, ApplyProductsOfDifferentialOperators (see eqs. (47) and (83)).

Besides the Maple Physics ability to handle noncommutative tensor operators and simplify such operators using commutator algebra rules, it is interesting to note: a) the ability of the system to factorize expressions involving products of noncommutative operands (see eqs. (90) and (108)) and b) the extension of the algorithms for simplifying tensorial expressions [5] to the noncommutativity domain, used throughout this presentation.

It is also worth mentioning how equation labels can reduce the whole computation to entering the main definitions, followed by applying a few commands to equation labels. That approach helps to reduce the chance of typographical errors to a very strict minimum. Likewise, the fact that commands and equations distribute over each other allows cumbersome manipulations to be performed in simple ways, as done, for instance, in eqs. (8), (9) and (13).

Finally, it was significantly helpful for us to have the typesetting of results using standard mathematical physics notation, as shown in the presentation above.

 

Appendix

 

 

In this presentation, the input lines are preceded by a prompt > and the commands used are of three kinds: some basic Maple manipulation commands, the main Physics package commands to set things and simplify expressions, and two commands of the Physics:-Library to perform specialized, convenient, operations in expressions.

 

The basic Maple commands used

 

• 

interface is used once at the beginning to set the letter used to represent the imaginary unit (default is I but we used i).

• 

isolate is used in several places to isolate a variable in an expression, for example isolating x in a*x+b = 0 results in x = -b/a

• 

lhs and rhs respectively get the left-hand side Aand right-hand side Bof an equation A = B

• 

subs substitutes the left-hand side of an equation by the righ-hand side in a given target, for example subs(A = B, A+C) results in B+C

• 

@ is used to compose commands. So(`@`(A, B))(x) is the same as A(B(x)). This command is useful to express an abstract combo of manipulations, for example as in (108) ≡ lhs = `@`(Factor, rhs).

 

The Physics commands used

 

• 

Setup is used to set algebra rules as well as the dimension of space, type of metric, and conventions as the kind of letter used to represent indices.

• 

Commutator computes the commutator between two objects using the algebra rules set using Setup. If no rules are known to the system, it outputs a representation for the commutator that the system understands.

• 

CompactDisplay is used to avoid redundant display of the functionality of a function.

• 

d_[n] represents the `∂`[n] tensorial differential operator.

• 

Define is used to define tensors, with or without specifying its components.

• 

Dagger  computes the Hermitian transpose of an expression.

• 

Normal, Expand, Factor respectively normalizes, expands and factorizes expressions that involve products of noncommutative operands.

• 

Simplify performs simplification of tensorial expressions involving products of noncommutative factors taking into account Einstein's sum rule for repeated indices, symmetries of the indices of tensorial subexpressions and custom commutator algebra rules.

• 

SortProducts uses the commutation rules set using Setup to sort the non-commutative operands of a product in an indicated ordering.

 

The Physics:-Library commands used

 

• 

Library:-ApplyProductsOfDifferentialOperators applies the differential operators found in a product to the product operands that appear to its right. For example, applying this command to  p*V(X)*m__e results in m__e*p(V(X))

• 

Library:-EqualizeRepeatedIndices equalizes the repeated indices in the terms of a sum, so for instance applying this command to L[a]^2+L[b]^2 results in 2*L[a]^2

 

References

 

[1] W. Pauli, "On the hydrogen spectrum from the standpoint of the new quantum mechanics,” Z. Phys. 36, 336–363 (1926)

[2] S. Weinberg, "Lectures on Quantum Mechanics, second edition, Cambridge University Press," 2015.

[3] Veronika Gáliková, Samuel Kováčik, and Peter Prešnajder, "Laplace-Runge-Lenz vector in quantum mechanics in noncommutative space", J. Math. Phys. 54, 122106 (2013)

[4] Castro, P.G., Kullock, R. "Physics of the so__q(4) hydrogen atom". Theor. Math. Phys. 185, 1678–1684 (2015).

[5] L. R. U. Manssur, R. Portugal, and B. F. Svaiter, "Group-Theoretic Approach for Symbolic Tensor Manipulation," International Journal of Modern Physics C, Vol. 13, No. 07, pp. 859-879 (2002).

 

Download Hidden_SO4_symmetry_of_the_hydrogen_atom.mw

Download Hidden_SO4_symmetry_submitted_to_CPC.pdf (all sections open)


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

hi guys....new image that I am very proud of....utilzing the number theoretic argument λ (n)

 

 

My co-author and I recently published the 3rd edition of our finite element book1 utilizing routines written with MAPLE. In this latest edition, we include a chapter on the meshless method. The meshless method is a unique numerical method for solving PDEs. The finite element method requires the establishment of a mesh associated with node points. Consideration must be given in establishing a good mesh (and minimizing the bandwidth associated with node numbering). The meshless method does not require a mesh to connect nodes. The following excerpt describes the application of the meshless method for a simple 1-D heat transfer simulation using six nodes.

Consider the 1-D expression for heat transfer in a bar defined by the relation2

     (1)

where the 1-D domain is bounded by 0 ≤ xL. The exact solution to this problem is

     (2)

with the exact derivative of the temperature given by

     (3)

 

In order to solve the 1-D problem, a multiquadric (MQ) radial basis function (RBF) is used 

     (4)

where r(x, xj) is the radial (Euclidean) distance from the expansion point (xj) to any point (x) , c is a shape parameter that controls the flatness of the RBF and is set by the user, and n is an integer. With n = 1, we retrieve the inverse multiquadric

     (5)

that will be used to solve Eq. (1). Other types of RBFs are available; the MQ is accurate and popular.

 

A global expansion for the 1-D temperature can be expressed as

     (6)

with the second derivative of the temperature given as

     (7)

Introducing the RBF expansion for the terms in the governing equation, and collocating at the interior points, we obtain

     (8)

At the boundaries, we collocate the RBF expansion to impose the boundary conditions

     (9)

Defining the operator

     (10)

we can now assemble into a fully populated matrix as,

     (11)

 

The solutions obtained using finite difference, finite volume, finite element, boundary element, and the meshless method are listed in Table 1 for 6 equally spaced nodes3 with To = 15 and TL = 25, and L = 1. The interior nodes do not have to be uniformly spaced.

Table 1. Comparison of errors for interior temperatures i = 2,3,…N-1

 

The Maple code listing follows:

> restart:
   with(LinearAlgebra):with(plots):

# MESHLESS METHOD SOLUTION USING MULTIQUADRIC RADIAL BASIS
FUNCTIONS (RBF) il:=6:To:=15:TL:=25:L:=1:
>   x:=[0,1/5,2/5,3/5,4/5,1]:
>   S:=1000:n:=1:dx:=1/(il-1):
>   C:=Array(1..il,1..il):phi:=Array(1..il,1..il):d2phi:=Array(1..il,1..il):
b:=Vector(1..il):TM:=Vector(1..il):alpha:=Vector(1..il):
for i from 1 to il do
   for j from 1 to il do
      phi[i,j]:=(1+(x[i]-x[j])^2/(S*dx^2))^(n-3/2):
      d2phi[i,j]:=3*((x[j]-x[i])/20)^2/(4*((x[j]-x[i])^2/40+1)^(5/2) )-1/(40*((x[j]-x[i])^2/40+1)^(3/2)):
   end do:
end do:
>   for i from 2 to il-1 do    
>       for j from 1 to il do
>         C[i,j]:=d2phi[i,j]+phi[i,j];
            b[i]:=-x[i];
>         C[1,j]:=phi[1,j];
           C[il,j]:=phi[il,j];
         end do:
      end do:
      b[1]:=To:b[il]:=TL:
> #ConditionNumber(C);
>   alpha:=LinearSolve(convert(C,Matrix),b):
TM[1]:=To:TM[6]:=TL:
for i from 2 to il-1 do
   for j from 1 to il do
>         TM[i]:=TM[i]+alpha[j]*(1+(x[i]-x[j])^2/(S*dx^2))^(n-3/2);
 >     end do:
>   end do:
     evalf(TM);

                     

> TE:=To*cos(xx)+(TL+L-To*cos(L))/sin(L)*sin(xx)-xx:
> TE:=subs(TE):
> TE:=plot(TE,xx=0..1,color=blue,legend="Exact",thickness=3):
> MEM:=[seq([subs(x[i]),subs(TM[i])],i=1..6)]:
> T:=plots[pointplot](MEM,style=line,color=red,legend="MEM", thickness=3): 
  MEM:=plots[pointplot](MEM,color=red,legend="MEM",symbol=box, symbolsize=15):

>  plots[display](TE,MEM,T,axes=BOXED,title="Solution - MEM");

              

 

Additional examples for two-dimensional domains are described in the text, along with a chapter on the boundary element method. The meshless method is an interesting numerical approach that belongs to the family of weighted residual techniques. The matrix condition number is on the order of 1010 and can give surprisingly good results – however, the solution can fluctuate when repeatedly executed, eventually returning to the nearly correct solution; this is not an issue when using local assembly instead of the global assemble performed here4. The method can be used to form hybrid schemes, e.g., a finite element method can easily be linked with a meshless method to solve a secondary system of equations for problems involving large domains. Results are not sensitive to the location of the nodes; a random placement of points gives qualitatively similar results as a uniform placement.

Just over a year ago, someone asked me if Maple could help them pick names for their family gift exchange, because they were fed up with trying to find a solution by hand that met all their requirements. I knew it could be done, of course, and I spent some time (and at least one family dinner conversation) talking about how to do it. There wasn’t enough time to help my friend last year, but I dusted off my ideas, and my somewhat rusty Maple programming skills, and put something together for this year.

The problem, as stated to me, was “Assign everyone in the group the name of someone else in the group (the person they will buy a present for), with the restriction that no one can be assigned their partner.”

I decided to generalize a bit so that you can specify more than one person in the “do not pick” list for each individual, and the restrictions do not have to be reciprocal. That way, you can use it with rules like “parents cannot pick their children”, or “Elizabeth got Martin two years running, so she can’t pick him again this year”.

Ultimately I went with a “guess and check” approach. For each person, pick a name from the pool of suitable candidates (excluding themselves, anyone on their “do not pick” list, and anyone who has been picked already). Keep assigning names until either everyone has a name, or you end up in a situation where you can’t give someone a name. This can happen, for instance, if Todd is the last name, and the only unmatched name is Catherine, and Todd cannot pick Catherine. If that happens, I tossed all the names back into the virtual hat, gave it a good shake (i.e. randomize()) and tried again. Not as elegant as I would have liked, but it seemed like an effective approach.

It does feel like there ought to be a “nicer” solution. Maybe using graph theory? I know that my code will get into trouble if the restrictions are such that no solution exists.  If anyone has any ideas on other/better ways to solve this problem I’d be happy to hear them (now that I’ve had the fun of solving it myself first!).  

The application can be found on the application center: Gift Exchange Helper. The name picker algorithm is in the start-up code.

Happy gift giving!

The computation of traces of products of Dirac matrices was implemented years ago - see Physics,Trace .

 

The simplification of products of Dirac matrices, however, was not. Now it is, and illustrating this new feature is the matter of this post. To reproduce the results below please update the Physics library with the one distributed at the Maplesoft R&D Physics webpage.

with(Physics)

 

First of all, when loading Physics, a frequent question is about the signature, the default is (- - - +)

Setup(signature)

[signature = `- - - +`]

(1)

This is important because the convention for the Algebra of Dirac Matrices depends on the signature. With the signatures (- - - +) as well as (+ - - -), the sign of the timelike component is 1

Library:-SignOfTimelikeComponent()

1

(2)

With the signatures (+ + + -) as well as (- + + +), the sign of the timelike component is of course -1

Library:-SignOfTimelikeComponent(`+ + + -`)

-1

(3)

The simplification of products of Dirac Matrices, illustrated below with the default signature, works fine with any of these signatures, and works without having to set a representation for the Dirac matrices -- all the results are representation-independent.

 

The examples below, however, also illustrate a new feature of Physics, for now implemented as a Library:-PerformMatrixOperations command (there is a related, also new, command, Library:-RewriteInMatrixForm, to just present the underlying matrix operations, without performing them). To illustrate this other new functionality , set a representation for the Dirac matrices, say the standard one

 

Setup(Dgamma = standard, math = true)

`* Partial match of  'Physics:-Dgamma' against keyword 'Dgammarepresentation'`

 

`* Partial match of  'math' against keyword 'mathematicalnotation'`

 

`Setting lowercaselatin letters to represent spinor indices `

 

`Defined Dirac gamma matrices (Dgamma) in standard representation`, gamma[1], gamma[2], gamma[3], gamma[4]

 

__________________________________________________

 

[Dgammarepresentation = standard, mathematicalnotation = true]

(4)

The four Dirac matrices are

TensorArray(Dgamma[`~mu`])

Array(%id = 18446744078360533342)

(5)

The definition of the Dirac matrices is implemented in Maple following the conventions of Landau books ([1] Quantum Electrodynamics, V4), and  does not depend on the signature, ie the form of these matrices is

"Library:-RewriteInMatrixForm(?)"

Array(%id = 18446744078360529726)

(6)

With the default signature, the space part components of  gamma[mu] change sign when compared with corresponding ones from gamma[`~mu`] while the timelike component remains unchanged

TensorArray(Dgamma[mu])

Array(%id = 18446744078565663678)

(7)

"Library:-RewriteInMatrixForm(?)"

Array(%id = 18446744078677131982)

(8)

For the default signature, the algebra of the Dirac Matrices, loaded by default when Physics is loaded, is (see page 80 of [1])

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

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

(9)

When the sign of the timelike component of the signature is -1, we have a -1 factor on the right-hand side of (9).

 

Note as well that in (9) the right-hand side has no matrix elements. This is standard in particle physics where the computations are performed algebraically, without performing the matrix operations. For the purpose of actually performing the underlying matrix operations, however, one may want to rewrite this algebra including a 4x4 identity matrix. For that purpose, see Algebra of Dirac Matrices with an identity matrix on the right-hand side. For the purpose of this illustration, below we proceed with the algebra as shown in (9), interpreting right-hand sides as if they involve an identity matrix.

 

Verify the algebra rule by performing all the involved matrix operations

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

Physics:-`*`(Physics:-Dgamma[`~mu`], Physics:-Dgamma[`~nu`])+Physics:-`*`(Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~mu`]) = 2*Physics:-g_[`~mu`, `~nu`]

(10)

Note that, regarding the spacetime indices, this is a 4x4 matrix, whose elements are in turn 4x4 matrices. Compute first the external 4x4 matrix related to mu and nu

TensorArray(Physics[`*`](Physics[Dgamma][`~mu`], Physics[Dgamma][`~nu`])+Physics[`*`](Physics[Dgamma][`~nu`], Physics[Dgamma][`~mu`]) = 2*Physics[g_][`~mu`, `~nu`])

Matrix(%id = 18446744078587020822)

(11)

Perform now all the matrix operations involved in each of the elements of this 4x4 matrix

"Library:-PerformMatrixOperations(?)"

Matrix(%id = 18446744078743243942)

(12)

By eye everything checks OK.NULL

 

Consider now the following five products of Dirac matrices

e0 := Dgamma[mu]^2

Physics:-`*`(Physics:-Dgamma[mu], Physics:-Dgamma[`~mu`])

(13)

e1 := Dgamma[mu]*Dgamma[`~nu`]*Dgamma[mu]

Physics:-`*`(Physics:-Dgamma[mu], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~mu`])

(14)

e2 := Dgamma[mu]*Dgamma[`~lambda`]*Dgamma[`~nu`]*Dgamma[mu]

Physics:-`*`(Physics:-Dgamma[mu], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~mu`])

(15)

e3 := Dgamma[mu]*Dgamma[`~lambda`]*Dgamma[`~nu`]*Dgamma[`~rho`]*Dgamma[mu]

Physics:-`*`(Physics:-Dgamma[mu], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~mu`])

(16)

e4 := Dgamma[mu]*Dgamma[`~lambda`]*Dgamma[`~nu`]*Dgamma[`~rho`]*Dgamma[`~sigma`]*Dgamma[mu]

Physics:-`*`(Physics:-Dgamma[mu], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~sigma`], Physics:-Dgamma[`~mu`])

(17)

New: the simplification of these products is now implemented

e0 = Simplify(e0)

Physics:-`*`(Physics:-Dgamma[mu], Physics:-Dgamma[`~mu`]) = 4

(18)

Verify this result performing the underlying matrix operations

T := SumOverRepeatedIndices(Physics[`*`](Physics[Dgamma][mu], Physics[Dgamma][`~mu`]) = 4)

Physics:-`*`(Physics:-Dgamma[1], Physics:-Dgamma[`~1`])+Physics:-`*`(Physics:-Dgamma[2], Physics:-Dgamma[`~2`])+Physics:-`*`(Physics:-Dgamma[3], Physics:-Dgamma[`~3`])+Physics:-`*`(Physics:-Dgamma[4], Physics:-Dgamma[`~4`]) = 4

(19)

Library:-PerformMatrixOperations(T)

Matrix(%id = 18446744078553169662) = 4

(20)

The same with the other expressions

e1 = Simplify(e1)

Physics:-`*`(Physics:-Dgamma[mu], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~mu`]) = -2*Physics:-Dgamma[`~nu`]

(21)

SumOverRepeatedIndices(Physics[`*`](Physics[Dgamma][mu], Physics[Dgamma][`~nu`], Physics[Dgamma][`~mu`]) = -2*Physics[Dgamma][`~nu`])

Physics:-`*`(Physics:-Dgamma[1], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~1`])+Physics:-`*`(Physics:-Dgamma[2], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~2`])+Physics:-`*`(Physics:-Dgamma[3], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~3`])+Physics:-`*`(Physics:-Dgamma[4], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~4`]) = -2*Physics:-Dgamma[`~nu`]

(22)

T := TensorArray(Physics[`*`](Physics[Dgamma][1], Physics[Dgamma][`~nu`], Physics[Dgamma][`~1`])+Physics[`*`](Physics[Dgamma][2], Physics[Dgamma][`~nu`], Physics[Dgamma][`~2`])+Physics[`*`](Physics[Dgamma][3], Physics[Dgamma][`~nu`], Physics[Dgamma][`~3`])+Physics[`*`](Physics[Dgamma][4], Physics[Dgamma][`~nu`], Physics[Dgamma][`~4`]) = -2*Physics[Dgamma][`~nu`])

Array(%id = 18446744078695012102)

(23)

Library:-PerformMatrixOperations(T)

Array(%id = 18446744078701714238)

(24)

For e2

e2 = Simplify(e2)

Physics:-`*`(Physics:-Dgamma[mu], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~mu`]) = 4*Physics:-g_[`~lambda`, `~nu`]

(25)

SumOverRepeatedIndices(Physics[`*`](Physics[Dgamma][mu], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~mu`]) = 4*Physics[g_][`~lambda`, `~nu`])

Physics:-`*`(Physics:-Dgamma[1], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~1`])+Physics:-`*`(Physics:-Dgamma[2], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~2`])+Physics:-`*`(Physics:-Dgamma[3], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~3`])+Physics:-`*`(Physics:-Dgamma[4], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~4`]) = 4*Physics:-g_[`~lambda`, `~nu`]

(26)

T := TensorArray(Physics[`*`](Physics[Dgamma][1], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~1`])+Physics[`*`](Physics[Dgamma][2], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~2`])+Physics[`*`](Physics[Dgamma][3], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~3`])+Physics[`*`](Physics[Dgamma][4], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~4`]) = 4*Physics[g_][`~lambda`, `~nu`])

Matrix(%id = 18446744078470204942)

(27)

Library:-PerformMatrixOperations(T)

Matrix(%id = 18446744078550068870)

(28)

For e3 we have

e3 = Simplify(e3)

Physics:-`*`(Physics:-Dgamma[mu], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~mu`]) = -2*Physics:-`*`(Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~lambda`])

(29)

Verify this result,

SumOverRepeatedIndices(Physics[`*`](Physics[Dgamma][mu], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~rho`], Physics[Dgamma][`~mu`]) = -2*Physics[`*`](Physics[Dgamma][`~rho`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~lambda`]))

Physics:-`*`(Physics:-Dgamma[1], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~1`])+Physics:-`*`(Physics:-Dgamma[2], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~2`])+Physics:-`*`(Physics:-Dgamma[3], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~3`])+Physics:-`*`(Physics:-Dgamma[4], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~4`]) = -2*Physics:-`*`(Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~lambda`])

(30)

In this case, with three free spacetime indices lambda, nu, rho, the spacetime components form an array 4x4x4 of 64 components, each of which is a matrix equation

T := TensorArray(Physics[`*`](Physics[Dgamma][1], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~rho`], Physics[Dgamma][`~1`])+Physics[`*`](Physics[Dgamma][2], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~rho`], Physics[Dgamma][`~2`])+Physics[`*`](Physics[Dgamma][3], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~rho`], Physics[Dgamma][`~3`])+Physics[`*`](Physics[Dgamma][4], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~rho`], Physics[Dgamma][`~4`]) = -2*Physics[`*`](Physics[Dgamma][`~rho`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~lambda`]))

Array(%id = 18446744078647326830)

(31)

For instance, the first element is

T[1, 1, 1]

Physics:-`*`(Physics:-Dgamma[1], Physics:-`^`(Physics:-Dgamma[`~1`], 4))+Physics:-`*`(Physics:-Dgamma[2], Physics:-`^`(Physics:-Dgamma[`~1`], 3), Physics:-Dgamma[`~2`])+Physics:-`*`(Physics:-Dgamma[3], Physics:-`^`(Physics:-Dgamma[`~1`], 3), Physics:-Dgamma[`~3`])+Physics:-`*`(Physics:-Dgamma[4], Physics:-`^`(Physics:-Dgamma[`~1`], 3), Physics:-Dgamma[`~4`]) = -2*Physics:-`^`(Physics:-Dgamma[`~1`], 3)

(32)

and it checks OK:

Library:-PerformMatrixOperations(T[1, 1, 1])

Matrix(%id = 18446744078647302614) = Matrix(%id = 18446744078647302974)

(33)

How can you test the 64 components of T all at once?

1. Compute the matrices, without displaying the whole thing, take the elements of the array and remove the indices (ie take the right-hand side); call it M

 

M := map(rhs, ArrayElems(Library:-PerformMatrixOperations(T)))

 

For instance,

M[1]

Matrix(%id = 18446744078629635726) = Matrix(%id = 18446744078629636206)

(34)

Now verify all these matrix equations at once: take the elements of the arrays on each side of the equations and verify that the are the same: we expect for output just {true}

 

map(proc (u) options operator, arrow; evalb(map(ArrayElems, u)) end proc, M)

{true}

(35)

The same for e4

e4 = Simplify(e4)

Physics:-`*`(Physics:-Dgamma[mu], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~sigma`], Physics:-Dgamma[`~mu`]) = 2*Physics:-`*`(Physics:-Dgamma[`~sigma`], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`])+2*Physics:-`*`(Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~sigma`])

(36)

SumOverRepeatedIndices(Physics[`*`](Physics[Dgamma][mu], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~rho`], Physics[Dgamma][`~sigma`], Physics[Dgamma][`~mu`]) = 2*Physics[`*`](Physics[Dgamma][`~sigma`], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~rho`])+2*Physics[`*`](Physics[Dgamma][`~rho`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~sigma`]))

Physics:-`*`(Physics:-Dgamma[1], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~sigma`], Physics:-Dgamma[`~1`])+Physics:-`*`(Physics:-Dgamma[2], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~sigma`], Physics:-Dgamma[`~2`])+Physics:-`*`(Physics:-Dgamma[3], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~sigma`], Physics:-Dgamma[`~3`])+Physics:-`*`(Physics:-Dgamma[4], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~sigma`], Physics:-Dgamma[`~4`]) = 2*Physics:-`*`(Physics:-Dgamma[`~sigma`], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`])+2*Physics:-`*`(Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~sigma`])

(37)

Regarding the spacetime indices this is now an array 4x4x4x4

T := TensorArray(Physics[`*`](Physics[Dgamma][1], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~rho`], Physics[Dgamma][`~sigma`], Physics[Dgamma][`~1`])+Physics[`*`](Physics[Dgamma][2], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~rho`], Physics[Dgamma][`~sigma`], Physics[Dgamma][`~2`])+Physics[`*`](Physics[Dgamma][3], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~rho`], Physics[Dgamma][`~sigma`], Physics[Dgamma][`~3`])+Physics[`*`](Physics[Dgamma][4], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~rho`], Physics[Dgamma][`~sigma`], Physics[Dgamma][`~4`]) = 2*Physics[`*`](Physics[Dgamma][`~sigma`], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~rho`])+2*Physics[`*`](Physics[Dgamma][`~rho`], Physics[Dgamma][`~nu`], Physics[Dgamma][`~lambda`], Physics[Dgamma][`~sigma`]))

Array(%id = 18446744078730196382)

(38)

For instance the first of these 256 matrix equations

T[1, 1, 1, 1]

Physics:-`*`(Physics:-Dgamma[1], Physics:-`^`(Physics:-Dgamma[`~1`], 5))+Physics:-`*`(Physics:-Dgamma[2], Physics:-`^`(Physics:-Dgamma[`~1`], 4), Physics:-Dgamma[`~2`])+Physics:-`*`(Physics:-Dgamma[3], Physics:-`^`(Physics:-Dgamma[`~1`], 4), Physics:-Dgamma[`~3`])+Physics:-`*`(Physics:-Dgamma[4], Physics:-`^`(Physics:-Dgamma[`~1`], 4), Physics:-Dgamma[`~4`]) = 4*Physics:-`^`(Physics:-Dgamma[`~1`], 4)

(39)

verifies OK:

Library:-PerformMatrixOperations(Physics[`*`](Physics[Dgamma][1], Physics[`^`](Physics[Dgamma][`~1`], 5))+Physics[`*`](Physics[Dgamma][2], Physics[`^`](Physics[Dgamma][`~1`], 4), Physics[Dgamma][`~2`])+Physics[`*`](Physics[Dgamma][3], Physics[`^`](Physics[Dgamma][`~1`], 4), Physics[Dgamma][`~3`])+Physics[`*`](Physics[Dgamma][4], Physics[`^`](Physics[Dgamma][`~1`], 4), Physics[Dgamma][`~4`]) = 4*Physics[`^`](Physics[Dgamma][`~1`], 4))

Matrix(%id = 18446744078727227382) = Matrix(%id = 18446744078727227862)

(40)

Now all the 256 matrix equations verified at once as done for e3

 

M := map(rhs, ArrayElems(Library:-PerformMatrixOperations(T)))

map(proc (u) options operator, arrow; evalb(map(ArrayElems, u)) end proc, M)

{true}

(41)

Finally, although there is more work to be done here, let's define some tensors and contract their product with these expressions involving products of Dirac matrices.

 

For example,

Define(A, B)

`Defined as tensors`

 

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

(42)

Contract with e1 and e2 and simplify

A[nu]*e1; % = Simplify(%)

A[nu]*Physics:-`*`(Physics:-Dgamma[mu], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~mu`]) = -2*A[`~nu`]*Physics:-Dgamma[nu]

(43)

A[nu]*B[lambda]*e2; % = Simplify(%)

A[nu]*B[lambda]*Physics:-`*`(Physics:-Dgamma[mu], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~mu`]) = 4*B[`~nu`]*A[nu]

(44)

 


 

Download DiracMatricesAndPerformMatrixOperation.mw

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

Are you a broke university or college student, living day to day, sustaining yourself on a package of ramen noodles every day? Scared that you’re missing necessary nutrients? Well this is the Maple application for you!

This app, refurbished and modified to fit the budget and nutritional needs of a student, based on this older app: https://www.mapleprimes.com/maplesoftblog/34983-The-Diet-Problem, takes a list of more than 20 foods, and based on nutritional and budget constraints that you choose and can change, will give you the cheapest option for foods that fits your needs!

Maybe you’re not a broke student, but rather someone who wants to keep track of calorie intake or macronutrient intake on a daily basis, well then this app still works for you!

Find it here:
https://www.maplesoft.com/applications/view.aspx?SID=154368

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

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 28 29 30 31 32 33 34 Last Page 30 of 77