Applications, Examples and Libraries

Simulation of the DNA with Maple

by: Lenin Araujo Castillo 1790 Product: Maple 2020

May 01 2020

6 0

This app shows the modeling and simulation of DNA carried out entirely in Maple. The mathematical model is inserted through the combination of trigonometric functions. It shows the graphs of the curvature vs time for its interpretation. Made for engineering and health science students.

MV_AC_R3_UNI_2020.mw

Lenin Araujo Castillo

Ambassador of Maple

Maple for Beginners

by: ecterrab 14440 Product: Maple

April 15 2020

11 0

Maple_for_Beginners.mw

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

Vectors in Spherical Coordinates using Tensor...

by: ecterrab 14440 Product: Maple

April 07 2020

6 3

Vectors in Spherical Coordinates using Tensor Notation

Edgardo S. Cheb-Terrab¹ and Pascal Szriftgiser²

(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

>

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

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

(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 , which in Cartesian coordinates has the simple form

>

(1.1)

To compute the line element in spherical coordinates, the starting point is the transformation

>

(1.2)

>

(1.3)

Since in are just symbols with no relationship to start transforming these differentials using the chain rule, computing the Jacobian of the transformation (1.2). In this Jacobian J, the first line is , ,

>

So in matrix notation,

>

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

(1.4)

To complete the computation of in spherical coordinates we can now use ChangeBasis , provided that next we substitute (1.4) in the result, expressing the abstract objects in terms of .

In two steps:

>

(1.5)

The line element

>

(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 , then taking the coefficients of the differentials (just substitute them by the number 1)

>

(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 , we see that corresponding vector, the line element in spherical coordinates , is not constructed by directly equating its components to the components of , so

The vector is constructed multiplying these contravariant components by the scaling factors, as

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 the same way we did in the line above with the line element, using the scale-factors , so that

$`#mover(mi("A"),mo("→"))` = 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 and the scale-factors . Because the system is orthogonal, each vector component satisfies

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

it follows that,

A__j = A__j/h__j

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

•	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: will represent A in Cartesian coordinates, while will represent A in spherical coordinates.

>

${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("μ",fontstyle = "normal")),mo("⁡"),mfenced(mi("y")),mo("="),mfrac(mrow(mo("&PartialD;"),msup(mi("y"),mi("μ",fontstyle = "normal"))),mrow(mo("&PartialD;"),msup(mi("x"),mi("ν",fontstyle = "normal"))),linethickness = "1"),mo("⁢"),mo("⁢"),msup(mi("A"),mi("ν",fontstyle = "normal")),mfenced(mi("x")))`$ , that is the same as the rule for the differential of the coordinates. Then, the transformation rule from to 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

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 , from spherical (4^th argument) to Cartesian (3^rd argument), which according to the 5^th bullet of TransformCoordinates will result in a transformation expressed in terms of the old coordinates (here the spherical ). Expand things to make the comparison with (1.4) possible by eye

>

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

(2.2)

We see that the transformation rule for a contravariant vector 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 , we know, by definition, that it is $`#mrow(msub(mi("A"),mi("μ",fontstyle = "normal")),mo("⁡"),mfenced(mi("y")),mo("="),mfrac(mrow(mo("&PartialD;"),msup(mi("x"),mi("ν",fontstyle = "normal"))),mrow(mo("&PartialD;"),msup(mi("y"),mi("μ",fontstyle = "normal"))),linethickness = "1"),mo("⁢"),mo("⁢"),msub(mi("A"),mi("ν",fontstyle = "normal")),mfenced(mi("x")))`$ , so the same transformation rule for the gradient , where and so on. We can experiment this by directly changing variables in the differential operators , for example

>

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 )

>

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

(2.4)

The corresponding transformation equations relating the tensors and 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](%id = 18446744078557373854) = Vector[column](%id = 18446744078557374334)

(2.5)

We see that the transformation rule for a covariant vector 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](%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:

>

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

>

(3.1)

>

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

>

%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 are the three lines in the right-hand side of (2.6) after replacing the covariant components by the derivatives of f with respect to the j^th coordinate, here displayed using indexed notation due to using CompactDisplay

>

(3.4)

>

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

>

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

(3.6)

>

(3.7)

Simplifying, we arrive at (3.2)

>

(3.8)

>

(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

>

(4.1)

>

(4.2)

>

%Divergence(%A__s_) = ((diff(A__r(S), r))*r+2*A__r(S))/r+((diff(`A__θ`(S), theta))*sin(theta)+`A__θ`(S)*cos(theta))/(r*sin(theta))+(diff(`A__φ`(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 . For that purpose, we first transform the coordinates and the metric introducing nonzero Christoffel symbols

>

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

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

(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

>

(4.5)

Consider now the divergence of the contravariant tensor, computed in tensor notation

>

(4.6)

>

(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

>

(4.8)

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

>

%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("→"))`$ in (4.3) ? The answer is in the relationship established at the end of Sec I between the components of the tensor and the components of the vector $`#mover(mi("\`A__s\`"),mo("→"))`$ , namely that the vector components are obtained multiplying the contravariant tensor components by the scale-factors . So, in the above we need to substitute the contravariant by the vector components divided by the scale-factors

>

(4.10)

>

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

(4.11)

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

>

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

(4.12)

•

The Curl

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

>

%Curl(%A__s_) = ((diff(`A__φ`(S), theta))*sin(theta)+`A__φ`(S)*cos(theta)-(diff(`A__θ`(S), phi)))*_r/(r*sin(theta))+(diff(A__r(S), phi)-(diff(`A__φ`(S), r))*r*sin(theta)-`A__φ`(S)*sin(theta))*_theta/(r*sin(theta))+((diff(`A__θ`(S), r))*r+`A__θ`(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

But in a curvilinear system is not a tensor, we need to use the non-Galilean form , where is the determinant of the metric. Moreover, since the expression 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

>

(4.14)

>

%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 by the vector components 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

>

(4.16)

Replace now the components of the tensor by the components of the 3D vector $`#mover(mi("\`A__s\`"),mo("→"))`$ using (4.10)

>

(4.17)

>

%Curl(%A__s_) = [((diff(`A__φ`(S), theta))*sin(theta)+`A__φ`(S)*cos(theta)-(diff(`A__θ`(S), phi)))/(r*sin(theta)), (diff(A__r(S), phi)-(diff(`A__φ`(S), r))*r*sin(theta)-`A__φ`(S)*sin(theta))/(r*sin(theta)), ((diff(`A__θ`(S), r))*r+`A__θ`(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__φ`(S), theta))*sin(theta)+`A__φ`(S)*cos(theta)-(diff(`A__θ`(S), phi)))*_r/(r*sin(theta))+(diff(A__r(S), phi)-(diff(`A__φ`(S), r))*r*sin(theta)-`A__φ`(S)*sin(theta))*_theta/(r*sin(theta))+((diff(`A__θ`(S), r))*r+`A__θ`(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 is a scalar, the covariant derivative becomes the standard derivative .

>

(4.21)

The above is the expected result (3.2)

>

(4.22)

•	The Laplacian

Likewise we can compute the Laplacian directly as

>

(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,

>

%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,

>

((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 , in curvilinear coordinates; rewrite it using standard diff derivatives and expand both sides of the equation for comparison

>

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:

>

(4.27)

Download Vectors_and_Spherical_coordinates_in_tensor_notation.mw

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

Нoles on an implicit plot

by: one man 1335 Product: Maple 17

April 01 2020

2 0

A way of cutting holes on an implicit plot. This is from the field of numerical parameterization of surfaces. On the example of the surface x3 = 0.01*exp (x1) / (0.01 + x1^4 + x2^4 + x3^4) consider the approach to producing holes. The surface is locally parameterized in some suitable way and the place for the hole and its size are selected. In the first example, the parametrization is performed on the basis of the section of the initial surface by perpendicular planes. In the second example, "round" parametrization. It is made on the basis of the cylinder and the planes passing through its axis. Holes can be of any size and any shape. In the figures, the cut out surface sections are colored green and are located above their own holes at an equidistant to the original surface.
HOLE_1.mw HOLE_2.mw

How high is your chance to die for CoViD19 once...

by: mmcdara 7691 Product: Maple 2015

March 29 2020

2 4

Hi,

The present work is aimed to show how bayesian inference methods can be used to infer (= to assess) the probabilility that a person detected infected by the SARS-Cov2 has to die (remark I did not write "has to die if it" because one never be sure of the reason of the death).
A lot of details are avaliable in the attached pdf file (I tried to be pedagogic enough so that the people not familiar with bayesian inference can get a global understanding of the subject, many links are provided for quick access to the different notions).

In particular, I explain why simple mathematics cannot provide a reliable estimate of this probability of death (sometimes referred to as the "death rate") as long as the epidemic continues to spread.

Even if the approach presented here is rather original, this is not the purpose of this post.
Since a long time I had in mind to post here an application concerning bayesian methods. The CoVid19 outbreak has only provided me with the most high-profile topic to do so.
I will say no more about the inference procedure itself (all the material is given in the attached pdf file) and I will only concentrate on the MAPLE implementation of the solution algorithm.

Bayesian Inference uses generally simple algorithms such as MCMC (Markov Chain Monte Carlo) or ABC (Approximate Bayesian Computation) to mention a few, and their corresponding pseudo code writes generally upon a few tens of lines.
This is something I already done with other languages but I found the task comparatively more difficult with Maple. Probably I was to obsess not to code in Maple as you code in Matlab or R for instance.
At the very end the code I wrote is rather slow, this because of the allocated memory size it uses.
In a question I posed weeks ago (How can I prevent the creation of random variables...) Preben gave a solution to limit the burst of the memory: the trick works well but I'm still stuked with memory size problems (Acer also poposed a solution but I wasn't capable to make it works... maybe I was too lazzy to modify deeply my code).

Anyway, the code is there, in case anyone would like to take up the challenge to make it more efficient (in which case I'll take it).

Note 1: this code contains a small "Maplet" to help you choose any country in the data file on which you would like to run the inference.
Note 2: Be careful: doing statistics, even bayesian statistics, needs enough data: some countries have history records ranging over a few days , or no recorded death at all; infering something from so loos date will probably be disappointing

The attached files:

The pdf file is the "companion document" where all or most of it is explained.It has been written a few days ago for another purpose and the results it presents were not ontained from the lattest data (march 21, 2020 coronavirus)
xls files are data files, they were loaded yesterday (march 28, 2020) from here coronavirus
the mw file... well, I guess you know what it is.

Bayesian_inference.pdf

total-cases-covid-19_NF.xls

total-deaths-covid-19_NF.xls

Bayesian_Inference_ABC+MCMC_NF_2.mw

Puzzle - cut it into 2 equal parts

by: Kitonum 21370

March 23 2020

5 0

The following puzzle prompted me to write this post: "A figure is drawn on checkered paper that needs to be cut into 2 equal parts (the cuts must pass along the sides of the squares.)" (parts are called equal if, after cutting, they can be superimposed on one another, that is, if one of them can be moved, rotated and (if need to) flip so that they completely coincide) (see the first picture below).
I could not solve it manually and wrote a procedure called CutTwoParts that does this automatically (of course, this procedure applies to other similar puzzles). This procedure uses my procedure AreIsometric published earlier https://www.mapleprimes.com/posts/200157-Testing-Of-Two-Plane-Sets-For-Isometry (for convenience, I have included its text here). In the procedure CutTwoParts the figure is specified by the coordinates of the centers of the squares of which it consists).

I advise everyone to first try to solve this puzzle manually in order to feel its non-triviality, and only then load the worksheet with the procedure for automatic solution.

For some reason, the worksheet did not load and I was only able to insert the link.

Cuttings.mw

First equilibrium condition with Maple

by: Lenin Araujo Castillo 1790 Product: Maple 2020

March 20 2020

2 0

With this application our students of science and engineering in the areas of physics will check the first condition of balance using Maple technology. Only with entering mass and angles we obtain graphs and data for a better interpretation.

First_equilibrium_condition.zip

Lenin AC

Ambassador of Maple

Add Reverb to Audio

by: Samir Khan 2091 Product: Maple 2020

March 18 2020

4 0

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

Exploring the CoVid19 outbreak

by: mmcdara 7691 Product: Maple

March 14 2020

6 14

Hi,

Two weeks ago, I started loading data on the CoVid19 outbreak in order to understand, out of any official communication from any country, what is really going on.

From february 29 to march 9 these data come from https://bnonews.com/index.php/2020/02/the-latest-coronavirus-cases/ and from 10 march until now from https://www.worldometers.info/coronavirus/#repro.In all cases the loading is done manually (copy-paste onto a LibreOffice spreadsheet plus correction and save into a xls file)Â for I wasn't capable to find csv data (csv data do exist here https://github.com/CSSEGISandData/COVID-19, by they end febreuary 15th).
So I copied-pasted the results from the two sources above into a LibreOffice spreadsheet, adjusted the names of some countries for they appeared differently (for instance "United States" instead of "USA"), removed the unnessary commas and saved the result in a xls file.

I also used data from https://www.worldometers.info/world-population/population-by-country/ to get the populations of more than 260 countries around the world and, finally, csv data from https://ourworldindata.org/coronavirus#covid-19-tests to get synthetic histories of confirmed and death cases (I have discovered this site only yesterday evening and I think it could replace all the data I initially loaded).

The two worksheet here are aimed to exploratory and visualization only.
An other one is in progress whose goal is to infer the true death rate (also known as CFR, Case Fatality Rate).

No analysis is presented, if for no other reason than that the available data (except the numbers of deaths) are extremely dependent on the testing policies in place. But some features can be drawn from the data used here.
For instance, if you select country = "China" in file Covid19_Evolution_bis.mw, you will observe very well known behaviour which is that the "Apparent Death Rate", I defined as the ratio of the cumulated number of death at time t by the cumulatibe number of confirmed cases at the same time, is always an underestimation of the death rate one can only known once the outbreak has ended. With this in mind, changing the country in this worksheet from China to Italy seems to lead to frightening scary interpolations... But here again, without knowing the test policy no solid conclusion can be drawn: maybe Italy tests mainly elder people with accute symptoms, thus the huge "Apparent Death Rate" Italy seems to have?

The work has been done with Maple 2015 and some graphics can be improved if a newer version is used (for instance, as Maple 2015 doesn't allow to change the direction of tickmarks, I overcome this limitation by assigning the date to the vertical axis on some plots).
The second Explore plot could probably be improved by using newer versions or Maplets or Embeded components.

Explore data from https://bnonews.com/index.php/2020/02/the-latest-coronavirus-cases/ and https://www.worldometers.info/coronavirus/#repro
Files to use
Covid19_Evolution.mw
Covid19_Data.m.zip
Population.xls

Explore data from https://ourworldindata.org/coronavirus#covid-19-tests
Files to use
Covid19_Evolution_bis.mw
daily-deaths-covid-19-who.xls
total-cases-covid-19-who.xls
Population.xls

I would be interested by any open collaboration with people interested by this post (it's not in my intention to write papers on the subject, my only motivation is scientific curiosity).

Optimizing Mini Golf in MapleSim

by: TechnicalSupport 750 Products: Global Optimization , MapleSim 2019

February 12 2020

13 0

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

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

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

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

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

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

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

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

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

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

Mini_golf_Global_Optimization.msim

A riddle !

by: vv 13625 Product: Maple 2018

February 11 2020

7 3

Can you guess what P() produces, without executing it?

P:=proc(N:=infinity) local q,r,t,k,n,l,h, f;
q,r,t,k,n,l,h := 1,0,1,1,3,3,0:
while h<N do 
   if 4*q+r-t < n*t
   then f:=`if`(++h mod 50=0,"\n",`if`(h mod 10=0," ","")); printf("%d"||f,n);   
        q,r,t,k,n,l := 10*q,10*(r-n*t),t,k,iquo(10*(3*q+r),t)-10*n,l
   else q,r,t,k,n,l := q*k,(2*q+r)*l,t*l,k+1,iquo(q*(7*k+2)+r*l,t*l),l+2
   fi
od: NULL
end:

I hope you will like it (maybe after execution).

Feynman Diagrams - the scattering matrix in coordi...

by: ecterrab 14440 Product: Maple

January 25 2020

7 0

Feynman Diagrams
The scattering matrix in coordinates and momentum representation

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

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

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

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

Introduction

A scattering matrix relates the initial and final states, and , of an interacting system. In an 4-dimensional spacetime with coordinates , can be written as:

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

This exponential can be expanded as

where

$S[n] = `#mrow(mo("⁡"),mfrac(msup(mi("i"),mi("n")),mrow(mi("n"),mo("&excl;")),linethickness = "1"),mo("⁢"),mo("∫"),mi("…"),mo("⁢"),mo("∫"),mi("T"),mo("⁡"),mfenced(mrow(mi("L"),mo("⁡"),mfenced(mi("\`X__1\`")),mo(","),mi("…"),mo(","),mi("L"),mo("⁡"),mfenced(mi("\`X__n\`")))),mo("⁢"),mo("&DifferentialD;"),msup(mi("\`X__1\`"),mn("4")),mo("⁢"),mi("…"),mo("⁢"),mo("&DifferentialD;"),msup(mi("\`X__n\`"),mn("4")))`$

and is the time-ordered product of n interaction Lagrangians evaluated at different points. The S matrix formulation is at the core of perturbative approaches in relativistic Quantum Field Theory.

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

a.	The expansion in coordinates representation up to arbitrary order (the limitation is now only your hardware)

b.

The S-matrix element in momentum representation up to arbitrary order for given number of loops and initial and final particles (the contents of the and states); optionally, also the transition probability density, constructed using the square of the scattering matrix element , as shown in formula (13) of sec. 21.1 of ref.[1].

c.	The Feynman diagrams (drawings) related to the different terms of the expansion of S or of its matrix elements .

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

Examples

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

>

(1.1)

Let be the interaction Lagrangian

>

(1.2)

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

>

(1.3)

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

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

>

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

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

(1.4)

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

A scattering process with matrix element in momentum representation, corresponding to the term with 4 external legs (symmetry factor = 72), could be any process where the total number of incoming + outgoing parties is equal to 4. For example, one with 2 incoming and 2 outgoing particles. The transition probability for that process is given by

>

(1.5)

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

The transition matrix element is related to the transition probability density (formula (13) of sec. 21.1 of ref.[1]) by

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

where represent the particle densities of each of the s particles in the initial state , the (Dirac) is the expected singular factor due to the conservation of the energy-momentum and the amplitude is related to via

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

To directly get the probability density instead ofuse the option output = probabilitydensity

>

(1.6)

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

>

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

(1.7)

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

>

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

(1.8)

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

>

(1.9)

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

>

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

(1.10)

A corresponding S-matrix element in momentum representation:

>

(1.11)

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

>

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

(1.12)

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

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

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

(1.13)

The matrix indices of the Dirac matrices are written explicitly and use conjugate to represent the Dirac conjugate

>

(1.14)

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

>

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

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

(1.15)

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

>

(1.16)

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

>

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

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

(1.17)

The corresponding S-matrix elements in momentum representation

>

(1.18)

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

>

(1.19)

Now the propagator for is the one of a massless vector field:

>

(1.20)

The self-energy of the photon:

>

(1.21)

where is the corresponding polarization vector.

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

>

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

(1.22)

>

(1.23)

>

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

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

(1.24)

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

>

(1.25)

(1.26)

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

>

(1.27)

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

>

(1.28)

>

(1.29)

>

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

(1.30)

>

(1.31)

>

(1.32)

>

(1.33)

>

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

(1.34)

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

>

References

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

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

FeynmanDiagrams_and_the_Scattering_Matrix.PDF

FeynmanDiagrams_and_the_Scattering_Matrix.mw

FeynmanDiagrams_-_help_page.mw

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

Happy New Year!

by: vv 13625 Product: Maple

December 31 2019

8 6

with(plots):
S:=cat("Happy New Year 2020!   "$3):
N:=length(S): a:=0.77*Pi: h:=2*Pi/N:
display(seq(textplot([cos(a-k*h), sin(a-k*h),S[k+1]], 
        rotation=-Pi/2+a-k*h, 'font'=["times","roman",24]), k=0..N-4), axes=none);

A Parametric Happy New Year

by: Samir Khan 2091 Product: Maple

December 31 2019

2 4

When plotted, these parametric equations say "happy new year" (and were constructed with this worksheet)

x := piecewise(t <= 58, -15.0*sin(1.43 + 0.650*t) - 14.8*sin(-1.64 + 0.703*t) - 1.28*sin(-2.97 + 1.25*t) - 11.9*sin(-1.58 + 0.540*t) - 1.07*sin(-1.60 + 1.35*t) - 3.85*sin(-2.09 + 1.41*t) - 7.13*sin(1.13 + 1.73*t) - 4.40*sin(1.32 + 1.30*t) - 26.3*sin(1.53 + 0.380*t) - 9.42*sin(-4.65 + 0.433*t) - 3.43*sin(1.42 + 2.06*t) - 7.57*sin(-1.77 + 2.11*t) - 2.65*sin(-4.34 + 0.323*t) - 1.95*sin(-4.57 + 2.54*t) - 5.39*sin(-1.38 + 2.60*t) - 49.2*sin(1.52 + 0.487*t) - 0.754*sin(-4.38 + 2.87*t) - 9.67*sin(-1.58 + 2.65*t) - 7.88*sin(-4.59 + 1.95*t) - 2.39*sin(-1.67 + 2.71*t) - 15.1*sin(1.53 + 0.108*t) - 39.0*sin(1.47 + 0.757*t) - 1.80*sin(1.37 + 2.22*t) - 4.22*sin(-1.95 + 0.973*t) - 7.72*sin(-1.44 + 2.17*t) - 8.80*sin(-1.66 + 0.813*t) - 3.59*sin(1.13 + 1.57*t) - 15.4*sin(-1.64 + 1.62*t) - 6.70*sin(1.36 + 1.19*t) - 791.*sin(-1.57 + 0.0540*t) - 2.55*sin(-1.55 + 1.89*t) - 6.92*sin(-1.87 + 1.68*t) - 3.95*sin(1.17 + 1.08*t) - 44.1*sin(-1.67 + 1.14*t) - 25.8*sin(1.51 + 0.597*t) - 31.4*sin(1.42 + 1.46*t) - 96.8*sin(-1.59 + 0.162*t) - 18.7*sin(1.53 + 0.217*t) - 7.87*sin(-4.66 + 2.98*t) - 4.99*sin(1.22 + 3.03*t) - 6.92*sin(1.43 + 2.44*t) - 48.3*sin(1.47 + 1.03*t) - 24.2*sin(1.48 + 1.52*t) - 9.58*sin(1.43 + 2.49*t) - 4.29*sin(1.33 + 2.27*t) - 6.34*sin(1.22 + 2.33*t) - 12.0*sin(1.45 + 2.00*t) - 0.388*sin(-1.25 + 2.92*t) - 2.74*sin(-1.43 + 1.79*t) - 6.71*sin(-1.66 + 1.84*t) - 0.713*sin(-3.63 + 2.38*t) - 43.1*sin(-1.59 + 0.271*t) - 2.51*sin(1.12 + 2.76*t) - 1.29*sin(-3.92 + 2.82*t) - 21.3*sin(-1.70 + 0.867*t) - 12.4*sin(1.50 + 0.920*t), 58 < t and t <= 84, -500 - 321.*sin(-8.608 + 0.121*t) - 7.18*sin(-12.408 + 0.241*t) - 57.1*sin(-22.608 + 0.361*t) - 21.9*sin(-26.682 + 0.484*t) - 21.3*sin(-33.474 + 0.603*t) - 50.2*sin(-43.800 + 0.725*t) - 20.6*sin(-50.760 + 0.845*t) - 41.5*sin(-54.756 + 0.967*t) - 9.74*sin(-61.89 + 1.09*t) - 41.1*sin(-72.03 + 1.21*t) - 2.49*sin(-78.88 + 1.33*t) - 3.30*sin(-83.227 + 1.45*t) - 6.73*sin(-89.99 + 1.57*t) - 5.88*sin(-96.59 + 1.69*t) - 16.4*sin(-106.99 + 1.81*t) - 1.61*sin(-111.8982 + 1.93*t) - 1.84*sin(-117.970 + 2.05*t) - 0.464*sin(-127.83 + 2.17*t) - 1.64*sin(-134.90 + 2.30*t) - 3.94*sin(-142.37 + 2.41*t) - 2.35*sin(-149.22 + 2.54*t) - 2.72*sin(-154.3362 + 2.66*t) - 8.41*sin(-160.453 + 2.78*t) - 4.39*sin(-171.17 + 2.90*t), 84 < t, -300 - 2.66*sin(-205.04 + 2.41*t) - 1.26*sin(-207.397 + 2.46*t) - 2.21*sin(-196.59 + 2.31*t) - 2.31*sin(-166.83 + 1.96*t) - 48.9*sin(-39.688 + 0.452*t) - 0.697*sin(-252.158 + 3.01*t) - 2.51*sin(-179.22 + 2.11*t) - 1.57*sin(-222.14 + 2.66*t) - 0.745*sin(-226.24 + 2.71*t) - 49.4*sin(-10.020 + 0.100*t) - 0.289*sin(-159.628 + 1.91*t) - 95.9*sin(-32.358 + 0.402*t) - 60.0*sin(-43.928 + 0.502*t) - 3.76*sin(-73.736 + 0.854*t) - 3.05*sin(-183.97 + 2.16*t) - 0.629*sin(-158.50 + 1.86*t) - 9.25*sin(-49.272 + 0.603*t) - 4.46*sin(-74.716 + 0.904*t) - 10.4*sin(-79.040 + 0.955*t) - 2.65*sin(-103.67 + 1.21*t) - 1.99*sin(-145.57 + 1.71*t) - 1.52*sin(-197.315 + 2.36*t) - 0.685*sin(-258.12 + 3.06*t) - 1.04*sin(-247.58 + 2.96*t) - 64.8*sin(-18.514 + 0.201*t) - 68.5*sin(-31.278 + 0.352*t) - 579.*sin(-5.8068 + 0.0502*t) - 6.52*sin(-95.20 + 1.11*t) - 5.03*sin(-96.28 + 1.16*t) - 0.396*sin(-211.620 + 2.51*t) - 7.28*sin(-150.00 + 1.76*t) - 2.42*sin(-153.92 + 1.81*t) - 10.4*sin(-112.11 + 1.31*t) - 24.8*sin(-85.95 + 1.00*t) - 3.91*sin(-124.83 + 1.46*t) - 1.69*sin(-185.369 + 2.21*t) - 1.18*sin(-189.238 + 2.26*t) - 16.6*sin(-56.662 + 0.653*t) - 1.33*sin(-222.31 + 2.61*t) - 0.593*sin(-238.70 + 2.81*t) - 1.88*sin(-233.58 + 2.76*t) - 3.91*sin(-133.01 + 1.56*t) - 4.94*sin(-134.16 + 1.61*t) - 9.59*sin(-128.89 + 1.51*t) - 1.02*sin(-240.2714 + 2.86*t) - 2.15*sin(-247.83 + 2.91*t) - 5.52*sin(-90.85 + 1.06*t) - 3.83*sin(-171.25 + 2.01*t) - 0.523*sin(-171.66 + 2.06*t) - 0.284*sin(-141.80 + 1.66*t) - 23.2*sin(-11.174 + 0.151*t) - 1.58*sin(-114.615 + 1.36*t) - 2.67*sin(-120.75 + 1.41*t) - 5.83*sin(-19.524 + 0.251*t) - 13.7*sin(-23.774 + 0.301*t) - 14.8*sin(-107.89 + 1.26*t) - 15.5*sin(-60.842 + 0.703*t) - 37.7*sin(-65.176 + 0.754*t) - 2.02*sin(-217.95 + 2.56*t) - 13.2*sin(-69.466 + 0.804*t) - 37.7*sin(-45.052 + 0.553*t)):

y := piecewise(t <= 58, -28.1*sin(1.45 + 1.62*t) - 2.23*sin(-2.39 + 1.89*t) - 17.8*sin(-1.51 + 1.19*t) - 4.85*sin(-1.61 + 2.38*t) - 2.52*sin(1.55 + 1.95*t) - 20.0*sin(1.55 + 2.11*t) - 24.8*sin(-1.62 + 2.00*t) - 19.9*sin(-1.81 + 2.06*t) - 4.22*sin(-0.422 + 2.60*t) - 6.94*sin(1.47 + 2.87*t) - 61.1*sin(1.49 + 0.323*t) - 13.9*sin(-4.68 + 0.540*t) - 3.97*sin(0.00256 + 2.33*t) - 69.8*sin(1.53 + 0.487*t) - 59.6*sin(1.50 + 0.813*t) - 132.*sin(-1.65 + 0.867*t) - 26.7*sin(-1.76 + 1.52*t) - 53.1*sin(1.40 + 1.57*t) - 139.*sin(1.57 + 0.0540*t) - 3.75*sin(-2.34 + 3.03*t) - 8.03*sin(1.24 + 1.73*t) - 22.9*sin(-4.61 + 0.217*t) - 16.7*sin(-1.67 + 0.703*t) - 23.3*sin(-1.82 + 1.68*t) - 78.9*sin(-4.70 + 0.271*t) - 2.72*sin(-2.38 + 2.49*t) - 3.45*sin(1.10 + 2.54*t) - 2.07*sin(-0.489 + 2.22*t) - 13.1*sin(-1.82 + 2.27*t) - 60.6*sin(-1.62 + 1.08*t) - 5.27*sin(1.55 + 2.44*t) - 4.17*sin(1.46 + 2.82*t) - 33.1*sin(-1.80 + 1.46*t) - 2.15*sin(-1.58 + 0.757*t) - 3.94*sin(-3.86 + 2.65*t) - 8.88*sin(1.51 + 1.79*t) - 9.97*sin(1.52 + 1.84*t) - 105.*sin(1.48 + 1.03*t) - 15.2*sin(-4.67 + 1.25*t) - 101.*sin(1.51 + 0.380*t) - 11.0*sin(-4.59 + 0.433*t) - 86.7*sin(1.50 + 0.973*t) - 170.*sin(1.53 + 0.597*t) - 41.2*sin(1.51 + 0.650*t) - 20.4*sin(-1.67 + 1.30*t) - 47.9*sin(-1.70 + 1.35*t) - 15.8*sin(-1.66 + 2.71*t) - 8.61*sin(-1.71 + 2.76*t) - 25.7*sin(-1.64 + 0.108*t) - 70.9*sin(1.55 + 0.162*t) - 0.668*sin(-2.42 + 2.92*t) - 4.78*sin(-4.60 + 2.98*t) - 106.*sin(1.49 + 0.920*t) - 17.6*sin(1.53 + 1.41*t) - 8.82*sin(1.05 + 2.17*t) - 113.*sin(-1.67 + 1.14*t), t <= 84, -800 - 7.30*sin(-171.17 + 2.90*t) - 3.28*sin(-6.550 + 0.121*t) - 1.46*sin(-17.878 + 0.241*t) - 20.4*sin(-22.438 + 0.361*t) - 28.9*sin(-29.862 + 0.484*t) - 9.13*sin(-36.364 + 0.603*t) - 45.3*sin(-40.650 + 0.725*t) - 97.4*sin(-50.770 + 0.845*t) - 13.1*sin(-54.916 + 0.967*t) - 80.8*sin(-61.97 + 1.09*t) - 39.1*sin(-71.92 + 1.21*t) - 42.8*sin(-78.87 + 1.33*t) - 108.*sin(-85.97 + 1.45*t) - 10.6*sin(-92.80 + 1.57*t) - 49.8*sin(-99.94 + 1.69*t) - 15.4*sin(-103.75 + 1.81*t) - 24.2*sin(-113.90 + 1.93*t) - 8.96*sin(-123.18 + 2.05*t) - 1.59*sin(-127.14 + 2.17*t) - 14.1*sin(-137.59 + 2.30*t) - 6.51*sin(-142.35 + 2.41*t) - 7.98*sin(-145.83 + 2.54*t) - 6.40*sin(-153.721 + 2.66*t) - 1.23*sin(-164.36 + 2.78*t), 84 < t, -1400 - 128.*sin(-32.358 + 0.402*t) - 68.5*sin(-43.928 + 0.502*t) - 2.55*sin(-242.18 + 2.86*t) - 6.86*sin(-219.136 + 2.61*t) - 5.76*sin(-222.904 + 2.66*t) - 2.39*sin(-226.835 + 2.71*t) - 101.*sin(-11.164 + 0.151*t) - 8.69*sin(-231.548 + 2.76*t) - 146.*sin(-31.268 + 0.352*t) - 8.30*sin(-179.37 + 2.11*t) - 2.68*sin(-261.69 + 3.06*t) - 10.4*sin(-162.98 + 1.91*t) - 30.1*sin(-73.606 + 0.854*t) - 24.1*sin(-77.946 + 0.904*t) - 10.0*sin(-146.01 + 1.71*t) - 72.5*sin(-69.416 + 0.804*t) - 8.91*sin(-85.97 + 1.00*t) - 8.58*sin(-175.51 + 2.06*t) - 27.4*sin(-109.01 + 1.31*t) - 16.8*sin(-113.17 + 1.36*t) - 162.*sin(-5.7968 + 0.0502*t) - 3.69*sin(-205.52 + 2.41*t) - 7.62*sin(-207.006 + 2.46*t) - 131.*sin(-53.522 + 0.653*t) - 95.3*sin(-60.882 + 0.703*t) - 8.53*sin(-197.627 + 2.36*t) - 1.74*sin(-247.32 + 2.91*t) - 27.2*sin(-121.51 + 1.46*t) - 51.7*sin(-49.332 + 0.603*t) - 8.81*sin(-104.925 + 1.26*t) - 10.2*sin(-100.703 + 1.21*t) - 9.35*sin(-183.90 + 2.16*t) - 7.82*sin(-188.20 + 2.21*t) - 42.8*sin(-26.964 + 0.301*t) - 16.8*sin(-48.312 + 0.553*t) - 15.2*sin(-9.980 + 0.100*t) - 213.*sin(-18.524 + 0.201*t) - 39.4*sin(-19.584 + 0.251*t) - 6.28*sin(-87.85 + 1.06*t) - 3.71*sin(-117.623 + 1.41*t) - 4.92*sin(-196.77 + 2.31*t) - 1.25*sin(-255.21 + 3.01*t) - 5.13*sin(-248.529 + 2.96*t) - 8.69*sin(-141.43 + 1.66*t) - 11.5*sin(-167.26 + 1.96*t) - 13.0*sin(-171.19 + 2.01*t) - 4.12*sin(-159.23 + 1.86*t) - 3.66*sin(-212.23 + 2.51*t) - 0.810*sin(-83.380 + 0.955*t) - 3.11*sin(-65.516 + 0.754*t) - 1.38*sin(-139.34 + 1.61*t) - 9.07*sin(-188.885 + 2.26*t) - 52.6*sin(-39.678 + 0.452*t) - 6.81*sin(-125.917 + 1.51*t) - 24.7*sin(-130.128 + 1.56*t) - 4.16*sin(-215.362 + 2.56*t) - 11.8*sin(-92.283 + 1.11*t) - 16.6*sin(-96.32 + 1.16*t) - 6.39*sin(-147.108 + 1.76*t) - 7.61*sin(-154.46 + 1.81*t) - 4.28*sin(-235.566 + 2.81*t)):

plot( [ x, y, t = 0 .. 146 ], scaling = constrained, discont = [ usefdiscont ], axes = boxed, thickness = 5, size = [600, 600]);

Drawing fireworks with Maple

by: mmcdara 7691 Product: Maple 2015

December 25 2019

10 3

Here is a little animation to wish all of you a Merry Christmas

FireWorks.mw

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post