Advanced Search

ecterrab

14645 Reputation

24 Badges

20 years, 152 days
Contact ecterrab

MaplePrimes Activity


These are answers submitted by ecterrab

Not a bug....

ecterrab 14645
October 12 2015
5 4


Hi vv,

Your system

> 

sys := [diff(y(t), t) = y(t)^2-y(t)^3, y(0) = 1/10]

[diff(y(t), t) = y(t)^2-y(t)^3, y(0) = 1/10]

 (1)
> 

sol := dsolve(sys)

y(t) = RootOf(ln(_Z-1)*_Z-ln(_Z)*_Z-ln(9/10)*_Z-I*Pi*_Z-ln(10)*_Z+_Z*t-10*_Z+1)

 (2)

As mentioned by Preben, this RootOf comes from solve. The expression sent to solve is, mainly,

> 

sol__2 := DEtools[remove_RootOf](y(t) = RootOf(ln(_Z-1)*_Z-ln(_Z)*_Z-ln(9/10)*_Z-I*Pi*_Z-ln(10)*_Z+_Z*t-10*_Z+1))

ln(y(t)-1)*y(t)-ln(y(t))*y(t)-ln(9/10)*y(t)-I*Pi*y(t)-ln(10)*y(t)+y(t)*t-10*y(t)+1 = 0

 (3)

This expression satisfies the initial conditions, also mentioned by Preben

> 

eval(sol__2, [t = 0, y(t) = 1/10])

0 = 0

 (4)

Now, your argument to say that dsolve's solution is wrong is that

> 

eval(sol, [t = 0])

y(0) = RootOf(-ln(_Z-1)*_Z+ln(_Z)*_Z+ln(9/10)*_Z+I*Pi*_Z+ln(10)*_Z-1+10*_Z)

 (5)
> 

evalf(y(0) = RootOf(-ln(_Z-1)*_Z+ln(_Z)*_Z+ln(9/10)*_Z+I*Pi*_Z+ln(10)*_Z-1+10*_Z))

y(0) = 0.7601399434e-1-0.4409375242e-1*I

 (6)

and in the above, you argue, if the solution were correct, the right-hand side should be equal to 1/10. It sounds reasonable.

 

My point here, however, is that I do not trust the numerical evaluation of any RootOf unless its argument is polynomial in _Z. Otherwise (as it happens in your example), there may be branch cuts in the expression entering as argument to RootOf, or multiple roots for which we have no closed formulas, and so the numerical algorithm that evaluates RootOf can end up computing "over the branch cut" and choosing one of the two values without you being aware of this arbitrariness of the numerical result, or even worse: start computing a root through one path, to suddenly jump to a different path sufficiently close but belonging to another root, leading to unexpected results without you noticing (typically resuting in a jig-saw plot for instance). In part for these reasons is that dsolve has the option 'implicit' also for ODE-IVP problems, which is rather unusual in computer algebra systems.

So le's take your example, and consider the expression (3) sent to solve (for which solve returned a RootOf) and evaluat it at t = 0, but without telling the value of y(0)

> 

eval(sol__2, [t = 0])

ln(y(0)-1)*y(0)-ln(y(0))*y(0)-ln(9/10)*y(0)-I*Pi*y(0)-ln(10)*y(0)+1-10*y(0) = 0

 (7)

Evaluate this numerically, then solve for y(0) .. requesting*all*solutions

> 

evalf(ln(y(0)-1)*y(0)-ln(y(0))*y(0)-ln(9/10)*y(0)-I*Pi*y(0)-ln(10)*y(0)+1-10*y(0) = 0)

ln(y(0)-1.)*y(0)-1.*ln(y(0))*y(0)-12.19722458*y(0)-(3.141592654*I)*y(0)+1. = 0.

 (8)
> 

solve(ln(y(0)-1.)*y(0)-1.*ln(y(0))*y(0)-12.19722458*y(0)-(3.141592654*I)*y(0)+1. = 0., y(0), allsolutions)

0.7601399433e-1-0.4409375241e-1*I, 0.9999999998e-1-0.3691860874e-11*I

 (9)

And there you are, with two solutions, not just one. The first solution is, actually, your argument to say that dsolve's output is wrong. The second solution, however is, actually, "y(0)=1/10,"which, for the same reason one could take as an indication that dsolve's output is correct. Putting all together, the only thing you can say is that, symbolically (ie exact solutions) , the expression (5) involving a RootOf is correct, in that, as shown above, if you remove the RootOf you see that the solution (3) verifies OK, exactly, as shown in (4), and for these problems you cannot verify correctness by numerically evaluating a RootOf whose argument is not polynomial in _Z.

To see that there are branch cuts in this expression you only need to check whether there is more than one value of y(t) for a given t (ie the function is multivalued). By eye, (3) is linear in t so in this case it is easy:

> 

isolate(subs(y(t) = y, ln(y(t)-1)*y(t)-ln(y(t))*y(t)-ln(9/10)*y(t)-I*Pi*y(t)-ln(10)*y(t)+y(t)*t-10*y(t)+1 = 0), t)

t = (-ln(y-1)*y+ln(y)*y+ln(9/10)*y+I*Pi*y+ln(10)*y+10*y-1)/y

 (10)

Compute now the branch cuts (if any) of the right-hand side

> 

FunctionAdvisor(branch_cuts, rhs(t = (-ln(y-1)*y+ln(y)*y+ln(9/10)*y+I*Pi*y+ln(10)*y+10*y-1)/y))

[(-ln(y-1)*y+ln(y)*y+ln(9/10)*y+I*Pi*y+ln(10)*y+10*y-1)/y, And(0 <= y, y < 1)]

 (11)

Ie there is a branch cut for "0<=y<1, "so the function is multivalued, you cannot rely on the numerical evaluation of RootOf. In these cases the only verification that I think is meaningful is the exact symbolic one, shown above, based on removing the RootOf from the solution.

ADDED AFTER POSTING THIS ANSWER: You can get the solution with LambertW using the appropriate optional argument, not obvious at all but anyway: skipimplicit

> 

dsolve(sys, skipimplicit)

y(t) = 1/(LambertW(9*exp(9)*exp(1)*exp(-t-1))+1)

 (12)
> 

odetest(y(t) = 1/(LambertW(9*exp(9)*exp(1)*exp(-t-1))+1), sys)

[0, 0]

 (13)


Download dsolve_not_a_bug.mw

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

simplify(convert(eq2, exp)) = 0...

ecterrab 14645
October 12 2015
2 1

Hi
Instead of simplify(eq2), try simplify(convert(eq2, exp)) and you receive 0. So, no, pdetest is not returning a wrong result.

Have in mind, in general, that pdetest, as well as odetest, have all sorts of simplification manipulations, some of them very sophisticated, towards zero recognition, that go far beyond a simple call to a simplification command, and that are not easy to imagine immediately. These manipulations were collected analyzing different situations reported by users along the (19) years since pdetest entered the Maple library.

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

Nothing wrong, 'I' is expected in what y...

ecterrab 14645
October 11 2015
2 0

 

Michael

Your definition of RicciT (by the way, why don't you use the already defined Ricci tensor directly with tetrad indices?) involves d_[c] where c is a tetrad index. Then, when you invoke RicciT[1,2], that requires computing the tetrad component d_[1] which, in order to be mapped into derivatives with respect to the spacetime (global system of references) variable, it needs to be multiplied by the tetrad. It so happens that for the tetrad metric you set the corresponding tetrad has square roots and the imaginary unit, and so d_[1] where 1 is a tetrad index will appear with the imaginary unit and square roots around.

 

Here is how it happens

> 

restart; with(Physics); with(Tetrads)

`Setting lowercaselatin letters to represent tetrad indices `

  

0, "%1 is not a command in the %2 package", Tetrads, Physics

  

0, "%1 is not a command in the %2 package", Tetrads, Physics

 (1)

Please note that when you load Tetrads, tetrad indices are automatically set to lowercaselatin, so you do not need to set them again to lowercaselatin (you do that in your Setup command).

 

Copying now from your worksheet, the tetrad metric you want to use us

> 

Matrix([[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]], shape = symmetric)

Matrix([[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])

 (2)

NOTE: you can have this for tetrad metric in place by just telling Physics that you want to work with a null (instead of orthonormal) tetrad and that the signature is +--- (instead of the default in Physics that is ---+). So instead of using the Matrix as you did, here I use these Physics commands directly that result in the same

> 

Physics:-Setup(coordinates = (X = [t, x, y, z]), signature = `+---`, tetradmetric = null)

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

  

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

  

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

  

[coordinatesystems = {X}, signature = `+ - - -`, tetradmetric = {(1, 2) = 1, (3, 4) = -1}]

 (3)
> 

eta_[]

eta[a, b] = (Matrix(4, 4, {(1, 1) = 0, (1, 2) = 1, (1, 3) = 0, (1, 4) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 3) = 0, (3, 4) = -1, (4, 4) = 0}, storage = triangular[upper], shape = [symmetric]))

 (4)

Before proceeding, check the tetrad

> 

e_[]

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

 (5)

You see it contains sqrt(2) and the imaginary unit. It is all correct, according to the definition

> 

e_[definition]

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

 (6)

Check if the definition holds:

> 

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

Matrix(%id = 18446744078378088558)

 (7)

So all is correct. Next you define two tensors, GammaT and RicciT (I imagine that you are using 'T' to mean "in the tetrad local system of references", but for that it suffices to use the Physics Ricci tensor and Christoffel directly with tetrad indices, no need to define yet another tensor). Here I do things in both ways.

 

First your way, I am just copying from your worksheet:

> 

Physics:-Define(GammaT[a, b, c], RicciT[b, c] = -%d_[c](GammaT[`~a`, b, a])+%d_[a](GammaT[`~a`, b, c])-Physics:-`*`(GammaT[`~m`, b, a], GammaT[`~a`, m, c])+Physics:-`*`(GammaT[`~m`, b, c], GammaT[`~a`, m, a])-Physics:-`*`(GammaT[`~a`, b, m], GammaT[`~m`, c, a]-GammaT[`~m`, a, c]))

`Defined objects with tensor properties`

  

{Physics:-Dgamma[mu], GammaT[a, b, c], Physics:-Psigma[mu], RicciT[b, c], Physics:-d_[mu], Physics:-Tetrads:-eta_[a, b], Physics:-g_[mu, nu], Physics:-Tetrads:-l_[mu], Physics:-Tetrads:-m_[mu], Physics:-Tetrads:-mb_[mu], Physics:-Tetrads:-n_[mu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

 (8)

Check your definition of RicciT

> 

RicciT[definition]

RicciT[b, c] = -%d_[c](GammaT[`~a`, b, a])+%d_[a](GammaT[`~a`, b, c])-GammaT[`~m`, b, a]*GammaT[`~a`, m, c]+GammaT[`~m`, b, c]*GammaT[`~a`, m, a]-GammaT[`~a`, b, m]*(GammaT[`~m`, c, a]-GammaT[`~m`, a, c])

 (9)

You see it involves the non-covariant derivative with tetrad indices, but if you want this to become derivatives with respect to the spacetime coordinates (from the global system of references) you will need to multiply by the tetrad. For example, define a tensor in the tetrad system of references, as you did (for that, it suffices to define it with a tetrad index)

> 

Physics:-Define(A[a])

`Defined objects with tensor properties`

  

{A[a], Physics:-Dgamma[mu], GammaT[a, b, c], Physics:-Psigma[mu], RicciT[b, c], Physics:-d_[mu], Physics:-Tetrads:-eta_[a, b], Physics:-g_[mu, nu], Physics:-Tetrads:-l_[mu], Physics:-Tetrads:-m_[mu], Physics:-Tetrads:-mb_[mu], Physics:-Tetrads:-n_[mu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

 (10)

Take now the non-covariant derivative using tetrad indices and indicate a dependency of A[a] with respect to the spacetime variables

> 

PDEtools:-declare(A(X))

A(t, x, y, z)*`will now be displayed as`*A

 (11)
> 

Physics:-d_[b](A[a](X))

Physics:-d_[b](A[a](X), [X])

 (12)

Rewrite now the type of indices of this expression so that it can be computed in terms of derivatives with respect to the spacetime variables

> 

Physics:-d_[b](A[a](X), [X]) = Library:-RewriteTypeOfIndices(Physics:-d_[b](A[a](X), [X]))

Physics:-d_[b](A[a](X), [X]) = Physics:-Tetrads:-e_[b, `~mu`]*Physics:-d_[mu](A[a](X), [X])

 (13)

Recalling now from (5) that the tetrad "`&efr;`[b]^(mu)" involves the imaginary unit and sqrt(2), when you compute the components of this 2x2 matrix behind this tensorial expression (12) you can expect the imaginary unit and sqrt(2) around. To compute this matrix (array), use

> 

TensorArray(Physics[d_][b](A[a](X), [X]))

Matrix(%id = 18446744078408303190)

 (14)

Therefore you will see the imaginary unit in your definition of RicciT, and that is the answer to your question in this thread.

 

Let me show now that you do not need to define RicciT nor GammaT to compute with these tensors in the local system of references (the tetrad system). First note that your spacetime is flat, from where the same is true with respect to your tetrad system, and hence the components of the Ricci tensor in both systems of references are all equal to 0.

 

So to illustrate the use of tetrad indices with Ricci let's first set the curvature to something, and cannot be the standard Schwarzchild metric because in that case too the components of Ricci are 0. So how about the Tolman metric

> 

g_[tol]

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

  

`Parameters: `[R(t, r), E(r)]

  

g[mu, nu] = (Matrix(4, 4, {(1, 1) = -Typesetting:-mcomplete(Typesetting:-msub(Typesetting:-mi("R"), Typesetting:-mi("r")), Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mo("&PartialD;"), Typesetting:-mrow(Typesetting:-mo("&PartialD;"), Typesetting:-mi("r"))), Typesetting:-mspace(width = "0.4em"), Typesetting:-mi("R"), Typesetting:-mo("&ApplyFunction;"), Typesetting:-mfenced(Typesetting:-mi("t"), Typesetting:-mi("r"))))^2/(1+2*E(r)), (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 2) = -R(t, r)^2, (2, 3) = 0, (2, 4) = 0, (3, 3) = -R(t, r)^2*sin(theta)^2, (3, 4) = 0, (4, 4) = 1}, storage = triangular[upper], shape = [symmetric]))

 (15)
> 

PDEtools:-declare(Physics:-g_[mu, nu] = Matrix(%id = 18446744078543569374))

E(r)*`will now be displayed as`*E

  

R(t, r)*`will now be displayed as`*R

 (16)

So here are the components of the Ricci tensor in the global (spacetime) system of references

> 

"Ricci[mu,nu,matrix) "# Here you can use Ricci[] as a shortcut too

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

 (17)

And these are the components of the Ricci tensor in the local (tetrad) system of references

> 

Ricci[a, b, matrix]

Physics:-Ricci[a, b] = Matrix(%id = 18446744078587469206)

 (18)

To understand how one is obtained form the other one you can use the same Library routine used some lines above

> 

Ricci[b, c]:

Physics:-Ricci[b, c] = Physics:-Tetrads:-e_[b, `~mu`]*Physics:-Tetrads:-e_[c, `~nu`]*Physics:-Ricci[mu, nu]

 (19)

This gives you an array of equations, where all the lhs (the components of R[a, b]) are equal to the rhs (the components of "`&efr;`[a]^(mu) `&efr;`[b]^(nu) R[mu,nu]")

> 

TensorArray(Physics[Ricci][b, c] = Physics:-Tetrads:-e_[b, `~mu`]*Physics:-Tetrads:-e_[c, `~nu`]*Physics[Ricci][mu, nu], simplifier = normal)

Matrix(%id = 18446744078558356654)

 (20)

Get the equations:

> 

map(rhs, ArrayElems(Matrix(%id = 18446744078558356654)))

{0 = 0, ((diff(R(t, r), t))*(diff(diff(R(t, r), r), t))-(diff(diff(R(t, r), t), t))*(diff(R(t, r), r))-(diff(E(r), r)))/(R(t, r)*(diff(R(t, r), r))) = ((diff(R(t, r), t))*(diff(diff(R(t, r), r), t))-(diff(diff(R(t, r), t), t))*(diff(R(t, r), r))-(diff(E(r), r)))/(R(t, r)*(diff(R(t, r), r))), (R(t, r)*(diff(diff(R(t, r), t), t))*(diff(R(t, r), r))+(diff(R(t, r), t))^2*(diff(R(t, r), r))+(diff(diff(R(t, r), r), t))*R(t, r)*(diff(R(t, r), t))-(diff(E(r), r))*R(t, r)-2*(diff(R(t, r), r))*E(r))/(R(t, r)^2*(diff(R(t, r), r))) = (R(t, r)*(diff(diff(R(t, r), t), t))*(diff(R(t, r), r))+(diff(R(t, r), t))^2*(diff(R(t, r), r))+(diff(diff(R(t, r), r), t))*R(t, r)*(diff(R(t, r), t))-(diff(E(r), r))*R(t, r)-2*(diff(R(t, r), r))*E(r))/(R(t, r)^2*(diff(R(t, r), r))), -((diff(diff(diff(R(t, r), r), t), t))*R(t, r)+(diff(R(t, r), t))*(diff(diff(R(t, r), r), t))+(diff(diff(R(t, r), t), t))*(diff(R(t, r), r))-(diff(E(r), r)))/(R(t, r)*(diff(R(t, r), r))) = -((diff(diff(diff(R(t, r), r), t), t))*R(t, r)+(diff(R(t, r), t))*(diff(diff(R(t, r), r), t))+(diff(diff(R(t, r), t), t))*(diff(R(t, r), r))-(diff(E(r), r)))/(R(t, r)*(diff(R(t, r), r)))}

 (21)

Check whether they are all just identities of the form A = A

> 

map(evalb, {0 = 0, ((diff(R(t, r), t))*(diff(diff(R(t, r), r), t))-(diff(diff(R(t, r), t), t))*(diff(R(t, r), r))-(diff(E(r), r)))/(R(t, r)*(diff(R(t, r), r))) = ((diff(R(t, r), t))*(diff(diff(R(t, r), r), t))-(diff(diff(R(t, r), t), t))*(diff(R(t, r), r))-(diff(E(r), r)))/(R(t, r)*(diff(R(t, r), r))), (R(t, r)*(diff(diff(R(t, r), t), t))*(diff(R(t, r), r))+(diff(R(t, r), t))^2*(diff(R(t, r), r))+(diff(diff(R(t, r), r), t))*R(t, r)*(diff(R(t, r), t))-(diff(E(r), r))*R(t, r)-2*(diff(R(t, r), r))*E(r))/(R(t, r)^2*(diff(R(t, r), r))) = (R(t, r)*(diff(diff(R(t, r), t), t))*(diff(R(t, r), r))+(diff(R(t, r), t))^2*(diff(R(t, r), r))+(diff(diff(R(t, r), r), t))*R(t, r)*(diff(R(t, r), t))-(diff(E(r), r))*R(t, r)-2*(diff(R(t, r), r))*E(r))/(R(t, r)^2*(diff(R(t, r), r))), -((diff(diff(diff(R(t, r), r), t), t))*R(t, r)+(diff(R(t, r), t))*(diff(diff(R(t, r), r), t))+(diff(diff(R(t, r), t), t))*(diff(R(t, r), r))-(diff(E(r), r)))/(R(t, r)*(diff(R(t, r), r))) = -((diff(diff(diff(R(t, r), r), t), t))*R(t, r)+(diff(R(t, r), t))*(diff(diff(R(t, r), r), t))+(diff(diff(R(t, r), t), t))*(diff(R(t, r), r))-(diff(E(r), r)))/(R(t, r)*(diff(R(t, r), r)))})

{true}

 (22)

Summarizing: nothing is wrong here, and for the way you defined RicciT you can expect the imaginary unit.

 

One last thing: is Ricci[a,b] equal to your RicciT[a,b]? To verify that take the right-hand side of (19) and express it in terms of Christoffel symbols

> 

lhs(Physics:-Ricci[b, c] = Physics:-Tetrads:-e_[b, `~mu`]*Physics:-Tetrads:-e_[c, `~nu`]*Physics:-Ricci[mu, nu]) = convert(rhs(Physics:-Ricci[b, c] = Physics:-Tetrads:-e_[b, `~mu`]*Physics:-Tetrads:-e_[c, `~nu`]*Physics:-Ricci[mu, nu]), Christoffel)

Physics:-Ricci[b, c] = Physics:-Tetrads:-e_[b, `~mu`]*Physics:-Tetrads:-e_[c, `~nu`]*(Physics:-d_[alpha](Physics:-Christoffel[`~alpha`, mu, nu], [X])-Physics:-d_[nu](Physics:-Christoffel[`~alpha`, alpha, mu], [X])+Physics:-Christoffel[`~beta`, mu, nu]*Physics:-Christoffel[`~alpha`, alpha, beta]-Physics:-Christoffel[`~beta`, alpha, mu]*Physics:-Christoffel[`~alpha`, beta, nu])

 (23)

Simplify now the repeated indices

> 

Physics:-Simplify(Physics:-Ricci[b, c] = Physics:-Tetrads:-e_[b, `~mu`]*Physics:-Tetrads:-e_[c, `~nu`]*(Physics:-d_[alpha](Physics:-Christoffel[`~alpha`, mu, nu], [X])-Physics:-d_[nu](Physics:-Christoffel[`~alpha`, alpha, mu], [X])+Physics:-Christoffel[`~beta`, mu, nu]*Physics:-Christoffel[`~alpha`, alpha, beta]-Physics:-Christoffel[`~beta`, alpha, mu]*Physics:-Christoffel[`~alpha`, beta, nu]))

Physics:-Ricci[b, c] = Physics:-d_[alpha](Physics:-Christoffel[`~alpha`, mu, nu], [X])*Physics:-Tetrads:-e_[b, `~mu`]*Physics:-Tetrads:-e_[c, `~nu`]-Physics:-Tetrads:-e_[b, `~mu`]*Physics:-d_[c](Physics:-Christoffel[alpha, mu, `~alpha`], [X])+Physics:-Christoffel[`~alpha`, alpha, beta]*Physics:-Christoffel[`~beta`, b, c]-Physics:-Christoffel[`~alpha`, beta, c]*Physics:-Christoffel[`~beta`, alpha, b]

 (24)

Compare with your definition, and I will substitute here GammaT by Christoffel in the definition just to make it more readable for comparison purposes

> 

subs(GammaT = Christoffel, RicciT[definition])

RicciT[b, c] = -%d_[c](Physics:-Christoffel[`~a`, b, a])+%d_[a](Physics:-Christoffel[`~a`, b, c])-Physics:-Christoffel[`~m`, b, a]*Physics:-Christoffel[`~a`, m, c]+Physics:-Christoffel[`~m`, b, c]*Physics:-Christoffel[`~a`, m, a]-Physics:-Christoffel[`~a`, b, m]*(Physics:-Christoffel[`~m`, c, a]-Physics:-Christoffel[`~m`, a, c])

 (25)

Let's compare the first term of Ricci[b,c] and RicciT[b,c]; if these two are equal we would have

> 

op(2, rhs(Physics:-Ricci[b, c] = Physics:-d_[alpha](Physics:-Christoffel[`~alpha`, mu, nu], [X])*Physics:-Tetrads:-e_[b, `~mu`]*Physics:-Tetrads:-e_[c, `~nu`]-Physics:-Tetrads:-e_[b, `~mu`]*Physics:-d_[c](Physics:-Christoffel[alpha, mu, `~alpha`], [X])+Physics:-Christoffel[`~alpha`, alpha, beta]*Physics:-Christoffel[`~beta`, b, c]-Physics:-Christoffel[`~alpha`, beta, c]*Physics:-Christoffel[`~beta`, alpha, b])) = op(1, rhs(RicciT[b, c] = -%d_[c](Physics:-Christoffel[`~a`, b, a])+%d_[a](Physics:-Christoffel[`~a`, b, c])-Physics:-Christoffel[`~m`, b, a]*Physics:-Christoffel[`~a`, m, c]+Physics:-Christoffel[`~m`, b, c]*Physics:-Christoffel[`~a`, m, a]-Physics:-Christoffel[`~a`, b, m]*(Physics:-Christoffel[`~m`, c, a]-Physics:-Christoffel[`~m`, a, c])))

-Physics:-Tetrads:-e_[b, `~mu`]*Physics:-d_[c](Physics:-Christoffel[alpha, mu, `~alpha`], [X]) = -%d_[c](Physics:-Christoffel[`~a`, a, b])

 (26)

But these two are not equal: to be equal you would need to move "`&efr;`[b]^(mu)" within the `&PartialD;`[c] derivative but you cannot do that because "`&efr;`[b]^(mu)" depends on the coordinates, it does not commute with `&PartialD;`[c]. One could think that your GammaT is not Christoffel at tetrad indices, so let's check what Christoffel at tetrad indices is:

> 

Christoffel[a, b, c]; -1; % = Library:-RewriteTypeOfIndices(%)

Physics:-Christoffel[a, b, c] = Physics:-Tetrads:-e_[a, `~mu`]*Physics:-Tetrads:-e_[b, `~nu`]*Physics:-Tetrads:-e_[c, `~alpha`]*Physics:-Christoffel[mu, alpha, nu]

 (27)

Summarizing about these releationships: I suggest you to check carefully whether RicciT and GammaT are the objects you want them to be, and whether you really need to introduce them or perhaps it is simpler to work with the Ricci and Christoffel symbols just using tetrad indices.

> 

``

NOTE (after having posted): it might also be the case that what you want to use instead of RicciT is actually Tetrads:-gamma_, that is, the Ricci rotation coefficients in the (local) tetrad system of references.

Download _No_Problem_with_Defined_Tetrad_Function_(reviewed).mw

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

alias...

ecterrab 14645
October 10 2015
2 0

You can always use an alias, for instance as shown in this image:

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

PDEtools and its symmetry commands...

ecterrab 14645
October 09 2015
1 1

Hi

Just to let you know, that the symmetry routines of PDEtools are Maple native (official) and more thorough by large and accross the board, besides being used (and so tested) by people around the world. Please give a look at the help page ?PDEtools, the section on Symmetry commands. If after that there is still a computation that you think you cannot do it with these PDEtools / Symmetry routines, please post the problem and I will step in to show how you can do it.

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

The Physics and Physics:-Vectors package...

ecterrab 14645
October 07 2015
3 0

Hi

Please give a look at the help page ?Physics,Examples, where in the first section, on Vector Analysis, you will find material like the one you posted using vector notation, and in the one for Relativity you will find material more advanced now using tensor notation (covariant and contravariant, for Galilean or curved spaces or spacetimes).

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

Bug. Fixed....

ecterrab 14645
October 05 2015
2 9

This is a subtle bug; although rare to happen in the DE programs, these contain a formidably large amount of code ... Anyway it is fixed now. The fix is available to everybody as usual, in the Mapelsoft R&D web page for Differential Equations and Mathematical Functions. Note however that this fix only works in Maple 2015 so maybe you want to consider upgrading, there are good reasons for that independent of this fix. If that is not at reach for you please write me an email to physics@maplesoft.com and I will see if I can help perhaps guiding in some way for you to have this fixed in Maple 16.

Oct/6 update: Here is is a maple library, Maple16_DE_update.mla.zip, that you can drop in the lib directory of your Maple 16 installation and fixes the issue, so that now you have:

 

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

more details?...

ecterrab 14645
October 03 2015
1 7

@Dmitry Lyakhov 

There is no command in the Physics package to automatically find the maximal abelian subalgebra of a given Lie algebra. 

I worked on an implementation of derived subalgebras in the context of differential equations time ago, and the related command is DEtools[solve_group]. I didn't find the time to do this more generally in the more modern symmetry routines of PDEtools.

But then there is the DifferentialGeometry package that we introduce not many years ago, with commands for performing these computations, in either semi-automatic way (not "one go", but almost), or step by step, with commands for performing - I'd say - all the seps.

So why don't you post your problem? Giving a look at the details it may be possible to understand why is that the existing commands are not sufficient (are they really not?) or what is that you find non-algorithmic (as far as I recall,  polynomial factorization is fully algorithmic - see for instance The complete analysis of a polynomial factorization algorithm over finite fields).

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

missing context for your question...

ecterrab 14645
October 02 2015
2 6

 

Hi Michael

Your question misses context: the d_ Physics command is an indexed differential operator that is applied onto a "derivand", you know. But your question shows nothing of that. So without further details from you, this is the basics about how it works

> 

with(Physics); -1; Setup(coordinates = X, tensors = A)

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

  

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

  

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

  

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

 (1)
> 

Physics:-d_[mu](A[mu](X))

Physics:-d_[mu](A[`~mu`](X), [X])

 (2)

The above cannot be expressed in terms of diff or D because mu has not been given any value and, you know, the differentiation variable depends on the value of mu.

Then, if you perform the summation implicit in the repeated indices, the output already comes in diff format:

> 

Physics:-SumOverRepeatedIndices(Physics:-d_[mu](A[`~mu`](X), [X]))

diff(A[`~1`](X), x1)+diff(A[`~2`](X), x2)+diff(A[`~3`](X), x3)+diff(A[`~4`](X), x4)

 (3)

The same does not happen with the inert form %d_, but then you can convert to diff or D using convert (that you say it has no effect, but actually works fine as expected):

> 

subs(Physics:-d_ = %d_, Physics:-d_[mu](A[`~mu`](X), [X]))

%d_[mu](A[`~mu`](X), [X])

 (4)
> 

Physics:-SumOverRepeatedIndices(%d_[mu](A[`~mu`](X), [X]))

%d_[1](A[`~1`](X), [X])+%d_[2](A[`~2`](X), [X])+%d_[3](A[`~3`](X), [X])+%d_[4](A[`~4`](X), [X])

 (5)
> 

convert(%d_[1](A[`~1`](X), [X])+%d_[2](A[`~2`](X), [X])+%d_[3](A[`~3`](X), [X])+%d_[4](A[`~4`](X), [X]), diff)

diff(A[`~1`](X), x1)+diff(A[`~2`](X), x2)+diff(A[`~3`](X), x3)+diff(A[`~4`](X), x4)

 (6)
> 

convert(%d_[1](A[`~1`](X), [X])+%d_[2](A[`~2`](X), [X])+%d_[3](A[`~3`](X), [X])+%d_[4](A[`~4`](X), [X]), D)

(D[1](A[`~1`]))(X)+(D[2](A[`~2`]))(X)+(D[3](A[`~3`]))(X)+(D[4](A[`~4`]))(X)

 (7)

So what is the context for you question? What is the actual expression you want to convert that you say is not being converted?

 

Download convert_inert_d_to_D.mw

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

The problem and its workaround...

ecterrab 14645
September 30 2015
1 1

I actually like and use 2D Math input, but it has its problems, frequently invisible, and this one you mention looks like that kind of thing.

The simple workaround, that more than that it is a saving of your time all around, is to not rewrite the commands found in this site you mention, but instead use these commands as already written by Maple. I mean:

  1. Open the help page by typing ?Poincare,examples (you see I am entering the question mark "?" to invoke the help page) then pressing Enter.
  2. The page that opens is actually a worksheet with everything written, you can reproduce things, play around, change, etc. and more than nothing: you don't need to spend your time rewriting :) 

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

don't use spacetime contracted indices a...

ecterrab 14645
September 29 2015
3 9

It is not, at all, that a loop could change or corrupt in this way a tetrad definition.

The answer actually is that you shall not use spacetime contracted indices, that enter your definition of RicciT, I am talking of a (see the line that starts with 'Define'), also as a loop dummy variable, because in this way you corrupt your definition (because the repeated spacetime index a ends up having the value of the loop index, that actually refers to a free, not repeated index).

So the correction in this case is to use something else for the loop dummies, for instance i and j, that do not enter as repeated indices in the definition of your tensors or, to avoid having to read things by eye to discard letters, once you choose using lowercase latin to represent spacetime indices, then, in general, use some other kind of letter for for-do-od loops, for instance alpha and beta, or A and B.

 

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

You can use the Physics package...

ecterrab 14645
September 25 2015
2 1

Hi Damon,
This is also the kind of operation you do directly by using the Physics package; in this example using its Vectors subpackage; please give a look at the help page ?Physics,Vectors

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

The Physics package...

ecterrab 14645
September 25 2015
3 0

Tensor algebra using standard tensor indicial notation, including the LeviCivita and KroneckerDelta tensors (your epsilon and delta) is implemented in the Physics package, please give a look at the help page ?Physics.

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

Yes...

ecterrab 14645
September 22 2015
1 0

Yes the output is correct. Recalling, differentiation with respect to a 'function' is not an ambiguous operation but just differentiation with respect to a jet variable of a jet space. Skipping technical details you can think of a jet space as one where independent variables, dependent variables (what we call 'functions') and their partial derivatives are all "independent objects in equal footing". That is, q and dq/dt are independent objects, and so  d/dq (dq/dt) = 0. See the help page for PDEtools:-ToJet. Another page you can give a look for a concrete example of the form you are posting (a bit more general) is ?Physics,diff, in the Examples section, starting at equation number (6).

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

The error message...

ecterrab 14645
September 19 2015
3 1

Hi

The error message tells what the problem is. Change, in your definition of a (first equation), f(x) by f(x(t)) and everything will work as you expect.

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

First 40 41 42 43 44 45 46 Last Page 42 of 60