ecterrab

13538 Reputation

24 Badges

18 years, 286 days

MaplePrimes Activity


These are answers submitted by ecterrab

Hi Andriy

In the first reply it is shown how to perform to the end the computation you wanted to perform. In the second one I mentioned the implementation, in the downloadable update of Physics, of the simplification N^k = N, where N = ap . am and ap, am respectively represent criation, annihilation fermionic operators. Besides that, in the update of Physics posted today (dowloadable from the Maple Physics: Research & Development updates page) the simplification am^2 = 0 and ap^2 = 0 is also implemented, completing what has been discussed in this thread. Recalling, in connection with one of your questions, Andriy, these new simplifications are performed not using simplify but using Simplify, the command of the Physics package (see ?Physics,Simplify).

Edgardo S. Cheb-Terrab
Physics, Maplesoft


restart

diff(z(x, y), x, x)+2*(diff(z(x, y), x, y))+diff(z(x, y), y, y) = 0

diff(diff(z(x, y), x), x)+2*(diff(diff(z(x, y), x), y))+diff(diff(z(x, y), y), y) = 0

(1)

The transformation you posted:

itr := {u = x+y, v = x-y, omega(u, v) = x*y-z(x, y)}

{u = x+y, v = x-y, omega(u, v) = x*y-z(x, y)}

(2)

Rewrite your transformation with x, y and z(x, y) on the left-hand sides

tr := solve(itr, {x, y, z(x, y)})

{x = (1/2)*u+(1/2)*v, y = -(1/2)*v+(1/2)*u, z(x, y) = (1/4)*u^2-(1/4)*v^2-omega(u, v)}

(3)

Change variables

PDEtools:-dchange(tr, diff(diff(z(x, y), x), x)+2*(diff(diff(z(x, y), x), y))+diff(diff(z(x, y), y), y) = 0)

-4*(diff(diff(omega(u, v), u), u))+2 = 0

(4)

In the above you already see the answer to your question, so isolating omega[u, u]

isolate(-4*(diff(diff(omega(u, v), u), u))+2 = 0, diff(omega(u, v), u, u))

diff(diff(omega(u, v), u), u) = 1/2

(5)


Download PartialDerivative.mw

Edgardo S. Cheb-Terrab
Physics, Maplesoft

Hi Peterm,


restart; with(Physics)

Setup(noncommutativepre = S, quiet)

[noncommutativeprefix = {S}]

(1)

Set your algebra as an algebra for powers of noncommutative operators, only two rules: one for Commutators, the other one for AntiCommutators, both rules respectively mapping into two functions f and g

Setup(%Commutator(S[a]^n, S[b]^m) = f(S[a], S[b], n, m), %AntiCommutator(S[a]^n, S[b]^m) = g(S[a], S[b], n, m))

[algebrarules = {%AntiCommutator(Physics:-`^`(S[a], n), Physics:-`^`(S[b], m)) = g(S[a], S[b], n, m), %Commutator(Physics:-`^`(S[a], n), Physics:-`^`(S[b], m)) = f(S[a], S[b], n, m)}]

(2)

Check it

AntiCommutator(S[a]^r, S[b]^s)

g(S[a], S[b], r, s)

(3)

Commutator(S[a]^r, S[b]^s)

f(S[a], S[b], r, s)

(4)

That is the most important part. You can now write f and g to accomplish the task in different ways.

 

Here is one way, starting with g

g := proc (A, B, n, m) if [n, m]::'Not(list(posint))' or `minus`({A, B}, {S[x], S[y], S[z]}) <> {} then ('AntiCommutator')(A^n, B^m) elif n::even then if m::even then 2 else AntiCommutator(1, B) end if elif m::even then AntiCommutator(A, 1) elif A = B then 2 else Physics:-`*`(A, B)+Physics:-`*`(B, A) end if end proc:

 

For f you need also to take care of cyclic permutations of [x, y, z] to get the correct sign; express first f  in terms of a procedure s, then concentrate on s

f := proc (A, B, n, m) if [n, m]::'Not(list(posint))' or `minus`({A, B}, {S[x], S[y], S[z]}) <> {} then ('Commutator')(A^n, B^m) elif A = B or n::even or m::even then 0 else I*remove(member, [S[x], S[y], S[z]], [A, B])[1]*s(A, B) end if end proc:

 

This is one possible way to code s

s := proc (A, B) options operator, arrow; if verify([A, B], [S[x], S[y], S[z], S[x]], sublist) then 1 else -1 end if end proc:

 

Check it out

s(S[x], S[y]), s(S[y], S[x])

1, -1

(5)

s(S[x], S[z]), s(S[z], S[x])

-1, 1

(6)

s(S[y], S[z]), s(S[z], S[y])

1, -1

(7)

All is in place. So now:

AntiCommutator(S[x], S[y])

Physics:-`*`(S[x], S[y])+Physics:-`*`(S[y], S[x])

(8)

AntiCommutator(S[x]^2, S[y])

2*S[y]

(9)

AntiCommutator(S[x]^2, S[y]^3)

2*S[y]

(10)

AntiCommutator(S[x]^4, S[y]^4)

2

(11)

Commutator(S[x], S[y])

I*S[z]

(12)

Commutator(S[x]^2, S[y])

0

(13)

Commutator(S[y]^3, S[x]^3)

-I*S[z]

(14)

Commutator(S[x]^4, S[y])

0

(15)

All this is reproducible in Maple 17, but anyway it is a good idea to update your Physics library - you can download it from the Maplesoft "Maple Physics: Research & Development" updates page. I will see next week to have things in place so that you can achieve the same results without having to write these f and g functions - say tackling this problem succesfully in a more natural way, similar to the one you tried.



Download PauliSigmaAlgebraRu.mw

 

Edgardo S. Cheb-Terrab
Physics, Maplesoft


Hi*ml76

I just saw your post ... answering two years after is probably irrelevant, but anyway: things have changed since 2011. Using the Physics library up-to-date (you can download it from the the "Maple Physics: Research & Development" updates page), current version:

Physics:-Version()

"/Users/ecterrab/Maple/lib/Physics.mla", `2013, September 6, 13:30 hours`

(1)

Set some commutation rules to see how it works

restart; with(Physics)

Setup(op = {A, B, C, f, g}, %Commutator(A, B) = 1, %Commutator(A, C) = g, %Commutator(B(x), C(y)) = f(x, y))

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

 

[algebrarules = {%Commutator(A, B) = 1, %Commutator(A, C) = g, %Commutator(B(x), C(y)) = f(x, y)}, quantumoperators = {A, B, C, f, g}]

(2)

The first commutator has a constant in the rhs: in these cases you can now index the operators A and B or give functionality to them as desired and the result will always be the right-hand side of type constant  :

Commutator(A, B)

1

(3)

Commutator(A(x), B(y))

1

(4)

Commutator(A[i](x), B[j](y))

1

(5)

In the second commutator, the right-hand side is not constant, so this one works as in previous releases:

Commutator(A, C)

g

(6)

Commutator(A(x), C(y))

Physics:-Commutator(A(x), C(y))

(7)

Commutator(A[i], C[j])

Physics:-Commutator(A[i], C[j])

(8)

The third commutator has functionality in the left-hand and right-hand sides, the same one, so in this case you have more flexibility

Commutator(B(x), C(y))

f(x, y)

(9)

Commutator(B(s), C(t))

f(s, t)

(10)

You can also transform f into a procedure

f := proc (x, y) options operator, arrow; x^y end proc

proc (x, y) options operator, arrow; Physics:-`^`(x, y) end proc

(11)

Commutator(B(u), C(w))

u^w

(12)

f := proc (x, y) options operator, arrow; 0 end proc

proc (x, y) options operator, arrow; 0 end proc

(13)

Commutator(B(a), C(b))

0

(14)

The Physics update is developed assuming the Maple 17 library, but the last approach, with f(x, y) in the right-hand side of the commutator followed by defining f as a procedure, will also work in previous Maples.


Download MaplePrimesCommutato.mw

Edgardo S. Cheb-Terrab 
Physics, Maplesoft

I suggest you to give a look at the section on Quantum Mechanics, subsection on Angular Momentum, within the help page ?Physics,Examples. The situation there illustrated is similar. From there you can either obtain what you want with few changes, or perhaps refine your question a bit more to help you more concretely.

Edgardo S. Cheb-Terrab
Physics, Maplesoft


Hi

I only saw your post now .. anyway:

restart; with(Physics)


From the formula you posted I assume you are working in a 3D Euclidean space; set things accordingly and define a tensor A

Setup(dimension = [3, `+`], coordinates = X, spacetimeindices = lowercaselatin);

`* Partial match of  'coordinates' against keyword 'coordinatesystems'`

 

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

 

`Default differentiation variables for d_, D_ and dAlembertian are: `*{X = (x1, x2, x3)}

 

`Systems of spacetime Coordinates are: `*{X = (x1, x2, x3)}

 

[coordinatesystems = {X}, dimension = 3, signature = `+`, spacetimeindices = lowercaselatin]

 

`Defined objects with tensor properties`

 

{A, Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, mu, nu], Physics:-SpaceTimeVector[mu](X)}

(1)

Your formula is (using the abbreviation ep_ for LeviCivita)

Physics:-`*`(Physics:-`*`(d_[i](d_[m](A[n](X))), ep_[m, n, j]), ep_[i, j, k])

Physics:-d_[i](Physics:-d_[m](A[n](X), [X]), [X])*Physics:-LeviCivita[j, m, n]*Physics:-LeviCivita[i, j, k]

(2)

To obtain the right-hand side you show, you can use Simplify

Simplify(Physics:-d_[i](Physics:-d_[m](A[n](X), [X]), [X])*Physics:-LeviCivita[j, m, n]*Physics:-LeviCivita[i, j, k])

-Physics:-dAlembertian(A[k](X), [X])+Physics:-d_[k](Physics:-d_[n](A[n](X), [X]), [X])

(3)

In the above, the repeated index n corresponds to your index i, and the dAlembertian corresonds to your d_[i]*d_[i]. To understand how this result is formed you can selectively simplify first the product of LeviCivita tensors using the dot operator, as in

`index/Physics/LeviCivita`, "expected %1 indices for the Levi-Civita symbol in a %1-dimensional spacetime; received %2", 4, 3

Physics:-d_[i](Physics:-d_[m](A[n](X), [X]), [X])*(-Physics:-KroneckerDelta[i, m]*Physics:-KroneckerDelta[k, n]+Physics:-KroneckerDelta[i, n]*Physics:-KroneckerDelta[k, m])

(4)

Simplify(Physics:-d_[i](Physics:-d_[m](A[n](X), [X]), [X])*(-Physics:-KroneckerDelta[i, m]*Physics:-KroneckerDelta[k, n]+Physics:-KroneckerDelta[i, n]*Physics:-KroneckerDelta[k, m]))

-Physics:-dAlembertian(A[k](X), [X])+Physics:-d_[k](Physics:-d_[n](A[n](X), [X]), [X])

(5)

``

Download MaplePrimesLeviCivit.mw

Edgardo S. Cheb-Terrab
Physics, Maplesoft

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

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