desolve: wrong solution?!

asked 2020-05-31 18:52:19 +0200

Anonymous

Hello,

Consider the following differential equation

y=function('y')(x)
eq=6*diff(y,x)-2*y==x*y^4
sol=desolve(eq,y,show_method=True)

This code gives the following solution

[e^(1/3*x)/(-1/2*(x - 1)*e^x + _C)^(1/3), 'bernoulli']

When I solve this with Mathematica it gives me 3 solutions for y(x) (for each root, as expected).

Question #1: Why Sage gives one solution? When I solve this differential equation analytically, I would multiply both sides by y^(-4) (assuming y!=0) and then change the variable y^3=v. The differential equation becomes a 1st order linear differential equation in the dependent variable v and independent variable x. After solving for v, the cube-root brings the 3 solutions Mathematica is giving.

Question #2: Ok, this is really bugging me the most: Consider the same differential equation but this time let's also give the initial condition y(0)=-2

y=function('y')(x)
eq=6*diff(y,x)-2*y==x*y^4
sol=desolve(eq,y,show_method=True,ics=[0,-2])

This gives the following

[e^(1/3*x)/(-1/2*(x - 1)*e^x - 1/2)^(1/3), 'bernoulli']

However, this function blows up at x=0! How can I have a solution by calling desolve with the initial condition y(0)=-2 and end up with a solution that doesn't satisfy the initial condition ?!

Question #3: I tried to figure out the constant by myself with the following code:

y=function('y')(x)
eq=6*diff(y,x)-2*y==x*y^4
sol=desolve(eq,y,show_method=True)
c=sol[0].variables()[0] #<--- this refs _C
icsol=solve(sol[0].subs(x==0)==-2,c)

icsol gives me

[(_C + 1/2)^(1/3) == (-1/2)]

Why doesn't it give as [_C==...] ?

Question #4: For the code given in Question #3, if I try to substitute icsol back in the solution as

sol[0].subs(icsol)

I got sol[0] thrown back at me unchanged since there is no (_C + 1/2)^(1/3) in sol[0].

How can I even substitute back the constant C in the solution?

Thanks in advance for your help.

edit retag flag offensive close merge delete

Comments

I don't understand your "manual solution". Could you amplify ?

EDIT : Never mind : I've rebuilt it.

Emmanuel Charpentier ( 2020-06-06 15:59:16 +0200 )edit

add a comment

2 Answers

Sort by » oldest newest most voted

answered 2020-06-03 22:12:19 +0200

Emmanuel Charpentier
7057 ●6 ●46 ●137

Partial answer : Mathematica's and Sage's algorithms give the same solutions :

sage: y=function("y")(x)
sage: SSol=desolve(6*diff(y,x)-2*y==x*y^4,y) ; SSol
e^(1/3*x)/(-1/2*(x - 1)*e^x + _C)^(1/3)
sage: MSol=mathematica("DSolve[6*y'[x]-2*y[x]==x*y[x]^4, y[x], x]") ; MSol
{{y[x] -> ((-2)^(1/3)*E^(x/3))/(-E^x + E^x*x - 2*C[1])^(1/3)}, 
 {y[x] -> -((2^(1/3)*E^(x/3))/(-E^x + E^x*x - 2*C[1])^(1/3))}, 
 {y[x] -> -(((-1)^(2/3)*2^(1/3)*E^(x/3))/(-E^x + E^x*x - 2*C[1])^(1/3))}}

Mathematica's answer is not (yet) easily translatable to Sage. However, one can use Mathematica to show that the three solutions given by Mathematica differ only by a multiplicative constant :

sage: KK=mathematica("E^(x/3)/(-E^x + E^x*x - 2*C[1])^(1/3)")
sage: [MSol[u][1][2]/KK for u in (1..3)]
[(-2)^(1/3), -2^(1/3), -((-1)^(2/3)*2^(1/3))]

Those three values are Mathematica's way of deboting the three cube roots of -2 :

sage: [mathematica.N(MSol[u][1][2]/KK) for u in (1..3)]
[0.6299605249474367 + 1.0911236359717214*I,
 -1.2599210498948732,
 0.6299605249474363 - 1.0911236359717216*I]
sage: [(MSol[u][1][2]/KK)^3 for u in (1..3)]
[-2, -2, -2]

With a bit of notational changes, the factor common to Mathematica's three solutions appears very close to the Sage's solution :

sage: C=function("C")
sage: var("_C")
_C
sage: SSol.subs(_C==C(2))
e^(1/3*x)/(-1/2*(x - 1)*e^x + C(2))^(1/3)
sage: KK.sage(locals={"C":C,"E":e}).factor()
e^(1/3*x)/((x - 1)*e^x - 2*C(1))^(1/3)

up to a factor (-1/2)^(1/3) appearing in the denominator, which would be a factor (-2)^(1/3) in the numerator.

Sage appears to denote by (-1/2)^(1/3) any quantity whose cube is -1/2. Up to this notational quirk, Sane's and Mathematica's solutions are identical.

edit flag offensive delete link

add a comment

answered 2020-06-08 20:24:12 +0200

curios_mind
115 ●1 ●5 ●15

updated 2020-06-08 21:36:37 +0200

Thanks for the answer. I see why Sagemath gives one answer and Mathematica gives 3. I tried to solve this using "fricas", it gave me the 3 roots.

Though I still think when the initial conditions are invoked, desolve solves it wrong with maxima, since the denominator is zero. The initial condition can be satisfied by one of the cube roots (checked with Mathematica)

Is there a way to get the (in this example--) cube roots from desolve result? or perhaps I should open up another thread for this.

edit flag offensive delete link

Comments

Briefly :

desolve(6*diff(y(x),x)-2*y(x)==x*y(x)^4), y(x), ics=[0,-2]) definitely gives a wrong answer. It's a bug.
desolve(6*diff(y(x),x)-2*y(x)==x*y(x)^4), y(x)) gives a (correct) general solution, to be exploited with keeping in mind the ambiguity of (-2^(1/3), thus allowing to find the "right" constant (and the "right" branch) satisfying your boundary condition. see this post on sage-support.
FWIW, Sympy can solve (correctly) the general case, but fails fto find the solution for boundary condition.

I want to peruse the known ODE bugs of Sage/Maxima before filing a new ticket.

Emmanuel Charpentier ( 2020-06-09 00:07:01 +0200 )edit

add a comment

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

desolve: wrong solution?!

Comments

2 Answers

Comments

Your Answer

Question Tools

Stats

Related questions

desolve: wrong solution?! edit

Comments

2 Answers

Comments

Your Answer

Question Tools

Stats

Related questions

desolve: wrong solution?!