ecterrab

7374 Reputation

19 Badges

14 years, 77 days

MaplePrimes Activity


These are answers submitted by ecterrab

I may not be understanding your question precisely - what do you mean by "ignoring the Christoffel symbols with coefficients of order Omega(r)^2". Given the metric g[mu, nu], you cannot force the Christoffel symbols be something different than their definition

Christoffel[definition]

Physics:-Christoffel[alpha, mu, nu] = (1/2)*Physics:-d_[nu](Physics:-g_[alpha, mu], [X])+(1/2)*Physics:-d_[mu](Physics:-g_[alpha, nu], [X])-(1/2)*Physics:-d_[alpha](Physics:-g_[mu, nu], [X])

(1)

So assuming you meant to depart from GAMMA[alpha, mu, nu] and create another tensor, say `Ѓ`[alpha, mu, nu] with the same dimensions and components but where all the ones with degree > 1 in Omega(r) are equal to 0, below you see one way of doing that. Using the definition of other tensors that are functions of GAMMA[alpha, mu, nu] you can compute the equivalent ones based on `Ѓ`[alpha, mu, nu].

with(Physics)

[`*`, `.`, Annihilation, AntiCommutator, Antisymmetrize, Assume, Bra, Bracket, Check, Christoffel, Coefficients, Commutator, CompactDisplay, Coordinates, Creation, D_, Dagger, Decompose, Define, Dgamma, Einstein, EnergyMomentum, Expand, ExteriorDerivative, Factor, FeynmanDiagrams, Fundiff, Geodesics, GrassmannParity, Gtaylor, Intc, Inverse, Ket, KillingVectors, KroneckerDelta, LeviCivita, Library, LieBracket, LieDerivative, Normal, NumericalRelativity, Parameters, PerformOnAnticommutativeSystem, Projector, Psigma, Redefine, Ricci, Riemann, Setup, Simplify, SortProducts, SpaceTimeVector, StandardModel, SubstituteTensor, SubstituteTensorIndices, SumOverRepeatedIndices, Symmetrize, TensorArray, Tetrads, ThreePlusOne, ToFieldComponents, ToSuperfields, Trace, TransformCoordinates, Vectors, Weyl, `^`, dAlembertian, d_, diff, g_, gamma_]

(2)

In the worksheet, you posted you forgot to set the coordinates, and from trial an error I can see you also changed the signature to be (+ - - -). I am adding that line here:

 

Setup(signature = `+---`, coordinates = spherical)

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

(3)

This is the metric you indicated

Setup(metric = -exp(2*alpha(r))*%d_(t)^2+exp(2*beta(r))*%d_(r)^2+r^2*%d_(theta)^2+r^2*sin(theta)^2*(%d_(phi)-Omega(r)*%d_(t))^2)

[metric = {(1, 1) = -exp(2*alpha(r))+(-r^2*cos(theta)^2+r^2)*Omega(r)^2, (1, 4) = -r^2*sin(theta)^2*Omega(r), (2, 2) = exp(2*beta(r)), (3, 3) = r^2, (4, 4) = r^2*sin(theta)^2}]

(4)

g_[]

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

(5)

These are the nonzero components

Christoffel[nonzero]

Physics:-Christoffel[alpha, mu, nu] = {(1, 1, 2) = -(diff(alpha(r), r))*exp(2*alpha(r))-r*Omega(r)*(cos(theta)-1)*(cos(theta)+1)*(r*(diff(Omega(r), r))+Omega(r)), (1, 1, 3) = r^2*sin(theta)*Omega(r)^2*cos(theta), (1, 2, 1) = -(diff(alpha(r), r))*exp(2*alpha(r))-r*Omega(r)*(cos(theta)-1)*(cos(theta)+1)*(r*(diff(Omega(r), r))+Omega(r)), (1, 2, 4) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (1, 3, 1) = r^2*sin(theta)*Omega(r)^2*cos(theta), (1, 3, 4) = -r^2*sin(theta)*Omega(r)*cos(theta), (1, 4, 2) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (1, 4, 3) = -r^2*sin(theta)*Omega(r)*cos(theta), (2, 1, 1) = (diff(alpha(r), r))*exp(2*alpha(r))+r*Omega(r)*(cos(theta)-1)*(cos(theta)+1)*(r*(diff(Omega(r), r))+Omega(r)), (2, 1, 4) = (1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (2, 2, 2) = (diff(beta(r), r))*exp(2*beta(r)), (2, 3, 3) = -r, (2, 4, 1) = (1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (2, 4, 4) = -r*sin(theta)^2, (3, 1, 1) = -r^2*sin(theta)*Omega(r)^2*cos(theta), (3, 1, 4) = r^2*sin(theta)*Omega(r)*cos(theta), (3, 2, 3) = r, (3, 3, 2) = r, (3, 4, 1) = r^2*sin(theta)*Omega(r)*cos(theta), (3, 4, 4) = -r^2*sin(theta)*cos(theta), (4, 1, 2) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (4, 1, 3) = -r^2*sin(theta)*Omega(r)*cos(theta), (4, 2, 1) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (4, 2, 4) = r*sin(theta)^2, (4, 3, 1) = -r^2*sin(theta)*Omega(r)*cos(theta), (4, 3, 4) = r^2*sin(theta)*cos(theta), (4, 4, 2) = r*sin(theta)^2, (4, 4, 3) = r^2*sin(theta)*cos(theta)}

(6)

You can tell the system to make equal to zero the components having degree > 1 in Omega(r), and use that Array to define another tensor, `Ѓ`[alpha, mu, nu]. For that purpose, take a closer look at the structure (6): you want to operate on all algebraic expressions (note that equations `=` are not of type algebraic). Also, the degree command will get confused with diff(Omega(r), r) but not with (D(Omega))(r) which is the same derivative but represented using D  notation. You switch from the diff  and D notations using convert . The command is subsindets . Puting all together,

convert(subsindets(convert(Physics:-Christoffel[alpha, mu, nu] = {(1, 1, 2) = -(diff(alpha(r), r))*exp(2*alpha(r))-r*Omega(r)*(cos(theta)-1)*(cos(theta)+1)*(r*(diff(Omega(r), r))+Omega(r)), (1, 1, 3) = r^2*sin(theta)*Omega(r)^2*cos(theta), (1, 2, 1) = -(diff(alpha(r), r))*exp(2*alpha(r))-r*Omega(r)*(cos(theta)-1)*(cos(theta)+1)*(r*(diff(Omega(r), r))+Omega(r)), (1, 2, 4) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (1, 3, 1) = r^2*sin(theta)*Omega(r)^2*cos(theta), (1, 3, 4) = -r^2*sin(theta)*Omega(r)*cos(theta), (1, 4, 2) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (1, 4, 3) = -r^2*sin(theta)*Omega(r)*cos(theta), (2, 1, 1) = (diff(alpha(r), r))*exp(2*alpha(r))+r*Omega(r)*(cos(theta)-1)*(cos(theta)+1)*(r*(diff(Omega(r), r))+Omega(r)), (2, 1, 4) = (1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (2, 2, 2) = (diff(beta(r), r))*exp(2*beta(r)), (2, 3, 3) = -r, (2, 4, 1) = (1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (2, 4, 4) = -r*sin(theta)^2, (3, 1, 1) = -r^2*sin(theta)*Omega(r)^2*cos(theta), (3, 1, 4) = r^2*sin(theta)*Omega(r)*cos(theta), (3, 2, 3) = r, (3, 3, 2) = r, (3, 4, 1) = r^2*sin(theta)*Omega(r)*cos(theta), (3, 4, 4) = -r^2*sin(theta)*cos(theta), (4, 1, 2) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (4, 1, 3) = -r^2*sin(theta)*Omega(r)*cos(theta), (4, 2, 1) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (4, 2, 4) = r*sin(theta)^2, (4, 3, 1) = -r^2*sin(theta)*Omega(r)*cos(theta), (4, 3, 4) = r^2*sin(theta)*cos(theta), (4, 4, 2) = r*sin(theta)^2, (4, 4, 3) = r^2*sin(theta)*cos(theta)}, D), algebraic, proc (u) options operator, arrow; if 1 < degree(u, Omega(r)) then 0 else u end if end proc), diff)

Physics:-Christoffel[alpha, mu, nu] = {(1, 1, 2) = -(diff(alpha(r), r))*exp(2*alpha(r)), (1, 1, 3) = 0, (1, 2, 1) = -(diff(alpha(r), r))*exp(2*alpha(r)), (1, 2, 4) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (1, 3, 1) = 0, (1, 3, 4) = -r^2*sin(theta)*Omega(r)*cos(theta), (1, 4, 2) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (1, 4, 3) = -r^2*sin(theta)*Omega(r)*cos(theta), (2, 1, 1) = (diff(alpha(r), r))*exp(2*alpha(r)), (2, 1, 4) = (1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (2, 2, 2) = (diff(beta(r), r))*exp(2*beta(r)), (2, 3, 3) = -r, (2, 4, 1) = (1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (2, 4, 4) = -r*sin(theta)^2, (3, 1, 1) = 0, (3, 1, 4) = r^2*sin(theta)*Omega(r)*cos(theta), (3, 2, 3) = r, (3, 3, 2) = r, (3, 4, 1) = r^2*sin(theta)*Omega(r)*cos(theta), (3, 4, 4) = -r^2*sin(theta)*cos(theta), (4, 1, 2) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (4, 1, 3) = -r^2*sin(theta)*Omega(r)*cos(theta), (4, 2, 1) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (4, 2, 4) = r*sin(theta)^2, (4, 3, 1) = -r^2*sin(theta)*Omega(r)*cos(theta), (4, 3, 4) = r^2*sin(theta)*cos(theta), (4, 4, 2) = r*sin(theta)^2, (4, 4, 3) = r^2*sin(theta)*cos(theta)}

(7)

In the above you see, for instance, that the component (1,1,3) is now equal to 0. As said, you cannot define GAMMA[alpha, mu, nu] to be something different than its definition, but you can  Define  another tensor with these components.

Define(`&GJcy;`[alpha, mu, nu] = Array(`$`(1 .. 4, 3), rhs(Physics:-Christoffel[alpha, mu, nu] = {(1, 1, 2) = -(diff(alpha(r), r))*exp(2*alpha(r)), (1, 1, 3) = 0, (1, 2, 1) = -(diff(alpha(r), r))*exp(2*alpha(r)), (1, 2, 4) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (1, 3, 1) = 0, (1, 3, 4) = -r^2*sin(theta)*Omega(r)*cos(theta), (1, 4, 2) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (1, 4, 3) = -r^2*sin(theta)*Omega(r)*cos(theta), (2, 1, 1) = (diff(alpha(r), r))*exp(2*alpha(r)), (2, 1, 4) = (1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (2, 2, 2) = (diff(beta(r), r))*exp(2*beta(r)), (2, 3, 3) = -r, (2, 4, 1) = (1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (2, 4, 4) = -r*sin(theta)^2, (3, 1, 1) = 0, (3, 1, 4) = r^2*sin(theta)*Omega(r)*cos(theta), (3, 2, 3) = r, (3, 3, 2) = r, (3, 4, 1) = r^2*sin(theta)*Omega(r)*cos(theta), (3, 4, 4) = -r^2*sin(theta)*cos(theta), (4, 1, 2) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (4, 1, 3) = -r^2*sin(theta)*Omega(r)*cos(theta), (4, 2, 1) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (4, 2, 4) = r*sin(theta)^2, (4, 3, 1) = -r^2*sin(theta)*Omega(r)*cos(theta), (4, 3, 4) = r^2*sin(theta)*cos(theta), (4, 4, 2) = r*sin(theta)^2, (4, 4, 3) = r^2*sin(theta)*cos(theta)})))

{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], Ѓ[alpha, mu, nu], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

(8)

Set a macro for `&GJcy;` to avoid having to use the palette all the time

macro(G = `&GJcy;`)

Check it out

G[1, 1, 3]

0

(9)

G[nonzero]

Ѓ[mu, nu, alpha] = {(1, 1, 2) = -(diff(alpha(r), r))*exp(2*alpha(r)), (1, 2, 1) = -(diff(alpha(r), r))*exp(2*alpha(r)), (1, 2, 4) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (1, 3, 4) = -r^2*sin(theta)*Omega(r)*cos(theta), (1, 4, 2) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (1, 4, 3) = -r^2*sin(theta)*Omega(r)*cos(theta), (2, 1, 1) = (diff(alpha(r), r))*exp(2*alpha(r)), (2, 1, 4) = (1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (2, 2, 2) = (diff(beta(r), r))*exp(2*beta(r)), (2, 3, 3) = -r, (2, 4, 1) = (1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (2, 4, 4) = -r*sin(theta)^2, (3, 1, 4) = r^2*sin(theta)*Omega(r)*cos(theta), (3, 2, 3) = r, (3, 3, 2) = r, (3, 4, 1) = r^2*sin(theta)*Omega(r)*cos(theta), (3, 4, 4) = -r^2*sin(theta)*cos(theta), (4, 1, 2) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (4, 1, 3) = -r^2*sin(theta)*Omega(r)*cos(theta), (4, 2, 1) = -(1/2)*r*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r)), (4, 2, 4) = r*sin(theta)^2, (4, 3, 1) = -r^2*sin(theta)*Omega(r)*cos(theta), (4, 3, 4) = r^2*sin(theta)*cos(theta), (4, 4, 2) = r*sin(theta)^2, (4, 4, 3) = r^2*sin(theta)*cos(theta)}

(10)

Of course, you can do the same with Ricci, Riemann, etc. to have versions of them where all the terms with degree > 2 in Omega(r) are discarded. The following is also not something you asked but illustrates well how to work with `&GJcy;` and reusing definitions of other tensors; take for instance Ricci

Ricci[definition]

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

(11)

You can now define an alternative R[mu, nu], say `&Rscr;`[mu, nu]based on `&GJcy;`[mu, nu, alpha]:

subs(Christoffel = G, Ricci = `&Rscr;`, Physics:-Ricci[mu, nu] = Physics:-d_[alpha](Physics:-Christoffel[`~alpha`, mu, nu], [X])-Physics:-d_[nu](Physics:-Christoffel[`~alpha`, mu, alpha], [X])+Physics:-Christoffel[`~beta`, mu, nu]*Physics:-Christoffel[`~alpha`, beta, alpha]-Physics:-Christoffel[`~beta`, mu, alpha]*Physics:-Christoffel[`~alpha`, nu, beta])

`&Rscr;`[mu, nu] = Physics:-d_[alpha](Ѓ[`~alpha`, mu, nu], [X])-Physics:-d_[nu](Ѓ[`~alpha`, alpha, mu], [X])+Ѓ[`~beta`, mu, nu]*Ѓ[`~alpha`, alpha, beta]-Ѓ[`~beta`, alpha, mu]*Ѓ[`~alpha`, beta, nu]

(12)

Define(`&Rscr;`[mu, nu] = Physics:-d_[alpha](Ѓ[`~alpha`, mu, nu], [X])-Physics:-d_[nu](Ѓ[`~alpha`, alpha, mu], [X])+Ѓ[`~beta`, mu, nu]*Ѓ[`~alpha`, alpha, beta]-Ѓ[`~beta`, alpha, mu]*Ѓ[`~alpha`, beta, nu])

`Defined objects with tensor properties`

 

{`&Rscr;`[mu, nu], 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], Ѓ[alpha, mu, nu], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

(13)

Create a macro for R to skip the palette

macro(R = `&Rscr;`)

R[1, 1]

(1/2)*(2*(-2*(diff(beta(r), r))+2*(diff(alpha(r), r)))*exp(-2*beta(r)+2*alpha(r))*(diff(alpha(r), r))*r+2*exp(-2*beta(r)+2*alpha(r))*(diff(diff(alpha(r), r), r))*r+exp(-2*beta(r)+2*alpha(r))*(diff(alpha(r), r))*(-2*cos(theta)^2*Omega(r)*(diff(Omega(r), r))*r^3*exp(-2*alpha(r))-2*cos(theta)^2*Omega(r)^2*r^2*exp(-2*alpha(r))+2*Omega(r)*(diff(Omega(r), r))*r^3*exp(-2*alpha(r))+2*r^2*Omega(r)^2*exp(-2*alpha(r))+2*(diff(alpha(r), r))*r+2*(diff(beta(r), r))*r+4)-2*(-cos(theta)^2*(diff(Omega(r), r))*Omega(r)*r^2*exp(-2*alpha(r))-2*cos(theta)^2*Omega(r)^2*r*exp(-2*alpha(r))+r^2*Omega(r)*(diff(Omega(r), r))*exp(-2*alpha(r))+2*r*Omega(r)^2*exp(-2*alpha(r))+2*(diff(alpha(r), r)))*exp(-2*beta(r)+2*alpha(r))*(diff(alpha(r), r))*r+exp(-2*beta(r))*sin(theta)^2*(r*(diff(Omega(r), r))+2*Omega(r))*(cos(theta)^2*Omega(r)^2*(diff(Omega(r), r))*exp(-2*alpha(r))*r^3+2*cos(theta)^2*Omega(r)^3*exp(-2*alpha(r))*r^2-Omega(r)^2*(diff(Omega(r), r))*exp(-2*alpha(r))*r^3-2*Omega(r)^3*exp(-2*alpha(r))*r^2-2*Omega(r)*(diff(alpha(r), r))*r+r*(diff(Omega(r), r))+2*Omega(r))*r+4*Omega(r)^2*cos(theta)^2*(cos(theta)^2*Omega(r)^2*r^2*exp(-2*alpha(r))-r^2*Omega(r)^2*exp(-2*alpha(r))+1)*r)/r

(14)

From Riemann[definition]

Riemann[definition]

Physics:-Riemann[alpha, beta, mu, nu] = Physics:-g_[alpha, lambda]*(Physics:-d_[mu](Physics:-Christoffel[`~lambda`, beta, nu], [X])-Physics:-d_[nu](Physics:-Christoffel[`~lambda`, beta, mu], [X])+Physics:-Christoffel[`~lambda`, upsilon, mu]*Physics:-Christoffel[`~upsilon`, beta, nu]-Physics:-Christoffel[`~lambda`, upsilon, nu]*Physics:-Christoffel[`~upsilon`, beta, mu])

(15)

applying the same procedure you can get one defined in terms of `&GJcy;`, and using the approach to discard terms used for Christoffel you can discard terms in this one too.

 

``


 

Download Alternative_Christoffel_without_some_terms.mw

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

This is fixed in V.428 or higher of the Physics Updates; you get

 

restart; pde := diff(u(x, t), t) = diff(u(x, t), `$`(x, 2)); bc := u(-Pi, t) = u(Pi, t), (D[1](u))(-Pi, t) = (D[1](u))(Pi, t); pdsolve([pde, bc], u(x, t))

diff(u(x, t), t) = diff(diff(u(x, t), x), x)

 

u(-Pi, t) = u(Pi, t), (D[1](u))(-Pi, t) = (D[1](u))(Pi, t)

 

u(x, t) = exp(-t)*_C3*(exp(I*x)*_C1+exp(-I*x)*_C2), u(x, t) = _C3*(_C1+_C2)

(1)

map(pdetest, [u(x, t) = exp(-t)*_C3*(exp(I*x)*_C1+exp(-I*x)*_C2), u(x, t) = _C3*(_C1+_C2)], [pde, bc])

[[0, 0, 0], [0, 0, 0]]

(2)

``


 

Download fixed_in_v.428.mw

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

Hi

This is fixed in the Physics Updates version 427 or higher. In version 426, unfortunately a file of the version under development got loaded inadvertently, creating an out of sync issue. Thanks for posting the problem.

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

Below I intercalated comments italized, showing how to use casesplit to tackle the problem. Note also that through casesplit you can choose to use rifsimp, DifferentialAlgebra or DifferentialThomas as differential elimination engines. By default, for this problem it uses rifsimp (the non-integer powers are replaced by auxiliary functions satisfying differential equations).

restart; with(PDEtools); declare(u(t, x), A(t), B(t), F(t))

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

 

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

 

` B`(t)*`will now be displayed as`*B

 

` F`(t)*`will now be displayed as`*F

(1)

det_eqs := [2*(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx)) = 0, 2*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx), x))+2*uxx*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), ux), uxx), uxx))+2*ux*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx), uxx)) = 0, 2*(diff(Lambda1(t, x, u, ux, uxx), uxx))*B(t)+(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx))*B(t)*uxx+(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx))*F(t)*u+(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx))*A(t)*u^n*ux = 0, 2*(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx))*ux+2*(diff(diff(Lambda1(t, x, u, ux, uxx), ux), uxx))*uxx-2*(diff(Lambda1(t, x, u, ux, uxx), ux))+2*(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), x)) = 0, 3*(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx))*B(t)+(diff(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx), uxx))*B(t)*uxx+(diff(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx), uxx))*F(t)*u+(diff(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx), uxx))*A(t)*u^n*ux = 0, ux^2*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), u), uxx))+2*ux*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx), x))-ux*(diff(diff(Lambda1(t, x, u, ux, uxx), u), ux))-uxx*(diff(diff(Lambda1(t, x, u, ux, uxx), ux), ux))+uxx^2*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), ux), ux), uxx))+uxx*(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx))+2*uxx*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), ux), uxx), x))+2*uxx*ux*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), ux), uxx))-(diff(diff(Lambda1(t, x, u, ux, uxx), ux), x))+diff(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), x), x) = 0, 2*uxx*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), ux), uxx), uxx))*A(t)*u^n*ux+4*(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), x))*B(t)+4*(diff(diff(Lambda1(t, x, u, ux, uxx), ux), uxx))*B(t)*uxx+2*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), ux), uxx), uxx))*B(t)*uxx^2+4*ux*(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx))*B(t)+2*ux*(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx))*F(t)+2*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx), x))*B(t)*uxx+2*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx), x))*F(t)*u+2*(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx))*A(t)*u^n*n*ux^2/u+2*uxx*(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx))*A(t)*u^n+2*uxx*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), ux), uxx), uxx))*F(t)*u+2*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx), uxx))*A(t)*u^n*ux^2+2*ux*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx), uxx))*B(t)*uxx+2*ux*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx), uxx))*F(t)*u+2*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx), x))*A(t)*u^n*ux = 0, 2*ux*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx), x))*B(t)*uxx+2*ux*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx), x))*F(t)*u+(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), u), uxx))*A(t)*u^n*ux^3+ux^2*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), u), uxx))*B(t)*uxx+ux^2*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), u), uxx))*F(t)*u+2*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx), x))*A(t)*u^n*ux^2-(diff(diff(Lambda1(t, x, u, ux, uxx), u), ux))*A(t)*u^n*ux^2-ux*(diff(diff(Lambda1(t, x, u, ux, uxx), u), ux))*F(t)*u-2*uxx*(diff(Lambda1(t, x, u, ux, uxx), ux))*A(t)*u^n+2*ux^2*(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx))*F(t)+(diff(diff(diff(Lambda1(t, x, u, ux, uxx), ux), ux), uxx))*B(t)*uxx^3+ux^2*(diff(diff(Lambda1(t, x, u, ux, uxx), u), u))*B(t)-(diff(Lambda1(t, x, u, ux, uxx), x))*A(t)*u^n+2*ux*(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), x))*F(t)-ux*(diff(Lambda1(t, x, u, ux, uxx), ux))*F(t)+2*ux*(diff(diff(Lambda1(t, x, u, ux, uxx), u), x))*B(t)+(diff(diff(Lambda1(t, x, u, ux, uxx), ux), x))*B(t)*uxx+2*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), ux), uxx), x))*B(t)*uxx^2+(diff(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), x), x))*F(t)*u-(diff(diff(Lambda1(t, x, u, ux, uxx), ux), x))*F(t)*u+(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx))*B(t)*uxx^2+uxx*(diff(Lambda1(t, x, u, ux, uxx), uxx))*F(t)+(diff(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), x), x))*B(t)*uxx+(diff(Lambda1(t, x, u, ux, uxx), u))*F(t)*u+2*(diff(Lambda1(t, x, u, ux, uxx), u))*B(t)*uxx-uxx*(diff(diff(Lambda1(t, x, u, ux, uxx), ux), ux))*F(t)*u+uxx*ux*(diff(diff(Lambda1(t, x, u, ux, uxx), u), ux))*B(t)+uxx*(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx))*F(t)*u+2*ux*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), ux), uxx))*B(t)*uxx^2+2*uxx*ux*(diff(diff(Lambda1(t, x, u, ux, uxx), ux), uxx))*F(t)+2*uxx*(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), x))*A(t)*u^n+2*uxx*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), ux), uxx), x))*F(t)*u+2*uxx*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), ux), uxx))*A(t)*u^n*ux^2+2*uxx*ux*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), u), ux), uxx))*F(t)*u+3*uxx*(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx))*A(t)*u^n*ux+2*uxx*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), ux), uxx), x))*A(t)*u^n*ux+uxx^2*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), ux), ux), uxx))*A(t)*u^n*ux-uxx*(diff(diff(Lambda1(t, x, u, ux, uxx), ux), ux))*A(t)*u^n*ux+Lambda1(t, x, u, ux, uxx)*F(t)+2*uxx*(diff(diff(Lambda1(t, x, u, ux, uxx), ux), uxx))*A(t)*u^n*n*ux^2/u+3*uxx*(diff(Lambda1(t, x, u, ux, uxx), uxx))*A(t)*u^n*n*ux/u+(diff(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), x), x))*A(t)*u^n*ux-(diff(diff(Lambda1(t, x, u, ux, uxx), ux), x))*A(t)*u^n*ux+(diff(diff(Lambda1(t, x, u, ux, uxx), x), x))*B(t)-(diff(Lambda1(t, x, u, ux, uxx), t))+2*uxx^2*(diff(diff(Lambda1(t, x, u, ux, uxx), ux), uxx))*A(t)*u^n+uxx^2*(diff(diff(diff(Lambda1(t, x, u, ux, uxx), ux), ux), uxx))*F(t)*u+2*(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx))*A(t)*u^n*n*ux^3/u+(diff(Lambda1(t, x, u, ux, uxx), uxx))*A(t)*u^n*n^2*ux^3/u^2-(diff(Lambda1(t, x, u, ux, uxx), uxx))*A(t)*u^n*n*ux^3/u^2+2*(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), x))*A(t)*u^n*n*ux^2/u-(diff(Lambda1(t, x, u, ux, uxx), ux))*A(t)*u^n*n*ux^2/u = 0, diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx) = 0, diff(diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx), uxx) = 0]

The command PDEtools:-casesplit handles no only non-integer powers but also known functions and arbitrary compositions of all of those. Now, your system involves 4 unknown functions

indets(det_eqs, unknown); show

{A(t), B(t), F(t), Lambda1(t, x, u, ux, uxx)}

(2)

Your call to DEtools[rifsimp with CL_multipliers := [Lambda1(t, x, u, ux, uxx)] is equivalent to take A, B  and F as arbitrary functions. For arbitrary values of them, the only solution to your system is

NULL

casesplit(det_eqs, Lambda1, arbitrary = {A, B, F})

`casesplit/ans`([diff(Lambda1(t, x, u, ux, uxx), t) = Lambda1(t, x, u, ux, uxx)*F(t), diff(Lambda1(t, x, u, ux, uxx), x) = 0, diff(Lambda1(t, x, u, ux, uxx), u) = 0, diff(Lambda1(t, x, u, ux, uxx), ux) = 0, diff(Lambda1(t, x, u, ux, uxx), uxx) = 0], [])

(3)

More interesting, you can split this problem into cases if you take all of A, B and F as solving variables ranked lower than Lambda1 (see discussion on rankings in the help page for PDEtools:-casesplit )

casesplit(det_eqs, [Lambda1, {A, B, F}])

`casesplit/ans`([diff(Lambda1(t, x, u, ux, uxx), t) = Lambda1(t, x, u, ux, uxx)*F(t), diff(Lambda1(t, x, u, ux, uxx), x) = 0, diff(Lambda1(t, x, u, ux, uxx), u) = 0, diff(Lambda1(t, x, u, ux, uxx), ux) = 0, diff(Lambda1(t, x, u, ux, uxx), uxx) = 0], [A(t) <> 0, B(t) <> 0]), `casesplit/ans`([diff(diff(Lambda1(t, x, u, ux, uxx), x), x) = (diff(Lambda1(t, x, u, ux, uxx), t)-Lambda1(t, x, u, ux, uxx)*F(t))/B(t), diff(Lambda1(t, x, u, ux, uxx), u) = 0, diff(Lambda1(t, x, u, ux, uxx), ux) = 0, diff(Lambda1(t, x, u, ux, uxx), uxx) = 0, A(t) = 0], [B(t) <> 0]), `casesplit/ans`([diff(diff(Lambda1(t, x, u, ux, uxx), uxx), x) = -(diff(diff(Lambda1(t, x, u, ux, uxx), ux), uxx))*uxx-(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx))*ux+diff(Lambda1(t, x, u, ux, uxx), ux), diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx) = 0, diff(Lambda1(t, x, u, ux, uxx), t) = ((diff(Lambda1(t, x, u, ux, uxx), uxx))*A(t)*u^n*n^2*ux^3+(diff(Lambda1(t, x, u, ux, uxx), ux))*A(t)*u^n*n*ux^2*u+3*uxx*(diff(Lambda1(t, x, u, ux, uxx), uxx))*A(t)*u^n*n*ux*u-(diff(Lambda1(t, x, u, ux, uxx), uxx))*A(t)*u^n*n*ux^3-(diff(Lambda1(t, x, u, ux, uxx), x))*A(t)*u^n*u^2+ux*(diff(Lambda1(t, x, u, ux, uxx), ux))*F(t)*u^2+(diff(Lambda1(t, x, u, ux, uxx), u))*F(t)*u^3+uxx*(diff(Lambda1(t, x, u, ux, uxx), uxx))*F(t)*u^2+Lambda1(t, x, u, ux, uxx)*F(t)*u^2)/u^2, B(t) = 0], [])

(4)

The caseplot option gives a pictorial view of the splitting into cases

casesplit(det_eqs, [Lambda1, {A, B, F}], caseplot)

`========= Pivots Legend =========`

 

p1 = B(t)

 

p2 = A(t)

 

 

`casesplit/ans`([diff(Lambda1(t, x, u, ux, uxx), t) = Lambda1(t, x, u, ux, uxx)*F(t), diff(Lambda1(t, x, u, ux, uxx), x) = 0, diff(Lambda1(t, x, u, ux, uxx), u) = 0, diff(Lambda1(t, x, u, ux, uxx), ux) = 0, diff(Lambda1(t, x, u, ux, uxx), uxx) = 0], [A(t) <> 0, B(t) <> 0]), `casesplit/ans`([diff(diff(Lambda1(t, x, u, ux, uxx), x), x) = (diff(Lambda1(t, x, u, ux, uxx), t)-Lambda1(t, x, u, ux, uxx)*F(t))/B(t), diff(Lambda1(t, x, u, ux, uxx), u) = 0, diff(Lambda1(t, x, u, ux, uxx), ux) = 0, diff(Lambda1(t, x, u, ux, uxx), uxx) = 0, A(t) = 0], [B(t) <> 0]), `casesplit/ans`([diff(diff(Lambda1(t, x, u, ux, uxx), uxx), x) = -(diff(diff(Lambda1(t, x, u, ux, uxx), ux), uxx))*uxx-(diff(diff(Lambda1(t, x, u, ux, uxx), u), uxx))*ux+diff(Lambda1(t, x, u, ux, uxx), ux), diff(diff(Lambda1(t, x, u, ux, uxx), uxx), uxx) = 0, diff(Lambda1(t, x, u, ux, uxx), t) = ((diff(Lambda1(t, x, u, ux, uxx), uxx))*A(t)*u^n*n^2*ux^3+(diff(Lambda1(t, x, u, ux, uxx), ux))*A(t)*u^n*n*ux^2*u+3*uxx*(diff(Lambda1(t, x, u, ux, uxx), uxx))*A(t)*u^n*n*ux*u-(diff(Lambda1(t, x, u, ux, uxx), uxx))*A(t)*u^n*n*ux^3-(diff(Lambda1(t, x, u, ux, uxx), x))*A(t)*u^n*u^2+ux*(diff(Lambda1(t, x, u, ux, uxx), ux))*F(t)*u^2+(diff(Lambda1(t, x, u, ux, uxx), u))*F(t)*u^3+uxx*(diff(Lambda1(t, x, u, ux, uxx), uxx))*F(t)*u^2+Lambda1(t, x, u, ux, uxx)*F(t)*u^2)/u^2, B(t) = 0], [])

(5)

From the pivot information you see the problem split into three cases depending on whether A or B are equal to 0, and when both are different from zero only the first case (equivalent to all of A, B and F being arbitrary). By the way, if this is the determining system of a PDE problem you may be interested in giving a look at the help page for PDEtools:-DEtermininingPDE  and all the related symmetry commands.

 

NULL


 

Download question_for_rif_(reviewed).mw

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

@nm 

Your comments are valuable; it would be helpful however if you could please

  • Transform your reply into a post.
  • Be more specific regarding 3) and 4); itemize the three things you think are more important than anything else for each of them, be clear enough for people to understand, unambiguously, what you are talking about. Giving an example for each of your points always help.

Then, people may or not agree with you, as usual, that is natural, and Maplesoft may or may not have other plans, that is also natural, but people are listening, and chances are that your opinions (yours and others that may add material to your post) are taken into account, entirely or to some point.

Regarding the Heun functions, the original topic of this question, I agree with you nm when you say "given the small [amount of] resources ... I would prefer Maplesoft put its important resources on things that can impact many more users".

For example, this year is one with a heavy dedication to the integral transforms (the inttrans package) in connection with further developments on the exact solutions of PDE & Boundary Conditions, where these transforms are a key part of the strategies. The integral transforms numerical evaluation, differentiation rules, and several fixes and improvements (all that included some versions ago in the Physics Updates for Maple 2019 - and there is more coming) have the potential for significantly bumping up the utility of these valuable functions.

The Heun functions, by the way, are implemented - I wrote the Maple implementation many years ago, and I think it is a thorough implementation. True, their existing numerical evaluation routines can be improved at this point, and the Lame and spheroidal wave functions are not implemented. I know all that. Nobody forgot. What is happening is that other things pop up year after year that appear to me more relevant and of more impact than enhancinthe existing (and in my opinion already sophisticated) numerical evaluation routines of Heun functions, or implementing Lame. I realize others may not agree with me, and that is also natural.

And about that comment "... esoteric parts of the physics package", it looks to me highly subjective - despite the way the comment is written as if it were an indisputable truth. My opinion is that nothing currently found in the Physics package is esoteric. We all need to be respectful of others opinions. And that is all that I'll say in this thread about the original question.

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

The command for applying a change of variables (transformation) to an algebraic expression (differential or not) is PDEtools:-dchange

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

The convention explained in dsolve's help page is that the integration constants introduced by dsolve are of the form _Cn with n an integer. If you change that by something else, odetest has no way to know what is the integration constant that dsolve introduced. That information is relevant to test implicit solutions, as the one you show in your question. As explained in an answer to another of your questions, the technique is simple: solve for the integration constant introduced by dsolve (_C1, identified by odetest because of the _Cn convention) then differentiate with respect to the independent variable (in your example that is x) and you will get 0 modulo the ODE.

In summary: no bug here, and if you want the code to work as indicated in the help page, you cannot change its conventions. Alternatively, you change the conventions, and then test this result manually, not using odetest. The technique is simple, explained in the previous paragraph.

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


 

From the version of Physics you are using I assume you are using Maple 2018?

 

I'm unable to reproduce the problem in Maple 2019, this is what I get:

 

Physics:-Version()

`The "Physics Updates" version in the MapleCloud is 419 and is the same as the version installed in this computer, created 2019, August 24, 16:37 hours, found in the directory /Users/ecterrab/maple/toolbox/2019/Physics Updates/lib/`

(1)

with(Physics)

Setup(coordinates = {X = [t, x, y, z]}, dimension = 4, metric = {(1, 1) = 1, (2, 2) = -1, (3, 3) = -1, (4, 4) = -1}, signature = `+---`, mathematicalnotation = true)

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

 

_______________________________________________________

 

[coordinatesystems = {X}, dimension = 4, mathematicalnotation = true, metric = {(1, 1) = 1, (2, 2) = -1, (3, 3) = -1, (4, 4) = -1}, signature = `+ - - -`]

(2)

Define(H[mu, nu] = Matrix(4, {(1, 1) = 1+h, (1, 2) = h, (2, 2) = -1+h, (3, 3) = -1, (4, 4) = -1}, shape = symmetric))

`Defined objects with tensor properties`

 

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

(3)

H[]

H[mu, nu] = Matrix(%id = 18446744078474560566)

(4)

"LinearAlgebra:-Determinant(rhs(?))"

-1

(5)

H[determinant]

-1

(6)

Note that in Maple 2019 you can also indicate the indices to make clear what determinant are you computing:

H[`~mu`, nu, determinant]

1

(7)

And that is Maple 2019.

Regarding previous versions of Maple, unfortunately I cannot help you with that, but in case this information is of use for you, from reading the code in Maple 2018 I can see it is computing the determinant of H[mu, `~nu`]:

"g_[~mu,~nu] H[nu,rho]"

Physics:-g_[`~mu`, `~nu`]*H[nu, rho]

(8)

Simplify(Physics[g_][`~mu`, `~nu`]*H[nu, rho])

H[rho, `~mu`]

(9)

TensorArray(H[rho, `~mu`])

Matrix(%id = 18446744078374452990)

(10)

"LinearAlgebra:-Determinant(?)"

1

(11)

That is, in Maple 2018 you are getting "det(H) = det(H[nu]^(mu))=det(g[]^(mu,nu)) det(H[nu,rho])."


 

Download It_works_fine_in_2019..mw


Edgardo S. Cheb-Terrab
Physics, Differential Equations andMathematicalFunctions, Maplesoft

  • When odeavisor returns a list of types (or classifications) of an ODE. This must mean the ODE can be of any one of these types?

Yes, when odeavisor "returns a list of types (or classifications) of an ODE" it means the ODE is of all of those types.

  • When Maple is given an ODE to solve, and it can be of number of types, how does it select which type to use to solve the ODE?

See the help page ?dsolve,setup. That page is old, by now incomplete, for instance it does not mention the method of integrating factors for ODEs of differential order > 1, a method as powerful as the Lie symmetry method. See also ?dsolve,education and ?dsolve,algorithms.

  • Is there any detailed document that explains more how odeadvisor determines the type of the ODE? This is a very powerful and a useful command but I can't find in help any hints on how it determines the type of ODE.

I wrote the odeavisor to be a sort of interactive e-book for ODEs: you pass an ODE and it not only determines the ODE types but also pops-up help pages explaining the solving methods available for those types. The point is that I wrote dsolve and the odeavisor at the same time, structuring the code so that the latter could take advantage of the subroutines of the former. The odeadvisor is thus a sort of "dsolve half of the way" plus the set of explanatory help pages on the corresponding solving methods for the ODE you passed, and its ability to pop-up these pages automatically. The implementation details are too many, and also this is not public software, even when at Maplesoft we are incredibly open regarding that - you can always read the code.

In general, about possibly more questions on the structure of Maple commands or the help pages mentioned above: I don't intend to make the help pages mentioned more complete/detailed. The competition copied some of our algorithms, for dsolve, pdsolve and also other stuff, e.g. the assuming and FunctionAdvisor commands, albeit they did a poor job at all that in my opinion. 

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

As pointed out by you, Oliveira, this was an issue visible only when the Physics setting assumingusesAssume was equal to true . The problem was actually within assuming, not  Physics:-Assume, and not noticed before because the old assume command redefines assumed variables, something that the more modern Physics:-Assume does not do. The issue is now fixed at its root and so assuming works as expected with the Physics setting assumingusesAssume. The adjustment is distributed for everybody within the Maplesoft Physics Updates v.416 or higher.

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

Hi

First of all: no, there is no leakage with regards to use and any other package  (Physics included). What this is: most packages have initialization files that set a few things here and there. That is not always documented properly. When you invoke use Package that initialization is executed, and whatever got set remains set after end use. That accounts for what Rouben Rostamian mentioned, about the use of the dot to display derivatives with respect to t . Historically, this was a Maple default, which however changed one or two releases ago, and I decided to keep that default in place when you load Physics, for natural reasons.

Second, there is a bug, not a leakage, with regards to the Physics's default setting 'assumingusesAssume'. When you load (or use) Physics, that setting is set to true, also for obvious reasons (Physics:-Assume does not redefine assumed variables, representing a significant advantage all around). This bug was mentioned as such yesterday, in this other Mapleprimes post by you, Oliveira. But this has nothing to do with leakages or use .

It is important to recall here that the implementation of both old and more modern ideas (e.g. assuming using either assume or the more modern Physics:-Assume), as every portion of code, is only as good as the amount of hours tested / debugged. Finding an issue in the interaction of assuming and Physics:-Assume is rare, I'm glad that the problem got uncovered. I noticed, however, that were not for Carl Love extra comments, you (Oliveira) were not going to mention the origin of your problem (that dsolve was returning only one branch instead of three, in connection with loading Physics and its automatic setting assumingusesAsume = true). Having good mathematical software is a collective endeavour. Reports of things not working as expected, with the necessary details for reproducing them, is perhaps the most important thing / contribution to the project. 

Note after post: That problem in assuming when using Physics with the setting assumingusesAssume is fixed and the adjustment distributed for everybody within the Maplesoft Physics Updates v.416 or higher.

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

Hi Rouben

Good catch, tricky problem, this issue is fixed and the fix available to everybody within the Maplesoft Physics Updates version 411 or higher. What follows is your worksheet, your text, with any added text in italic style. and the new result and its verification.

 

We solve Laplace's equation in the domain a < r and r < b, c < t and t < d
in polar coordinates subject to prescribed Dirichlet data.

Maple produces a solution in the form of an infinite sum,
but that solution fails to satisfy the boundary condition
on the domain's outer arc.  Is this a bug or am I missing
something?

The issue is resolved, the solution, shown below, is correct.

restart;

# kernelopts(version);

Physics:-Version()[2]

`2019, August 7, 17:56 hours, version in the MapleCloud: 411, version installed in this computer: 411.`

(1)

with(plots):

pde := diff(u(r,t),r,r) + diff(u(r,t),r)/r + diff(u(r,t),t,t)/r^2 = 0;

diff(diff(u(r, t), r), r)+(diff(u(r, t), r))/r+(diff(diff(u(r, t), t), t))/r^2 = 0

(2)

a, b, c, d := 1, 2, Pi/6, Pi/2;

1, 2, (1/6)*Pi, (1/2)*Pi

(3)

bc := u(r,c)=c, u(r,d)=0, u(a,t)=0, u(b,t)=t;

u(r, (1/6)*Pi) = (1/6)*Pi, u(r, (1/2)*Pi) = 0, u(1, t) = 0, u(2, t) = t

(4)

We plot the boundary data on the domain's outer arc:

p1 := plots:-spacecurve([b*cos(t), b*sin(t), t], t=c..d, color=red, thickness=5);

 

Solve the PDE:

pdsol := pdsolve({pde, bc});

u(r, t) = Sum((1/3)*sin((3/2)*n1*(-2*t+Pi))*8^n1*(-r^(3*n1)+r^(-3*n1))*(-3+(-1)^n1)/((64^n1-1)*n1), n1 = 1 .. infinity)+Sum(-(1/3)*((-1)^n-1)*sin(n*Pi*ln(r)/ln(2))*(exp((1/6)*Pi*n*(7*Pi-6*t)/ln(2))-exp((1/6)*Pi*n*(Pi+6*t)/ln(2)))/(n*(exp(n*Pi^2/ln(2))-exp((1/3)*n*Pi^2/ln(2)))), n = 1 .. infinity)

(5)

Try pdetest, and verify that what remains is actually 0 in the region of interests (see boundary conditions passed in (4))

remains := pdetest(pdsol, [pde, bc])

[0, Sum(-((1/6)*I)*((-1)^n-1)*(r^(-I*n*Pi/ln(2))-r^(I*Pi*n/ln(2)))/n, n = 1 .. infinity)-(1/6)*Pi, 0, 0, Sum(-(1/3)*sin((3/2)*n*(-2*t+Pi))*(-3+(-1)^n)/n, n = 1 .. infinity)-t]

(6)

zero__r := remains[2]

Sum(-((1/6)*I)*((-1)^n-1)*(r^(-I*n*Pi/ln(2))-r^(I*Pi*n/ln(2)))/n, n = 1 .. infinity)-(1/6)*Pi

(7)

Digits := 20

20

(8)

The following verifies OK, the plot is sufficiently close to 0 in the region of interest, r = a .. b

plot(eval(zero__r, [infinity = 100, Sum = add]), r = a .. b)

 

Next 'remains' to be verified

zero__t := remains[-1]

Sum(-(1/3)*sin((3/2)*n*(-2*t+Pi))*(-3+(-1)^n)/n, n = 1 .. infinity)-t

(9)

The following also verifies OK, the plot is again sufficiently close to 0 in the region of interest, t = c .. d

plot(eval(zero__t, [infinity = 100, Sum = add]), t = c .. d)

 

Try your verification approach just changing 'value' by eval at Sum = add  and using 100 terms instead of just 20; also change the orientation a bit to make the comparison more evident

Truncate the infinite sum at [20] 200 terms, and plot the result:

eval(rhs(pdsol), [infinity = 100, Sum = add])

p2 := plot3d([r*cos(t), r*sin(t), %], r=a..b, t=c..d, orientation=[32, 40, 20]);

 

Here is the combined plot of the solution and the boundary condition.
We see that the proposed solution [after installing the Maplesoft Physics Updates v.411 or higher] does satisfy the boundary condition.

plots:-display([p1, p2], orientation=[32, 40, 20]);

 

``


 

Download mw-6_(reviewed).mw

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

The solution returned by pdsolve is correct, the starting term with n = 0 should be computed taking limits (see vv's reply). On the other hand, I can also see the issue you are pointing at: it is standard to compute the sum by just evaluating the summation index, and if you do so you receive division by 0.

A change in the code that resolves that issue, presenting the special cases computed as branches in a piecewise, is distributed within the Maplesoft Physics updates v.410 and higher, so that after installing that Update you get

f := piecewise(0 <= x and x <= 1, (1/10)*x, 1 < x and x <= 3, 1/10)

pde := diff(u(x, t), t, t)+(4*Pi*(1/3))*(diff(u(x, t), t)) = 16*(diff(u(x, t), x, x))
ic := (D[2](u))(x, 0) = 0, u(x, 0) = f

bc := u(0, t) = 0, (D[1](u))(3, t) = 0

 

pdsolve([pde, ic, bc])

u(x, t) = Sum(piecewise(n = 0, (4/5)*sin((1/6)*Pi*x)*(Pi*t+3/2)*exp(-(2/3)*Pi*t)/Pi^2, (3/10)*sin((1/6)*(1+2*n)*Pi*x)*((2*n^(1/2)*(n+1)^(1/2)+I)*exp(((2/3)*I)*Pi*(-2*n^(1/2)*(n+1)^(1/2)+I)*t)-(-2*n^(1/2)*(n+1)^(1/2)+I)*exp(((2/3)*I)*Pi*(2*n^(1/2)*(n+1)^(1/2)+I)*t))*(3^(1/2)*sin((1/3)*Pi*n)+cos((1/3)*Pi*n))/(n^(1/2)*(n+1)^(1/2)*(1+2*n)^2*Pi^2)), n = 0 .. infinity)

(1)

Another example

f := 'f'

pde := diff(u(x, t), t, t)+2*(diff(u(x, t), t)) = diff(u(x, t), x, x)
ic := (D[2](u))(x, 0) = 0, u(0, t) = 0, u(x, 0) = f(x)

bc := u(0, t) = 0, u(Pi, t) = 0

sol := `assuming`([pdsolve([pde, ic, bc])], [t > 0])

u(x, t) = Sum(piecewise(n = 1, 2*(int(sin(x)*f(x), x = 0 .. Pi))*exp(-t)*(t+1)*sin(x)/Pi, (Int(sin(x*n)*f(x), x = 0 .. Pi))*((-1+(-n^2+1)^(1/2))*exp(-((-n^2+1)^(1/2)+1)*t)+exp((-1+(-n^2+1)^(1/2))*t)*((-n^2+1)^(1/2)+1))*sin(x*n)/((-n^2+1)^(1/2)*Pi)), n = 1 .. infinity)

(2)

``


 

Download piecewise_solution_with_special_cases.mw

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

Hi
I see no bug here. In
> ee := (((Q(a)^3)^(5/4))^(15/7))^(6/8);

you have a four times nested anything^rational type object. So, for the call on the left-hand side of what you posted, for any P, the call

> subsindets(ee, anything^rational, P)

will be applied, from the inside to the outside, 4 times, resulting in

    P(P(P(P(Q(a)^3)^(5/4))^(15/7))^(3/4))

But then on the right-hand side of what you posted, I see subsindets applied to objects of type specfunc(Q)^rational, that is a different thing, and within ee I see only one object of that type, so the call 

> subsindets(ee, specfunc(Q)^rational, P);

results, again as expected, in

    ((P(Q(a)^3)^(5/4))^(15/7))^(3/4)

Make now P be

> P := proc(x) if type(x, specfunc(Q)^rational) then 'Q(x)' else 'x' fi end

and you see the two results you got, all this is as expected.

Note as well that if you don't want subsindets to be applied recursively, you can use its option 'flat'. Also, to see subsindets at work (after you have assigned the procedure P) it suffices to trace P, entering trace(P); then call subsindets again

> subsindets(ee, anything^rational, P);
{--> enter P, args = Q(a)^3
<-- exit P (now at top level) = Q(Q(a)^3)}
{--> enter P, args = Q(Q(a)^3)^(5/4)
<-- exit P (now at top level) = Q(Q(Q(a)^3)^(5/4))}
{--> enter P, args = Q(Q(Q(a)^3)^(5/4))^(15/7)
<-- exit P (now at top level) = Q(Q(Q(Q(a)^3)^(5/4))^(15/7))}
{--> enter P, args = Q(Q(Q(Q(a)^3)^(5/4))^(15/7))^(3/4)
<-- exit P (now at top level) = Q(Q(Q(Q(Q(a)^3)^(5/4))^(15/7))^(3/4))}

while for the other type you used you have

> subsindets(ee, specfunc(Q)^rational, P);
{--> enter P, args = Q(a)^3
<-- exit P (now at top level) = Q(Q(a)^3)}

In summary that I see no bug here.

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

 

 

Have in mind that pFq, more than most functions actually, satisfy a myriad of identities (try FunctionAdvisor(identities, hypergeom([a, b], [c], z)), each of which has an apparently different integral representation. This is the same simplification problem mentioned in your previous post, two things that are the same but look different, making zero recognition difficult, not algorithmically possible in all cases.

Besides, the solution returned by dsolve uses Intat, not Int. There is no conversion to Intat, a more intelligent comand (think of it as an integral evaluated at a point - the equivalent to a derivative evaluated at a point represented in Maple by the D command).

As for your comment about odetest "does not give zero", its help page explains that a solution could be correct but, because of zero recognition, odetest will not find zero. I thought it was clear in the replies you received for the previous post you made today. I also showed, in one of those replies, what you could do in those cases: test the solution with a plot (equivalent to saying "numerically"). You could copy the input lines I posted in reply to your previous post and try to apply them for this your hand_solution.

 

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