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

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

         (1)

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

         (2)

    with the exact derivative of the temperature given by

         (3)

     

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

         (4)

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

         (5)

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

     

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

         (6)

    with the second derivative of the temperature given as

         (7)

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

         (8)

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

         (9)

    Defining the operator

         (10)

    we can now assemble into a fully populated matrix as,

         (11)

     

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

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

     

    The Maple code listing follows:

    > restart:
       with(LinearAlgebra):with(plots):
    
    # MESHLESS METHOD SOLUTION USING MULTIQUADRIC RADIAL BASIS
    FUNCTIONS (RBF) il:=6:To:=15:TL:=25:L:=1:
    >   x:=[0,1/5,2/5,3/5,4/5,1]:
    >   S:=1000:n:=1:dx:=1/(il-1):
    >   C:=Array(1..il,1..il):phi:=Array(1..il,1..il):d2phi:=Array(1..il,1..il):
    b:=Vector(1..il):TM:=Vector(1..il):alpha:=Vector(1..il):
    for i from 1 to il do
       for j from 1 to il do
          phi[i,j]:=(1+(x[i]-x[j])^2/(S*dx^2))^(n-3/2):
          d2phi[i,j]:=3*((x[j]-x[i])/20)^2/(4*((x[j]-x[i])^2/40+1)^(5/2) )-1/(40*((x[j]-x[i])^2/40+1)^(3/2)):
       end do:
    end do:
    >   for i from 2 to il-1 do    
    >       for j from 1 to il do
    >         C[i,j]:=d2phi[i,j]+phi[i,j];
                b[i]:=-x[i];
    >         C[1,j]:=phi[1,j];
               C[il,j]:=phi[il,j];
             end do:
          end do:
          b[1]:=To:b[il]:=TL:
    > #ConditionNumber(C);
    >   alpha:=LinearSolve(convert(C,Matrix),b):
    TM[1]:=To:TM[6]:=TL:
    for i from 2 to il-1 do
       for j from 1 to il do
    >         TM[i]:=TM[i]+alpha[j]*(1+(x[i]-x[j])^2/(S*dx^2))^(n-3/2);
     >     end do:
    >   end do:
         evalf(TM);
    
                         
    
    > TE:=To*cos(xx)+(TL+L-To*cos(L))/sin(L)*sin(xx)-xx:
    > TE:=subs(TE):
    > TE:=plot(TE,xx=0..1,color=blue,legend="Exact",thickness=3):
    > MEM:=[seq([subs(x[i]),subs(TM[i])],i=1..6)]:
    > T:=plots[pointplot](MEM,style=line,color=red,legend="MEM", thickness=3): 
      MEM:=plots[pointplot](MEM,color=red,legend="MEM",symbol=box, symbolsize=15):
    
    >  plots[display](TE,MEM,T,axes=BOXED,title="Solution - MEM");
    
                  
    
    

     

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

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

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

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

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

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

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

    Happy gift giving!

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

     

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

    with(Physics)

     

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

    Setup(signature)

    [signature = `- - - +`]

    (1)

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

    Library:-SignOfTimelikeComponent()

    1

    (2)

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

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

    -1

    (3)

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

     

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

     

    Setup(Dgamma = standard, math = true)

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

     

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

     

    `Setting lowercaselatin letters to represent spinor indices `

     

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

     

    __________________________________________________

     

    [Dgammarepresentation = standard, mathematicalnotation = true]

    (4)

    The four Dirac matrices are

    TensorArray(Dgamma[`~mu`])

    Array(%id = 18446744078360533342)

    (5)

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

    "Library:-RewriteInMatrixForm(?)"

    Array(%id = 18446744078360529726)

    (6)

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

    TensorArray(Dgamma[mu])

    Array(%id = 18446744078565663678)

    (7)

    "Library:-RewriteInMatrixForm(?)"

    Array(%id = 18446744078677131982)

    (8)

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

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

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

    (9)

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

     

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

     

    Verify the algebra rule by performing all the involved matrix operations

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

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

    (10)

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

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

    Matrix(%id = 18446744078587020822)

    (11)

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

    "Library:-PerformMatrixOperations(?)"

    Matrix(%id = 18446744078743243942)

    (12)

    By eye everything checks OK.NULL

     

    Consider now the following five products of Dirac matrices

    e0 := Dgamma[mu]^2

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

    (13)

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

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

    (14)

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

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

    (15)

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

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

    (16)

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

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

    (17)

    New: the simplification of these products is now implemented

    e0 = Simplify(e0)

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

    (18)

    Verify this result performing the underlying matrix operations

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

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

    (19)

    Library:-PerformMatrixOperations(T)

    Matrix(%id = 18446744078553169662) = 4

    (20)

    The same with the other expressions

    e1 = Simplify(e1)

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

    (21)

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

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

    (22)

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

    Array(%id = 18446744078695012102)

    (23)

    Library:-PerformMatrixOperations(T)

    Array(%id = 18446744078701714238)

    (24)

    For e2

    e2 = Simplify(e2)

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

    (25)

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

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

    (26)

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

    Matrix(%id = 18446744078470204942)

    (27)

    Library:-PerformMatrixOperations(T)

    Matrix(%id = 18446744078550068870)

    (28)

    For e3 we have

    e3 = Simplify(e3)

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

    (29)

    Verify this result,

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

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

    (30)

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

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

    Array(%id = 18446744078647326830)

    (31)

    For instance, the first element is

    T[1, 1, 1]

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

    (32)

    and it checks OK:

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

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

    (33)

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

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

     

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

     

    For instance,

    M[1]

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

    (34)

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

     

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

    {true}

    (35)

    The same for e4

    e4 = Simplify(e4)

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

    (36)

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

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

    (37)

    Regarding the spacetime indices this is now an array 4x4x4x4

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

    Array(%id = 18446744078730196382)

    (38)

    For instance the first of these 256 matrix equations

    T[1, 1, 1, 1]

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

    (39)

    verifies OK:

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

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

    (40)

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

     

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

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

    {true}

    (41)

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

     

    For example,

    Define(A, B)

    `Defined as tensors`

     

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

    (42)

    Contract with e1 and e2 and simplify

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

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

    (43)

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

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

    (44)

     


     

    Download DiracMatricesAndPerformMatrixOperation.mw

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

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

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

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

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

    Application to generate vector exercises using very easy to use patterns. These exercises are for sum of vectors in calculations of magnitude, direction, graph and projections. Basically to evaluate our students on the board and thus not lose the results. For university professors and engineering students.

    Generator_of_exercises_with_vectors.mw

    Lenin Araujo Castillo

    Ambassador of Maple

    We added a small new feature to MaplePrimes that will make it easier to tag other members in your questions, posts and comments.

    To use the feature, type "@" in a message. If you wait a moment, a list of names will be populated, starting with members who are participating in the current thread, followed by other members ordered by reputation. You can also begin to type after the "@" to search for a particular member. Similar to the above, members who match your search who are in the current thread will appear first, followed by others.

    If you are tagged in this way, an indication will appear in your notifications pane.

    We hope you find this new feature useful!

     

    Many of you enjoyed our profile on one of our developers, Paulina Chin, so we’re happy to bring you another one!

    Today, we’ll be talking with John May, Senior Developer of Maple. Let’s get started.

    1. What do you do at Maplesoft?
      Until recently I was consulting on-site at the NASA Jet Propulsion Laboratory helping people there more effectively solve their engineering problems using Maplesoft products.  But my main job that I am back to full time now is the development and maintenance of various parts of the Maple library.
       
    2. What did you study in school?
      I studied both Pure and Applied Mathematics at the University of Oregon,  focusing a lot on Abstract Algebra.  In graduate school, I specialized more in computation mathematics like computer algebra and numerical analysis.  My Ph.D. work focused on effective numerical algorithms for problems in polynomial algebra – with implementations in Maple!
       
    3. What area(s) of Maple are you currently focusing on in your development?
      Right now I am focused on addressing complaints I’ve gotten from engineers about the usability of units with other parts of the math library.
       
    4. What’s the coolest feature of Maple that you’ve had a hand in developing?
      A lot of the cool things I’ve built live pretty deep in the internals of Maple.  I’ve done a lot of meta-heuristic tuning to seamlessly integrate high-performance libraries into top-level Maple commands.

      I had a lot of fun developing a lot of the stuff for manipulation and visualization of colors in the ColorTools package.
       
    5. What do you like most about working at Maplesoft? How long have you worked here?
      I started working at Maple in 2007, but I’ve been a Maple user since 1997.  I love being part of the magic that brings powerful algorithmic mathematics to everyone.  The R&D team is also full of eccentric nerds who are great fun to work with.
       
    6. Favourite hobby?
      It varies by the season, but right now it is prime for mountain biking in southern California.  I ride my local trails a couple times a week, and when I get I chance, I love to get away on epic bikepacking adventures (like this one: https://www.bikemag.com/features/two-wheeled-escape-one-hour-from-l-a/  this is me: https://cdn.bikemag.com/uploads/2016/05/16File.jpg ).
       
    7. What do you like on your pizza?
      Anything and everything. Something different every time. My all-time favorite pie my from grad school days is the “Rio Rancho” from the dearly departed That’s Amore Pizza (which was next to the comic book store and across the street from North Carolina State University).  It was an olive oil and mozzarella pizza with chopped bacon that was covered in sliced fresh roma tomatoes and drizzled with ranch dressing when it came out of the oven. 
       
    8. What’s your favourite movie?
      It’s really hard to pick just one.  So, I’ll go with the safe answer and say the greatest movie of all time, and “Weird Al” Yankovic’s only foray into movies, UHF, is my favorite.
      http://www.imdb.com/title/tt0098546/
       
    9. What skill would you love to learn? (That you haven’t already) Why?
      Another hard one.  I feel like I’ve dabbled in lots of things that I would like to get better at.  At the top of the list is probably unicycling.  I’d love to get good enough to play Unicyle Football or do Muni (mountain unicyling).
      https://en.wikipedia.org/wiki/Mountain_unicycling
      http://www.unicyclefootball.com/
       
    10. Who’s your favourite mathematician?
      Batman. https://youtu.be/AcMEckOyoaM

     

    The Perimeter Institute for Theoretical Physics (PI) is a place where bold ideas flourish. It brings together great minds in a shared effort to achieve scientific breakthroughs that will transform our future. PI is an independent research center in foundational theoretical physics.

    One of the key mission objectives of Perimeter is to provide scientific training and educational outreach activities to the general public. Maplesoft is proud to be part of this endeavor, as PI’s Educational Outreach Champion. Maplesoft’s contributions support Perimeter’s Teacher Network training activities, core educational resources, development of online course material, support of events such as the EinsteinPlus workshop for high school teachers, the International Summer School for Young Physicists (ISSYP), and other initiatives.
     

    The annual International Summer School for Young Physicists (ISSYP) is a two-week camp that brings together 40 exceptional physics-minded students from high schools across the globe to PI. The program covers many different topics in physics such as quantum mechanics, special relativity, general relativity, and cosmology. Each year students receive a complimentary copy of Maple, and use the product to practice and strengthen their math skills. The program receives an average of three-hundred applications from students in grades eleven and twelve from around the world. Competition is intense, and students who are chosen for the program are extremely bright and advanced for their age; however there is some variation in their level of math and physics knowledge. . Students are asked to review a “math primer” document to prepare them with the background needed for the program.The ISSYP program now uses Möbius, Maplesoft’s online courseware platform to administer this primer. With Möbius, PI has moved from a pdf document primer to fully online material, which has motivated more students to complete the material and be more engaged in their courses. The interactivity and engagement that technology provides has made the summer program more productive and dynamic.

    EinsteinPlus is a one-week intensive workshop for Canadian and international high school teachers that focuses on modern physics, including quantum physics, special relativity, and cosmology. EinsteinPlus also provides unique opportunities to learn some of the latest developments in physics from expert researchers at the forefront of their fields. Maplesoft proudly supports this workshop by giving teachers access to and training in Maple.

     

    Perimeter Institute also organizes a lively program of seminars, regularly exposing researchers and students to current ideas in the wider theoretical physics community. The talks provide content outside of, but related to, core disciplines.  Recently Maplesoft’s own physics expert Dr. Edgardo Cheb-Terrab conducted a lecture and training session at PI on Computer Algebra for Theoretical Physics.

    Dr. Edgardo Cheb-Terrab is the force behind the algorithms and Maple libraries of the ODE and PDE symbolic solvers, the MathematicalFunctions package (an expert system in special functions) and the Physics package, among other things.  

    In his talk at PI, the Physics project at Maplesoft was presented and the resulting Physics package was illustrated through simple problems in classical field theory, quantum mechanics and general relativity, and through tackling the computations of some recent Physical Review papers in those areas.   In addition there was a hands-on workshop where attendees were offered four choices of activity:  follow the mini-course; explore items of the worksheet of the morning presentation; or bring their own problems so that Dr. Cheb-Terrab could guide them on how to tackle it using the Physics package and Maple in general.

    As a company that strives to continuously improve student learning, and empower instructors and researchers with the tools necessary to compete in an ever changing and demanding educational environment, Maplesoft’s partnership with the Perimeter Institute allows us to do just that. We take great pride and joy in bringing our technology to outreach programs for students and teachers, making these opportunities a more productive and dynamic experience for all.

     

    In a previous Mapleprimes question related to Dirac Matrices, I was asked how to represent the algebra of Dirac matrices with an identity matrix on the right-hand side of  %AntiCommutator(Physics:-Dgamma[j], Physics:-Dgamma[k]) = 2*g[j, k]. Since this is a hot-topic in general, in that, making it work, involves easy and useful functionality however somewhat hidden, not known in general, it passed through my mind that this may be of interest in general. (To reproduce the computations below you need to update your Physics library with the one distributed at the Maplesoft R&D Physics webpage.)

     

    restart

    with(Physics)

     

    First of all, this shows the default algebra rules loaded when you load the Physics package, for the Pauli  and Dirac  matrices

    Library:-DefaultAlgebraRules()

    %Commutator(Physics:-Psigma[j], Physics:-Psigma[k]) = (2*I)*(Physics:-Psigma[1]*Physics:-LeviCivita[4, j, k, `~1`]+Physics:-Psigma[2]*Physics:-LeviCivita[4, j, k, `~2`]+Physics:-Psigma[3]*Physics:-LeviCivita[4, j, k, `~3`]), %AntiCommutator(Physics:-Psigma[j], Physics:-Psigma[k]) = 2*Physics:-KroneckerDelta[j, k], %AntiCommutator(Physics:-Dgamma[j], Physics:-Dgamma[k]) = 2*Physics:-g_[j, k]

    (1)

    Now, you can always overwrite these algebra rules.

     

    For instance, to represent the algebra of Dirac matrices with an identity matrix on the right-hand side, one can proceed as follows.

    First create the identity matrix. To emulate what we do with paper and pencil, where we write I to represent an identity matrix without having to see the actual table 2x2 with the number 1 in the diagonal and a bunch of 0, I will use the old matrix command, not the new Matrix (see more comments on this at the end). One way of entering this identity matrix is

    `𝕀` := matrix(4, 4, proc (i, j) options operator, arrow; KroneckerDelta[i, j] end proc)

    array( 1 .. 4, 1 .. 4, [( 4, 1 ) = (0), ( 1, 2 ) = (0), ( 2, 3 ) = (0), ( 1, 3 ) = (0), ( 2, 2 ) = (1), ( 4, 2 ) = (0), ( 3, 4 ) = (0), ( 1, 4 ) = (0), ( 3, 1 ) = (0), ( 4, 4 ) = (1), ( 3, 2 ) = (0), ( 1, 1 ) = (1), ( 2, 1 ) = (0), ( 4, 3 ) = (0), ( 3, 3 ) = (1), ( 2, 4 ) = (0)  ] )

    (2)

    The most important advantage of the old matrix command is that I is of type algebraic and, consequently, this is the important thing, one can operate with it algebraically and its contents are not displayed:

    type(`𝕀`, algebraic)

    true

    (3)

    `𝕀`

    `𝕀`

    (4)

    And so, most commands of the Maple library, that only work with objects of type algebraic, will handle the task. The contents are displayed only on demand, for instance using eval

    eval(`𝕀`)

    array( 1 .. 4, 1 .. 4, [( 4, 1 ) = (0), ( 1, 2 ) = (0), ( 2, 3 ) = (0), ( 1, 3 ) = (0), ( 2, 2 ) = (1), ( 4, 2 ) = (0), ( 3, 4 ) = (0), ( 1, 4 ) = (0), ( 3, 1 ) = (0), ( 4, 4 ) = (1), ( 3, 2 ) = (0), ( 1, 1 ) = (1), ( 2, 1 ) = (0), ( 4, 3 ) = (0), ( 3, 3 ) = (1), ( 2, 4 ) = (0)  ] )

    (5)

    Returning to the topic at hands: set now the algebra the way you want, with an I matrix on the right-hand side, and without seeing a bunch of 0 and 1

    %AntiCommutator(Dgamma[mu], Dgamma[nu]) = 2*g_[mu, nu]*`𝕀`

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

    (6)

    Setup(algebrarules = (%AntiCommutator(Physics[Dgamma][mu], Physics[Dgamma][nu]) = 2*Physics[g_][mu, nu]*`𝕀`))

    [algebrarules = {%AntiCommutator(Physics:-Dgamma[mu], Physics:-Dgamma[nu]) = 2*Physics:-g_[mu, nu]*`𝕀`}]

    (7)

    And that is all.

     

    Check it out

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

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

    (8)

    Set now a Dirac spinor; this is how you could do that, step-by-step.

     

    Again you can use {vector, matrix, array} or {Vector, Matrix, Array}, and again, if you use the Upper case commands, you always have the components visible, and cannot compute with the object normally since they are not of type algebraic. So I use matrix, not Matrix, and matrix instead of vector so that the Dirac spinor that is both algebraic and matrix, is also displayed in the usual display as a "column vector"

     

    _local(Psi)

    Setup(anticommutativeprefix = {Psi, psi})

    [anticommutativeprefix = {_lambda, psi, :-Psi}]

    (9)

    In addition, following your question, in this example I explicitly specify the components of the spinor, in any preferred way, for example here I use psi[j]

    Psi := matrix(4, 1, [psi[1], psi[2], psi[3], psi[4]])

    array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] )

    (10)

    Check it out:

    Psi

    Psi

    (11)

    type(Psi, algebraic)

    true

    (12)

    Let's see all this working together by multiplying the anticommutator equation by Psi

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

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

    (13)

    Suppose now that you want to see the matrix form of this equation

    Library:-RewriteInMatrixForm(Physics[`*`](%AntiCommutator(Physics[Dgamma][mu], Physics[Dgamma][nu]), Psi) = 2*Physics[g_][mu, nu]*Physics[`*`](`𝕀`, Psi))

    Physics:-`.`(%AntiCommutator(Physics:-Dgamma[mu], Physics:-Dgamma[nu]), array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] )) = 2*Physics:-g_[mu, nu]*Physics:-`.`(array( 1 .. 4, 1 .. 4, [( 4, 1 ) = (0), ( 1, 2 ) = (0), ( 2, 3 ) = (0), ( 1, 3 ) = (0), ( 2, 2 ) = (1), ( 4, 2 ) = (0), ( 3, 4 ) = (0), ( 1, 4 ) = (0), ( 3, 1 ) = (0), ( 4, 4 ) = (1), ( 3, 2 ) = (0), ( 1, 1 ) = (1), ( 2, 1 ) = (0), ( 4, 3 ) = (0), ( 3, 3 ) = (1), ( 2, 4 ) = (0)  ] ), array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] ))

    (14)

    The above has the matricial operations delayed; unleash them

    %

    Physics:-`.`(%AntiCommutator(Physics:-Dgamma[mu], Physics:-Dgamma[nu]), array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] )) = 2*Physics:-g_[mu, nu]*(array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] ))

    (15)

    Or directly perform in one go the matrix operations behind (13)

    Library:-PerformMatrixOperations(Physics[`*`](%AntiCommutator(Physics[Dgamma][mu], Physics[Dgamma][nu]), Psi) = 2*Physics[g_][mu, nu]*Physics[`*`](`𝕀`, Psi))

    Physics:-`.`(%AntiCommutator(Physics:-Dgamma[mu], Physics:-Dgamma[nu]), array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] )) = 2*Physics:-g_[mu, nu]*(array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] ))

    (16)

    REMARK: As shown above, in general, the representation using lowercase commands allows you to use `*` or `.` depending on whether you want to represent the operation or perform the operation. For example this represents the operation, as an exact mimicry of what we do with paper and pencil, both regarding input and output

    `𝕀`*Psi

    Physics:-`*`(`𝕀`, Psi)

    (17)

    And this performs the operation

    `𝕀`.Psi

    array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] )

    (18)

    Or to only displaying the operation

    Library:-RewriteInMatrixForm(Physics[`*`](`𝕀`, Psi))

    Physics:-`.`(array( 1 .. 4, 1 .. 4, [( 4, 1 ) = (0), ( 1, 2 ) = (0), ( 2, 3 ) = (0), ( 1, 3 ) = (0), ( 2, 2 ) = (1), ( 4, 2 ) = (0), ( 3, 4 ) = (0), ( 1, 4 ) = (0), ( 3, 1 ) = (0), ( 4, 4 ) = (1), ( 3, 2 ) = (0), ( 1, 1 ) = (1), ( 2, 1 ) = (0), ( 4, 3 ) = (0), ( 3, 3 ) = (1), ( 2, 4 ) = (0)  ] ), array( 1 .. 4, 1 .. 1, [( 4, 1 ) = (psi[4]), ( 3, 1 ) = (psi[3]), ( 1, 1 ) = (psi[1]), ( 2, 1 ) = (psi[2])  ] ))

    (19)

    And to perform all these matricial operations within an algebraic expression,

    Library:-PerformMatrixOperations(Physics[`*`](`𝕀`, Psi))

    Matrix(%id = 18446744079185513758)

    (20)

    ``

     


     

    Download DiracAlgebraWithIdentityMatrix.mw

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

    Hey guys...

    One of my favorite new images rendered in Maple.....

    integrates mod congruency arguments like the following, into an 8 stage animated 3D cylindrical coordinates context::

     

    u*mod(x^(1/3)cos(x),23)

     

     

    As a powerlifter, I constantly find myself calculating my total between my competition lifts, bench, squat, and deadlift. Following that, I always end up calculating my wilks score. The wilks score is a score used to compare lifters between weight classes (https://en.wikipedia.org/wiki/Wilks_Coefficient ), as comparing someone who weighs 59kg in competition like me, to someone who weighs 120kg+, the other end of the spectrum; obviously the 120kg+ lifter is going to massively out-lift me.

    So I decided to program a wilks calculator for quick use, rather than needing to go search for one on the internet. For anyone curious about specific scoring, a score of 300+ is very strong for the average gym goer, and is about normal for a powerlifter. A score of 400+ makes you strong for a powerlifter, putting you in the top 250 powerlifters, while 500+ is the very top, as far as unequipped powerlifting goes, including the top 30. For anyone wondering, my best score at my best meet was 390, although given the progress I’ve made in the gym, should be above 400 by my next meet.

    Hope you all enjoy!

     Find it here: https://maple.cloud#doc=5687076260413440&key=301A440EFD2C4EDD8480D60B5E3147BF40CA460F842942449C939AB8D2E7D679

    The problem is to analyze the behavior of an involute of the cubical parabola near the singular points and at infinity.

    To do this efficiently, it is best to work with partially evaluated expressions.

    We want to investigate the singular points of an involute of the cubical parabola y = x^3.

    If the curve is given parametrically by "conjugate(r)(t)", an involute is given by

    "conjugate(E)(t)=conjugate(r)(t)-(conjugate(r)'(t))/(|conjugate(r)'(t)|)*(∫)[a]^(t)|conjugate(r)'(tau)| ⅆtau"

    What we want to do is to leave the integral unevaluated and compute only its derivatives. The simplest way is to define S := proc (a, t) options operator, arrow; Int(s(tau), tau = a .. t) end proc. Both diff and evalf will be handled automatically by Maple. We'll do it in a slightly more complicated way, leaving S(a, t) completely inert with the exception of the condition S(a, a)=0 and implementing a more efficient numerical evaluator.

    forget(diff); forget(series); r := proc (t) options operator, arrow; [t, t^3] end proc; s := unapply(sqrt((diff(r(t)[1], t))^2+(diff(r(t)[2], t))^2), t); S := proc (a, t) options operator, arrow; `if`(t = a, 0, 'S(a, t)') end proc; `diff/S` := proc (at, ft, t) options operator, arrow; s(ft)*(diff(ft, t))-s(at)*(diff(at, t)) end proc; E := unapply(r(t)-`~`[`*`](diff(r(t), t), S(a, t)/s(t)), [a, t])

    forget(evalf); `evalf/S` := (proc () local acur, Ssol; acur := undefined; Ssol := subs(dsolve({diff(f(t), t) = s(t), f(a) = 0}, f(t), numeric, parameters = [a], output = listprocedure), f(t)); proc (a, t) if [a, t]::(list(numeric)) then if a <> acur then Ssol(parameters = [a]); acur := a end if; Ssol(t) else 'S(a, t)' end if end proc end proc)()

    plot([[op(r(t)), t = -1 .. 3], [op(E(1, t)), t = -1.5 .. 4]], view = -1 .. 3, scaling = constrained)

     

    simplify(diff(E(a, t), t))

    [18*S(a, t)*t^3/(9*t^4+1)^(3/2), -6*t*S(a, t)/(9*t^4+1)^(3/2)]

    (1)

    The singular points are easily determined now: we need either S(a, t)=0, which implies t=a, or t=0. We assume a>0.

    Compute the Taylor series around t=a first.

    `~`[series](E(a, t+a), t = 0, 4); ser := collect(convert(%, polynom), t, normal)

    [-18*a^2*(6*a^4-1)*t^3/(9*a^4+1)^2+9*a^3*t^2/(9*a^4+1)+a, 2*(36*a^4-1)*t^3/(9*a^4+1)^2-3*a*t^2/(9*a^4+1)+a^3]

    (2)

    The center of the expansion E(a, a) is the point [a, a^3] on the cubical parabola.

    There are two quadratic terms, which give us the slope of the tangent line. The tangent coincides with the normal to the cubical parabola at this point.

    Rotate the axes to make the tangent vertical.

    rot := proc (e, a) options operator, arrow; convert(Student:-LinearAlgebra:-RotationMatrix(a).Vector(e), list) end proc

    rot(ser-E(a, a), (`@`(arctan, op))(`~`[coeff](ser, t, 2))); ser := `assuming`([collect(%, t, normal)], [a > 0])

    [-12*a^2*t^3/(9*a^4+1)^(3/2), -2*(162*a^8+9*a^4-1)*t^3/(9*a^4+1)^(5/2)+3*a*t^2/(9*a^4+1)^(1/2)]

    (3)

    In the new coordinates, the leading terms are t^3 and t^2, and the involute has the shape of a semicubical parabola. When t increases, the involute moves from the first to the second quadrant. In terms of x and y, `~`[x]*alpha*y^(3/2), and the coefficient alpha is found from the two leading terms. For t>0:

    `assuming`([simplify(coeff(ser[1], t, 3)/coeff(ser[2], t, 2)^(3/2))], [a > 0])

    -(4/3)*a^(1/2)*3^(1/2)/(9*a^4+1)^(3/4)

    (4)

    Now consider the second singular point t=0 in the original coordinate system.

    ser := convert(`~`[series](E(a, t), t = 0, 6), polynom)

    [-S(a, 0)+(9/2)*S(a, 0)*t^4+(18/5)*t^5, -3*S(a, 0)*t^2-2*t^3]

    (5)

    The center is the point [-S(a, 0), 0], which corresponds to the point [0, 0] on the cubical parabola. The leading terms t^4 and t^2 indicate that the tangent line is vertical and that this is a return point, with both branches lying in the second quadrant in the coordinate system shifted by [-S(a, 0), 0].

    In terms of x and y, `~`[x]*alpha*y^2, and to determine the next term, we need to analyze the two series expansions. t^4 gives a t^4 term and fixes alpha. t^2*(t^3) gives a t^5 term, but the coefficient doesn't match the coefficient at t^5 in x. So the next term has to be y^(5/2) in order to match t^5. For t>0, (b*t^3+a*t^2)^(5/2) = t^5*(b*t+a)^(5/2), which gives a^(5/2)*t^5 plus higher order terms. No other terms contribute, and we can equate the coefficients at t^4 and t^5 in x and y:

    collect(alpha*ser[2]^2+beta*coeff(ser[2], t, 2)^(5/2)*t^5, t); (`@`(simplify, solve))(`~`[coeff](%-ser[1], t, [4, 5]), {alpha, beta})

    4*alpha*t^6+(beta*(-3*S(a, 0))^(5/2)+12*alpha*S(a, 0))*t^5+9*alpha*S(a, 0)^2*t^4

     

    {alpha = (1/2)/S(a, 0), beta = -(4/45)*3^(1/2)/(-S(a, 0))^(5/2)}

    (6)

    We have obtained the coefficients in `~`[x]*alpha*y^2+beta*y^(5/2). alpha and beta are negative, while for t<0, beta has the opposite sign, and the branch corresponding to t<0 is the one which is closer to the vertical tangent.

    This could have been done by using Groebner:-Basis to eliminate t from the two series and then using algcurves:-puiseux. But it is still necessary to determine the truncation order.

    There is also a separate case a=0.

    ser := `~`[series](E(0, t), t = 0, 6)

    [series((18/5)*t^5+O(t^6),t,6), series(-2*t^3+O(t^7),t,7)]

    (7)

    Now the singularity is an inflection point, and the involute has a vertical tangent and moves from the second to the fourth quadrant. In the fourth quadrant, `~`[x]*alpha*(-y)^(5/3), and we find alpha in the same way as before:

    coeff(ser[1], t, 5)/(-coeff(ser[2], t, 3))^(5/3)

    (9/10)*2^(1/3)

    (8)

    Finally let's consider the asymptotic behavior of an involute at infinity. Now we do have to investigate S(a, t). One possible way is to apply the binomial formula to s(t) = 3*t^2*sqrt(1+1/(9*t^4))  and integrate term by term:

    `assuming`([int(3*tau^2*binomial(1/2, k)*(9*tau^4)^(-k), tau = a .. t)], [0 < a and a < t])

    3*binomial(1/2, k)*(9^(-k)*a^(-4*k+3)-9^(-k)*t^(-4*k+3))/(4*k-3)

    (9)

    This is the power series expansion of S(a, t) for large t, valid for 9*a^4 > 1. Substituting it into E(a, t), we obtain the required asymptotics.

    E(a, t)

    [-S(a, t)/(9*t^4+1)^(1/2)+t, -3*t^2*S(a, t)/(9*t^4+1)^(1/2)+t^3]

    (10)

    In the x component, we get -(1/3)*t+o(1)+t. In the y component, we get -t^3-C+o(1)+t^3. Thus the involute approaches a horizontal asymptote when t goes to +infinity, and the constant term -C gives the position of the asymptote:

    `C__+` := `assuming`([unapply(sum(3*binomial(1/2, k)*9^(-k)*a^(-4*k+3)/(4*k-3), k = 0 .. infinity), a)], [a > 1/sqrt(3)])

    proc (a) options operator, arrow; -a^3*hypergeom([-3/4, -1/2], [1/4], -(1/9)/a^4) end proc

    (11)

    That in fact extends to all positive a. For the asymptotics at -infinity, first we compute the integral over [a, -a], this time using the power series for small t. For the integral from -a to -t, S(-a, -t) = -S(a, t), and we need to subtract the value found above:

    `C__-` := `assuming`([unapply(sum(int(binomial(1/2, k)*(9*tau^4)^k, tau = a .. -a), k = 0 .. infinity)-`C__+`(a), a)], [0 < a and a < 1/sqrt(3)])

    proc (a) options operator, arrow; -2*a*hypergeom([-1/2, 1/4], [5/4], -9*a^4)+a^3*hypergeom([-3/4, -1/2], [1/4], -(1/9)/a^4) end proc

    (12)

    For negative a, `C__+` and `C__-` are swapped.

    This could have been done by converting S(a, t) into a difference of two integrals of hypergeometric type, for which Maple is able to find the asymptotics automatically.

    inv := proc (a, T) plots:-display(plot([op(r(t)), t = -1 .. 2], color = black), plot([op(E(a, t)), t = -1.5 .. T], color = red), plottools:-line(r(T), E(a, T), linestyle = dot), map(proc (y) options operator, arrow; plottools:-line([-1, y], [3, y], linestyle = spacedot) end proc, [-`C__+`(a), -`C__-`(a)]), scaling = constrained, view = [-1 .. 3, -1 .. 3]) end proc

    Explore(inv(a, t), a = -.3 .. 1.3, t = -1. .. 4., initialvalues = [a = 1., t = 3.])

    ``

    Download cubic.mw

    I submit a bug through MaplePrimes because I can't do it as usually (Hope some people understand me.). Let us consider

    with(LinearAlgebra):
    M := Matrix(5, 5,  (i, j) -> (10*i+j)*sin((1/180)*Pi*(10*i+j))):
    MatrixInverse(M);
     #One sees a long and wrong output instead of the warning "Matrix M is singular"
    

    Indeed,

    Digits := 500; evalf(Determinant(M), 495);
                                   
                               1.3 10 ^(-488)   
    

    Bug_in_MatrixInverse.mw

    I submit a bug through MaplePrimes because I can't do it as usually (Hope some people understand me.). Let us consider

    restart; pdsolve([diff(u(t, x), t, t) = diff(u(t, x), x, x), u(t, 0) = 0, u(t, Pi) = 0]);
    pdsolve([diff(u(t, x), t, t) = diff(u(t, x), x, x), u(t, 0) = 0, u(t, Pi) = 0], generalsolution);
    u(t, x) = Sum(sin(n*x)*(_C5(n)*cos(n*t)+_C1(n)*sin(n*t)), n = 1 .. infinity)
    u(t, x) = Sum(sin(n*x)*(_C5(n)*cos(n*t)+_C1(n)*sin(n*t)), n = 1 .. infinity)
    

    The question arises: what do these outputs mean? I don't see any explanation in ?pdsolve and ?examples,pdsolve_boundaryconditions. What are _C1(n) and _C5(n)? Under which conditions does the above series converge?

    Moreover,

    pdetest(%, [diff(u(t, x), t, t) = diff(u(t, x), x, x), u(t, 0) = 0, u(t, Pi) = 0]);
                               [0, 0, 0]
    

    I think the above is simply a fake: it is possible to differentiate  a series only under certain conditions.

    Bug_in_pdsolve.mw

    Please, don't convert my post to a question. This is not correct and fair. Hope some people understand me.

    Inspired by this question. There is a simple iterated procedure that can generate the Puiseux expansion.

    Note the use of convert(..., rational, exact) to preserve the digits of the floating-point numbers.

    (proc () global COcrit, P; Digits := 30; COcrit := 1/convert(2*0.115e-12, rational, exact); P := (`@`(`@`(rcurry(unapply, [NO2, O3]), numer), rcurry(convert, rational, exact)))(1.8027350279544425625*O3^8/10^49+(2.982942852010948125*CO/10^49+2.27611508925183215625*NO2/10^47+3.754849807690565625/10^37)*O3^7+(1.2339511797887015625*CO^2/10^49+2.4836920553034140625*CO/10^37+(-1)*5.9759150186390484375*NO2/10^35+6.3862221928648528125*NO2^2/10^46+1.88311680046014890625*CO*NO2/10^47+(-1)*9.69704144499346875/10^24)*O3^6+((-1)*1.71193039685859375*CO^2/10^37+(-1)*5.7098636544065625*CO/10^24+(-1)*1.277325786277575*NO2/10^21+(-1)*1.0801570017944671875*NO2^2/10^32+(-1)*2.9081778565815421875*CO*NO2/10^34+(-1)*3.66453227489203125/10^11)*O3^5+(1.9152220505625*CO^2/10^25+8.1035862984375*CO/10^13+1.40846609345625*NO2/10^10+(-1)*5.19285353257125*NO2^2/10^21+(-1)*1.55036925507375*CO*NO2/10^21+(-1)*1.844759695120875*CO^2*NO2/10^34+(-1)*1.876842472427325*CO*NO2^2/10^32-7.201634275625)*O3^4+(1.793091274970625*NO2^2/10^7+8618.14231275*NO2+(-1)*2.298266460675*CO^2*NO2/10^22+(-1)*9.5902239009375*CO*NO2/10^10+(-1)*1.685705248305*CO*NO2^2/10^20+9.2666503797075*NO2^3/10^19)*O3^3+((-1)*2.5638555726*10^6*NO2^2+6.894799382025*NO2^2*CO^2/10^20+2.7563954788125*CO*NO2^2/10^7+3.5073544682475*NO2^3*CO/10^18+(-1)*0.604340578881e-4*NO2^3+(-1)*3.5683519372605*NO2^4/10^16)*O3^2+((-1)*8.75499651*10^6*NO2^3+0.482686765875e-5*NO2^3*CO+(-1)*0.98216263665e-4*NO2^4)*O3+5.4066609*10^7*NO2^4) end proc)()

    We're interested in the asymptotics of RootOf(P(NO2, _Z)) for small NO2.

    plot3d(RootOf(P(NO2, _Z)), CO = .5*COcrit .. 2*COcrit, NO2 = 10^(-3) .. 10^(-6))

    Start with P(NO2, Z) and look for the expansion for Z when NO2 is small.

    Expand P(NO2, Z) into monomials and convert each monomial a*Typesetting:-mi("NO2", italic = "true", mathvariant = "italic")^Typesetting:-mi("p", italic = "true", mathvariant = "italic")*Typesetting:-mi("Z", italic = "true", mathvariant = "italic")^Typesetting:-mi("q", italic = "true", mathvariant = "italic") into the point [p, q].

    ([op])(collect(P(NO2, Z), [NO2, Z], distributed, normal)); map(proc (e) options operator, arrow; `~`[degree](e, [NO2, Z]) end proc, %)

    [[4, 2], [4, 1], [4, 0], [3, 3], [3, 2], [3, 1], [2, 6], [2, 5], [2, 4], [2, 3], [2, 2], [1, 7], [1, 6], [1, 5], [1, 4], [1, 3], [0, 8], [0, 7], [0, 6], [0, 5], [0, 4]]

    (1)

    The Newton polygon is the convex hull of the points.

    newtonPoly := proc (P, xy) local terms, pts; terms := ([op])(collect(P, xy, distributed, normal)); pts := map(proc (e) options operator, arrow; `~`[degree](e, xy) end proc, terms); plots:-display(plottools:-polygon(simplex:-convexhull(pts)), plottools:-point(pts), symbol = solidcircle, symbolsize = 15, style = line, thickness = 3, color = [magenta, black], axis = [gridlines = spacing(1)]) end proc

    newtonPoly(P(NO2, Z), [NO2, Z])

     

    The side closest to the origin will correspond to the asymptotics for small NO2.

    Sum the corresponding monomials and solve for Z.

    sideAsympt := proc (P, xy, pts) local lead; lead := add(coeff(coeff(P, xy[1], p[1]), xy[2], p[2])*xy[1]^p[1]*xy[2]^p[2], `in`(p, pts)); ([solve])(lead, xy[2]) end proc

    sideAsympt(P(NO2, Z), [NO2, Z], [[0, 4], [1, 3], [2, 2], [3, 1], [4, 0]])

    [600*NO2, 600*NO2, -57000000000000000*NO2/(4071*CO-17795000000000000), -228000000000000000*NO2/(4071*CO+35020000000000000)]

    (2)

    We have obtained the first term in the expansion. If CO>COcrit, the smallest positive root is the one asymptotic to 600*NO2. That will be the value of RootOf(P(NO2, _Z)).

    The expansion can be continued in the same manner.

    newtonPoly(P(NO2, NO2*(600+Z)), [NO2, Z])

     

    In this case we don't have to choose between sides with different slopes. (Compare to "P(NO2,600*NO2+Z).)"

    sideAsympt(P(NO2, NO2*(600+Z)), [NO2, Z], [[4, 2], [5, 0]]); `assuming`([simplify(`assuming`([simplify(%)], [NO2 > 0]))], [CO > COcrit])

    [(531/430000)*354^(1/2)*NO2^(1/2)*(200000000000000+23*CO)^(1/2)/(23*CO-100000000000000)^(1/2), -(531/430000)*354^(1/2)*NO2^(1/2)*(200000000000000+23*CO)^(1/2)/(23*CO-100000000000000)^(1/2)]

    (3)

    We get `~`[Z]*NO2^(1/2), so the next order in the expansion is NO2^(3/2).

    Again, if we want the solution that corresponds to RootOf(P(NO2, _Z)), we should take the negative term, as the principal root will be the smaller one.

    Find one more term. To avoid fractional powers, take NO2r = NO2^(1/2).

    approx := 600*NO2r^2+NO2r^3*(-(531/430000)*sqrt(354)*sqrt(200000000000000+23*CO)/sqrt(23*CO-100000000000000)+Z)

    newtonPoly(P(NO2r^2, approx), [NO2r, Z])

     

    sideAsympt(P(NO2r^2, approx), [NO2r, Z], [[10, 1], [11, 0]])

    [(4779/36980000000000)*(4341503*CO^2+29548445000000000000*CO-209981000000000000000000000000000)*NO2r/(529*CO^2-4600000000000000*CO+10000000000000000000000000000)]

    (4)

    We get `~`[Z]*NO2r, so we have obtained the coefficient at NO2r^3*NO2r = NO2^2.

    All of this can be done using the algcurves:-puiseux command. The only difficulty is the unwieldy expressions that algcurves:-puiseux will generate for P.

    Let's also find the NO2^2 term by the method of dominant balance. Suppose that we don't know p yet and are looking for the term alpha*NO2^p = alpha*NO2p.

    approx := 600*NO2r^2-(531/430000)*sqrt(354)*sqrt(200000000000000+23*CO)*NO2r^3/sqrt(23*CO-100000000000000)+alpha*NO2p

    terms := (`@`([op], collect))(P(NO2r^2, approx), [NO2r, NO2p], distributed, normal)

    Find all exponents of NO2, coming from NO2r and from NO2p.

    pows := map(proc (t) options operator, arrow; (1/2)*degree(t, NO2r)+p*degree(t, NO2p) end proc, terms)

    plot(pows, p = 1 .. 2, view = 4 .. 6, annotation = pows)

     

    At p=2, the two dominant terms NO2^(11/2) and NO2^(7/2+p) can balance each other.

    (`@`(normal, coeff))(subs(NO2p = NO2r^4, P(NO2r^2, approx)), NO2r, 11); ([solve])(%, alpha)

    [(4779/36980000000000)*(4341503*CO^2+29548445000000000000*CO-209981000000000000000000000000000)/(529*CO^2-4600000000000000*CO+10000000000000000000000000000)]

    (5)

    ``

    Download np.mw

    First 46 47 48 49 50 51 52 Last Page 48 of 302