ecterrab

13431 Reputation

24 Badges

19 years, 356 days

MaplePrimes Activity


These are Posts that have been published by ecterrab


In a recent question in Mapleprimes, a spacetime (metric) solution to Einstein's equations, from chapter 27 of the book of Exact Solutions to Einstein's equations [1] was discussed. One of the issues was about computing a tetrad for that solution [27, 37, 1] such that the corresponding Weyl scalars are in canonical form. This post illustrates how to do that, with precisely that spacetime metric solution, in two different ways: 1) automatically, all in one go, and 2) step-by-step. The step-by-step computation is useful to verify results and also to compute different forms of the tetrads or Weyl scalars. The computation below is performed using the latest version of the Maplesoft Physics Updates.

 

with(Physics)

Physics:-Version()

`The "Physics Updates" version in the MapleCloud is 851 and is the same as the version installed in this computer, created 2020, October 19, 13:47 hours Pacific Time.`

(1)

The starting point is this image of page 421 of the book of Exact Solutions to Einstein's equations, formulas (27.37)

 

Load the solution [27, 37, 1] from Maple's database of solutions to Einstein's equations

g_[[27, 37, 1]]

_______________________________________________________

 

`Systems of spacetime coordinates are:`*{X = (z, zb, r, u)}

 

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

 

`The `*`Robinson and Trautman (1962)`*` metric in coordinates `*[z, zb, r, u]

 

`Parameters: `*[P(z, zb, u), H(X)]

 

"`Comments: ` admits geodesic, shearfree, twistfree null congruence, rho=-1/r=rho_b"

 

`Resetting the signature of spacetime from `*`- - - +`*` to `*`+ + + -`*` in order to match the signature in the database of metrics`

 

_______________________________________________________

 

`Setting `*lowercaselatin_is*` letters to represent `*space*` indices`

 

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

(2)

"CompactDisplay(?)"

H(X)*`will now be displayed as`*H

 

P(z, zb, u)*`will now be displayed as`*P

(3)

The assumptions on the metric's parameters are

Assume(P(z, zb, u) > 0, (H(X))::real, r >= 0)

 

The line element is as shown in the second line of the image above

g_[lineelement]

2*r^2*Physics:-d_(z)*Physics:-d_(zb)/P(z, zb, u)^2-2*Physics:-d_(r)*Physics:-d_(u)-2*H(X)*Physics:-d_(u)^2

(4)

Load Tetrads

with(Tetrads)

_______________________________________________________

 

`Setting `*lowercaselatin_ah*` letters to represent `*tetrad*` indices`

 

((`Defined as tetrad tensors `*`see <a href='http://www.maplesoft.com/support/help/search.aspx?term=Physics,tetrads`*`,' target='_new'>?Physics,tetrads`*`,</a> `*`&efr;`[a, mu]*`, `)*eta[a, b]*`, `*gamma[a, b, c]*`, `)*lambda[a, b, c]

 

((`Defined as spacetime tensors representing the NP null vectors of the tetrad formalism `*`see <a href='http://www.maplesoft.com/support/help/search.aspx?term=Physics,tetrads`*`,' target='_new'>?Physics,tetrads`*`,</a> `*l[mu]*`, `)*n[mu]*`, `*m[mu]*`, `)*conjugate(m[mu])

 

_______________________________________________________

(5)

The Petrov type of this spacetime solution is

PetrovType()

"II"

(6)

The null tetrad computed by the Maple system using a general algorithms is

Setup(tetrad = null)

e_[]

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

(7)

 

According to the help page TransformTetrad , the canonical form of the Weyl scalars for each different Petrov type is

 

So for type II, when the tetrad is in canonical form, we expect only `&Psi;__2` and `&Psi;__3` different from 0. For the tetrad computed automatically, however, the scalars are

Weyl[scalars]

psi__0 = -P(z, zb, u)*(2*(diff(P(z, zb, u), z))*(diff(H(X), z))+P(z, zb, u)*(diff(diff(H(X), z), z)))/(r^2*(H(X)^2+1)^(1/2)), psi__1 = ((1/2)*I)*(-(diff(diff(H(X), r), z))*P(z, zb, u)^2*r+2*P(z, zb, u)^2*(diff(H(X), z))-(diff(P(z, zb, u), u))*(diff(P(z, zb, u), z))*r+(diff(diff(P(z, zb, u), u), z))*r*P(z, zb, u))/(P(z, zb, u)*r^2*(H(X)^2+1)^(1/4)), psi__2 = (1/6)*((diff(diff(H(X), r), r))*r^2+2*(diff(P(z, zb, u), z))*(diff(P(z, zb, u), zb))-2*(diff(H(X), r))*r-2*(diff(diff(P(z, zb, u), z), zb))*P(z, zb, u)+2*H(X))/r^2, psi__3 = 0, psi__4 = 0

(8)

The question is, how to bring the tetrad `&efr;`[a, mu] (equation (7)) into canonical form. The plan for that is outlined in Chapter 7, by Chandrasekhar, page 388, of the book "General Relativity, an Einstein centenary survey", edited by S.W. Hawking and W.Israel. In brief, for Petrov type II, use a transformation ofClass[2] to make Psi[0] = `&Psi;__1` and `&Psi;__1` = 0, then a transformation of Class[1] making Psi[4] = 0, finally use a transformation of Class[3] making Psi[3] = 1. For an explanation of these transformations see the help page for TransformTetrad . This plan, however, is applicable if and only if the starting tetrad results in `&psi;__4` <> 0, which we see in (8) it is not the case, so we need, in addition, before applying this plan, to perform a transformation of Class[1] making `&psi;__4` <> 0.

 

In what follows, the transformations mentioned are first performed automatically, in one go, letting the computer deduce each intermediate transformation, by passing to TransformTetrad the optional argument canonicalform. Then, the same result is obtained by transforming the starting tetrad  one step at at time, arriving at the same Weyl scalars. That illustrates well both how to get the result exploiting advanced functionality but also how to verify the result performing each step, and also how to get any desired different form of the Weyl scalars.

 

Although it is possible to perform both computations, automatically and step-by-step, departing from the tetrad (7), that tetrad and the corresponding Weyl scalars (8) have radicals, making the readability of the formulas at each step less clear. Both computations, can be presented in more readable form without radicals departing from the tetrad shown in the book, that is

e_[a, mu] = (Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = -1, (2, 1) = 0, (2, 2) = r/P(z, zb, u), (2, 3) = 0, (2, 4) = 0, (3, 1) = r/P(z, zb, u), (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = -1, (4, 4) = -H(X)}))

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

(9)

"IsTetrad(?)"

`Type of tetrad: `*null

 

true

(10)

The corresponding Weyl scalars free of radicals are

"WeylScalars(?)"

psi__0 = P(z, zb, u)*(2*(diff(P(z, zb, u), z))*(diff(H(X), z))+P(z, zb, u)*(diff(diff(H(X), z), z)))/r^2, psi__1 = -(1/2)*(-(diff(diff(H(X), r), z))*P(z, zb, u)^2*r+2*P(z, zb, u)^2*(diff(H(X), z))-(diff(P(z, zb, u), u))*(diff(P(z, zb, u), z))*r+(diff(diff(P(z, zb, u), u), z))*r*P(z, zb, u))/(r^2*P(z, zb, u)), psi__2 = -(1/6)*(-(diff(diff(H(X), r), r))*r^2+2*(diff(H(X), r))*r-2*(diff(P(z, zb, u), z))*(diff(P(z, zb, u), zb))+2*(diff(diff(P(z, zb, u), z), zb))*P(z, zb, u)-2*H(X))/r^2, psi__3 = 0, psi__4 = 0

(11)

So set this tetrad as the starting point

"Setup(?)"

[tetrad = {(1, 4) = -1, (2, 2) = r/P(z, zb, u), (3, 1) = r/P(z, zb, u), (4, 3) = -1, (4, 4) = -H(X)}]

(12)


All the transformations performed automatically, in one go

 

To arrive in one go, automatically, to a tetrad whose Weyl scalars are in canonical form as in (31), use the optional argument canonicalform:

T__5 := TransformTetrad(canonicalform)

WeylScalars(T__5)

psi__0 = 0, psi__1 = 0, psi__2 = -(1/6)*(-(diff(diff(H(X), r), r))*r^2+2*(diff(H(X), r))*r-2*(diff(P(z, zb, u), z))*(diff(P(z, zb, u), zb))+2*(diff(diff(P(z, zb, u), z), zb))*P(z, zb, u)-2*H(X))/r^2, psi__3 = 1, psi__4 = 0

(13)

Note the length of T__5

length(T__5)

58242

(14)

That length corresponds to several pages long. That happens frequently, you get Weyl scalars with a minimum of residual invariance, at the cost of a more complicated tetrad.

 

The transformations step-by-step leading to the same canonical form of the Weyl scalars

 

Step 0

 

As mentioned above, to apply the plan outlined by Chandrasekhar, the starting point needs to be a tetrad with `&Psi;__4` <> 0, not the case of (9), so in this step 0 we use a transformation of Class[1] making `&psi;__4` <> 0. This transformation introduces a complex parameter E and to get `&psi;__4` <> 0 any value of E suffices. We use E = 1:

TransformTetrad(nullrotationwithfixedl_)

Matrix(%id = 18446744078634914990)

(15)

"`T__0` := eval(?,E=1)"

Matrix(%id = 18446744078634940646)

(16)

Indeed, for this tetrad, `&Psi;__4` <> 0:

WeylScalars(T__0)[-1]

psi__4 = ((diff(diff(H(X), r), r))*r^2*P(z, zb, u)+P(z, zb, u)^3*(diff(diff(H(X), z), z))+2*(diff(P(z, zb, u), z))*(diff(H(X), z))*P(z, zb, u)^2+2*(diff(diff(H(X), r), z))*P(z, zb, u)^2*r+2*(diff(P(z, zb, u), z))*(diff(P(z, zb, u), zb))*P(z, zb, u)+2*(diff(P(z, zb, u), u))*(diff(P(z, zb, u), z))*r-4*P(z, zb, u)^2*(diff(H(X), z))-2*(diff(diff(P(z, zb, u), z), zb))*P(z, zb, u)^2-2*(diff(H(X), r))*P(z, zb, u)*r-2*(diff(diff(P(z, zb, u), u), z))*r*P(z, zb, u)+2*H(X)*P(z, zb, u))/(r^2*P(z, zb, u))

(17)

Step 1

Next is a transformation of Class__2 to make `&Psi;__0` = 0, that in the case of Petrov type II also implies on `&Psi;__1` = 0.According to the the help page TransformTetrad , this transformation introduces a parameter B that, according to the plan outlined by Chandrasekhar in Chapter 7 page 388, is one of the two identical roots (out of the four roots) of the principalpolynomial. To see the principal polynomial, or, directly, its roots you can use the PetrovType  command:

PetrovType(principalroots = 'R')

"II"

(18)

The first two are the same and equal to -1

R[1 .. 2]

[-1, -1]

(19)

So the transformed tetrad T__1 is

T__1 := eval(TransformTetrad(T__0, nullrotationwithfixedn_), B = -1)

Matrix(%id = 18446744078641721462)

(20)

Check this result and the corresponding Weyl scalars to verify that we now have `&Psi;__0` = 0 and `&Psi;__1` = 0

IsTetrad(T__1)

`Type of tetrad: `*null

 

true

(21)

WeylScalars(T__1)[1 .. 2]

psi__0 = 0, psi__1 = 0

(22)

Step 2

Next is a transformation of Class__1 that makes `&Psi;__4` = 0. This transformation introduces a parameter E, that according to Chandrasekhar's plan can be taken equal to one of the roots of Weyl scalar `&Psi;__4`that corresponds to the transformed tetrad. So we need to proceed in three steps:

a. 

transform the tetrad introducing a parameter E in the tetrad's components

b. 

compute the Weyl scalars for that transformed tetrad

c. 

take `&Psi;__4` = 0 and solve for E

d. 

apply the resulting value of E to the transformed tetrad obtained in step a.

 

a.Transform the tetrad and for simplicity take E real

T__2 := eval(TransformTetrad(T__1, nullrotationwithfixedl_), conjugate(E) = E)

Matrix(%id = 18446744078624751238)

(23)

"IsTetrad(?)"

`Type of tetrad: `*null

 

true

(24)

b. Compute `&Psi;__4` for this tetrad

simplify(WeylScalars(T__2)[-1])

psi__4 = (r^2*P(z, zb, u)*(E-1)^2*(diff(diff(H(X), r), r))-2*r*P(z, zb, u)^2*(E-1)*(diff(diff(H(X), r), z))+P(z, zb, u)^3*(diff(diff(H(X), z), z))-2*P(z, zb, u)^2*(E-1)^2*(diff(diff(P(z, zb, u), z), zb))+2*r*P(z, zb, u)*(E-1)*(diff(diff(P(z, zb, u), u), z))-2*r*P(z, zb, u)*(E-1)^2*(diff(H(X), r))+4*P(z, zb, u)^2*(E+(1/2)*(diff(P(z, zb, u), z))-1)*(diff(H(X), z))+2*((P(z, zb, u)*(E-1)*(diff(P(z, zb, u), zb))-(diff(P(z, zb, u), u))*r)*(diff(P(z, zb, u), z))+H(X)*P(z, zb, u)*(E-1))*(E-1))/(r^2*P(z, zb, u))

(25)

c. Solve `&Psi;__4` = 0 discarding the case E = 0 which implies on no transformation

simplify(solve({rhs(psi__4 = (r^2*P(z, zb, u)*(E-1)^2*(diff(diff(H(X), r), r))-2*r*P(z, zb, u)^2*(E-1)*(diff(diff(H(X), r), z))+P(z, zb, u)^3*(diff(diff(H(X), z), z))-2*P(z, zb, u)^2*(E-1)^2*(diff(diff(P(z, zb, u), z), zb))+2*r*P(z, zb, u)*(E-1)*(diff(diff(P(z, zb, u), u), z))-2*r*P(z, zb, u)*(E-1)^2*(diff(H(X), r))+4*P(z, zb, u)^2*(E+(1/2)*(diff(P(z, zb, u), z))-1)*(diff(H(X), z))+2*((P(z, zb, u)*(E-1)*(diff(P(z, zb, u), zb))-(diff(P(z, zb, u), u))*r)*(diff(P(z, zb, u), z))+H(X)*P(z, zb, u)*(E-1))*(E-1))/(r^2*P(z, zb, u))) = 0, E <> 0}, {E}, explicit)[1])

{E = ((diff(diff(H(X), r), r))*r^2*P(z, zb, u)+(diff(diff(H(X), r), z))*P(z, zb, u)^2*r-2*(diff(H(X), r))*P(z, zb, u)*r-(diff(diff(P(z, zb, u), u), z))*r*P(z, zb, u)+(diff(P(z, zb, u), u))*(diff(P(z, zb, u), z))*r-2*P(z, zb, u)^2*(diff(H(X), z))-2*(diff(diff(P(z, zb, u), z), zb))*P(z, zb, u)^2+2*(diff(P(z, zb, u), z))*(diff(P(z, zb, u), zb))*P(z, zb, u)+2*H(X)*P(z, zb, u)+(-P(z, zb, u)^4*((diff(diff(H(X), r), r))*r^2+2*(diff(P(z, zb, u), z))*(diff(P(z, zb, u), zb))-2*(diff(H(X), r))*r-2*(diff(diff(P(z, zb, u), z), zb))*P(z, zb, u)+2*H(X))*(diff(diff(H(X), z), z))+P(z, zb, u)^4*(diff(diff(H(X), r), z))^2*r^2+(-2*r^2*(diff(diff(P(z, zb, u), u), z))*P(z, zb, u)^3+2*r^2*P(z, zb, u)^2*(diff(P(z, zb, u), u))*(diff(P(z, zb, u), z))-4*r*(diff(H(X), z))*P(z, zb, u)^4)*(diff(diff(H(X), r), z))-2*(diff(P(z, zb, u), z))*(diff(H(X), z))*P(z, zb, u)^3*(diff(diff(H(X), r), r))*r^2+P(z, zb, u)^2*(diff(diff(P(z, zb, u), u), z))^2*r^2+(-2*r^2*P(z, zb, u)*(diff(P(z, zb, u), u))*(diff(P(z, zb, u), z))+4*r*(diff(H(X), z))*P(z, zb, u)^3)*(diff(diff(P(z, zb, u), u), z))+4*(diff(P(z, zb, u), z))*(diff(H(X), z))*P(z, zb, u)^4*(diff(diff(P(z, zb, u), z), zb))+4*(diff(H(X), z))^2*P(z, zb, u)^4+4*P(z, zb, u)^2*(diff(P(z, zb, u), z))*((diff(H(X), r))*P(z, zb, u)*r-(diff(P(z, zb, u), z))*(diff(P(z, zb, u), zb))*P(z, zb, u)-(diff(P(z, zb, u), u))*r-H(X)*P(z, zb, u))*(diff(H(X), z))+(diff(P(z, zb, u), u))^2*(diff(P(z, zb, u), z))^2*r^2)^(1/2))/(P(z, zb, u)*((diff(diff(H(X), r), r))*r^2+2*(diff(P(z, zb, u), z))*(diff(P(z, zb, u), zb))-2*(diff(H(X), r))*r-2*(diff(diff(P(z, zb, u), z), zb))*P(z, zb, u)+2*H(X)))}

(26)

d. Apply this result to the tetrad (23). In doing so, do not display the result, just measure its length (corresponds to two+ pages)

T__3 := simplify(eval(T__2, {E = ((diff(diff(H(X), r), r))*r^2*P(z, zb, u)+(diff(diff(H(X), r), z))*P(z, zb, u)^2*r-2*(diff(H(X), r))*P(z, zb, u)*r-(diff(diff(P(z, zb, u), u), z))*r*P(z, zb, u)+(diff(P(z, zb, u), u))*(diff(P(z, zb, u), z))*r-2*P(z, zb, u)^2*(diff(H(X), z))-2*(diff(diff(P(z, zb, u), z), zb))*P(z, zb, u)^2+2*(diff(P(z, zb, u), z))*(diff(P(z, zb, u), zb))*P(z, zb, u)+2*H(X)*P(z, zb, u)+(-P(z, zb, u)^4*((diff(diff(H(X), r), r))*r^2+2*(diff(P(z, zb, u), z))*(diff(P(z, zb, u), zb))-2*(diff(H(X), r))*r-2*(diff(diff(P(z, zb, u), z), zb))*P(z, zb, u)+2*H(X))*(diff(diff(H(X), z), z))+P(z, zb, u)^4*(diff(diff(H(X), r), z))^2*r^2+(-2*r^2*(diff(diff(P(z, zb, u), u), z))*P(z, zb, u)^3+2*r^2*P(z, zb, u)^2*(diff(P(z, zb, u), u))*(diff(P(z, zb, u), z))-4*r*(diff(H(X), z))*P(z, zb, u)^4)*(diff(diff(H(X), r), z))-2*(diff(P(z, zb, u), z))*(diff(H(X), z))*P(z, zb, u)^3*(diff(diff(H(X), r), r))*r^2+P(z, zb, u)^2*(diff(diff(P(z, zb, u), u), z))^2*r^2+(-2*r^2*P(z, zb, u)*(diff(P(z, zb, u), u))*(diff(P(z, zb, u), z))+4*r*(diff(H(X), z))*P(z, zb, u)^3)*(diff(diff(P(z, zb, u), u), z))+4*(diff(P(z, zb, u), z))*(diff(H(X), z))*P(z, zb, u)^4*(diff(diff(P(z, zb, u), z), zb))+4*(diff(H(X), z))^2*P(z, zb, u)^4+4*P(z, zb, u)^2*(diff(P(z, zb, u), z))*((diff(H(X), r))*P(z, zb, u)*r-(diff(P(z, zb, u), z))*(diff(P(z, zb, u), zb))*P(z, zb, u)-(diff(P(z, zb, u), u))*r-H(X)*P(z, zb, u))*(diff(H(X), z))+(diff(P(z, zb, u), u))^2*(diff(P(z, zb, u), z))^2*r^2)^(1/2))/(P(z, zb, u)*((diff(diff(H(X), r), r))*r^2+2*(diff(P(z, zb, u), z))*(diff(P(z, zb, u), zb))-2*(diff(H(X), r))*r-2*(diff(diff(P(z, zb, u), z), zb))*P(z, zb, u)+2*H(X)))}[1]))

length(T__3)

12589

(27)

Check the scalars, we expect `&Psi;__0` = `&Psi;__1` and `&Psi;__1` = `&Psi;__4` and `&Psi;__4` = 0

WeylScalars(T__3); %[1 .. 2], %[-1]

psi__0 = 0, psi__1 = 0, psi__4 = 0

(28)

Step 3

Use a transformation of Class[3] making Psi[3] = 1. Such a transformation changes Psi[3]^` '` = A*exp(-I*Omega)*Psi[3], where we need to take A*exp(-I*Omega) = 1/`&Psi;__3`, and without loss of generality we can take Omega = 0.

Check first the value of `&Psi;__3` in the last tetrad computed

WeylScalars(T__3)[4]

psi__3 = (1/2)*(-2*(diff(P(z, zb, u), z))*(diff(H(X), z))*P(z, zb, u)^3*(diff(diff(H(X), r), r))*r^2-P(z, zb, u)^4*(diff(diff(H(X), z), z))*(diff(diff(H(X), r), r))*r^2+P(z, zb, u)^4*(diff(diff(H(X), r), z))^2*r^2+4*(diff(P(z, zb, u), z))*(diff(H(X), r))*(diff(H(X), z))*P(z, zb, u)^3*r+2*(diff(H(X), r))*P(z, zb, u)^4*(diff(diff(H(X), z), z))*r-4*(diff(P(z, zb, u), zb))*(diff(P(z, zb, u), z))^2*(diff(H(X), z))*P(z, zb, u)^3+4*(diff(P(z, zb, u), z))*(diff(H(X), z))*P(z, zb, u)^4*(diff(diff(P(z, zb, u), z), zb))-4*(diff(H(X), z))*P(z, zb, u)^4*(diff(diff(H(X), r), z))*r+2*(diff(P(z, zb, u), u))*(diff(P(z, zb, u), z))*P(z, zb, u)^2*(diff(diff(H(X), r), z))*r^2-2*(diff(P(z, zb, u), zb))*(diff(P(z, zb, u), z))*P(z, zb, u)^4*(diff(diff(H(X), z), z))+2*P(z, zb, u)^5*(diff(diff(H(X), z), z))*(diff(diff(P(z, zb, u), z), zb))-2*P(z, zb, u)^3*(diff(diff(P(z, zb, u), u), z))*(diff(diff(H(X), r), z))*r^2+4*(diff(H(X), z))^2*P(z, zb, u)^4-4*(diff(P(z, zb, u), u))*(diff(P(z, zb, u), z))*(diff(H(X), z))*P(z, zb, u)^2*r-4*H(X)*(diff(P(z, zb, u), z))*(diff(H(X), z))*P(z, zb, u)^3+4*(diff(H(X), z))*P(z, zb, u)^3*(diff(diff(P(z, zb, u), u), z))*r+(diff(P(z, zb, u), u))^2*(diff(P(z, zb, u), z))^2*r^2-2*(diff(P(z, zb, u), u))*(diff(P(z, zb, u), z))*P(z, zb, u)*(diff(diff(P(z, zb, u), u), z))*r^2-2*H(X)*P(z, zb, u)^4*(diff(diff(H(X), z), z))+P(z, zb, u)^2*(diff(diff(P(z, zb, u), u), z))^2*r^2)^(1/2)/(r^2*P(z, zb, u))

(29)

So, the transformed tetrad T__4 to which corresponds Weyl scalars in canonical form, with `&Psi;__0` = `&Psi;__1` and `&Psi;__1` = `&Psi;__4` and `&Psi;__4` = 0 and `&Psi;__3` = 1, is

T__4 := simplify(eval(TransformTetrad(T__3, boostsn_l_plane), A = 1/rhs(psi__3 = (1/2)*(-2*(diff(P(z, zb, u), z))*(diff(H(X), z))*P(z, zb, u)^3*(diff(diff(H(X), r), r))*r^2-P(z, zb, u)^4*(diff(diff(H(X), z), z))*(diff(diff(H(X), r), r))*r^2+P(z, zb, u)^4*(diff(diff(H(X), r), z))^2*r^2+4*(diff(P(z, zb, u), z))*(diff(H(X), r))*(diff(H(X), z))*P(z, zb, u)^3*r+2*(diff(H(X), r))*P(z, zb, u)^4*(diff(diff(H(X), z), z))*r-4*(diff(P(z, zb, u), zb))*(diff(P(z, zb, u), z))^2*(diff(H(X), z))*P(z, zb, u)^3+4*(diff(P(z, zb, u), z))*(diff(H(X), z))*P(z, zb, u)^4*(diff(diff(P(z, zb, u), z), zb))-4*(diff(H(X), z))*P(z, zb, u)^4*(diff(diff(H(X), r), z))*r+2*(diff(P(z, zb, u), u))*(diff(P(z, zb, u), z))*P(z, zb, u)^2*(diff(diff(H(X), r), z))*r^2-2*(diff(P(z, zb, u), zb))*(diff(P(z, zb, u), z))*P(z, zb, u)^4*(diff(diff(H(X), z), z))+2*P(z, zb, u)^5*(diff(diff(H(X), z), z))*(diff(diff(P(z, zb, u), z), zb))-2*P(z, zb, u)^3*(diff(diff(P(z, zb, u), u), z))*(diff(diff(H(X), r), z))*r^2+4*(diff(H(X), z))^2*P(z, zb, u)^4-4*(diff(P(z, zb, u), u))*(diff(P(z, zb, u), z))*(diff(H(X), z))*P(z, zb, u)^2*r-4*H(X)*(diff(P(z, zb, u), z))*(diff(H(X), z))*P(z, zb, u)^3+4*(diff(H(X), z))*P(z, zb, u)^3*(diff(diff(P(z, zb, u), u), z))*r+(diff(P(z, zb, u), u))^2*(diff(P(z, zb, u), z))^2*r^2-2*(diff(P(z, zb, u), u))*(diff(P(z, zb, u), z))*P(z, zb, u)*(diff(diff(P(z, zb, u), u), z))*r^2-2*H(X)*P(z, zb, u)^4*(diff(diff(H(X), z), z))+P(z, zb, u)^2*(diff(diff(P(z, zb, u), u), z))^2*r^2)^(1/2)/(r^2*P(z, zb, u)))))

IsTetrad(T__4)

`Type of tetrad: `*null

 

true

(30)

WeylScalars(T__4)

psi__0 = 0, psi__1 = 0, psi__2 = -(1/6)*(-(diff(diff(H(X), r), r))*r^2+2*(diff(H(X), r))*r+2*(diff(diff(P(z, zb, u), z), zb))*P(z, zb, u)-2*(diff(P(z, zb, u), z))*(diff(P(z, zb, u), zb))-2*H(X))/r^2, psi__3 = 1, psi__4 = 0

(31)

These are the same scalars computed in one go in (13)

psi__0 = 0, psi__1 = 0, psi__2 = -(1/6)*(-(diff(diff(H(X), r), r))*r^2+2*(diff(H(X), r))*r-2*(diff(P(z, zb, u), z))*(diff(P(z, zb, u), zb))+2*(diff(diff(P(z, zb, u), z), zb))*P(z, zb, u)-2*H(X))/r^2, psi__3 = 1, psi__4 = 0

psi__0 = 0, psi__1 = 0, psi__2 = -(1/6)*(-(diff(diff(H(X), r), r))*r^2+2*(diff(H(X), r))*r-2*(diff(P(z, zb, u), z))*(diff(P(z, zb, u), zb))+2*(diff(diff(P(z, zb, u), z), zb))*P(z, zb, u)-2*H(X))/r^2, psi__3 = 1, psi__4 = 0

(32)

``


 

Download The_metric_[27_37_1]_in_canonical_form.mw

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

Hi,
Some people using the Windows platform have had problems installing MapleCloud packages, including the Maplesoft Physics Updates. This problem does not happen in Macintosh or Linux/Unix, also does not happen with all Windows computers but with some of them, and is not a problem of the MapleCloud packages themselves, but a problem of the installer of packages.

I understand that a solution to this problem will be presented within an upcoming Maple dot release.

Meantime, there is a solution by installing a helper library; after that, MapleCloud packages install without problems in all Windows machines. So whoever is having trouble installing MapleCloud packages in Windows and prefers not to wait until that dot release, instead wants to install this helper library, please email me at physics@maplesoft.com

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

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

``

 

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

 

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

 

Vector notation

 

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

with(Physics); with(Vectors)

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

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

(1)

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

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

(2)

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

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

(3)

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

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

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

(4)

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

 

Tensor notation

 

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

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

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

 

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

 

_______________________________________________________

 

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

 

_______________________________________________________

 

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

 

_______________________________________________________

(5)

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

X[j](t)

(X)[j](t)

(6)

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

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

(7)

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

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

(8)

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

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

`Defined objects with tensor properties`

 

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

(9)

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

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

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

(10)

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

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

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

(11)

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

TensorArray(d_(X[j]))

Array(%id = 18446744078354237310)

(12)

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

 

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

 

Geodesics(t)

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

(13)

 

Curvilinear coordinates

 

Vector notation

 

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

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

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

(14)

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

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

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

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

(15)

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

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

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

(16)

The second term is due to the dependency of `#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))` on the coordinate phi together with the chain rule diff(`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`(t), t) = (diff(`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`, phi))*(diff(phi(t), t)) and (diff(`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`, phi))*(diff(phi(t), t)) = `#mover(mi("&phi;",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,

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

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

(17)

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

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

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

(18)

 

Newton's 2nd law becomes

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

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

(19)

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

`~`[`=`]({coeffs(rhs(diff(diff(r_(t), t), t) = _rho(t)*(diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2)+_phi(t)*(2*(diff(rho(t), t))*(diff(phi(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k), [`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`(t), `#mover(mi("&phi;",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 `#mover(mi("F"),mo("&rarr;"))`(t) = 0, although diff(z(t), t, t) = 0, no acceleration in the `#mover(mi("k"),mo("&and;"))` direction, is naturally expected, the same cannot be said about the other two equations for diff(phi(t), t, t) and diff(rho(t), t, t). Those two equations are discussed below under Coriolis and Centripetal forces. The key observation at this point, however, is that the right-hand sides of both unexpected equations involve diff(phi(t), t), rotation around the z axis.

 

Tensor notation

 

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

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

_______________________________________________________

 

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

 

_______________________________________________________

 

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

 

_______________________________________________________

(22)

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

Geodesics(t)

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

(23)

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

 

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

1. 

Make rho be a function of t directly in the metric

2. 

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

3. 

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

 

For step 1, we have

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

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

(24)

Set this metric where `&equiv;`(rho, rho(t))

"Setup(?):"

_______________________________________________________

 

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

 

_______________________________________________________

 

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

 

_______________________________________________________

(25)

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

Geodesics(t, tensornotation)

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

(26)

Step 3: compute the components of this tensorial equation

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

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

(27)

These are the same equations (23).

 

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

Define(F[`~j`] = [`F__&rho;`(t), `F__&phi;`(t), F__z(t)])

`Defined objects with tensor properties`

 

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

(28)

 

Newton's 2nd law (19)

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

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

(29)

now using tensorial notation, becomes

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

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

(30)

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

Array(%id = 18446744078329063774)

(31)

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

 

Coriolis force and centripetal force

 

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

 

`#mover(mi("r"),mo("&rarr;"))`(t) = z(t)*`#mover(mi("k"),mo("&and;"))`+`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`(t)*rho(t)

 

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

 

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

 

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

 

To see that the equations

 

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

 

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

`#mover(mi("&omega;",fontstyle = "normal"),mo("&rarr;"))` = -(diff(phi(t), t))*`#mover(mi("k"),mo("&and;"))`

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

(32)

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

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

`#mover(mi("a"),mo("&rarr;"))`(t) = (diff(rho(t), t, t))*`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`(t)+rho(t)*(diff(phi(t), t, t))*`#mover(mi("&phi;",fontstyle = "normal"),mo("&and;"))`(t)+(diff(z(t), t, t))*`#mover(mi("k"),mo("&and;"))`

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

(33)

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

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

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

(34)

Third the centripetal force divided by the mass, defined by

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

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

(35)

Adding these three terms,

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

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

(36)

So that

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

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

(37)

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

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

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

(38)

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

0

(39)

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

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

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

(40)

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

 

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

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

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

(41)

resulting in the three equations

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

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

(42)

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

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

(43)

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

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

(44)

which are the equations returned by Geodesics in (23)

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

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

(45)

``

References

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


 

Download The_equations_of_motion_in_curvilinear_coordinates.mw

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



Maple_for_Beginners.mw

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


Vectors in Spherical Coordinates using Tensor Notation

Edgardo S. Cheb-Terrab1 and Pascal Szriftgiser2

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

(1) Maplesoft

 

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

 

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

 

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

 

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

with(Physics); with(Vectors)

Setup(dimension = 3, coordinates = cartesian, g_ = `+`, spacetimeindices = lowercaselatin)

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

 

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

 

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

 

_______________________________________________________

 

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

 

_______________________________________________________

 

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

 

(`Defined Pauli sigma matrices (Psigma): `*sigma[1]*`, `*sigma[2]*`, `)*sigma[3]

 

__________________________________________________

 

_______________________________________________________

(1)

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

 

 

In vector calculus, at the root of everything there is the line element `#mover(mi("dr"),mo("&rarr;"))`, which in Cartesian coordinates has the simple form

dr_ = _i*dx+_j*dy+_k*dz

dr_ = _i*dx+_j*dy+_k*dz

(1.1)

To compute the line element  `#mover(mi("dr"),mo("&rarr;"))` in spherical coordinates, the starting point is the transformation

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

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

(1.2)

Coordinates(S = [r, theta, phi])

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

(1.3)

Since in (dr_ = _i*dx+_j*dy+_k*dz)*[dx, dy, dz] are just symbols with no relationship to "[x,y,z],"start transforming these differentials using the chain rule, computing the Jacobian of the transformation (1.2). In this Jacobian J, the first line is "[(&PartialD;x)/(&PartialD;r)dr", "(&PartialD;x)/(&PartialD;theta)"`d&theta;`, "(&PartialD;x)/(&PartialD;phi)dphi]"

J := VectorCalculus:-Jacobian(map(rhs, [x = r*sin(theta)*cos(phi), y = r*sin(theta)*sin(phi), z = r*cos(theta)]), [S])

 

So in matrix notation,

Vector([dx, dy, dz]) = J.Vector([dr, dtheta, dphi])

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

(1.4)

To complete the computation of  `#mover(mi("dr"),mo("&rarr;"))` in spherical coordinates we can now use ChangeBasis , provided that next we substitute (1.4) in the result, expressing the abstract objects [dx, dy, dz] in terms of [dr, `d&theta;`, `d&phi;`].

 

In two steps:

lhs(dr_ = _i*dx+_j*dy+_k*dz) = ChangeBasis(rhs(dr_ = _i*dx+_j*dy+_k*dz), spherical)

dr_ = (dx*sin(theta)*cos(phi)+dy*sin(theta)*sin(phi)+dz*cos(theta))*_r+(dx*cos(phi)*cos(theta)+dy*sin(phi)*cos(theta)-dz*sin(theta))*_theta+(cos(phi)*dy-sin(phi)*dx)*_phi

(1.5)

The line element

"simplify(subs(convert(lhs(?) =~ rhs(?),set),dr_ = (dx*sin(theta)*cos(phi)+dy*sin(theta)*sin(phi)+dz*cos(theta))*_r+(dx*cos(phi)*cos(theta)+dy*sin(phi)*cos(theta)-dz*sin(theta))*_theta+(cos(phi)*dy-sin(phi)*dx)*_phi))"

dr_ = _phi*dphi*r*sin(theta)+_theta*dtheta*r+_r*dr

(1.6)

This result is important: it gives us the so-called scale factors, the key that connect 3D vectors with the related covariant and contravariant tensors in curvilinear coordinates. The scale factors are computed from (1.6) by taking the scalar product with each of the unit vectors [`#mover(mi("r"),mo("&and;"))`, `#mover(mi("&theta;",fontstyle = "normal"),mo("&and;"))`, `#mover(mi("&phi;",fontstyle = "normal"),mo("&and;"))`], then taking the coefficients of the differentials [dr, `d&theta;`, `d&phi;`] (just substitute them by the number 1)

h := subs(`~`[`=`]([dr, `d&theta;`, `d&phi;`], 1), [seq(rhs(dr_ = _phi*dphi*r*sin(theta)+_theta*dtheta*r+_r*dr).q, q = [`#mover(mi("r"),mo("&and;"))`, `#mover(mi("&theta;",fontstyle = "normal"),mo("&and;"))`, `#mover(mi("&phi;",fontstyle = "normal"),mo("&and;"))`])])

[1, r, r*sin(theta)]

(1.7)

The scale factors are relevant because the components of the 3D vector and the corresponding tensor are not the same in curvilinear coordinates. For instance, representing the differential of the coordinates as the tensor dS^j = [dr, `d&theta;`, `d&phi;`], we see that corresponding vector, the line element in spherical coordinates `#mover(mi("dS"),mo("&rarr;"))`, is not  constructed by directly equating its components to the components of dS^j = [dr, `d&theta;`, `d&phi;`], so  

 

 `#mover(mi("dS"),mo("&rarr;"))` <> `d&phi;`*`#mover(mi("&phi;",fontstyle = "normal"),mo("&and;"))`+dr*`#mover(mi("r"),mo("&and;"))`+`d&theta;`*`#mover(mi("&theta;",fontstyle = "normal"),mo("&and;"))` 

 

The vector `#mover(mi("dS"),mo("&rarr;"))` is constructed multiplying these contravariant components [dr, `d&theta;`, `d&phi;`] by the scaling factors, as

 

 `#mover(mi("dS"),mo("&rarr;"))` = `d&phi;`*`h__&phi;`*`#mover(mi("&phi;",fontstyle = "normal"),mo("&and;"))`+dr*h__r*`#mover(mi("r"),mo("&and;"))`+`d&theta;`*`h__&theta;`*`#mover(mi("&theta;",fontstyle = "normal"),mo("&and;"))` 

 

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

 

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

 

where on the right-hand side we see the contravariant components "A[]^(j)" and the scale-factors h[j]. Because the system is orthogonal, each vector component `#msub(mi("A",fontstyle = "normal"),mfenced(mi("j")))`satisfies

A__j = h[j]*A[`~j`]

 

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

 

LinearAlgebra[Norm](`#mover(mi("A"),mo("&rarr;"))`) = A[j]*A[`~j`]

it follows that,

A__j = A__j/h__j

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

 

• 

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

 

II. Transformation of contravariant and covariant tensors

 

 

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

Define(A__c[j], A__s[j])

`Defined objects with tensor properties`

 

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

(2.1)

Transformation rule for a contravariant tensor

 

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

 

"TensorArray(`A__c`[]^(j))=TransformCoordinates(tr,`A__s`[]^(j),[X],[S])"

where on the left-hand side we get, isolated, the three components of A in Cartesian coordinates, and on the right-hand side we transform the spherical components "`A__c`[]^(j)", from spherical S = (r, theta, phi) (4th argument) to Cartesian X = (x, y, z) (3rd argument), which according to the 5th bullet of TransformCoordinates  will result in a transformation expressed in terms of the old coordinates (here the spherical [S]). Expand things to make the comparison with (1.4) possible by eye

 

Vector[column](TensorArray(A__c[`~j`])) = TransformCoordinates(tr, A__s[`~j`], [X], [S], simplifier = expand)

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

(2.2)

We see that the transformation rule for a contravariant vector "`A__c`[]^(j)"is, indeed, as the transformation (1.4) for the differential of the coordinates.

Transformation rule for a covariant tensor

 

For the transformation rule for the components of a covariant tensor A__c[j], we know, by definition, that it is `#mrow(msub(mi("A"),mi("&mu;",fontstyle = "normal")),mo("&ApplyFunction;"),mfenced(mi("y")),mo("&equals;"),mfrac(mrow(mo("&PartialD;"),msup(mi("x"),mi("&nu;",fontstyle = "normal"))),mrow(mo("&PartialD;"),msup(mi("y"),mi("&mu;",fontstyle = "normal"))),linethickness = "1"),mo("&InvisibleTimes;"),mo("&InvisibleTimes;"),msub(mi("A"),mi("&nu;",fontstyle = "normal")),mfenced(mi("x")))`, so the same transformation rule for the gradient [`&PartialD;`[x], `&PartialD;`[y], `&PartialD;`[z]], where `&PartialD;`[x] = (proc (u) options operator, arrow; diff(u, x) end proc) and so on. We can experiment this by directly changing variables in the differential operators [`&PartialD;`[x], `&PartialD;`[y], `&PartialD;`[z]], for example

d_[x] = PDEtools:-dchange(tr, proc (u) options operator, arrow; diff(u, x) end proc, simplify)

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

(2.3)

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

H := simplify(LinearAlgebra:-Transpose(1/J))

Vector([d_[x], d_[y], d_[z]]) = H.Vector([d_[r], d_[theta], d_[phi]])

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

(2.4)

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

Vector[column](TensorArray(A__c[j])) = TransformCoordinates(tr, A__s[j], [X], [r, theta, phi], simplifier = expand)

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

(2.5)

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

 

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

Vector[column](TensorArray(A__s[j])) = TransformCoordinates(tr, A__c[j], [S], [X], simplifier = expand)

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

(2.6)

III. Deriving the transformation rule for the Gradient using TransformCoordinates

 

 

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

ON


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

(%Nabla = Nabla)(f(X))

%Nabla(f(X)) = (diff(f(X), x))*_i+(diff(f(X), y))*_j+(diff(f(X), z))*_k

(3.1)

(%Nabla = Nabla)(f(S))

%Nabla(f(S)) = (diff(f(S), r))*_r+(diff(f(S), theta))*_theta/r+(diff(f(S), phi))*_phi/(r*sin(theta))

(3.2)

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

 

Start changing the vector basis in the gradient (3.1)

lhs(%Nabla(f(X)) = (diff(f(X), x))*_i+(diff(f(X), y))*_j+(diff(f(X), z))*_k) = ChangeBasis(rhs(%Nabla(f(X)) = (diff(f(X), x))*_i+(diff(f(X), y))*_j+(diff(f(X), z))*_k), spherical)

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

(3.3)

By eye, we see that in this result the coefficients of [`#mover(mi("r"),mo("&and;"))`, `#mover(mi("&theta;",fontstyle = "normal"),mo("&and;"))`, `#mover(mi("&phi;",fontstyle = "normal"),mo("&and;"))`] are the three lines in the right-hand side of (2.6) after replacing the covariant components A__j by the derivatives of f with respect to the jth coordinate, here displayed using indexed notation due to using CompactDisplay

`~`[`=`]([A__s[1], A__s[2], A__s[3]], [diff(f(S), r), diff(f(S), theta), diff(f(S), phi)])

[A__s[1] = Physics:-Vectors:-diff(f(S), r), A__s[2] = Physics:-Vectors:-diff(f(S), theta), A__s[3] = Physics:-Vectors:-diff(f(S), phi)]

(3.4)

`~`[`=`]([A__c[1], A__c[2], A__c[3]], [diff(f(X), x), diff(f(X), y), diff(f(X), z)])

[A__c[1] = Physics:-Vectors:-diff(f(X), x), A__c[2] = Physics:-Vectors:-diff(f(X), y), A__c[3] = Physics:-Vectors:-diff(f(X), z)]

(3.5)

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

"subs([A__s[1] = Physics:-Vectors:-diff(f(S),r), A__s[2] = Physics:-Vectors:-diff(f(S),theta), A__s[3] = Physics:-Vectors:-diff(f(S),phi)],[A__c[1] = Physics:-Vectors:-diff(f(X),x), A__c[2] = Physics:-Vectors:-diff(f(X),y), A__c[3] = Physics:-Vectors:-diff(f(X),z)],?)"

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

(3.6)

"subs(convert(lhs(?) =~ rhs(?),set),%Nabla(f(X)) = (diff(f(X),x)*sin(theta)*cos(phi)+diff(f(X),y)*sin(theta)*sin(phi)+diff(f(X),z)*cos(theta))*_r+(diff(f(X),x)*cos(phi)*cos(theta)+diff(f(X),y)*sin(phi)*cos(theta)-diff(f(X),z)*sin(theta))*_theta+(-diff(f(X),x)*sin(phi)+cos(phi)*diff(f(X),y))*_phi)"

%Nabla(f(X)) = ((sin(theta)*cos(phi)*(diff(f(S), r))+cos(theta)*cos(phi)*(diff(f(S), theta))/r-sin(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(theta)*cos(phi)+(sin(theta)*sin(phi)*(diff(f(S), r))+cos(theta)*sin(phi)*(diff(f(S), theta))/r+cos(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(theta)*sin(phi)+(cos(theta)*(diff(f(S), r))-sin(theta)*(diff(f(S), theta))/r)*cos(theta))*_r+((sin(theta)*cos(phi)*(diff(f(S), r))+cos(theta)*cos(phi)*(diff(f(S), theta))/r-sin(phi)*(diff(f(S), phi))/(r*sin(theta)))*cos(phi)*cos(theta)+(sin(theta)*sin(phi)*(diff(f(S), r))+cos(theta)*sin(phi)*(diff(f(S), theta))/r+cos(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(phi)*cos(theta)-(cos(theta)*(diff(f(S), r))-sin(theta)*(diff(f(S), theta))/r)*sin(theta))*_theta+(-(sin(theta)*cos(phi)*(diff(f(S), r))+cos(theta)*cos(phi)*(diff(f(S), theta))/r-sin(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(phi)+cos(phi)*(sin(theta)*sin(phi)*(diff(f(S), r))+cos(theta)*sin(phi)*(diff(f(S), theta))/r+cos(phi)*(diff(f(S), phi))/(r*sin(theta))))*_phi

(3.7)

Simplifying, we arrive at (3.2)

(lhs = `@`(`@`(expand, simplify), rhs))(%Nabla(f(X)) = ((sin(theta)*cos(phi)*(diff(f(S), r))+cos(theta)*cos(phi)*(diff(f(S), theta))/r-sin(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(theta)*cos(phi)+(sin(theta)*sin(phi)*(diff(f(S), r))+cos(theta)*sin(phi)*(diff(f(S), theta))/r+cos(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(theta)*sin(phi)+(cos(theta)*(diff(f(S), r))-sin(theta)*(diff(f(S), theta))/r)*cos(theta))*_r+((sin(theta)*cos(phi)*(diff(f(S), r))+cos(theta)*cos(phi)*(diff(f(S), theta))/r-sin(phi)*(diff(f(S), phi))/(r*sin(theta)))*cos(phi)*cos(theta)+(sin(theta)*sin(phi)*(diff(f(S), r))+cos(theta)*sin(phi)*(diff(f(S), theta))/r+cos(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(phi)*cos(theta)-(cos(theta)*(diff(f(S), r))-sin(theta)*(diff(f(S), theta))/r)*sin(theta))*_theta+(-(sin(theta)*cos(phi)*(diff(f(S), r))+cos(theta)*cos(phi)*(diff(f(S), theta))/r-sin(phi)*(diff(f(S), phi))/(r*sin(theta)))*sin(phi)+cos(phi)*(sin(theta)*sin(phi)*(diff(f(S), r))+cos(theta)*sin(phi)*(diff(f(S), theta))/r+cos(phi)*(diff(f(S), phi))/(r*sin(theta))))*_phi)

%Nabla(f(X)) = (diff(f(S), r))*_r+(diff(f(S), theta))*_theta/r+(diff(f(S), phi))*_phi/(r*sin(theta))

(3.8)

%Nabla(f(S)) = (diff(f(S), r))*_r+(diff(f(S), theta))*_theta/r+(diff(f(S), phi))*_phi/(r*sin(theta))

%Nabla(f(S)) = (diff(f(S), r))*_r+(diff(f(S), theta))*_theta/r+(diff(f(S), phi))*_phi/(r*sin(theta))

(3.9)

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

 

 

• 

The Divergence

 

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

A__s_ := A__r(S)*_r+`A__&theta;`(S)*_theta+`A__&phi;`(S)*_phi

A__r(S)*_r+`A__&theta;`(S)*_theta+`A__&phi;`(S)*_phi

(4.1)

CompactDisplay(%)

` A__r`(S)*`will now be displayed as`*A__r

 

` A__&phi;`(S)*`will now be displayed as`*`A__&phi;`

 

` A__&theta;`(S)*`will now be displayed as`*`A__&theta;`

(4.2)

%Divergence(%A__s_) = Divergence(A__s_)

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

(4.3)

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

TransformCoordinates(tr, g_[j, k], [S], setmetric)

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

 

`Changing the differentiation variables used to compute the Christoffel symbols from `[x, y, z]*` to `[r, theta, phi]*` while the spacetime metric depends on `[r, theta]

 

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

 

_______________________________________________________

 

`Coordinates: `[r, theta, phi]*`. Signature: `(`+ + -`)

 

_______________________________________________________

 

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

 

_______________________________________________________

 

`Setting `*greek*` letters to represent `*space*` indices`

(4.4)

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

Riemann[nonzero]

Physics:-Riemann[a, b, c, d] = {}

(4.5)

Consider now the divergence of the contravariant "`A__s`[]^(j)"tensor, computed in tensor notation

CompactDisplay(A__s(S))

` A__s`(S)*`will now be displayed as`*A__s

(4.6)

D_[j](A__s[`~j`](S))

Physics:-D_[j](A__s[`~j`](S), [S])

(4.7)

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

D_[j](A__s[`~j`](S), [S]) = convert(D_[j](A__s[`~j`](S), [S]), d_)

Physics:-D_[j](A__s[`~j`](S), [S]) = Physics:-d_[j](A__s[`~j`](S), [S])+Physics:-Christoffel[`~j`, a, j]*A__s[`~a`](S)

(4.8)

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

%D_[j](%A__s[`~j`]) = SumOverRepeatedIndices(D_[j](A__s[`~j`](S), [S]))

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

(4.9)

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

[seq(A__s[Library:-Contravariant(j)](S) = Component(A__s_, j)/h[j], j = 1 .. 3)]

[A__s[`~1`](S) = A__r(S), A__s[`~2`](S) = `A__&theta;`(S)/r, A__s[`~3`](S) = `A__&phi;`(S)/(r*sin(theta))]

(4.10)

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

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

(4.11)

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

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

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

(4.12)
• 

The Curl

 

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

%Curl(%A__s_) = Curl(A__s_)

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

(4.13)

 

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

 

"`&epsilon;`[i,j,k] `&dtri;`[]^(j)(`A__s`[]^(k))"

 

But in a curvilinear system `&epsilon;`[i, j, k] is not a tensor, we need to use the non-Galilean form Epsilon[i, j, k] = sqrt(%g_[determinant])*`&epsilon;`[i, j, k], where %g_[determinant] is the determinant of the metric. Moreover, since the expression "Epsilon[i,j,k] `&dtri;`[]^(j)(`A__s`[]^(k))"has one free covariant index (the first one), to compare with the vectorial formula (4.12) this index also needs to be rewritten as a vector component as discussed at the end of Sec. I, using

A__j = A__j/h__j

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

Setup(levicivita = nongalilean)

[levicivita = nongalilean]

(4.14)

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

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

(4.15)

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

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

%Curl(%A__s_) = [(sin(theta)^3*(diff(A__s[`~3`](S), theta))*r^2+2*sin(theta)^2*cos(theta)*A__s[`~3`](S)*r^2-(diff(A__s[`~2`](S), phi))*sin(theta)*r^2)/(%h[1]*sin(theta)^2*r^2), (-sin(theta)^3*(diff(A__s[`~3`](S), r))*r^4-2*sin(theta)^3*A__s[`~3`](S)*r^3+(diff(A__s[`~1`](S), phi))*sin(theta)*r^2)/(%h[2]*sin(theta)^2*r^2), (sin(theta)^3*(diff(A__s[`~2`](S), r))*r^4+2*sin(theta)^3*A__s[`~2`](S)*r^3-sin(theta)^3*(diff(A__s[`~1`](S), theta))*r^2)/(%h[3]*sin(theta)^2*r^2)]

(4.16)

Replace now the components of the tensor "`A__s`[]^(j)" by the components of the 3D vector `#mover(mi("\`A__s\`"),mo("&rarr;"))` using (4.10)

lhs(%Curl(%A__s_) = [(sin(theta)^3*(diff(A__s[`~3`](S), theta))*r^2+2*sin(theta)^2*cos(theta)*A__s[`~3`](S)*r^2-(diff(A__s[`~2`](S), phi))*sin(theta)*r^2)/(%h[1]*sin(theta)^2*r^2), (-sin(theta)^3*(diff(A__s[`~3`](S), r))*r^4-2*sin(theta)^3*A__s[`~3`](S)*r^3+(diff(A__s[`~1`](S), phi))*sin(theta)*r^2)/(%h[2]*sin(theta)^2*r^2), (sin(theta)^3*(diff(A__s[`~2`](S), r))*r^4+2*sin(theta)^3*A__s[`~2`](S)*r^3-sin(theta)^3*(diff(A__s[`~1`](S), theta))*r^2)/(%h[3]*sin(theta)^2*r^2)]) = value(subs[eval]([A__s[`~1`](S) = A__r(S), A__s[`~2`](S) = `A__&theta;`(S)/r, A__s[`~3`](S) = `A__&phi;`(S)/(r*sin(theta))], rhs(%Curl(%A__s_) = [(sin(theta)^3*(diff(A__s[`~3`](S), theta))*r^2+2*sin(theta)^2*cos(theta)*A__s[`~3`](S)*r^2-(diff(A__s[`~2`](S), phi))*sin(theta)*r^2)/(%h[1]*sin(theta)^2*r^2), (-sin(theta)^3*(diff(A__s[`~3`](S), r))*r^4-2*sin(theta)^3*A__s[`~3`](S)*r^3+(diff(A__s[`~1`](S), phi))*sin(theta)*r^2)/(%h[2]*sin(theta)^2*r^2), (sin(theta)^3*(diff(A__s[`~2`](S), r))*r^4+2*sin(theta)^3*A__s[`~2`](S)*r^3-sin(theta)^3*(diff(A__s[`~1`](S), theta))*r^2)/(%h[3]*sin(theta)^2*r^2)])))

%Curl(%A__s_) = [(sin(theta)^3*((diff(`A__&phi;`(S), theta))/(r*sin(theta))-`A__&phi;`(S)*cos(theta)/(r*sin(theta)^2))*r^2+2*sin(theta)*cos(theta)*`A__&phi;`(S)*r-(diff(`A__&theta;`(S), phi))*r*sin(theta))/(h[1]*sin(theta)^2*r^2), (-sin(theta)^3*((diff(`A__&phi;`(S), r))/(r*sin(theta))-`A__&phi;`(S)/(r^2*sin(theta)))*r^4-2*sin(theta)^2*`A__&phi;`(S)*r^2+(diff(A__r(S), phi))*sin(theta)*r^2)/(h[2]*sin(theta)^2*r^2), (sin(theta)^3*((diff(`A__&theta;`(S), r))/r-`A__&theta;`(S)/r^2)*r^4+2*sin(theta)^3*`A__&theta;`(S)*r^2-sin(theta)^3*(diff(A__r(S), theta))*r^2)/(h[3]*sin(theta)^2*r^2)]

(4.17)

(lhs = `@`(simplify, rhs))(%Curl(%A__s_) = [(sin(theta)^3*((diff(`A__&phi;`(S), theta))/(r*sin(theta))-`A__&phi;`(S)*cos(theta)/(r*sin(theta)^2))*r^2+2*sin(theta)*cos(theta)*`A__&phi;`(S)*r-(diff(`A__&theta;`(S), phi))*r*sin(theta))/(h[1]*sin(theta)^2*r^2), (-sin(theta)^3*((diff(`A__&phi;`(S), r))/(r*sin(theta))-`A__&phi;`(S)/(r^2*sin(theta)))*r^4-2*sin(theta)^2*`A__&phi;`(S)*r^2+(diff(A__r(S), phi))*sin(theta)*r^2)/(h[2]*sin(theta)^2*r^2), (sin(theta)^3*((diff(`A__&theta;`(S), r))/r-`A__&theta;`(S)/r^2)*r^4+2*sin(theta)^3*`A__&theta;`(S)*r^2-sin(theta)^3*(diff(A__r(S), theta))*r^2)/(h[3]*sin(theta)^2*r^2)])

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

(4.18)

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

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

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

(4.19)
• 

The Gradient

 

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

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

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

(4.20)

where on the left-hand side we have the vectorial Nabla  differential operator and on the right-hand side, since f(S) is a scalar, the covariant derivative `&dtri;`[j](f) becomes the standard derivative `&PartialD;`[j](f).

lhs(%Nabla(f(S)) = Physics[d_][j](f(S), [S])/%h[j]) = eval(value(Library:-TensorComponents(rhs(%Nabla(f(S)) = Physics[d_][j](f(S), [S])/%h[j]))))

%Nabla(f(S)) = [Physics:-Vectors:-diff(f(S), r), (diff(f(S), theta))/r, (diff(f(S), phi))/(r*sin(theta))]

(4.21)

The above is the expected result (3.2)

%Nabla(f(S)) = (diff(f(S), r))*_r+(diff(f(S), theta))*_theta/r+(diff(f(S), phi))*_phi/(r*sin(theta))

%Nabla(f(S)) = (diff(f(S), r))*_r+(diff(f(S), theta))*_theta/r+(diff(f(S), phi))*_phi/(r*sin(theta))

(4.22)
• 

The Laplacian

 

Likewise we can compute the Laplacian directly as

%Laplacian(f(S)) = D_[j](D_[j](f(S)))

%Laplacian(f(S)) = Physics:-D_[j](Physics:-d_[`~j`](f(S), [S]), [S])

(4.23)

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

SumOverRepeatedIndices(%Laplacian(f(S)) = D_[j](Physics[d_][`~j`](f(S), [S]), [S]))

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

(4.24)

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

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

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

(4.25)

On the right-hand side we see the dAlembertian , "`&square;`(f(S)),"in curvilinear coordinates; rewrite it using standard diff  derivatives and expand both sides of the equation for comparison

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

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

(4.26)

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

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

true

(4.27)


 

Download Vectors_and_Spherical_coordinates_in_tensor_notation.mw

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

2 3 4 5 6 7 8 Last Page 4 of 17