Product Tips & Techniques

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

MacDude posted a worksheet for creating help files using the makehelp command that I found very useful.  However, his worksheet created a table of contents that organized the command help pages under a single folder.  I wanted to create a help file table of contents with the commands organized into sub-folders under the main application folder. ie. Package Name,Folder Name, command. The help page for makehelp is not very informative; in particular the example that shows (purportedly) how to override the existing help pages with your own help file is very misleading.  Also, the description of the parameters to the browser option of the makehelp command is too vague. I needed an error message to tell me that the browser option expects parameters in the form List(Name,String).  In the end, I was able to get a folder/sub-folder structure using a structure [`Folder Name`,`Subfolder Name`,"Command"].  Sub-folders are sorted alphabetically. If anyone has a need, the attached worksheet shows how I created the table of contents structure. 

helpcode.mw

I like tweaking plots to get the look and feel I want, and luckily Maple has many plotting options that I often play with. Here, I visualize the same data several times, but each time with different styling.

First, some data.

restart:
data_1 := [[0,0],[1,2],[2,1.3],[3,6]]:
data_2 := [[0.5,3],[1,1],[2,5],[3,2]]:
data_3 := [[-0.5,3],[1.3,1],[2.5,5],[4.5,2]]:

This is the default look.

plot([data_1, data_2, data_3])

I think the darker background on this plot makes it easier to look at.

gray_grid :=
 background      = "LightGrey"
,color           = [ ColorTools:-Color("RGB",[150/255, 40 /255, 27 /255])
                    ,ColorTools:-Color("RGB",[0  /255, 0  /255, 0  /255])
                    ,ColorTools:-Color("RGB",[68 /255, 108/255, 179/255]) ]
,axes            = frame
,axis[2]         = [color = black, gridlines = [10, thickness = 1, color = ColorTools:-Color("RGB", [1, 1, 1])]]
,axis[1]         = [color = black, gridlines = [10, thickness = 1, color = ColorTools:-Color("RGB", [1, 1, 1])]]
,axesfont        = [Arial]
,labelfont       = [Arial]
,size            = [400*1.78, 400]
,labeldirections = [horizontal, vertical]
,filled          = false
,transparency    = 0
,thickness       = 5
,style           = line:

plot([data_1, data_2, data_3], gray_grid);

I call the next style Excel, for obvious reasons.

excel :=
 background      = white
,color           = [ ColorTools:-Color("RGB",[79/255,  129/255, 189/255])
                    ,ColorTools:-Color("RGB",[192/255, 80/255,   77/255])
                    ,ColorTools:-Color("RGB",[155/255, 187/255,  89/255])]
,axes            = frame
,axis[2]         = [gridlines = [10, thickness = 0, color = ColorTools:-Color("RGB",[134/255,134/255,134/255])]]
,font            = [Calibri]
,labelfont       = [Calibri]
,size            = [400*1.78, 400]
,labeldirections = [horizontal, vertical]
,transparency    = 0
,thickness       = 3
,style           = point
,symbol          = [soliddiamond, solidbox, solidcircle]
,symbolsize      = 15:

plot([data_1, data_2, data_3], excel)

This style makes the plot look a bit like the oscilloscope I have in my garage.

dark_gridlines :=
 background      = ColorTools:-Color("RGB",[0,0,0])
,color           = white
,axes            = frame
,linestyle       = [solid, dash, dashdot]
,axis            = [gridlines = [10, linestyle = dot, color = ColorTools:-Color("RGB",[0.5, 0.5, 0.5])]]
,font            = [Arial]
,size            = [400*1.78, 400]:

plot([data_1, data_2, data_3], dark_gridlines);

The colors in the next style remind me of an Autumn morning.

autumnal :=
 background      =  ColorTools:-Color("RGB",[236/255, 240/255, 241/255])
,color           = [  ColorTools:-Color("RGB",[144/255, 54/255, 24/255])
                     ,ColorTools:-Color("RGB",[105/255, 108/255, 51/255])
                     ,ColorTools:-Color("RGB",[131/255, 112/255, 82/255]) ]
,axes            = frame
,font            = [Arial]
,size            = [400*1.78, 400]
,filled          = true
,axis[2]         = [gridlines = [10, thickness = 1, color = white]]
,axis[1]         = [gridlines = [10, thickness = 1, color = white]]
,symbol          = solidcircle
,style           = point
,transparency    = [0.6, 0.4, 0.2]:

plot([data_1, data_2, data_3], autumnal);

In honor of a friend and ex-colleague, I call this style "The Swedish".

swedish :=
 background      = ColorTools:-Color("RGB", [0/255, 107/255, 168/255])
,color           = [ ColorTools:-Color("RGB",[169/255, 158/255, 112/255])
                    ,ColorTools:-Color("RGB",[126/255,  24/255,   9/255])
                    ,ColorTools:-Color("RGB",[254/255, 205/255,   0/255])]
,axes            = frame
,axis            = [gridlines = [10, color = ColorTools:-Color("RGB",[134/255,134/255,134/255])]]
,font            = [Arial]
,size            = [400*1.78, 400]
,labeldirections = [horizontal, vertical]
,filled          = false
,thickness       = 10:

plot([data_1, data_2, data_3], swedish);

This looks like a plot from a journal article.

experimental_data_mono :=

background       = white
,color           = black
,axes            = box
,axis            = [gridlines = [linestyle = dot, color = ColorTools:-Color("RGB",[0.5, 0.5, 0.5])]]
,font            = [Arial, 11]
,legendstyle     = [font = [Arial, 11]]
,size            = [400, 400]
,labeldirections = [horizontal, vertical]
,style           = point
,symbol          = [solidcircle, solidbox, soliddiamond]
,symbolsize      = [15,15,20]:

plot([data_1, data_2, data_3], experimental_data_mono, legend = ["Annihilation", "Authority", "Acceptance"]);

If you're willing to tinker a little bit, you can add some real character and personality to your visualizations. Try it!

I'd also be very interested to learn what you find attractive in a plot - please do let me know.

We have just released an update to Maple, Maple 2020.1.

Maple 2020.1 includes corrections and improvement to the mathematics engine, export to PDF, MATLAB connectivity, support for Ubuntu 20.04, and more.  We recommend that all Maple 2020 users install these updates.

This update is available through Tools>Check for Updates in Maple, and is also available from our website on the Maple 2020.1 download page, where you can also find more details.

In particular, please note that this update includes a fix to the SMTLIB problem reported on MaplePrimes. Thanks for the feedback!

 

We’re excited to announce a new version of MapleSim! The MapleSim 2020 family of products lets you build and test models faster than ever, including faster simulations, powerful new tools for machine builders, and expanded modeling capabilities. Improvements include:

  • Faster results, with more efficient models, faster simulations, and more powerful design tools.
  • Powerful new features for machine builders, with new components, improved visualizations, and automation-focused connectivity tools that make it faster than ever to build and test digital twins.
  • Improved modeling capabilities, with an extensive collection of updates to components, libraries, and analysis tools.
  • More realistic machine visualizations with an expanded Kinematic Cam Generation App.
  • New product: MapleSim Insight, giving machine builders powerful, simulation-based debugging and 3-D visualization capabilities that connect directly to common automation platforms.
  • New add-on library:  MapleSim Ropes and Pulleys Library for the easy creation of winch and pulley systems as part of your machine development. 

See What’s New in MapleSim 2020 for more information about these and other improvements!
 

This is my attempt to produce a subplot within a larger plot for magnifying data in a small region, and putting that subplot into the white space of the figure.
Based on the questions: How to insert a plot into another plot? and Inset figure in Maple, I wrote a couple of procedures that create sub-plots and allow the user to place the subplot window as he/she chooses. This avoids the graininess issues mentioned by acer in the second link (and experienced by me).

So far, I only have this completed for point plots, but using acer's method of piecewise functions posted in the plotin2b.mw of the second article, with the subplot function being defined only if it satisfies your conditions, would allow the subplot generating procedure to be generalized easily enough. But the data I'm working with all point plots, so that's the example here.

The basic idea  is to use one procedure to create boxes, make tickmarks on the expanded region, and make tickmark labels, combine all of those into one graph. Then create scaled and shifted versions of the data series, then make graphs of those. Lastly, combine them all into one picture.

Hope this helps someone who has to do the same.

Mapleprimes isn't inserting the contents, but here is the download of the file: SubPlotBoxesandVectorDataSeries.mw

 

The equations of motion in curvilinear coordinates, tensor notation and Coriolis force

``

 

The formulation of the equations of motion of a particle is simple in Cartesian coordinates using vector notation. However, depending on the problem, for example when describing the motion of a particle as seen from a non-inertial system of references (e.g. a rotating planet, like earth), there is advantage in using curvilinear coordinates and also tensor notation. When the particle's movement is observed from such a rotating referential, we also see "acceleration" that is not due to any force but to the fact that the referential itself is accelerated. The material below discusses and formulates these topics, and derives the expression for the Coriolis and centripetal force in cylindrical coordinates as seen from a rotating system of references.

 

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

 

Vector notation

 

Generally speaking the equations of motion of a particle are easy to formulate: the position vector is a function of time, the velocity is its first derivative and the acceleration is its second derivative. For instance, in Cartesian coordinates

with(Physics); with(Vectors)

r_(t) = x(t)*_i+y(t)*_j+z(t)*_k

r_(t) = x(t)*_i+y(t)*_j+z(t)*_k

(1)

diff(r_(t) = x(t)*_i+y(t)*_j+z(t)*_k, t)

diff(r_(t), t) = (diff(x(t), t))*_i+(diff(y(t), t))*_j+(diff(z(t), t))*_k

(2)

diff(diff(r_(t), t) = (diff(x(t), t))*_i+(diff(y(t), t))*_j+(diff(z(t), t))*_k, t)

diff(diff(r_(t), t), t) = (diff(diff(x(t), t), t))*_i+(diff(diff(y(t), t), t))*_j+(diff(diff(z(t), t), t))*_k

(3)

Newton's 2nd law, that in an inertial system of references when there is force there is acceleration and viceversa, is then given by

F_(t) = m*lhs(diff(diff(r_(t), t), t) = (diff(diff(x(t), t), t))*_i+(diff(diff(y(t), t), t))*_j+(diff(diff(z(t), t), t))*_k)

F_(t) = m*(diff(diff(r_(t), t), t))

(4)

where `#mover(mi("F"),mo("→"))`(t) = F__x(t)*`#mover(mi("i"),mo("∧"))`+F__y(t)*`#mover(mi("j"),mo("∧"))`+F__z(t)*`#mover(mi("k"),mo("∧"))` represents the total force acting on the particle. This vectorial equation represents three second order differential equations which, for given initial conditions, can be integrated to arrive at a closed form expression for `#mover(mi("r"),mo("→"))`(t) as a function of t.

 

Tensor notation

 

In Cartesian coordinates, the tensorial form of the equations (4) is also straightforward. In a flat spacetime - Galilean system of references - the three space coordinates x, y, z form a 3D tensor, and so does its first derivate and the second one. Set the spacetime to be 3-dimensional and Euclidean and use lowercaselatin indices for the corresponding tensors

Setup(coordinates = cartesian, metric = Euclidean, dimension = 3, spacetimeindices = lowercaselatin)

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

 

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

 

_______________________________________________________

 

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

 

_______________________________________________________

 

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

 

_______________________________________________________

(5)

The position, velocity and acceleration vectors are expressed in tensor notation as done in (1), (2) and (3)

X[j](t)

(X)[j](t)

(6)

diff((X)[j](t), t)

Physics:-Vectors:-diff((Physics:-SpaceTimeVector[j](X))(t), t)

(7)

diff(Physics[Vectors]:-diff((Physics[SpaceTimeVector][j](X))(t), t), t)

Physics:-Vectors:-diff(Physics:-Vectors:-diff((Physics:-SpaceTimeVector[j](X))(t), t), t)

(8)

Setting a tensor F__j(t) to represent the three Cartesian components of the force

Define(F[j] = [F__x(t), F__y(t), F__z(t)])

`Defined objects with tensor properties`

 

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

(9)

Newton's 2nd law (4), now expressed in tensorial notation, is given by

F[j] = m*Physics[Vectors]:-diff(Physics[Vectors]:-diff((Physics[SpaceTimeVector][j](X))(t), t), t)

F[j] = m*(diff(diff((Physics:-SpaceTimeVector[j](X))(t), t), t))

(10)

The three differential equations behind this tensorial form of (4) are as expected

TensorArray(F[j] = m*(diff(diff((Physics[SpaceTimeVector][j](X))(t), t), t)), output = setofequations)

{F__x(t) = m*(diff(diff(x(t), t), t)), F__y(t) = m*(diff(diff(y(t), t), t)), F__z(t) = m*(diff(diff(z(t), t), t))}

(11)

Things are straightforward in Cartesian coordinates because the components of the line element `#mover(mi("dr"),mo("→"))` = dx*`#mover(mi("i"),mo("∧"))`+dy*`#mover(mi("j"),mo("∧"))`+dz*`#mover(mi("k"),mo("∧"))` are exactly the components of the tensor d(X[j])

TensorArray(d_(X[j]))

Array(%id = 18446744078354237310)

(12)

and so, the form factors (see related Mapleprimes post) are all equal to 1.

 

In the case of no external forces, `#mover(mi("F"),mo("→"))`(t) = 0 and 0 = F[j] and the equations of motion, whose solution are the trajectory, can be formulated as the path of minimal length between two points, a geodesic. In the case under consideration, because the spacetime is flat (Galilean) those two points lie on a plane, we are talking about Euclidean geometry, that information is encoded in the metric (the 3x3 identity matrix (5)), and the geodesic is a straight line. The differential equations of this geodesic are thus the equations of motion (11) with  `#mover(mi("F"),mo("→"))`(t) = 0, and can be computed using Geodesics

 

Geodesics(t)

[diff(diff(z(t), t), t) = 0, diff(diff(y(t), t), t) = 0, diff(diff(x(t), t), t) = 0]

(13)

 

Curvilinear coordinates

 

Vector notation

 

The form of these equations in the case of curvilinear coordinates, for example in cylindrical or spherical variables, is obtained performing a change of coordinates.

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

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

(14)

This change keeps the z axis unchanged, so the corresponding unit vector `#mover(mi("k"),mo("∧"))` remains unchanged.

Changing the basis and coordinates used to represent the position vector `#mover(mi("r"),mo("→"))`(t) = x(t)*`#mover(mi("i"),mo("∧"))`+y(t)*`#mover(mi("j"),mo("∧"))`+z(t)*`#mover(mi("k"),mo("∧"))`, it becomes

r_(t) = ChangeBasis(rhs(r_(t) = x(t)*_i+y(t)*_j+z(t)*_k), cylindrical, alsocomponents)

r_(t) = z(t)*_k+rho(t)*_rho(t)

(15)

where since in (1) the coordinates (x, y, z) depend on t, in (15), not just rho(t) and z(t) but also the unit vector `#mover(mi("ρ",fontstyle = "normal"),mo("∧"))`(t)depends on t. The velocity is computed as usual, differentiating

diff(r_(t) = z(t)*_k+rho(t)*_rho(t), t)

diff(r_(t), t) = (diff(z(t), t))*_k+(diff(rho(t), t))*_rho(t)+rho(t)*(diff(phi(t), t))*_phi(t)

(16)

The second term is due to the dependency of `#mover(mi("ρ",fontstyle = "normal"),mo("∧"))` on the coordinate phi together with the chain rule diff(`#mover(mi("ρ",fontstyle = "normal"),mo("∧"))`(t), t) = (diff(`#mover(mi("ρ",fontstyle = "normal"),mo("∧"))`, phi))*(diff(phi(t), t)) and (diff(`#mover(mi("ρ",fontstyle = "normal"),mo("∧"))`, phi))*(diff(phi(t), t)) = `#mover(mi("φ",fontstyle = "normal"),mo("∧"))`(t)*(diff(phi(t), t)). The dependency of curvilinear unit vectors on the coordinates is automatically taken into account when differentiating due to the Setup setting geometricdifferentiation = true.

 

For the acceleration,

diff(diff(r_(t), t) = (diff(z(t), t))*_k+(diff(rho(t), t))*_rho(t)+rho(t)*(diff(phi(t), t))*_phi(t), t)

diff(diff(r_(t), t), t) = _rho(t)*(diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2)+_phi(t)*(2*(diff(rho(t), t))*(diff(phi(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k

(17)

where the term involving (diff(phi(t), t))^2 comes from differentiating `#mover(mi("φ",fontstyle = "normal"),mo("∧"))`(t) in (16) taking into account the dependency of `#mover(mi("φ",fontstyle = "normal"),mo("∧"))` on the coordinate "phi." This result can also be obtained by directly changing variables in the acceleration diff(`#mover(mi("r"),mo("→"))`(t), t, t), in equation (3)

lhs(diff(diff(r_(t), t), t) = (diff(diff(x(t), t), t))*_i+(diff(diff(y(t), t), t))*_j+(diff(diff(z(t), t), t))*_k) = ChangeBasis(rhs(diff(diff(r_(t), t), t) = (diff(diff(x(t), t), t))*_i+(diff(diff(y(t), t), t))*_j+(diff(diff(z(t), t), t))*_k), cylindrical, alsocomponents)

diff(diff(r_(t), t), t) = _rho(t)*(diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2)+_phi(t)*(2*(diff(rho(t), t))*(diff(phi(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k

(18)

 

Newton's 2nd law becomes

F_(t) = m*rhs(diff(diff(r_(t), t), t) = _rho(t)*(diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2)+_phi(t)*(2*(diff(rho(t), t))*(diff(phi(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k)

F_(t) = m*(_rho(t)*(diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2)+_phi(t)*(2*(diff(rho(t), t))*(diff(phi(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k)

(19)

In the absence of external forces, equating to 0 the vector components (coefficients of the unit vectors) of the acceleration diff(`#mover(mi("r"),mo("→"))`(t), t, t)we get the system of differential equations in cylindrical coordinates whose solution is the trajectory of the particle expressed in the ("rho(t),phi(t),z(t))"

`~`[`=`]({coeffs(rhs(diff(diff(r_(t), t), t) = _rho(t)*(diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2)+_phi(t)*(2*(diff(rho(t), t))*(diff(phi(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k), [`#mover(mi("ρ",fontstyle = "normal"),mo("∧"))`(t), `#mover(mi("φ",fontstyle = "normal"),mo("∧"))`(t), `#mover(mi("k"),mo("∧"))`])}, 0)

{2*(diff(rho(t), t))*(diff(phi(t), t))+rho(t)*(diff(diff(phi(t), t), t)) = 0, diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2 = 0, diff(diff(z(t), t), t) = 0}

(20)

solve({2*(diff(rho(t), t))*(diff(phi(t), t))+rho(t)*(diff(diff(phi(t), t), t)) = 0, diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2 = 0, diff(diff(z(t), t), t) = 0}, {diff(phi(t), t, t), diff(rho(t), t, t), diff(z(t), t, t)})

{diff(diff(phi(t), t), t) = -2*(diff(rho(t), t))*(diff(phi(t), t))/rho(t), diff(diff(rho(t), t), t) = rho(t)*(diff(phi(t), t))^2, diff(diff(z(t), t), t) = 0}

(21)

In this formulation (21) with `#mover(mi("F"),mo("→"))`(t) = 0, although diff(z(t), t, t) = 0, no acceleration in the `#mover(mi("k"),mo("∧"))` direction, is naturally expected, the same cannot be said about the other two equations for diff(phi(t), t, t) and diff(rho(t), t, t). Those two equations are discussed below under Coriolis and Centripetal forces. The key observation at this point, however, is that the right-hand sides of both unexpected equations involve diff(phi(t), t), rotation around the z axis.

 

Tensor notation

 

The same equations (19) and (21) result when using tensor notation. For that purpose, one can transform the position, velocity and acceleration tensors (6), (7), (8), but since they are expressed as functions of a parameter (the time), it is simpler to transform only the underlying metric using TransformCoordinates. That has the advantage that all the geometrical subtleties of curvilinear coordinates, like scale factors and dependency of unit vectors on curvilinear coordinates, get automatically, very succinctly, encoded in the metric:

TransformCoordinates(tr, g_[j, k], [rho, phi, z], setmetric)

_______________________________________________________

 

`Coordinates: `[rho, phi, z]*`. Signature: `(`+ + +`)

 

_______________________________________________________

 

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

 

_______________________________________________________

(22)

The computation of geodesics assumes that the coordinates (rho, phi, z) depend on a parameter. That parameter is passed as the first argument to Geodesics

Geodesics(t)

[diff(diff(rho(t), t), t) = rho(t)*(diff(phi(t), t))^2, diff(diff(phi(t), t), t) = -2*(diff(rho(t), t))*(diff(phi(t), t))/rho(t), diff(diff(z(t), t), t) = 0]

(23)

These equations of motion (23) are the same as the equations (21) computed using standard vector notation, differentiating and taking into account the dependency of curvilinear unit vectors on the curvilinear coordinates in (16) and (17).  One of the interesting features of computing with tensors is, as said, that all those geometrical algebraic subtleties of curvilinear coordinates are automatically encoded in the metric (22).

 

To understand how are the geodesic equations computed in one go in (23), one can perform the calculation in steps:

1. 

Make rho be a function of t directly in the metric

2. 

Compute - not the final form of the equations (23) - but the intermediate form expressing the geodesic equation using tensor notation, in terms of Christoffel symbols

3. 

Compute the components of that tensorial equation for the geodesic (using TensorArray)

 

For step 1, we have

subs(rho = rho(t), g_[])

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

(24)

Set this metric where `≡`(rho, rho(t))

"Setup(?):"

_______________________________________________________

 

`Coordinates: `[rho, phi, z]*`. Signature: `(`+ + +`)

 

_______________________________________________________

 

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

 

_______________________________________________________

(25)

Step 2, the geodesic equations in tensor notation with the coordinates depending on the time t are computed passing the optional argument tensornotation

Geodesics(t, tensornotation)

diff(diff((Physics:-SpaceTimeVector[`~a`](X))(t), t), t)+Physics:-Christoffel[`~a`, b, c]*(diff((Physics:-SpaceTimeVector[`~b`](X))(t), t))*(diff((Physics:-SpaceTimeVector[`~c`](X))(t), t)) = 0

(26)

Step 3: compute the components of this tensorial equation

TensorArray(diff(diff((Physics[SpaceTimeVector][`~a`](X))(t), t), t)+Physics[Christoffel][`~a`, b, c]*(diff((Physics[SpaceTimeVector][`~b`](X))(t), t))*(diff((Physics[SpaceTimeVector][`~c`](X))(t), t)) = 0, output = listofequations)

[diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2 = 0, diff(diff(phi(t), t), t)+2*(diff(rho(t), t))*(diff(phi(t), t))/rho(t) = 0, diff(diff(z(t), t), t) = 0]

(27)

These are the same equations (23).

 

Having the tensorial equation (26) is also useful to formulate the equations of motion in tensorial form in the presence of force. For that purpose, redefine the contravariant tensor F^j to represent the force in the cylindrical basis

Define(F[`~j`] = [`F__ρ`(t), `F__φ`(t), F__z(t)])

`Defined objects with tensor properties`

 

{Physics:-D_[a], Physics:-Dgamma[a], F[j], Physics:-Psigma[a], Physics:-Ricci[a, b], Physics:-Riemann[a, b, c, d], Physics:-Weyl[a, b, c, d], Physics:-d_[a], Physics:-g_[a, b], Physics:-Christoffel[a, b, c], Physics:-Einstein[a, b], Physics:-LeviCivita[a, b, c], Physics:-SpaceTimeVector[a](X)}

(28)

 

Newton's 2nd law (19)

F_(t) = m*(_rho(t)*(diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2)+_phi(t)*(2*(diff(rho(t), t))*(diff(phi(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k)

F_(t) = m*(_rho(t)*(diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2)+_phi(t)*(2*(diff(rho(t), t))*(diff(phi(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k)

(29)

now using tensorial notation, becomes

F[`~a`] = m*lhs(diff(diff((Physics[SpaceTimeVector][`~a`](X))(t), t), t)+Physics[Christoffel][`~a`, b, c]*(diff((Physics[SpaceTimeVector][`~b`](X))(t), t))*(diff((Physics[SpaceTimeVector][`~c`](X))(t), t)) = 0)

F[`~a`] = m*(diff(diff((Physics:-SpaceTimeVector[`~a`](X))(t), t), t)+Physics:-Christoffel[`~a`, b, c]*(diff((Physics:-SpaceTimeVector[`~b`](X))(t), t))*(diff((Physics:-SpaceTimeVector[`~c`](X))(t), t)))

(30)

TensorArray(F[`~a`] = m*(diff(diff((Physics[SpaceTimeVector][`~a`](X))(t), t), t)+Physics[Christoffel][`~a`, b, c]*(diff((Physics[SpaceTimeVector][`~b`](X))(t), t))*(diff((Physics[SpaceTimeVector][`~c`](X))(t), t))))

Array(%id = 18446744078329063774)

(31)

where we recall (see related Mapleprimes post) that to obtain the vector components entering `#mover(mi("F"),mo("→"))`(t) from these tensor components F[`~a`]we need to multiply the latter by the scale factors (1, rho(t), 1), the component of `#mover(mi("F"),mo("→"))`(t) in the direction of `#mover(mi("φ",fontstyle = "normal"),mo("∧"))` is given by rho(t)*m*(diff(phi(t), t, t)+2*(diff(rho(t), t))*(diff(phi(t), t))/rho(t)).

 

Coriolis force and centripetal force

 

After changing variables the position vector of the particle got expressed in (15) as

 

`#mover(mi("r"),mo("→"))`(t) = z(t)*`#mover(mi("k"),mo("∧"))`+`#mover(mi("ρ",fontstyle = "normal"),mo("∧"))`(t)*rho(t)

 

A distinction needs to be made here, according to whether the unit vector `#mover(mi("ρ",fontstyle = "normal"),mo("∧"))` depends or not on the time t, the former being the general case. When `#mover(mi("ρ",fontstyle = "normal"),mo("∧"))` is a constant, the value of the coordinate phi - the angle between `#mover(mi("ρ",fontstyle = "normal"),mo("∧"))` and the x axis - does not change, there is no rotation around the z axis. On the other hand, when `≡`(`#mover(mi("ρ",fontstyle = "normal"),mo("∧"))`, `#mover(mi("ρ",fontstyle = "normal"),mo("∧"))`(t)), the orientation of `#mover(mi("ρ",fontstyle = "normal"),mo("∧"))` and so the coordinate phi changes with time, so either the force `#mover(mi("F"),mo("→"))`(t)acting on the particle has a component in the `#mover(mi("φ",fontstyle = "normal"),mo("∧"))` direction that produces rotation around the z axis, or the system of references - itself - is rotating around the z axis.

 

Likewise, the expression (15)  can represent the position vector measured in the original Galilean (inertial) system of references, where a force `#mover(mi("F"),mo("→"))`(t)is producing rotation around the z axis, or it can represent the position of the particle measured in a rotating, non-inertial system references. Hence the transformation (14) can also be interpreted in two different ways, as representing a choice of different functions (generalized coordinates) to represent the position of the particle in the original inertial system of references, or it can represent a transformation from an inertial to another rotating, non-inertial, system of references.

 

This equivalence between the trajectory of a particle subject to an external force, as observed in an inertial system of references, and the same trajectory observed from a non-inertial accelerated system of references where there is no external force, actually at the root of the formulation of general relativity, is also well known in classical mechanics. The (apparent) forces observed in the rotating non-inertial system of references, due only to its acceleration, are called Coriolis and centripetal forces.

 

To see that the equations

 

diff(rho(t), t, t) = (diff(phi(t), t))^2*rho(t), diff(phi(t), t, t) = -2*(diff(phi(t), t))*(diff(rho(t), t))/rho(t)

 

that appeared in (27) when in the inertial system of references `#mover(mi("F"),mo("→"))`(t) = m*(diff(`#mover(mi("r"),mo("→"))`(t), t, t)) and m*(diff(`#mover(mi("r"),mo("→"))`(t), t, t)) = 0, are related to the Coriolis and centripetal forces in the non-inertial referencial, following [1] introduce a vector `#mover(mi("ω",fontstyle = "normal"),mo("→"))`representing the rotation of that referencial around the z axis (when, in the inertial system of references, the non-inertial system rotates clockwise, in the non-inertial system phi increases in value in the anti-clockwise direction)

`#mover(mi("ω",fontstyle = "normal"),mo("→"))` = -(diff(phi(t), t))*`#mover(mi("k"),mo("∧"))`

omega_ = -(diff(phi(t), t))*_k

(32)

According to [1], (39.7), the acceleration diff(`#mover(mi("r"),mo("→"))`(t), t, t)in the inertial system is expressed in terms of the quantities in the non-inertial rotating system by the sum of the following three vectorial terms.

First, the components of the acceleration `#mover(mi("a"),mo("→"))`(t)measured in the non-inertial system are given by the second derivatives of the coordinates (rho(t), phi(t), z(t)) multiplied by the scale factors, which in this case are given by (1, rho(t), 1) (see this post in Mapleprimes)

`#mover(mi("a"),mo("→"))`(t) = (diff(rho(t), t, t))*`#mover(mi("ρ",fontstyle = "normal"),mo("∧"))`(t)+rho(t)*(diff(phi(t), t, t))*`#mover(mi("φ",fontstyle = "normal"),mo("∧"))`(t)+(diff(z(t), t, t))*`#mover(mi("k"),mo("∧"))`

a_(t) = (diff(diff(rho(t), t), t))*_rho(t)+rho(t)*(diff(diff(phi(t), t), t))*_phi(t)+(diff(diff(z(t), t), t))*_k

(33)

Second, the Coriolis force divided by the mass, by definition given by

2*`&x`(diff(r_(t), t) = (diff(z(t), t))*_k+(diff(rho(t), t))*_rho(t)+rho(t)*(diff(phi(t), t))*_phi(t), omega_ = -(diff(phi(t), t))*_k)

2*Physics:-Vectors:-`&x`(diff(r_(t), t), omega_) = -2*rho(t)*(diff(phi(t), t))^2*_rho(t)+2*(diff(rho(t), t))*(diff(phi(t), t))*_phi(t)

(34)

Third the centripetal force divided by the mass, defined by

`&x`(omega_ = -(diff(phi(t), t))*_k, `&x`(r_(t) = z(t)*_k+rho(t)*_rho(t), omega_ = -(diff(phi(t), t))*_k))

Physics:-Vectors:-`&x`(omega_, Physics:-Vectors:-`&x`(r_(t), omega_)) = rho(t)*(diff(phi(t), t))^2*_rho(t)

(35)

Adding these three terms,

(a_(t) = (diff(diff(rho(t), t), t))*_rho(t)+rho(t)*(diff(diff(phi(t), t), t))*_phi(t)+(diff(diff(z(t), t), t))*_k)+(2*Physics[Vectors][`&x`](diff(r_(t), t), omega_) = -2*rho(t)*(diff(phi(t), t))^2*_rho(t)+2*(diff(rho(t), t))*(diff(phi(t), t))*_phi(t))+(Physics[Vectors][`&x`](omega_, Physics[Vectors][`&x`](r_(t), omega_)) = rho(t)*(diff(phi(t), t))^2*_rho(t))

a_(t)+2*Physics:-Vectors:-`&x`(diff(r_(t), t), omega_)+Physics:-Vectors:-`&x`(omega_, Physics:-Vectors:-`&x`(r_(t), omega_)) = _rho(t)*(diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2)+_phi(t)*(2*(diff(rho(t), t))*(diff(phi(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k

(36)

So that

diff(`#mover(mi("r"),mo("→"))`(t), t, t) = lhs(a_(t)+2*Physics[Vectors][`&x`](diff(r_(t), t), omega_)+Physics[Vectors][`&x`](omega_, Physics[Vectors][`&x`](r_(t), omega_)) = _rho(t)*(diff(diff(rho(t), t), t)-(diff(phi(t), t))^2*rho(t))+_phi(t)*(2*(diff(phi(t), t))*(diff(rho(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k)

diff(diff(r_(t), t), t) = a_(t)+2*Physics:-Vectors:-`&x`(diff(r_(t), t), omega_)+Physics:-Vectors:-`&x`(omega_, Physics:-Vectors:-`&x`(r_(t), omega_))

(37)

and where the right-hand side of (36) is, actually, the result computed lines above in (18)

diff(diff(r_(t), t), t) = _rho(t)*(diff(diff(rho(t), t), t)-(diff(phi(t), t))^2*rho(t))+_phi(t)*(2*(diff(phi(t), t))*(diff(rho(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k

diff(diff(r_(t), t), t) = _rho(t)*(diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2)+_phi(t)*(2*(diff(rho(t), t))*(diff(phi(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k

(38)

rhs(a_(t)+2*Physics[Vectors][`&x`](diff(r_(t), t), omega_)+Physics[Vectors][`&x`](omega_, Physics[Vectors][`&x`](r_(t), omega_)) = _rho(t)*(diff(diff(rho(t), t), t)-(diff(phi(t), t))^2*rho(t))+_phi(t)*(2*(diff(phi(t), t))*(diff(rho(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k)-rhs(diff(diff(r_(t), t), t) = _rho(t)*(diff(diff(rho(t), t), t)-(diff(phi(t), t))^2*rho(t))+_phi(t)*(2*(diff(phi(t), t))*(diff(rho(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k)

0

(39)

From (37), in the absence of external forces diff(`#mover(mi("r"),mo("→"))`(t), t, t) = 0 and so the acceleration `#mover(mi("a"),mo("→"))`(t) measured in the rotating system is given by the sum of the Coriolis and centripetal accelerations

isolate(rhs(diff(diff(r_(t), t), t) = a_(t)+2*Physics[Vectors][`&x`](diff(r_(t), t), omega_)+Physics[Vectors][`&x`](omega_, Physics[Vectors][`&x`](r_(t), omega_))), `#mover(mi("a"),mo("→"))`(t))

a_(t) = -2*Physics:-Vectors:-`&x`(diff(r_(t), t), omega_)-Physics:-Vectors:-`&x`(omega_, Physics:-Vectors:-`&x`(r_(t), omega_))

(40)

In other words: in the absence of external forces, the acceleration of a particle observed in a rotating (non-inertial) system of references is not zero.

 

Expressing this equation (40) in terms of (rho(t), phi(t), z(t)) we get

subs(a_(t) = (diff(diff(rho(t), t), t))*_rho(t)+rho(t)*(diff(diff(phi(t), t), t))*_phi(t)+(diff(diff(z(t), t), t))*_k, -(2*Physics[Vectors][`&x`](diff(r_(t), t), omega_) = -2*rho(t)*(diff(phi(t), t))^2*_rho(t)+2*(diff(rho(t), t))*(diff(phi(t), t))*_phi(t)), Physics[Vectors][`&x`](omega_, Physics[Vectors][`&x`](r_(t), omega_)) = rho(t)*(diff(phi(t), t))^2*_rho(t), a_(t) = -2*Physics[Vectors][`&x`](diff(r_(t), t), omega_)-Physics[Vectors][`&x`](omega_, Physics[Vectors][`&x`](r_(t), omega_)))

(diff(diff(rho(t), t), t))*_rho(t)+rho(t)*(diff(diff(phi(t), t), t))*_phi(t)+(diff(diff(z(t), t), t))*_k = rho(t)*(diff(phi(t), t))^2*_rho(t)-2*(diff(rho(t), t))*(diff(phi(t), t))*_phi(t)

(41)

resulting in the three equations

((diff(diff(rho(t), t), t))*_rho(t)+rho(t)*(diff(diff(phi(t), t), t))*_phi(t)+(diff(diff(z(t), t), t))*_k = rho(t)*(diff(phi(t), t))^2*_rho(t)-2*(diff(rho(t), t))*(diff(phi(t), t))*_phi(t)).`#mover(mi("ρ",fontstyle = "normal"),mo("∧"))`(t)

diff(diff(rho(t), t), t) = rho(t)*(diff(phi(t), t))^2

(42)

((diff(diff(rho(t), t), t))*_rho(t)+rho(t)*(diff(diff(phi(t), t), t))*_phi(t)+(diff(diff(z(t), t), t))*_k = rho(t)*(diff(phi(t), t))^2*_rho(t)-2*(diff(rho(t), t))*(diff(phi(t), t))*_phi(t)).`#mover(mi("φ",fontstyle = "normal"),mo("∧"))`(t)

rho(t)*(diff(diff(phi(t), t), t)) = -2*(diff(rho(t), t))*(diff(phi(t), t))

(43)

((diff(diff(rho(t), t), t))*_rho(t)+rho(t)*(diff(diff(phi(t), t), t))*_phi(t)+(diff(diff(z(t), t), t))*_k = rho(t)*(diff(phi(t), t))^2*_rho(t)-2*(diff(rho(t), t))*(diff(phi(t), t))*_phi(t)).`#mover(mi("k"),mo("∧"))`

diff(diff(z(t), t), t) = 0

(44)

which are the equations returned by Geodesics in (23)

[diff(diff(rho(t), t), t) = rho(t)*(diff(phi(t), t))^2, diff(diff(phi(t), t), t) = -2*(diff(rho(t), t))*(diff(phi(t), t))/rho(t), diff(diff(z(t), t), t) = 0]

[diff(diff(rho(t), t), t) = rho(t)*(diff(phi(t), t))^2, diff(diff(phi(t), t), t) = -2*(diff(rho(t), t))*(diff(phi(t), t))/rho(t), diff(diff(z(t), t), t) = 0]

(45)

``

References

[1] L.D. Landau, E.M. Lifchitz, Mechanics, Course of Theoretical Physics, Volume 1, third edition, Elsevier.


 

Download The_equations_of_motion_in_curvilinear_coordinates.mw

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

Since we are getting many questions on how to create Math apps to add to the Maple Cloud. I wanted to go over the different GUI aspects of how you go about creating a Math App in Maple. The following Document also includes some code examples that are used in the the Math App but doesn't go into them in detail. For more details on the type of coding you do in a Math App see the DocumentTools package help page.

Some of the graphical features of the Math app don't display on Maple Primes so I'd recommend downloading this worksheet from here: HowToMathApp.mw to follow along.


 

NULL

How to make a Math App (An example of using the Document Tools).

 

This Document will provide a beginners guide on one way to make a Math app in Maple.

It will contain some coding examples as well as where to find different options in the user interface.

Step 1 Insert a Table

 

 

• 

When making a Math App in Maple I often start with a table. You can enter a table by going to Insert > Table...

  

 

• 

I often make the table 1 x 2 to start with as this gives an area for input and an area for the output (such as plots).

NULL

 

Add a plot component to one of the cells of the table

 

 

• 

From the Components  Palette you can add a Plot Component . Add it to the table by clicking and dragging it over.

 

 

NULL

NULL

Add another table inside the other cell

 

 

• 

In the other cell of the table I'll add another table to organize my use of buttons, sliders, and other components.
NULL

NULL

Add some components to the new table

 

 

• 

From the Components Palette I'll add a slider, or dial, or something else for interaction.

 

• 

You may also want a Math region for an area to enter functions and a button to tell Maple to do something with it.

 

NULL

NULL

Arrange the Components to look nice

 

 

• 

You can change how the components are placed either by resizing the tables or changing the text orientation of the contents of the cells.

 

NULL

Write some code for the interaction of the buttons.

 

 

• 

Using the DocumentTools  package there are lots of ways you can use the components. I often will start writing my code using a code edit region  as it provides better visualization for syntax. On MaplePrimes these display as collapsed so I will also include code blocks for the code.

 

NULL

NULL

Let's write something that takes the value of the slider and applies it to the dial

 

 

• 

Note that the names of the components will change in each section as they are copies of the previous section.

 

with(DocumentTools):

14

with(DocumentTools):
sv:=GetProperty('Slider2',value);
SetProperty('Dial2',value,sv);
• 

This code will only execute when run using the  button. Change the value of the slider below then run the code above to see what happens.

 

NULL

NULL

Move the code 'inside' the slider

 

 

• 

Instead of putting the code inside the code edit region where it needs to be executed, we'll next add the code to the value changed code of the slider.

 

• 

Right click the Slider then select "Edit Value Changed Code".

 

 

• 

This will open the code editor for the Slider

 

 

• 

Enter your code (ensuring you're using the correct name for the slider and dial).

 

• 

Notice that you don't need to use the with(DocumentTools): command as "use DocumentTools in ... end use;" is already filled in for you.

 

• 

Save the code in the Slider and hit the  button inside it once.

• 

Now move the slider.

 

• 

On future uses of the App you won't need to hit  as the code will be run on startup.

``

NULL

NULL

Add some more details to your App

 

 

• 

Let's make this app do something a bit more interesting than change the contents of a dial when a slider moves.

 

• 

The plan in the next few steps is to make this app allow a user to explore parameters changing in a sinusoidal expression.

 

• 

I'm going to add a second Math Component, put the expression A*sin(t*theta+phi)into both then uncheck the box in the context panel that says "Editable".

 

• 

To make the Math containers fit nicely I'll check the Auto-fit container box and set the Minimum Width Pixels to 200.

 

``

Add code to change the value of phi in the second Math Container when the Slider changes

 

 

Note: Maple uses Radians for trigonometric functions so we should convert the value of phi to Radians.

use DocumentTools in

 

use DocumentTools in 
phi_s:=GetProperty(Slider5,value);
expr:= GetProperty(MathContainer6,expression);
new_expr:=algsubs(phi=phi_s*Pi/180,expr);

SetProperty(MathContainer7,expression,new_expr);
end use:

``

``

Make the Dial go from 0 to 360°

 

 

• 

Click the Dial and look at the options in the context panel on the right.

 

• 

Update the values in the Dial so that the highest position is 360 and the spacing makes sense for the app.

  NULL

``

Have the Dial update the theta value of the expression

 

 

• 

Add the following code to the Dial

 

use DocumentTools in
use DocumentTools in 
theta_d:=GetProperty(Dial7,value);
phi_s:=GetProperty(Slider7,value); #This is added so that phi also has the value updated

expr:= GetProperty(MathContainer10,expression);
new_expr0:=algsubs(theta=theta_d*Pi/180,expr);
new_expr:=algsubs(phi=phi_s*Pi/180,new_expr0);  #This is added so that phi also has the value updated

SetProperty(MathContainer11,expression,new_expr);
end use:

 

• 

Update the value in the slider to include the value from the dial

 

use DocumentTools in

 

use DocumentTools in 

theta_d:=GetProperty(Dial7,value); #This is added so that theta also has the value updated
phi_s:=GetProperty(Slider7,value); 

expr:= GetProperty(MathContainer10,expression);
new_expr0:=algsubs(theta=theta_d*Pi/180,expr); #This is added so that theta also has the value updated
new_expr:=algsubs(phi=phi_s*Pi/180,new_expr0);  

SetProperty(MathContainer11,expression,new_expr);

end use:

 

``

``

Notice that the code in the Dial and Slider are the same

 

 

• 

Since the code in the Dial and Slider are the same it makes sense to put the code into a procedure that can be called from multiple places.

 

Note: The changes in the code such as local and the single quotes are not needed but make the code easier to read and less likely to run into errors if edited in the future (for example if you create a variable called dial8 it won't interfere now that the names are in quotes).

 

 

UpdateMath:=proc() 

UpdateMath:=proc()
local theta_d, phi_s, expr, new_expr, new_expr0;
use DocumentTools in 
theta_d:=GetProperty('Dial8','value'); #Get value of theta from Dial
phi_s:=GetProperty('Slider8','value'); #Get value of phi from slider

expr:= GetProperty('MathContainer12','expression');
new_expr0:=algsubs('theta'=theta_d*Pi/180,expr);  # Put value of theta in expression
new_expr:=algsubs('phi'=phi_s*Pi/180,new_expr0);  # Put value of phi in expression
SetProperty('MathContainer13','expression',new_expr); # Update expression
end use:
end proc:

 

• 

Now change the code in the components to call the function using UpdateMath().

 

• 

Since the code above is only defined there it will need to be run once (but only once) before moving the components. Instead of leaving it here you can add it to the Startup code by clicking  or going to Edit > Startup code.  This code will run every time you open the Math App ensuring that it works right away.

 

• 

The startup code isn't defined in this document to allow progression of these steps.

 

``

Make the button initialize the app

 

 

• 

Since the startup code isn't defined in this document we are going to move this function into the button.

 

UpdateMath:=proc()

 

UpdateMath:=proc()
local theta_d, phi_s, expr, new_expr, new_expr0;
use DocumentTools in 
theta_d:=GetProperty('Dial9','value'); #Get value of theta from Dial
phi_s:=GetProperty('Slider9','value'); #Get value of phi from slider

expr:= GetProperty('MathContainer14','expression');
new_expr0:=algsubs('theta'=theta_d*Pi/180,expr);  # Put value of theta in expression
new_expr:=algsubs('phi'=phi_s*Pi/180,new_expr0);  # Put value of phi in expression
SetProperty('MathContainer15','expression',new_expr); # Update expression
end use:
end proc:
• 

First click the button to rename it, you'll see the  option in the context panel on the right. Then add the code above to the button in the same way as the Slider an Dial (Right click and select Edit Click Code).

 

``

``

Now it is easy to add new components

 

 

• 

Now if we want to add new components we just have to change the one procedure.  Let's add a Volume Gauge to change the value of A.

 

• 

Click in the cell containing the Dial, the context panel will show the option to Insert a row below the Dial.

• 

Now drag a Volume Gauge into the new cell.

 

• 

Click in the cell and choose the alignment (from the context panel) that looks best to you. In this case I chose center:

 

``

 

NULL

``

Update the procedure code for the Gauge

 

 

• 

Add two lines for the volume gauge to get the value and sub it into the expression

UpdateMath:=proc()

UpdateMath:=proc()
local theta_d, phi_s, expr, new_expr, new_expr0;
use DocumentTools in 
theta_d:=GetProperty('Dial11','value'); #Get value of theta from the Dial
phi_s:=GetProperty('Slider11','value'); #Get value of phi from the Slider
A_g:=GetProperty('VolumeGauge1','value'); #Get value of A from the Guage

expr:= GetProperty('MathContainer18','expression');
new_expr0:=algsubs('theta'=theta_d*Pi/180,expr);  # Put value of theta in expression
new_expr1:=algsubs('phi'=phi_s*Pi/180,new_expr0);  # Put value of phi in expression
new_expr:=algsubs('A'=A_g,new_expr1);  # Put value of A in expression

SetProperty('MathContainer19','expression',new_expr); # Update expression
end use:
end proc:
• 

Now add

UpdateMath();

  to the Gauge.

  ``

``

Plot the changing expression

 

 

• 

Make a procedure to get the value in the second Math Container and plot it

 

PlotMath:=proc()

PlotMath:=proc()
	local expr, p;
	use DocumentTools in 

	expr:=GetProperty('MathContainer21','expression'); 

	p:=plot(expr,'t'=-Pi/2..Pi/2,'view'=[-Pi/2..Pi/2,-100..100]):

	SetProperty('Plot14','value',p)
	end use:
end proc:
• 

Put this procedure in the Initialize button and the call to it in the components.

 

NULL

``

Tidy up the app

 

 

• 

Now that we have an interactive app let's tidy it up a bit.

 

• 

The first thing I'd recommend in your own app is moving the code from the initialize button to startup code. In this document we choose to use the button instead to preserve earlier versions.

 

• 

You can also remove the borders around the components by clicking in the table and selecting "Interior Borders" > "None" and "Exterior Borders" > "None" from the context panel.

NULL

``

``

Now you have a Math App

 

 

• 

You can upload your Math App to the Maple Cloud to share with others by going to "File" > "Save to Cloud".

 

• 

I'd recommend also including a description of what your app does. You can do this nicely using another table and Text mode.

 

 

 

``

``

NULL

HowToMathApp.mw

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.

First 10 11 12 13 14 15 16 Last Page 12 of 66