Does Mathematica evaluate Sqrt[1] to +1 or -1 in this differential equation?

Question

Bug introduced in 8.0 or earlier and persisting through 14.0 or later

---

This is known ode that DSolve generates a wrong extra solution.

Been there since 2010. I do not know if this was posted about here before or not. If it was, will delete this question. I did a search here and did not find one.

DSolve gives two solutions to this ode. But only one is correct.

I've been trying to guess how the wrong solution could have been generated and at what point.

ClearAll[y,x]
DSolve[{y'[x]==2*y[x]* (x *Sqrt[y[x]]-1),y[0]==1},y[x],x]

The first solution above should not be there. To find how DSolve came up with this, I solved this by hand to see where it could have gone wrong.

It seems this first solution could be generated if DSolve evaluated or replaced $\sqrt 1$ to $\pm 1$ near the end, when finding the constant of integration. But this would be very strange thing to do. As we always take the positive root (principal root) of square root of a positive number.

Here is my solution

Solve
\begin{gather*} \boxed{y^{\prime}-2 y \left(x \sqrt{y}-1\right)=0} \end{gather*} With initial conditions $$ y \left(0\right) = 1 $$

This is a Bernoulli ODE.

\begin{align*} y' = -2 y +2 x y^{\frac{3}{2}} \tag{1} \end{align*} Dividing both sides of ODE (1) by $y^{\frac{3}{2}}$ gives \begin{align*} y'\frac{1}{y^{\frac{3}{2}}} &= -\frac{2}{\sqrt{y}} +2 x \tag{4} \end{align*} Let \begin{align*} w &=\frac{1}{\sqrt{y}} \tag{5} \end{align*} Taking derivative of equation (5) w.r.t $x$ gives \begin{align*} w' &= -\frac{1}{2 y^{\frac{3}{2}}}y' \tag{6} \end{align*} Substituting equations (5) and (6) into equation (4) gives \begin{align*} -2 w^{\prime}\left(x \right)&= -2 w \left(x \right)+2 x\\ w' &= w -x \tag{7} \end{align*} The above now is a linear ODE in $w \left(x \right)$. The integrating factor $\mu$ is \begin{align*} \mu &= {\mathrm e}^{-\int d x}\\ &= {\mathrm e}^{-x} \end{align*} Eq (7) becomes \begin{align*} \frac{\mathop{\mathrm{d}}}{ \mathop{\mathrm{d}x}}\left( \mu w\right) &= \left(\mu \right) \left(-x\right) \\ \frac{\mathop{\mathrm{d}}}{ \mathop{\mathrm{d}x}} \left(w \,{\mathrm e}^{-x}\right) &= -x {\mathrm e}^{-x} \end{align*} Integrating gives \begin{align*} w \,{\mathrm e}^{-x} &= \int{- x {\mathrm e}^{-x} \,\mathrm{d} x}\\ w \,{\mathrm e}^{-x} &= \left(x +1\right) {\mathrm e}^{-x} + c_1 \end{align*} Dividing both sides by the integrating factor $\mu={\mathrm e}^{-x}$ results in \begin{align*} w &= x +1+c_{1} {\mathrm e}^{x} \end{align*}

Replacing $w$ in the above by $\frac{1}{\sqrt{y}}$ using equation (5) gives the final solution.
$$ \frac{1}{\sqrt{y}} = x +1+c_{1} {\mathrm e}^{x} $$ Initial conditions are now used to solve for $c_{1}$.

Substituting $x=0$ and $y=1$ in the above solution gives an equation to solve for the constant of integration. This is by taking $\sqrt{1}=1$. This gives \begin{align*} c_{1} = 0 \end{align*} Substituting $c_{1}$ found above in the general solution found above gives $$ \frac{1}{\sqrt{y}} = x +1 $$ Solving for $y \left(x \right)$ from the $\frac{1}{\sqrt{y \left(x \right)}} = x +1$ gives the following solution \begin{align*} y \left(x \right) &= \frac{1}{\left(x +1\right)^{2}} \end{align*}

Now, if DSolve also used $\sqrt{1}=-1$, then $c_1=-2$ and this would result in that extra first solution DSolve gave.

Question is: Do you think this is what happened? Or do you see another path that could have caused DSolve to generate the first solution?

V 13.0 on windows 10

Answer 1

The procedure of DSolve

Where the OP solves for C[1] in the by-hand solution ($ {1}/{\sqrt{y}} = x +1+c_{1} e^{x}$), DSolve first solves for $y$ -- its process is get the general solution if possible; then specialize to the BCs. Solving for $y$ squares both sides of the equation, introducing an extraneous solution. You could say DSolve flunks 8th grade algebra. :)

@user64494's 1st & 4th steps are basically what happens (+1). Here are the steps I refer to, which I repeat to define some variables used later, although I think DSolve calls Solve instead of Reduce:

ode = y'[x] == 2*y[x]*(x*Sqrt[y[x]] - 1);
gensol = DSolve[ode, y[x], x]
(*  {{y[x] -> 1/(1 + x + E^x C[1])^2}}  *)

const = Solve[y[0] == 1 /. DSolve`DSolveToPureFunction[gensol]]
(*  {{C[1] -> -2}, {C[1] -> 0}}  *)

ivpsol = Flatten[gensol /. const, 1]
(*  {{y[x] -> 1/(1 - 2 E^x + x)^2}, {y[x] -> 1/(1 + x)^2}}  *)

Checking solutions

Checking functional solutions can be hard when we're dealing with branch cuts in equations. Here's one way it could be done in this case (and could often be expensive in other cases). At least this works in V13; I don't remember Simplify returning an inequality given an equation to simplify before, but maybe I never tried. We get the domains of each solution. Each domain should be an interval containing the boundary condition. The code below checks for an interval but it does not check that 0 belongs to it. (You'd have to parse the locations of the BCs, and the code is going to get longer. In this case DSolve returns solutions all of whose domains contain 0. Perhaps this always happens and there's no need to check it.)

ivpdomains = Simplify[
  ode /. DSolve`DSolveToPureFunction[ivpsol],
  x \[Element] Reals]
truesol = Pick[ivpsol,
  RegionDimension@ImplicitRegion[#, {x}] & /@ ivpdomains,
  1]
(*
  {x == 0, x/(1 + x) == 0 || x >= -1}
  {{y[x] -> 1/(1 + x)^2}}
*)

Alternate approach

You can also use my trick from Accessing Reduce from DSolve to use Reduce in Solve, which is more fastidious about 8th grade algebra:

Internal`InheritedBlock[{Solve}, Unprotect[Solve];
 Solve[eq_, v_, opts___] /; ! TrueQ[$in] := 
  Block[{$in = True, $res1, $res2}, 
   Solve[eq, v, Method -> Reduce, opts]];
 Protect[Solve];
 DSolve[{ode, y[0] == 1}, y[x], x] // 
  Simplify[#, x \[Element] Reals] &
 ]

Solve::useq: The answer found by Solve contains equational condition(s)...

{{y[x] -> ConditionalExpression[1/(1 + x)^2, x >= -1]}}

Answer 2

Sometimes DSolve has problems with Sqrt.

If you substitute Sqrt[y[x]]->z[x] the new ode

ode = Simplify[Derivative[1][y][x] == 2 (-1 + x Sqrt[y[x]]) y[x] /.y -> Function[x, z[x]^2] , z[x] > 0] 
(*x z[x]^2 == z[x] + Derivative[1][z][x]*)

is easily and correctly solved:

DSolve[{ode,z[0]==1}, z,x]
(*{{z -> Function[{x}, 1/(1 + x)]}}*)

Hope it helps!

Answer 3

Here is my explanation. Mathematica 13.0.0 correctly solves

DSolve[{y'[x] == 2*y[x]*(x*Sqrt[y[x]] - 1)}, y, x]

{{y -> Function[{x}, 1/(1 + x + E^x C[1])^2]}}

, but does not take into account that this is a true solution under the condition 1 + x + E^x C[1]>=0 only:

y'[x] - 2*y[x]*(x*Sqrt[y[x]] - 1) /. 
{{y -> Function[{x}, 1/(1 + x + E^x C[1])^2]}}

{-((2 (1 + E^x C[1]))/(1 + x + E^x C[1])^3) - ( 2 (-1 + x Sqrt[1/(1 + x + E^x C[1])^2]))/(1 + x + E^x C[1])^2}

, and then

FullSimplify[{-((2 (1 + E^x C[1]))/(1 + x + E^x C[1])^3) -
( 2 (-1 + x Sqrt[1/(1 + x + E^x C[1])^2]))/(1 + x + E^x C[1])^2}, 
 Assumptions ->1 +x + E^x C[1] >= 0]

{0}

FullSimplify[{-((2 (1 + E^x C[1]))/(1 + x + E^x C[1])^3) - (2 (-1 + 
x Sqrt[1/(1 + x + E^x C[1])^2]))/(1 + x + E^x C[1])^2}, 
Assumptions ->1 +x + E^x C[1] < 0]

{-((2 x (-1 + Sign[1 + x + E^x C[1]]^3))/(1 + x + E^x C[1])^3)}

As I understand it, she does

Reduce[1 == 1/(1 + x + E^x C[1])^2 /. x -> 0, C[1]]

C[1] == -2 || C[1] == 0

, not selecting only C[1]==0 (1+0+(-2)==-1 so C[1]==-2 does not produce a solution).

The solution of the problem is unique because of

D[2*y*(x*Sqrt[y] - 1), y] /. {x -> 0, y -> 1}

-2

We see a bug, however, Mathematica does not think Sqrt[1]==-1.

Answer 4

To me it looks like under 12.3, DSolve, when given initial conditions, is returning two solutions because Mathematica is treating $\sqrt{y}$ in the DE as multi-valued. This to me is the more general solution. The solution $y_1(x)=\frac{1}{(1+x)^2}$ satisfies the IVP: $$y'=2y(x\sqrt{y}-1); y(0)=1 $$ and the solution $y_2(x)=\frac{1}{(-1+2e^x-x)^2}$ satisfies: $$ y'=2y(-x\sqrt{y}-1); y(0)=1 $$ when $\sqrt{y}$ is interpreted as principal-valued.

Consider the IVP: $$ y'=2y(xy^{1/3}-1),y(0)=1 $$ In this case DSolve returns three solutions apparently corresponding to the three values of $y^{1/3}=|y|^{1/3}e^{\frac{i}{3}(\theta+2n\pi)}$. I have not checked it but I suspect in this case, the three solutions satisfy the DE:

$$ y'=2y(x|y|^{1/3}e^{\frac{i}{3}(\theta+2 n\pi)}-1),y(0)=1, n=0,1,2 $$