ecterrab

12714 Reputation

24 Badges

17 years, 263 days

MaplePrimes Activity


These are answers submitted by ecterrab

For vector calculus problems in Cartesian, cylindrical or spherical coordinates/basis, mostly always Physics:-Vectors provides a simpler approach.

with(Physics:-Vectors)

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

(1)

Your vector field - you don't need, but it is more visual if you give it a name:

A_ := f(r, theta, phi)*_theta

f(r, theta, phi)*_theta

(2)

This is your PDE: note that diff (actually, Physics:-Vectors:-diff) knows about the coordinates dependency of the spherical unit vectors, so differentiating you get

pde := diff(A_, theta) = 0

-f(r, theta, phi)*_r+(diff(f(r, theta, phi), theta))*_theta = 0

(3)

If you want to solve this as a system of equations for f(r, theta, phi), you need to rewrite the vector basis as constant unit vectors, i.e. Cartesian, since otherwise, e.g. the spherical unit vectors depend on theta. So,

ChangeBasis(-f(r, theta, phi)*_r+(diff(f(r, theta, phi), theta))*_theta = 0, Cartesian)

cos(phi)*((diff(f(r, theta, phi), theta))*cos(theta)-f(r, theta, phi)*sin(theta))*_i+sin(phi)*((diff(f(r, theta, phi), theta))*cos(theta)-f(r, theta, phi)*sin(theta))*_j+(-f(r, theta, phi)*cos(theta)-sin(theta)*(diff(f(r, theta, phi), theta)))*_k = 0

(4)

The PDE system here is obtained by taking the components; the simplest logical way to get them is to take the scalar produt with the unit vectors `#mover(mi("i"),mo("∧"))`, `#mover(mi("j"),mo("∧"))`, `#mover(mi("k"),mo("∧"))`

pde_system := [_i.(cos(phi)*((diff(f(r, theta, phi), theta))*cos(theta)-f(r, theta, phi)*sin(theta))*_i+sin(phi)*((diff(f(r, theta, phi), theta))*cos(theta)-f(r, theta, phi)*sin(theta))*_j+(-f(r, theta, phi)*cos(theta)-sin(theta)*(diff(f(r, theta, phi), theta)))*_k = 0), _j.(cos(phi)*((diff(f(r, theta, phi), theta))*cos(theta)-f(r, theta, phi)*sin(theta))*_i+sin(phi)*((diff(f(r, theta, phi), theta))*cos(theta)-f(r, theta, phi)*sin(theta))*_j+(-f(r, theta, phi)*cos(theta)-sin(theta)*(diff(f(r, theta, phi), theta)))*_k = 0), _k.(cos(phi)*((diff(f(r, theta, phi), theta))*cos(theta)-f(r, theta, phi)*sin(theta))*_i+sin(phi)*((diff(f(r, theta, phi), theta))*cos(theta)-f(r, theta, phi)*sin(theta))*_j+(-f(r, theta, phi)*cos(theta)-sin(theta)*(diff(f(r, theta, phi), theta)))*_k = 0)]

[cos(phi)*((diff(f(r, theta, phi), theta))*cos(theta)-f(r, theta, phi)*sin(theta)) = 0, sin(phi)*((diff(f(r, theta, phi), theta))*cos(theta)-f(r, theta, phi)*sin(theta)) = 0, -f(r, theta, phi)*cos(theta)-sin(theta)*(diff(f(r, theta, phi), theta)) = 0]

(5)

pdsolve(pde_system)

{f(r, theta, phi) = 0}

(6)

In brief, unless I am missing something here, the solution to your equation diff(pde, theta) = 0 is f(r, theta, phi) = 0.

 

Download spherical_derivative_using_Physics_Vectors.mw

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

The answer to this question is the same one I gave to a question you asked yesterday about the same thing. Recalling the gist of that answer, "If your ODE is solvable using the method for class XXX of ODEs, then the ode is classified as of class XXX." From your two questions about the same, I suppose this sentence should be written explicitly on the odeadvisor help page.

The issue underlying your question here: there is, of course, a distinction between "computer syntactic" versus "mathematical normal form in a solving context." You seem to be disregarding this distinction, while all solving commands make this distinction and take the second one as representing your input expression. The ODE you are presenting now to an ODE solving (classifying) command,  (x+y(x))*diff(y(x),x) = 0, as an ODE, it is normalized to diff(y(x),x) = 0, and therefore of course it is solved by dsolve using the method for linear ODEs (try dsolve(ODE, [linear]) to see that), thus classified as linear by the odeadvisor, no bug here.

A not-equal, however similar situation you have when algebraically solving ee := (a*y^2 + b*y)/y = 0 with respect to y. Is this an expression linear in y? From the computer syntactic point of view, of course not. Input type(ee, linear(y)) , the answer is false. But in the solving context, when you pass it to solve, it is first normalized; as such, ee represents (a*y+b) = 0 and therefore is solved as a linear equation.

Moving then to the other question of your post: you could check whether the "ODE input expression" is syntactically linear in y' and if so, isolate y' before sending it to DEtools:-convertAlg. There are, of course, many other possible ways of addressing your computational need, but that one will match what dsolve does (and therefore what the odeadvisor does, given the first paragraph above).

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

Hi
First, the references for most of the classifications are at the end of the help page ?odeadvisor, and in particular, for this "classify ODEs" project, which is a subproduct of the new dsolve for Maple of 1997, I used Kamke's book. I say "most of the classifications" because the odeadvisor contains all those found in the textbook references but also several other ones - e.g. using symmetry patterns, or for Abel ODEs - that we identified ourselves as more general - entirely solvable - classes, that were (still are) not mentioned in textbooks, only in our scientific published papers - see ?dsolve,references.

And then there is what in Education is called "instrumental knowledge" which I summarize here - e.g., as: if the method for systematically solving dAlembert ODES can solve your example, then your example is of dAlembert type. That doesn't mean "Maple uses its own definition of dAlembert ODEs." The definition, as said, is the one shown in Kamke's book. In this context, to suggest - e.g. - that f = 0 or g = 0 are impediments to classify the ODE as dAlembert is not appropriate. If those values of f or g make the ODE solvable by other methods (e.g. missing x), that changes nothing: the ODE is solvable by more than one method, belongs to more than one class.

And how is the odeavisor used "instrumentally" in Maple? When writing dsolve (1997), I organized the computational flow around a table of methods (see ?dsolve,setup), where new methods could be plugged, removed, or their ordering changed at will, even by a user. Each method consists of three subroutines:

  • The manager routine (receives the problem in some established format, returns NULL or a solution)
  • The odeavisor routine (tells whether the problem passed is or not solvable by the method, and if so, on the way it collects the relevant solving information)
  • The integrate routine (given the OK by the odeavisor routine of the method and relevant information, it proceeds to integrate the ODE using instrumental formulas).

For example, for dAlembert ODEs, these routines are: `odsolve/dAlembert`, `odeadv/dAlembert` and `odsolve/dAlembert/integrate`. You can display these routines on your screen. The odeadvisor only uses the odeadv routine to tell you whether an ODE is or not of that type.

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

 

Hi @Rouben Rostamian@vv;

Download SpecializeArbitraryFunctions_-_fixed.mw
 

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

This was a subtle bug in the routines deciding a convenient ranking to tackle the system. The problem is fixed, and the fix is distributed as usual for everybody using the latest Maple version (2022) within the Maplesoft Physics Updates v.1312 or newer. With the fix in place, pdsolve returns a solution. I also see you are using Maple 2018. Unfortunately, I have no way to make the fix retroactive to previous Maple versions, so in your case I think Carl's solution is the way to go.

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

Download Physics_tensors.mw

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

Hi

The package you mentioned is not a Maplesoft package but one coded by Maple users. I am not sure how well know that package is, but I am not familiar with it. On the other hand, Maple comes, right from the box, with three packages for these purposes you mention: CalculusOfVariations, DifferentialGeometry and Physics. I suggest you try your calculations with the existing Maple package, and if you find something missing, please post that as a question about the existing Maple packages. Most probably, you will receive answers right away.

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

The transformation y(x) = arccos(u(t)) (see ?PDEtools:-dchange) removes the cos(y(x)) from the problem. Still, it remains sqrt(1 - u(t)^2); you'd need to experiment further to remove that to obtain a purely rational equation. In any case, after removing cos(y(x)) you see this equation admits point symmetries (see ?DEtools:-gensym). That doesn't solve the problem but is a good prognostic: if you happen to compute a solution to the output of gensym, that is, if you find a symmetry, you can use it to construct a solution, either using dsolve with one of its symmetry options (see ?dsolve,details) or directly PDEtools:-Invariant solutions (check its help page). You may also want to consider the determination of an integrating factor (here too, see ?gensym), and if you succeed, see ?dsolve,options to use it to compute a solution.

In general, you 1) always try to find a transformation that makes the ODE simpler. Then either you find 2) a symmetry or an integrating factor, and if you succeed, 3) you can use it to integrate the problem. You can use the computer for all 1), 2) and 3).

 

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

Use dsolve's option useint and you get a solution. However, instead of combine, use simplify(solution, size). The solution is large in size. The reason for using this option: some integrals appear in the solution without the option useint, and those integrals make it difficult to solve for the integration constants to match the initial conditions.

I note as well that a solution involving the integrals (so without the option useint) is at reach if you pass only 4 of the 6 initial conditions, then you may want to compute the integral using value, finally adjusting the two remaining integration constants interactively.

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

This is as @acer said; for this problem, solve returns something that, unfortunately, results in RootOf(0) when simplified. The problem is caught now, and PDEtools:-Solve returns NULL. This change is distributed to everybody using Maple 2022 within the Maplesoft Physics Udpates v.1288 or newer.

About this other issue, PDEtools:-Solve(sin(t)^2 + cos(t)^2 - 1, t);, mentioned by @Carl Love, that is an entirely different problem, with no RootOf involved. This other problem is also fixed. The fix is distributed in the same Physics Updates v.1288 or newer.

By the way, PDEtoosl:-Solve does not have another solver within. What PDEtools:-Solve has is several additional strategies extending the applicability of solve and some other solving Maple commands. Recalling from it's help page,

The PDEtools:-Solve command computes the value of solving_variables that solves a system of equations sys. The system being solved can involve algebraic or differential equations, or both. You can request an exact (default), numeric, or series solution (respectively use the option numeric or series). In this sense, Solve is a unified command that understands when to call solve, fsolve, dsolve, or pdsolve according to your input, also facilitating the analysis of different types of solutions by just adding the keywords series or numeric.

In addition ....[etc].

 

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

Hi Rouben,

After with(Vectors), type(A_, PhysicsVectors) works as you expect. I need to document all the Physics:-Library:-PhysicsTypes ... for which nops(with(Physics:-Library:-PhysicsTypes)) -> 114. These are the types used internally by the Physics package routines and are made available for programming purposes together with the Physics:-Library (this one is documented), which contains the internal routines of the Physics package, also documented for programming with them.

By the way, this type should also be mentioned in the Physics:-Vectors help pages - at this moment, it is mentioned on the page ?convert,PhysicsVectors.

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

 

You know the form of Dirac matrices in a Minkowski space with metric  (+++-) , either in a flat spacetime or in the local system of references implemented as the tetrad system.

So take the transformation that transforms the eta_ Minkowski metric (+++-) into (++--), you asked about the same in a previous question I answered a week or so ago, and apply it to the standard form of each Dirac matrix; then use Physics:-Define to define a new tensor (with spacetime or tetrad indices, depending on your problem) whose components are those transformed Dirac matrices. See the page ?Physics,Define and ?Physics,Tensors for the different ways of defining such a tensor.

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

 

This is your example, showing the ODE explicitly

 

ode := sqrt(y(x))+(1+x)*(diff(y(x), x)) = 0; ic := y(0) = 1

y(x)^(1/2)+(1+x)*(diff(y(x), x)) = 0

(1)

To start with, whenever you see something like sqrt(y(x)) you could expect issues: the implicit solution is basically always correct and the explicit solution has problems, because isolating y(x) may require ignoring branch cuts and the like.

 

So let's first see the implicit solution

implicit_sol := dsolve([ode, ic], implicit)

y(x)^(1/2)+(1/2)*ln(1+x)-1 = 0

(2)

odetest(implicit_sol, [ode, ic])

[0, 0]

(3)

 

Good. Try isolating y(x)

solve(implicit_sol, {y(x)})

{y(x) = (1/4)*ln(1+x)^2-ln(1+x)+1}

(4)

odetest({y(x) = (1/4)*ln(1+x)^2-ln(1+x)+1}, [ode, ic])

[(1/2)*csgn(ln(1+x)-2)*ln(1+x)-csgn(ln(1+x)-2)+(1/2)*ln(1+x)-1, 0]

(5)

You see the problem. In fact, if you insert (4) into (2) you see (4) is not really a solution of (2) but for some regions only

simplify(eval(implicit_sol, {y(x) = (1/4)*ln(1+x)^2-ln(1+x)+1}))

(1/2)*(csgn(ln(1+x)-2)+1)*(ln(1+x)-2) = 0

(6)

In the above I would expect 0, but the above is zero only when csgn(ln(1+x)-2) = -1

eval((1/2)*(csgn(ln(1+x)-2)+1)*(ln(1+x)-2) = 0, csgn = -1)

0 = 0

(7)

Now that the problem is understood, let's see your solution``

your_sol := y(x) = (1/4)*(ln(1+x)-2)^2

y(x) = (1/4)*(ln(1+x)-2)^2

(8)

simplify(eval(implicit_sol, your_sol))

(1/2)*(csgn(ln(1+x)-2)+1)*(ln(1+x)-2) = 0

(9)

As expected, it has the same problem, valid only when csgn(ln(1+x)-2)+1 = -1, More important: the "implicit form" of your_solution that you construct taking lhs-rhs is also not-a-solution, let's show it

(lhs-rhs)(your_sol) = 0

y(x)-(1/4)*(ln(1+x)-2)^2 = 0

(10)

Compare with the actual solution computed by dsolve

implicit_sol

y(x)^(1/2)+(1/2)*ln(1+x)-1 = 0

(11)

You see that you can go from (11) to (10) by squaring, but not the other way around taking square roots

sqrt(y(x)) = -(ln(1 + x)/2 - 1)

y(x)^(1/2) = -(1/2)*ln(1+x)+1

(12)

map(`^`,(12), 2);  # going from (11) to (10)

 

y(x) = (-(1/2)*ln(1+x)+1)^2

(13)

map(sqrt, (13));   # going back taking sqrt you do not get (11)

y(x)^(1/2) = ((-(1/2)*ln(1+x)+1)^2)^(1/2)

(14)

The reason being that x <> sqrt(x^2), as we know.

 

So your question could be reformulated as "when testing not-a-solution to an ODE + IC, why am I receiving different not-zero output from odetest depending on how I sent the not-a-solution ?"

 

The literal answer is the one you gave, of course, different sectors of the program are used when you pass an explicit instead of an implicit solution. The rationale for that: it is frequently the case that implicit solutions test OK, while explicit forms of them obtained e.g. thinking that x = sqrt(x^2) of course do not test OK. We want to be able to distinguish between these two cases. And for that we use different computational strategies: when testing implicit solutions we do not isolate the unknown, in that way avoid hitting wrong assumptions like NULLx = sqrt(x^2).

 

In summary: in cases like this one you post, when testing not-a-solution, you can expect different not-zero output from odetest depending on whether you pass the not-a-solution in explicit form or taking lhs-rhs

 

One could further ask: "Why is solve returning a solution that assumes NULLx = sqrt(x^2)? " The answer, I imagine, is historical reasons. Check for instance Mathematica: it also returns a solution that assumes  NULLx = sqrt(x^2) (although it presents a Warning message indicating the solution may be not valid for some values of x ). Maple's solve could present a warning too

solve(sqrt(x^2) = X^2, {x});  # this is ok

{x = X^2}, {x = -X^2}

(17)

This next is ok too in that it is better than NULL, but could present a warning: the solution is not valid for all values of X2.

solve(sqrt(x2) = X2, {x2});

{x2 = X2^2}

(18)

NULL


 

Download Your_problem.mw

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

The default signatures implemented in the Physics package have only one sign different, corresponding to the timelike component of the metric. However, it is possible to work with an arbitrary signature, where by that I mean an arbitrary (symmetric) form of the tetrad metric eta[a, b] of a local flat system of references. In particular, you mentioned the signature "(++--)", so below I derive that case but the procedure shown is the way to go in general.

 

with(Physics); with(Tetrads)

_______________________________________________________

(1)

Not necessary, but to simplify steps, set the signature as close as possible to the one you want

Setup(signature = "+++-")

 

At this point I could do all what follows with a Minkowski metric - the procedure is the same - but it would look artificially easy, so let's complicate only one bit by introducing some curvature, the Schwarzschild

metric, so that the metric (setting it and showing it at the same time via "g_[sc])", tetrad and tetrad metric have for components

g_[sc]

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

(2)

   

The tetrad

e_[]

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

(3)

The orthonormal tetrad metric, a.k.a "the signature"

eta_[]

Physics:-Tetrads:-eta_[a, b] = Matrix(%id = 36893488153385350796)

(4)

We want the latter, (4), to be "something else". In particular, you asked for `&eta;__a,b` = (Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -1, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = -1})). This case is easy, everything is diagonal. If not, the form of A that follows needs to be more general. So for this case search for a transformation matrixA of this form

A := Matrix(4, symbol = a, shape = diagonal)

Matrix(%id = 36893488153385343204)

(5)

The "new" `&eta;__a,b` is computed via

A.rhs(eta_[]).LinearAlgebra:-Transpose(A)

Matrix(%id = 36893488153385329836)

(6)

By eye, the solution A you are looking for involves the imaginary unit I and is

A := Matrix(4, {(1, 1) = 1, (2, 2) = 1, (3, 3) = I, (4, 4) = 1})

Matrix(%id = 36893488153385314052)

(7)

Check it out: this is the new eta

Eta[a, b] = A.rhs(eta_[]).LinearAlgebra:-Transpose(A)

Eta[a, b] = Matrix(%id = 36893488153385305380)

(8)

Not necessary, but to facilitate a posterior verification, define this Eta[a, b], as a tensor

"Define(?)"

{Physics:-D_[mu], Physics:-Dgamma[mu], Eta[a, b], 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:-gamma_[i, j], Physics:-Tetrads:-gamma_[a, b, c], Physics:-Tetrads:-l_[mu], Physics:-Tetrads:-lambda_[a, b, c], Physics:-Tetrads:-m_[mu], Physics:-Tetrads:-mb_[mu], Physics:-Tetrads:-n_[mu], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

(9)

Now, from the relationship that defines the tetrad `&efr;`[a, mu] in terms of the tetrad metric eta[a, b]

e_[definition]

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

(10)

you see that the transformation defining your new tetrad is

E[a, mu] = A.rhs(e_[])

E[a, mu] = Matrix(%id = 36893488151943554100)

(11)

Again, not necessary, but to facilitate manipulations, define it

"Define(?)"

{Physics:-D_[mu], Physics:-Dgamma[mu], E[a, mu], Eta[a, b], 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:-gamma_[i, j], Physics:-Tetrads:-gamma_[a, b, c], Physics:-Tetrads:-l_[mu], Physics:-Tetrads:-lambda_[a, b, c], Physics:-Tetrads:-m_[mu], Physics:-Tetrads:-mb_[mu], Physics:-Tetrads:-n_[mu], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

(12)

So this is the relationship that must be satisfied by the new tetrad and new tetrad metric

"E[a, mu] * E[b,~mu] = Eta[a,b]"

E[a, mu]*E[b, `~mu`] = Eta[a, b]

(13)

TensorArray(E[a, mu]*E[b, `~mu`] = Eta[a, b])

Matrix(%id = 36893488153401475188)

(14)

Now that everything matches, just set the tetrad and the tetrad metric as these ones derived above, equations (8) and (11)

"Setup(tetradmetric= rhs(?), tetrad = rhs(?))"

[tetrad = {(1, 1) = -I*r^(1/2)/(2*m-r)^(1/2), (2, 2) = r, (3, 3) = I*r*sin(theta), (4, 4) = -I*(2*m-r)^(1/2)/r^(1/2)}, tetradmetric = {(1, 1) = 1, (2, 2) = 1, (3, 3) = -1, (4, 4) = -1}]

(15)

Now your tetrad and tetrad metric are set as you want, you are working with a signature "(++--)". You can check it out as follows

e_[]

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

(16)

eta_[]

Physics:-Tetrads:-eta_[a, b] = Matrix(%id = 36893488153385358508)

(17)

e_[definition]

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

(18)

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

Matrix(%id = 36893488153377682244)

(19)

Likewise,

e_[a, mu]*e_[a, nu] = g_[mu, nu]

Physics:-Tetrads:-e_[a, mu]*Physics:-Tetrads:-e_[`~a`, nu] = Physics:-g_[mu, nu]

(20)

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

Matrix(%id = 36893488153356963580)

(21)

NULL


 

Download arbitrary_signature.mw

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

dsolve(ode, useInt);  # tests OK, the output involves uncomputed Int(1/(y^2 - 1)^(2/3), y)
value(%);                     # the solution you show, it is problematic

Summarizing, the issue is with the value of Int(1/(y^2 - 1)^(2/3), y).

By the way, dsolve's options useInt and implicit are the ones that help disentangling an issue in dsolve from issues in int or solve.

2 3 4 5 6 7 8 Last Page 4 of 52