MaplePrimes Posts

MaplePrimes Posts are for sharing your experiences, techniques and opinions about Maple, MapleSim and related products, as well as general interests in math and computing.

Latest Post
  • Latest Posts Feed
  • Formulating and solving the equivalence problem for Schwarzschild metric in a simple case

     

    In connection with the digitizing in Maple 2016 of the database of solutions to Einstein's equations of the book Exact Solutions to Einstein Field Equations. I was recently asked about a statement found in the "What is new in Physics in Maple 2016" page:

      

    In the Maple PDEtools package, you have the mathematical tools - including a complete symmetry approach - to work with the underlying [Einstein’s] partial differential equations. [By combining that functionality with the one in the Physics and Physics:-Tetrads package] you can also formulate and, depending on the metrics also resolve, the equivalence problem; that is: to answer whether or not, given two metrics, they can be obtained from each other by a transformation of coordinates, as well as compute the transformation.

    This question posed is a reasonable one: "could you please provide one example?" This post provides that example.

     

    First of all the existing science behind: in my opinion, the main reference regarding the equivalence problem is at the paper "A Review of the Geometrical Equivalence of Metrics in General Relativity", General Relativity and Gravitation, Vol. 12, No. 9, 1980, by A. Karlhede (University of Stockholm). This approach got refined later by others and, generally speaking, it is currently know as the Cartan-Karlhede method, summarized in chapter 9.2 of the book Exact Solutions to Einstein Field Equations. whose solutions were all digitized within the Physics and DifferentialGeometry packages for Maple 2016. This method of Chapter 9.2 (see also Tetrads and Weyl scalars in canonical form, Mapleprimes post), however, is not the only approach to the problem, and sometimes simpler methods can handle the problem faster, or just in simpler forms.

     

    The example worked out below is actually the example from Karlhede's paper just mentioned, on pages 704 - 706: "Show that the Schwarzschild metric and its form written in terms of isotropic spherical coordinates are equivalent, and derive the transformation that relates them". Because this problem happens to be simple for nowadays computer algebra, below I also tackle it modified, slightly more difficult variants of it. The approach shown works for more complicated cases as well.

     

    Below we tackle Karlhede's paper-problem using: one PDEtools command, the Physics:-TransformCoordinates, the Physics:-Weyl command to compute the Weyl scalars and the Physics:-Tetrads:-PetrovType to see the Petrov type of the metrics involved. The transformation resolving the equivalence is explicitly derived.

     

    Start loading the Physics and Tetrads package. To reproduce the computations below, as usual, update your Physics library with the one available for download at the Maplesoft R&D Physics webpage

    with(Physics); with(Tetrads); Setup(auto = true, tetradmetric = null, signature = `+---`)

    `Setting lowercaselatin letters to represent tetrad indices `

     

    0, "%1 is not a command in the %2 package", Tetrads, Physics

     

    0, "%1 is not a command in the %2 package", Tetrads, Physics

     

    `* Partial match of  'auto' against keyword 'automaticsimplification'`

     

    [automaticsimplification = true, signature = `+ - - -`, tetradmetric = {(1, 2) = 1, (3, 4) = -1}]

    (1)

    To formulate the problem, set first some symbols to represent the changed metric, changed mass and changed coordinates - no mathematics at this point

    gt, mt, tt, rt, thetat, phit := `𝔤`, `𝔪`, `𝔱`, `𝔯`, `ϑ`, `ϕ`

    `𝔤`, `𝔪`, `𝔱`, `𝔯`, vartheta, varphi

    (2)

    Set now a new coordinates system, call it Y, involving the new coordinates (in the paper they are represented with a tilde on top of the letters)

    Coordinates(Y = [tt, rt, thetat, phit])

    `Default differentiation variables for d_, D_ and dAlembertian are: `*{Y = (`𝔱`, `𝔯`, `ϑ`, `ϕ`)}

     

    `Systems of spacetime Coordinates are: `*{Y = (`𝔱`, `𝔯`, `ϑ`, `ϕ`)}

     

    {Y}

    (3)

    According to eq.(7.6) of the paper, the line element of Schwarzschild solution in isotropic spherical coordinates is given by

    `#msup(mi("ds"),mn("2"))` := ((1-mt/(2*rt))/(1+mt/(2*rt)))^2*d_(tt)^2-(1+mt/(2*rt))^4*(d_(rt)^2+rt^2*d_(thetat)^2+rt^2*sin(thetat)^2*d_(phit)^2)

    (-2*`𝔯`+`𝔪`)^2*Physics:-d_(`𝔱`)^2/(2*`𝔯`+`𝔪`)^2-(1/16)*(2*`𝔯`+`𝔪`)^4*(Physics:-d_(`𝔯`)^2+`𝔯`^2*Physics:-d_(vartheta)^2+`𝔯`^2*sin(vartheta)^2*Physics:-d_(varphi)^2)/`𝔯`^4

    (4)

    Set this to be the metric

    Setup(metric = `#msup(mi("ds"),mn("2"))`)

    Check it out

    g_[]

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

    (5)

    In connection with the transformation used further below, compute now the Petrov type and the Weyl scalars for this metric, just to have an idea of what is behind this metric.

    PetrovType()

    "D"

    (6)

    Weyl[scalars]

    psi__0 = 0, psi__1 = 0, psi__2 = -64*`𝔯`^3*`𝔪`/(2*`𝔯`+`𝔪`)^6, psi__3 = 0, psi__4 = 0

    (7)

    We see that the Weyl scalars are already in canonical form (see post in Mapleprimes about canonical forms): only `&Psi;__2` <> 0 and the important thing: it depends on only one coordinate, `&rfr;` .

     

    Now: we want to see if this metric (5) is equivalent to Schwarzschild metric in standard spherical coordinates

    g_[sc]

    `Systems of spacetime Coordinates are: `*{X = (t, r, theta, phi), Y = (`&tfr;`, `&rfr;`, `&vartheta;`, `&varphi;`)}

     

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

     

    `The Schwarzschild metric in coordinates `[t, r, theta, phi]

     

    `Parameters: `[m]

     

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

    (8)

    The equivalence we want to resolve is regarding an arbitrary relationship `&mfr;`(m)between the masses used in (5) and (8) and a generic change of variables from X to Y

    TR := {phi = Phi(Y), r = R(Y), t = Tau(Y), theta = Theta(Y)}

    {phi = Phi(Y), r = R(Y), t = Tau(Y), theta = Theta(Y)}

    (9)

    Using a differential equation mindset, the formulation of the equivalence between (8) and (5) under the transformation (9) is actually simple: change variables in (8), using (9) and the Physics:-TransformCoordinates command (this is the command that changes variables in tensorial expressions), then equate the result to (5), then try to solve the problem for the unknowns `&mfr;`(m), Phi(Y), R(Y), Theta(Y) and Tau(Y).

     

    We note at this point, however, that the Weyl scalars for Schwarzschild metric in this standard form (8) are also in canonical form of Petrov type D and also depend on only one variable, r 

    PetrovType()

    "D"

    (10)

    Weyl[scalars]

    psi__0 = 0, psi__1 = 0, psi__2 = -m/r^3, psi__3 = 0, psi__4 = 0

    (11)

    The fact that the Weyl scalars in both cases ((7) and (11)) are in canonical form (only `&Psi;__2` <> 0 ) and in both cases this scalar depends on only one coordinate is already an indicator that the transformation involved changes only one variable in terms of the other one. So one could just search for a transformation of the form r = R(`&rfr;`) and resolve the problem instantly. Still, to make the problem slightly more general, consider instead a generic transformation for r in terms of all of Y = (`&tfr;`, `&rfr;`, `&vartheta;`, `&varphi;`)

    tr := r = R(Y)

    r = R(Y)

    (12)

    PDEtools:-declare(r = R(Y))

    R(`&tfr;`, `&rfr;`, vartheta, varphi)*`will now be displayed as`*R

    (13)

    Transform the  coordinates in the metric (because of having used PDEtools:-declare, derivatives of the unknowns R are displayed indexed, for compact notation)

    TransformCoordinates(tr, g_[mu, nu])

    Matrix(%id = 18446744078873927542)

    (14)

    Proceed equating (14) to (5) to obtain a set of equations that entirely formulates the problem

    "convert(rhs(?)=? ,setofequations)"

    {0 = (diff(R(Y), `&rfr;`))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), vartheta))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), (-2*`&rfr;`+`&mfr;`)^2/(2*`&rfr;`+`&mfr;`)^2 = ((diff(R(Y), `&tfr;`))^2*R(Y)^2-4*(-(1/2)*R(Y)+m)^2)/(R(Y)*(-R(Y)+2*m)), -(1/16)*(2*`&rfr;`+`&mfr;`)^4/`&rfr;`^4 = (diff(R(Y), `&rfr;`))^2*R(Y)/(-R(Y)+2*m), -(1/16)*(2*`&rfr;`+`&mfr;`)^4/`&rfr;`^2 = -(diff(R(Y), vartheta))^2*R(Y)/(R(Y)-2*m)-R(Y)^2, -(1/16)*(2*`&rfr;`+`&mfr;`)^4*sin(vartheta)^2/`&rfr;`^2 = 2*((1/2)*(diff(R(Y), varphi))^2+(cos(vartheta)-1)*R(Y)*(cos(vartheta)+1)*(-(1/2)*R(Y)+m))*R(Y)/(-R(Y)+2*m)}

    (15)

    This problem, shown in Karlhede's paper as the example of the approach he summarized, is solvable using the differential equation commands of PDEtools (in this case casesplit) in one go and no time, obtaining the same solution shown in the paper with equation number (7.10), the problem actually admits two solutions

    PDEtools:-casesplit({0 = (diff(R(Y), `&rfr;`))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), vartheta))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), (-2*`&rfr;`+`&mfr;`)^2/(2*`&rfr;`+`&mfr;`)^2 = ((diff(R(Y), `&tfr;`))^2*R(Y)^2-4*(-(1/2)*R(Y)+m)^2)/(R(Y)*(-R(Y)+2*m)), -(1/16)*(2*`&rfr;`+`&mfr;`)^4/`&rfr;`^4 = (diff(R(Y), `&rfr;`))^2*R(Y)/(-R(Y)+2*m), -(1/16)*(2*`&rfr;`+`&mfr;`)^4/`&rfr;`^2 = -(diff(R(Y), vartheta))^2*R(Y)/(R(Y)-2*m)-R(Y)^2, -(1/16)*(2*`&rfr;`+`&mfr;`)^4*sin(vartheta)^2/`&rfr;`^2 = 2*((1/2)*(diff(R(Y), varphi))^2+(cos(vartheta)-1)*R(Y)*(cos(vartheta)+1)*(-(1/2)*R(Y)+m))*R(Y)/(-R(Y)+2*m)}, [R, mt])

    `casesplit/ans`([R(Y) = -(1/4)*(m-2*`&rfr;`)^2/`&rfr;`, `&mfr;` = -m], []), `casesplit/ans`([R(Y) = (1/4)*(2*`&rfr;`+m)^2/`&rfr;`, `&mfr;` = m], [])

    (16)

    By all means this does not mean this differential equation approach is better than the general approach mentioned in the paper (also in section 9.2 of the Exact Solutions book). This presentation above only makes the point of the paragraph mentioned at the beginning of this worksheet "... [in Maple 2016] you can also formulate and, depending on the the metrics also resolve, the equivalence problem; that is: to answer whether or not, given two metrics, they can be obtained from each other by a transformation of coordinates, as well as compute the transformation." 

     

    In any case this problem above is rather easy for the computer. Consider a slightly more difficult problem, where `&mfr;` <> m. For example:

    "subs(mt = 1/(mt^(2)),?)"

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

    (17)

    Tackle now the same problem

    "convert(rhs(?)=? ,setofequations)"

    {0 = (diff(R(Y), `&rfr;`))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), vartheta))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), (-2*`&rfr;`+1/`&mfr;`^2)^2/(2*`&rfr;`+1/`&mfr;`^2)^2 = ((diff(R(Y), `&tfr;`))^2*R(Y)^2-4*(-(1/2)*R(Y)+m)^2)/(R(Y)*(-R(Y)+2*m)), -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4/`&rfr;`^4 = (diff(R(Y), `&rfr;`))^2*R(Y)/(-R(Y)+2*m), -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4/`&rfr;`^2 = -(diff(R(Y), vartheta))^2*R(Y)/(R(Y)-2*m)-R(Y)^2, -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4*sin(vartheta)^2/`&rfr;`^2 = 2*((1/2)*(diff(R(Y), varphi))^2+(cos(vartheta)-1)*R(Y)*(cos(vartheta)+1)*(-(1/2)*R(Y)+m))*R(Y)/(-R(Y)+2*m)}

    (18)

    The solutions to the equivalence between (17) and (5) are then given by

    PDEtools:-casesplit({0 = (diff(R(Y), `&rfr;`))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), vartheta))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), (-2*`&rfr;`+1/`&mfr;`^2)^2/(2*`&rfr;`+1/`&mfr;`^2)^2 = ((diff(R(Y), `&tfr;`))^2*R(Y)^2-4*(-(1/2)*R(Y)+m)^2)/(R(Y)*(-R(Y)+2*m)), -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4/`&rfr;`^4 = (diff(R(Y), `&rfr;`))^2*R(Y)/(-R(Y)+2*m), -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4/`&rfr;`^2 = -(diff(R(Y), vartheta))^2*R(Y)/(R(Y)-2*m)-R(Y)^2, -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4*sin(vartheta)^2/`&rfr;`^2 = 2*((1/2)*(diff(R(Y), varphi))^2+(cos(vartheta)-1)*R(Y)*(cos(vartheta)+1)*(-(1/2)*R(Y)+m))*R(Y)/(-R(Y)+2*m)}, [R, mt])

    `casesplit/ans`([R(Y) = -(1/4)*(m-2*`&rfr;`)^2/`&rfr;`, `&mfr;`^2 = -1/m], [`&mfr;` <> 0]), `casesplit/ans`([R(Y) = (1/4)*(2*`&rfr;`+m)^2/`&rfr;`, `&mfr;`^2 = 1/m], [`&mfr;` <> 0])

    (19)

    Moreover, despite that the Weyl scalars suggest that a transformation of only one variable is sufficient to solve the problem, one could also consider a more general transformation, of more variables. Provided we exclude theta (because there is cos(theta) around and that would take us to solve differential equations for Theta(theta), that involve things like cos(Theta(theta))), and also to speed up matters let's remove the change in phi, consider an arbitrary change in r and t

    TR := select(has, {phi = Phi(Y), r = R(Y), t = Tau(Y), theta = Theta(Y)}, {r, t})

    {r = R(Y), t = Tau(Y)}

    (20)

    PDEtools:-declare({r = R(Y), t = Tau(Y)})

    R(`&tfr;`, `&rfr;`, vartheta, varphi)*`will now be displayed as`*R

     

    Tau(`&tfr;`, `&rfr;`, vartheta, varphi)*`will now be displayed as`*Tau

    (21)

    So our transformation now involve two arbitrary variables, each one depending on all the four coordinates, and a more complicated function `&mfr;`(m). Change variables (because of having used PDEtools:-declare, derivatives of the unknowns R and Tau are displayed indexed, for compact notation)

    TransformCoordinates(TR, g_[mu, nu])

    Matrix(%id = 18446744078309268046)

    (22)

    Construct the set of Partial Differential Equations to be tackled

    "convert(rhs(?)=?,setofequations)"

    {0 = (-4*(diff(Tau(Y), `&rfr;`))*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), `&tfr;`))+(diff(R(Y), `&rfr;`))*(diff(R(Y), `&tfr;`))*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), 0 = (-4*(diff(Tau(Y), varphi))*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), `&rfr;`))+(diff(R(Y), varphi))*(diff(R(Y), `&rfr;`))*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), 0 = (-4*(diff(Tau(Y), varphi))*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), `&tfr;`))+(diff(R(Y), varphi))*(diff(R(Y), `&tfr;`))*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), 0 = (-4*(diff(Tau(Y), varphi))*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), vartheta))+(diff(R(Y), varphi))*(diff(R(Y), vartheta))*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), 0 = (-4*(diff(Tau(Y), vartheta))*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), `&rfr;`))+(diff(R(Y), vartheta))*(diff(R(Y), `&rfr;`))*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), 0 = (-4*(diff(Tau(Y), vartheta))*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), `&tfr;`))+(diff(R(Y), vartheta))*(diff(R(Y), `&tfr;`))*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), (-2*`&rfr;`+1/`&mfr;`^2)^2/(2*`&rfr;`+1/`&mfr;`^2)^2 = (-4*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), `&tfr;`))^2+(diff(R(Y), `&tfr;`))^2*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4/`&rfr;`^4 = (-4*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), `&rfr;`))^2+(diff(R(Y), `&rfr;`))^2*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4/`&rfr;`^2 = (diff(Tau(Y), vartheta))^2*(R(Y)-2*m)/R(Y)-(diff(R(Y), vartheta))^2*R(Y)/(R(Y)-2*m)-R(Y)^2, -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4*sin(vartheta)^2/`&rfr;`^2 = (-4*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), varphi))^2+2*((1/2)*(diff(R(Y), varphi))^2+(cos(vartheta)-1)*R(Y)*(cos(vartheta)+1)*(-(1/2)*R(Y)+m))*R(Y)^2)/(R(Y)*(-R(Y)+2*m))}

    (23)

    Solve the problem running a differential elimination (actually without solving any differential equations): there are more than two solutions

    sol := PDEtools:-casesplit({0 = (-4*(diff(Tau(Y), `&rfr;`))*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), `&tfr;`))+(diff(R(Y), `&rfr;`))*(diff(R(Y), `&tfr;`))*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), 0 = (-4*(diff(Tau(Y), varphi))*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), `&rfr;`))+(diff(R(Y), varphi))*(diff(R(Y), `&rfr;`))*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), 0 = (-4*(diff(Tau(Y), varphi))*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), `&tfr;`))+(diff(R(Y), varphi))*(diff(R(Y), `&tfr;`))*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), 0 = (-4*(diff(Tau(Y), varphi))*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), vartheta))+(diff(R(Y), varphi))*(diff(R(Y), vartheta))*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), 0 = (-4*(diff(Tau(Y), vartheta))*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), `&rfr;`))+(diff(R(Y), vartheta))*(diff(R(Y), `&rfr;`))*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), 0 = (-4*(diff(Tau(Y), vartheta))*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), `&tfr;`))+(diff(R(Y), vartheta))*(diff(R(Y), `&tfr;`))*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), (-2*`&rfr;`+1/`&mfr;`^2)^2/(2*`&rfr;`+1/`&mfr;`^2)^2 = (-4*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), `&tfr;`))^2+(diff(R(Y), `&tfr;`))^2*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4/`&rfr;`^4 = (-4*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), `&rfr;`))^2+(diff(R(Y), `&rfr;`))^2*R(Y)^2)/(R(Y)*(-R(Y)+2*m)), -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4/`&rfr;`^2 = (diff(Tau(Y), vartheta))^2*(R(Y)-2*m)/R(Y)-(diff(R(Y), vartheta))^2*R(Y)/(R(Y)-2*m)-R(Y)^2, -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4*sin(vartheta)^2/`&rfr;`^2 = (-4*(-(1/2)*R(Y)+m)^2*(diff(Tau(Y), varphi))^2+2*((1/2)*(diff(R(Y), varphi))^2+(cos(vartheta)-1)*R(Y)*(cos(vartheta)+1)*(-(1/2)*R(Y)+m))*R(Y)^2)/(R(Y)*(-R(Y)+2*m))}, [R, mt])

    `casesplit/ans`([R(Y) = (1/4)*(2*`&rfr;`+m)^2/`&rfr;`, `&mfr;`^2 = 1/m, diff(Tau(Y), `&tfr;`) = -1, diff(Tau(Y), `&rfr;`) = 0, diff(Tau(Y), vartheta) = 0, diff(Tau(Y), varphi) = 0], [`&mfr;` <> 0]), `casesplit/ans`([R(Y) = (1/4)*(2*`&rfr;`+m)^2/`&rfr;`, `&mfr;`^2 = 1/m, diff(Tau(Y), `&tfr;`) = 1, diff(Tau(Y), `&rfr;`) = 0, diff(Tau(Y), vartheta) = 0, diff(Tau(Y), varphi) = 0], [`&mfr;` <> 0]), `casesplit/ans`([R(Y) = -(1/4)*(m-2*`&rfr;`)^2/`&rfr;`, `&mfr;`^2 = -1/m, diff(Tau(Y), `&tfr;`) = -1, diff(Tau(Y), `&rfr;`) = 0, diff(Tau(Y), vartheta) = 0, diff(Tau(Y), varphi) = 0], [`&mfr;` <> 0]), `casesplit/ans`([R(Y) = -(1/4)*(m-2*`&rfr;`)^2/`&rfr;`, `&mfr;`^2 = -1/m, diff(Tau(Y), `&tfr;`) = 1, diff(Tau(Y), `&rfr;`) = 0, diff(Tau(Y), vartheta) = 0, diff(Tau(Y), varphi) = 0], [`&mfr;` <> 0])

    (24)

    Consider for instance the first one

    sol[1]

    `casesplit/ans`([R(Y) = (1/4)*(2*`&rfr;`+m)^2/`&rfr;`, `&mfr;`^2 = 1/m, diff(Tau(Y), `&tfr;`) = -1, diff(Tau(Y), `&rfr;`) = 0, diff(Tau(Y), vartheta) = 0, diff(Tau(Y), varphi) = 0], [`&mfr;` <> 0])

    (25)

    Compute the actual solution behind this case :

    pdsolve(`casesplit/ans`([R(Y) = (1/4)*(2*`&rfr;`+m)^2/`&rfr;`, `&mfr;`^2 = 1/m, diff(Tau(Y), `&tfr;`) = -1, diff(Tau(Y), `&rfr;`) = 0, diff(Tau(Y), vartheta) = 0, diff(Tau(Y), varphi) = 0], [`&mfr;` <> 0]), {R, Tau, mt})

    {`&mfr;` = -1/m^(1/2), R(Y) = (1/4)*(2*`&rfr;`+m)^2/`&rfr;`, Tau(Y) = -`&tfr;`+_C1}, {`&mfr;` = 1/m^(1/2), R(Y) = (1/4)*(2*`&rfr;`+m)^2/`&rfr;`, Tau(Y) = -`&tfr;`+_C1}

    (26)

    The fact that the time t appears defined in terms of the transformed time Tau(Y) = -`&tfr;`+_C1 involving an arbitrary constant is expected: the time does not enter the metric, it only enters through derivatives of Tau(Y) entering the Jacobian of the transformation used to change variables in tensorial expressions (the metric) in (22).

     

    Summary: the approach shown above, based on formulating the problem for the transformation functions of the equivalence and solving for them the differential equations using the commands in PDEtools, after restricting the generality of the transformation functions by looking at the form of the Weyl scalars, works well for other cases too, specially now that, in Maple 2016, the Weyl scalars can be expressed also in canonical form in one go (see previous Mapleprimes post on "Tetrads and Weyl scalars in canonical form").  Also important: in Maple 2016 it is present the functionality necessary to implement the approach of section 9.2 of the Exact solutions book as well.

      

     

     

    Download Equivalence_-_Schwarzschild.mw

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

    Tetrads and Weyl scalars in canonical form

     

    The material below is about a new development that didn't arrive in time for the launch of Maple 2016 (March) and that complements in a relevant way the ones introduced in Physics in Maple 2016. It is at topic in general relativity, the computation of a canonical form of a tetrad, so that, generally speaking (skipping a technical description) the Weyl scalars are fixed as much as possible (either equal to 0 or to 1) regarding transformations that leave invariant the tetrad metric in a tetrad system of references. Bringing a tetrad in canonical form is a relevant step in the tackling of the equivalence problem between two spacetime metrics (Mapleprimes post), and it is relevant in connection with the digitizing in Maple 2016 of the database of solutions to Einstein's equations of the book Exact Solutions to Einstein Field Equations.

    The reference for this development is the book "General Relativity, an Einstein century survey", edited by S.W. Hawking (Cambridge) and W. Israel (U. Alberta, Canada), specifically Chapter 7 written by S. Chandrasekhar, and more specifically exploring what is said in page 388 about the Petrov classification.


    A canonical form for the tetrad and Weyl scalars admits alternate forms; the implementation is as implicit in page 388:

     

    `&Psi;__0`

    `&Psi;__1`

    `&Psi;__2`

    `&Psi;__3`

    `&Psi;__4`

    Residual invariance

    Petrov type I

    0

    "<>0"

    "<>0"

    1

    0

    none

    Petrov type II

    0

    0

    "<>0"

    1

    0

    none

    Petrov type III

    0

    0

    0

    1

    0

    none

    Petrov type D

    0

    0

    "<>0"

    0

    0

    `&Psi;__2`  remains invariant under rotations of Class III

    Petrov type N

    0

    0

    0

    0

    1

    `&Psi;__4` remains invariant under rotations of Class II

     

    The transformations (rotations of the tetrad system of references) used are of Class I, II and III as defined in Chandrasekar's chapter - equations (7.79) in page 384, (7.83) and (7.84) in page 385. Transformations of Class I can be performed with the command Physics:-Tetrads:-TransformTetrad using the optional argument nullrotationwithfixedl_, of Class II using nullrotationwithfixedn_ and of Class III by calling TransformTetrad(spatialrotationsm_mb_plan, boostsn_l_plane), so with the two optional arguments simultaneously.

     

    In this development, a new optional argument, canonicalform got implemented to TransformTetrad so that the whole sequence of three transformations of Classes I, II and III is performed automatically, in one go. Regarding the canonical form of the tetrad, the main idea is that from the change in the Weyl scalars one can derive the parameters entering tetrad transformations that result in a canonical form of the tetrad. 

     

    with(Physics); with(Tetrads)

    `Setting lowercaselatin letters to represent tetrad indices `

     

    0, "%1 is not a command in the %2 package", Tetrads, Physics

     

    0, "%1 is not a command in the %2 package", Tetrads, Physics

     

    [IsTetrad, NullTetrad, OrthonormalTetrad, PetrovType, SimplifyTetrad, TransformTetrad, e_, eta_, gamma_, l_, lambda_, m_, mb_, n_]

    (1)

    (Note the Tetrads:-PetrovType command, unfinished in the first release of Maple 2016.) To run the following computations you need to update your Physics library to the latest version from the Maplesoft R&D Physics webpage, so with this datestamp or newer:

    Physics:-Version()

    "/Users/ecterrab/Maple/lib/Physics2016.mla", `2016, April 20, 12:56 hours`

    (2)

    An Example of Petrov type I

    There are six Petrov types: I, II, III, D, N and O. Start with a spacetime metric of Petrov type "I"  (the numbers always refer to the equation number in the "Exact solutions to Einstein's field equations" textbook)

    g_[[12, 21, 1]]

    `Systems of spacetime Coordinates are: `*{X = (t, x, y, phi)}

     

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

     

    `The McLenaghan, Tariq (1975), Tupper (1976) metric in coordinates `[t, x, y, phi]

     

    `Parameters: `[a, k, kappa0]

     

    "`Comments: `_k parametrizes the most general electromagnetic invariant with respect to the last 3 Killing vectors"

     

    `Resetting the signature of spacetime from "+ - - -" to \`- + + +\` in order to match the signature in the database of metrics:`

     

    g[mu, nu] = (Matrix(4, 4, {(1, 1) = -1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 2*y, (2, 1) = 0, (2, 2) = a^2/x^2, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = a^2/x^2, (3, 4) = 0, (4, 1) = 2*y, (4, 2) = 0, (4, 3) = 0, (4, 4) = x^2-4*y^2}))

    (3)

    The Weyl scalars

    Weyl[scalars]

    psi__0 = (1/4)*((4*I)*x^3*abs(x)^3-abs(x)^6+abs(x)^4*x^2+abs(x)^2*x^4-x^6)/(a^2*abs(x)^4*x^2), psi__1 = 0, psi__2 = -(1/4)*(x^2+abs(x)^2)*(x^4+abs(x)^4)/(a^2*abs(x)^4*x^2), psi__3 = 0, psi__4 = (1/4)*((4*I)*x^3*abs(x)^3-abs(x)^6+abs(x)^4*x^2+abs(x)^2*x^4-x^6)/(a^2*abs(x)^4*x^2)

    (4)

    ... there is abs around. Let's assume everything is positive to simplify formulas, use Capital Physics:-Assume  (the lower case assume  command redefines the assumed variables, so it is not compatible with Physics, DifferentialGeometry and VectorCalculus among others).

    Assume(x > 0, y > 0, a > 0)

    {a::(RealRange(Open(0), infinity))}, {x::(RealRange(Open(0), infinity))}, {y::(RealRange(Open(0), infinity))}

    (5)

    The scalars are now simpler, although still not in "canonical form" because `&Psi;__4` <> 0 and `&Psi;__3` <> 1.

    Weyl[scalars]

    psi__0 = I/a^2, psi__1 = 0, psi__2 = -1/a^2, psi__3 = 0, psi__4 = I/a^2

    (6)

    The Petrov type

    PetrovType()

    "I"

    (7)

    The  call to Tetrads:-TransformTetrad two lines below transforms the current tetrad ,

    e_[]

    Physics:-Tetrads:-e_[a, mu] = Matrix(%id = 18446744078512745638)

    (8)

    into another tetrad such that the Weyl scalars are in canonical form, which for Petrov "I" type happens when `&Psi;__0` = 0, `&Psi;__4` = 0 and `&Psi;__3` = 1.

    TransformTetrad(canonicalform)

    Matrix(%id = 18446744078500192254)

    (9)

    Despite the fact that the result is a much more complicated tetrad, this is an amazing result in that the resulting Weyl scalars are all fixed (see below).  Let's first verify that this is indeed a tetrad, and that now the Weyl scalars are in canonical form

    "IsTetrad(?)"

    `Type of tetrad: null `

     

    true

    (10)

    Set (9) to be the tetrad in use and recompute the Weyl scalars

    "Setup(tetrad = ?):"

    Inded we now have `&Psi;__0` = 0, `&Psi;__4` = 0 and `&Psi;__3` = 1 

    simplify([Weyl[scalars]])

    [psi__0 = 0, psi__1 = (-1/2-(3/2)*I)/a^4, psi__2 = (-1+I)/a^2, psi__3 = 1, psi__4 = 0]

    (11)

    So Weyl scalars computed after setting the canonical tetrad (9) to be the tetrad in use are in canonical form. Great! NOTE: computing the canonicalWeyl scalars is not really the difficult part, and within the code, these scalars (11) are computed before arriving at the tetrad (9). What is really difficult (from the point of view of computational complexity and simplifications) is to compute the actual canonical form of the tetrad (9).

     

    An Example of Petrov type II

    Consider this other solution to Einstein's equation (again, the numbers in g_[[24,37,7]] always refer to the equation number in the "Exact solutions to Einstein's field equations" textbook)

    g_[[24, 37, 7]]

    `Systems of spacetime Coordinates are: `*{X = (u, v, x, y)}

     

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

     

    `The Stephani metric in coordinates `[u, v, x, y]

     

    `Parameters: `[f(x), a, Psi1(u, x, y)]

     

    "`Comments: `Case 6 from Table 24.1:_Psi1(u,x,y): diff(_Psi1(u,x,y),x,x)+diff(_Psi1(u,x,y),y,y)=0, diff(x*diff(_M(u,x,y),x),x)+x*diff(_M(u,x,y),y,y)=_kappa0*(diff(_Psi(u,x,y),x)^2+diff(_Psi(u,x,y),y)^2)"

     

    g[mu, nu] = (Matrix(4, 4, {(1, 1) = -2*x*(f(x)+y*a), (1, 2) = -x, (1, 3) = 0, (1, 4) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 3) = 1/x^(1/2), (3, 4) = 0, (4, 4) = 1/x^(1/2)}, storage = triangular[upper], shape = [symmetric]))

    (12)

    Check the Petrov type

    PetrovType()

    "II"

    (13)

    The starting tetrad

    e_[]

    Physics:-Tetrads:-e_[a, mu] = Matrix(%id = 18446744078835577550)

    (14)

    results in Weyl scalars not in canonical form:

    Weyl[scalars]

    psi__0 = 0, psi__1 = 0, psi__2 = (1/8)/x^(3/2), psi__3 = 0, psi__4 = -((3*I)*a-2*x*(diff(diff(f(x), x), x))-3*(diff(f(x), x)))/(x^(1/2)*(4*y*a+4*f(x)))

    (15)

    For Petrov type "II", the canonical form is as for type "I" but in addition `&Psi;__1` = 0. Again let's assume positive, not necessary, but to get simpler formulas around

    Assume(f(x) > 0, x > 0, y > 0, a > 0)

    {a::(RealRange(Open(0), infinity))}, {x::(RealRange(Open(0), infinity)), (-f(x))::(RealRange(-infinity, Open(0))), (f(x))::(RealRange(Open(0), infinity))}, {y::(RealRange(Open(0), infinity))}

    (16)

    Compute now a canonical form for the tetrad, to be used instead of (14)

    TransformTetrad(canonicalform)

    Matrix(%id = 18446744078835949430)

    (17)

    Set this tetrad and check the Weyl scalars again

    "Setup(tetrad = ?):"

    Weyl[scalars]

    psi__0 = 0, psi__1 = 0, psi__2 = (1/8)/x^(3/2), psi__3 = 1, psi__4 = 0

    (18)

    This result (18) is fantastic. Compare these Weyl scalars with the ones (15) before transforming the tetrad.

     

    An Example of Petrov type III

    g_[[12, 35, 1]]

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

     

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

     

    `The Kaigorodov (1962), Cahen (1964), Siklos (1981), Ozsvath (1987) metric in coordinates `[u, x, y, z]

     

    `Parameters: `[Lambda]

     

    g[mu, nu] = (Matrix(4, 4, {(1, 1) = 0, (1, 2) = exp(-2*z), (1, 3) = 0, (1, 4) = 0, (2, 2) = exp(4*z), (2, 3) = 2*exp(z), (2, 4) = 0, (3, 3) = 2*exp(-2*z), (3, 4) = 0, (4, 4) = 3/abs(Lambda)}, storage = triangular[upper], shape = [symmetric]))

    (19)

    Assume(z > 0, Lambda > 0)

    {Lambda::(RealRange(Open(0), infinity))}, {z::(RealRange(Open(0), infinity))}

    (20)

    The Petrov type and the original tetrad

    PetrovType()

    "III"

    (21)

    e_[]

    Physics:-Tetrads:-e_[a, mu] = Matrix(%id = 18446744078349449926)

    (22)

    This tetrad results in the following scalars

    Weyl[scalars]

    psi__0 = -2*Lambda*2^(1/2)+(11/4)*Lambda, psi__1 = -(1/2)*Lambda*2^(1/2)+(3/4)*Lambda, psi__2 = (1/4)*Lambda, psi__3 = -(1/2)*Lambda*2^(1/2)-(3/4)*Lambda, psi__4 = 2*Lambda*2^(1/2)+(11/4)*Lambda

    (23)

    that are not in canonical form, which for Petrov type III is as in Petrov type II but in addition we should have `&Psi;__2` = 0.

    Compute now a canonical form for the tetrad

    TransformTetrad(canonicalform)

    Matrix(%id = 18446744078500057566)

    (24)

    Set this one to be the tetrad in use and recompute the Weyl scalars

    "Setup(tetrad = ?):"

    Weyl[scalars]

    psi__0 = 0, psi__1 = 0, psi__2 = 0, psi__3 = 1, psi__4 = 0

    (25)

    Great!``

    An Example of Petrov type N

    g_[[12, 6, 1]]

    `Systems of spacetime Coordinates are: `*{X = (u, v, y, z)}

     

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

     

    `The Defrise (1969) metric in coordinates `[u, v, y, z]

     

    `Parameters: `[Lambda, kappa0]

     

    "`Comments: `_Lambda < 0 required for a pure radiation solution"

     

    g[mu, nu] = (Matrix(4, 4, {(1, 1) = 0, (1, 2) = -(3/2)/(y^2*Lambda), (1, 3) = 0, (1, 4) = 0, (2, 2) = -3/(y^4*Lambda), (2, 3) = 0, (2, 4) = 0, (3, 3) = 3/(y^2*Lambda), (3, 4) = 0, (4, 4) = 3/(y^2*Lambda)}, storage = triangular[upper], shape = [symmetric]))

    (26)

    Assume(y > 0, Lambda > 0)

    {Lambda::(RealRange(Open(0), infinity))}, {y::(RealRange(Open(0), infinity))}

    (27)

    PetrovType()

    "N"

    (28)

    The original tetrad and related Weyl scalars are not in canonical form:

    e_[]

    Physics:-Tetrads:-e_[a, mu] = Matrix(%id = 18446744078404437406)

    (29)

    Weyl[scalars]

    psi__0 = -(1/4)*Lambda, psi__1 = -((1/4)*I)*Lambda, psi__2 = (1/4)*Lambda, psi__3 = ((1/4)*I)*Lambda, psi__4 = -(1/4)*Lambda

    (30)

    For Petrov type "N", the canonical form has `&Psi;__4` <> 0 and all the other `&Psi;__n` = 0.

    Compute a canonical form, set it to be the tetrad in use and recompute the Weyl scalars

    TransformTetrad(canonicalform)

    Matrix(%id = 18446744078518486190)

    (31)

    "Setup(tetrad = ?):"

    Weyl[scalars]

    psi__0 = 0, psi__1 = 0, psi__2 = 0, psi__3 = 0, psi__4 = 1

    (32)

    All as expected.

    An Example of Petrov type D

     

    g_[[12, 8, 4]]

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

     

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

     

    `The  metric in coordinates `[t, x, y, z]

     

    `Parameters: `[A, B]

     

    "`Comments: `k = 0, kprime = 1, not an Einstein metric"

     

    g[mu, nu] = (Matrix(4, 4, {(1, 1) = -B^2*sin(z)^2, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 2) = A^2, (2, 3) = 0, (2, 4) = 0, (3, 3) = A^2*x^2, (3, 4) = 0, (4, 4) = B^2}, storage = triangular[upper], shape = [symmetric]))

    (33)

    Assume(A > 0, B > 0, x > 0, 0 <= z and z <= (1/4)*Pi)

    {A::(RealRange(Open(0), infinity))}, {B::(RealRange(Open(0), infinity))}, {x::(RealRange(Open(0), infinity))}, {z::(RealRange(0, (1/4)*Pi))}

    (34)

    PetrovType()

    "D"

    (35)

    The default tetrad and related Weyl scalars are not in canonical form, which for Petrov type "D" is with `&Psi;__2` <> 0 and all the other `&Psi;__n` = 0

    e_[]

    Physics:-Tetrads:-e_[a, mu] = Matrix(%id = 18446744078503920694)

    (36)

    Weyl[scalars]

    psi__0 = (1/4)/B^2, psi__1 = 0, psi__2 = (1/12)/B^2, psi__3 = 0, psi__4 = (1/4)/B^2

    (37)

    Transform the  tetrad, set it and recompute the Weyl scalars

    TransformTetrad(canonicalform)

    Matrix(%id = 18446744078814996830)

    (38)

    "Setup(tetrad=?):"

    Weyl[scalars]

    psi__0 = 0, psi__1 = 0, psi__2 = -(1/6)/B^2, psi__3 = 0, psi__4 = 0

    (39)

    Again the expected canonical form of the Weyl scalars, and `&Psi;__2` <> 0 remains invariant under transformations of Class III.

     

    An Example of Petrov type O

     

    Finally an example of type "O". This corresponds to a conformally flat spacetime, for which the Weyl tensor (and with it all the Weyl scalars) vanishes. So the code just interrupts with "not implemented for conformally flat spactimes of Petrov type O"

    g_[[8, 33, 1]]

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

     

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

     

    `The  metric in coordinates `[t, x, y, z]

     

    `Parameters: `[K]

     

    "`Comments: `_K=3*_Lambda, _K>0 de Sitter, _K<0 anti-de Sitte"

     

    g[mu, nu] = z

    (40)

    PetrovType()

    "O"

    (41)

    The Weyl tensor and its scalars all vanish:

    Weyl[nonzero]

    Physics:-Weyl[mu, nu, alpha, beta] = {}

    (42)

    simplify(evala([Weyl[scalars]]))

    [psi__0 = 0, psi__1 = 0, psi__2 = 0, psi__3 = 0, psi__4 = 0]

    (43)

    TransformTetrad(canonicalform)

    Error, (in Tetrads:-CanonicalForm) canonical form is not implemented for flat or conformally flat spacetimes of Petrov type "O"

     

    NULL

     

    Download TetradsAndWeylScalarsInCanonicalForm.mw

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

    An update to Maple 2016 is now available. Maple 2016.1 provides:

    • Updated translations for Simplified and Traditional Chinese,  French, Greek, Japanese, Brazilian Portuguese, and Spanish
    • Updates to the new Maple Workbook
    • Enhancements to Maple’s context-sensitive menus
    • A variety of improvements to the math engine and interface

     

    To get this update, use Tools>Check for Updates from within Maple, or visit the Maple 2016.1 downloads page.

     

    eithne

     

    MapleSim 2016 is here!

    MapleSim 2016 provides variety of improvements to streamline the user experience, expand modeling scope, and enhance connectivity with other tools. Here are some highlights:

    • Collapsible task panes provide a larger model workspace, so you can see more of your model at once.
    • Improved layout ensures the tools you need for your current task are available at your fingertips.
    • The expanded Multibody component library now supports contact modeling.
    • A new add-on library, the MapleSim Pneumatics Library from Modelon, supports the modeling and simulation of pneumatic systems.
    • The MapleSim CAD Toolbox has been extended to support the latest versions of Inventor®, NX®, SOLIDWORKS®, CATIA® V5, Solid Edge®, PTC® Creo Parametric™, and more.
    • The MapleSim Connector, which provides connectivity to Simulink®, now supports single precision export of S-functions so you can run your MapleSim models on hardware that only supports single precision.

    See What’s New in MapleSim 2016 for more information about these and other improvements.

     

    eithne

    New color schemes for plotting have been added to Maple 2016. They are summarized here. You can also see more details on the plot/colorscheme help pages, which have been reworked (hopefully in a way that makes them more useful to you).

    We released a new video a few weeks ago describing one of these features: coloring by value. The worksheet I used for the video is available here: ColorByValueWebinar.mw

     

    According to Sphere Packing Solved in Higher Dimensions, the best way, i.e., most compact way, to pack spheres in dimensions 8 and 24 are done with the E8 lattice and Leech lattice, respectively. According to the Wikipedia article Leech lattice, the number of spheres that can be packed around any one sphere is 240 and 196,560 (!), respectively, the latter number of spheres counter-intuitively large. It inspired me to try to check that there is indeed room in these lattices for (at least) this number of spheres.

    Starting with the E8 lattice: It is generated by the sum (over the integers) of all the 240 roots of E8. Following the prescription given in the subsection 'Construction' in the Wikipedia article E8 (mathematics), these roots may be constructed as follows:

    ROOTS := map(x -> Vector(x),[
       # Coordinates all integers: 112 roots
       combinat[permute]([+1,+1,0,0,0,0,0,0])[],
       combinat[permute]([+1,-1,0,0,0,0,0,0])[],
       combinat[permute]([-1,-1,0,0,0,0,0,0])[],
       # Coordinates all half-integers: 128 roots
       seq(combinat[permute]([
          (+1/2)$(  2*n),   # Even number of +1/2
          (-1/2)$(8-2*n)    # Even number of -1/2
       ])[],n = 0..4)
    ]):

    This Maple code gives a list of 240 eight-dimensional vectors. All these roots have the same length (the lattice thus being simply laced):

    convert(map(x -> Norm(x,2),ROOTS),set)[];

    If the distance between any pair of different roots is at least this length, then there will be room for 240 spheres of radius equal to this length around any one single sphere. And that is indeed the case:

    DIST_ROOTS := Matrix(nops(ROOTS)$2,(i,j) ->
       Norm(ROOTS[i] - ROOTS[j],2)
    ,shape = symmetric):
    min(convert(DIST_ROOTS,set) minus {0});

    Using the above method for the Leech lattice will fail on grounds of hopeless performance, not the least because DIST_ROOTS will take ages to calculate, if at all possible. So any reader is welcome to weigh in with ideas on how to check the Leech lattice case.

    PS: By the way, I was surprised to find that the three exceptional Lie algebras E6, E7, and E8 are seemingly not accessible through the Maple command SimpleLieAlgebraData, see its help page. Only the four infinite families A,B,C,D, as well as the two exceptional Lie algebras G2 and F4 are. Using Maple 17, I would like to know if that has been changed in Maple 17+, and if not, why not.

    In a recent post, the following inequality was proved with Maple:



    (a,b,c,d >= 0).

    Here is another direct proof attempt.

    f:=(a+b+c+d)^2*(a*b+a*c+a*d+b*c+b*d+c*d)^2-144*(a^2+b^2+c^2+d^2)*a*b*c*d:
    g:=expand(eval(f,d=1)):
    s:=minimize(g, [a=0..infinity,b=0..infinity,c=0..infinity]):
    length(s);   # huge
            304856
    map(evalf@evalf[500],s);

    map(evalf@evalf[1000],s);


    So, Maple returns the expression

    min(0,r1,r2,r3)


    where r1,r2,r3 are huge expressions containing RootOfs. In order to evaluate them, several hundreds of digits are needed.
    The solution seems to be correct, but the question is: may we (mathematically) accept it? What do you think?

     

    I’m pleased to announce the release of Maple T.A. 2016, our online assessment system.

    For this release, we put a lot of effort into streamlining the authoring experience. We worked closely with customers to find out how they authored content, the places where they found the interface awkward, the tasks that took longer than they should have, and what they’d like to see changed. Then we made it better.

    Right away you’ll notice that questions and assignments are no longer in separate places in Maple T.A. All your content is stored in a convenient location that makes it simple to browse your content. Contextual navigation, filtering options, sorting tools, question details, drag and drop organization, combined import feature, and more make it easier than ever to find and organize your content. The Maple T.A. Cloud also sees improvements. Not only can questions be shared, but assignments and entire course modules can be as well.

    For question creation, we consolidated all question authoring into the question designer, so you have a single starting point no matter what kind of question you want to create. We have also refined the text editor to help authors find the tools they need to modify their questions. This includes embedding external content, importing questions from the repository, text formatting options, and more.

    Of course, once you have questions, you’ll want to put them into an assignment, and assignment creation is now easier than ever. A key change is that you can now create and modify questions while you are creating an assignment, without having to leave the assignment editor. There are also changes to how you preview questions, set properties, and even save your assignments, all of which contribute to making assignment creation simpler and faster.

    Of course, there’s more than just a significantly improved author workflow. Here are some highlights:

    • Assignment groups for efficient organization, both in the content repository and on the class homepage.
    • Easy-to-create sorting questions – no coding required!
    • HTML questions, which can be authored directly in the question designer.
    • Clickable image questions are Java-free and easier to author.
    • Maximum word counts and other improvements to the essay question type.
    • A new scanned document feature lets instructors upload and even grade scanned documents.
    • Officially certified LTI integration for connectivity with a wide range of course management systems

    See What’s New in Maple T.A. 2016 for more information on these and other new features.

    Jonny Zivku
    Maplesoft Product Manager, Online Education Products

     

    How to prove the inequality 12*sqrt((a^2+b^2+c^2+d^2)*a*b*c*d) <= (a+b+c+d)*(a*b+a*c+a*d+b*c+b*d+c*d) , assuming that the  variables are nonnegative? That hard question  was asked by arqady in dxdy and answered  by himself  in a complicated way. Maple proves the inequality by the LagrangeMultipliers command which is strong. I think these calculations cannot be done by hand at all. Without loss of generality one may assume a+b+c+d = 1. Then

     restart:with(Student[MultivariateCalculus]):

    ans := [LagrangeMultipliers((a+b+c+d)*(a*b+a*c+a*d+b*c+b*d+c*d)-12*sqrt((a^2+b^2+c^2+d^2)*a*b*c*d), [a+b+c+d-1], [a, b, c, d], output = detailed)]:

    We have to remove complex solutions by
    ans1:=remove(c -> has(evalf(c), I),ans):

    The next big output is  only partly seen in the post (look in the attached file for the whole one).

    ans2:=simplify(ans1,radical);

    [[a = 1/6, b = 1/2, c = 1/6, d = 1/6, lambda[1] = 0, -12*sqrt((a^2+b^2+c^2+d^2)*a*b*c*d)+(b+c+d)*a^2+(b^2+(3*c+3*d)*b+c^2+3*c*d+d^2)*a+(d+c)*b^2+(c^2+3*c*d+d^2)*b+c^2*d+c*d^2 = 0],[a = 1/4, b = 1/4, c = 1/4, d = 1/4, lambda[1] = 0, -12*sqrt((a^2+b^2+c^2+d^2)*a*b*c*d)+(b+c+d)*a^2+(b^2+(3*c+3*d)*b+c^2+3*c*d+d^2)*a+(d+c)*b^2+(c^2+3*c*d+d^2)*b+c^2*d+c*d^2 = 0],[a = 13/72-(1/216)*sqrt(3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))+(1/216)*sqrt(3)*sqrt(2)*sqrt((11*(11548+108*sqrt(13397))^(1/3)-(11548+108*sqrt(13397))^(2/3)+284)/(11548+108*sqrt(13397))^(1/3)+273*sqrt(3)/sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))), b = 11/24+(1/72)*sqrt(3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))-(1/72)*sqrt(3)*sqrt(2)*sqrt((11*(11548+108*sqrt(13397))^(1/3)-(11548+108*sqrt(13397))^(2/3)+284)/(11548+108*sqrt(13397))^(1/3)+273*sqrt(3)/sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))), c = 13/72-(1/216)*sqrt(3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))+(1/216)*sqrt(3)*sqrt(2)*sqrt((11*(11548+108*sqrt(13397))^(1/3)-(11548+108*sqrt(13397))^(2/3)+284)/(11548+108*sqrt(13397))^(1/3)+273*sqrt(3)/sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))), d = 13/72-(1/216)*sqrt(3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))+(1/216)*sqrt(3)*sqrt(2)*sqrt((11*(11548+108*sqrt(13397))^(1/3)-(11548+108*sqrt(13397))^(2/3)+284)/(11548+108*sqrt(13397))^(1/3)+273*sqrt(3)/sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))), lambda[1] = -(5/36)*(sqrt(2)*(sqrt(3)*(sqrt(13397)-(71/27)*(11548+108*sqrt(13397))^(1/3)-(103/540)*(11548+108*sqrt(13397))^(2/3)+2887/27)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))-15*sqrt(13397)+(355/9)*(11548+108*sqrt(13397))^(1/3)+(109/36)*(11548+108*sqrt(13397))^(2/3)-14435/9)*sqrt((11*(11548+108*sqrt(13397))^(1/3)-(11548+108*sqrt(13397))^(2/3)+284)/(11548+108*sqrt(13397))^(1/3)+273*sqrt(3)/sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3)))-(133/15)*(11548+108*sqrt(13397))^(2/3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))+(2*((sqrt(13397)+2374/45)*(11548+108*sqrt(13397))^(1/3)+(103/5)*sqrt(13397)+(449/90)*(11548+108*sqrt(13397))^(2/3)+132727/45))*sqrt(3))/((11548+108*sqrt(13397))^(2/3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))), -12*sqrt((a^2+b^2+c^2+d^2)*a*b*c*d)+(b+c+d)*a^2+(b^2+(3*c+3*d)*b+c^2+3*c*d+d^2)*a+(d+c)*b^2+(c^2+3*c*d+d^2)*b+c^2*d+c*d^2 = -(13/46656)*(((2/13)*sqrt((11*(11548+108*sqrt(13397))^(1/3)-(11548+108*sqrt(13397))^(2/3)+284)/(11548+108*sqrt(13397))^(1/3)+273*sqrt(3)/sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3)))*(11548+108*sqrt(13397))^(1/3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))+sqrt(2)*(sqrt(3)*(11548+108*sqrt(13397))^(1/3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))-(11/13)*(11548+108*sqrt(13397))^(1/3)-(2/13)*(11548+108*sqrt(13397))^(2/3)+568/13))*sqrt(5)*sqrt((sqrt(3)*sqrt(2)*sqrt((11*(11548+108*sqrt(13397))^(1/3)-(11548+108*sqrt(13397))^(2/3)+284)/(11548+108*sqrt(13397))^(1/3)+273*sqrt(3)/sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3)))-sqrt(3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))-33)*(sqrt(3)*sqrt(2)*sqrt((11*(11548+108*sqrt(13397))^(1/3)-(11548+108*sqrt(13397))^(2/3)+284)/(11548+108*sqrt(13397))^(1/3)+273*sqrt(3)/sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3)))-sqrt(3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))+39)*(sqrt(2)*(sqrt(3)*(11548+108*sqrt(13397))^(1/3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))+(11/5)*(11548+108*sqrt(13397))^(1/3)+(2/5)*(11548+108*sqrt(13397))^(2/3)-568/5)*sqrt((11*(11548+108*sqrt(13397))^(1/3)-(11548+108*sqrt(13397))^(2/3)+284)/(11548+108*sqrt(13397))^(1/3)+273*sqrt(3)/sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3)))-(216/5)*(11548+108*sqrt(13397))^(1/3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))-(328/5*((11548+108*sqrt(13397))^(1/3)+(5/164)*(11548+108*sqrt(13397))^(2/3)-355/41))*sqrt(3))/((11548+108*sqrt(13397))^(1/3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))))-(180/13)*sqrt(2)*(sqrt(3)*(11548+108*sqrt(13397))^(1/3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))+(11/5)*(11548+108*sqrt(13397))^(1/3)+(2/5)*(11548+108*sqrt(13397))^(2/3)-568/5)*sqrt((11*(11548+108*sqrt(13397))^(1/3)-(11548+108*sqrt(13397))^(2/3)+284)/(11548+108*sqrt(13397))^(1/3)+273*sqrt(3)/sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3)))-(15552/13)*(11548+108*sqrt(13397))^(1/3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3))+(11808/13*((11548+108*sqrt(13397))^(1/3)+(5/164)*(11548+108*sqrt(13397))^(2/3)-355/41))*sqrt(3))/((11548+108*sqrt(13397))^(1/3)*sqrt((2*(11548+108*sqrt(13397))^(2/3)+11*(11548+108*sqrt(13397))^(1/3)-568)/(11548+108*sqrt(13397))^(1/3)))]

    (1)

    evalf(ans2);

    [[a = .1666666667, b = .5000000000, c = .1666666667, d = .1666666667, lambda[1] = 0., -12.*((a^2+b^2+c^2+d^2)*a*b*c*d)^(1/2)+(b+c+d)*a^2+(b^2+(3.*c+3.*d)*b+c^2+3.*c*d+d^2)*a+(d+c)*b^2+(c^2+3.*c*d+d^2)*b+c^2*d+c*d^2 = 0.], [a = .2500000000, b = .2500000000, c = .2500000000, d = .2500000000, lambda[1] = 0., -12.*((a^2+b^2+c^2+d^2)*a*b*c*d)^(1/2)+(b+c+d)*a^2+(b^2+(3.*c+3.*d)*b+c^2+3.*c*d+d^2)*a+(d+c)*b^2+(c^2+3.*c*d+d^2)*b+c^2*d+c*d^2 = 0.], [a = .1666666667, b = .1666666667, c = .5000000000, d = .1666666667, lambda[1] = 0., -12.*((a^2+b^2+c^2+d^2)*a*b*c*d)^(1/2)+(b+c+d)*a^2+(b^2+(3.*c+3.*d)*b+c^2+3.*c*d+d^2)*a+(d+c)*b^2+(c^2+3.*c*d+d^2)*b+c^2*d+c*d^2 = 0.], [a = .1666666667, b = .1666666667, c = .1666666667, d = .5000000000, lambda[1] = 0., -12.*((a^2+b^2+c^2+d^2)*a*b*c*d)^(1/2)+(b+c+d)*a^2+(b^2+(3.*c+3.*d)*b+c^2+3.*c*d+d^2)*a+(d+c)*b^2+(c^2+3.*c*d+d^2)*b+c^2*d+c*d^2 = 0.], [a = .5000000000, b = .1666666667, c = .1666666667, d = .1666666667, lambda[1] = 0., -12.*((a^2+b^2+c^2+d^2)*a*b*c*d)^(1/2)+(b+c+d)*a^2+(b^2+(3.*c+3.*d)*b+c^2+3.*c*d+d^2)*a+(d+c)*b^2+(c^2+3.*c*d+d^2)*b+c^2*d+c*d^2 = 0.], [a = .2118620934, b = .3644137199, c = .2118620934, d = .2118620934, lambda[1] = 0.2834790478e-2, -12.*((a^2+b^2+c^2+d^2)*a*b*c*d)^(1/2)+(b+c+d)*a^2+(b^2+(3.*c+3.*d)*b+c^2+3.*c*d+d^2)*a+(d+c)*b^2+(c^2+3.*c*d+d^2)*b+c^2*d+c*d^2 = 0.9449303017e-3], [a = 0.3692850681e-1, b = .8892144797, c = 0.3692850681e-1, d = 0.3692850681e-1, lambda[1] = 0.9303874297e-1, -12.*((a^2+b^2+c^2+d^2)*a*b*c*d)^(1/2)+(b+c+d)*a^2+(b^2+(3.*c+3.*d)*b+c^2+3.*c*d+d^2)*a+(d+c)*b^2+(c^2+3.*c*d+d^2)*b+c^2*d+c*d^2 = 0.3101291407e-1], [a = .8892144797, b = 0.3692850681e-1, c = 0.3692850681e-1, d = 0.3692850681e-1, lambda[1] = 0.9303874297e-1, -12.*((a^2+b^2+c^2+d^2)*a*b*c*d)^(1/2)+(b+c+d)*a^2+(b^2+(3.*c+3.*d)*b+c^2+3.*c*d+d^2)*a+(d+c)*b^2+(c^2+3.*c*d+d^2)*b+c^2*d+c*d^2 = 0.3101291407e-1], [a = .3644137199, b = .2118620934, c = .2118620934, d = .2118620934, lambda[1] = 0.2834790478e-2, -12.*((a^2+b^2+c^2+d^2)*a*b*c*d)^(1/2)+(b+c+d)*a^2+(b^2+(3.*c+3.*d)*b+c^2+3.*c*d+d^2)*a+(d+c)*b^2+(c^2+3.*c*d+d^2)*b+c^2*d+c*d^2 = 0.9449303017e-3], [a = 0.3692850681e-1, b = 0.3692850681e-1, c = 0.3692850681e-1, d = .8892144797, lambda[1] = 0.9303874297e-1, -12.*((a^2+b^2+c^2+d^2)*a*b*c*d)^(1/2)+(b+c+d)*a^2+(b^2+(3.*c+3.*d)*b+c^2+3.*c*d+d^2)*a+(d+c)*b^2+(c^2+3.*c*d+d^2)*b+c^2*d+c*d^2 = 0.3101291407e-1], [a = .2118620934, b = .2118620934, c = .2118620934, d = .3644137199, lambda[1] = 0.2834790478e-2, -12.*((a^2+b^2+c^2+d^2)*a*b*c*d)^(1/2)+(b+c+d)*a^2+(b^2+(3.*c+3.*d)*b+c^2+3.*c*d+d^2)*a+(d+c)*b^2+(c^2+3.*c*d+d^2)*b+c^2*d+c*d^2 = 0.9449303017e-3], [a = 0.3692850681e-1, b = 0.3692850681e-1, c = .8892144797, d = 0.3692850681e-1, lambda[1] = 0.9303874297e-1, -12.*((a^2+b^2+c^2+d^2)*a*b*c*d)^(1/2)+(b+c+d)*a^2+(b^2+(3.*c+3.*d)*b+c^2+3.*c*d+d^2)*a+(d+c)*b^2+(c^2+3.*c*d+d^2)*b+c^2*d+c*d^2 = 0.3101291407e-1], [a = .2118620934, b = .2118620934, c = .3644137199, d = .2118620934, lambda[1] = 0.2834790478e-2, -12.*((a^2+b^2+c^2+d^2)*a*b*c*d)^(1/2)+(b+c+d)*a^2+(b^2+(3.*c+3.*d)*b+c^2+3.*c*d+d^2)*a+(d+c)*b^2+(c^2+3.*c*d+d^2)*b+c^2*d+c*d^2 = 0.9449303017e-3]]

    (2)

    Indeed, the minimum value of the target function is exactly 0. Quod erat demonstrantum.

    NULL

     inequality.mw

     

    A wealth of knowledge is on display in MaplePrimes as our contributors share their expertise and step up to answer others’ queries. This post picks out one such response and further elucidates the answers to the posted question. I hope these explanations appeal to those of our readers who might not be familiar with the techniques embedded in the original responses.

    The Question: Transforming functions to names

    Bendesarts wanted to know how to make programmatic changes to characters in a list. He wrote:

    I have this list :

    T:=[alpha(t),beta(t)]

    I would like to create this list automatically:

    Tmod:=[alpha_,beta_]

    In other words, how can I remove the 3 characters "(t)" and replace it by "_"

    Do you have ideas to do so ?

    Thanks a lot for your help

    Joe Riel provided a complete answer that had three different approaches. He wrote:

    Map, over the list, a procedure that extracts the name of the function and catenates an underscore to it. The function name is extracted via the op procedure, e.g. op(0,f(a,b,c)) evaluates to f. Thus 

    map(f->cat(op(0,f),_),T);

    Note that this fails if the name of a function is indexed, e.g. f[2](a). No error is generated but the catenation doesn't evaluate to a symbol. Presumably that isn't an issue here.  One way to handle that case is to first convert the indexed name to a symbol, then catenate the underscore.  So a more robust version is

    map(f->cat(convert(op(0,f),'symbol'),_),T);

    However, if you are actually dealing with indexed names you might want a different result. Another way to do the conversion, and combine it with the catenation, is to use nprintf, which generates a name (symbol). Thus

    map(f -> nprintf("%a_", op(0,f)),T);

     

    Let’s discuss each approach by understanding the definitions and functionalities of the commands used. 

    The map command, map(fcn, expr, arg1, ..., argN) applies a procedure or name, fcn, to the operands or elements of an expression, expr. The result of a call to map is a copy of expr with the ith operand of expr replaced by the result of applying fcn to the ith operand.  This concept is easier to grasp by looking at a few examples related to the usage of map in this question.

    Example 1.  map(x-> x2,x+y)         returns     x2+y2                    

    Example 2. map(a -> a-b, sin(x))    returns     sin(x-b)

     

    The cat function, cat(a,b,c,…), is commonly used to concatenate (or join) string and names together. This function’s parameters: a,b,c…, can be any expressions.

    Example 1. cat(a,2)                      returns     a2

    Example 2.  cat(“a”,3,sin(x))          returns    “a3sin(x)”

     

    The op function, op(i..j,e), extracts operands from an expression. The parameters i and j are the integers indicating positions of the operands and e is the expression. For functions, as in this example, op(0,e) is the name of the function.

    Example 1.  op(0,alpha(t))            returns   the symbol alpha

    Example 2.  op(0, sin(x))              returns    sin

     

    Now analyzing Joe Riel's code will be easier.

    1. map(f->cat(op(0,f),_),T);

    In this approach Joe is extracting the name of the functions, alpha and beta, and then concatenating it to the underscore symbol. Then using the mapping function he applies the previous procedure to the list T.

    1. map(f->cat(convert(op(0,f),'symbol'),_),T);

    This approach is a lot similar to the previous one, but he added the convert function in case the function inside of map was indexed. Convert(expr, form, arg3,..), is used to change an expression from one form to another. In this example op(0,f) has been changed from type name to type symbol.

    1. map(f -> nprintf("%a_", op(0,f)),T);

    Again this is a similar approach but it uses nprintf. This command, nprintf(fmt,x1,..xn), is based on a C standard library command of the same name. It uses the format specifications in the fmt string to format and writes the expression into a Maple symbol, which is returned. In this example the format specified is the algebraic format “%a”.

     

    This blog was written by Maplesoft’s intern Pia under the supervision of Dr. Robert Lopez. We both hope that you find this useful. If there is a particular question on MaplePrimes that you would like further explained, please let us know. 

    The mechanism of transport of the material of the sewing machine M 1022 class: mathematical animation.   BELORUS.mw 




     

    How to prove the inequality x^(4*y)+y^(4*x) <= 2 provided x^2+y^2 = 2, 0 <= x, 0 <= y? That problem was posed  by Israeli mathematician nicked by himself as arqady in Russian math forum and was not answered there.I know how to prove that with Maple and don't know how to prove that without Maple. Neither LagrangeMultipliers nor extrema work here. The difficulty consists in the nonlinearity both the target function and the main constraint. The first step is to linearize the main constraint and the second step is to reduce the number of variables to one.

    restart; A := eval(x^(4*y)+y^(4*x), [x = sqrt(u), y = sqrt(v)]);

    (u^(1/2))^(4*v^(1/2))+(v^(1/2))^(4*u^(1/2))

    (1)

     

    B := expand(A);

    u^(2*v^(1/2))+v^(2*u^(1/2))

    (2)

    C := eval(B, u = 2-v);

    (2-v)^(2*v^(1/2))+v^(2*(2-v)^(1/2))

    (3)

    It is more or less clear that the plot of F is symmetric wrt  the straight line v=1. This motivates the following change of variable  to obtain an even function.

    F := simplify(expand(eval(C, v = z+1)), symbolic, power);

    (1-z)^(2*(z+1)^(1/2))+(z+1)^(2*(1-z)^(1/2))

    (4)

    NULL

    The plots suggest the only maximim of F at z=0 and its concavity.

    Student[Calculus1]:-FunctionPlot(F, z = -1 .. 1);

     

    Student[Calculus1]:-FunctionPlot(diff(F, z, z), z = -1 .. 1);

     

    As usually, numeric global solvers cannot prove certain inequalities. However, the GlobalSearch command of the DirectSearch package indicates the only local maximum of  F and F''.NULL

    Digits := 25; DirectSearch:-GlobalSearch(F, {z = -1 .. 1}, maximize, solutions = 3, tolerances = 10^(-15)); DirectSearch:-GlobalSearch(diff(F, z, z), {z = -1 .. 1}, maximize, solutions = 3, tolerances = 10^(-15));

    Array([[0.8e-23, [z = -0.1980181305884928531875965e-12], 36]])

    (5)

    The series command confirms a local maximum of F at z=0.

    series(F, z, 6);

    series(2-(2/3)*z^4+O(z^6),z,6)

    (6)

    The extrema command indicates only the value of F at a critical point, not outputting its position.

    extrema(F, z); extrema(F, z, 's');

    {2}

    (7)

    solve(F = 2);

    RootOf((1-_Z)^(2*(_Z+1)^(1/2))+(_Z+1)^(2*(1-_Z)^(1/2))-2)

    (8)

    DirectSearch:-SolveEquations(F = 2, {z = -1 .. 1}, AllSolutions, solutions = 3);

    Matrix(1, 4, {(1, 1) = 0., (1, 2) = Vector(1, {(1) = 0.}), (1, 3) = [z = -0.5463886313e-6], (1, 4) = 27})

    (9)

    DirectSearch:-SolveEquations(F = 2, {z = -1 .. 1}, AllSolutions, solutions = 3, assume = integer);

    Matrix(1, 4, {(1, 1) = 0., (1, 2) = Vector(1, {(1) = 0.}), (1, 3) = [z = 0], (1, 4) = 30})

    (10)

    NULL

     PS. I see my proof needs an additional explanation. The DirectSearch command establishes the only both local and global  maximum of F is located at z= -1.98*10^(-13) up to default error 10^(-9). After that  the series command confirms a local maximum at z=0. Combining these, one draws the conclusion that the global maximum is placed exactly at z=0 and equals 2. In order to confirm that the only real root of F=2 at z=0  is found approximately and exactly by the DirectSearch.

    Download maxi.mw

    I thought I would share some code for computing sparse matrix products in Maple.  For floating point matrices this is done quickly, but for algebraic datatypes there is a performance problem with the builtin routines, as noted in http://www.mapleprimes.com/questions/205739-How-Do-I-Improve-The-Performance-Of

    Download spm.txt

    The code is fairly straightforward in that it uses op(1,A) to extract the dimensions and op(2,A) to extract the non-zero elements of a Matrix or Vector, and then loops over those elements.  I included a sparse map function for cases where you want to map a function (like expand) over non-zero elements only.

    # sparse matrix vector product
    spmv := proc(A::Matrix,V::Vector)
    local m,n,Ae,Ve,Vi,R,e;
    n, m := op(1,A);
    if op(1,V) <> m then error "incompatible dimensions"; end if;
    Ae := op(2,A);
    Ve := op(2,V);
    Vi := map2(op,1,Ve);
    R := Vector(n, storage=sparse);
    for e in Ae do
    n, m := op(1,e);
    if member(m, Vi) then R[n] := R[n] + A[n,m]*V[m]; end if;
    end do;
    return R;
    end proc:
    # sparse matrix product
    spmm := proc(A::Matrix, B::Matrix)
    local m,n,Ae,Be,Bi,R,l,e,i;
    n, m := op(1,A);
    i, l := op(1,B);
    if i <> m then error "incompatible dimensions"; end if;
    Ae := op(2,A);
    Be := op(2,B);
    R := Matrix(n,l,storage=sparse);
    for i from 1 to l do
    Bi, Be := selectremove(type, Be, (anything,i)=anything);
    Bi := map2(op,[1,1],Bi);
    for e in Ae do
    n, m := op(1,e);
    if member(m, Bi) then R[n,i] := R[n,i] + A[n,m]*B[m,i]; end if;
    end do;
    end do;
    return R;
    end proc:
    # sparse map
    smap := proc(f, A::{Matrix,Vector})
    local B, Ae, e;
    if A::Vector then
    B := Vector(op(1,A),storage=sparse):
    else
    B := Matrix(op(1,A),storage=sparse):
    end if;
    Ae := op(2,A);
    for e in Ae do
    B[op(1,e)] := f(op(2,e),args[3..nargs]);
    end do;
    return B;
    end proc:


    As for how it performs, here is a demo inspired by the original post.

    n := 674;
    k := 6;
    A := Matrix(n,n,storage=sparse):
    for i to n do
      for j to k do
        A[i,irem(rand(),n)+1] := randpoly(x):
      end do:
    end do:
    V := Vector(n):
    for i to k do
      V[irem(rand(),n)+1] := randpoly(x):
    end do:
    C := CodeTools:-Usage( spmv(A,V) ):  # 7ms, 25x faster
    CodeTools:-Usage( A.V ):  # 174 ms
    B := Matrix(n,n,storage=sparse):
    for i to n do
      for j to k do
        B[i,irem(rand(),n)+1] := randpoly(x):
      end do:
    end do:
    C := CodeTools:-Usage( spmm(A,B) ):  # 2.74 sec, 50x faster
    CodeTools:-Usage( A.B ):  # 2.44 min
    # expand and collect like terms
    C := CodeTools:-Usage( smap(expand, C) ):

    CameraProfiler.mw

    D700_Profile.xlsx

    This Application uses the xrite color checker to obtain Forward and Color Matrices as defined by Adobe.

    It compares the camera gamut based on the above matrices to the 1931 Standard Observer using LEDs.

     

     

     

      Continuation.
      One way to get rolling without slipping animation in 3d. The trajectory and circle are divided into segments of equal length. In the next segment of the trajectory we construct circle, taking into account the fact that it turned on one segment. Rolling sphere or cylinder can be simulated, if we take plottools templates of the same radius, and replace them on the site of our circle.

    ROLLING_WITHOUT_3d.mw













    First 53 54 55 56 57 58 59 Last Page 55 of 291