Why does DSolve return two solutions for my ODE?

Question

I wanted to solve the differential equation: $y’ = (1+2x)\sqrt{y}$ with $y(0) = 1$. It can be done by hand, and to check my answer I typed the following in Mathematica

DSolve[{y'[x] == (1 + 2 x)Sqrt[y[x]], y[0] == 1}, y[x] ,x]// FullSimplify

Surprisingly Mathematica gave me two solutions:

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

The second solution is what I got by hand. However, I don’t see how Mathematica got the first solution.

If I use it, I get y’(0) = -1, but according to the differential equation $y’(0)= (1+2*0)\sqrt{y(0)} = 1$.

I can see that it is a solution of $y'(x)^2 = ((1+2x)\sqrt{y})^2$, but that ODE is not equivalent to $y’ = (1+2x)\sqrt{y}$.

And when I run the following Mathematica code:

s1[x_] = 1/4 (-2 + x + x^2)^2; 
s2[x_] = 1/4 (2 + x + x^2)^2;
check1[x_] = s1'[x] - (1 + 2 x) Sqrt[s1[x]];
check2[x_] = s2'[x] - (1 + 2 x) Sqrt[s2[x]];
Plot[{check1[x]}, {x, -3, 1.5}]
Plot[{check2[x]}, {x, -3, 1.5}]

I can see that check2[x] is identically zero, but not check1[x]. What am I missing here?

Answer 1

EDIT: added handling special case of x == -1/2

Expanding on comment by@murray

eqns = {y'[x] == (1 + 2 x) Sqrt[y[x]], y[0] == 1};

sol = DSolve[eqns, y, x]

{{y -> Function[{x}, (1/4)*(4 - 4*x - 3*x^2 + 2*x^3 + x^4)]}, {y -> Function[{x}, (1/4)*(4 + 4*x + 5*x^2 + 2*x^3 + x^4)]}}

eqns /. sol // Simplify

{{(1 + 2*x)*(-2 + x + x^2 - Sqrt[(-2 + x + x^2)^2]) == 0, True}, {(1 + 2*x)*(2 + x + x^2 - Sqrt[(2 + x + x^2)^2]) == 0, True}}

The second solution is valid for

Reduce[2 + x + x^2 >= 0 || x == -1/2, x]

Element[x, Reals] || x == -(1/2)

However, Element[x, Reals] encompasses the special case of x == -1/2

eqns /. sol[[2]] // Simplify[#, Element[x, Reals]] &

{True, True}

The first solution requires that

Reduce[-2 + x + x^2 >= 0 || x == -1/2, x]

x <= -2 || x >= 1 || x == -(1/2)

So the second solution is valid for all real x and both solutions are valid for x <= -2 or x > 1 or x == -1/2

eqns /. sol // Simplify[#, x < -2 || x >= 1 || x == -1/2] &

{{True, True}, {True, True}}

Answer 2

First: Getting the right solution

Using Internal`WithLocalSettings to set the Method -> Reduce option on Solve, and adding the assumption x ∈ Reals, we get the correct answer:

sys = {y'[x] == (1 + 2 x) Sqrt[y[x]], y[0] == 1};

With[{opts = Options[Solve]},
 Internal`WithLocalSettings[
  SetOptions[Solve, Method -> Reduce],
  s = Assuming[x ∈ Reals,
    y -> Function @@ {x, y[x] /. #} & /@    (* turn the expression into a Function *)
     Simplify@DSolve[sys, y[x], x]          (* Simplify the solution *)
    ],
  SetOptions[Solve, opts]
  ]]
sys /. s // Simplify[#, x ∈ Reals] &

> Solve::useq: The answer found by Solve contains equational condition(s) {0==1/2 (x+x^2+Reduce`ReduceVar[2]-Sqrt[Plus[<<3>>]^2])}. A likely reason for this is that the solution set depends on branch cuts of Wolfram Language functions. >>

(*
  {y -> Function[x, 1/4 (2 + x + x^2)^2]}
  {True, True}
*)

Without the assumption x ∈ Reals , we also get the right answer but in a less convenient form, with an unsimplified ConditionalExpression alluded to in the warning message:

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

Note: I prefer the Function form of the solution s because the derivative(s) in a differential equation sys are automatically computed in sys /. s. This seems convenient to me.

Second: What DSolve is doing

I think DSolve is doing what one normally does. It solves for the general solution and then solves the initial condition for the constant of integration. The problem is that it carries the general solution a little too far, which Chip Hurst remarked on.

s = DSolve[{y'[x] == (1 + 2 x) Sqrt[y[x]]}, y, x]
(*  {{y -> Function[{x}, 1/4 (x^2 + 2 x^3 + x^4 + 2 x C[1] + 2 x^2 C[1] + C[1]^2)]}}  *)

Solve[y[0] == 1 /. s, C[1]]
(*  {{C[1] -> -2}, {C[1] -> 2}}  *)

Note the C[1]^2 in s. One problem at this stage is that one cannot determine which value of C[1] is correct from the initial condition alone. You have to go back and check the ODE. So Method -> Reduce is ineffective at this step:

Solve[y[0] == 1 /. s, C[1], Method -> Reduce]
(*  {{C[1] -> -2}, {C[1] -> 2}}  *)

Apparently DSolve does not go back and check. I suppose there are reasons for this "laziness." It might take an unknown and long time to verify a potentially complicated solution, for instance. But that could be handled with TimeConstrained and a warning message on time-out.

One can implement this strategy oneself, with or without TimeConstrained, as follows:

s = DSolve[{y'[x] == (1 + 2 x) Sqrt[y[x]], y[0] == 1}, y, x];
Pick[s, TrueQ /@ And @@@ Simplify[sys /. s, x ∈ Reals]]
(* {{y -> Function[{x}, 1/4 (4 + 4 x + 5 x^2 + 2 x^3 + x^4)]}}  *)

Answer 3

If we run WolframAlpha["y'[x] == (1 + 2 x)Sqrt[y[x]], y[0] == 1"], we can see the step by step solution, which exposes the issue at hand.

WolframAlpha["y'[x] == (1 + 2 x)Sqrt[y[x]], y[0] == 1"]
(* click the Step-by-step solution *)

Now, I cut off the image here on purpose.

Notice to solve for y[x], we had to square both sides. This could give extraneous solutions, dependent on x, so Solve just gives the generic solution. This seems acceptable IMHO since Solve tries to return a list of rules. Reduce on the other hand does a better job:

Reduce[2 Sqrt[y[x]] == x^2 + x + C[1], y[x]]

During evaluation of In[58]:= Reduce::useq: The answer found by Reduce contains unsolved equation(s) {0==1/2 (x+x^2+C[1]-Sqrt[Plus[<<3>>]^2])}. A likely reason for this is that the solution set depends on branch cuts of Wolfram Language functions. >>

0 == 1/2 (x + x^2 + C[1] - Sqrt[(x + x^2 + C[1])^2]) && 
 y[x] == 1/4 (x + x^2 + C[1])^2

Here's a post asking how to force DSolve to use Reduce instead of Solve. I'm not sure how well the solution works...

If you want to see a simpler example of why Sqrt's can give extraneous solutions, consider this example:

WolframAlpha["Sqrt[x] == -1"]
(* click the Step-by-step solution *)