Product Tips & Techniques

Tips and Tricks on how to get the most about Maple and MapleSim

When re-initializing a module (for example, while executing its ModuleApply), you'd likely want to re-initialize all of the remember tables and caches of its procedures. Unfortunately, forget(thismodule) (a direct self reference) doesn't work (I wonder why not?), and it doesn't even tell you that it didn't work. You could do forget(external name of module(an indirect self reference), but I consider that not robust because you'd need to change it if you change the external name. Even less robust is forget~([procedure 1, procedure 2, submodule 3, ...]). Here's something that works and is robust:

forget~([seq(x, x= thismodule)])[]

Or, using Maple 2019 or later syntax,

forget~([for local x in thismodule do x od])[]

Both of these handle submodules and ignore non-procedures. Both handle both locals and exports.

With the new features added to the Student[LinearAlgebra] package I wanted to go over some of the basics on how someone can do Linear Algebra in Maple without require them to do any programming.  I was recently asked about this and thought that the information may be useful to others.
 

This post will be focussed towards new Maple users. I hope that this will be helpful to students using Maple for the first time and professors who want their students to use Maple without needing to spend time learning the language.
 

In addition to the following post you can find a detailed video on using Maple to do Linear Algebra without programming here. You can also find some of the tools that are new to Maple 2020 for Linear Algebra here.

The biggest tools you'll be using are the Matrix palette on the left of Maple, and the Context Panel on the right of Maple.

First you should load the Student[LinearAlgebra] package by entering:

with(Student[LinearAlgebra]);

at the beginning of your document. If you end it with a colon rather than a semi colon it won't display the commands in the package.

Use the Matrix Palette on the left to input Matrices:

 


Once you have a Matrix you can use the context panel on the right to apply a variety of operations to it:

 


The Student Linear Algebra Menu will give you many linear algebra commands.

 


You can also access Maple's Tutors from the Tools Menu

Tools > Tutors > Linear Algebra



If you're interested in using the commands for Student[LinearAlegbra] in Maple you can view the help pages here or by entering:

?Student[LinearAlegbra]

into Maple.

I hope that this helps you get started using Maple for Linear Algebra.



Maple_for_Beginners.mw

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


Vectors in Spherical Coordinates using Tensor Notation

Edgardo S. Cheb-Terrab1 and Pascal Szriftgiser2

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

(1) Maplesoft

 

The following is a topic that appears frequently in formulations: given a 3D vector in spherical (or any curvilinear) coordinates, how do you represent and relate, in simple terms, the vector and the corresponding vectorial operations Gradient, Divergence, Curl and Laplacian using tensor notation?

 

The core of the answer is in the relation between the - say physical - vector components and the more abstract tensor covariant and contravariant components. Focusing the case of a transformation from Cartesian to spherical coordinates, the presentation below starts establishing that relationship between 3D vector and tensor components in Sec.I. In Sec.II, we verify the transformation formulas for covariant and contravariant components on the computer using TransformCoordinates. In Sec.III, those tensor transformation formulas are used to derive the vectorial form of the Gradient in spherical coordinates. In Sec.IV, we switch to using full tensor notation, a curvilinear metric and covariant derivatives to derive the 3D vector analysis traditional formulas in spherical coordinates for the Divergence, Curl, Gradient and Laplacian. On the way, some useful technics, like changing variables in 3D vectorial expressions, differential operators, using Jacobians, and shortcut notations are shown.

 

The computation below is reproducible in Maple 2020 using the Maplesoft Physics Updates v.640 or newer.

 

Start setting the spacetime to be 3-dimensional, Euclidean, and use Cartesian coordinates

with(Physics); with(Vectors)

Setup(dimension = 3, coordinates = cartesian, g_ = `+`, spacetimeindices = lowercaselatin)

`The dimension and signature of the tensor space are set to `[3, `+ + +`]

 

`Default differentiation variables for d_, D_ and dAlembertian are:`*{X = (x, y, z)}

 

`Systems of spacetime coordinates are:`*{X = (x, y, z)}

 

_______________________________________________________

 

`The Euclidean metric in coordinates `*[x, y, z]

 

_______________________________________________________

 

Physics:-g_[mu, nu] = Matrix(%id = 18446744078312229334)

 

(`Defined Pauli sigma matrices (Psigma): `*sigma[1]*`, `*sigma[2]*`, `)*sigma[3]

 

__________________________________________________

 

_______________________________________________________

(1)

I. The line element in spherical coordinates and the scale-factors

 

 

In vector calculus, at the root of everything there is the line element `#mover(mi("dr"),mo("→"))`, which in Cartesian coordinates has the simple form

dr_ = _i*dx+_j*dy+_k*dz

dr_ = _i*dx+_j*dy+_k*dz

(1.1)

To compute the line element  `#mover(mi("dr"),mo("→"))` in spherical coordinates, the starting point is the transformation

tr := `~`[`=`]([X], ChangeCoordinates([X], spherical))

[x = r*sin(theta)*cos(phi), y = r*sin(theta)*sin(phi), z = r*cos(theta)]

(1.2)

Coordinates(S = [r, theta, phi])

`Systems of spacetime coordinates are:`*{S = (r, theta, phi), X = (x, y, z)}

(1.3)

Since in (dr_ = _i*dx+_j*dy+_k*dz)*[dx, dy, dz] are just symbols with no relationship to "[x,y,z],"start transforming these differentials using the chain rule, computing the Jacobian of the transformation (1.2). In this Jacobian J, the first line is "[(∂x)/(∂r)dr", "(∂x)/(∂theta)"`dθ`, "(∂x)/(∂phi)dphi]"

J := VectorCalculus:-Jacobian(map(rhs, [x = r*sin(theta)*cos(phi), y = r*sin(theta)*sin(phi), z = r*cos(theta)]), [S])

 

So in matrix notation,

Vector([dx, dy, dz]) = J.Vector([dr, dtheta, dphi])

Vector[column](%id = 18446744078518652550) = Vector[column](%id = 18446744078518652790)

(1.4)

To complete the computation of  `#mover(mi("dr"),mo("→"))` in spherical coordinates we can now use ChangeBasis , provided that next we substitute (1.4) in the result, expressing the abstract objects [dx, dy, dz] in terms of [dr, `dθ`, `dφ`].

 

In two steps:

lhs(dr_ = _i*dx+_j*dy+_k*dz) = ChangeBasis(rhs(dr_ = _i*dx+_j*dy+_k*dz), spherical)

dr_ = (dx*sin(theta)*cos(phi)+dy*sin(theta)*sin(phi)+dz*cos(theta))*_r+(dx*cos(phi)*cos(theta)+dy*sin(phi)*cos(theta)-dz*sin(theta))*_theta+(cos(phi)*dy-sin(phi)*dx)*_phi

(1.5)

The line element

"simplify(subs(convert(lhs(?) =~ rhs(?),set),dr_ = (dx*sin(theta)*cos(phi)+dy*sin(theta)*sin(phi)+dz*cos(theta))*_r+(dx*cos(phi)*cos(theta)+dy*sin(phi)*cos(theta)-dz*sin(theta))*_theta+(cos(phi)*dy-sin(phi)*dx)*_phi))"

dr_ = _phi*dphi*r*sin(theta)+_theta*dtheta*r+_r*dr

(1.6)

This result is important: it gives us the so-called scale factors, the key that connect 3D vectors with the related covariant and contravariant tensors in curvilinear coordinates. The scale factors are computed from (1.6) by taking the scalar product with each of the unit vectors [`#mover(mi("r"),mo("∧"))`, `#mover(mi("θ",fontstyle = "normal"),mo("∧"))`, `#mover(mi("φ",fontstyle = "normal"),mo("∧"))`], then taking the coefficients of the differentials [dr, `dθ`, `dφ`] (just substitute them by the number 1)

h := subs(`~`[`=`]([dr, `dθ`, `dφ`], 1), [seq(rhs(dr_ = _phi*dphi*r*sin(theta)+_theta*dtheta*r+_r*dr).q, q = [`#mover(mi("r"),mo("∧"))`, `#mover(mi("θ",fontstyle = "normal"),mo("∧"))`, `#mover(mi("φ",fontstyle = "normal"),mo("∧"))`])])

[1, r, r*sin(theta)]

(1.7)

The scale factors are relevant because the components of the 3D vector and the corresponding tensor are not the same in curvilinear coordinates. For instance, representing the differential of the coordinates as the tensor dS^j = [dr, `dθ`, `dφ`], we see that corresponding vector, the line element in spherical coordinates `#mover(mi("dS"),mo("→"))`, is not  constructed by directly equating its components to the components of dS^j = [dr, `dθ`, `dφ`], so  

 

 `#mover(mi("dS"),mo("&rarr;"))` <> `d&phi;`*`#mover(mi("&phi;",fontstyle = "normal"),mo("&and;"))`+dr*`#mover(mi("r"),mo("&and;"))`+`d&theta;`*`#mover(mi("&theta;",fontstyle = "normal"),mo("&and;"))` 

 

The vector `#mover(mi("dS"),mo("&rarr;"))` is constructed multiplying these contravariant components [dr, `d&theta;`, `d&phi;`] by the scaling factors, as

 

 `#mover(mi("dS"),mo("&rarr;"))` = `d&phi;`*`h__&phi;`*`#mover(mi("&phi;",fontstyle = "normal"),mo("&and;"))`+dr*h__r*`#mover(mi("r"),mo("&and;"))`+`d&theta;`*`h__&theta;`*`#mover(mi("&theta;",fontstyle = "normal"),mo("&and;"))` 

 

This rule applies in general. The vectorial components of a 3D vector in an orthogonal system (curvilinear or not) are always expressed in terms of the contravariant components A^j the same way we did in the line above with the line element, using the scale-factors h__j, so that

 

 `#mover(mi("A"),mo("&rarr;"))` = Sum(h[j]*A^j*`#mover(mi("\`e__j\`"),mo("&circ;"))`, j = 1 .. 3)

 

where on the right-hand side we see the contravariant components "A[]^(j)" and the scale-factors h[j]. Because the system is orthogonal, each vector component `#msub(mi("A",fontstyle = "normal"),mfenced(mi("j")))`satisfies

A__j = h[j]*A[`~j`]

 

The scale-factors h[j] do not constitute a tensor, so on the right-hand side we do not sum over j.  Also, from

 

LinearAlgebra[Norm](`#mover(mi("A"),mo("&rarr;"))`) = A[j]*A[`~j`]

it follows that,

A__j = A__j/h__j

where on the right-hand side we now have the covariant tensor components A__j.

 

• 

This relationship between the components of a 3D vector and the contravariant and covariant components of a tensor representing the vector is key to translate vector-component to corresponding tensor-component formulas.

 

II. Transformation of contravariant and covariant tensors

 

 

Define here two representations for one and the same tensor: A__c will represent A in Cartesian coordinates, while A__s will represent A in spherical coordinates.

Define(A__c[j], A__s[j])

`Defined objects with tensor properties`

 

{A__c[j], A__s[j], Physics:-Dgamma[a], Physics:-Psigma[a], Physics:-d_[a], Physics:-g_[a, b], Physics:-LeviCivita[a, b, c], Physics:-SpaceTimeVector[a](S), Physics:-SpaceTimeVector[a](X)}

(2.1)

Transformation rule for a contravariant tensor

 

We know, by definition, that the transformation rule for the components of a contravariant tensor is `#mrow(msup(mi("A"),mi("&mu;",fontstyle = "normal")),mo("&ApplyFunction;"),mfenced(mi("y")),mo("&equals;"),mfrac(mrow(mo("&PartialD;"),msup(mi("y"),mi("&mu;",fontstyle = "normal"))),mrow(mo("&PartialD;"),msup(mi("x"),mi("&nu;",fontstyle = "normal"))),linethickness = "1"),mo("&InvisibleTimes;"),mo("&InvisibleTimes;"),msup(mi("A"),mi("&nu;",fontstyle = "normal")),mfenced(mi("x")))`, that is the same as the rule for the differential of the coordinates. Then, the transformation rule from "`A__c`[]^(j)" to "`A__s`[]^(j)"computed using TransformCoordinates should give the same relation (1.4). The application of the command, however, requires attention, because, as in (1.4), we want the Cartesian (not the spherical) components isolated. That is like performing a reversed transformation. So we will use

 

"TensorArray(`A__c`[]^(j))=TransformCoordinates(tr,`A__s`[]^(j),[X],[S])"

where on the left-hand side we get, isolated, the three components of A in Cartesian coordinates, and on the right-hand side we transform the spherical components "`A__c`[]^(j)", from spherical S = (r, theta, phi) (4th argument) to Cartesian X = (x, y, z) (3rd argument), which according to the 5th bullet of TransformCoordinates  will result in a transformation expressed in terms of the old coordinates (here the spherical [S]). Expand things to make the comparison with (1.4) possible by eye

 

Vector[column](TensorArray(A__c[`~j`])) = TransformCoordinates(tr, A__s[`~j`], [X], [S], simplifier = expand)

Vector[column](%id = 18446744078459463070) = Vector[column](%id = 18446744078459463550)

(2.2)

We see that the transformation rule for a contravariant vector "`A__c`[]^(j)"is, indeed, as the transformation (1.4) for the differential of the coordinates.

Transformation rule for a covariant tensor

 

For the transformation rule for the components of a covariant tensor A__c[j], we know, by definition, that it is `#mrow(msub(mi("A"),mi("&mu;",fontstyle = "normal")),mo("&ApplyFunction;"),mfenced(mi("y")),mo("&equals;"),mfrac(mrow(mo("&PartialD;"),msup(mi("x"),mi("&nu;",fontstyle = "normal"))),mrow(mo("&PartialD;"),msup(mi("y"),mi("&mu;",fontstyle = "normal"))),linethickness = "1"),mo("&InvisibleTimes;"),mo("&InvisibleTimes;"),msub(mi("A"),mi("&nu;",fontstyle = "normal")),mfenced(mi("x")))`, so the same transformation rule for the gradient [`&PartialD;`[x], `&PartialD;`[y], `&PartialD;`[z]], where `&PartialD;`[x] = (proc (u) options operator, arrow; diff(u, x) end proc) and so on. We can experiment this by directly changing variables in the differential operators [`&PartialD;`[x], `&PartialD;`[y], `&PartialD;`[z]], for example

d_[x] = PDEtools:-dchange(tr, proc (u) options operator, arrow; diff(u, x) end proc, simplify)

Physics:-d_[x] = (proc (u) options operator, arrow; ((-r*cos(theta)^2+r)*cos(phi)*(diff(u, r))+sin(theta)*cos(phi)*cos(theta)*(diff(u, theta))-(diff(u, phi))*sin(phi))/(r*sin(theta)) end proc)

(2.3)

This result, and the equivalent ones replacing x by y or z in the input above can be computed in one go, in matricial and simplified form, using the Jacobian of the transformation computed in . We need to take the transpose of the inverse of J (because now we are transforming the components of the gradient   [`&PartialD;`[x], `&PartialD;`[y], `&PartialD;`[z]])

H := simplify(LinearAlgebra:-Transpose(1/J))

Vector([d_[x], d_[y], d_[z]]) = H.Vector([d_[r], d_[theta], d_[phi]])

Vector[column](%id = 18446744078518933014) = Vector[column](%id = 18446744078518933254)

(2.4)

The corresponding transformation equations relating the tensors A__c and A__s in Cartesian and spherical coordinates is computed with TransformCoordinates  as in (2.2), just lowering the indices on the left and right hand sides (i.e., remove the tilde ~)

Vector[column](TensorArray(A__c[j])) = TransformCoordinates(tr, A__s[j], [X], [r, theta, phi], simplifier = expand)

Vector[column](%id = 18446744078557373854) = Vector[column](%id = 18446744078557374334)

(2.5)

We see that the transformation rule for a covariant vector A__c[j] is, indeed, as the transformation rule (2.4) for the gradient.

 

To the side: once it is understood how to compute these transformation rules, we can have the inverse of (2.5) as follows

Vector[column](TensorArray(A__s[j])) = TransformCoordinates(tr, A__c[j], [S], [X], simplifier = expand)

Vector[column](%id = 18446744078557355894) = Vector[column](%id = 18446744078557348198)

(2.6)

III. Deriving the transformation rule for the Gradient using TransformCoordinates

 

 

Turn ON the CompactDisplay  notation for derivatives, so that the differentiation variable is displayed as an index:

ON


The gradient of a function f in Cartesian coordinates and spherical coordinates is respectively given by

(%Nabla = Nabla)(f(X))

%Nabla(f(X)) = (diff(f(X), x))*_i+(diff(f(X), y))*_j+(diff(f(X), z))*_k

(3.1)

(%Nabla = Nabla)(f(S))

%Nabla(f(S)) = (diff(f(S), r))*_r+(diff(f(S), theta))*_theta/r+(diff(f(S), phi))*_phi/(r*sin(theta))

(3.2)

What we want now is to depart from (3.1) in Cartesian coordinates and obtain (3.2) in spherical coordinates using the transformation rule for a covariant tensor computed with TransformCoordinates in (2.5). (An equivalent derivation, simpler and with less steps is done in Sec. IV.)

 

Start changing the vector basis in the gradient (3.1)

lhs(%Nabla(f(X)) = (diff(f(X), x))*_i+(diff(f(X), y))*_j+(diff(f(X), z))*_k) = ChangeBasis(rhs(%Nabla(f(X)) = (diff(f(X), x))*_i+(diff(f(X), y))*_j+(diff(f(X), z))*_k), spherical)

%Nabla(f(X)) = ((diff(f(X), x))*sin(theta)*cos(phi)+(diff(f(X), y))*sin(theta)*sin(phi)+(diff(f(X), z))*cos(theta))*_r+((diff(f(X), x))*cos(phi)*cos(theta)+(diff(f(X), y))*sin(phi)*cos(theta)-(diff(f(X), z))*sin(theta))*_theta+(-(diff(f(X), x))*sin(phi)+cos(phi)*(diff(f(X), y)))*_phi

(3.3)

By eye, we see that in this result the coefficients of [`#mover(mi("r"),mo("&and;"))`, `#mover(mi("&theta;",fontstyle = "normal"),mo("&and;"))`, `#mover(mi("&phi;",fontstyle = "normal"),mo("&and;"))`] are the three lines in the right-hand side of (2.6) after replacing the covariant components A__j by the derivatives of f with respect to the jth coordinate, here displayed using indexed notation due to using CompactDisplay

`~`[`=`]([A__s[1], A__s[2], A__s[3]], [diff(f(S), r), diff(f(S), theta), diff(f(S), phi)])

[A__s[1] = Physics:-Vectors:-diff(f(S), r), A__s[2] = Physics:-Vectors:-diff(f(S), theta), A__s[3] = Physics:-Vectors:-diff(f(S), phi)]

(3.4)

`~`[`=`]([A__c[1], A__c[2], A__c[3]], [diff(f(X), x), diff(f(X), y), diff(f(X), z)])

[A__c[1] = Physics:-Vectors:-diff(f(X), x), A__c[2] = Physics:-Vectors:-diff(f(X), y), A__c[3] = Physics:-Vectors:-diff(f(X), z)]

(3.5)

So since (2.5) is the inverse of (2.6), replace A by ∂ f in (2.5), the formula computed using TransformCoordinates, then insert the result in (3.3) to relate the gradient in Cartesian and spherical coordinates. We expect to arrive at the formula for the gradient in spherical coordinates (3.2) .

"subs([A__s[1] = Physics:-Vectors:-diff(f(S),r), A__s[2] = Physics:-Vectors:-diff(f(S),theta), A__s[3] = Physics:-Vectors:-diff(f(S),phi)],[A__c[1] = Physics:-Vectors:-diff(f(X),x), A__c[2] = Physics:-Vectors:-diff(f(X),y), A__c[3] = Physics:-Vectors:-diff(f(X),z)],?)"

Vector[column](%id = 18446744078344866862) = Vector[column](%id = 18446744078344866742)

(3.6)

"subs(convert(lhs(?) =~ rhs(?),set),%Nabla(f(X)) = (diff(f(X),x)*sin(theta)*cos(phi)+diff(f(X),y)*sin(theta)*sin(phi)+diff(f(X),z)*cos(theta))*_r+(diff(f(X),x)*cos(phi)*cos(theta)+diff(f(X),y)*sin(phi)*cos(theta)-diff(f(X),z)*sin(theta))*_theta+(-diff(f(X),x)*sin(phi)+cos(phi)*diff(f(X),y))*_phi)"

%Nabla(f(X)) = ((sin(theta)*cos(phi)*(diff(f(S), r))+cos(theta)*cos(phi)*(diff(f(S), theta))/r-sin(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(theta)*cos(phi)+(sin(theta)*sin(phi)*(diff(f(S), r))+cos(theta)*sin(phi)*(diff(f(S), theta))/r+cos(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(theta)*sin(phi)+(cos(theta)*(diff(f(S), r))-sin(theta)*(diff(f(S), theta))/r)*cos(theta))*_r+((sin(theta)*cos(phi)*(diff(f(S), r))+cos(theta)*cos(phi)*(diff(f(S), theta))/r-sin(phi)*(diff(f(S), phi))/(r*sin(theta)))*cos(phi)*cos(theta)+(sin(theta)*sin(phi)*(diff(f(S), r))+cos(theta)*sin(phi)*(diff(f(S), theta))/r+cos(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(phi)*cos(theta)-(cos(theta)*(diff(f(S), r))-sin(theta)*(diff(f(S), theta))/r)*sin(theta))*_theta+(-(sin(theta)*cos(phi)*(diff(f(S), r))+cos(theta)*cos(phi)*(diff(f(S), theta))/r-sin(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(phi)+cos(phi)*(sin(theta)*sin(phi)*(diff(f(S), r))+cos(theta)*sin(phi)*(diff(f(S), theta))/r+cos(phi)*(diff(f(S), phi))/(r*sin(theta))))*_phi

(3.7)

Simplifying, we arrive at (3.2)

(lhs = `@`(`@`(expand, simplify), rhs))(%Nabla(f(X)) = ((sin(theta)*cos(phi)*(diff(f(S), r))+cos(theta)*cos(phi)*(diff(f(S), theta))/r-sin(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(theta)*cos(phi)+(sin(theta)*sin(phi)*(diff(f(S), r))+cos(theta)*sin(phi)*(diff(f(S), theta))/r+cos(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(theta)*sin(phi)+(cos(theta)*(diff(f(S), r))-sin(theta)*(diff(f(S), theta))/r)*cos(theta))*_r+((sin(theta)*cos(phi)*(diff(f(S), r))+cos(theta)*cos(phi)*(diff(f(S), theta))/r-sin(phi)*(diff(f(S), phi))/(r*sin(theta)))*cos(phi)*cos(theta)+(sin(theta)*sin(phi)*(diff(f(S), r))+cos(theta)*sin(phi)*(diff(f(S), theta))/r+cos(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(phi)*cos(theta)-(cos(theta)*(diff(f(S), r))-sin(theta)*(diff(f(S), theta))/r)*sin(theta))*_theta+(-(sin(theta)*cos(phi)*(diff(f(S), r))+cos(theta)*cos(phi)*(diff(f(S), theta))/r-sin(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(phi)+cos(phi)*(sin(theta)*sin(phi)*(diff(f(S), r))+cos(theta)*sin(phi)*(diff(f(S), theta))/r+cos(phi)*(diff(f(S), phi))/(r*sin(theta))))*_phi)

%Nabla(f(X)) = (diff(f(S), r))*_r+(diff(f(S), theta))*_theta/r+(diff(f(S), phi))*_phi/(r*sin(theta))

(3.8)

%Nabla(f(S)) = (diff(f(S), r))*_r+(diff(f(S), theta))*_theta/r+(diff(f(S), phi))*_phi/(r*sin(theta))

%Nabla(f(S)) = (diff(f(S), r))*_r+(diff(f(S), theta))*_theta/r+(diff(f(S), phi))*_phi/(r*sin(theta))

(3.9)

IV. Deriving the transformation rule for the Divergence, Curl, Gradient and Laplacian, using TransformCoordinates and Covariant derivatives

 

 

• 

The Divergence

 

Introducing the vector A in spherical coordinates, its Divergence is given by

A__s_ := A__r(S)*_r+`A__&theta;`(S)*_theta+`A__&phi;`(S)*_phi

A__r(S)*_r+`A__&theta;`(S)*_theta+`A__&phi;`(S)*_phi

(4.1)

CompactDisplay(%)

` A__r`(S)*`will now be displayed as`*A__r

 

` A__&phi;`(S)*`will now be displayed as`*`A__&phi;`

 

` A__&theta;`(S)*`will now be displayed as`*`A__&theta;`

(4.2)

%Divergence(%A__s_) = Divergence(A__s_)

%Divergence(%A__s_) = ((diff(A__r(S), r))*r+2*A__r(S))/r+((diff(`A__&theta;`(S), theta))*sin(theta)+`A__&theta;`(S)*cos(theta))/(r*sin(theta))+(diff(`A__&phi;`(S), phi))/(r*sin(theta))

(4.3)

We want to see how this result, (4.3), can be obtained using TransformCoordinates and departing from a tensorial representation of the object, this time the covariant derivative "`&dtri;`[j](`A__s`[]^(j))". For that purpose, we first transform the coordinates and the metric introducing nonzero Christoffel symbols

TransformCoordinates(tr, g_[j, k], [S], setmetric)

`Systems of spacetime coordinates are:`*{S = (r, theta, phi), X = (x, y, z)}

 

`Changing the differentiation variables used to compute the Christoffel symbols from `[x, y, z]*` to `[r, theta, phi]*` while the spacetime metric depends on `[r, theta]

 

`Default differentiation variables for d_, D_ and dAlembertian are:`*{S = (r, theta, phi)}

 

_______________________________________________________

 

`Coordinates: `[r, theta, phi]*`. Signature: `(`+ + -`)

 

_______________________________________________________

 

Physics:-g_[a, b] = Matrix(%id = 18446744078312216446)

 

_______________________________________________________

 

`Setting `*greek*` letters to represent `*space*` indices`

(4.4)

To the side: despite having nonzero Christoffel symbols, the space still has no curvature, all the components of the Riemann tensor are equal to zero

Riemann[nonzero]

Physics:-Riemann[a, b, c, d] = {}

(4.5)

Consider now the divergence of the contravariant "`A__s`[]^(j)"tensor, computed in tensor notation

CompactDisplay(A__s(S))

` A__s`(S)*`will now be displayed as`*A__s

(4.6)

D_[j](A__s[`~j`](S))

Physics:-D_[j](A__s[`~j`](S), [S])

(4.7)

To the side: the covariant derivative  expressed using the D_  operator can be rewritten in terms of the non-covariant d_  and Christoffel  symbols as follows

D_[j](A__s[`~j`](S), [S]) = convert(D_[j](A__s[`~j`](S), [S]), d_)

Physics:-D_[j](A__s[`~j`](S), [S]) = Physics:-d_[j](A__s[`~j`](S), [S])+Physics:-Christoffel[`~j`, a, j]*A__s[`~a`](S)

(4.8)

Summing over the repeated indices in (4.7), we have

%D_[j](%A__s[`~j`]) = SumOverRepeatedIndices(D_[j](A__s[`~j`](S), [S]))

%D_[j](%A__s[`~j`]) = diff(A__s[`~1`](S), r)+diff(A__s[`~2`](S), theta)+diff(A__s[`~3`](S), phi)+2*A__s[`~1`](S)/r+cos(theta)*A__s[`~2`](S)/sin(theta)

(4.9)

How is this related to the expression of the VectorCalculus[Nabla].`#mover(mi("\`A__s\`"),mo("&rarr;"))` in (4.3) ? The answer is in the relationship established at the end of Sec I between the components of the tensor "`A__s`[]^(j)"and the components of the vector `#mover(mi("\`A__s\`"),mo("&rarr;"))`, namely that the vector components are obtained multiplying the contravariant tensor components by the scale-factors h__j. So, in the above we need to substitute the contravariant "`A__s`[]^(j)" by the vector components A__j divided by the scale-factors

[seq(A__s[Library:-Contravariant(j)](S) = Component(A__s_, j)/h[j], j = 1 .. 3)]

[A__s[`~1`](S) = A__r(S), A__s[`~2`](S) = `A__&theta;`(S)/r, A__s[`~3`](S) = `A__&phi;`(S)/(r*sin(theta))]

(4.10)

subs[eval]([A__s[`~1`](S) = A__r(S), A__s[`~2`](S) = `A__&theta;`(S)/r, A__s[`~3`](S) = `A__&phi;`(S)/(r*sin(theta))], %D_[j](%A__s[`~j`]) = diff(A__s[`~1`](S), r)+diff(A__s[`~2`](S), theta)+diff(A__s[`~3`](S), phi)+2*A__s[`~1`](S)/r+cos(theta)*A__s[`~2`](S)/sin(theta))

%D_[j](%A__s[`~j`]) = diff(A__r(S), r)+(diff(`A__&theta;`(S), theta))/r+(diff(`A__&phi;`(S), phi))/(r*sin(theta))+2*A__r(S)/r+cos(theta)*`A__&theta;`(S)/(sin(theta)*r)

(4.11)

Comparing with (4.3), we see these two expressions are the same:

expand(%Divergence(%A__s_) = ((diff(A__r(S), r))*r+2*A__r(S))/r+((diff(`A__&theta;`(S), theta))*sin(theta)+`A__&theta;`(S)*cos(theta))/(r*sin(theta))+(diff(`A__&phi;`(S), phi))/(r*sin(theta)))

%Divergence(%A__s_) = diff(A__r(S), r)+(diff(`A__&theta;`(S), theta))/r+(diff(`A__&phi;`(S), phi))/(r*sin(theta))+2*A__r(S)/r+cos(theta)*`A__&theta;`(S)/(sin(theta)*r)

(4.12)
• 

The Curl

 

The Curl of the the vector `#mover(mi("\`A__s\`"),mo("&rarr;"))` in spherical coordinates is given by

%Curl(%A__s_) = Curl(A__s_)

%Curl(%A__s_) = ((diff(`A__&phi;`(S), theta))*sin(theta)+`A__&phi;`(S)*cos(theta)-(diff(`A__&theta;`(S), phi)))*_r/(r*sin(theta))+(diff(A__r(S), phi)-(diff(`A__&phi;`(S), r))*r*sin(theta)-`A__&phi;`(S)*sin(theta))*_theta/(r*sin(theta))+((diff(`A__&theta;`(S), r))*r+`A__&theta;`(S)-(diff(A__r(S), theta)))*_phi/r

(4.13)

 

One could think that the expression for the Curl in tensor notation is as in a non-curvilinear system

 

"`&epsilon;`[i,j,k] `&dtri;`[]^(j)(`A__s`[]^(k))"

 

But in a curvilinear system `&epsilon;`[i, j, k] is not a tensor, we need to use the non-Galilean form Epsilon[i, j, k] = sqrt(%g_[determinant])*`&epsilon;`[i, j, k], where %g_[determinant] is the determinant of the metric. Moreover, since the expression "Epsilon[i,j,k] `&dtri;`[]^(j)(`A__s`[]^(k))"has one free covariant index (the first one), to compare with the vectorial formula (4.12) this index also needs to be rewritten as a vector component as discussed at the end of Sec. I, using

A__j = A__j/h__j

The formula (4.13) for the vectorial Curl is thus expressed using tensor notation as

Setup(levicivita = nongalilean)

[levicivita = nongalilean]

(4.14)

%Curl(%A__s_) = LeviCivita[i, j, k]*D_[`~j`](A__s[`~k`](S))/%h[i]

%Curl(%A__s_) = Physics:-LeviCivita[i, j, k]*Physics:-D_[`~j`](A__s[`~k`](S), [S])/%h[i]

(4.15)

followed by replacing the contravariant tensor components "`A__s`[]^(k)" by the vector components A__k/h__k using (4.10). Proceeding the same way we did with the Divergence, expand this expression. We could use TensorArray , but Library:-TensorComponents places a comma between components making things more readable in this case

lhs(%Curl(%A__s_) = Physics[LeviCivita][i, j, k]*D_[`~j`](A__s[`~k`](S), [S])/%h[i]) = Library:-TensorComponents(rhs(%Curl(%A__s_) = Physics[LeviCivita][i, j, k]*D_[`~j`](A__s[`~k`](S), [S])/%h[i]))

%Curl(%A__s_) = [(sin(theta)^3*(diff(A__s[`~3`](S), theta))*r^2+2*sin(theta)^2*cos(theta)*A__s[`~3`](S)*r^2-(diff(A__s[`~2`](S), phi))*sin(theta)*r^2)/(%h[1]*sin(theta)^2*r^2), (-sin(theta)^3*(diff(A__s[`~3`](S), r))*r^4-2*sin(theta)^3*A__s[`~3`](S)*r^3+(diff(A__s[`~1`](S), phi))*sin(theta)*r^2)/(%h[2]*sin(theta)^2*r^2), (sin(theta)^3*(diff(A__s[`~2`](S), r))*r^4+2*sin(theta)^3*A__s[`~2`](S)*r^3-sin(theta)^3*(diff(A__s[`~1`](S), theta))*r^2)/(%h[3]*sin(theta)^2*r^2)]

(4.16)

Replace now the components of the tensor "`A__s`[]^(j)" by the components of the 3D vector `#mover(mi("\`A__s\`"),mo("&rarr;"))` using (4.10)

lhs(%Curl(%A__s_) = [(sin(theta)^3*(diff(A__s[`~3`](S), theta))*r^2+2*sin(theta)^2*cos(theta)*A__s[`~3`](S)*r^2-(diff(A__s[`~2`](S), phi))*sin(theta)*r^2)/(%h[1]*sin(theta)^2*r^2), (-sin(theta)^3*(diff(A__s[`~3`](S), r))*r^4-2*sin(theta)^3*A__s[`~3`](S)*r^3+(diff(A__s[`~1`](S), phi))*sin(theta)*r^2)/(%h[2]*sin(theta)^2*r^2), (sin(theta)^3*(diff(A__s[`~2`](S), r))*r^4+2*sin(theta)^3*A__s[`~2`](S)*r^3-sin(theta)^3*(diff(A__s[`~1`](S), theta))*r^2)/(%h[3]*sin(theta)^2*r^2)]) = value(subs[eval]([A__s[`~1`](S) = A__r(S), A__s[`~2`](S) = `A__&theta;`(S)/r, A__s[`~3`](S) = `A__&phi;`(S)/(r*sin(theta))], rhs(%Curl(%A__s_) = [(sin(theta)^3*(diff(A__s[`~3`](S), theta))*r^2+2*sin(theta)^2*cos(theta)*A__s[`~3`](S)*r^2-(diff(A__s[`~2`](S), phi))*sin(theta)*r^2)/(%h[1]*sin(theta)^2*r^2), (-sin(theta)^3*(diff(A__s[`~3`](S), r))*r^4-2*sin(theta)^3*A__s[`~3`](S)*r^3+(diff(A__s[`~1`](S), phi))*sin(theta)*r^2)/(%h[2]*sin(theta)^2*r^2), (sin(theta)^3*(diff(A__s[`~2`](S), r))*r^4+2*sin(theta)^3*A__s[`~2`](S)*r^3-sin(theta)^3*(diff(A__s[`~1`](S), theta))*r^2)/(%h[3]*sin(theta)^2*r^2)])))

%Curl(%A__s_) = [(sin(theta)^3*((diff(`A__&phi;`(S), theta))/(r*sin(theta))-`A__&phi;`(S)*cos(theta)/(r*sin(theta)^2))*r^2+2*sin(theta)*cos(theta)*`A__&phi;`(S)*r-(diff(`A__&theta;`(S), phi))*r*sin(theta))/(h[1]*sin(theta)^2*r^2), (-sin(theta)^3*((diff(`A__&phi;`(S), r))/(r*sin(theta))-`A__&phi;`(S)/(r^2*sin(theta)))*r^4-2*sin(theta)^2*`A__&phi;`(S)*r^2+(diff(A__r(S), phi))*sin(theta)*r^2)/(h[2]*sin(theta)^2*r^2), (sin(theta)^3*((diff(`A__&theta;`(S), r))/r-`A__&theta;`(S)/r^2)*r^4+2*sin(theta)^3*`A__&theta;`(S)*r^2-sin(theta)^3*(diff(A__r(S), theta))*r^2)/(h[3]*sin(theta)^2*r^2)]

(4.17)

(lhs = `@`(simplify, rhs))(%Curl(%A__s_) = [(sin(theta)^3*((diff(`A__&phi;`(S), theta))/(r*sin(theta))-`A__&phi;`(S)*cos(theta)/(r*sin(theta)^2))*r^2+2*sin(theta)*cos(theta)*`A__&phi;`(S)*r-(diff(`A__&theta;`(S), phi))*r*sin(theta))/(h[1]*sin(theta)^2*r^2), (-sin(theta)^3*((diff(`A__&phi;`(S), r))/(r*sin(theta))-`A__&phi;`(S)/(r^2*sin(theta)))*r^4-2*sin(theta)^2*`A__&phi;`(S)*r^2+(diff(A__r(S), phi))*sin(theta)*r^2)/(h[2]*sin(theta)^2*r^2), (sin(theta)^3*((diff(`A__&theta;`(S), r))/r-`A__&theta;`(S)/r^2)*r^4+2*sin(theta)^3*`A__&theta;`(S)*r^2-sin(theta)^3*(diff(A__r(S), theta))*r^2)/(h[3]*sin(theta)^2*r^2)])

%Curl(%A__s_) = [((diff(`A__&phi;`(S), theta))*sin(theta)+`A__&phi;`(S)*cos(theta)-(diff(`A__&theta;`(S), phi)))/(r*sin(theta)), (diff(A__r(S), phi)-(diff(`A__&phi;`(S), r))*r*sin(theta)-`A__&phi;`(S)*sin(theta))/(r*sin(theta)), ((diff(`A__&theta;`(S), r))*r+`A__&theta;`(S)-(diff(A__r(S), theta)))/r]

(4.18)

We see these are exactly the components of the Curl (4.13)

%Curl(%A__s_) = ((diff(`A__&phi;`(S), theta))*sin(theta)+`A__&phi;`(S)*cos(theta)-(diff(`A__&theta;`(S), phi)))*_r/(r*sin(theta))+(diff(A__r(S), phi)-(diff(`A__&phi;`(S), r))*r*sin(theta)-`A__&phi;`(S)*sin(theta))*_theta/(r*sin(theta))+((diff(`A__&theta;`(S), r))*r+`A__&theta;`(S)-(diff(A__r(S), theta)))*_phi/r

%Curl(%A__s_) = ((diff(`A__&phi;`(S), theta))*sin(theta)+`A__&phi;`(S)*cos(theta)-(diff(`A__&theta;`(S), phi)))*_r/(r*sin(theta))+(diff(A__r(S), phi)-(diff(`A__&phi;`(S), r))*r*sin(theta)-`A__&phi;`(S)*sin(theta))*_theta/(r*sin(theta))+((diff(`A__&theta;`(S), r))*r+`A__&theta;`(S)-(diff(A__r(S), theta)))*_phi/r

(4.19)
• 

The Gradient

 

Once the problem is fully understood, it is easy to redo the computations of Sec.III for the Gradient, this time using tensor notation and the covariant derivative. In tensor notation, the components of the Gradient are given by the components of the right-hand side

%Nabla(f(S)) = `&dtri;`[j](f(S))/%h[j]

%Nabla(f(S)) = Physics:-d_[j](f(S), [S])/%h[j]

(4.20)

where on the left-hand side we have the vectorial Nabla  differential operator and on the right-hand side, since f(S) is a scalar, the covariant derivative `&dtri;`[j](f) becomes the standard derivative `&PartialD;`[j](f).

lhs(%Nabla(f(S)) = Physics[d_][j](f(S), [S])/%h[j]) = eval(value(Library:-TensorComponents(rhs(%Nabla(f(S)) = Physics[d_][j](f(S), [S])/%h[j]))))

%Nabla(f(S)) = [Physics:-Vectors:-diff(f(S), r), (diff(f(S), theta))/r, (diff(f(S), phi))/(r*sin(theta))]

(4.21)

The above is the expected result (3.2)

%Nabla(f(S)) = (diff(f(S), r))*_r+(diff(f(S), theta))*_theta/r+(diff(f(S), phi))*_phi/(r*sin(theta))

%Nabla(f(S)) = (diff(f(S), r))*_r+(diff(f(S), theta))*_theta/r+(diff(f(S), phi))*_phi/(r*sin(theta))

(4.22)
• 

The Laplacian

 

Likewise we can compute the Laplacian directly as

%Laplacian(f(S)) = D_[j](D_[j](f(S)))

%Laplacian(f(S)) = Physics:-D_[j](Physics:-d_[`~j`](f(S), [S]), [S])

(4.23)

In this case there are no free indices nor tensor components to be rewritten as vector components, so there is no need for scale-factors. Summing over the repeated indices,

SumOverRepeatedIndices(%Laplacian(f(S)) = D_[j](Physics[d_][`~j`](f(S), [S]), [S]))

%Laplacian(f(S)) = Physics:-dAlembertian(f(S), [S])+2*(diff(f(S), r))/r+cos(theta)*(diff(f(S), theta))/(sin(theta)*r^2)

(4.24)

Evaluating the  Vectors:-Laplacian on the left-hand side,

value(%Laplacian(f(S)) = Physics[dAlembertian](f(S), [S])+2*(diff(f(S), r))/r+cos(theta)*(diff(f(S), theta))/(sin(theta)*r^2))

((diff(diff(f(S), r), r))*r+2*(diff(f(S), r)))/r+((diff(diff(f(S), theta), theta))*sin(theta)+cos(theta)*(diff(f(S), theta)))/(r^2*sin(theta))+(diff(diff(f(S), phi), phi))/(r^2*sin(theta)^2) = Physics:-dAlembertian(f(S), [S])+2*(diff(f(S), r))/r+cos(theta)*(diff(f(S), theta))/(sin(theta)*r^2)

(4.25)

On the right-hand side we see the dAlembertian , "`&square;`(f(S)),"in curvilinear coordinates; rewrite it using standard diff  derivatives and expand both sides of the equation for comparison

expand(convert(((diff(diff(f(S), r), r))*r+2*(diff(f(S), r)))/r+((diff(diff(f(S), theta), theta))*sin(theta)+cos(theta)*(diff(f(S), theta)))/(r^2*sin(theta))+(diff(diff(f(S), phi), phi))/(r^2*sin(theta)^2) = Physics[dAlembertian](f(S), [S])+2*(diff(f(S), r))/r+cos(theta)*(diff(f(S), theta))/(sin(theta)*r^2), diff))

diff(diff(f(S), r), r)+(diff(diff(f(S), theta), theta))/r^2+(diff(diff(f(S), phi), phi))/(r^2*sin(theta)^2)+2*(diff(f(S), r))/r+cos(theta)*(diff(f(S), theta))/(sin(theta)*r^2) = diff(diff(f(S), r), r)+(diff(diff(f(S), theta), theta))/r^2+(diff(diff(f(S), phi), phi))/(r^2*sin(theta)^2)+2*(diff(f(S), r))/r+cos(theta)*(diff(f(S), theta))/(sin(theta)*r^2)

(4.26)

This is an identity, the left and right hand sides are equal:

evalb(diff(diff(f(S), r), r)+(diff(diff(f(S), theta), theta))/r^2+(diff(diff(f(S), phi), phi))/(r^2*sin(theta)^2)+2*(diff(f(S), r))/r+cos(theta)*(diff(f(S), theta))/(sin(theta)*r^2) = diff(diff(f(S), r), r)+(diff(diff(f(S), theta), theta))/r^2+(diff(diff(f(S), phi), phi))/(r^2*sin(theta)^2)+2*(diff(f(S), r))/r+cos(theta)*(diff(f(S), theta))/(sin(theta)*r^2))

true

(4.27)


 

Download Vectors_and_Spherical_coordinates_in_tensor_notation.mw

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

Until now I have been reading Maple Help files on the MAPLE website.  For convenience and mark-up, I have often dowloaded the help files and printed them on paper, only to find that the text over-runs the margins, and is therefore annoyingly incomplete.  On reflection, this is not surprising as the content is formatted for internet/web display!

TIP:  Instead, go to the bottom of the MAPLE help webpage of interest, click on the "Download Help Document" link and so download the Maple file:  e.g. pdsolve-numeric.mw.  The helpfile can then be read (& executed) using MAPLE.

Melvin

So here's something silly but cool you can do with Maple while you're "working" from home.

  • Record a few seconds of your voice on a microphone that's close to your mouth (probably using a headset). This is your dry audio.
  • On your phone, record a single clap of your hands in an enclosed space, like your shower cubicle or a closet. Trim this audio to the clap, and the reverb created by your enclosed space. This is your impulse response.
  • Send both sound files to whatever computer you have Maple on.
  • Using AudioTools:-Convolution, convolve the dry audio with the impulse response . This your wet audio and should sound a little bit like your voice was recorded in your enclosed space.

Here's some code. I've also attached my dry audio, an impulse response recorded in my shower (yes, I stood inside my shower, closed the door, and recorded a single clap of my hands on my phone), and the resulting wet audio.

with( AudioTools ):
dry_audio := Read( "MaryHadALittleLamb_sc.wav" ):
impulse_response := Read( "clap_sc.wav" ):
wet_audio := Normalize( Convolution( dry_audio, impulse_response ) ):
Write("wet_audio.wav", wet_audio );

A full Maple worksheet is here.

AudioSamplesForReverb.zip

An expression sequence is the underlying data structure for lists, sets, and function call arguments in Maple. Conceptually, a list is just a sequence enclosed in "[" and "]", a set is a sequence (with no duplicate elements) enclosed in "{" and "}", and a function call is a sequence enclosed in "(" and ")". A sequence can also be used as a data structure itself:

> Q := x, 42, "string", 42;
                           Q := x, 42, "string", 42

> L := [ Q ];
                          L := [x, 42, "string", 42]

> S := { Q };
                            S := {42, "string", x}

> F := f( Q );
                          F := f(x, 42, "string", 42)

A sequence, like most data structures in Maple, is immutable. Once created, it cannot be changed. This means the same sequence can be shared by multiple data structures. In the example above, the list assigned to and the function call assigned to both share the same instance of the sequence assigned to . The set assigned to refers to a different sequence, one with the duplicate 42 removed, and sorted into a canonical order.

Appending an element to a sequence creates a new sequence. The original remains unaltered, and still referenced by the list and function call:

> Q := Q, a+b;
                        Q := x, 42, "string", 42, a + b

> L;
                             [x, 42, "string", 42]

> S;
                               {42, "string", x}

> F;
                            f(x, 42, "string", 42)

Because appending to a sequence creates a new sequence, building a long sequence by appending one element at a time is very inefficient in both time and space. Building a sequence of length this way creates sequences of lengths 1, 2, ..., -1, . The extra space used will eventually be reclaimed by Maple's garbage collector, but this takes time.

This leads to the subject of this article, which is how to create long sequences efficiently. For the remainder of this article, the sequence we will use is the Fibonacci numbers, which are defined as follows:

  • Fib(0) = 0
  • Fib(1) = 1
  • Fib() = Fib(-1) + Fib(-2) for all > 1

In a computer algebra system like Maple, the simplest way to generate individual members of this sequence is with a recursive function. This is also very efficient if option is used (and very inefficient if it is not; computing Fib() requires 2 Fib() - 1 calls, and Fib() grows exponentially):

> Fib := proc(N)
>     option remember;
>     if N = 0 then
>         0
>     elif N = 1 then
>         1
>     else
>         Fib(N-1) + Fib(N-2)
>     end if
> end proc:
> Fib(1);
                                       1

> Fib(2);
                                       1

> Fib(5);
                                       5

> Fib(10);
                                      55

> Fib(20);
                                     6765

> Fib(50);
                                  12586269025

> Fib(100);
                             354224848179261915075

> Fib(200);
                  280571172992510140037611932413038677189525

Let's start with the most straightforward, and most inefficient way to generate a sequence of the first 100 Fibonacci numbers, starting with an empty sequence and using a for-loop to append one member at a time. Part of the output has been elided below in the interests of saving space:

> Q := ();
                                     Q :=

> for i from 0 to 99 do
>     Q := Q, Fib(i)
> end do:
> Q;
0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584,

    4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811,

    ...

    51680708854858323072, 83621143489848422977, 135301852344706746049,

    218922995834555169026

As mentioned previously, this actually produces 100 sequences of lengths 1 to 100, of which 99 will (eventually) be recovered by the garbage collector. This method is O(2) (Big O Notation) in time and space, meaning that producing a sequence of 200 values this way will take 4 times the time and memory as a sequence of 100 values.

The traditional Maple wisdom is to use the seq function instead, which produces only the requested sequence, and no intermediate ones:

> Q := seq(Fib(i),i=0..99);
Q := 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597,

    2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418,

    ...

    51680708854858323072, 83621143489848422977, 135301852344706746049,

    218922995834555169026

This is O() in time and space; generating a sequence of 200 elements takes twice the time and memory required for a sequence of 100 elements.

As of Maple 2019, it is also possible to achieve O() performance by constructing a sequence directly using a for-expression, without the cost of constructing the intermediate sequences that a for-statement would incur:

> Q := (for i from 0 to 99 do Fib(i) end do);
Q := 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597,

    2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418,

    ...

    51680708854858323072, 83621143489848422977, 135301852344706746049,

    218922995834555169026

This method is especially useful when you wish to add a condition to the elements selected for the sequence, since the full capabilities of Maple loops can be used (see The Two Kinds of Loops in Maple). The following two examples produce a sequence containing only the odd members of the first 100 Fibonacci numbers, and the first 100 odd Fibonacci numbers respectively:

> Q := (for i from 0 to 99 do
>           f := Fib(i);
>           if f :: odd then
>               f
>           else
>               NULL
>           end if
>       end do);
Q := 1, 1, 3, 5, 13, 21, 55, 89, 233, 377, 987, 1597, 4181, 6765, 17711, 28657,

    75025, 121393, 317811, 514229, 1346269, 2178309, 5702887, 9227465,

    ...

    19740274219868223167, 31940434634990099905, 83621143489848422977,

    135301852344706746049

> count := 0:
> Q := (for i from 0 while count < 100 do
>           f := Fib(i);
>           if f :: odd then
>               count += 1;
>               f
>           else
>               NULL
>           end if
>       end do);
Q := 1, 1, 3, 5, 13, 21, 55, 89, 233, 377, 987, 1597, 4181, 6765, 17711, 28657,

    75025, 121393, 317811, 514229, 1346269, 2178309, 5702887, 9227465,

    ...

    898923707008479989274290850145, 1454489111232772683678306641953,

    3807901929474025356630904134051, 6161314747715278029583501626149

> i;
                                      150

A for-loop used as an expression generates a sequence, producing one member for each iteration of the loop. The value of that member is the last expression computed during the iteration. If the last expression in an iteration is NULL, no value is produced for that iteration.

Examining after the second loop completes, we can see that 149 Fibonacci numbers were generated to find the first 100 odd ones. (The loop control variable is incremented before the while condition is checked, hence is one more than the number of completed iterations.)

Until now, we've been using calls to the Fib function to generate the individual Fibonacci numbers. These numbers can of course also be generated by a simple loop which, together with assignment of its initial conditions, can be written as a single sequence:

> Q := ((f0 := 0),
>       (f1 := 1),
>       (for i from 2 to 99 do
>            f0, f1 := f1, f0 + f1;
>            f1
>        end do));
Q := 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597,

    2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418,

    ...

    51680708854858323072, 83621143489848422977, 135301852344706746049,

    218922995834555169026

A Maple Array is a mutable data structure. Changing an element of an Array modifies the Array in-place; no new copy is generated:

> A := Array([a,b,c]);
                                A := [a, b, c]

> A[2] := d;
                                   A[2] := d

> A;
                                   [a, d, c]

It is also possible to append elements to an array, either by using programmer indexing, or the recently introduced ,= operator:

> A(numelems(A)+1) := e; # () instead of [] denotes "programmer indexing"
                               A := [a, d, c, e]

> A;
                                 [a, d, c, e]

Like appending to a sequence, this sometimes causes the existing data to be discarded and new data to be allocated, but this is done in chunks proportional to the current size of the Array, resulting in time and memory usage that is still O(). This can be used to advantage to generate sequences efficiently:

> A := Array(0..1,[0,1]);
                              [ 0..1 1-D Array       ]
                         A := [ Data Type: anything  ]
                              [ Storage: rectangular ]
                              [ Order: Fortran_order ]

> for i from 2 to 99 do
>     A ,= A[i-1] + A[i-2]
> end do:
> A;
                           [ 0..99 1-D Array      ]
                           [ Data Type: anything  ]
                           [ Storage: rectangular ]
                           [ Order: Fortran_order ]

> Q := seq(A);
Q := 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597,

    2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418,

    ...

    51680708854858323072, 83621143489848422977, 135301852344706746049,

    218922995834555169026

Although unrelated specifically to the goal of producing sequences, the same techniques can be used to construct Maple strings efficiently:

> A := Array("0");
                                   A := [48]

> for i from 1 to 99 do
>    A ,= " ", String(Fib(i))
> end do:
> A;
                           [ 1..1150 1-D Array     ]
                           [ Data Type: integer[1] ]
                           [ Storage: rectangular  ]
                           [ Order: Fortran_order  ]

> A[1..10];
                   [48, 32, 49, 32, 49, 32, 50, 32, 51, 32]

> S := String(A);
S := "0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 \
    10946 17711 28657 46368 75025 121393 196418 317811 514229 832040 134626\
    9 2178309 3524578 5702887 9227465 14930352 24157817 39088169 63245986 1\
    02334155 165580141 267914296 433494437 701408733 1134903170 1836311903 \
    2971215073 4807526976 7778742049 12586269025 20365011074 32951280099 53\
    316291173 86267571272 139583862445 225851433717 365435296162 5912867298\
    79 956722026041 1548008755920 2504730781961 4052739537881 6557470319842\
     10610209857723 17167680177565 27777890035288 44945570212853 7272346024\
    8141 117669030460994 190392490709135 308061521170129 498454011879264 80\
    6515533049393 1304969544928657 2111485077978050 3416454622906707 552793\
    9700884757 8944394323791464 14472334024676221 23416728348467685 3788906\
    2373143906 61305790721611591 99194853094755497 160500643816367088 25969\
    5496911122585 420196140727489673 679891637638612258 1100087778366101931\
     1779979416004714189 2880067194370816120 4660046610375530309 7540113804\
    746346429 12200160415121876738 19740274219868223167 3194043463499009990\
    5 51680708854858323072 83621143489848422977 135301852344706746049 21892\
    2995834555169026"

A call to the Array constructor with a string as an argument produces an array of bytes (Maple data type integer[1]). The ,= operator can then be used to append additional characters or strings, with O() efficiency. Finally, the Array can be converted back into a Maple string.

Constructing sequences in Maple is a common operation when writing Maple programs. Maple gives you many ways to do this, and it's worthwhile taking the time to choose a method that is efficient, and suitable to the task at hand.

Maple 2020 offers many improvements motivated and driven by our users.

Every single update in a new release has a story behind it. It might be a new function that a customer wants, a response to some feedback about usability, or an itch that a developer needs to scratch.

I’ll end this post with a story about acoustic guitars and how they drove improvements in signal and audio processing. But first, here are some of my personal favorites from Maple 2020.

Graph theory is a big focus of Maple 2020. The new features include more control over visualization, additional special graphs, new analysis functions, and even an interactive layout tool.

I’m particularly enamoured by these:

  • We’ve introduced new centrality measures - these help you determine the most influential vertices, based on their connections to other vertices
  • You now have more control over the styling of graphs – for example, you can vary the size or color of a nodebased on its centrality

I’ve used these two new features to identify the most influential MaplePrimes users. Get the worksheet here.

@Carl Love – looks like you’re the biggest mover and shaker on MaplePrimes (well, according to the eigenvector centrality of the MaplePrimes interaction graph).

We’ve also started using graph theory elsewhere in Maple. For example, you can generate static call graph to visualize dependencies between procedures calls in a procedure

You now get smoother edges for 3d surfaces with non-numeric values. Just look at the difference between Maple 2019 and 2020 for this plot.

Printing and PDF export has gotten a whole lot better.  We’ve put a lot of work into the proper handling of plots, tables, and interactive components, so the results look better than before.

For example, plots now maintain their aspect ratio when printed. So your carefully constructed psychrometric chart will not be squashed and stretched when exported to a PDF.

We’ve overhauled the start page to give it a cleaner, less cluttered look – this is much more digestible for new users (experienced users might find the new look attractive as well!). There’s a link to the Maple Portal, and an updated Maple Fundamentals guide that helps new users learn the product.

We’ve also linked to a guide that helps you choose between Document and Worksheet, and a link to a new movie.

New messages also guide new users away from some very common mistakes. For example, students often type “e” when referring to the exponential constant – a warning now appears if that is detected

We’re always tweaking existing functions to make them faster. For example, you can now compute the natural logarithm of large integers much more quickly and with less memory.

This calculation is about 50 times faster in Maple 2020 than in prior versions:

Many of our educators have asked for this – the linear algebra tutorials now return step by step solutions to the main document, so you have a record of what you did after the tutor is closed.

Continuing with this theme, the Student:-LinearAlgebra context menu features several new linear algebra visualizations to the Student:-LinearAlgebra Context Menu. This, for example, is an eigenvector plot.

Maple can now numerically evaluate various integral transforms.

The numerical inversion of integral transforms has application in many branches of science and engineering.

Maple is the world’s best tool for the symbolic solution of ODEs and PDEs, and in each release we push the boundary back further.

For example, Maple 2020 has improved tools for find hypergeometric solutions for linear PDEs.

This might seem like a minor improvement that’s barely worth mentions, but it’s one I now use all the time! You can now reorder worksheet tabs just by clicking and dragging.

The Hough transform lets you detect straight lines and line segments in images.

Hough transforms are widely used in automatic lane detection systems for autonomous driving. You can even detect the straight lines on a Sudoku grid!

The Physics package is always a pleasure to write about because it's something we do far better than the competition.

The new explore option in TensorArray combines two themes in Maple - Physics and interactive components. It's an intuitive solution to the real problem of viewing the contents of higher dimensional tensorial expressions.

There are many more updates to Physics in Maple 2020, including a completely rewritten FeynmanDiagrams command.

The Quantum Chemistry Toolbox has been updated with more analysis tools and curriculum material.

There’s more teaching content for general chemistry.

Among the many new analysis functions, you can now visualize transition orbitals.

I promised you a story about acoustic guitars and Maple 2020, didn’t I?

I often start a perfectly innocuous conversation about Maple that descends into several weeks of intense, feverish work.

The work is partly for me, but mostly for my colleagues. They don’t like me for that.

That conversation usually happens on a Friday afternoon, when we’re least prepared for it. On the plus side, this often means a user has planted a germ of an idea for a new feature or improvement, and we just have to will it into existence.

One Friday afternoon last year, I was speaking to a user about acoustic guitars. He wanted to synthetically generate guitar chords with reverb, and export the sound to a 32-bit Wave file. All of this, in Maple.

This started a chain of events that that involved least-square filters, frequency response curves, convolution, Karplus-Strong string synthesis and more. We’ll package up the results of this work, and hand it over to you – our users – over the next one or two releases.

Let me tell you what made it into Maple 2020.

Start by listening to this:

It’s a guitar chord played twice, the second time with reverb, both generated with Maple.

The reverb was simulated with convolving the artificially generated guitar chord with an impulse response. I had a choice of convolution functions in the SignalProcessing and AudioTools packages.

Both gave the same results, but we found that SignalProcessing:-Convolution was much faster than its AudioTools counterpart.

There’s no reason for the speed difference, so R&D modified AudioTools:-Convolution to leverage SignalProcessing:-Convolution for the instances for which their options are compatible. In this application, AudioTools:-Convolution is 25 times faster in Maple 2020 than Maple 2019!

We also discovered that the underlying library we use for the SignalProcessing package (the Intel IPP) gives two options for convolution that we were previously not using; a method which use an explicit formula and a “fast” method that uses FFTs. We modified SignalProcessing:-Convolution to accept both options (previously, we used just one of the methods),

That’s the story behind two new features in Maple 2020. Look at the entirety of what’s new in this release – there’s a tale for each new feature. I’d love to tell you more, but I’d run out of ink before I finish.

To read about everything that’s new in Maple 2020, go to the new features page.

When discussing Maple programming, we often refer to for-loops, while-loops, until-loops, and do-loops (the latter being an infinite loop). But under the hood, Maple has only two kinds of loop, albeit very flexible and powerful ones that can combine the capabilities of any or all of the above, making it possible to write very concise code in a natural way.

Before looking at some actual examples, here is the formal definition of the loops' syntax, expressed in Wirth Syntax Notation, where "|" denotes alternatives, "[...]" denotes an optional part, "(...)" denotes grouping, and Maple keywords are in boldface:

[ for  ] [ from  ] [ by  ] [ to  ]
    [ while  ]
do
    
( end do | until  )
[ for  [ , variable ] ] in 
    [ while  ]
do
    
( end do | until  )

In the first form, every part of the loop syntax is optional, except the do keyword before the body of the loop, and either end do or an until clause after the body. (For those who prefer it, end do can also be written as od.) In the second form, only the in clause is required.

The simplest loop is just:

do
    
end do

This will repeat the forever, unless a break or return statement is executed, or an error occurs.

One or two loop termination conditions can be added:

  • A while clause can be written before the do, specifying a condition that is tested before each iteration begins. If the condition evaluates to false, the loop ends.
  • An until clause can be written instead of the end do, specifying a condition that is tested after each iteration finishes. If the condition evaluates to true, the loop ends.

A so-called for-loop is just a loop to which iteration clauses have been added. These can take one of two forms:

  • Any combination of for (with a single variable), from, by, and to clauses. The last three can appear in any order.
  • A for clause with one or two variables, followed by an in clause.

The following for-loop executes 10 times:

for  from 1 to 10 do
    
end do

However, if the doesn't depend on the value of , both the for and from clauses can be omitted:

to 10 do
    
end do

In this case, Maple supplies an implicit for clause (with an inaccessible internal variable), as well as an implicit "from 1" clause. In fact, all of the clauses are optional, and the infinite loop shown earlier is understood by Maple in exactly the same way as:

for  from 1 by 1 to infinity while true do
    
until false

When looping over the contents of a container, such as a one-dimensional array A, there are several possible approaches. The one closest to how it would be done in most other programming languages is (this example and those that follow can be copied and pasted into a Maple session):

 := Array([,"foo",42]);
for  from lowerbound() to upperbound() do
    print([],[])
end do;

If only the entries in the container are of interest, it is not necessary to loop over the indices. Instead, one can write:

 := Array([,"foo",42]);
for  in  do
    print()
end do;

If both the indices and values are needed, one can write:

 := Array([,"foo",42]);
for ,  in  do
    print([],)
end do;

For a numerically indexed container such as an Array, this is equivalent to the for-from-to example. However, this method also works with arbitrarily indexed containers such as a Matrix or table:

 := LinearAlgebra:-RandomMatrix(2,3);
for ,  in  do
    print([],)
end do;
 := table({1="one","hello"="world",=42});
for ,  in eval() do
    print([],)
end do;

(The second example requires the call to eval due to last-name evaluation of tables in Maple, a topic for another post.)

As with a simple do-loop, a while and/or until clause can be added. For example, the following finds the first negative entry, if any, in a Matrix (traversing the Matrix in storage order):

 := LinearAlgebra:-RandomMatrix(2,3);
for ,  in  do
    # nothing to do here
until  < 0;
if  < 0 then
    print([],)
end if;

Notice that the test, < 0, is written twice, since it is possible that the Matrix has no negative entry. Another way to write the same loop but only perform the test once is as follows:

 := LinearAlgebra:-RandomMatrix(2,3);
for ,  in  do
    if  < 0 then
	print([],);
	break
    end if;
end do;

Here, we perform the test within the loop, perform the desired processing on the found value (just printing in this case), and use a break statement to terminate the loop.

Sometimes, it is useful to abort the current iteration of the loop and move on to the next one. The next statement does exactly that. The following loop prints all the indices but only the positive values in a Matrix:

 := LinearAlgebra:-RandomMatrix(2,3);
for ,  in  do
    print(=[]);
    if  < 0 then
	next
    end if;
    print(=);
end do;

(Note that a simple example like this would be better written by enclosing the printing of the value in an if-statement instead of using next. The latter is generally only used if the former is not possible.)

Maple's loop statements are very flexible and powerful, making it possible to write loops with complex combinations of termination conditions in a concise yet readable way. The ability to use while and/or until in conjunction with for means that break statements are often unnecessary, further improving clarity.

Playing mini-golf recently, I realized that my protractor can only help me so far since it can't calculate the speed of the swing needed.  I decided a more sophisticated tool was needed and modeled a trick-shot in MapleSim.

To start, I laid out the obstacles, the ball and club, the ground, and some additional visualizations in the MapleSim environment.

 

When running the simulation, my first result wasn't even close to the hole (similar to when I play in real life!).

 

The model clearly needed to be optimized. I went to the Optimization app in MapleSim (this can be found under Add Apps or Templates  on the left hand side).

 

Inside the app I clicked "Load System" then selected the parameters I wanted to optimize.

 

For this case, I'm optimizing 's' (the speed of the club) and 'theta' (the angle of the club). For the Objective Function I added a Relative Translation Sensor to the model and attached a probe to the Vector Norm of the output.

 

Inside the app, I switched to the Objective Function section.  Selecting Probes, I added the new probe as the Objective Function by giving it a weight of 1.

 

 

Scrolling down to "Execute Parameter Optimization", I checked the "Use Global Optimization Toolbox" checkbox, and clicked Run Parameter Optimization.

 

Following a run time of 120 seconds, the app returns the graph of the objective function. 

 

Below the plot, optimal values for the parameters are given. Plugging these back into the parameter block for the simulation we see that the ball does in fact go into the hole. Success!

 

 

Mini_golf_Global_Optimization.msim

 

Feynman Diagrams
The scattering matrix in coordinates and momentum representation

  

Mathematical methods for particle physics was one of the weak spots in the Physics package. There existed a FeynmanDiagrams command, but its capabilities were too minimal. People working in the area asked for more functionality. These diagrams are the cornerstone of calculations in particle physics (collisions involving from the electron to the Higgs boson), for example at the CERN. As an introduction for people curious, not working in the area, see "Why Feynman Diagrams are so important".

  

This post is thus about a new development in Physics: a full rewriting of the FeynmanDiagrams command, now including a myriad of new capabilities (mainly a. b. and c. in the Introduction), reversing the previous status of things entirely. This is work in collaboration with Davide Polvara from Durham University, Centre for Particle Theory.

  


The complexity of this material is high, so the introduction to the presentation below is as brief as it can get, emphasizing the examples instead. This material is reproducible in Maple 2019.2 after installing the Physics Updates, v.598 or higher.

  

 

  

At the end they are attached the worksheet corresponding to this presentation and a PDF version of it, as well as the new FeynmanDiagrams help page with all the explanatory details.

Introduction

  

A scattering matrix S relates the initial and final states, `#mfenced(mrow(mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&verbar;",close = "&rang;")` and `#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;")),open = "&verbar;",close = "&rang;")`, of an interacting system. In an 4-dimensional spacetime with coordinates X, S can be written as:

S = T(exp(i*`#mrow(mo("&int;"),mi("L"),mo("&ApplyFunction;"),mfenced(mi("X")),mo("&DifferentialD;"),msup(mi("X"),mn("4")))`))

  

where i is the imaginary unit  and L is the interaction Lagrangian, written in terms of quantum fields  depending on the spacetime coordinates  X. The T symbol means time-ordered. For the terminology used in this page, see for instance chapter IV, "The Scattering Matrix", of ref.[1] Bogoliubov, N.N., and Shirkov, D.V. Quantum Fields.

  

This exponential can be expanded as

S = 1+S[1]+S[2]+S[3]+`...`

  

where

S[n] = `#mrow(mo("&ApplyFunction;"),mfrac(msup(mi("i"),mi("n")),mrow(mi("n"),mo("&excl;")),linethickness = "1"),mo("&InvisibleTimes;"),mo("&int;"),mi("&hellip;"),mo("&InvisibleTimes;"),mo("&int;"),mi("T"),mo("&ApplyFunction;"),mfenced(mrow(mi("L"),mo("&ApplyFunction;"),mfenced(mi("\`X__1\`")),mo("&comma;"),mi("&hellip;"),mo("&comma;"),mi("L"),mo("&ApplyFunction;"),mfenced(mi("\`X__n\`")))),mo("&InvisibleTimes;"),mo("&DifferentialD;"),msup(mi("\`X__1\`"),mn("4")),mo("&InvisibleTimes;"),mi("&hellip;"),mo("&InvisibleTimes;"),mo("&DifferentialD;"),msup(mi("\`X__n\`"),mn("4")))`

  

and T(L(X[1]), `...`, L(X[n])) is the time-ordered product of n interaction Lagrangians evaluated at different points. The S matrix formulation is at the core of perturbative approaches in relativistic Quantum Field Theory.

  

In connection, the FeynmanDiagrams  command has been rewritten entirely for Maple 2020. In brief, the new functionality includes computing:

a. 

The expansion S = 1+S[1]+S[2]+S[3]+`...` in coordinates representation up to arbitrary order (the limitation is now only your hardware)

b. 

The S-matrix element `#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` in momentum representation up to arbitrary order for given number of loops and initial and final particles (the contents of the `#mfenced(mrow(mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&verbar;",close = "&rang;")` and `#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;")),open = "&verbar;",close = "&rang;")` states); optionally, also the transition probability density, constructed using the square of the scattering matrix element abs(`#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")`)^2, as shown in formula (13) of sec. 21.1 of ref.[1].

c. 

The Feynman diagrams (drawings) related to the different terms of the expansion of S or of its matrix elements `#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")`.

  

Interaction Lagrangians involving derivatives of fields, typically appearing in non-Abelian gauge theories, are also handled, and several options are provided enabling restricting the outcome in different ways, regarding the incoming and outgoing particles, the number of loops, vertices or external legs, the propagators and normal products, or whether to compute tadpoles and 1-particle reducible terms.

 

Examples

 

For illustration purposes set three coordinate systems , and set phi to represent a quantum operator

with(Physics)

Setup(mathematicalnotation = true, coordinates = [X, Y, Z], quantumoperators = phi)

`Systems of spacetime coordinates are:`*{X = (x1, x2, x3, x4), Y = (y1, y2, y3, y4), Z = (z1, z2, z3, z4)}

 

_______________________________________________________

 

[coordinatesystems = {X, Y, Z}, mathematicalnotation = true, quantumoperators = {phi}]

(1.1)

Let L be the interaction Lagrangian

L := lambda*phi(X)^4

lambda*Physics:-`^`(phi(X), 4)

(1.2)

The expansion of S in coordinates representation, computed by default up to order = 3 (you can change that using the option order = n), by definition containing all possible configurations of external legs, displaying the related Feynman Diagrams, is given by

%eval(S, `=`(order, 3)) = FeynmanDiagrams(L, diagrams)

 

 

 

%eval(S, order = 3) = 1+%FeynmanIntegral(lambda*_GF(_NP(phi(X), phi(X), phi(X), phi(X))), [[X]])+%FeynmanIntegral(16*lambda^2*_GF(_NP(phi(X), phi(X), phi(X), phi(Y), phi(Y), phi(Y)), [[phi(X), phi(Y)]])+96*lambda^2*_GF(_NP(phi(X), phi(Y)), [[phi(X), phi(Y)], [phi(X), phi(Y)], [phi(X), phi(Y)]])+72*lambda^2*_GF(_NP(phi(X), phi(X), phi(Y), phi(Y)), [[phi(X), phi(Y)], [phi(X), phi(Y)]]), [[X], [Y]])+%FeynmanIntegral(1728*lambda^3*_GF(_NP(phi(X), phi(X), phi(Y), phi(Y), phi(Z), phi(Z)), [[phi(X), phi(Z)], [phi(X), phi(Y)], [phi(Z), phi(Y)]])+2592*lambda^3*_GF(_NP(phi(X), phi(X), phi(Y), phi(Y)), [[phi(X), phi(Z)], [phi(X), phi(Z)], [phi(Z), phi(Y)], [phi(Z), phi(Y)]])+10368*lambda^3*_GF(_NP(phi(X), phi(Y), phi(Z), phi(Z)), [[phi(X), phi(Y)], [phi(X), phi(Y)], [phi(X), phi(Z)], [phi(Y), phi(Z)]])+10368*lambda^3*_GF(_NP(phi(X), phi(Y)), [[phi(X), phi(Y)], [phi(X), phi(Z)], [phi(X), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)]])+3456*lambda^3*_GF(_NP(phi(X), phi(X)), [[phi(X), phi(Y)], [phi(X), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)]])+576*lambda^3*_GF(_NP(phi(X), phi(X), phi(X), phi(Y), phi(Y), phi(Z), phi(Z), phi(Z)), [[phi(X), phi(Y)], [phi(Y), phi(Z)]]), [[X], [Y], [Z]])

(1.3)


The expansion of S  in coordinates representation to a specific order shows in a compact way the topology of the underlying Feynman diagrams. Each integral is represented with a new command, FeynmanIntegral , that works both in coordinates and momentum representation. To each term of the integrands corresponds a diagram, and the correspondence is always clear from the symmetry factors.

In a typical situation, one wants to compute a specific term, or scattering process, instead of the S matrix up to some order with all possible configurations of external legs. For example, to compute only the terms of this result that correspond to diagrams with 1 loop use numberofloops = 1 (for tree-level, use numerofloops = 0)

%eval(S, [`=`(order, 3), `=`(loops, 1)]) = FeynmanDiagrams(L, numberofloops = 1, diagrams)

%eval(S, [order = 3, loops = 1]) = %FeynmanIntegral(72*lambda^2*_GF(_NP(phi(X), phi(X), phi(Y), phi(Y)), [[phi(X), phi(Y)], [phi(X), phi(Y)]]), [[X], [Y]])+%FeynmanIntegral(1728*lambda^3*_GF(_NP(phi(X), phi(X), phi(Y), phi(Y), phi(Z), phi(Z)), [[phi(X), phi(Z)], [phi(X), phi(Y)], [phi(Z), phi(Y)]]), [[X], [Y], [Z]])

(1.4)


In the result above there are two terms, with 4 and 6 external legs respectively.

A scattering process with matrix element `#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` in momentum representation, corresponding to the term with 4 external legs (symmetry factor = 72), could be any process where the total number of incoming + outgoing parties is equal to 4. For example, one with 2 incoming and 2 outgoing particles. The transition probability for that process is given by

`#mfenced(mrow(mo("&InvisibleTimes;"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&comma;",mathcolor = "olive"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&comma;",mathcolor = "olive"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&InvisibleTimes;",mathcolor = "olive")),open = "&lang;",close = "&rang;")` = FeynmanDiagrams(L, incomingparticles = [phi, phi], outgoingparticles = [phi, phi], numberofloops = 1, diagrams)

 

`#mfenced(mrow(mo("&InvisibleTimes;"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&comma;",mathcolor = "olive"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&comma;",mathcolor = "olive"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&InvisibleTimes;",mathcolor = "olive")),open = "&lang;",close = "&rang;")` = %FeynmanIntegral((9/8)*lambda^2*Dirac(-P__3-P__4+P__1+P__2)/(Pi^6*(E__1*E__2*E__3*E__4)^(1/2)*(p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1-P__2-p__2)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])+%FeynmanIntegral((9/8)*lambda^2*Dirac(-P__3-P__4+P__1+P__2)/(Pi^6*(E__1*E__2*E__3*E__4)^(1/2)*(p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1+P__3-p__2)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])+%FeynmanIntegral((9/8)*lambda^2*Dirac(-P__3-P__4+P__1+P__2)/(Pi^6*(E__1*E__2*E__3*E__4)^(1/2)*(p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1+P__4-p__2)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])

(1.5)

When computing in momentum representation, only the topology of the corresponding Feynman diagrams is shown (i.e. the diagrams associated to the corresponding Feynman integral in coordinates representation).

The transition matrix element `#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` is related to the transition probability density dw (formula (13) of sec. 21.1 of ref.[1]) by

dw = (2*Pi)^(3*s-4)*n__1*`...`*n__s*abs(F(p[i], p[f]))^2*delta(sum(p[i], i = 1 .. s)-(sum(p[f], f = 1 .. r)))*` d `^3*p[1]*` ...`*`d `^3*p[r]

where n__1*`...`*n__s represent the particle densities of each of the s particles in the initial state `#mfenced(mrow(mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&verbar;",close = "&rang;")`, the delta (Dirac) is the expected singular factor due to the conservation of the energy-momentum and the amplitude F(p[i], p[f])is related to `#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` via

`#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` = F(p[i], p[f])*delta(sum(p[i], i = 1 .. s)-(sum(p[f], f = 1 .. r)))

To directly get the probability density dw instead of`#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")`use the option output = probabilitydensity

FeynmanDiagrams(L, incomingparticles = [phi, phi], outgoingparticles = [phi, phi], numberofloops = 1, output = probabilitydensity)

Physics:-FeynmanDiagrams:-ProbabilityDensity(4*Pi^2*%mul(n[i], i = 1 .. 2)*abs(F)^2*Dirac(-P__3-P__4+P__1+P__2)*%mul(dP_[f]^3, f = 1 .. 2), F = %FeynmanIntegral((9/8)*lambda^2/(Pi^6*(E__1*E__2*E__3*E__4)^(1/2)*(p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1-P__2-p__2)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])+%FeynmanIntegral((9/8)*lambda^2/(Pi^6*(E__1*E__2*E__3*E__4)^(1/2)*(p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1+P__3-p__2)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])+%FeynmanIntegral((9/8)*lambda^2/(Pi^6*(E__1*E__2*E__3*E__4)^(1/2)*(p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1+P__4-p__2)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]]))

(1.6)

In practice, the most common computations involve processes with 2 or 4 external legs. To restrict the expansion of the scattering matrix in coordinates representation to that kind of processes use the numberofexternallegs option. For example, the following computes the expansion of S up to order = 3, restricting the outcome to the terms corresponding to diagrams with only 2 external legs

%eval(S, [`=`(order, 3), `=`(legs, 2)]) = FeynmanDiagrams(L, numberofexternallegs = 2, diagrams)

%eval(S, [order = 3, legs = 2]) = %FeynmanIntegral(96*lambda^2*_GF(_NP(phi(X), phi(Y)), [[phi(X), phi(Y)], [phi(X), phi(Y)], [phi(X), phi(Y)]]), [[X], [Y]])+%FeynmanIntegral(3456*lambda^3*_GF(_NP(phi(X), phi(X)), [[phi(X), phi(Y)], [phi(X), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)]])+10368*lambda^3*_GF(_NP(phi(X), phi(Y)), [[phi(X), phi(Y)], [phi(X), phi(Z)], [phi(X), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)]]), [[X], [Y], [Z]])

(1.7)


This result shows two Feynman integrals, with respectively 2 and 3 loops, the second integral with two terms. The transition probability density in momentum representation for a process related to the first integral (1 term with symmetry factor = 96) is then

FeynmanDiagrams(L, incomingparticles = [phi], outgoingparticles = [phi], numberofloops = 2, diagrams, output = probabilitydensity)

Physics:-FeynmanDiagrams:-ProbabilityDensity((1/2)*%mul(n[i], i = 1 .. 1)*abs(F)^2*Dirac(-P__2+P__1)*%mul(dP_[f]^3, f = 1 .. 1)/Pi, F = %FeynmanIntegral(%FeynmanIntegral(((3/8)*I)*lambda^2/(Pi^7*(E__1*E__2)^(1/2)*(p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*(p__3^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1-p__2-p__3)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]]), [[p__3]]))

(1.8)

In the above, for readability, the contracted spacetime indices in the square of momenta entering the amplitude F (as denominators of propagators) are implicit. To make those indices explicit, use the option putindicesinsquareofmomentum

F = FeynmanDiagrams(L, incoming = [phi], outgoing = [phi], numberofloops = 2, indices)

`* Partial match of  '`*indices*`' against keyword '`*putindicesinsquareofmomentum*`' `

 

F = %FeynmanIntegral(%FeynmanIntegral(((3/8)*I)*lambda^2*Dirac(-P__2[`~kappa`]+P__1[`~kappa`])/(Pi^7*(E__1*E__2)^(1/2)*(p__2[mu]*p__2[`~mu`]-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*(p__3[nu]*p__3[`~nu`]-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1[beta]-p__2[beta]-p__3[beta])*(-P__1[`~beta`]-p__2[`~beta`]-p__3[`~beta`])-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]]), [[p__3]])

(1.9)

This computation can also be performed to higher orders. For example, with 3 loops, in coordinates and momentum representations, corresponding to the other two terms and diagrams in (1.7)

%eval(S[3], [`=`(legs, 2), `=`(loops, 3)]) = FeynmanDiagrams(L, legs = 2, loops = 3)

`* Partial match of  '`*legs*`' against keyword '`*numberoflegs*`' `

 

`* Partial match of  '`*loops*`' against keyword '`*numberofloops*`' `

 

%eval(S[3], [legs = 2, loops = 3]) = %FeynmanIntegral(3456*lambda^3*_GF(_NP(phi(X), phi(X)), [[phi(X), phi(Y)], [phi(X), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)]])+10368*lambda^3*_GF(_NP(phi(X), phi(Y)), [[phi(X), phi(Y)], [phi(X), phi(Z)], [phi(X), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)]]), [[X], [Y], [Z]])

(1.10)

A corresponding S-matrix element in momentum representation:

%eval(%Bracket(phi, S[3], phi), `=`(loops, 3)) = FeynmanDiagrams(L, incomingparticles = [phi], outgoingparticles = [phi], numberofloops = 3)

%eval(%Bracket(phi, S[3], phi), loops = 3) = %FeynmanIntegral(%FeynmanIntegral(%FeynmanIntegral((9/32)*lambda^3*Dirac(-P__2+P__1)/(Pi^11*(E__1*E__2)^(1/2)*(p__3^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*(p__4^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*(p__5^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-p__3-p__4-p__5)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1+P__2+p__3+p__4+p__5)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__3]]), [[p__4]]), [[p__5]])+2*%FeynmanIntegral(%FeynmanIntegral(%FeynmanIntegral((9/32)*lambda^3*Dirac(-P__2+P__1)/(Pi^11*(E__1*E__2)^(1/2)*(p__3^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*(p__4^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*(p__5^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-p__3-p__4-p__5)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1+p__4+p__5)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__3]]), [[p__4]]), [[p__5]])+%FeynmanIntegral(%FeynmanIntegral((1/2048)*lambda*Dirac(-P__2+P__1)*%FeynmanIntegral(576*lambda^2/((p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-p__2-p__4-p__5)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])/(Pi^11*(E__1*E__2)^(1/2)*(p__4^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*(p__5^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1+p__4+p__5)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__4]]), [[p__5]])

(1.11)

Consider the interaction Lagrangian of Quantum Electrodynamics (QED). To formulate this problem on the worksheet, start defining the vector field A[mu].

Define(A[mu])

`Defined objects with tensor properties`

 

{A[mu], Physics:-Dgamma[mu], P__1[mu], P__2[mu], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], p__1[mu], p__2[mu], p__3[mu], p__4[mu], p__5[mu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X), Physics:-SpaceTimeVector[mu](Y), Physics:-SpaceTimeVector[mu](Z)}

(1.12)

Set lowercase Latin letters from i to s to represent spinor indices (you can change this setting according to your preference, see Setup ), also the (anticommutative) spinor field will be represented by psi, so set psi as an anticommutativeprefix, and set A and psi as quantum operators

Setup(spinorindices = lowercaselatin_is, anticommutativeprefix = psi, op = {A, psi})

`* Partial match of  '`*op*`' against keyword '`*quantumoperators*`' `

 

_______________________________________________________

 

[anticommutativeprefix = {psi}, quantumoperators = {A, phi, psi}, spinorindices = lowercaselatin_is]

(1.13)

The matrix indices of the Dirac matrices  are written explicitly and use conjugate  to represent the Dirac conjugate conjugate(psi[j])

L__QED := alpha*conjugate(psi[j](X))*Dgamma[mu][j, k]*psi[k](X)*A[mu](X)

alpha*Physics:-`*`(conjugate(psi[j](X)), psi[k](X), A[mu](X))*Physics:-Dgamma[`~mu`][j, k]

(1.14)

Compute S[2], only the terms with 4 external legs, and display the diagrams: all the corresponding graphs have no loops

%eval(S[2], `=`(legs, 4)) = FeynmanDiagrams(L__QED, numberofvertices = 2, numberoflegs = 4, diagrams)

%eval(S[2], legs = 4) = %FeynmanIntegral(-2*alpha^2*Physics:-Dgamma[`~mu`][j, k]*Physics:-Dgamma[`~alpha`][i, l]*_GF(_NP(psi[k](X), A[mu](X), conjugate(psi[i](Y)), A[alpha](Y)), [[psi[l](Y), conjugate(psi[j](X))]])+alpha^2*Physics:-Dgamma[`~mu`][j, k]*Physics:-Dgamma[`~alpha`][i, l]*_GF(_NP(conjugate(psi[j](X)), psi[k](X), conjugate(psi[i](Y)), psi[l](Y)), [[A[mu](X), A[alpha](Y)]]), [[X], [Y]])

(1.15)

The same computation but with only 2 external legs results in the diagrams with 1 loop that correspond to the self-energy of the electron and the photon (page 218 of ref.[1])

%eval(S[2], `=`(legs, 2)) = FeynmanDiagrams(L__QED, numberofvertices = 2, numberoflegs = 2, diagrams)

 

 

%eval(S[2], legs = 2) = %FeynmanIntegral(-2*alpha^2*Physics:-Dgamma[`~mu`][j, k]*Physics:-Dgamma[`~alpha`][i, l]*_GF(_NP(psi[k](X), conjugate(psi[i](Y))), [[A[mu](X), A[alpha](Y)], [psi[l](Y), conjugate(psi[j](X))]])-alpha^2*Physics:-Dgamma[`~mu`][j, k]*Physics:-Dgamma[`~alpha`][i, l]*_GF(_NP(A[mu](X), A[alpha](Y)), [[psi[l](Y), conjugate(psi[j](X))], [psi[k](X), conjugate(psi[i](Y))]]), [[X], [Y]])

(1.16)

where the diagram with two spinor legs is the electron self-energy. To restrict the output furthermore, for example getting only the self-energy of the photon, you can specify the normal products you want:

%eval(S[2], [`=`(legs, 2), `=`(products, _NP(A, A))]) = FeynmanDiagrams(L__QED, numberofvertices = 2, numberoflegs = 2, normalproduct = _NP(A, A))

`* Partial match of  '`*normalproduct*`' against keyword '`*normalproducts*`' `

 

%eval(S[2], [legs = 2, products = _NP(A, A)]) = %FeynmanIntegral(alpha^2*Physics:-Dgamma[`~mu`][j, k]*Physics:-Dgamma[`~alpha`][i, l]*_GF(_NP(A[mu](X), A[alpha](Y)), [[conjugate(psi[j](X)), psi[l](Y)], [psi[k](X), conjugate(psi[i](Y))]]), [[X], [Y]])

(1.17)

The corresponding S-matrix elements in momentum representation

`#mfenced(mrow(mo("&InvisibleTimes;"),mi("&psi;",fontstyle = "normal"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("&psi;",fontstyle = "normal"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` = FeynmanDiagrams(L__QED, incomingparticles = [psi], outgoing = [psi], numberofloops = 1, diagrams)

 

`#mfenced(mrow(mo("&InvisibleTimes;"),mi("&psi;",fontstyle = "normal"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("&psi;",fontstyle = "normal"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` = -%FeynmanIntegral((1/8)*Physics:-FeynmanDiagrams:-Uspinor[psi][i](P__1_)*conjugate(Physics:-FeynmanDiagrams:-Uspinor[psi][l](P__2_))*(-Physics:-g_[alpha, nu]+p__2[nu]*p__2[alpha]/m__A^2)*alpha^2*Physics:-Dgamma[`~alpha`][l, m]*Physics:-Dgamma[`~nu`][n, i]*((P__1[beta]+p__2[beta])*Physics:-Dgamma[`~beta`][m, n]+m__psi*Physics:-KroneckerDelta[m, n])*Dirac(-P__2+P__1)/(Pi^3*(p__2^2-m__A^2+I*Physics:-FeynmanDiagrams:-epsilon)*((P__1+p__2)^2-m__psi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])

(1.18)


In this result we see u[psi] spinor (see ref.[2]), and the propagator of the field A[mu] with a mass m[A]. To indicate that this field is massless use

Setup(massless = A)

`* Partial match of  '`*massless*`' against keyword '`*masslessfields*`' `

 

_______________________________________________________

 

[masslessfields = {A}]

(1.19)

Now the propagator for A[mu] is the one of a massless vector field:

FeynmanDiagrams(L__QED, incoming = [psi], outgoing = [psi], numberofloops = 1)

-%FeynmanIntegral(-(1/8)*Physics:-FeynmanDiagrams:-Uspinor[psi][i](P__1_)*conjugate(Physics:-FeynmanDiagrams:-Uspinor[psi][l](P__2_))*Physics:-g_[alpha, nu]*alpha^2*Physics:-Dgamma[`~alpha`][l, m]*Physics:-Dgamma[`~nu`][n, i]*((P__1[beta]+p__2[beta])*Physics:-Dgamma[`~beta`][m, n]+m__psi*Physics:-KroneckerDelta[m, n])*Dirac(-P__2+P__1)/(Pi^3*(p__2^2+I*Physics:-FeynmanDiagrams:-epsilon)*((P__1+p__2)^2-m__psi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])

(1.20)

The self-energy of the photon:

`#mfenced(mrow(mo("&InvisibleTimes;"),mi("A"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("A"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` = FeynmanDiagrams(L__QED, incomingparticles = [A], outgoing = [A], numberofloops = 1)

`#mfenced(mrow(mo("&InvisibleTimes;"),mi("A"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("A"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` = -%FeynmanIntegral((1/16)*Physics:-FeynmanDiagrams:-PolarizationVector[A][nu](P__1_)*conjugate(Physics:-FeynmanDiagrams:-PolarizationVector[A][alpha](P__2_))*(m__psi*Physics:-KroneckerDelta[l, n]+p__2[beta]*Physics:-Dgamma[`~beta`][l, n])*alpha^2*Physics:-Dgamma[`~alpha`][n, i]*Physics:-Dgamma[`~nu`][m, l]*((P__1[tau]+p__2[tau])*Physics:-Dgamma[`~tau`][i, m]+m__psi*Physics:-KroneckerDelta[i, m])*Dirac(-P__2+P__1)/(Pi^3*(E__1*E__2)^(1/2)*(p__2^2-m__psi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((P__1+p__2)^2-m__psi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])

(1.21)

where epsilon[A] is the corresponding polarization vector.

When working with non-Abelian gauge fields, the interaction Lagrangian involves derivatives. FeynmanDiagrams  can handle that kind of interaction in momentum representation. Consider for instance a Yang-Mills theory with a massless field B[mu, a] where a is a SU2 index (see eq.(12) of sec. 19.4 of ref.[1]). The interaction Lagrangian can be entered as follows

Setup(su2indices = lowercaselatin_ah, massless = B, op = B)

`* Partial match of  '`*massless*`' against keyword '`*masslessfields*`' `

 

`* Partial match of  '`*op*`' against keyword '`*quantumoperators*`' `

 

_______________________________________________________

 

[masslessfields = {A, B}, quantumoperators = {A, B, phi, psi, psi1}, su2indices = lowercaselatin_ah]

(1.22)

Define(B[mu, a], quiet)

F__B[mu, nu, a] := d_[mu](B[nu, a](X))-d_[nu](B[mu, a](X))

Physics:-d_[mu](B[nu, a](X), [X])-Physics:-d_[nu](B[mu, a](X), [X])

(1.23)

L := (1/2)*g*LeviCivita[a, b, c]*F__B[mu, nu, a]*B[mu, b](X)*B[nu, c](X)+(1/4)*g^2*LeviCivita[a, b, c]*LeviCivita[a, e, f]*B[mu, b](X)*B[nu, c](X)*B[mu, e](X)*B[nu, f](X)

(1/2)*g*Physics:-LeviCivita[a, b, c]*Physics:-`*`(Physics:-d_[mu](B[nu, a](X), [X])-Physics:-d_[nu](B[mu, a](X), [X]), B[`~mu`, b](X), B[`~nu`, c](X))+(1/4)*g^2*Physics:-LeviCivita[a, b, c]*Physics:-LeviCivita[a, e, f]*Physics:-`*`(B[mu, b](X), B[nu, c](X), B[`~mu`, e](X), B[`~nu`, f](X))

(1.24)

The transition probability density at tree-level for a process with two incoming and two outgoing B particles is given by

FeynmanDiagrams(L, incomingparticles = [B, B], outgoingparticles = [B, B], numberofloops = 0, output = probabilitydensity, factor, diagrams)

`* Partial match of  '`*factor*`' against keyword '`*factortreelevel*`' `

(1.25)

 

 

Physics:-FeynmanDiagrams:-ProbabilityDensity(4*Pi^2*%mul(n[i], i = 1 .. 2)*abs(F)^2*Dirac(-P__3[`~sigma`]-P__4[`~sigma`]+P__1[`~sigma`]+P__2[`~sigma`])*%mul(dP_[f]^3, f = 1 .. 2), F = (((1/8)*I)*Physics:-LeviCivita[a1, a3, h]*((-P__1[`~kappa`]-P__2[`~kappa`]-P__4[`~kappa`])*Physics:-g_[`~lambda`, `~tau`]+(P__1[`~lambda`]+P__2[`~lambda`]+P__3[`~lambda`])*Physics:-g_[`~kappa`, `~tau`]-Physics:-g_[`~kappa`, `~lambda`]*(P__3[`~tau`]-P__4[`~tau`]))*Physics:-LeviCivita[a2, d, g]*((P__1[`~beta`]+(1/2)*P__2[`~beta`])*Physics:-g_[`~alpha`, `~sigma`]+(-(1/2)*P__1[`~sigma`]+(1/2)*P__2[`~sigma`])*Physics:-g_[`~alpha`, `~beta`]-(1/2)*Physics:-g_[`~beta`, `~sigma`]*(P__1[`~alpha`]+2*P__2[`~alpha`]))*Physics:-g_[sigma, tau]*Physics:-KroneckerDelta[a2, a3]/((-P__1[chi]-P__2[chi])*(-P__1[`~chi`]-P__2[`~chi`])+I*Physics:-FeynmanDiagrams:-epsilon)-((1/16)*I)*((-P__1[`~beta`]+P__3[`~beta`]-P__4[`~beta`])*Physics:-g_[`~lambda`, `~tau`]+(P__1[`~lambda`]-P__2[`~lambda`]-P__3[`~lambda`])*Physics:-g_[`~beta`, `~tau`]+Physics:-g_[`~beta`, `~lambda`]*(P__2[`~tau`]+P__4[`~tau`]))*Physics:-LeviCivita[a1, a3, g]*((P__1[`~sigma`]+P__3[`~sigma`])*Physics:-g_[`~alpha`, `~kappa`]+(-2*P__1[`~kappa`]+P__3[`~kappa`])*Physics:-g_[`~alpha`, `~sigma`]+Physics:-g_[`~kappa`, `~sigma`]*(P__1[`~alpha`]-2*P__3[`~alpha`]))*Physics:-LeviCivita[a2, d, h]*Physics:-g_[sigma, tau]*Physics:-KroneckerDelta[a2, a3]/((-P__1[chi]+P__3[chi])*(-P__1[`~chi`]+P__3[`~chi`])+I*Physics:-FeynmanDiagrams:-epsilon)-((1/16)*I)*((-P__1[`~beta`]-P__3[`~beta`]+P__4[`~beta`])*Physics:-g_[`~kappa`, `~tau`]+(P__1[`~kappa`]-P__2[`~kappa`]-P__4[`~kappa`])*Physics:-g_[`~beta`, `~tau`]+Physics:-g_[`~beta`, `~kappa`]*(P__2[`~tau`]+P__3[`~tau`]))*Physics:-LeviCivita[a3, g, h]*((P__1[`~sigma`]+P__4[`~sigma`])*Physics:-g_[`~alpha`, `~lambda`]+(P__1[`~alpha`]-2*P__4[`~alpha`])*Physics:-g_[`~lambda`, `~sigma`]-2*Physics:-g_[`~alpha`, `~sigma`]*(P__1[`~lambda`]-(1/2)*P__4[`~lambda`]))*Physics:-LeviCivita[a1, a2, d]*Physics:-g_[sigma, tau]*Physics:-KroneckerDelta[a2, a3]/((-P__1[chi]+P__4[chi])*(-P__1[`~chi`]+P__4[`~chi`])+I*Physics:-FeynmanDiagrams:-epsilon)-((1/16)*I)*(Physics:-KroneckerDelta[g, h]*Physics:-KroneckerDelta[a1, d]*(Physics:-g_[`~alpha`, `~beta`]*Physics:-g_[`~kappa`, `~lambda`]+Physics:-g_[`~alpha`, `~kappa`]*Physics:-g_[`~beta`, `~lambda`]-2*Physics:-g_[`~alpha`, `~lambda`]*Physics:-g_[`~beta`, `~kappa`])+Physics:-KroneckerDelta[d, h]*(Physics:-g_[`~alpha`, `~beta`]*Physics:-g_[`~kappa`, `~lambda`]-2*Physics:-g_[`~alpha`, `~kappa`]*Physics:-g_[`~beta`, `~lambda`]+Physics:-g_[`~alpha`, `~lambda`]*Physics:-g_[`~beta`, `~kappa`])*Physics:-KroneckerDelta[a1, g]-2*(Physics:-g_[`~alpha`, `~beta`]*Physics:-g_[`~kappa`, `~lambda`]-(1/2)*Physics:-g_[`~beta`, `~kappa`]*Physics:-g_[`~alpha`, `~lambda`]-(1/2)*Physics:-g_[`~alpha`, `~kappa`]*Physics:-g_[`~beta`, `~lambda`])*Physics:-KroneckerDelta[d, g]*Physics:-KroneckerDelta[a1, h]))*g^2*conjugate(Physics:-FeynmanDiagrams:-PolarizationVector[B][kappa, h](P__3_))*conjugate(Physics:-FeynmanDiagrams:-PolarizationVector[B][lambda, a1](P__4_))*Physics:-FeynmanDiagrams:-PolarizationVector[B][alpha, d](P__1_)*Physics:-FeynmanDiagrams:-PolarizationVector[B][beta, g](P__2_)/(Pi^2*(E__1*E__2*E__3*E__4)^(1/2)))

(1.26)

To simplify the repeated indices, us the option simplifytensorindices. To check the indices entering a result like this one use Check ; there are no free indices, and regarding the repeated indices:

Check(Physics[FeynmanDiagrams]:-ProbabilityDensity(4*Pi^2*%mul(n[i], i = 1 .. 2)*abs(F)^2*Dirac(-P__3[`~sigma`]-P__4[`~sigma`]+P__1[`~sigma`]+P__2[`~sigma`])*%mul(dP_[f]^3, f = 1 .. 2), F = (((1/8)*I)*Physics[LeviCivita][a1, a3, h]*((-P__1[`~kappa`]-P__2[`~kappa`]-P__4[`~kappa`])*Physics[g_][`~lambda`, `~tau`]+(P__1[`~lambda`]+P__2[`~lambda`]+P__3[`~lambda`])*Physics[g_][`~kappa`, `~tau`]-Physics[g_][`~kappa`, `~lambda`]*(P__3[`~tau`]-P__4[`~tau`]))*Physics[LeviCivita][a2, d, g]*((P__1[`~beta`]+(1/2)*P__2[`~beta`])*Physics[g_][`~alpha`, `~sigma`]+(-(1/2)*P__1[`~sigma`]+(1/2)*P__2[`~sigma`])*Physics[g_][`~alpha`, `~beta`]-(1/2)*Physics[g_][`~beta`, `~sigma`]*(P__1[`~alpha`]+2*P__2[`~alpha`]))*Physics[g_][sigma, tau]*Physics[KroneckerDelta][a2, a3]/((-P__1[chi]-P__2[chi])*(-P__1[`~chi`]-P__2[`~chi`])+I*Physics[FeynmanDiagrams]:-epsilon)-((1/16)*I)*((-P__1[`~beta`]+P__3[`~beta`]-P__4[`~beta`])*Physics[g_][`~lambda`, `~tau`]+(P__1[`~lambda`]-P__2[`~lambda`]-P__3[`~lambda`])*Physics[g_][`~beta`, `~tau`]+Physics[g_][`~beta`, `~lambda`]*(P__2[`~tau`]+P__4[`~tau`]))*Physics[LeviCivita][a1, a3, g]*((P__1[`~sigma`]+P__3[`~sigma`])*Physics[g_][`~alpha`, `~kappa`]+(-2*P__1[`~kappa`]+P__3[`~kappa`])*Physics[g_][`~alpha`, `~sigma`]+Physics[g_][`~kappa`, `~sigma`]*(P__1[`~alpha`]-2*P__3[`~alpha`]))*Physics[LeviCivita][a2, d, h]*Physics[g_][sigma, tau]*Physics[KroneckerDelta][a2, a3]/((-P__1[chi]+P__3[chi])*(-P__1[`~chi`]+P__3[`~chi`])+I*Physics[FeynmanDiagrams]:-epsilon)-((1/16)*I)*((-P__1[`~beta`]-P__3[`~beta`]+P__4[`~beta`])*Physics[g_][`~kappa`, `~tau`]+(P__1[`~kappa`]-P__2[`~kappa`]-P__4[`~kappa`])*Physics[g_][`~beta`, `~tau`]+Physics[g_][`~beta`, `~kappa`]*(P__2[`~tau`]+P__3[`~tau`]))*Physics[LeviCivita][a3, g, h]*((P__1[`~sigma`]+P__4[`~sigma`])*Physics[g_][`~alpha`, `~lambda`]+(P__1[`~alpha`]-2*P__4[`~alpha`])*Physics[g_][`~lambda`, `~sigma`]-2*Physics[g_][`~alpha`, `~sigma`]*(P__1[`~lambda`]-(1/2)*P__4[`~lambda`]))*Physics[LeviCivita][a1, a2, d]*Physics[g_][sigma, tau]*Physics[KroneckerDelta][a2, a3]/((-P__1[chi]+P__4[chi])*(-P__1[`~chi`]+P__4[`~chi`])+I*Physics[FeynmanDiagrams]:-epsilon)-((1/16)*I)*(Physics[KroneckerDelta][g, h]*Physics[KroneckerDelta][a1, d]*(Physics[g_][`~alpha`, `~beta`]*Physics[g_][`~kappa`, `~lambda`]+Physics[g_][`~alpha`, `~kappa`]*Physics[g_][`~beta`, `~lambda`]-2*Physics[g_][`~alpha`, `~lambda`]*Physics[g_][`~beta`, `~kappa`])+Physics[KroneckerDelta][d, h]*(Physics[g_][`~alpha`, `~beta`]*Physics[g_][`~kappa`, `~lambda`]-2*Physics[g_][`~alpha`, `~kappa`]*Physics[g_][`~beta`, `~lambda`]+Physics[g_][`~alpha`, `~lambda`]*Physics[g_][`~beta`, `~kappa`])*Physics[KroneckerDelta][a1, g]-2*(Physics[g_][`~alpha`, `~beta`]*Physics[g_][`~kappa`, `~lambda`]-(1/2)*Physics[g_][`~alpha`, `~lambda`]*Physics[g_][`~beta`, `~kappa`]-(1/2)*Physics[g_][`~alpha`, `~kappa`]*Physics[g_][`~beta`, `~lambda`])*Physics[KroneckerDelta][d, g]*Physics[KroneckerDelta][a1, h]))*g^2*conjugate(Physics[FeynmanDiagrams]:-PolarizationVector[B][kappa, h](P__3_))*conjugate(Physics[FeynmanDiagrams]:-PolarizationVector[B][lambda, a1](P__4_))*Physics[FeynmanDiagrams]:-PolarizationVector[B][alpha, d](P__1_)*Physics[FeynmanDiagrams]:-PolarizationVector[B][beta, g](P__2_)/(Pi^2*(E__1*E__2*E__3*E__4)^(1/2))), all)

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

 

[{a1, a2, a3, alpha, beta, chi, d, g, h, kappa, lambda, sigma, tau}], {}

(1.27)


This process can be computed with 1 or more loops, in which case the number of terms increases significantly. As another interesting non-Abelian model, consider the interaction Lagrangian of the electro-weak part of the Standard Model

Coordinates(clear, Z)

`Unaliasing `*{Z}*` previously defined as a system of spacetime coordinates`

(1.28)

Setup(quantumoperators = {W, Z})

[quantumoperators = {A, B, W, Z, phi, psi, psi1}]

(1.29)

Define(W[mu], Z[mu])

`Defined objects with tensor properties`

 

{A[mu], B[mu, a], Physics:-Dgamma[mu], P__1[mu], P__2[mu], P__3[alpha], P__4[alpha], Physics:-Psigma[mu], W[mu], Z[mu], Physics:-d_[mu], Physics:-g_[mu, nu], p__1[mu], p__2[mu], p__3[mu], p__4[mu], p__5[mu], psi[j], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X), Physics:-SpaceTimeVector[mu](Y)}

(1.30)

CompactDisplay((W, Z)(X))

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

 

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

(1.31)

F__W[mu, nu] := d_[mu](W[nu](X))-d_[nu](W[mu](X))

Physics:-d_[mu](W[nu](X), [X])-Physics:-d_[nu](W[mu](X), [X])

(1.32)

F__Z[mu, nu] := d_[mu](Z[nu](X))-d_[nu](Z[mu](X))

Physics:-d_[mu](Z[nu](X), [X])-Physics:-d_[nu](Z[mu](X), [X])

(1.33)

L__WZ := I*g*cos(`&theta;__w`)*((Dagger(F__W[mu, nu])*W[mu](X)-Dagger(W[mu](X))*F__W[mu, nu])*Z[nu](X)+W[nu](X)*Dagger(W[mu](X))*F__Z[mu, nu])

I*g*cos(theta__w)*(Physics:-`*`(Physics:-`*`(Physics:-d_[mu](Physics:-Dagger(W[nu](X)), [X])-Physics:-d_[nu](Physics:-Dagger(W[mu](X)), [X]), W[`~mu`](X))-Physics:-`*`(Physics:-Dagger(W[mu](X)), Physics:-d_[`~mu`](W[nu](X), [X])-Physics:-d_[nu](W[`~mu`](X), [X])), Z[`~nu`](X))+Physics:-`*`(W[nu](X), Physics:-Dagger(W[mu](X)), Physics:-d_[`~mu`](Z[`~nu`](X), [X])-Physics:-d_[`~nu`](Z[`~mu`](X), [X])))

(1.34)

This interaction Lagrangian contains six different terms. The S-matrix element for the tree-level process with two incoming and two outgoing W particles is shown in the help page for FeynmanDiagrams .

NULL

References

 

[1] Bogoliubov, N.N., and Shirkov, D.V. Quantum Fields. Benjamin Cummings, 1982.

[2] Weinberg, S., The Quantum Theory Of Fields. Cambridge University Press, 2005.

 

FeynmanDiagrams_and_the_Scattering_Matrix.PDF

FeynmanDiagrams_and_the_Scattering_Matrix.mw

FeynmanDiagrams_-_help_page.mw


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

The ideas here are to allow 3D plotting commands such as plot3d to handle a `size` option similarly to how 2D plotting commands do so, and for the plots:-display command to also handle it for 3D plots.

The size denotes the dimensions of the inlined plotting window, and not the relative lengths of the three axes.

I'd be interested in any new problems introduced with this, eg. export, etc.

restart;

#
# Using ToInert/FromInert
#
# This might go in an initialzation file.
#
try
  __ver:=(u->:-parse(u[7..:-StringTools:-Search(",",u)-1]))(:-sprintf("%s",:-kernelopts(':-version')));
  if __ver>=18.0 and __ver<=2019.2 then
    __KO:=:-kernelopts(':-opaquemodules'=false);
    :-unprotect(:-Plot:-Options:-Verify:-ProcessParameters);
    __KK:=ToInert(eval(:-Plot:-Options:-Verify:-ProcessParameters)):
    __LL:=[:-op([1,2,1,..],__KK)]:
    __NN:=:-nops(:-remove(:-type,[:-op([1,..],__KK)],
                          ':-specfunc(:-_Inert_SET)'))
          +:-select(:-has,[:-seq([__i,__LL[__i]],
                                 __i=1..:-nops(__LL))],
                    "size")[1][1];
    if :-has(:-op([5,2,2,2,1],__KK),:-_Inert_PARAM(__NN)) then
      __KK:=:-subsop([5,2,2,2,1]
                     =:-subs([:-_Inert_PARAM(__NN)=:-NULL],
                              :-op([5,2,2,2,1],__KK)),__KK);
      :-Plot:-Options:-Verify:-ProcessParameters:=
      :-FromInert(:-subsop([5,2,2,3,1]
                  =:-subs([:-_Inert_STRING("size")=:-NULL],
                          :-op([5,2,2,3,1],__KK)),__KK));
      :-print("3D size patch done");
    else
      :-print("3D size patch not appropriate; possibly already done");
    end if;
  else
    :-print(sprintf("3D size patch not appropriate for version %a"),__ver);
  end if;
catch:
  :-print("3D size patch failed");
finally
  :-protect(:-Plot:-Options:-Verify:-ProcessParameters);
  :-kernelopts(':-opaquemodules'=__KO);
end try:

"3D size patch done"

 

P := plot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1, size=[150,150],
            font=[Times,5], labels=["","",""]):
P;

plots:-display(P, size=[300,300], font=[Times,10]);

#
# inherited from the contourplot3d (the plot3d is unset).
#
plots:-display(
  plots:-contourplot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1,
                       thickness=3, contours=20, size=[800,800]),
  plot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1, color="Gray",
         transparency=0.1, style=surface)
);

# Some options should still act as 2D-plot-specific.
#
try plot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1, legend="Q");
    print("Not OK");
catch:
if StringTools:-FormatMessage(lastexception[2..-1])
   ="the legend option is not available for 3-D plots"
then print("OK"); else print("Not OK"); error; end if; end try;

"OK"

 

Download 3Dsize_hotedit.mw

If this works fine then it might be a candidate for inclusion in an initialization file, so that it's
automatically available.

I was trying to display a Physics[Vectors] vector name in a 3dplot with an up arrow
on it. I found that this old 2008 trick still works in MAPLE 2018.

 


 

restart;

with(plots):
with(Physics[Vectors]):

# Using MAPLE 2018.2

a:=arrow([-1,1,1],view=[-1.5..1.5,-1.5..1.5,-1.5..1.5]):

v_;
t:= textplot3d([-1.1,1.1,1,v_]):
display(a,t);

v_

 

 

# I found this on an old 2008 post
t:= textplot3d([-1.1,1.1,1,typeset(`#mover(mi(` || v ||  `),mo("→"))`)]):
display(a,t);

 


 

Download VectorTypeSetting.mw

This update fixes the problems inadvertently introduced in Maple 2019.2, namely:

  • Maple failed to run the code in the maple.ini/.mapleinit initialization files when loading existing worksheets containing a restart() command
  • Installing some packages from the MapleCloud was unsuccessful

For anyone who installed the 2019.2 update, installing 2019.2.1 will fix these problems.

If you are at Maple 2019.1 or earlier, installing this update will bring you straight to Maple 2019.2.1.

This update is available through Tools>Check for Updates in Maple, and is also available from our website on the Maple 2019.2.1 download page.

If you are a MapleSim user, please note that these problems do not affect your use of MapleSim. If you use Maple on its own, and if you use Maple command initialization files and/or you need to install a package from the MapleCloud that does not work, please contact Maplesoft Technical Support for assistance.

We sincerely apologize for the inconvenience and thank you for your patience as we worked through this issue.

I'm only just hearing (haven't experienced) about some serious issues with the 2019.2 updates.  I would recommend waiting for Maplesoft to release an emergency 2019.3 fix update - Maplesoft can NOT leave the last update of 2019 in this state.

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