Ask Your Question
1

desolve: wrong solution?!

asked 2020-05-31 11:52:19 -0500

anonymous user

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 gravatar imageEmmanuel Charpentier ( 2020-06-06 08:59:16 -0500 )edit

2 answers

Sort by ยป oldest newest most voted
1

answered 2020-06-03 15:12:19 -0500

Emmanuel Charpentier gravatar image

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 more
0

answered 2020-06-08 13:24:12 -0500

curios_mind gravatar image

updated 2020-06-08 14:36:37 -0500

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 more

Comments

Briefly :

  1. 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.

  2. 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.

  3. 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 gravatar imageEmmanuel Charpentier ( 2020-06-08 17:07:01 -0500 )edit

Your Answer

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

Add Answer

Question Tools

1 follower

Stats

Asked: 2020-05-31 11:52:19 -0500

Seen: 107 times

Last updated: Jun 08