ecterrab

7080 Reputation

19 Badges

14 years, 11 days

MaplePrimes Activity


These are answers submitted by ecterrab

The form f[a](b) is useful if you want to manipulate the index within the procedure assigned to f. For example:
> f := proc(x) [procname, [args]] end;

and now execute it you see that procname appears indexed. So for instance

> f := proc(x) if procname::indexed then idx := op(procname) else idx := NULL fi; "idx" = idx end;

upon execute will show you the value of idx, and of course, you can then do whatever you want with idx, within a procedure. In summary, t is a useful feature that allows you to work with indexed functions using the indexation as a variable.

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

Besides PDEtools:-difforder, typically, you find other tools in ?PDEtools,Library. In this particular case (degree) what you are looking for is in DifferentialAlgebra[Tools][LeadingDerivative]

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

 

Understanding that you do have a copy of Maple (otherwise, let me know), all you need to do is to input ?PDEtools to open the help page and, on that page, take a look under "Symmetry solution related". The infinitesimals you show are computed with the command PDEtools:-Infinitesimals. For the infinitesimal generators in operator form use InfinitesimalGenerator. The rest is pretty straightforward, I'd say more like an exercise to familiarize with the package, and can be performed with the commands listed and linked under "Symmetry solution related", plus possibly some of the ones in the PDEtools:-Library (see ?PDEtools:-Library).

If you find this suggestion too dense, I suggest you to try to put on a Maple worksheet your attempt to formulate and resolve the problem, and we take that as the departure point. Or let me know here in case you do not have a copy of Maple.

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

Write your equations on a Maple worksheet, then you can use PDEtools:-casesplit to uncouple the system.

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

Check the corresponding help page (?Physics:-Fundiff). That command computes functional derivatives of a given functional.

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

With small touches to your worksheet you see there is no formal_series solution for this problem:
 

restart

NULL

eq1 := diff(A(r), r, r)+(diff(A(r), r))/r+A(r)/r^2-a*r*A(r)+b*r^2*f*B(r)

diff(diff(A(r), r), r)+(diff(A(r), r))/r+A(r)/r^2-a*r*A(r)+b*r^2*f*B(r)

(1)

eq2 := diff(B(r), r, r)+(diff(B(r), r))/r+B(r)/r^2-c*r*A(r)+d*r^2*B(r)

diff(diff(B(r), r), r)+(diff(B(r), r))/r+B(r)/r^2-c*r*A(r)+d*r^2*B(r)

(2)

NULL

NULL

dsolve({eq1, eq2}, {A(r), B(r)})

{A(r) = DESol({-(-b*c*f*r^7+a*d*r^7-d*r^4+2*a*r^3-17)*_Y(r)/r^4-(-d*r^5-a*r^4-3*r)*(diff(_Y(r), r))/r^4-(-d*r^6+a*r^5-r^2)*(diff(diff(_Y(r), r), r))/r^4-2*(diff(diff(diff(_Y(r), r), r), r))/r+diff(diff(diff(diff(_Y(r), r), r), r), r)}, {_Y(r)}), B(r) = (a*r^3*A(r)-(diff(diff(A(r), r), r))*r^2-(diff(A(r), r))*r-A(r))/(b*f*r^4)}

(3)

Instead of tackling the whole system at once, decouple it first

{eq1, eq2}

{diff(diff(A(r), r), r)+(diff(A(r), r))/r+A(r)/r^2-a*r*A(r)+b*r^2*f*B(r), diff(diff(B(r), r), r)+(diff(B(r), r))/r+B(r)/r^2-c*r*A(r)+d*r^2*B(r)}

(4)

PDEtools:-casesplit({diff(diff(A(r), r), r)+(diff(A(r), r))/r+A(r)/r^2-a*r*A(r)+b*r^2*f*B(r), diff(diff(B(r), r), r)+(diff(B(r), r))/r+B(r)/r^2-c*r*A(r)+d*r^2*B(r)})

`casesplit/ans`([B(r) = (a*r^3*A(r)-(diff(diff(A(r), r), r))*r^2-(diff(A(r), r))*r-A(r))/(b*f*r^4), diff(diff(diff(diff(A(r), r), r), r), r) = (-c*A(r)*b*f*r^7+d*A(r)*a*r^7-d*(diff(diff(A(r), r), r))*r^6-d*(diff(A(r), r))*r^5+(diff(diff(A(r), r), r))*a*r^5-d*A(r)*r^4-(diff(A(r), r))*a*r^4+2*a*r^3*A(r)+2*(diff(diff(diff(A(r), r), r), r))*r^3-(diff(diff(A(r), r), r))*r^2-3*(diff(A(r), r))*r-17*A(r))/r^4], [])

(5)

You see: solve for A(r) a 4th order equation, then plug the solution in the 1st equation above to get B(r)

op([1, -1], `casesplit/ans`([B(r) = (a*r^3*A(r)-(diff(diff(A(r), r), r))*r^2-(diff(A(r), r))*r-A(r))/(b*f*r^4), diff(diff(diff(diff(A(r), r), r), r), r) = (-c*A(r)*b*f*r^7+d*A(r)*a*r^7-d*(diff(diff(A(r), r), r))*r^6-d*(diff(A(r), r))*r^5+(diff(diff(A(r), r), r))*a*r^5-d*A(r)*r^4-(diff(A(r), r))*a*r^4+2*a*r^3*A(r)+2*(diff(diff(diff(A(r), r), r), r))*r^3-(diff(diff(A(r), r), r))*r^2-3*(diff(A(r), r))*r-17*A(r))/r^4], []))

diff(diff(diff(diff(A(r), r), r), r), r) = (-c*A(r)*b*f*r^7+d*A(r)*a*r^7-d*(diff(diff(A(r), r), r))*r^6-d*(diff(A(r), r))*r^5+(diff(diff(A(r), r), r))*a*r^5-d*A(r)*r^4-(diff(A(r), r))*a*r^4+2*a*r^3*A(r)+2*(diff(diff(diff(A(r), r), r), r))*r^3-(diff(diff(A(r), r), r))*r^2-3*(diff(A(r), r))*r-17*A(r))/r^4

(6)

Now you see there is no formal_series solution with polynomial coeffs as you want

dsolve(diff(diff(diff(diff(A(r), r), r), r), r) = (-c*A(r)*b*f*r^7+d*A(r)*a*r^7-d*(diff(diff(A(r), r), r))*r^6-d*(diff(A(r), r))*r^5+(diff(diff(A(r), r), r))*a*r^5-d*A(r)*r^4-(diff(A(r), r))*a*r^4+2*a*r^3*A(r)+2*(diff(diff(diff(A(r), r), r), r))*r^3-(diff(diff(A(r), r), r))*r^2-3*(diff(A(r), r))*r-17*A(r))/r^4, A(r), 'formal_series', 'coeffs' = 'polynomial')


There is actually no formal series solution available even without restricting coeffs

NULL

NULL


 

Download dsolve_(reviewed).mw

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

Hi

I can't read well your expected output in LaTeX as you posted, but, in general, this type of computation is concretely simpler using the Physics package.

with(Physics)

 

Your line element (just type * instead of &t in what you posted)

ds2 := -(1-2*M*mu/r)^(1/mu)*dt^2+(1-2*M*mu/r)^(-1/mu)*dr^2+r^2*(1-2*M*mu/r)^(1-1/mu)*(dtheta^2+sin(theta)^2*dphi^2)

 

Set the coordinates, and the metric using your line element

Setup(metric = ds2, coordinates = {X = [t, r, theta, phi]}, quiet)

 

The Ricci trace is

Ricci[trace]

2*M^2*(mu^2-1)*(-(2*M*mu-r)/r)^(1/mu)/(r^2*(2*M*mu-r)^2)

(1)

 

The Kretchmann scalar is the full contraction of the Riemann tensor, so

Riemann[alpha, beta, mu, nu]^2

Physics:-Riemann[alpha, beta, mu, nu]*Physics:-Riemann[`~alpha`, `~beta`, `~mu`, `~nu`]

(2)

SumOverRepeatedIndices(Physics[Riemann][alpha, beta, mu, nu]*Physics[Riemann][`~alpha`, `~beta`, `~mu`, `~nu`], simplifier = simplify)

3*((M*mu-(1/2)*r)^2*((mu+1)^2*(mu^2+(2/3)*mu+1)*M^2-(8/3)*r*(mu+1)^2*M+(8/3)*r^2)*((-2*M*mu+r)/r)^((-2*mu+2)/mu)+(1/3)*r^2*((-2*M*mu+r)/r)^(2/mu)*((mu+1)*M-r)^2)*M^2/((M*mu-(1/2)*r)^4*r^6)

(3)

 

But all the symbols entering here are positive, so

`assuming`([simplify(3*((M*mu-(1/2)*r)^2*((mu+1)^2*(mu^2+(2/3)*mu+1)*M^2-(8/3)*r*(mu+1)^2*M+(8/3)*r^2)*((-2*M*mu+r)/r)^((-2*mu+2)/mu)+(1/3)*r^2*((-2*M*mu+r)/r)^(2/mu)*((mu+1)*M-r)^2)*M^2/((M*mu-(1/2)*r)^4*r^6))], [positive])

12*r^((-4*mu-2)/mu)*((mu+1)^2*(mu^2+(2/3)*mu+7/3)*M^2-(8/3)*r*(mu+2)*(mu+1)*M+4*r^2)*(-2*M*mu+r)^((-4*mu+2)/mu)*M^2

(4)

 

I leave below the input lines you presented using DifferentialGeometry, so that you can compare how much simpler are the Physics input lines above.

 

restart

M > 

# Obtaining Ricci and Kretchmann;
with(DifferentialGeometry):with(Tensor):

DGsetup([t, r, theta, phi], M);
g := evalDG(-(1-2*M*mu/r)^(1/mu)*dt &t dt+(1-2*M*mu/r)^(-1/mu)*`&t`(dr, dr)+r^2*(1-2*M*mu/r)^(1-1/mu)*(`&t`(dtheta, dtheta)+sin(theta)^2*`&t`(dphi, dphi)));
C := Christoffel(g):

`frame name: M`

 

_DG([["tensor", M, [["cov_bas", "cov_bas"], []]], [[[1, 1], -(-(2*M*mu-r)/r)^(1/mu)], [[2, 2], (-(2*M*mu-r)/r)^(-1/mu)], [[3, 3], r^2*(-(2*M*mu-r)/r)^((mu-1)/mu)], [[4, 4], r^2*(-(2*M*mu-r)/r)^((mu-1)/mu)*sin(theta)^2]]])

(5)

Rie := CurvatureTensor(C):
R := RicciScalar(g,Rie);
h := InverseMetric(g):
kretchmann := ContractIndices(RaiseLowerIndices(g, Rie, [1]), RaiseLowerIndices(h, Rie, [2, 3, 4]), [[1, 1], [2, 2], [3, 3], [4, 4]]);

2*(-(2*M*mu-r)/r)^(1/mu)*M^2*(mu^2-1)/(r^2*(2*M*mu-r)^2)

 

4*(-(2*M*mu-r)/r)^(-2*(mu-1)/mu)*M^2*(M*mu^2+2*M*mu+M-2*r)^2/(r^6*(2*M*mu-r)^2)+4*(-(2*M*mu-r)/r)^(2/mu)*M^2*(M*mu^2+M*mu-r)^2/((2*M*mu-r)^4*r^4)+20*(-(2*M*mu-r)/r)^(2/mu)*(M*mu+M-r)^2*M^2/((2*M*mu-r)^4*r^4)+4*(-(2*M*mu-r)/r)^(-2*(mu-1)/mu)*M^2*(M*mu^2+M*mu-r)^2/(r^6*(2*M*mu-r)^2)+4*(-(2*M*mu-r)/r)^(-2*(mu-1)/mu)*M^2*(M*mu+M-r)^2/(r^6*(2*M*mu-r)^2)

(6)
M > 

# simplification

M > 

simplify(normal(R),symbolic)

2*(-1)^(1/mu)*(2*M*mu-r)^((1-2*mu)/mu)*r^((-1-2*mu)/mu)*M^2*(mu^2-1)

(7)
M > 

simplify(kretchmann,size,symbolic)

4*(-(2*M*mu-r)/r)^(-2*(mu-1)/mu)*M^2*(M*mu^2+2*M*mu+M-2*r)^2/(r^6*(2*M*mu-r)^2)+4*(-(2*M*mu-r)/r)^(2/mu)*M^2*(M*mu^2+M*mu-r)^2/((2*M*mu-r)^4*r^4)+20*(-(2*M*mu-r)/r)^(2/mu)*(M*mu+M-r)^2*M^2/((2*M*mu-r)^4*r^4)+4*(-(2*M*mu-r)/r)^(-2*(mu-1)/mu)*M^2*(M*mu^2+M*mu-r)^2/(r^6*(2*M*mu-r)^2)+4*(-(2*M*mu-r)/r)^(-2*(mu-1)/mu)*M^2*(M*mu+M-r)^2/(r^6*(2*M*mu-r)^2)

(8)
M > 

 


 

Download RicciScalarKretchmann.mw

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

Hi nm
Good catch, the unexpected error interruption had to do with an incomplete detection of singularities in the summand. The issue is fixed and the fix distributed within the Maplesoft Physics Updates v.368 (or higher). Both examples you posted are solved normally.

 

lap := diff(u(r, z, t), `$`(r, 2))+(diff(u(r, z, t), r))/r+diff(u(r, z, t), `$`(z, 2)); bc := u(r, 0, t) = 0, u(r, 1, t) = 0, u(1, z, t) = 0; ic := u(r, z, 0) = f(r, z); `assuming`([pdsolve([diff(u(r, z, t), t) = lap, bc, ic], u(r, z, t))], [t > 0])

diff(diff(u(r, z, t), r), r)+(diff(u(r, z, t), r))/r+diff(diff(u(r, z, t), z), z)

 

u(r, 0, t) = 0, u(r, 1, t) = 0, u(1, z, t) = 0

 

u(r, z, 0) = f(r, z)

 

u(r, z, t) = `casesplit/ans`(Sum(Sum(4*BesselJ(0, lambda[n1]*r)*sin(n*Pi*z)*exp(-t*(Pi^2*n^2+lambda[n1]^2))*(Int(BesselJ(0, lambda[n1]*r)*r*(Int(sin(n*Pi*z)*f(r, z), z = 0 .. 1, AllSolutions)), r = 0 .. 1, AllSolutions))/hypergeom([1/2], [1, 2], -lambda[n1]^2), n = 1 .. infinity), n1 = 1 .. infinity), {And(lambda[n1] = BesselJZeros(0, n1), 0 <= lambda[n1])})

(1)

lap := VectorCalculus:-Laplacian(u(r, z, t), cylindrical[r, theta, z]); bc := u(r, 0, t) = 0, u(r, 1, t) = 0, u(1, z, t) = 0; ic := u(r, z, 0) = f(r, z); pdsolve([diff(u(r, z, t), t) = lap, bc, ic], u(r, z, t), HINT = boundedseries(r = 0))

(diff(u(r, z, t), r)+r*(diff(diff(u(r, z, t), r), r))+r*(diff(diff(u(r, z, t), z), z)))/r

 

u(r, 0, t) = 0, u(r, 1, t) = 0, u(1, z, t) = 0

 

u(r, z, 0) = f(r, z)

 

u(r, z, t) = `casesplit/ans`(Sum(Sum(4*BesselJ(0, lambda[n1]*r)*sin(n*Pi*z)*exp(-t*(Pi^2*n^2+lambda[n1]^2))*(Int(BesselJ(0, lambda[n1]*r)*r*(Int(sin(n*Pi*z)*f(r, z), z = 0 .. 1, AllSolutions)), r = 0 .. 1, AllSolutions))/hypergeom([1/2], [1, 2], -lambda[n1]^2), n = 1 .. infinity), n1 = 1 .. infinity), {And(lambda[n1] = BesselJZeros(0, n1), 0 <= lambda[n1])})

(2)

``


 

Download fixed_in_v368.mw

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

Hi nm,

The unexpected error interruption is fixed; as usual, the fix is distributed for everybody within the Maplesoft Physics Updates v.366 or higher. A couple of other ones you posted here the last month are also fixed there. I note however that for this example, pdsolve returns NULL instead of the solution you get using the workaround in this thread by Thomas Richard. By all means, the problem is easy and adjusting the code should be easy too, but I'm kinda overloaded this time ... until the end of June, so maybe this will need to wait.

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

Hi
Below I answer your questions, intercalating some comments in italics, compacting a few input lines when it seemed to me more clear - for the rest it is just your worksheet reviewed.

restart

with(Physics)

Setup(signature = `-+++`, coordinates = (X = [t, x, sigma, phi]))

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

 

`Systems of spacetime coordinates are:`*{X = (t, x, sigma, phi)}

 

_______________________________________________________

 

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

(1)

Setup(g_ = -(c(t, x, sigma)^2-v[1](t, x, sigma)^2-v[2](t, x, sigma)^2)*dt^2+2*v[1](t, x, sigma)*dt*dx+2*v[2](t, x, sigma)*dt*dsigma+dx^2+dsigma^2+sigma^2*dphi^2)

[metric = {(1, 1) = -c(t, x, sigma)^2+v[1](t, x, sigma)^2+v[2](t, x, sigma)^2, (1, 2) = v[1](t, x, sigma), (1, 3) = v[2](t, x, sigma), (2, 2) = 1, (3, 3) = 1, (4, 4) = sigma^2}]

(2)

g_[]

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

(3)

"CompactDisplay(?)"

` c`(t, x, sigma)*`will now be displayed as`*c

 

` v`(t, x, sigma)*`will now be displayed as`*v

(4)

g_[]

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

(5)

Here is an important correction: when you input the functions c and v, you need to specify their functional dependency, that is not displayed only in the output. I am correcting this following input line with regards to that:

Define(l[`~mu`] = [1/sqrt(c(t, x, sigma)^2-v[2](t, x, sigma)^2), (-v[1](t, x, sigma)+sqrt(c(t, x, sigma)^2-v[2](t, x, sigma)^2))/sqrt(c(t, x, sigma)^2-v[2](t, x, sigma)^2), 0, 0])

`Defined objects with tensor properties`

 

{Physics:-D_[mu], Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], l[`~mu`], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

(6)

Check the definition just entered, also with more mathematical format

l[definition]

l[`~mu`] = [1/(c(t, x, sigma)^2-v[2](t, x, sigma)^2)^(1/2), (-v[1](t, x, sigma)+(c(t, x, sigma)^2-v[2](t, x, sigma)^2)^(1/2))/(c(t, x, sigma)^2-v[2](t, x, sigma)^2)^(1/2), 0, 0]

(7)

Only now, that the dependency is there in the input (6), is that this vector is null:

l[]

l[mu] = Array(%id = 18446744078192522654)

(8)

l[mu]^2

l[mu]*l[`~mu`]

(9)

The following is zero as expected

SumOverRepeatedIndices(l[mu]*l[`~mu`])

0

(10)

Likewise, this other vector definition also requires a touch indicating the dependency of the functions c and v:

Define(n[`~mu`] = [1/sqrt(c(t, x, sigma)^2-v[2](t, x, sigma)^2), (-v[1](t, x, sigma)-sqrt(c(t, x, sigma)^2-v[2](t, x, sigma)^2))/sqrt(c(t, x, sigma)^2-v[2](t, x, sigma)^2), 0, 0])

`Defined objects with tensor properties`

 

{Physics:-D_[mu], Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], l[`~mu`], n[`~mu`], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

(11)

n[definition]

n[`~mu`] = [1/(c(t, x, sigma)^2-v[2](t, x, sigma)^2)^(1/2), (-v[1](t, x, sigma)-(c(t, x, sigma)^2-v[2](t, x, sigma)^2)^(1/2))/(c(t, x, sigma)^2-v[2](t, x, sigma)^2)^(1/2), 0, 0]

(12)

n[mu]^2

n[mu]*n[`~mu`]

(13)

SumOverRepeatedIndices(n[mu]*n[`~mu`])

0

(14)

Note as well that there are the null vectors  that come with the Physics:-Tetrads package (they are defined on the fly when you define the metric)

with(Tetrads)

_______________________________________________________

 

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

 

((`Defined as tetrad tensors `(`see ?Physics,tetrads`)*`, `*`&efr;`[a, mu]*`, `)*eta[a, b]*`, `*gamma[a, b, c]*`, `)*lambda[a, b, c]

 

((`Defined as spacetime tensors representing the NP null vectors of the tetrad formalism `(`see ?Physics,tetrads`)*`: `(`see ?Physics,tetrads`)*`, `*l[mu]*`, `)*n[mu]*`, `*m[mu]*`, `)*conjugate(m[mu])

 

_______________________________________________________

(15)

"l_[~]"

Physics:-Tetrads:-l_[`~mu`] = Array(%id = 18446744078238865518)

(16)

l_[mu]^2

Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-l_[`~mu`]

(17)

SumOverRepeatedIndices(Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-l_[`~mu`])

0

(18)

"n_[~]"

Physics:-Tetrads:-n_[`~mu`] = Array(%id = 18446744078378704038)

(19)

n_[mu]^2

Physics:-Tetrads:-n_[mu]*Physics:-Tetrads:-n_[`~mu`]

(20)

SumOverRepeatedIndices(Physics:-Tetrads:-n_[mu]*Physics:-Tetrads:-n_[`~mu`])

0

(21)

The definition of `l__&mu;` and `n__&mu;` in the Tetrads  package is according to standards in the Newman-Penrose formalism

l_[definition]

Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-l_[`~mu`] = 0, Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-n_[`~mu`] = -1, Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-m_[`~mu`] = 0, Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-mb_[`~mu`] = 0, Physics:-g_[mu, nu] = -Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-n_[nu]-Physics:-Tetrads:-l_[nu]*Physics:-Tetrads:-n_[mu]+Physics:-Tetrads:-m_[mu]*Physics:-Tetrads:-mb_[nu]+Physics:-Tetrads:-m_[nu]*Physics:-Tetrads:-mb_[mu]

(22)

With these null vectors you construct a null tetrad. Of course you can rotate the tetrad system of references to obtain different forms of these four null vectors - also match the ones you defined - for that purpose see TransformTetrad .

 

Regarding your other question, a directional derivative, there are two comments.
 

1) Within the Physics package, there is the LieDerivative  command. As you know, in the case of a tensor "derivand", the Lie derivative is the difference between two directional derivatives, but in the case of a scalar, it is a single directional derivative (this is 1/2 of the answer to your question). So for an arbitrary function Phi(X), the directional derivative in the direction of l[mu] or n[mu]is given by, for example

LieDerivative[l](Phi(X))

l[`~mu`]*Physics:-d_[mu](Phi(X), [X])

(23)

2) A directional derivative command can be defined in straightforwardly, no need of a command for that. For example in the direction of "l[]^(mu)"(your vector, not the one that comes with Tetrads, although of course you can do the same with any vector)

L := proc (T) options operator, arrow; l[`~mu`]*`&dtri;`[mu](T) end proc

Note that you are already using the index mu in the definition of L, so do not use it in the argument you pass to L.

For example, define a tensor A

Define(A[nu])

`Defined as tensors`

 

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

(24)

The directional derivative of A^nu in the direction of l is given by

L(A[`~nu`](X))

l[`~mu`]*Physics:-D_[mu](A[`~nu`](X), [X])

(25)

CompactDisplay(l[`~mu`]*D_[mu](A[`~nu`](X), [X]))

` A`(X)*`will now be displayed as`*A

(26)

Not necessary, but useful to be sure what you are computing: check the contents of (25)

convert(l[`~mu`]*D_[mu](A[`~nu`](X), [X]), d_)

l[`~mu`]*(Physics:-d_[mu](A[`~nu`](X), [X])+Physics:-Christoffel[`~nu`, alpha, mu]*A[`~alpha`](X))

(27)

To compute now the components of this vectorial equation, i.e. the four components of this directional derivative. you can use TensorArray , as in TensorArray(l[`~mu`]*(Physics[d_][mu](A[`~nu`](X), [X])+Physics[Christoffel][`~nu`, alpha, mu]*A[`~alpha`](X))). The same applies to the directional derivative of any other tensor. For example, the directional derivative of l^muwith respect to itself  ( from your post, it seems this is what you wanted to compute? ) is given by

L(l[`~nu`])

l[`~mu`]*Physics:-D_[mu](l[`~nu`], [X])

(28)

convert(l[`~mu`]*D_[mu](l[`~nu`], [X]), d_)

l[`~mu`]*(Physics:-d_[mu](l[`~nu`], [X])+Physics:-Christoffel[`~nu`, alpha, mu]*l[`~alpha`])

(29)

To compute the components, again you can use TensorArray

NULL


 

Download try2_(reviewed).mw


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


Hi umbli,

What you are doing is correct, i.e. you can set the coordinates and the line element the way you posted:

restart

with(Physics)

Note the signature:

Setup(sign)

`* Partial match of  '`*sign*`' against keyword '`*signature*`' `

 

_______________________________________________________

 

[signature = `- - - +`]

(1)

So place t in position 4 in the list of coordinates

Setup(coordinates = {X = [x, sigma, theta, t]}, g_ = a(t, x, sigma)*dt^2+b(t, x, sigma)*dt*dx+dx^2+dsigma^2+sigma^2*dtheta^2)

`Systems of spacetime coordinates are:`*{X = (x, sigma, theta, t)}

 

_______________________________________________________

 

[coordinatesystems = {X}, metric = {(1, 1) = 1, (1, 4) = (1/2)*b(t, x, sigma), (2, 2) = 1, (3, 3) = sigma^2, (4, 4) = a(t, x, sigma)}]

(2)

This is the metric

g_[]

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

(3)

And that automatically indicated the dependence of the a and b functions. Independent of that, it is frequently convenient to avoid redundant display of functionality (you know a and b are function, no need to show their functionality all around ...). For that purpose,

"CompactDisplay(?)"

` a`(t, x, sigma)*`will now be displayed as`*a

 

` b`(t, x, sigma)*`will now be displayed as`*b

(4)

g_[]

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

(5)

The above only hides the functionality, but it is there

show

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

(6)

See CompactDisplay. By the way, one thing I think is of the utmost importance when you use computer algebra: always give a look at the help page of a command before using it (even if it is a very brief look). That avoids frustrations and saves a lot of your time.

 

Regarding the Christoffel symbols, they are computed automatically as soon as you enter the metric, e.g:

Christoffel[nonzero]

Physics:-Christoffel[alpha, mu, nu] = {(1, 2, 4) = (1/4)*(diff(b(t, x, sigma), sigma)), (1, 4, 2) = (1/4)*(diff(b(t, x, sigma), sigma)), (1, 4, 4) = (1/2)*(diff(b(t, x, sigma), t))-(1/2)*(diff(a(t, x, sigma), x)), (2, 1, 4) = -(1/4)*(diff(b(t, x, sigma), sigma)), (2, 3, 3) = -sigma, (2, 4, 1) = -(1/4)*(diff(b(t, x, sigma), sigma)), (2, 4, 4) = -(1/2)*(diff(a(t, x, sigma), sigma)), (3, 2, 3) = sigma, (3, 3, 2) = sigma, (4, 1, 1) = (1/2)*(diff(b(t, x, sigma), x)), (4, 1, 2) = (1/4)*(diff(b(t, x, sigma), sigma)), (4, 1, 4) = (1/2)*(diff(a(t, x, sigma), x)), (4, 2, 1) = (1/4)*(diff(b(t, x, sigma), sigma)), (4, 2, 4) = (1/2)*(diff(a(t, x, sigma), sigma)), (4, 4, 1) = (1/2)*(diff(a(t, x, sigma), x)), (4, 4, 2) = (1/2)*(diff(a(t, x, sigma), sigma)), (4, 4, 4) = (1/2)*(diff(a(t, x, sigma), t))}

(7)

"seq(Christoffel[~j,mu,nu, matrix], ~j = [~0, ~1, ~2, ~3])"

Physics:-Christoffel[`~4`, mu, nu] = Matrix(%id = 18446744078247888406), Physics:-Christoffel[`~1`, mu, nu] = Matrix(%id = 18446744078247884190), Physics:-Christoffel[`~2`, mu, nu] = Matrix(%id = 18446744078247885630), Physics:-Christoffel[`~3`, mu, nu] = Matrix(%id = 18446744078318387558)

(8)

I noticed in another post you made that you prefer to work with a signature with time in position 1. That is all OK:

Redefine(setall, tosignature = `-+++`)

[X], Matrix(%id = 18446744078366166782)

(9)

So now you have

Coordinates()

`Systems of spacetime coordinates are:`*{X = (t, x, sigma, theta)}

 

{X}

(10)

g_[]

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

(11)

Just be careful with using 0 as an index, since it always points to the timelike position,which now is 1

"Christoffel[~0,mu,nu, matrix]"

Physics:-Christoffel[`~1`, mu, nu] = Matrix(%id = 18446744078318380086)

(12)

NULL


 

Download signature_and_line_element.mw

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

Hi,
Instead of assume, you can use Physics:-Assume, which is free of this issue of redefininig the variable each time you place an assumption on it, or you add to the assumptions placed. The syntax / use of Physics:-Assume is exactly the same as that of the older lowercase assume; Physics:-Assume also combines the functionality of assume and additionally into one (check the help page for details).

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

 

These first three work fine, as you say in your post

NULL

NULL

NULL

 

The next was the problematic one. A fix for this issue is available to everybody, provided within the Maplesoft Physics Updates v.358.

Physics:-Version()[2]

`2019, May 21, 15:22 hours, version in the MapleCloud: 358, version installed in this computer: 358.`

(1)

This is the system

sys := [diff(u(t, x), t) = diff(u(t, x), x, x), u(0, x) = 1, u(t, -1) = 0, u(t, 1) = 0]

[diff(u(t, x), t) = diff(diff(u(t, x), x), x), u(0, x) = 1, u(t, -1) = 0, u(t, 1) = 0]

(2)

The solution

sol := pdsolve(sys)

u(t, x) = Sum(-4*exp(-(1/4)*Pi^2*(2*n-1)^2*t)*cos((1/2)*(2*n-1)*Pi*x)*(-1)^n/((2*n-1)*Pi), n = 1 .. infinity)

(3)

pdetest(sol, sys)

[0, Sum(-4*cos((1/2)*(2*n-1)*Pi*x)*(-1)^n/((2*n-1)*Pi), n = 1 .. infinity)-1, 0, 0]

(4)

So we need to test by hand the first condition, u(0, x) = 1. For that: take the solution at t=0, change infinity by - say - 1000 terms

u(0, x) = eval(rhs(u(t, x) = Sum(-4*exp(-(1/4)*Pi^2*(2*n-1)^2*t)*cos((1/2)*(2*n-1)*Pi*x)*(-1)^n/((2*n-1)*Pi), n = 1 .. infinity)), [t = 0, infinity = 1000])

u(0, x) = Sum(-4*cos((1/2)*(2*n-1)*Pi*x)*(-1)^n/((2*n-1)*Pi), n = 1 .. 1000)

(5)

Plot from -1 to 1 (the values of x in the 2nd and 3rd boundary conditions), expected: a constant segment at 1 in the vertical axis:

plot(rhs(u(0, x) = Sum(-4*cos((1/2)*(2*n-1)*Pi*x)*(-1)^n/((2*n-1)*Pi), n = 1 .. 1000)), x = -1 .. 1)

 


 

Download PDE_4th_working_fine.mw

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

 

  

pde := (a*y+b*x+c)*(diff(w(x, y), x))-(b*y+k*x+s)*(diff(w(x, y), y)) = 0

(a*y+b*x+c)*(diff(w(x, y), x))-(b*y+k*x+s)*(diff(w(x, y), y)) = 0

(1)

This pde admits several differential invariants. The pdsolve command finds one and proceeds ahead constructing the solution. It is unclear to me whether you can predict that there exists a differential invariant that will lead to a simpler solution without just doing trial and error, and in doing so running the risk of significantly complicating the process while trying to search for a "simpler" solution. In other words, if another program using the same (characteristic strip) method (as Mathematica does) arrives at a simplier solution it is just by chance.

 

Having said that, here is a human-guided way to make the "simpler" solution be discovered before the other ones: change proc (x) options operator, arrow; 1/x end proc, solve the transformed pde, then chage variables back.

PDEtools:-dchange({x = 1/xi, w(x, y) = W(xi, y)}, pde, [xi, W])

-(a*y+b/xi+c)*(diff(W(xi, y), xi))*xi^2-(b*y+k/xi+s)*(diff(W(xi, y), y)) = 0

(2)

pdsolve(-(a*y+b/xi+c)*(diff(W(xi, y), xi))*xi^2-(b*y+k/xi+s)*(diff(W(xi, y), y)) = 0)

W(xi, y) = _F1(-(1/2)*(a*xi^2*y^2+2*c*xi^2*y+2*b*xi*y+2*s*xi+k)/xi^2)

(3)

PDEtools:-dchange({xi = 1/x, W(xi, y) = w(x, y)}, W(xi, y) = _F1(-(1/2)*(a*xi^2*y^2+2*c*xi^2*y+2*b*xi*y+2*s*xi+k)/xi^2), known = all, [x, w], normal)

w(x, y) = _F1(-(1/2)*a*y^2-b*y*x-(1/2)*k*x^2-c*y-s*x)

(4)

pdetest(w(x, y) = _F1(-(1/2)*a*y^2-b*y*x-(1/2)*k*x^2-c*y-s*x), pde)

0

(5)

``


 

Download simplier_solution.mw

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


Maybe I am missing something ... what you suggest, already exists, it is called  macro

macro(v^2 = a^2+b^2)

v^2

a^2+b^2

(1)

v := 5

5

(2)

v^2

a^2+b^2

(3)

v^3

125

(4)

macro(a*b = 3*x+5*y2)

a*b

3*x+5*y2

(5)

a, b := 1, 2

1, 2

(6)

a*b

3*x+5*y2

(7)


 

Download macro.mw

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

1 2 3 4 5 6 7 Last Page 1 of 34