ecterrab

14540 Reputation

24 Badges

20 years, 22 days

MaplePrimes Activity


These are answers submitted by ecterrab

Hi

Unfortunately, you use the MTM package and it redefines the 'dsolve' command in a way that it stops working properly. You do that after equation (2) of your worksheet.

So first you need to undo that redefinition (or otherwise do not load MTM). You can undo this wrong redefinition as shown in the reviewed worksheet (linked at the end) entering uprotect(dsolve): dsolve := eval(:-dsolve); protect(dsolve);

After doing that there is still another source of error interruption in your formulation: some letters x1... x10 depend on t, as in x10(t), others do not, as x8, and some appear both with and without dependency at the same time, in the system sys1, that you want to solve - that is not correct syntax. So second you need to revise who depends on t and who doesn't. In the reviewed worksheet linked at the end of this answer I made x9 and x10 depend on t always (as oppposed to some times yes and some other times no), but left x8 untouched for instance. In this way you can run dsolve without error interruptions. But then the output is "NULL" meaning the algorithms installed cannot solve this problem.

In summary: the question you mentioned in this post (what is the source of the error interruption) is resolved (the MTM issue) but you need to revise your formulation anyway in order to tell whether the system is or not solvable by Maple's dsolve.

M_(reviewed).mw

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

Hi

I enhanced one bit the design of the types in Physics so that from now on there are no more error interruptions in type commutative, anticommutative or noncommutative. So, the error interruptiion you noticed doesn't happen anymore. In order to install the new version you need to update your Physics.mla library with the one distributed on the Maplesoft R&D Physics webpage.

Independent of that, there is an issue with the generation of Matlab code: after equation (26), within the double loop seen in your worksheet, the code interrupts as shown below, but this has not to do with Physics and only with the way you are formulating your problem, I think. Note as well that in the worksheet you posted, from Physics, you only use diff and `.`, so you could as well, instead of "with(Physics)", use "with(Physics, diff, `.`)".

transverse_linearization_(reviewed).mw

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

The idea of providing these updates literally around the clock is relatively new. So there are updates only for Maple 17 and ahead till the current Maple 2015. By the way the update covers not only Physics, but also the Differential Equations and Mathematical Functions libraries.

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


Just saw your post today. A symbolic matrix is just an object for which the product is not commutative. Therefore you can use any noncommutative objects to produce the result you want. The Physics package is the tool for computing with noncommutative objects, including differentiation among other things.

For example:

with(Physics):

Set three operators (there is a one to one correspondence between a operator and a matrix; you can, if you want, but you do not need to indicate the dimension of the matrices)

Physics:-Setup(operator = {A, B, C})

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

 

[quantumoperators = {A, B, C}]

(1)

Define now your matricial equation (use exponentiation with the product operator as exponent to represent the Dagger, which in the case of real matrices is the transpose)

C(t) = Typesetting:-delayDotProduct(Physics:-Dagger(A(t)), B(t))

C(t) = Physics:-`*`(Physics:-Dagger(A(t)), B(t))

(2)

Differentiate:

diff(C(t) = Physics:-`*`(Physics:-Dagger(A(t)), B(t)), t)

diff(C(t), t) = Physics:-`*`(diff(Physics:-Dagger(A(t)), t), B(t))+Physics:-`*`(Physics:-Dagger(A(t)), diff(B(t), t))

(3)

And that is essentially the result you wanted, all symbolic, it is 100% doable in Maple.

Another way of doing the same: set a noncommutative prefix:

Physics:-Setup(noncommutativeprefix = Z)

[noncommutativeprefix = {Z}]

(4)

Z__1(t) = Typesetting:-delayDotProduct(Physics:-Dagger(Z__2(t)), Z__3(t))

Z__1(t) = Physics:-`*`(Physics:-Dagger(Z__2(t)), Z__3(t))

(5)

diff(Z__1(t) = Physics:-`*`(Physics:-Dagger(Z__2(t)), Z__3(t)), t)

diff(Z__1(t), t) = Physics:-`*`(diff(Physics:-Dagger(Z__2(t)), t), Z__3(t))+Physics:-`*`(Physics:-Dagger(Z__2(t)), diff(Z__3(t), t))

(6)

Or you can also set many prefixes and use each one to represent one of the matrices of your problem:

Physics:-Setup(noncommutativeprefix = {R, S, T})

[noncommutativeprefix = {R, S, T}]

(7)

R(t) = Typesetting:-delayDotProduct(Physics:-Dagger(S(t)), T(t))

R(t) = Physics:-`*`(Physics:-Dagger(S(t)), T(t))

(8)

diff(R(t) = Physics:-`*`(Physics:-Dagger(S(t)), T(t)), t)

diff(R(t), t) = Physics:-`*`(diff(Physics:-Dagger(S(t)), t), T(t))+Physics:-`*`(Physics:-Dagger(S(t)), diff(T(t), t))

(9)

``


Download DifferentiatingSymbolicMatrices.mw

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

 

with(PDEtools):

declare(u(x))

u(x)*`will now be displayed as`*u

(1)

PDE := diff(u(x), x, x)+A*u(x)^2

diff(diff(u(x), x), x)+A*u(x)^2

(2)

The call you are doing to DeterminingPDE automatically enforces integrability conditions assuming the parameters - in this case A -  are arbitrary, and you get as you show:

DetSys := DeterminingPDE([PDE])

{diff(diff(_xi[x](x, u), x), x) = 0, diff(_xi[x](x, u), u) = 0, _eta[u](x, u) = -2*(diff(_xi[x](x, u), x))*u}

(3)

Leading to the infinitesimals

pdsolve(DetSys)

{_eta[u](x, u) = -2*_C1*u, _xi[x](x, u) = _C1*x+_C2}

(4)

Library:-Specialize_Cn({_eta[u](x, u) = -2*_C1*u, _xi[x](x, u) = _C1*x+_C2})

{_eta[u](x, u) = -2*u, _xi[x](x, u) = x}, {_eta[u](x, u) = 0, _xi[x](x, u) = 1}

(5)

And that is the same you obtain in one go if you directly call Infinitesimals

Infinitesimals(PDE)

[_xi[x](x, u) = 1, _eta[u](x, u) = 0], [_xi[x](x, u) = x, _eta[u](x, u) = -2*u]

(6)

If you want to split into cases according to the values of the parameters involved, call DeterminingPDE with the option to not aplpy integrability conditions so that the system you obtain is not simplified and in this example includes A

DetSys := DeterminingPDE([PDE], integrability = false)

`* Partial match of  'integrability' against keyword 'integrabilityconditions'`

 

{diff(diff(_eta[u](x, u), u), u)-2*(diff(diff(_xi[x](x, u), u), x)), 3*(diff(_xi[x](x, u), u))*A*u^2+2*(diff(diff(_eta[u](x, u), u), x))-(diff(diff(_xi[x](x, u), x), x)), -(diff(_eta[u](x, u), u))*A*u^2+2*(diff(_xi[x](x, u), x))*A*u^2+diff(diff(_eta[u](x, u), x), x)+2*_eta[u](x, u)*A*u, diff(diff(_xi[x](x, u), u), u)}

(7)

Now split this system into cases taking A as a parameter

PDEtools:-casesplit({diff(diff(_eta[u](x, u), u), u)-2*(diff(diff(_xi[x](x, u), u), x)), 3*(diff(_xi[x](x, u), u))*A*u^2+2*(diff(diff(_eta[u](x, u), u), x))-(diff(diff(_xi[x](x, u), x), x)), -(diff(_eta[u](x, u), u))*A*u^2+2*(diff(_xi[x](x, u), x))*A*u^2+diff(diff(_eta[u](x, u), x), x)+2*_eta[u](x, u)*A*u, diff(diff(_xi[x](x, u), u), u)}, parameters = A)

`casesplit/ans`([diff(_eta[u](x, u), x) = 0, diff(_eta[u](x, u), u) = _eta[u](x, u)/u, diff(_xi[x](x, u), x) = -(1/2)*_eta[u](x, u)/u, diff(_xi[x](x, u), u) = 0], []), `casesplit/ans`([diff(diff(diff(_eta[u](x, u), u), u), x) = 0, diff(diff(diff(_eta[u](x, u), u), u), u) = 0, diff(diff(_eta[u](x, u), x), x) = 0, diff(diff(_xi[x](x, u), x), x) = 2*(diff(diff(_eta[u](x, u), u), x)), diff(diff(_eta[u](x, u), u), u) = 2*(diff(diff(_xi[x](x, u), u), x)), diff(diff(_xi[x](x, u), u), u) = 0, A = 0], [])

(8)

Leading to the following infinitesimals

map(pdsolve, [%], parameters = A)

Warning, the names [A] indicated as parameters were not found in the system.

 

[{_eta[u](x, u) = _C1*u, _xi[x](x, u) = -(1/2)*_C1*x+_C2}, {A = 0, _eta[u](x, u) = _C1*u*x+_C2*x+(1/2)*_C3*u^2+_C4*u+_C5, _xi[x](x, u) = _C1*x^2+(1/2)*(_C3*u+2*_C7)*x+_C6*u+_C8}]

(9)

map(Library:-Specialize_Cn, [{_eta[u](x, u) = _C1*u, _xi[x](x, u) = -(1/2)*_C1*x+_C2}, {A = 0, _eta[u](x, u) = _C1*u*x+_C2*x+(1/2)*_C3*u^2+_C4*u+_C5, _xi[x](x, u) = _C1*x^2+(1/2)*(_C3*u+2*_C7)*x+_C6*u+_C8}])

[{_eta[u](x, u) = u, _xi[x](x, u) = -(1/2)*x}, {_eta[u](x, u) = 0, _xi[x](x, u) = 1}, {A = 0, _eta[u](x, u) = u*x, _xi[x](x, u) = x^2}, {A = 0, _eta[u](x, u) = x, _xi[x](x, u) = 0}, {A = 0, _eta[u](x, u) = (1/2)*u^2, _xi[x](x, u) = (1/2)*u*x}, {A = 0, _eta[u](x, u) = u, _xi[x](x, u) = 0}, {A = 0, _eta[u](x, u) = 1, _xi[x](x, u) = 0}, {A = 0, _eta[u](x, u) = 0, _xi[x](x, u) = u}, {A = 0, _eta[u](x, u) = 0, _xi[x](x, u) = x}, {A = 0, _eta[u](x, u) = 0, _xi[x](x, u) = 1}]

(10)

NULL

 

 

Download PDETools_(reviewed).mw

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


To have this working the way you want either use p_ . p_, as you said, or in addition to Vectors also load Physics:-`^` and Physics:-`*`, as shown below, and note please that this only works this way in Maple 2015 after updating the Physics library with the one available in the Maplesoft R&D Physics webpage.

restart;with(Physics,`^`,`*`):with(Physics[Vectors]);
Setup(mathematicalnotation=true);
Setup(geometricdifferentiation=true);

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

 

[mathematicalnotation = true]

 

[geometricdifferentiation = true]

(1)

(* Derivation of the standard Hamiltonian from moderately-first principles.
 *)
Physics:-Version();

"/Users/ecterrab/Maple/lib/Physics2015.mla", `2015, November 8, 4:36 hours`

(2)

# We start with the relativistic Hamiltonian H=T+V = E:

H := sqrt(p_^2*c^2+m^2*c^4);

(Physics:-`^`(p_, 2)*c^2+m^2*c^4)^(1/2)

(3)

Note first that the expansion of a power of a vector is implemented, returning the square of the norm:

expand(H)

(Physics:-Vectors:-Norm(p_)^2*c^2+m^2*c^4)^(1/2)

(4)

Introduce now the components as you did

p_ := p1*_i+p2*_j+p3*_k;

_i*p1+_j*p2+_k*p3

(5)

So now H, or its expanded form, are given by

H;

(Physics:-`^`(_i*p1+_j*p2+_k*p3, 2)*c^2+m^2*c^4)^(1/2)

(6)

expand(H)

(c^4*m^2+c^2*p1^2+c^2*p2^2+c^2*p3^2)^(1/2)

(7)

You can differentiate directly, without having to indicate that c is real:

dH := diff(H, p1)

p1*c^2/(Physics:-`^`(_i*p1+_j*p2+_k*p3, 2)*c^2+m^2*c^4)^(1/2)

(8)

``


Download Derivation_of_H_(reviewed).mw


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


with(Physics:-Vectors);

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

(1)

If all what you want is to use _x, _y and _z with the following meaning:

_x = _i, _y = _j, _z = _k;

_x = _i, _y = _j, _z = _k

(2)

Then you can use macro (check its help page), as in

macro(_x = _i, _y = _j, _z = _k):

So that now you can use _x, _y, _z meaning

_x, _y, _z

_i, _j, _k

(3)

If, however, what you want is for _i, _j, and _k to be displayed as _x, _y and _z, then you can use alias (again give a look at its help page please), as in

alias(_x = _i, _y = _j, _z = _k):

_x, _z, _z

_x, _z, _z

(4)

``


Download UnitVectors.mw

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


with(Physics):

This is a vector (postfix is an underscore - see the ?Physics,Vectors page)

f_(x, y, z)

f_(x, y, z)

(1)

Physics:-Vectors:-`+`(Physics:-Vectors:-Gradient(Physics:-Vectors:-Divergence(f_(x, y, z))), Physics:-`*`(K, Physics:-Vectors:-Laplacian(f_(x, y, z))))

Physics:-diff(Physics:-Vectors:-Divergence(f_(x, y, z)), x)*_i+Physics:-diff(Physics:-Vectors:-Divergence(f_(x, y, z)), y)*_j+Physics:-diff(Physics:-Vectors:-Divergence(f_(x, y, z)), z)*_k+K*(diff(diff(f_(x, y, z), y), y))+K*(diff(diff(f_(x, y, z), z), z))+K*(diff(diff(f_(x, y, z), x), x))

(2)

eval(Physics:-diff(Physics:-Vectors:-Divergence(f_(x, y, z)), x)*_i+Physics:-diff(Physics:-Vectors:-Divergence(f_(x, y, z)), y)*_j+Physics:-diff(Physics:-Vectors:-Divergence(f_(x, y, z)), z)*_k+K*(diff(diff(f_(x, y, z), y), y))+K*(diff(diff(f_(x, y, z), z), z))+K*(diff(diff(f_(x, y, z), x), x)), f_(x, y, z) = Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(f[1](x, y, z), _i), Physics:-`*`(f[2](x, y, z), _j)), Physics:-`*`(f[3](x, y, z), _k)));

(diff(diff(f[1](x, y, z), x), x)+diff(diff(f[2](x, y, z), x), y)+diff(diff(f[3](x, y, z), x), z))*_i+(diff(diff(f[1](x, y, z), x), y)+diff(diff(f[2](x, y, z), y), y)+diff(diff(f[3](x, y, z), y), z))*_j+(diff(diff(f[1](x, y, z), x), z)+diff(diff(f[2](x, y, z), y), z)+diff(diff(f[3](x, y, z), z), z))*_k+K*((diff(diff(f[1](x, y, z), y), y))*_i+(diff(diff(f[2](x, y, z), y), y))*_j+(diff(diff(f[3](x, y, z), y), y))*_k)+K*((diff(diff(f[1](x, y, z), z), z))*_i+(diff(diff(f[2](x, y, z), z), z))*_j+(diff(diff(f[3](x, y, z), z), z))*_k)+K*((diff(diff(f[1](x, y, z), x), x))*_i+(diff(diff(f[2](x, y, z), x), x))*_j+(diff(diff(f[3](x, y, z), x), x))*_k)

(3)

convert((diff(diff(f[1](x, y, z), x), x)+diff(diff(f[2](x, y, z), x), y)+diff(diff(f[3](x, y, z), x), z))*_i+(diff(diff(f[1](x, y, z), x), y)+diff(diff(f[2](x, y, z), y), y)+diff(diff(f[3](x, y, z), y), z))*_j+(diff(diff(f[1](x, y, z), x), z)+diff(diff(f[2](x, y, z), y), z)+diff(diff(f[3](x, y, z), z), z))*_k+K*((diff(diff(f[1](x, y, z), y), y))*_i+(diff(diff(f[2](x, y, z), y), y))*_j+(diff(diff(f[3](x, y, z), y), y))*_k)+K*((diff(diff(f[1](x, y, z), z), z))*_i+(diff(diff(f[2](x, y, z), z), z))*_j+(diff(diff(f[3](x, y, z), z), z))*_k)+K*((diff(diff(f[1](x, y, z), x), x))*_i+(diff(diff(f[2](x, y, z), x), x))*_j+(diff(diff(f[3](x, y, z), x), x))*_k), setofequations)

{diff(diff(f[1](x, y, z), x), x)+diff(diff(f[2](x, y, z), x), y)+diff(diff(f[3](x, y, z), x), z)+K*(diff(diff(f[1](x, y, z), y), y))+K*(diff(diff(f[1](x, y, z), z), z))+K*(diff(diff(f[1](x, y, z), x), x)) = 0, diff(diff(f[1](x, y, z), x), y)+diff(diff(f[2](x, y, z), y), y)+diff(diff(f[3](x, y, z), y), z)+K*(diff(diff(f[2](x, y, z), y), y))+K*(diff(diff(f[2](x, y, z), z), z))+K*(diff(diff(f[2](x, y, z), x), x)) = 0, diff(diff(f[1](x, y, z), x), z)+diff(diff(f[2](x, y, z), y), z)+diff(diff(f[3](x, y, z), z), z)+K*(diff(diff(f[3](x, y, z), y), y))+K*(diff(diff(f[3](x, y, z), z), z))+K*(diff(diff(f[3](x, y, z), x), x)) = 0}

(4)

 

Alternatively, you can also enter (1) directly specifying the function components and using compact display, as in

PDEtools:-declare(f(x, y, z))

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

(5)

f_(x, y, z) := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(f[1](x, y, z), _i), Physics:-`*`(f[2](x, y, z), _j)), Physics:-`*`(f[3](x, y, z), _k))

f[1](x, y, z)*_i+f[2](x, y, z)*_j+f[3](x, y, z)*_k

(6)

Physics:-Vectors:-`+`(Physics:-Vectors:-Gradient(Physics:-Vectors:-Divergence(f_(x, y, z))), Physics:-`*`(K, Physics:-Vectors:-Laplacian(f_(x, y, z))))

_i*(diff(diff(f[1](x, y, z), x), x)+diff(diff(f[2](x, y, z), x), y)+diff(diff(f[3](x, y, z), x), z)+K*(diff(diff(f[1](x, y, z), y), y))+K*(diff(diff(f[1](x, y, z), z), z))+K*(diff(diff(f[1](x, y, z), x), x)))+_j*(diff(diff(f[1](x, y, z), x), y)+diff(diff(f[2](x, y, z), y), y)+diff(diff(f[3](x, y, z), y), z)+K*(diff(diff(f[2](x, y, z), y), y))+K*(diff(diff(f[2](x, y, z), z), z))+K*(diff(diff(f[2](x, y, z), x), x)))+_k*(diff(diff(f[1](x, y, z), x), z)+diff(diff(f[2](x, y, z), y), z)+diff(diff(f[3](x, y, z), z), z)+K*(diff(diff(f[3](x, y, z), y), y))+K*(diff(diff(f[3](x, y, z), z), z))+K*(diff(diff(f[3](x, y, z), x), x)))

(7)

convert(_i*(diff(diff(f[1](x, y, z), x), x)+diff(diff(f[2](x, y, z), x), y)+diff(diff(f[3](x, y, z), x), z)+K*(diff(diff(f[1](x, y, z), y), y))+K*(diff(diff(f[1](x, y, z), z), z))+K*(diff(diff(f[1](x, y, z), x), x)))+_j*(diff(diff(f[1](x, y, z), x), y)+diff(diff(f[2](x, y, z), y), y)+diff(diff(f[3](x, y, z), y), z)+K*(diff(diff(f[2](x, y, z), y), y))+K*(diff(diff(f[2](x, y, z), z), z))+K*(diff(diff(f[2](x, y, z), x), x)))+_k*(diff(diff(f[1](x, y, z), x), z)+diff(diff(f[2](x, y, z), y), z)+diff(diff(f[3](x, y, z), z), z)+K*(diff(diff(f[3](x, y, z), y), y))+K*(diff(diff(f[3](x, y, z), z), z))+K*(diff(diff(f[3](x, y, z), x), x))), setofequations)

{diff(diff(f[1](x, y, z), x), x)+diff(diff(f[2](x, y, z), x), y)+diff(diff(f[3](x, y, z), x), z)+K*(diff(diff(f[1](x, y, z), y), y))+K*(diff(diff(f[1](x, y, z), z), z))+K*(diff(diff(f[1](x, y, z), x), x)) = 0, diff(diff(f[1](x, y, z), x), y)+diff(diff(f[2](x, y, z), y), y)+diff(diff(f[3](x, y, z), y), z)+K*(diff(diff(f[2](x, y, z), y), y))+K*(diff(diff(f[2](x, y, z), z), z))+K*(diff(diff(f[2](x, y, z), x), x)) = 0, diff(diff(f[1](x, y, z), x), z)+diff(diff(f[2](x, y, z), y), z)+diff(diff(f[3](x, y, z), z), z)+K*(diff(diff(f[3](x, y, z), y), y))+K*(diff(diff(f[3](x, y, z), z), z))+K*(diff(diff(f[3](x, y, z), x), x)) = 0}

(8)

``


Download VectorialEquation.mw

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

 

This is your input:

plot3d([r, theta, -3.3203*HeunT(.4995036958*3^(2/3)/([1.1233, 6.8291, 12.2089, 17.4216, 32.9018, 53.2699, 64.3282]*[1])^(4/3), 0, (.3138423830*[1.1233, 6.8291, 12.2089, 17.4216, 32.9018, 53.2699, 64.3282])*[1]*3^(1/3), -.3258398511*3^(2/3)*r)*exp(0.5743565187e-1*r*(2.710463448*r^2+3))+4.9407*exp(-(.1148713037*(1.355231724*r^2+3/2))*r)*HeunT(.4277706929*3^(2/3), 0, .3525391488*3^(1/3), .3258398511*3^(2/3)*r)], r = 0 .. 1, theta = 0 .. 2*Pi, coords = cylindrical)

And as you noticed, it produces nothing.  This is the problem:

([1.1233, 6.8291, 12.2089, 17.4216, 32.9018, 53.2699, 64.3282]*[1])^(4/3)

([1.1233, 6.8291, 12.2089, 17.4216, 32.9018, 53.2699, 64.3282]*[1])^(4/3)

(1)

lprint(%)

([1.1233, 6.8291, 12.2089, 17.4216, 32.9018, 53.2699, 64.3282]*[1])^(4/3)

 

The `*` is incorrect of course, not your fault but an issue with 2D, that you can avoid by going to Preferences -> Interface -> Smart Operators: uncheck it so that an invisible multiplication is not inserted automatically between ][ or )(

 

To make the problem evident, when you find something strange using 2D input, the first thing to do is to place the same input in a 1D input line: just open a prompt below and press F5 to go from 2D to 1D input mode, now paste:

plot3d([r, theta, -3.3203*HeunT(.4995036958*3^(2/3)/([1.1233, 6.8291, 12.2089, 17.4216, 32.9018, 53.2699, 64.3282]*[1])^(4/3), 0, (.3138423830*[1.1233, 6.8291, 12.2089, 17.4216, 32.9018, 53.2699, 64.3282])*[1]*3^(1/3), -.3258398511*3^(2/3)*r)*exp(0.5743565187e-1*r*(2.710463448*r^2+3))+4.9407*exp(-(.1148713037*(1.355231724*r^2+3/2))*r)*HeunT(.4277706929*3^(2/3), 0, .3525391488*3^(1/3), .3258398511*3^(2/3)*r)], r = 0 .. 1, theta = 0 .. 2*Pi, coords = cylindrical);

 

Then look closer and you see the problem aforementioned. OK. Fix it now to obtain correct input, and hence the plot

plot3d([r, theta, -3.3203*HeunT(.4995036958*3^(2/3)/([1.1233, 6.8291, 12.2089, 17.4216, 32.9018, 53.2699, 64.3282][1])^(4/3), 0, (.3138423830*[1.1233, 6.8291, 12.2089, 17.4216, 32.9018, 53.2699, 64.3282])[1]*3^(1/3), -.3258398511*3^(2/3)*r)*exp(0.5743565187e-1*r*(2.710463448*r^2+3))+4.9407*exp(-(.1148713037*(1.355231724*r^2+3/2))*r)*HeunT(.4277706929*3^(2/3), 0, .3525391488*3^(1/3), .3258398511*3^(2/3)*r)], r = 0 .. 1, theta = 0 .. 2*Pi, coords = cylindrical)

 

Finally, answering the question you actually posted: by all means Heun functions can be plotted.

 

Download Heun_(reviewed).mw

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


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

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

 

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

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

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

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

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

First 39 40 41 42 43 44 45 Last Page 41 of 59