# does DSolve solution violate Picard-Lindelof theorem (existence and uniqueness of solution of first order ode) in this problem?

I searched this site looking for similar question but could not find one. Feel free to close this if it is duplicate.

Given the Bernoulli ode

\begin{align} y^{\prime} & =y+x\sqrt{y}\tag{1}\\ y\left( 0\right) & =4\tag{2} \end{align}

I have two questions about Mathematica solution to the above ode.

The first is: I solved this by hand and obtained the general solution (before applying the initial conditions) as (I can post the hand solution if needed)

$$\begin{equation} \sqrt{y}=-2-x+4e^{\frac{x}{2}}\tag{3} \end{equation}$$

Also Maple gives same solution as above: But Mathematica gives the general solution as

ClearAll[x, y];
ode = y'[x] - y[x] == x*y[x]^(1/2);
ic = y == 4;
solNoIc = DSolve[ode, y[x], x] Mathematica automatically squared both side of the solution. But this causes a problem later on when finding the initial conditions, since it introduces extra root.

Question: Should Mathematica have left the general solution with the sqrt in the general solution as in (3)?

The second question on this solution: Adding now initial conditions, Mathematica gives two solutions. This is as side effect of squaring both side of (3). It gives

solWithIc = DSolve[{ode, ic}, y[x], x] But these two solutions are the same only at $$x=0$$. Any interval around $$x=0$$ these two solutions are not the same. But this violates existence and uniqueness theory of first order ode (Picard-Lindelof theorem). Writing (1) as $$y^{\prime}=f\left( x,y\right) =$$ $$y+x\sqrt{y}$$, this shows $$f$$ is continuous for all $$x,y$$. And $$\frac{\partial f}{\partial y}=1+\frac{x}{2\sqrt{y}}$$. This shows that $$\frac{\partial f}{\partial y}$$ is continuous for all $$x,y$$ except at $$y=0$$.

So the solution can not cross $$y=0$$. So there exists solution which is unique in some rectangle either in upper half of the plane or lower half of the plane. But initial conditions says $$y=4$$ at $$x=0$$.

Since initial conditions is in upper half of plane. This means there exists a unique solution in some rectangle around the initial conditions point and exists only in upper half of plane.

But Mathematica gives 2 solutions. These are the same only at $$x=0$$. So there is no interval (no matter how small) around $$x=0$$ where the solution of the ode is unique. This violates Picard-Lindelof. The ode should have only one solution.

## Verify Mathematic's solutions

Here is attempt to verify Mathematica's own solutions. It verifies the first solution but not the second:

ClearAll[x,y];
ode=y'[x]-y[x]==x*y[x]^(1/2);
ic=y==4;
solWithIc=DSolve[{ode,ic},y,x]
{ode,ic}/.solWithIc//FullSimplify Reduce[x (2 - 4 E^(x/2) + x + Sqrt[(2 - 4 E^(x/2) + x)^2]) ==   0, Reals] Reduce[ x (2 + x + Sqrt[(2 + x)^2]) == 0, Reals] So Mathematica says its second solution is only valid at initial conditions $$x=0$$ and for $$x\leq -2$$. But it should be an interval centered at $$x=0$$.

The question is: Is Mathematica solution mathematically correct? Does Mathematica solution for the IVP violate Picard-Lindelof theorem?

• The equation $y' = |y|^{1/2}$ is a standard textbook example where Picard-Lindelof does not apply, see here for example. Jul 16, 2022 at 8:09
• Suppose we interpret the square root in the conventional way as mapping $[0,\infty)$ to $[0,\infty)$. Given $y(0) = 4$, then yes in a small $x$-intervall around $x=0$ there is a unique solution. That unique solution is $y(x) = (-2-x+4e^{x/2})^2$. The function $y(x) = (x+2)^2$ is not a solution. If we interpret the square root in non-conventional way as mapping $[0,\infty)$ to $(-\infty,0]$ then $y(x)=(x+2)^2$ is the unique solution. Why Mathematica includes this one, I do not know. Note: Both functions solve the equation $(y'-y)^2 = x^2 y$ obtained by squaring the original equation. Jul 16, 2022 at 9:26
• Mathematica cannot violate a theorem, unless the theorem is not a theorem; however, it can be wrong. I'm not sure how good Mma is at solving/verifying functional equations. Solve solves equations algebraically. I don't know how easy it is to verify that a solution is valid in a neighborhood and not just at a point. Further, I'm not sure DSolve checks the results of Solve. The solution that''s "valid" only at a point is not a valid solution; otherwise, one could accept any expression whose value at x=0 is y=4. — BTW, Picard-Lindelöf does apply here since the IC is at y=4, not y=0. Jul 16, 2022 at 14:17
• The solutions are perhaps valid in the 17/18th cent. Eulerische sense, that the same equation satisfies the ODE over some interval and at the IC. (That's neither Euler's view per se nor the modern view. It's my irreverent view of pre-modern analysis, in which Sqrt[] was multivalued and convergence & continuity weren't always a concern, leastways not in the modern sense.) Anyway, the spurious solution does not violate P-L because it does not satisfy the ode in an interval containing x=0. Only one of the solutions satisfies the ode in such an interval. Jul 16, 2022 at 14:32
• Thanks to all the great answers. I wish I can accept both of them. Both are useful to learn from. I think Mathematica should return only one solution and not two. The second solution results, because it squared both sides of the general solution to obtain explicit solution. Jul 17, 2022 at 3:08

This is in response to your last question: You are right that one solution returned by Mathematica is wrong, in the sense that it is not a solution to the initial value problem for $$x$$ close to zero.

Simpler problem. Let me first consider the following initial value problem:

DSolve[{y'[x] == + Sqrt[y[x]], y == 1}, y[x], x] // Simplify
DSolve[{y'[x] == - Sqrt[y[x]], y == 1}, y[x], x] // Simplify


The two differential equations differ only in the sign in front of the square root. In this case, the solutions returned by Mathematica V12.3 are the correct unique solutions for $$x$$ close to $$0$$ at least: Note that these are not actually solutions for all $$x$$: Plugging $$y(x) = \frac{1}{4}(2+x)^2$$ into $$y' = \sqrt{y}$$ gives $$\frac{1}{2}(2+x) = \frac{1}{2}|2+x|$$ which only holds for $$x \geq -2$$. But let me ignore this here and just accept that for $$x$$ close to $$0$$ at least, Mathematica gave the unique correct solution.

Your problem (a Bernoulli differential equation). Let me now take your initial value problem, as well as the one with the other sign in front of the square root:

DSolve[{y'[x] == y[x] + x*Sqrt[y[x]], y == 4}, y[x], x]
DSolve[{y'[x] == y[x] - x*Sqrt[y[x]], y == 4}, y[x], x]


In this case Mathematica V12.3 returns, without any warning, two solutions for the first and two solutions for the second equation, and actually it returns exactly the same two solutions in both cases: So effectively the sign in the differential equation is ignored. You are right that in each case there is actually only one solution to the initial value problem for $$x$$ close to $$0$$, not two. Mathematica returns too many solutions. I will leave the calculation for the reader, and show a plot instead.

Slope field plot. Here is the slope field for the equation $$y' = y + x\sqrt{y}$$ The correct solution $$(-2-x+4e^{x/2})^2$$ is in green, the wrong solution $$(2+x)^2$$ that is also returned by Mathematica is in orange, and is clearly not a solution: The plot was generated using this code:

With[{a=0.5},Show[{
VectorPlot[{1,y+x*Sqrt[y]},{x,-a,a},{y,2,6},VectorColorFunction->None,VectorStyle->{Blue,Opacity[0.4]},FrameLabel->{"x","y"}],
Plot[(-2+4*Exp[x/2]-x)^2,{x,-a,a},PlotStyle->Green],
Plot[(2+x)^2,{x,-a,a},PlotStyle->Orange]
}]]


Should Mathematica have left the general solution with the sqrt in the general solution as in (3)?

This was connected to the spurious solution $$y=(x+2)^2$$ to the IVP. It's not clear that those have to be connected. If they can be separated, then I'd rather get an explicit function from DSolve than an implicit solution.

Algebra is generally easier with rational function than algebraic ones, and I surmise that DSolve rationalizes the ODE. As I alluded to in one of my comments, at one point people thought the square root should be multivalued, before the functional approach took over analysis. I like the (algebraic) geometry of ODEs so much that the spurious solutions from DSolve are interesting to me. To someone solving a practical problem, spurious solutions are distracting and extra work and a darn nuisance. Personally, I think it would be nice if DSolve came with a way to verify solutions. Oh wait, there's this from it doc page:

With a pure function output, eqn/.{{u->f},[Ellipsis]} can be used to verify the solution.  »

There's even a section in the tutorial Working with DSolve: A User's Guide with the title "Verification of the Solution". One of the methods seems to be conclusive about the spurious solution (which yields the second pair in the output below):

{ode, ic} /. Equal -> Subtract /.
solWithIc /. {x -> RandomReal[{-1, 1}]}
(*  {{0., 0}, {1.73578, 0}}  <-- Zero means a solution *)


Miraculously, none of the examples in the tutorial show how to reject a spurious solution. I think it would be nice if they added a PossibleDSolveSolutionQ that would work like PossibleZeroQ, doing a short symbolic check followed by a numeric check. It could even be an option to DSolve. It would also be nice to check the solution over an interval containing the initial condition. The OP shows a way to do that (there are others), but that sort of check would have eliminated the spurious solution.

To summarize the answer to the first question: I do not think that returning an implicit solution to the general ODE is the way to fix the spurious solution to the IVP. I would rather have a function returned.

For fun, here's a picture of what's going on. The ode defines a manifold in {x, y[x], y'[x]} space (the greenish surface below). The conjugate ODE (ode /. Sqrt[r_] :> - Sqrt[r]) is the reddish surface. Together they form the surface of the rationalized ODE:

(-y[x] + y'[x])^2 == -x^2 y[x]


A true solution to ode defines a curve on the green manifold. Some solutions may be extended to a solution of the conjugate ODE. The true solution to the OP's IVP is shown in thick blue, and the spurious solution is shown in thick golden; the white dot represents the initial condition. The spurious solution is a valid solution when it's on the green manifold; strictly speaking, it's not a solution in a neighborhood of the initial condition. Theoretically, one can resolve the singularity (where the surface intersects itself), by differentiating the ODE, but DSolve fails to solve it, whether ode or its rationalized form:

DSolve[{D[ode, x], ic, ode /. x -> 0 /. First@Solve[ic]}, y[x], x]

oderat = (-y[x] + Derivative[y][x])^2 == -x^2 y[x];
DSolve[{D[oderat, x], ic, ode /. x -> 0 /. First@Solve[ic]}, y[x], x]


Code for plot

solWithIc2 = MapThread[
DSolve[{ode, y[#] == #2}, y, x] &,
Transpose@
MeshCoordinates@
DiscretizeRegion[Rectangle[{-5, 1}, {3, 6}],
MaxCellMeasure -> 3/4]
] // Apply@Join;

Show[
ParametricPlot3D[{x, u^2, u^2 - u x},
{x, -5, 3}, {u, -2.5, 2.5},
MeshFunctions -> {(*#1(#2-#3)*)#5 &}, Mesh -> {{0}},
Directive[Opacity[0.5], ColorData]},
PlotPoints -> {15, 40}],
ParametricPlot3D[{x, y[x], y'[x]} /. solWithIc2 // Evaluate,
{x, -5, 3},
PlotStyle ->
Directive[AbsoluteThickness[0.4], Opacity[0.7],
Darker@ColorData]],
ParametricPlot3D[{x, y[x], y'[x]} /. solWithIc // Evaluate,
{x, -5, 3}],
Graphics3D[{
Sphere[{0, 4,
y' /. First@Solve[ode /. {y[x] -> 4, x -> 0}, y']}, 0.15]}],
PlotRange -> {{-5, 3}, {-2, 5}, {-5, 5}},
BoxRatios -> Automatic, AxesLabel -> {x, y[x], y'[x]}
]

• I do not think that returning an implicit solution to the general ODE is the way to fix the spurious solution to the IVP If there was no initial conditions present, then yes, I agree. But when there are initial conditions like this problem, then the explicit general solution now results in two constants of integrations when substituting the initial conditions. This gives rise to the two solution given by Mathematica. If Mathematica had kept the general solution implicit, then only one constant of integration will result. This is what I meant. Maple also gives one solution to this ode. Jul 17, 2022 at 2:59
• @Nasser I don't see why the implicit solution has to be returned as the general solution in order to solve the IVP correctly. Solve the IVP correctly, then solve for y[x] and return that. Whichever problem is being solved, it could certainly solve for y[x] at the end, right? Jul 17, 2022 at 3:25
• When I solved this by hand, I see that two constant of integration will result, due to forcing the general solution to be explicit beforehand (side effect of squaring both sides). But when keeping the general solution implicit, only one constant results in this case. If you like I can post the step by step solution to this ode. May be I made an error but I do not think so, since it agrees with Maple's solution. If there was no I.C., given then it is safe to square both sides and return an explicit general solution then (which will have unknown constant of integration in it)... Jul 17, 2022 at 3:34
• ... at the very end, i.e. after solving for constant of integration, it is then save to solve for explicit $y(x)$. But not before finding the constant of integration. I find this way avoids introducing spurious solutions. Jul 17, 2022 at 3:37
• @Nasser DSolve could return this: DSolveDSolveBoundaryValueProblem[DSolve, Unevaluated@Solve[Sqrt[y[x]] == 1/2 (-4 - 2 x) + E^(x/2) C, y[x]], {y == 4}, {C}, {x}] Jul 17, 2022 at 3:38

I know I'm a little bit late...

Mathematica often shows problems concerning Sqrt- expressions inside NDSolve, because it tries to rationalize the ode.

Substituting Sqrt[y[x]] -> z[x] as workaround avoids this problem and gives the unique solution of the problem:

odez = ode /. y -> (z[#]^2 &) // Simplify[#, z[x] > 0] &
(*x + z[x] == 2 Derivative[z][x]*)

solz = DSolve[odez, z[x] , x][]
(*z[x] -> 1/2 (-4 - 2 x) + E^(x/2) C*)


Resubstitution gives y[x]

solz /. z[x] -> Sqrt[y[x]]
(*{Sqrt[y[x]] -> 1/2 (-4 - 2 x) + E^(x/2) C}*)


I tried something like this here just now, but in this case, one has to add initial conditions up to the second order.

DSolve[{ode, Reduce[{D[ode, x], ode, ic} /. x -> 0]}, y[x], x]

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


where the initial conditions are calculated to be

Reduce[{D[ode, x], ode, ic} /. x -> 0]

y'' == 6 && y' == 4 && y == 4


Actually, I'm surprised DSolve will use them above the order of ode`. In any case, it has worked twice now. It almost seems like a method.