We all know that for real $x,y$ that $\sqrt{x y} = \sqrt{x} \sqrt{y}$ only when $x,y$ are nonnegative. But Mathematica DSolve automatically does this when solving an ode. The question is why?

Here is an example

DSolve[y'[x] == Sqrt[x y[x]], y[x], x]

(if you solve it by hand, how would you do this?)

Mathematica gives

Mathematica graphics

How did it do it? It clearly rewrote Sqrt[x y[x]] as Sqrt[x]* Sqrt[y[x]] and made it separable and just integrated. This can be seen by calling WolframAlpha and looking at the step-by-step


enter image description here

Lets ask Reduce what it thinks about all of this

Reduce[Sqrt[x*y] == Sqrt[x]*Sqrt[y], Reals]

Mathematica graphics

Even if do not use Reals above

Reduce[Sqrt[x*y] == Sqrt[x]*Sqrt[y]]

Mathematica graphics

The point is, we can't just replace Sqrt[x*y] by Sqrt[x]*Sqrt[y]] without assumptions. Right?

This becomes more clear when adding assumptions to DSolve itself. Now it gives different answers depending if we tell it $x$ is positive or negative

DSolve[y'[x] == Sqrt[x y[x]], y[x], x, Assumptions -> Element[{x, y}, Reals]]

Mathematica graphics

DSolve[y'[x] == Sqrt[x y[x]], y[x], x, Assumptions -> {x > 0, y > 0}]

Mathematica graphics

DSolve[y'[x] == Sqrt[x y[x]], y[x], x, Assumptions -> {x < 0, y > 0}]

Mathematica graphics

This shows that the answer depends on the sign.

I am just trying to understand if Mathematica's answer is valid mathematically wise and why it breaks the sqrt like this without using assumptions.

One possibility I was thinking about why the above can be done when solving the ode, is that the solution has an arbitrary constant of integration and that the constant of integration can be complex itself. But I am not sure that this will give one permission to break the sqrt as done here at the start. But that is what DSolve does.

When I solve this ode myself by hand, I will at least add assumptions that $x>0,y>$ before breaking up the sqrt. So may be if DSolve generated the answer with this assumption included, then that would be better.

V 13.2 on windows 10

  • $\begingroup$ This is not an MMA answer, so I am leaving it here as a comment: (if you solve it by hand, how would you do this?). Well, the first observation is that $y(x)=0$ is a solution $\forall x$. So, let's find $y(x) \neq 0$. Divide LHS and RHS by $\sqrt{y}$ to bring it in the form: $(2 \sqrt{y})^{\prime}=\sqrt{x}$ which can be readily integrated to give $y = \left(\frac{1}{3} x^{3/2} + \text{constant} \right)^2$. $\endgroup$
    – bmf
    Commented Jan 26, 2023 at 5:20
  • 1
    $\begingroup$ Divide LHS and RHS by $\sqrt y$ @bmf thanks. But this is basically what Mathematica did. made it separable. But one can't do $\frac{\sqrt{x y}}{\sqrt y}$ and simplify it to $\sqrt x$ unless both are $x,y$ positive. And this is the main question I have. Screen shot !Mathematica graphics I think if I do this in the exam, the teacher might take off few points unless I explicitly say assuming $x,y>0$ before. $\endgroup$
    – Nasser
    Commented Jan 26, 2023 at 7:07
  • 2
    $\begingroup$ @Nasser Mathematica solves incorrectly the squared ode y'[x]^2 == x y[x] $\endgroup$ Commented Jan 26, 2023 at 7:14
  • $\begingroup$ @Nasser sorry for the delayed response, but I have a project that picked up and has lots of computations. We have established that we are looking for solution $y(x) \neq 0$. What I meant in my comment is that you can do this. Without assumptions on $x$ and $y$, right? $\endgroup$
    – bmf
    Commented Jan 26, 2023 at 10:52
  • $\begingroup$ @bmf but PowerExpand assumes all quantities are positive. From help, The transformations made by PowerExpand are correct in general only if c is an integer or a and b are positive real numbers. so it is assuming $x,y$ are positive. We are back to square one. One can't simplify as DSolve did and break the sqrt without this assumption. $\endgroup$
    – Nasser
    Commented Jan 26, 2023 at 11:36

1 Answer 1

eqn = y'[x] == Sqrt[x y[x]];

sol = DSolve[y'[x] == Sqrt[x y[x]], y, x][[1]]

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

Verifying the solution will provide the required constraints.

eqn2 = eqn /. sol // Simplify

(* Sqrt[x (2 x^(3/2) + 3 C[1])^2] == 2 x^2 + 3 Sqrt[x] C[1] *)

Reduce[eqn2, x, Reals] // ToRadicals

(* (C[1] < 0 && (x == 0 || x >= (-(3/2))^(2/3) C[1]^(2/3))) || 
  (C[1] >= 0 && x >= 0) *)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.