Advanced Search

ecterrab

14702 Reputation

24 Badges

20 years, 240 days
Contact ecterrab

MaplePrimes Activity


These are replies submitted by ecterrab

  Hi Rouben, Having a ...

ecterrab 14702
July 26 2022
0 0

 

Hi Rouben,

Having a unique zero is an intentional design, not a mistake, not a bug nor something of the like.

 

First, note that the error message you show is due to the way you programmed your procedure f enforcing a type checking of its argument that only accepts a vector (what you'd call a non-zero vector). For example, if instead of that you accept input of type algebraic and use Vectors:-Norm, you have the f you seem to have tried to program

 

with(Physics:-Vectors)

 

f := proc (q::algebraic) options operator, arrow; Norm(q)^2 end proc

proc (q::algebraic) options operator, arrow; Physics:-Vectors:-Norm(q)^2 end proc

 (1)

f(a_)

Physics:-Vectors:-Norm(a_)^2

 (2)

f(k*a_)

Physics:-Vectors:-Norm(k*a_)^2

 (3)

expand(%)

k^2*Physics:-Vectors:-Norm(a_)^2

 (4)

f(0*a_)

0

 (5)

All the product and sum operations between objects of type algebraic are well defined with 0

0+a_

a_

 (6)

0.a_

0

 (7)

`&x`(0, a_)

0

 (8)

In this context, distinguishing between the algebraic 0 you see in (7) and the algebraic 0 you see in (8) is of no use, not for the purpose of performing algebraic computations (computing with objects of type algebraic). Such a distinction would then be artificial; an unnecessary complication.

 

Regarding your comment about other packages, VectorCalculus and LinearAlgebra compute with matrix representations of vectors, and there is indeed a computational distinction between the number 0 and a matrix all of whose components are equal to 0; so these two packages need to distinguish between the number 0 and a "0 vector" (represented as a column matrix of zeros). One of the main disadvantage of the matrix representation of vectors is that you cannot compute with the basic object A_, a non-projected vector, as you do when computing with Physics:-Vectors.

 

By the way, one of the key design differences between Physics and other packages is that, in Physics, instead of restricting the representations of mathematical objects to scalars or matrices, the domains of the operations were extended to handle the objects we use when computing with paper and pencil. That includes abstract vectors, tensors, Einstein's sum rule for repeated indices, anti-commutative and non-commutative objects, spinors, etc. All these objects, not present in the standard domain of computer algebra systems before Physics, are implemented with type algebraic. The number of advantages of having all these objects of type algebraic (when that is possible) is too large to list here, but for me the main one is that it allows for the closest match between what we do when computing with paper and pencil and what we can now do on a computer algebra worksheet.

 

If regardless of the above you still want to have a type that includes 0 as a vector you can use Or(0, PhysicsVectors). Then your procedure f will work as you planned, and in addition you can assign a name to this combo, as in TPV := Or(0, PhysicsVectors).

Download unique_zero.mw

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

 

Same procedure solves the problem...

ecterrab 14702
July 15 2022
0 0

@bikramphy 

You have now added information, the metric, OK. Below is basically the same worksheet, with small changes, taking into account the form of the metric.

 

with(Physics); with(Tetrads)

_______________________________________________________

 (1)

Simplify steps setting the signature as close as possible to the one you want, and set coordinates as you now indicated they are

Setup(signature = "+++-", coordinates = (X = [Q, q, u, t]))

[coordinatesystems = {X}, signature = `+ + + -`]

 (2)

 

Input the metric that you are now indicating

ds2 := dQ*dq+dQ*dq+du*dt+du*dt-Typesetting[delayDotProduct](2*Q/m.(diff(V(q, t), q)), dt, true)*dt

2*dQ*dq+2*du*dt-2*Q*(diff(V(q, t), q))*dt^2/m

 (3)

Set compact display to make things more readable 

CompactDisplay(2*dQ*dq+2*du*dt-2*Q*(diff(V(q, t), q))*dt^2/m)

V(q, t)*`will now be displayed as`*V

 (4)

Setup(metric = ds2)

[metric = {(1, 2) = 1, (3, 4) = 1, (4, 4) = -2*Q*(diff(V(q, t), q))/m}]

 (5)

The tetrad computed by default is orthonormal

e_[]

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

 (6)

"IsTetrad(?)"

true

 (7)

From

e_[definition]

Physics:-Tetrads:-e_[a, mu]*Physics:-Tetrads:-e_[b, `~mu`] = Physics:-Tetrads:-eta_[a, b]

 (8)

You see the form of the tetrad metric eta[a, b] that indicates the signature, is "(+++-)"

e_[a, mu]*e_[b, mu]; % = TensorArray(%, simplifier = `@`(simplify, evala))

Physics:-Tetrads:-e_[a, mu]*Physics:-Tetrads:-e_[b, `~mu`] = Matrix(%id = 36893488153619389964)

 (9)

You need to find a matrix A such that when redefining the tetrad as A * e_ = E, and substituting into the above, you get "E[a,mu] E[b]^(mu)=([[[1,0,0,0],[0,1,0,0],[0,0,-1,0],[0,0,0,-1]]])", where in the right-hand side you see the signature "(++--)".


I'd search for a matrix of this form, adding one element out of the diagonal, mostly for illustration purposes (you see in the result that was actually not necessary)

A := subs(`~`[`=`]([a[1, 2], a[1, 3], a[1, 4], a[2, 4], a[3, 4]], 0), Matrix(4, symbol = a, shape = symmetric))

Matrix(%id = 36893488151942187836)

 (10)

E[a, mu] = A.rhs(e_[])

E[a, mu] = Matrix(%id = 36893488151939865164)

 (11)

Define this as a tensor, to use the tensor functionality to build a set of equations for the `a__i,j`

"Define(?)"

{Physics:-D_[mu], Physics:-Dgamma[mu], E[a, mu], Physics:-Psigma[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-Tetrads:-e_[a, mu], Physics:-Tetrads:-eta_[a, b], Physics:-g_[mu, nu], Physics:-gamma_[i, j], Physics:-Tetrads:-gamma_[a, b, c], Physics:-Tetrads:-l_[mu], Physics:-Tetrads:-lambda_[a, b, c], Physics:-Tetrads:-m_[mu], Physics:-Tetrads:-mb_[mu], Physics:-Tetrads:-n_[mu], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

 (12)

Introduce the expected form of the tetrad metric to build the set of equations to be solved for the `a__i,j`

Eta[a, b] = Matrix(4, 4, [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, -1, 0], [0, 0, 0, -1]])

Eta[a, b] = Matrix(%id = 36893488153673893876)

 (13)

"Define(?)"

{Physics:-D_[mu], Physics:-Dgamma[mu], E[a, mu], Eta[a, b], Physics:-Psigma[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-Tetrads:-e_[a, mu], Physics:-Tetrads:-eta_[a, b], Physics:-g_[mu, nu], Physics:-gamma_[i, j], Physics:-Tetrads:-gamma_[a, b, c], Physics:-Tetrads:-l_[mu], Physics:-Tetrads:-lambda_[a, b, c], Physics:-Tetrads:-m_[mu], Physics:-Tetrads:-mb_[mu], Physics:-Tetrads:-n_[mu], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

 (14)

This is the identity you use to adjust the `a__i,j`

E[a, mu]*E[b, `~mu`] = Eta[a, b]

E[a, mu]*E[b, `~mu`] = Eta[a, b]

 (15)

TensorArray(E[a, mu]*E[b, `~mu`] = Eta[a, b], simplifier = simplify)

Matrix(%id = 36893488153624419132)

 (16)

"convert(?,setofequations)"

{0 = 0, a[1, 1]^2 = 1, a[2, 3]*(a[2, 2]+a[3, 3]) = 0, -a[4, 4]^2 = -1, a[2, 2]^2+a[2, 3]^2 = 1, a[2, 3]^2+a[3, 3]^2 = -1}

 (17)

solve({0 = 0, a[1, 1]^2 = 1, a[2, 3]*(a[2, 2]+a[3, 3]) = 0, -a[4, 4]^2 = -1, a[2, 2]^2+a[2, 3]^2 = 1, a[2, 3]^2+a[3, 3]^2 = -1})

{a[1, 1] = 1, a[2, 2] = 1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = 1}, {a[1, 1] = 1, a[2, 2] = 1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = -1}, {a[1, 1] = 1, a[2, 2] = -1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = 1}, {a[1, 1] = 1, a[2, 2] = -1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = -1}, {a[1, 1] = -1, a[2, 2] = 1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = 1}, {a[1, 1] = -1, a[2, 2] = 1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = -1}, {a[1, 1] = -1, a[2, 2] = -1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = 1}, {a[1, 1] = -1, a[2, 2] = -1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = -1}

 (18)

Any of these solutions serve your purpose, take the first one, remove that RootOf

sol := convert(({a[1, 1] = 1, a[2, 2] = 1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = 1}, {a[1, 1] = 1, a[2, 2] = 1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = -1}, {a[1, 1] = 1, a[2, 2] = -1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = 1}, {a[1, 1] = 1, a[2, 2] = -1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = -1}, {a[1, 1] = -1, a[2, 2] = 1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = 1}, {a[1, 1] = -1, a[2, 2] = 1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = -1}, {a[1, 1] = -1, a[2, 2] = -1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = 1}, {a[1, 1] = -1, a[2, 2] = -1, a[2, 3] = 0, a[3, 3] = RootOf(_Z^2+1), a[4, 4] = -1})[1], radical)

{a[1, 1] = 1, a[2, 2] = 1, a[2, 3] = 0, a[3, 3] = I, a[4, 4] = 1}

 (19)

The tetrad that results in the tetrad metric "(++--)"is then given by

subs(sol, E[])

E[a, mu] = Matrix(%id = 36893488153642495860)

 (20)

Set now this matrix as the tetrad

"Setup(tetrad = rhs(?))"

[tetrad = {(1, 1) = ((1/2)*I)*2^(1/2), (1, 2) = -((1/2)*I)*2^(1/2), (2, 3) = (1/2)*2^(1/2)/((Q^2*(diff(V(q, t), q))^2+m^2)^(1/2)/m)^(1/2), (2, 4) = -(1/2)*2^(1/2)*(Q*(diff(V(q, t), q))-(Q^2*(diff(V(q, t), q))^2+m^2)^(1/2))/(((Q^2*(diff(V(q, t), q))^2+m^2)^(1/2)/m)^(1/2)*m), (3, 3) = ((1/2)*I)*2^(1/2)/(-(Q^2*(diff(V(q, t), q))^2+m^2)^(1/2)/m)^(1/2), (3, 4) = -((1/2)*I)*2^(1/2)*(Q*(diff(V(q, t), q))+(Q^2*(diff(V(q, t), q))^2+m^2)^(1/2))/((-(Q^2*(diff(V(q, t), q))^2+m^2)^(1/2)/m)^(1/2)*m), (4, 1) = -((1/2)*I)*2^(1/2), (4, 2) = -((1/2)*I)*2^(1/2)}]

 (21)

And that is all.

 

Verify now:

e_[]

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

 (22)

This is the tetrad metric you were looking for, associated to a signature of the form "(++--)"

eta_[]

Physics:-Tetrads:-eta_[a, b] = Matrix(%id = 36893488151947672860)

 (23)

This is the tetrad definition

e_[definition]

Physics:-Tetrads:-e_[a, mu]*Physics:-Tetrads:-e_[b, `~mu`] = Physics:-Tetrads:-eta_[a, b]

 (24)

This definition is satisfied (by construction)

TensorArray(Physics:-Tetrads:-e_[a, mu]*Physics:-Tetrads:-e_[b, `~mu`] = Physics:-Tetrads:-eta_[a, b])

Matrix(%id = 36893488153642493340)

 (25)

This is the other way of the definition of the tetrad, relating it to the spacetime metric

e_[a, mu]*e_[a, nu] = g_[mu, nu]

Physics:-Tetrads:-e_[a, mu]*Physics:-Tetrads:-e_[`~a`, nu] = Physics:-g_[mu, nu]

 (26)

TensorArray(Physics:-Tetrads:-e_[a, mu]*Physics:-Tetrads:-e_[`~a`, nu] = Physics[g_][mu, nu])

Matrix(%id = 36893488151942188916)

 (27)

The equation components are all identities by eye, also by computer (useful when the expressions are difficult to read by eye).

"map(evalb,?)"

Matrix(%id = 36893488151918155164)

 (28)

Summarizing: the procedure is the same one shown in the first reply, just that instead of solving the system by eye as I did from (6) to (7) of the previous workshet-reply (something correct in that simple case of Schwarzschild solution) you solve it using the computer as done above from (15) to (18), also arriving at a simple solution.

 

NULL


 

Download arbitrary_signature_(II).mw
 

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

Hi, could you please upload a worksheet ...

ecterrab 14702
July 13 2022
0 0

Hi, could you please upload a worksheet with the metric? You can use the green arrow for that.

DiracComb...

ecterrab 14702
June 03 2022
0 0

Hi
Wikipedia mentions the DiracComb function. Anyone has an idea of how frequently used is this function? Maybe it is worth implementing it in Maple. Opinions?

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

May 27: Answer and add a worksheet ...

ecterrab 14702
May 27 2022
0 0

May 27: Answer and add a worksheet showing how to do this with Physics

First the answer to your question: change you input 

Oneform := RaiseLowerIndices(g, ON[1], [1])

by

Oneform := convert(RaiseLowerIndices(g, ON[1], [1]), DGform)

And that resolves the problem you encountered. Next is how you do the same computation using Physics

> 

NULL
> 

Delta := -2*M*R+R^2+a^2
> 

Sigma2 := (R^2+a^2)^2-Delta*a^2*sin(Theta)^2
> 

rho2 := R^2+a^2*cos(Theta)^2
> 

z := 2*M*R/rho2

  

Load Physics

> 

with(Physics)

  

Set the coordinates

> 

Coordinates(X = [T, R, Theta, Phi])

`Systems of spacetime coordinates are:`*{X = (T, R, Theta, Phi)}

  

{X}

 (1)

Set the metric, I am doing this in two steps so that you see how it is done: just remove the "&t" from around

> 

`#msup(mi("ds"),mn("2"))` := (-1+z)*dT^2+dT*dR+dT*dR+rho2*dTheta^2-z*a*sin(Theta)^2*(dPhi*dT+dPhi*dT)-a*sin(Theta)^2*(dR*dPhi+dR*dPhi)+Sigma2*sin(Theta)^2*dPhi^2/rho2

(-1+2*M*R/(R^2+a^2*cos(Theta)^2))*dT^2+2*dT*dR+(R^2+a^2*cos(Theta)^2)*dTheta^2-4*M*R*a*sin(Theta)^2*dPhi*dT/(R^2+a^2*cos(Theta)^2)-2*a*sin(Theta)^2*dR*dPhi+((R^2+a^2)^2-(-2*M*R+R^2+a^2)*a^2*sin(Theta)^2)*sin(Theta)^2*dPhi^2/(R^2+a^2*cos(Theta)^2)

 (2)

Set the metric and visually check it out (if its not right, adjust the line element above accordingly)

> 

Setup(metric = (-1+2*M*R/(R^2+a^2*cos(Theta)^2))*dT^2+2*dT*dR+(R^2+a^2*cos(Theta)^2)*dTheta^2-4*M*R*a*sin(Theta)^2*dPhi*dT/(R^2+a^2*cos(Theta)^2)-2*a*sin(Theta)^2*dR*dPhi+((R^2+a^2)^2-(-2*M*R+R^2+a^2)*a^2*sin(Theta)^2)*sin(Theta)^2*dPhi^2/(R^2+a^2*cos(Theta)^2))

_______________________________________________________

  

`Coordinates: `[T, R, Theta, Phi]*`. Signature: `(`- - - +`)

  

_______________________________________________________

  

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

  

_______________________________________________________

  

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

 (3)

To work with tetrads and everything related, there is a specialized package:

> 

with(Tetrads)

_______________________________________________________

  

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

  

(`Defined as tetrad tensors `(`see ?Physics,tetrads`)*`, `*`𝔢`[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 ?Physics,tetrads`)*`, `*l[mu]*` , `*n[mu]*` , `)*m[mu]*` , `*conjugate(m[mu])

  

_______________________________________________________

  

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

 (4)

For a brief description of each command, see the overview help page Tetrads

Set the signature, it follows the tetrad, by default it is of orthonormal type

> 

Setup(signature = "-+++")

[signature = `- + + +`]

 (5)
> 

e_[]

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

 (6)
> 

"IsTetrad(?)"

`Type of tetrad: `*orthonormal

  

true

 (7)

Check its definition

> 

e_[definition]

Physics:-Tetrads:-e_[a, mu]*Physics:-Tetrads:-e_[b, `~mu`] = Physics:-Tetrads:-eta_[a, b]

 (8)

Compute the components of this definition

> 

TensorArray(Physics:-Tetrads:-e_[a, mu]*Physics:-Tetrads:-e_[b, `~mu`] = Physics:-Tetrads:-eta_[a, b], simplifier = simplify)

Matrix(%id = 36893488159798381732)

 (9)

In the context of Physics you don't need to use a command to raise or lower an index, instead: just index with a covariant or contravariant index (the latter, by convention, are preceded by tilde ~). So to take one component of the tetrad, with either index covariant or contravariant, just input the tetrad with the desired indicesFor example, "`𝔢`[2]^(1) "

> 

"e_[~1,2]"

I*(2*a^2*cos(Theta)^2+2*R^2)^(1/2)/(-2*a^2*cos(Theta)^2+4*M*R-2*R^2)^(1/2)

 (10)

To compute the exterior derivative of an expression (see Physics:-ExteriorDerivative ), e.g. of  "`𝔢`[2]^(1) "

> 

"ExteriorDerivative[mu](e_[~1,2])"

I*(-(1/2)*(2*a^2*cos(Theta)^2+2*R^2)^(1/2)*(4*a^2*cos(Theta)*Physics:-d_[mu](Theta, [X])*sin(Theta)+4*M*Physics:-d_[mu](R, [X])-4*R*Physics:-d_[mu](R, [X]))/(-2*a^2*cos(Theta)^2+4*M*R-2*R^2)^(3/2)+(1/2)*(-4*a^2*cos(Theta)*Physics:-d_[mu](Theta, [X])*sin(Theta)+4*R*Physics:-d_[mu](R, [X]))/((-2*a^2*cos(Theta)^2+4*M*R-2*R^2)^(1/2)*(2*a^2*cos(Theta)^2+2*R^2)^(1/2)))

 (11)

To compute the ExteriorDerivative of a column of the tetrad, not just one component, do the same

> 

ExteriorDerivative[mu](e_[1, nu])

(1/2)*Physics:-d_[mu](Physics:-Tetrads:-e_[1, nu], [X])-(1/2)*Physics:-d_[nu](Physics:-Tetrads:-e_[1, mu], [X])

 (12)

To see the components of this expression always use TensorArray

> 

TensorArray((1/2)*Physics[d_][mu](Physics:-Tetrads:-e_[1, nu], [X])-(1/2)*Physics[d_][nu](Physics:-Tetrads:-e_[1, mu], [X]))

Matrix(%id = 36893488159805085084)

 (13)

You can also play around with these things, e.g. set a tensor whose components are given by (12)

> 

T[mu, nu] = (1/2)*Physics[d_][mu](Physics:-Tetrads:-e_[1, nu], [X])-(1/2)*Physics[d_][nu](Physics:-Tetrads:-e_[1, mu], [X])

T[mu, nu] = (1/2)*Physics:-d_[mu](Physics:-Tetrads:-e_[1, nu], [X])-(1/2)*Physics:-d_[nu](Physics:-Tetrads:-e_[1, mu], [X])

 (14)
> 

Define(T[mu, nu] = (1/2)*Physics[d_][mu](Physics:-Tetrads:-e_[1, nu], [X])-(1/2)*Physics[d_][nu](Physics:-Tetrads:-e_[1, mu], [X]))

`Defined objects with tensor properties`

  

{Physics:-D_[mu], Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], T[mu, nu], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-Tetrads:-e_[a, mu], Physics:-Tetrads:-eta_[a, b], Physics:-g_[mu, nu], Physics:-Tetrads:-gamma_[a, b, c], Physics:-gamma_[i, j], Physics:-Tetrads:-l_[mu], Physics:-Tetrads:-lambda_[a, b, c], Physics:-Tetrads:-m_[mu], Physics:-Tetrads:-mb_[mu], Physics:-Tetrads:-n_[mu], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

 (15)

Compute now any component of this tensor, e.g. T[`~1`, 2]

> 

"T[~1, 2]"

(4*I)*M*a^2*(-sin(Theta)^2*cos(Theta)^4*a^4-cos(2*Theta)*cos(Theta)^4*a^4+2*M*sin(Theta)^2*cos(Theta)^2*R*a^2+a^4*cos(Theta)^4+M*cos(2*Theta)*cos(Theta)^2*R*a^2-2*M*sin(Theta)^2*R^3-M*R*a^2*cos(Theta)^2+R^4*sin(Theta)^2-M*cos(2*Theta)*R^3+cos(2*Theta)*R^4+M*R^3-R^4)/((R^2+a^2*cos(Theta)^2)*(2*a^2*cos(Theta)^2+2*R^2)^(3/2)*(-2*a^2*cos(Theta)^2+4*M*R-2*R^2)^(3/2))

 (16)

 

Download WithPhysicsTetrads.mw

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

A comment and Maplesoft Physics Updates ...

ecterrab 14702
May 24 2022
0 0

In the latest Maplesoft Physics Updates (v.1235) there is a fix to "not placing the default nonnegint assumption" I mentioned a couple of replies above. So with the latest Physics Updates for Maple 2022, we now get

> D[1 $ j, 2 $ k](f)(x, y);  # computes both derivatives, x$j and y$k

                  /   k                                                 /

                  | -----                                               |

                  |  \                                      (2 j - _k1) |

pochhammer(-j, j) |   )    binomial(k, _k1) GAMMA(1 + 2 j) y            | ...

                  |  /                                                  |

                  | -----                                               |

                  |_k1 = 0                                              |

                  \                                                     \

                                     ...  

> eval(value(eval(%, [ j = 2, k = 4 ])), y = 0);

                                       6

 

PS: @Carl Love, although I see the motivation for your comment, I don't foresee adding to the diff command an option like "return a piecewise with all the cases".  The design of diff is to return a default value, and handle assumptions when possible - only that. It is similar to how the conversion network for mathematical functions works (see ?convert, to_special_function).

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

@vv I pointed out that the symbolic...

ecterrab 14702
May 23 2022
0 0

@vv 
I pointed out that the symbolic differentiation result, for symbols j and k, is correct. You objected to that. I only replied showing that, indeed, the result is correct. That closes the question.

Next you pointed out that the symbolic differentiation result is fragile. That concept is debatable. I do agree though, with a previous comment by you, that the most common situation is j, k nonnegint; so if choosing a way to express the result using GAMMA functions, do that for the most common case. The symbolic differentiation code is already doing that. But not in one branch of a subroutine, and the problem posted hits that branch. Although the result is correct, it can be better. It is something to adjust.

But then you come up with applying blindly the manipulation I used for the default output, after (you) using the opposite assumptions. With due respect, @vv, I don't see your point. You realize that if the result contains GAMMA as in this example, under direct evaluation it will either work with j negint, or for j posint, and not with both at the same time.

Add to that my opinion, that returning piecewise functions all around is not a convenient idea. Instead, I prefer presenting results with a convenient default and leave to the user to use assumptions for the other cases. Also relevant: understand that having Sums, and especially in the higher-order symbolic differentiation case, the results may require additional algebraic manipulations when evaluating the symbols involved at some values; typically, cases where correctness would require using limits or removing apparent singularities.

Returning to your example assuming j::nonegint, k::nonegint (above this reply), the manipulation after ex1 is as follows:

ex2 := simplify(eval(ex1, [j=2,k = 4, Sum = add]), GAMMA);
simplify(eval(ex2, [y=0]));
                                              6


So, the result by the symbolic differentiation code is also correct when assuming j::nonegint, k::nonegint.

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

@vv  It is about removable/...

ecterrab 14702
May 21 2022
0 0

@vv 

It is about removable/apparent singularities.

f := proc (x, y) options operator, arrow; 1/(2+x*y^2) end proc

(D[`$`(1, j), `$`(2, k)](f))(x, y)

pochhammer(-j, j)*(Sum(binomial(k, _k1)*(-1)^_k1*y^(-_k1)*(y^2)^j*GAMMA(-2*j+_k1)*(diff((x*y^2+2)^(-1-j), [`$`(y, k-_k1)]))/GAMMA(-2*j), _k1 = 0 .. k))

 (1)

We know that, for x = 0, y = 0, j = 2, k = 4 the value of (D[`$`(1, j), `$`(2, k)](f))(x, y)is 6, the issue is whether the answer above is or is not correct:

 

simplify(eval(pochhammer(-j, j)*(Sum(binomial(k, _k1)*(-1)^_k1*y^(-_k1)*(y^2)^j*GAMMA(-2*j+_k1)*(diff((x*y^2+2)^(-1-j), [`$`(y, k-_k1)]))/GAMMA(-2*j), _k1 = 0 .. k)), [k = 4, Sum = add]), GAMMA)

8*(x*y^2+2)^(-1-j)*(32*j^4+(-208*x*y^2-96)*j^3+(366*x^2*y^4+216*x*y^2+88)*j^2+(-180*x^3*y^6+18*x^2*y^4-56*x*y^2-24)*j+15*x^4*y^8-60*x^3*y^6+12*x^2*y^4)*pochhammer(-j, j)*(y^2)^j/((x*y^2+2)^4*y^4)

 (2)


eval(8*(x*y^2+2)^(-1-j)*(32*j^4+(-208*x*y^2-96)*j^3+(366*x^2*y^4+216*x*y^2+88)*j^2+(-180*x^3*y^6+18*x^2*y^4-56*x*y^2-24)*j+15*x^4*y^8-60*x^3*y^6+12*x^2*y^4)*pochhammer(-j, j)*(y^2)^j/((x*y^2+2)^4*y^4), [j = 2, x = 0])

6``

 (3)

NULL

Download Symbolic_differentation_is_correct.mw

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

in `eval/Sum`, not sum...

ecterrab 14702
May 21 2022
0 0

Hi @vv 

The problem you point at is in `eval/Sum`. Not in sum and not in D, both work correctly, that is why D[1$2,2$4](f)(0,0) returns 6, correct. To see the problem in `eval/Sum` try

f := (x, y) -> 1/(2 + x*y^2):
D[1 $ j, 2 $ k](f)(x, y);  # correct result for the jth and kth symbolic derivative

                                ...

eval(%, [x = 0, y = 0]);  # problem in `eval/Sum`
                                   /   k     \
                                   | -----   |
                                   |  \      |
                 pochhammer(-j, j) |   )    0|
                                   |  /      |
                                   | -----   |
                                   \_k1 = 0  /

The other questions...

ecterrab 14702
April 19 2022
0 0

 

Regarding your first question, the command to use is Library:-SubstituteOperator. Your reply indirectly introduces two more questions, one about reordering noncommuting operators, the command is SortProducts, and the other one is about commutativity of indexed quantum operators - use noncommutative variables instead, or quantum tensorial operators.

 

Finally, when you present a power as in A[1]*B[2]*B[2] = A[1]*B[2]^2, the code was requiring an extra argument to act on powers. I changed that default behavior so it now acts on powers automatically (unless you explicitly specify otherwise) - so for that you need to update to the latest Maplesoft Physics Updates version 1126. The rest of what follows doesn't require that.

 

Formulating the problem,

> 

with(Physics); Setup(quantumoperators = {A, B, C})

[quantumoperators = {A, B, C}]

 (1)
> 

Setup(algebrarules = {%Commutator(A[m], B[n]) = C[m+n]})

[algebrarules = {%Commutator(A[m], B[n]) = C[m+n]}]

 (2)

The answer to your original question is Library:-SubstituteOperator

> 

Library:-SubstituteOperator(A[1]*B[2] = 5, A[1]*B[2]*B[3]+C[3])

5*B[3]+C[3]

 (3)

Next question, when you are substituting into noncommutative operators, the order of the operands is relevant. For example, this doesn't work:

> 

('Library:-SubstituteOperator')(A[1]*B[2] = 5, A[1]*B[2]*B[3]+C[3])

Physics:-Library:-SubstituteOperator(Physics:-`*`(A[1], B[2]) = 5, Physics:-`*`(B[2], A[1], B[3])+C[3])

 (4)
> 

Physics[Library]:-SubstituteOperator(Physics[`*`](A[1], B[2]) = 5, Physics[`*`](B[2], A[1], B[3])+C[3])

Physics:-`*`(B[2], A[1], B[3])+C[3]

 (5)

Why? Because A[1]*B[2] <> A[1]*B[2]. So first you'd need to sort products, and in doing so any commutator rule is automatically taken into account. For example in view of

> 

(%Commutator = Commutator)(A[2], B[1])

%Commutator(A[2], B[1]) = C[3]

 (6)

The following two-step operation works

> 

SortProducts(Physics[`*`](B[2], A[1], B[3])+C[3], [A[1], B[2]])

Physics:-`*`(A[1], B[2], B[3])-Physics:-`*`(C[3], B[3])+C[3]

 (7)
> 

Library:-SubstituteOperator(A[1]*B[2] = 5, Physics[`*`](A[1], B[2], B[3])-Physics[`*`](C[3], B[3])+C[3])

5*B[3]-Physics:-`*`(C[3], B[3])+C[3]

 (8)

One could make these two steps be automatically be only one step, at the cost of having less control; something to think about.

 

Next, about commutation:

> 

Library:-Commute(B[1], B[2])

true

 (9)

Why is that? Because setting B as a quantum operator results in this action on kets.

> 

B[1].Ket(B, m, n)

m*Physics:-Ket(B, m, n)

 (10)
> 

B[2].Ket(B, m, n)

n*Physics:-Ket(B, m, n)

 (11)

That is, B[1] and B[2] are operators that act on disjointed spaces, the indexation is used to select the corresponding eigenvalue in a tensor product of states Ket(B, m, n), therefore B[1] and B[2] commute.

 

You could avoid this subtlety if you work with noncommutative variables as opposed to quantum operators. For example

> 

Setup(noncommutativeprefix = Z)

[noncommutativeprefix = {Z}]

 (12)
> 

Library:-Commute(Z[1], Z[2])

false

 (13)

You can also give a different meaning to quantum operators B[1] and B[2] and, in doing so, indicate they do not commute. Suppose you know the dimension of the index, for example 3; e.g. angular momentum. Define a tensor with its components here I use a different letter but any letter is as good

> 

Setup(spaceindices = lowercaselatin_is, op = L, tensors = L[m])

`* Partial match of '`*op*`' against keyword '`*quantumoperators*`' `

  

_______________________________________________________

  

[quantumoperators = {A, B, C, L}, spaceindices = lowercaselatin_is, tensors = {Physics:-Dgamma[mu], L[m], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-gamma_[i, j], Physics:-LeviCivita[alpha, beta, mu, nu]}]

 (14)

So now

> 

Library:-Commute(L[1], L[2])

false

 (15)

You see the dimension around:

> 

L[m]^2

Physics:-`*`(L[m], L[`~m`])

 (16)
> 

SumOverRepeatedIndices(Physics[`*`](L[m], L[`~m`]))

Physics:-`*`(L[1], L[`~1`])+Physics:-`*`(L[2], L[`~2`])+Physics:-`*`(L[3], L[`~3`])

 (17)

Now suppose you do not know the dimension. You can still define tensor of a space with unknown dimension, here generally called "gauge". So define corresponding indices and a tensorial quantum operator F[a] of that space

> 

Setup(gaugeindices = lowercase_ah, op = F, tensors = F[a])

`* Partial match of '`*op*`' against keyword '`*quantumoperators*`' `

  

_______________________________________________________

  

[gaugeindices = lowercaselatin_ah, quantumoperators = {A, B, C, F, L}, tensors = {Physics:-Dgamma[mu], F[a], L[m], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-gamma_[i, j], Physics:-LeviCivita[alpha, beta, mu, nu]}]

 (18)

Again, the components don't commute

> 

Library:-Commute(F[1], F[2])

false

 (19)

This time the dimension is not known, so for example the sum over the repeated indices cannot be performed:

> 

F[a]^2

Physics:-`^`(F[a], 2)

 (20)
> 

SumOverRepeatedIndices(Physics[`^`](F[a], 2))

Physics:-`^`(F[a], 2)

 (21)

Still, tensorial simplification operations continue working

> 

F[a]^2-F[b]^2

Physics:-`^`(F[a], 2)-Physics:-`^`(F[b], 2)

 (22)
> 

Simplify(Physics[`^`](F[a], 2)-Physics[`^`](F[b], 2))

0

 (23)
> 

KroneckerDelta[a, b]*F[b]

Physics:-KroneckerDelta[a, b]*F[b]

 (24)
> 

Simplify(Physics[KroneckerDelta][a, b]*F[b])

F[a]

 (25)

Now to handle powers

> 

A[1]*B[2]*B[2]-C[3]*B[3]+C[3]

Physics:-`*`(A[1], Physics:-`^`(B[2], 2))-Physics:-`*`(C[3], B[3])+C[3]

 (26)

without updating to the latest Physics Updates, you'd need to add the word 'true' at the end of the arguments passed. After installing the latest Physics updates that is not necessary:

> 

Library:-SubstituteOperator(A[1]*B[2] = 5, Physics[`*`](A[1], Physics[`^`](B[2], 2))-Physics[`*`](C[3], B[3])+C[3])

5*B[2]-Physics:-`*`(C[3], B[3])+C[3]

 (27)


PS: when using the Green Arrow to insert a worksheet, you can also ask for the contents to be shown, and filename.maple are not recommended for upload; use filename.mw instead.

Download SubstituteOperator_(II).mw

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

@etian2 then you don't have the...

ecterrab 14702
April 16 2022
0 0

@etian2 

But then you need to write the preamble properly - and that is not what you are doing. One workaround is for you to convert any worksheet, get the right preamble, then use it as a template. If that doesn't work AND you keep not posting a worksheet (do you expect people to start reading your image and trying to reproduce by typing etc.?), there is no chance.

Anyway, this is my last reply to this thread.

@etian2  I don"t see your work...

ecterrab 14702
April 16 2022
0 0

@etian2  I don"t see your worksheet. Also, the heading of the file you show is not the one you get ia File >Export> Latex. 

Post your worksheet please...

ecterrab 14702
April 15 2022
0 0

@etian2 

From your answers, guesswork is not helping. Let's move away from guessing. Could you please post the worksheet you are translating to LaTeX? I will give it a look.

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

maple.sty...

ecterrab 14702
April 15 2022
0 0

@etian2 

In usepackage, it is maple(.sty). Not maple2e(.sty). Yes you need to copy the sty files to the directory where you put the tex file, that is the standard LaTeX procedure for sty files.

@etian2  type a_ + b_ anywhere, sa...

ecterrab 14702
April 15 2022
0 0

@etian2 

type a_ + b_ anywhere, say within the text. Then mark the equation with the mouse, right-click, convert to > 2D-math non-executable, and you will see it with the arrow on top, as you say, symbols in nonexecutable math.

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

First 10 11 12 13 14 15 16 Last Page 12 of 65