ecterrab

The equations of motion in curvilinear c...

Posted: ecterrab 14540 Product: Maple

May 21 2020

5 0

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

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

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

Vector notation

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

>

(1)

>

(2)

>

(3)

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

>

(4)

where represents the total force acting on the particle. This vectorial equation represents three second order differential equations which, for given initial conditions, can be integrated to arrive at a closed form expression for as a function of .

Tensor notation

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

>

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

(5)

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

>

(6)

>

(7)

>

(8)

Setting a tensor to represent the three Cartesian components of the force

>

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

(9)

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

>

(10)

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

>

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

(11)

Things are straightforward in Cartesian coordinates because the components of the line element are exactly the components of the tensor

>

(12)

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

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

>

(13)

Curvilinear coordinates

Vector notation

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

>

(14)

This change keeps the z axis unchanged, so the corresponding unit vector remains unchanged.

Changing the basis and coordinates used to represent the position vector , it becomes

>

(15)

where since in (1) the coordinates () depend on t, in (15), not just and but also the unit vector depends on t. The velocity is computed as usual, differentiating

>

(16)

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

For the acceleration,

>

(17)

where the term involving comes from differentiating in (16) taking into account the dependency of on the coordinate This result can also be obtained by directly changing variables in the acceleration , in equation (3)

>

(18)

Newton's 2^nd law becomes

>

(19)

In the absence of external forces, equating to 0 the vector components (coefficients of the unit vectors) of the acceleration we get the system of differential equations in cylindrical coordinates whose solution is the trajectory of the particle expressed in the (

>

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

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

(20)

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

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

(21)

In this formulation (21) with , although , no acceleration in the direction, is naturally expected, the same cannot be said about the other two equations for and . Those two equations are discussed below under Coriolis and Centripetal forces. The key observation at this point, however, is that the right-hand sides of both unexpected equations involve , rotation around the z axis.

Tensor notation

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

>

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

(22)

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

>

(23)

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

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

1.	Make be a function of directly in the metric

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

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

For step 1, we have

>

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

(24)

Set this metric where

>

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

(25)

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

>

(26)

Step 3: compute the components of this tensorial equation

>

(27)

These are the same equations (23).

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

>

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

(28)

Newton's 2^nd law (19)

>

(29)

now using tensorial notation, becomes

>

(30)

>

Array(%id = 18446744078329063774)

(31)

where we recall (see related Mapleprimes post) that to obtain the vector components entering from these tensor components we need to multiply the latter by the scale factors (), the component of in the direction of is given by .

Coriolis force and centripetal force

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

A distinction needs to be made here, according to whether the unit vector depends or not on the time , the former being the general case. When is a constant, the value of the coordinate - the angle between and the x axis - does not change, there is no rotation around the z axis. On the other hand, when , the orientation of and so the coordinate changes with time, so either the force acting on the particle has a component in the direction that produces rotation around the z axis, or the system of references - itself - is rotating around the z axis.

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

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

To see that the equations

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

that appeared in (27) when in the inertial system of references , are related to the Coriolis and centripetal forces in the non-inertial referencial, following [1] introduce a vector representing the rotation of that referencial around the z axis (when, in the inertial system of references, the non-inertial system rotates clockwise, in the non-inertial system increases in value in the anti-clockwise direction)

>

(32)

According to [1], (39.7), the acceleration in the inertial system is expressed in terms of the quantities in the non-inertial rotating system by the sum of the following three vectorial terms.

First, the components of the acceleration measured in the non-inertial system are given by the second derivatives of the coordinates () multiplied by the scale factors, which in this case are given by () (see this post in Mapleprimes)

>

(33)

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

>

(34)

Third the centripetal force divided by the mass, defined by

>

(35)

Adding these three terms,

>

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

(36)

So that

>

(37)

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

>

(38)

>

(39)

From (37), in the absence of external forces and so the acceleration measured in the rotating system is given by the sum of the Coriolis and centripetal accelerations

>

(40)

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

Expressing this equation (40) in terms of () we get

>

(41)

resulting in the three equations

>

(42)

>

(43)

>

(44)

which are the equations returned by Geodesics in (23)

>

(45)

>

References

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

Download The_equations_of_motion_in_curvilinear_coordinates.mw

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

Maple for Beginners...

Posted: ecterrab 14540 Product: Maple

April 15 2020

11 0

Maple_for_Beginners.mw

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

Vectors in Spherical Coordinates using T...

Posted: ecterrab 14540 Product: Maple

April 07 2020

6 3

Vectors in Spherical Coordinates using Tensor Notation

Edgardo S. Cheb-Terrab¹ and Pascal Szriftgiser²

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

(1) Maplesoft

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

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

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

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

>

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

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

(1)

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

In vector calculus, at the root of everything there is the line element , which in Cartesian coordinates has the simple form

>

(1.1)

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

>

(1.2)

>

(1.3)

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

>

So in matrix notation,

>

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

(1.4)

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

In two steps:

>

(1.5)

The line element

>

(1.6)

This result is important: it gives us the so-called scale factors, the key that connect 3D vectors with the related covariant and contravariant tensors in curvilinear coordinates. The scale factors are computed from (1.6) by taking the scalar product with each of the unit vectors , then taking the coefficients of the differentials (just substitute them by the number 1)

>

(1.7)

The scale factors are relevant because the components of the 3D vector and the corresponding tensor are not the same in curvilinear coordinates. For instance, representing the differential of the coordinates as the tensor , we see that corresponding vector, the line element in spherical coordinates , is not constructed by directly equating its components to the components of , so

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

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

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

where on the right-hand side we see the contravariant components and the scale-factors . Because the system is orthogonal, each vector component satisfies

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

it follows that,

A__j = A__j/h__j

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

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

II. Transformation of contravariant and covariant tensors

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

>

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

(2.1)

Transformation rule for a contravariant tensor

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

where on the left-hand side we get, isolated, the three components of A in Cartesian coordinates, and on the right-hand side we transform the spherical components , from spherical (4^th argument) to Cartesian (3^rd argument), which according to the 5^th bullet of TransformCoordinates will result in a transformation expressed in terms of the old coordinates (here the spherical ). Expand things to make the comparison with (1.4) possible by eye

>

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

(2.2)

We see that the transformation rule for a contravariant vector is, indeed, as the transformation (1.4) for the differential of the coordinates.

Transformation rule for a covariant tensor

For the transformation rule for the components of a covariant tensor , we know, by definition, that it is $`#mrow(msub(mi("A"),mi("μ",fontstyle = "normal")),mo("⁡"),mfenced(mi("y")),mo("="),mfrac(mrow(mo("&PartialD;"),msup(mi("x"),mi("ν",fontstyle = "normal"))),mrow(mo("&PartialD;"),msup(mi("y"),mi("μ",fontstyle = "normal"))),linethickness = "1"),mo("⁢"),mo("⁢"),msub(mi("A"),mi("ν",fontstyle = "normal")),mfenced(mi("x")))`$ , so the same transformation rule for the gradient , where and so on. We can experiment this by directly changing variables in the differential operators , for example

>

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

(2.3)

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

>

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

(2.4)

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

>

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

(2.5)

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

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

>

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

(2.6)

III. Deriving the transformation rule for the Gradient using TransformCoordinates

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

>

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

>

(3.1)

>

(3.2)

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

Start changing the vector basis in the gradient (3.1)

>

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

(3.3)

By eye, we see that in this result the coefficients of are the three lines in the right-hand side of (2.6) after replacing the covariant components by the derivatives of f with respect to the j^th coordinate, here displayed using indexed notation due to using CompactDisplay

>

(3.4)

>

(3.5)

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

>

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

(3.6)

>

(3.7)

Simplifying, we arrive at (3.2)

>

(3.8)

>

(3.9)

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

•	The Divergence

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

>

(4.1)

>

(4.2)

>

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

(4.3)

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

>

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

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

(4.4)

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

>

(4.5)

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

>

(4.6)

>

(4.7)

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

>

(4.8)

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

>

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

(4.9)

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

>

(4.10)

>

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

(4.11)

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

>

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

(4.12)

•

The Curl

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

>

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

(4.13)

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

But in a curvilinear system is not a tensor, we need to use the non-Galilean form , where is the determinant of the metric. Moreover, since the expression has one free covariant index (the first one), to compare with the vectorial formula (4.12) this index also needs to be rewritten as a vector component as discussed at the end of Sec. I, using

A__j = A__j/h__j

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

>

(4.14)

>

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

(4.15)

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

>

(4.16)

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

>

(4.17)

>

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

(4.18)

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

>

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

(4.19)

•	The Gradient

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

>

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

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

(4.20)

where on the left-hand side we have the vectorial Nabla differential operator and on the right-hand side, since is a scalar, the covariant derivative becomes the standard derivative .

>

(4.21)

The above is the expected result (3.2)

>

(4.22)

•	The Laplacian

Likewise we can compute the Laplacian directly as

>

(4.23)

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

>

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

(4.24)

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

>

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

(4.25)

On the right-hand side we see the dAlembertian , in curvilinear coordinates; rewrite it using standard diff derivatives and expand both sides of the equation for comparison

>

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

(4.26)

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

>

(4.27)

Download Vectors_and_Spherical_coordinates_in_tensor_notation.mw

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

Feynman Diagrams - the scattering matrix...

Posted: ecterrab 14540 Product: Maple

January 25 2020

7 0

Feynman Diagrams
The scattering matrix in coordinates and momentum representation

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

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

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

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

Introduction

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

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

This exponential can be expanded as

where

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

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

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

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

b.

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

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

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

Examples

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

>

(1.1)

Let be the interaction Lagrangian

>

(1.2)

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

>

(1.3)

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

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

>

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

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

(1.4)

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

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

>

(1.5)

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

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

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

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

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

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

>

(1.6)

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

>

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

(1.7)

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

>

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

(1.8)

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

>

(1.9)

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

>

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

(1.10)

A corresponding S-matrix element in momentum representation:

>

(1.11)

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

>

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

(1.12)

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

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

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

(1.13)

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

>

(1.14)

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

>

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

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

(1.15)

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

>

(1.16)

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

>

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

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

(1.17)

The corresponding S-matrix elements in momentum representation

>

(1.18)

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

>

(1.19)

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

>

(1.20)

The self-energy of the photon:

>

(1.21)

where is the corresponding polarization vector.

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

>

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

(1.22)

>

(1.23)

>

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

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

(1.24)

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

>

(1.25)

(1.26)

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

>

(1.27)

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

>

(1.28)

>

(1.29)

>

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

(1.30)

>

(1.31)

>

(1.32)

>

(1.33)

>

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

(1.34)

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

>

References

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

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

FeynmanDiagrams_and_the_Scattering_Matrix.PDF

FeynmanDiagrams_and_the_Scattering_Matrix.mw

FeynmanDiagrams_-_help_page.mw

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

Splitting PDE parameterized symmetries a...

Posted: ecterrab 14540 Product: Maple

October 18 2019

5 6

Splitting PDE parameterized symmetries

and Parameter-continuous symmetry transformations

The determination of symmetries for partial differential equation systems (PDE) is relevant in several contexts, the most obvious of which is of course the determination of the PDE solutions. For instance, generally speaking, the knowledge of a N-dimensional Lie symmetry group can be used to reduce the number of independent variables of PDE by N. So if PDE depends only on N independent variables, that amounts to completely solving it. If only N-1 symmetries are known or can be successfully used then PDE becomes and ODE; etc., all advantageous situations. In Maple, a complete set of symmetry commands, to perform each step of the symmetry approach or several of them in one go, is part of the PDEtools package.

Besides the dependent and independent variables, PDE frequently depends on some constant parameters, and besides the PDE symmetries for arbitrary values of those parameters, for some particular values of them, PDE transforms into a completely different problem, admitting different symmetries. The question then is: how can you determine those particular values of the parameters and the corresponding different symmetries? That was the underlying subject of a recent question in Mapleprimes. The answer to those questions is relatively simple and yet not entirely obvious for most of us, motivating this post, organized briefly around one example.

To reproduce the input/output below you need Maple 2019 and to have installed the Physics Updates v.449 or higher.

Consider the family of Korteweg-de Vries equation for involving three constant parameters . For convenience (simpler input and more readable output) use the diff_table and declare commands

>

(1)

>

(2)

This pde admits a 4-dimensional symmetry group, whose infinitesimals - for arbitrary values of the parameters - are given by

>

(3)

Looking at pde (1) as a nonlinear problem in u, a, b and q, it splits into four cases for some particular values of the parameter:

>	$pde__cases := casesplit(bu(x, t)(diff(u(x, t), x))+a(diff(u(x, t), x))+q(diff(diff(diff(u(x, t), x), x), x))+diff(u(x, t), t) = 0, parameters = {a, b, q}, caseplot)$

`casesplit/ans`([diff(diff(diff(u(x, t), x), x), x) = -(b*u(x, t)*(diff(u(x, t), x))+a*(diff(u(x, t), x))+diff(u(x, t), t))/q], [q <> 0]), `casesplit/ans`([diff(u(x, t), x) = -(diff(u(x, t), t))/(b*u(x, t)+a), q = 0], [b*u(x, t)+a <> 0]), `casesplit/ans`([u(x, t) = -a/b, q = 0], [b <> 0]), `casesplit/ans`([diff(u(x, t), t) = 0, a = 0, b = 0, q = 0], [])

(4)

The legend above indicates the pivots and the tree of cases, depending on whether each pivot is equal or different from 0. At the end there is the algebraic sequence of cases. The first case is the general case, for which the symmetry infinitesimals were computed as above, but clearly the other three cases admit more general symmetries. Consider for instance the second case, pass the to ignore the parameterizing equation , and you get

>

[_xi[x](x, t, u) = _F3(x, t, u), _xi[t](x, t, u) = Intat(((b*u+a)*(D[1](_F3))(_a, ((b*u+a)*t-x+_a)/(b*u+a), u)-_F1(u, ((b*u+a)*t-x)/(b*u+a))*b+(D[2](_F3))(_a, ((b*u+a)*t-x+_a)/(b*u+a), u))/(b*u+a)^2, _a = x)+_F2(u, ((b*u+a)*t-x)/(b*u+a)), _eta[u](x, t, u) = _F1(u, ((b*u+a)*t-x)/(b*u+a))]

(5)

These infinitesimals are indeed much more general than , in fact so general that (5) is almost unreadable ... Specialize the three arbitrary functions into something "easy" just to be able follow - e.g. take _F1 to be just the + operator, _F2 the * operator and

>

[_xi[x](x, t, u) = 1, _xi[t](x, t, u) = Intat(-(u+((b*u+a)*t-x)/(b*u+a))*b/(b*u+a)^2, _a = x)+u*((b*u+a)*t-x)/(b*u+a), _eta[u](x, t, u) = u+((b*u+a)*t-x)/(b*u+a)]

(6)

>

[_xi[x](x, t, u) = 1, _xi[t](x, t, u) = (b^3*t*u^4+((3*a*t-x)*u^3-u^2*x-t*x*u)*b^2+((3*a^2*t-2*a*x)*u^2-a*u*x-a*t*x+x^2)*b+a^2*u*(a*t-x))/(b*u+a)^3, _eta[u](x, t, u) = (b*u^2+(b*t+a)*u+a*t-x)/(b*u+a)]

(7)

This symmetry is of course completely different than computed for the general case.

The symmetry (7) can be verified against or directly against pde after substituting .

(8)

>

(9)

>

(10)

Summarizing: "to split PDE symmetries into cases according to the values of the PDE parameters, split the PDE into cases with respect to these parameters (command PDEtools:-casesplit ) then compute the symmetries for each case"

Parameter continuous symmetry transformations

A different, however closely related question, is whether pde admits "symmetries with respect to the parameters a, b and q", so whether exists continuous transformations of the parameters a, b and q that leave pde invariant in form.

Beforehand, note that since the parameters are constants with regards to the dependent and independent variables (here ), such continuous symmetry transformations cannot be used directly to compute a solution for pde. They can, however, be used to reduce the number of parameters. And in some contexts, that is exactly what we need, for example to entirely remove the splitting into cases due to their presence, or to proceed applying a solving method that is valid only when there are no parameters (frequently the case when computing exact solutions to "PDE & Boundary Conditions").

To compute such "continuous symmetry transformations of the parameters" that leave pde invariant one can always think of these parameters as "additional independent variables of pde". In terms of formulation, that amounts to replacing the dependency in the dependent variable, i.e. replace by

>

(11)

Compute now the infinitesimals: note there are now three additional ones, related to continuous transformations of and - for readability, avoid displaying the redundant functionality on the left-hand sides of these infinitesimals

>

[_xi[x] = (1/3)*(_F4(a, b, q)*q+_F3(a, b, q))*x/q+_F6(a, b, q)*t+_F7(a, b, q), _xi[t] = _F4(a, b, q)*t+_F5(a, b, q), _xi[a] = _F1(a, b, q), _xi[b] = _F2(a, b, q), _xi[q] = _F3(a, b, q), _eta[u] = (1/3)*((b*u+a)*_F3(a, b, q)-2*((b*u+a)*_F4(a, b, q)+(3/2)*u*_F2(a, b, q)+(3/2)*_F1(a, b, q)-(3/2)*_F6(a, b, q))*q)/(b*q)]

(12)

This result is more general than what is convenient for algebraic manipulations, so specialize the seven arbitrary functions of and keep only the first symmetry that result from this specialization: that suffices to illustrate the removal of any of the three parameters a, b, or q

>

(13)

To remove the parameters, as it is standard in the symmetry approach, compute a transformation to canonical coordinates, with respect to the parameter a. That means a transformation that changes the list of infinitesimals, or likewise its infinitesimal generator representation,

>

proc (f) options operator, arrow; diff(f, a)-(diff(f, u))/b end proc

(14)

into or its equivalent generator representation

That same transformation, when applied to , entirely removes the parameter a.

The transformation is computed using CanonicalCoordinates and the last argument indicates the "independent variable" (in our case a parameter) that the transformation should remove. We choose to remove a

>

(15)

>

(16)

Invert this transformation in order to apply it

>

(17)

The next step is not necessary, but just to understand how all this works, verify its action over the infinitesimal generator proc (f) options operator, arrow; diff(f, a)-(diff(f, u))/b end proc

>

ChangeSymmetry({a = alpha, b = beta, q = chi, t = tau, x = xi, u(x, t, a, b, q) = (upsilon(xi, tau, alpha, beta, chi)*beta-alpha)/beta}, proc (f) options operator, arrow; diff(f, a)-(diff(f, u))/b end proc, [upsilon(xi, tau, alpha, beta, chi), xi, tau, alpha, beta, chi])

(18)

Now that we see the transformation (17) is the one we want, just use it to change variables in

>	$PDEtools:-dchange({a = alpha, b = beta, q = chi, t = tau, x = xi, u(x, t, a, b, q) = (upsilon(xi, tau, alpha, beta, chi)*beta-alpha)/beta}, pde__xtabq, [upsilon(xi, tau, alpha, beta, chi), xi, tau, alpha, beta, chi], simplify)$

(19)

As expected, this result depends only on two parameters, , and , and the one equivalent to a (that is , see the transformation used (17)), is not present anymore.

To remove b or q we use the same steps, (15), (17) and (19), just changing the parameter to be removed, indicated as the last argument in the call to CanonicalCoordinates . For example, to eliminate b (represented in the new variables by ), input

>

(20)

>

(21)

>	$PDEtools:-dchange({a = beta, b = alpha, q = chi, t = tau, x = xi, u(x, t, a, b, q) = (upsilon(xi, tau, alpha, beta, chi)*alpha-beta)/alpha}, pde__xtabq, [upsilon(xi, tau, alpha, beta, chi), xi, tau, alpha, beta, chi], simplify)$

(22)

and as expected this result does not contain To remove a second parameter, the whole cycle is repeated starting with computing infinitesimals, for instance for (22). Finally, the case of function parameters is treated analogously, by considering the function parameters as additional dependent variables instead of independent ones.

Download How_to_split_symmetries_into_cases_(II).mw

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

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

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

Ask a Question

Create a Post

14540 Reputation

24 Badges

MaplePrimes Activity

The equations of motion in curvilinear c...

Maple for Beginners...

Vectors in Spherical Coordinates using T...

Feynman Diagrams - the scattering matrix...

Splitting PDE parameterized symmetries a...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again