ecterrab

13316 Reputation

24 Badges

18 years, 77 days

MaplePrimes Activity


These are answers submitted by ecterrab

Hi
When entering z1*z2*z1*z1; I receive z1 z2 z1^2 which is correct. What version of Maple are you using?

Edgardo S. Cheb-Terrab
Physics, Maplesoft


Hi Andriy,

restart; with(Physics)

Note first the use of the up-to-date Physics library (dowloadable from the Maple Physics: Research & Development updates page):

Physics:-Version()

"/Users/ecterrab/Maple/lib/Physics.mla", `2013, September 4, 15:20 hours`

(1)

 

Now, Psi is a mathematical function (check ?Psi) so, although possible, it's better to not use that name to avoid potential ambiguous meaning. I all what follows below I use psi instead.

Setup(mathematicalnotation = true, anticommutativeprefix = psi);

[anticommutativeprefix = {_lambda, psi}, mathematicalnotation = true]

(2)

Your annihilation operators: as mentioned by Carl, you can delay evaluation of these procedures am, ap, nn until i, j and sigma are all nonnegative integers. For that purpose, add the line you see below after "local k"

am := proc (i, j, sigma) local k; if [i, j, sigma]::(Not(list(nonnegint))) then return 'procname(args)' end if; k := 6*i-8+2*j+sigma; Annihilation(psi, k, notation = explicit) end proc;

proc (i, j, sigma) local k; if [i, j, sigma]::(Not(list(nonnegint))) then return 'procname(args)' end if; k := Physics:-`*`(6, i)-8+Physics:-`*`(2, j)+sigma; Physics:-Annihilation(psi, k, notation = explicit) end proc

(3)

Insert the same line for ap that generats your creation operators, and also for nn, the occupation number at state i, j, sigma

ap := proc (i, j, sigma) local k; if [i, j, sigma]::(Not(list(nonnegint))) then return 'procname(args)' end if; k := 6*i-8+2*j+sigma; Creation(psi, k, notation = explicit) end proc;

proc (i, j, sigma) local k; if [i, j, sigma]::(Not(list(nonnegint))) then return 'procname(args)' end if; k := Physics:-`*`(6, i)-8+Physics:-`*`(2, j)+sigma; Physics:-Creation(psi, k, notation = explicit) end proc

(4)

nn := proc (i, j, sigma) if [i, j, sigma]::(Not(list(nonnegint))) then return 'procname(args)' end if; ap(i, j, sigma).am(i, j, sigma) end proc

proc (i, j, sigma) if [i, j, sigma]::(Not(list(nonnegint))) then return 'procname(args)' end if; Typesetting:-delayDotProduct(ap(i, j, sigma), am(i, j, sigma)) end proc

(5)

IMPORTANT: the internal routines know that if A is a quantum operator (and all annihilation and creation operators generated by Annihilation and Creation are automatically quantum operators), a sum of A objects is also a quantum operator. However, because these am, ap and nn defined above are returning unevaluated when any of i, j, sigma are not numbers, the internal routines have no way to know these am, ap, nn programs are quantum operators, restricting the use of these three procedures when i, j, sigma are not numbers. You resolve this by directly passing this piece of information, as in

Setup(op = {am, ap, nn})

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

 

[quantumoperators = {am, ap, nn}]

(6)

And now you can operate with am, ap and nn as quantum operators even when i, j, sigma are not numbers. For instance, you can do sum, or Sum - but then you cannot us add for the summation where  f is a summation limit and is not defined. So for instance use sum for the addition from 1 to f, and use add for the rest

N := sum(add(add(nn(i, j, sigma), sigma = 1 .. 2), j = 1 .. 3), i = 1 .. f)

sum(nn(i, 1, 1)+nn(i, 1, 2)+nn(i, 2, 1)+nn(i, 2, 2)+nn(i, 3, 1)+nn(i, 3, 2), i = 1 .. f)

(7)

As said, because you indicated that nn is a quantum operator, the system now knows that this sum 'N ' is also a quantum operator:

type(N, Library:-PhysicsType:-ExtendedQuantumOperator)

true

(8)

You can now apply N to a Ket even not saying the value of f, although of course the operation cannot be performed to the end until you tell that value (note also that `.` automatically distributes over the sum)

K := Ket(psi, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0);

Physics:-Ket(psi, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0)

(9)

NK := Typesetting:-delayDotProduct(N, K);

sum(Physics:-`.`(nn(i, 1, 1), Physics:-Ket(psi, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0))+Physics:-`.`(nn(i, 1, 2), Physics:-Ket(psi, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0))+Physics:-`.`(nn(i, 2, 1), Physics:-Ket(psi, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0))+Physics:-`.`(nn(i, 2, 2), Physics:-Ket(psi, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0))+Physics:-`.`(nn(i, 3, 1), Physics:-Ket(psi, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0))+Physics:-`.`(nn(i, 3, 2), Physics:-Ket(psi, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0)), i = 1 .. f)

(10)

Evaluate for instance at  f = 2

eval(NK, f = 2)

5*Physics:-Ket(psi, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0)

(11)

This result is correct and I understand the answer to your question on how to do this, which is pretty close to what you have done in fact.

 

For curiosity, let's anyway verify this result. You need to check that


a) N is an operator that indeed returns the total number of particles.

b) This total number of particles of the state represented by the ket K is indeed 5.

 

Start with b) just because it is a 1 step thing: the total number of particles of a Ket is equal to the sum of its occupation numbers.  So, discarding the label of K (first operand), the occupation numbers are

op(2 .. -1, K)

1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0

(12)

and you can add them via

`+`(%)

5

(13)

So b) checks OK.

 

For a) you need to show that the summation NK indeed adds all these occupation numbers of K. For that, get the summand

summand := op(1, NK)

Physics:-`.`(nn(i, 1, 1), Physics:-Ket(psi, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0))+Physics:-`.`(nn(i, 1, 2), Physics:-Ket(psi, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0))+Physics:-`.`(nn(i, 2, 1), Physics:-Ket(psi, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0))+Physics:-`.`(nn(i, 2, 2), Physics:-Ket(psi, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0))+Physics:-`.`(nn(i, 3, 1), Physics:-Ket(psi, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0))+Physics:-`.`(nn(i, 3, 2), Physics:-Ket(psi, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0))

(14)

Get the nn operands in this summand

nn_operands := map2(op, 1, summand)

nn(i, 1, 1)+nn(i, 1, 2)+nn(i, 2, 1)+nn(i, 2, 2)+nn(i, 3, 1)+nn(i, 3, 2)

(15)

Check now the quantum number over which each of these Inn operators act. For that purpose note Library:-GetFingerprint (a section describing it is found within the help page for the Physics, Library )

am(1, 1, 4)

`a-`[psi[4]]

(16)

Library:-GetFingerprint(%)

"annihilation", "label" = psi, "quantum_numbers" = [4], "fermion"

(17)

The quantum number returned by nn(1, 1, 1).K is then the third element returned by GetFingerprint above. So get all the products am.ap being added and get the occupation number each of them is returning

nn_operands := map(op, [seq(nn_operands, i = 1 .. 2)])

[Physics:-`.`(`a+`[psi[1]], `a-`[psi[1]]), Physics:-`.`(`a+`[psi[2]], `a-`[psi[2]]), Physics:-`.`(`a+`[psi[3]], `a-`[psi[3]]), Physics:-`.`(`a+`[psi[4]], `a-`[psi[4]]), Physics:-`.`(`a+`[psi[5]], `a-`[psi[5]]), Physics:-`.`(`a+`[psi[6]], `a-`[psi[6]]), Physics:-`.`(`a+`[psi[7]], `a-`[psi[7]]), Physics:-`.`(`a+`[psi[8]], `a-`[psi[8]]), Physics:-`.`(`a+`[psi[9]], `a-`[psi[9]]), Physics:-`.`(`a+`[psi[10]], `a-`[psi[10]]), Physics:-`.`(`a+`[psi[11]], `a-`[psi[11]]), Physics:-`.`(`a+`[psi[12]], `a-`[psi[12]])]

(18)

map(proc (u) options operator, arrow; rhs(Library:-GetFingerprint(op(1, u))[3]) end proc, nn_operands)

[[1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12]]

(19)

In summary: the sum NK is returning a sum of each of the 12 occupation numbers and, indeed, twelve is the number of "occupation numbers" of K, so this sum is adding all the twelve occupation numbers:

nops(K)-1

12

(20)

and because the sum of these occupation numbers is also the total number of particles of the state represented by the Ket then NK is indeed returning the total number of particles.

 

``

Download MaplePrimesNK.mw

Edgardo S. Cheb-Terrab
Physics, Maplesoft

Hi Peter,

restart; with(Physics); with(Vectors)

[`&x`, `+`, `.`, ChangeBasis, Component, Curl, DirectionalDiff, Divergence, Gradient, Identify, Laplacian, Nabla, Norm, Setup, diff]

(1)

Note the use of the up-to-date Physics library (dowloadable from the Maple Physics: Research & Development updates page):

Physics:-Version()

"/Users/ecterrab/Maple/lib/Physics.mla", `2013, September 3, 19:11 hours`

(2)

Now regarding your post: to work in 3D, you can set the dimension to 3 as you did, yes, or perhaps simpler just set the metric of a the default 4D spacetime to have dt^2 = 1 and set the letters to use for spaceindices (not spacetimeindices), for instance as in

Setup(mathematicalnotation = true, coordinates = cylindrical, spaceindices = lowercaselatin_is, metric = drho^2+rho(t)^2*dphi^2+dz^2+dt^2, quiet)

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

 

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

 

[coordinatesystems = {X}, mathematicalnotation = true, metric = {(1, 1) = 1, (2, 2) = rho(t)^2, (3, 3) = 1, (4, 4) = 1}, spaceindices = lowercaselatin_is]

(3)

This is the metric (you do not need to use print as in print(g_[]))

g_[]

g[mu, nu] = (Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = rho(t)^2, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}))

(4)

You can now work with the 3D part by referring to it using spaceindices (not spacetimeindices).

 

Regarding your use of `.` instead of `*` in the Lagrangian below (I assume you did that to produce tensor contractions), note that the identity (X[]^(j))(t)=g[]^(i,j) (X[i])(t), valid in Euclidean or Minkowski spaces, is not valid in this case where diff and g[mu, nu] do not commute because the Christoffel symbols have nonzero components

Christoffel[nonzero]

Physics:-Christoffel[mu, nu, alpha] = {(1, 2, 2) = -rho(t), (2, 1, 2) = rho(t), (2, 2, 1) = rho(t)}

(5)

In short, in this system of coordinates (diff(X(t), t))^i is not a tensor

diff(X[`~mu`](t), t)

diff((Physics:-SpaceTimeVector[`~mu`](X))(t), t)

(6)

Library:-IsTensorInCurrentCoordinates(diff((Physics:-SpaceTimeVector[`~mu`](X))(t), t))

false

(7)

The use of `.` to contract the indices then produces no contraction (it would be wrong, equivalent to commuting diff and g[mu, nu]) and here is then the Lagrangian - the 3D part using spaceindices

L := (1/2)*m*g_[i, j].(diff(X[`~i`](t), t)).(diff(X[`~j`](t), t))

(1/2)*m*Physics:-g_[i, j]*(diff((Physics:-SpaceTimeVector[`~i`](X))(t), t))*(diff((Physics:-SpaceTimeVector[`~j`](X))(t), t))

(8)

Regarding your question: the command to perform a summation over the repeated indices is:

SumOverRepeatedIndices(L)

(1/2)*m*(diff(rho(t), t))^2+(1/2)*m*rho(t)^2*(diff(phi(t), t))^2+(1/2)*m*(diff(z(t), t))^2

(9)

Below you have a section showing the same but setting the dimention of 'spacetime' to 3 as you did originally.

 

Regarding your other question, on how to expand the vector expression

B0_.`&x`(v_, B0_+B1_)

Physics:-Vectors:-`.`(B0_, Physics:-Vectors:-`&x`(v_, B0_+B1_))

(10)

recalling I am working with the updated Physics library, you can use expand, Expand, or - in the case of Vectors expandable commands as this one - also with Simplify

expand(Physics:-Vectors:-`.`(B0_, Physics:-Vectors:-`&x`(v_, B0_+B1_)))

Physics:-Vectors:-`.`(Physics:-Vectors:-`&x`(v_, B1_), B0_)

(11)

Expand(Physics:-Vectors:-`.`(B0_, Physics:-Vectors:-`&x`(v_, B0_+B1_)))

Physics:-Vectors:-`.`(Physics:-Vectors:-`&x`(v_, B1_), B0_)

(12)

Simplify(Physics:-Vectors:-`.`(B0_, Physics:-Vectors:-`&x`(v_, B0_+B1_)))

Physics:-Vectors:-`.`(Physics:-Vectors:-`&x`(v_, B1_), B0_)

(13)

Alternative way, setting the dimension of "spacetime" to 3

 

Download MaplePrimesPhysicsVe.mw

 
Edgardo S. Cheb-Terrab
Physics, Maplesoft 

In brief: dsolve is a general solver, able to solve different kinds of equations, and with a strategy of taking advantage of different methods when more than one can solve a given equation - nothing of that is present in DEtools[bernoullisols] that can only solve Bernoulli equations and always using the same method.

 

Edgardo S. Cheb-Terrab
Physics, Maplesoft



Hi,

In what follows you have the symbolic exact solution to the ODE system you posted, matching the ics you posted, and without specializing any of the symbolic parameters a1, a2, a3, b1, b2 or b3.

 

Your ODE system:

dC1 := diff(C1(t), t) = -a1*C1(t)+b1*C2(t);

diff(C1(t), t) = -a1*C1(t)+b1*C2(t)

 

diff(C2(t), t) = a1*C1(t)-(b1+a2)*C2(t)+b2*C3(t)

 

diff(C3(t), t) = a2*C2(t)-(b2+a3)*C3(t)+b3*O1(t)

 

diff(O1(t), t) = a3*C3(t)-b3*O1(t)

 

diff(C1(t), t) = -a1*C1(t)+b1*C2(t), diff(C2(t), t) = a1*C1(t)-(b1+a2)*C2(t)+b2*C3(t), diff(C3(t), t) = a2*C2(t)-(b2+a3)*C3(t)+b3*O1(t), diff(O1(t), t) = a3*C3(t)-b3*O1(t)

 

C1(0) = 1, C2(0) = 0, C3(0) = 0, C1(t)+C2(t)+C3(t)+O1(t) = 1

 

[C1(t), C2(t), C3(t), O1(t)]

(1)

Strategy:

1) rewrite this system as an ODE system for three unknowns C1, C2, C3, plus an expression giving O1 in terms of them

2) compute the exact solution for C1, C2, C3 in terms of three arbitrary constants

3) adjust the three arbitrary constants to match your initial conditions ics for C1, C2 and C3

4) substitute this solution for C1, C2, C3 matching your initial conditions into the expression giving O1 in terms of them.

 

While performing these steps, do not display expressions that are too large.

 

Step 1), one of the equations you present as ics is actually a constraint between the functions:

ics[-1]

C1(t)+C2(t)+C3(t)+O1(t) = 1

(2)

The actual system to be solved then involves 4 DEs and one constraint

sys__0 := [syst, ics[-1]]:

nops(sys__0)

5

(3)

indets(sys__0, unknown)

{C1(t), C2(t), C3(t), O1(t)}

(4)

Rewrite now this as "one constraint and three DEs involving only C1, C2 and C3" (see casesplit )

PDEtools:-casesplit(sys__0, [O1, {C1, C2, C3}]);

`casesplit/ans`([O1(t) = -C1(t)-C2(t)-C3(t)+1, diff(C1(t), t) = -a1*C1(t)+b1*C2(t), diff(C2(t), t) = a1*C1(t)-b1*C2(t)-a2*C2(t)+b2*C3(t), diff(C3(t), t) = -b3*C1(t)-b3*C2(t)+a2*C2(t)-a3*C3(t)-b3*C3(t)-b2*C3(t)+b3], [])

(5)

In the above you see a system entirely equivalent to yours, but now involving 3 DEs instead of four. Remove then the constraint (that has O1 for lhs) to work just with the ODE system, and assign it to sys__1

sol__O1, sys__1 := selectremove(has, {op(op(1, `casesplit/ans`([O1(t) = -C1(t)-C2(t)-C3(t)+1, diff(C1(t), t) = -a1*C1(t)+b1*C2(t), diff(C2(t), t) = a1*C1(t)-b1*C2(t)-a2*C2(t)+b2*C3(t), diff(C3(t), t) = -b3*C1(t)-b3*C2(t)+a2*C2(t)-a3*C3(t)-b3*C3(t)-b2*C3(t)+b3], [])))}, O1)

{O1(t) = -C1(t)-C2(t)-C3(t)+1}, {diff(C1(t), t) = -a1*C1(t)+b1*C2(t), diff(C2(t), t) = a1*C1(t)-b1*C2(t)-a2*C2(t)+b2*C3(t), diff(C3(t), t) = -b3*C1(t)-b3*C2(t)+a2*C2(t)-a3*C3(t)-b3*C3(t)-b2*C3(t)+b3}

(6)

nops(sys__1)

3

(7)

Step 2) Compute now the exact solution to this system; note that dsolve will solve and simplify the solution - most of the time is spent in simplifying

_t0 := time():

timeconsumed = 583.306

(8)

nops(sol)

3

(9)

As mentioned by others, these solutions involve large algebraic expressions

map(length, sol)

{42200, 79352, 124604}

(10)

You can transform these into slightly shorter expressions using simplify/size

sol__1 := simplify(sol, size):

map(length, sol__1)

{41094, 64662, 84377}

(12)

Before proceeding further, let's verify this solution

odetest(sol__1, sys__1)

{0}

(13)

Step 3) Now you need to adjust the integration constants to match your ics for C1, C2 and C3

ics__1 := remove(has, [ics], O1)

[C1(0) = 1, C2(0) = 0, C3(0) = 0]

(14)

 are:

Cn := indets(`ics__1 _at_0`, suffixed(_C))

{_C1, _C2, _C3}

(15)

So: evaluate sol__1 at t = 0 and substitute in ics__1 in order to have a system for the Cn

`ics__1 _at_0` := subs(eval(sol__1, t = 0), ics__1):

Check lengths

map(length, `ics__1 _at_0`)

[350, 23917, 43632]

(16)

The symbols in these ics are

indets(`ics__1 _at_0`, symbol)

{_C1, _C2, _C3, a1, a2, a3, b1, b2, b3}

(17)

Solve now for the Cn and simplify the size; check lengths and left-hand sides

sol__Cn := simplify(solve(`ics__1 _at_0`, Cn), size):

map(length, sol__Cn)

{39639, 116164, 142373}

(18)

map(lhs, sol__Cn)

{_C1, _C2, _C3}

(19)

At this point we have solved the problem: sol__1 is the general solution in terms of three arbitrary constants Cn and sol__Cn gives the value of these three arbitrary constants that match your ics__1 for C1, C2 and C3. So evaluate sol__1 at sol__Cn to get the answer you are looking for; check lengths

answer__0 := eval(sol__1, sol__Cn):

map(lhs, answer__0)

{C1(t), C2(t), C3(t)}

(20)

map(length, answer__0)

{339223, 362791, 382506}

(21)

Step 4) To get the answer to O1, substitute these solutions for C1, C2 and C3 in the sol__O1

sol__O1

{O1(t) = -C1(t)-C2(t)-C3(t)+1}

(22)

answer__O1 := subs(answer__0, sol__O1):

length(answer__O1)

1084508

(23)

In addition, the following three would help in compacting the expressions, but they are rather large and simplifying takes time - I am keeping the lines here in comments in case you want to try them

NULL

NULL



Download MaplePrimesODESyste.mw

Edgardo S. Cheb-Terrab
Physics, Maplesoft 

Hi,

Suppose a set of N vectors where there may be one of them parallel to one another: their cross product is zero. So, in Maple, for j to N-1, for k from j+1 to N, if V[j] &x V[k] = 0 then V[j] is parallel to V[k]. You can perform this loop operation either using VectorCalculus or Physics:-Vectors.

Edgardo S. Cheb-Terrab 
Physics, Maplesoft

 

Hi,
Yes, in Maple 17 you can do substitutions of tensors of the form you mention, taking into account free and dummy indices the same way we do when computing with paper and pencil. You need to update your Physics library with the latest one, available for download at the Maplesoft Physics: Research & Development webpage. The command you use to do these substitutions, SubstituteTensor, is currently a Physics:-Library command, part of an upcomming more general Physics:-Substitute command, to handle also tensorial subexpressions within noncommutative products similar to what algsubs and simplify/siderels do with commutative nontensorial products.

 

with(Physics); with(Library)

Define(A, B, C, F, G)

`Defined objects with tensor properties`

 

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

(1)

A[mu] = Physics:-`*`(Physics:-`*`(G[nu, alpha], A[`~alpha`]), F[mu, `~nu`])

A[mu] = G[nu, alpha]*A[`~alpha`]*F[mu, `~nu`]

(2)

The repeated and free indices of the lhs and rhs of (2)

Check(A[mu] = G[nu, alpha]*A[`~alpha`]*F[mu, `~nu`], all)

`The repeated indices per term are: `[{`...`}, {`...`}, `...`]*`; the free indices are: `*{`...`}

 

([{}], {mu}) = ([{alpha, nu}], {mu})

(3)

The easy case

SubstituteTensor(A[mu] = G[nu, alpha]*A[`~alpha`]*F[mu, `~nu`], A[mu])

G[nu, alpha]*A[`~alpha`]*F[mu, `~nu`]

(4)

The free index in the target expression is not mu but nu

SubstituteTensor(A[mu] = G[nu, alpha]*A[`~alpha`]*F[mu, `~nu`], A[nu])

G[beta, alpha]*A[`~alpha`]*F[nu, `~beta`]

(5)

Distinction between covariant and contravariant indices

SubstituteTensor(A[mu] = G[nu, alpha]*A[`~alpha`]*F[mu, `~nu`], A[`~nu`])

G[beta, alpha]*A[`~alpha`]*F[`~nu`, `~beta`]

(6)

The index nu found repeated in the rhs of the substitution equation (2) apears also repeated in the following target expression

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

A[nu]*A[`~nu`]

(7)

SubstituteTensor(A[mu] = G[nu, alpha]*A[`~alpha`]*F[mu, `~nu`], A[nu]*A[`~nu`])

G[beta, alpha]*A[`~alpha`]*F[nu, `~beta`]*G[lambda, kappa]*A[`~kappa`]*F[`~nu`, `~lambda`]

(8)

Check(%, all)

`The repeated indices per term are: `[{`...`}, {`...`}, `...`]*`; the free indices are: `*{`...`}

 

[{alpha, beta, kappa, lambda, nu}], {}

(9)

Substitute (2) in the rhs of (2) itself

SubstituteTensor(A[mu] = G[nu, alpha]*A[`~alpha`]*F[mu, `~nu`], rhs(A[mu] = G[nu, alpha]*A[`~alpha`]*F[mu, `~nu`]))

G[nu, alpha]*G[kappa, beta]*A[`~beta`]*F[`~alpha`, `~kappa`]*F[mu, `~nu`]

(10)

Check(%, all)

`The repeated indices per term are: `[{`...`}, {`...`}, `...`]*`; the free indices are: `*{`...`}

 

[{alpha, beta, kappa, nu}], {mu}

(11)

The example you posted as a question

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

A[mu]*B[`~mu`]

(12)

A[alpha] = Physics:-`*`(Physics:-`*`(C[alpha], C[mu]), C[`~mu`])

A[alpha] = C[alpha]*C[mu]*C[`~mu`]

(13)

SubstituteTensor(A[alpha] = C[alpha]*C[mu]*C[`~mu`], A[mu]*B[`~mu`])

C[mu]*C[nu]*C[`~nu`]*B[`~mu`]

(14)

``


Download SubstituteTensor.mw

Edgardo S. Cheb-Terrab
Physics, Maplesoft

Hi
allvalues is the natural command to use in that - when solve can rewrite the RootOf expression free of RootOf - you will get rid of it.

There are other situations though, where you would like to understand the RootOf expression without RootOfs around, or allvalues is of no help because solve cannot rewrite the expression without them. This is a typical situation when working with differential equations. A command helpful in these cases, not so well known, is DEtools[remove_RootOf].

Edgardo S. Cheb-Terrab
Physics, Maplesoft

You can try other commands of PDEtools as PolynomialSolutions, FunctionFieldSolutions and mainly InvariantSolutions combined with its optional arguments. It would also help seeing the PDE you want to solve - sometimes one can see a change of variables that helps.

Edgardo S. Cheb-Terrab
Physics, Maplesoft


Hi. My suggestion was for you to do these computations using the Physics package, not the old 'tensor' package. Here it is how.

with(Physics):

Example with the Schwarzschild metric in spherical coordinates: instead of using Setup (you can, but you do not need to), set the spacetime metric directly from the the metric command g_

g_[sc]

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

 

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

 

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

 

`Parameters: `[m]

 

g[mu, nu] = (Matrix(4, 4, {(1, 1) = r/(-r+2*m), (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = -r^2, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -r^2*sin(theta)^2, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = (r-2*m)/r}))

(1)

The Kretschmann scalar is given by

K := Riemann[alpha, beta, gamma, delta]^2

Physics:-Riemann[alpha, beta, delta, gamma]*Physics:-Riemann[`~alpha`, `~beta`, `~delta`, `~gamma`]

(2)

Simplify(K);

48*m^2/r^6

(3)

The above concerns your fist post, and to see the Ricci tensor you can always

Ricci[]

Physics:-Ricci[mu, nu] = Matrix(%id = 18446744078267576438)

(4)

Now in your second post, however, you do not use the Physics package as suggested but the old 'tensor' package. The computations with the old 'tensor' are overly complicated and that package is by now obsolete.

 

My suggestion in the previous replies were computing this using Physics, it is really simpler.

 

Copying the metric info from the worksheet attached in your last reply and pasting here:

g00 := (1-2*m/r)^(1+q)

(1-2*m/r)^(1+q)

(5)

g11 := -(1-2*m/r)^(-q)*(1-sin(theta)^2*m^2)^(-q*(2+q))/(1-2*m/r)

-(1-2*m/r)^(-q)*(1-sin(theta)^2*m^2)^(-q*(2+q))/(1-2*m/r)

(6)

g22 := -(1-2*m/r)^(-q)*(1-sin(theta)^2*m^2)^(-q*(2+q))*r^2

-(1-2*m/r)^(-q)*(1-sin(theta)^2*m^2)^(-q*(2+q))*r^2

(7)

g33 := -(1-2*m/r)^(-q)*(1-sin(theta)^2*m^2)^(-q*(2+q))*r^2*sin(theta)^2

-(1-2*m/r)^(-q)*(1-sin(theta)^2*m^2)^(-q*(2+q))*r^2*sin(theta)^2

(8)

Noticing now that your coordinates [t, r, theta, phi] are the same introduced above with the Schwarzschild metric where x[0] = x[4] = t, you can now re-enter the metric in different ways, see Physics:-Setup  .For instance, to avoid disquisitions about the ordering of the coordinates, enter the metric directly as the square of the line element

ds2 := g00*d_(t)^2+g11*d_(r)^2+g22*d_(theta)^2+g33*d_(phi)^2

Setup(metric = ds2)

 

That is all.

 

To compute the Kretschmann scalar for this metric, just Simplify(K) again, and no need to re-enter, re-construct or re-define the Riemann tensor or anything else.

KScalar := Simplify(K)

 

Regarding the Schwarzschild limit when q = 0 mentioned in your worksheet, this is the metric with q <> 0

g_[]

g[mu, nu] = (Matrix(4, 4, {(1, 1) = ((r-2*m)/r)^(-q)*(cos(theta)^2*m^2-m^2+1)^(-q*(2+q))*r/(-r+2*m), (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = -((r-2*m)/r)^(-q)*(cos(theta)^2*m^2-m^2+1)^(-q*(2+q))*r^2, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -((r-2*m)/r)^(-q)*(cos(theta)^2*m^2-m^2+1)^(-q*(2+q))*r^2*sin(theta)^2, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = (r-2*m)*((r-2*m)/r)^q/r}))

(9)

(If you want to see the inversse of the metric, just enter g_[~mu, ~nu, matrix]. )

This is the value of g[mu, nu] at q = 0 :

eval(g[mu, nu] = (Matrix(4, 4, {(1, 1) = ((r-2*m)/r)^(-q)*(cos(theta)^2*m^2-m^2+1)^(-q*(2+q))*r/(-r+2*m), (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = -((r-2*m)/r)^(-q)*(cos(theta)^2*m^2-m^2+1)^(-q*(2+q))*r^2, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -((r-2*m)/r)^(-q)*(cos(theta)^2*m^2-m^2+1)^(-q*(2+q))*r^2*sin(theta)^2, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = (r-2*m)*((r-2*m)/r)^q/r})), q = 0)

g[mu, nu] = (Matrix(4, 4, {(1, 1) = r/(-r+2*m), (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = -r^2, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -r^2*sin(theta)^2, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = (r-2*m)/r}))

(10)

This is the value of the K scalar when q = 0

simplify(eval(KScalar, q = 0))

48*m^2/r^6

(11)

And this is the KScalar for q <> 0 (a bit of a wallpaper even after simplifying in size)

simplify(KScalar, size)

4*((r-2*m)/r)^(2*q)*((cos(theta)+1)*((2+q)*m-r)^2*(cos(theta)^2*m^2-m^2+1)^4*(-(1/2)*r+m)*(cos(theta)-1)*(2+q)^2*cos(theta)^2*m^2*q^2*r*(cos(theta)^2*m^2-m^2+1)^(2*q^2+4*q-2)+(cos(theta)+1)*(cos(theta)^2*m^2-m^2+1)^3*(1+q)^2*(-(1/2)*r+m)*(cos(theta)-1)*(2+q)^2*cos(theta)^2*m^4*q^2*r*(cos(theta)^2*m^2-m^2+1)^(2*q^2+4*q-1)+(7/4)*(((2+q)^2*(q^2+(16/7)*q+12/7)*m^4-(4/7)*r*(2+q)*(q^2+2*q+2)*(q^2+3*q+3)*m^3+(4/7)*(q^8+8*q^7+22*q^6+(43/2)*q^5+(5/2)*q^4+(19/2)*q^3+26*q^2+10*q+3)*r^2*m^2-(4/7)*(q^6+6*q^5+10*q^4+(1/2)*q^3-(9/2)*q^2+(11/2)*q+1/2)*(2+q)*q*r^3*m+(1/7)*q^2*r^4*(q^4+4*q^3+2*q^2-4*q+2)*(2+q)^2)*m^6*cos(theta)^8-4*((2+q)^2*(q^2+(16/7)*q+12/7)*m^6-(3/7)*(q^4+(16/3)*q^3+(40/3)*q^2+16*q+8)*(2+q)*r*m^5+(((89/14)*q^5+(48/7)*q^3+(16/7)*q^4+(88/7)*q^2+12/7+(16/7)*q^7+(2/7)*q^8+(44/7)*q^6+(40/7)*q)*r^2-(2+q)^2*(q^2+(16/7)*q+12/7))*m^4-(2/7)*((q^6+6*q^5+10*q^4+(3/4)*q^3-(11/4)*q^2+(17/2)*q+1)*q*r^2-q^4-25*q-(13/2)*q^3-(39/2)*q^2-12)*(2+q)*r*m^3+(1/14)*(q^2*(q^4+4*q^3+2*q^2-4*q+3)*(2+q)^2*r^2+4*q^6+18*q^5-5*q^4-123*q^3-190*q^2-88*q-24)*r^2*m^2-(2/7)*(2+q)*(q^4+(7/2)*q^3-(3/2)*q^2-(41/4)*q-5/4)*q*r^3*m+(1/14)*q^2*r^4*(q^2+2*q-4)*(2+q)^2)*m^4*cos(theta)^6+6*((2+q)^2*(q^2+(16/7)*q+12/7)*m^8-(2/7)*r*(2+q)*(q^4+6*q^3+18*q^2+24*q+12)*m^7+(((232/21)*q^2+(16/21)*q^7+(158/21)*q^3+(44/21)*q^6+(18/7)*q^4+(40/7)*q+(7/3)*q^5+(2/21)*q^8+12/7)*r^2-2*(2+q)^2*(q^2+(16/7)*q+12/7))*m^6-(2/21)*((q^6+6*q^5+10*q^4+(3/2)*q^3+(9/2)*q^2+22*q+3)*q*r^2-4*q^4-29*q^3-103*q^2-146*q-72)*(2+q)*r*m^5+((1/42)*q^2*(q^4+4*q^3+2*q^2-4*q+8)*(2+q)^2*r^4+(-24/7+(4/7)*q^5+(4/21)*q^6-(331/21)*q^3-(478/21)*q^2-(248/21)*q-(23/7)*q^4)*r^2+(2+q)^2*(q^2+(16/7)*q+12/7))*m^4-(4/21)*(q*(q^4+3*q^3-9*q^2-(101/4)*q-13/4)*r^2+37*q+(11/2)*q^3+(1/2)*q^4+18+(49/2)*q^2)*(2+q)*r*m^3+(1/21)*r^2*(q^2*(q^2+2*q-10)*(2+q)^2*r^2+3*q^5+51*q^4+189*q^3+262*q^2+128*q+36)*m^2-(1/21)*q*r^3*(2+q)*(q^3+31*q^2+65*q+7)*m+(1/3)*q^2*r^4*(2+q)^2)*m^2*cos(theta)^4-4*((2+q)^2*(q^2+(16/7)*q+12/7)*m^7-(1/7)*r*(2+q)*(q^4+8*q^3+32*q^2+48*q+24)*m^6+((12/7+(3/14)*q^5+(16/7)*q^4+(52/7)*q^3+(72/7)*q^2+(40/7)*q)*r^2-2*(2+q)^2*(q^2+(16/7)*q+12/7))*m^5-(1/14)*((q^4+15*q^3+30*q^2+4*q)*r^2-2*q^4-22*q^3-114*q^2-188*q-96)*(2+q)*r*m^4+((3/14)*q^2*r^4*(2+q)^2+(-(76/7)*q-(43/14)*q^4-(121/7)*q^2-24/7-(3/14)*q^5-(157/14)*q^3)*r^2+(2+q)^2*(q^2+(16/7)*q+12/7))*m^3+(1/14)*((q^4+23*q^3+49*q^2+7*q)*r^2-6*q^3-50*q^2-92*q-48)*(2+q)*r*m^2-(5/14)*(q^2*(2+q)^2*r^2-(3/5)*(1+q)*(q^3+6*q^2+16*q+8))*r^2*m-(3/14)*q*r^3*(2+q)*(1+q))*(m-1)*(m+1)*m*cos(theta)^2+((2+q)^2*(q^2+(16/7)*q+12/7)*m^6-(4/7)*r*(2+q)*(1+q)*(q^2+6*q+6)*m^5+(((10/7)*q^4+(40/7)*q+12/7+(46/7)*q^3+(72/7)*q^2)*r^2-2*(2+q)^2*(q^2+(16/7)*q+12/7))*m^4-(8/7)*(((1/4)*q+q^3+(9/4)*q^2)*r^2-(11/2)*q^2-11*q-(1/2)*q^3-6)*(2+q)*r*m^3+((2/7)*q^2*r^4*(2+q)^2-(2/7)*(1+q)*(q^3+6*q^2+20*q+12)*r^2+(2+q)^2*(q^2+(16/7)*q+12/7))*m^2+(2/7)*r*(2+q)*(1+q)*(q*r^2-8*q-12)*m+(12/7)*r^2*(1+q)^2)*(m-1)^2*(m+1)^2)*(cos(theta)^2*m^2-m^2+1)^(2*q*(2+q)))*m^2/((cos(theta)^2*m^2-m^2+1)^4*(-(1/2)*r+m)^2*r^6)

(12)

For more information please check the help pages for g_, Ricci, Riemann and Setup.

``

 

Download KretschmannScalar.mw


Edgardo S. Cheb-Terrab
Physics, Maplesoft

Hi
Note that in Maple 16 you have the same functionality (and furthermore) within the distributed Physics package.

Edgardo S. Cheb-Terrab
Physics, Maplesoft

Hi
Note also there is now access to the version of Physics under development, available for download at http://www.maplesoft.com/products/maple/features/physicsresearch.aspx, with adjustments and fixes as they are ready, including the one I mentioned in my previous reply for this problem you posted.

Edgardo S. Cheb-Terrab
Physics, Maplesoft

Hi
Note also tthere is now access to the version of Physics under development, available for download at http://www.maplesoft.com/products/maple/features/physicsresearch.aspx, including adjustments and fixes as they are ready.

Edgardo S. Cheb-Terrab
Physics, Maplesoft

Hi Peter137
A fix to the problem you posted is now in place, and new: now there is access to the version of Physics under development, available for download at http://www.maplesoft.com/products/maple/features/physicsresearch.aspx. This version includes post 17.01 material, adjustments, fixes, and novelties as they are ready. In the page you can also present your suggestions, etc

Edgardo S. Cheb-Terrab
Physics, Maplesoft


dsolve's solution is correct:

ode := [diff(u(t),t)-f(t)=0, diff(f(t),t)+t*f(t)=0, diff(f(t),t$2)+(t*sin(a)*sin(b))^2*f(t)=0];

[diff(u(t), t)-f(t) = 0, diff(f(t), t)+t*f(t) = 0, diff(diff(f(t), t), t)+t^2*sin(a)^2*sin(b)^2*f(t) = 0]

(1)

dsolve(ode);

{f(t) = 0, u(t) = _C1}

(2)

This simplification (see ?PDEtools/casesplit) explains dsolve's solution:

PDEtools:-casesplit(ode);

`casesplit/ans`([diff(u(t), t) = 0, f(t) = 0], [])

(3)

You can arrive at the same result by hand. Take the first equation, and you have the solution for f(t) as a function of u(t)

isolate(ode[1], f(t));

f(t) = diff(u(t), t)

(4)

Remove now f(t) from the system

eval(ode, %);

[0 = 0, diff(diff(u(t), t), t)+t*(diff(u(t), t)) = 0, diff(diff(diff(u(t), t), t), t)+t^2*sin(a)^2*sin(b)^2*(diff(u(t), t)) = 0]

(5)

Take the second equation

isolate(%[2], diff(u(t),t,t));

diff(diff(u(t), t), t) = -t*(diff(u(t), t))

(6)

Substitute into the third equation

factor(PDEtools:-dsubs(%, %%[3]));

(diff(u(t), t))*(t^2*sin(a)^2*sin(b)^2+t^2-1) = 0

(7)

So, the solution for u(t) is

dsolve(%);

u(t) = _C1

(8)

Then from (4) the solution for f(t) is 0.

 

Download dsolve_solution_is_c.mw

Edgardo S. Cheb-Terrab
Physics, Maplesoft

First 49 50 51 52 53 54 55 Page 51 of 55