Why does DSolve fail in this first order differential equation?

I am trying to solve the following first order differential equation: $$g'(R)=-2 \sqrt{R^2-g(R)}+2 R$$

By direct substitution it can be verified that an obvious solution is $$g(R)=R^2$$

ode=2R-2Sqrt[R^2-g[R]]-Derivative[1][g][R]
g[R_]:=R^2
ode==0


Giving True as expected.

However, using DSolve (Mathematica 11.2 for Windows 10):

Clear[g]
DSolve[ode==0, g[R], R]//FullSimplify//Expand


we find

{{g[R]->-1+E^(2 C[1])-2 R-2 Sqrt[-E^(2 C[1]) (1+R)^2]},
{g[R]->-1+E^(2 C[1])-2 R+2 Sqrt[-E^(2 C[1]) (1+R)^2]}}


As seen, Mathematica completely missed the correct solution!

Edit 1: For more details on the mathematics aspects of the solution (and the two branches of the solution), see my post at the Mathematics stack.

Edit 2: Adding initial conditions like $$g(0)=0$$ or $$g(1)=1$$, as suggested in the comments, unfortunately doesn't help.

Edit 3: Just for future reference I add the way I found to solve the ode and recover the $$R^2$$ solution. To do that define a new function $$y(R)$$ that should be zero when $$g(R)=R^2$$, in the following manner

g[R_] := R^2 (1 - y[R]^2)
ode1 = ode // FullSimplify // PowerExpand // FullSimplify;
DSolve[ode1 == 0, y[R], R] // Expand


which gives

{{y[R] -> 0}, {y[R] -> 1 + C[1]/R}}


The solution $$y[R] \rightarrow 0$$ implies $$g(R)=R^2$$ as expected.

After massaging the ODE as shown above, I can get Mathematica to give the correct solution, however I still don't understand why it failed in the first place.

So, my question is, why did DSolve fail in this simple example?

• If you integrating equation ,then where is integrating constant C[1] ? Feb 5, 2020 at 17:27
• The fact that you're getting two solutions here makes me wonder whether Mathematica is increasing the degree of the ODE (i.e., internally rewriting it in terms of ${g'}^2$), which then causes spurious results. Feb 5, 2020 at 19:44
• Also, how long does MM 11 take to find results for this equation? I'm running MM 12 on a Mac, and when I run your code it doesn't provide any results before I give up and abort the calculation (several minutes). Feb 5, 2020 at 19:45
• This is an interesting question. I've looked at numerical solutions given by NDSolve and I can't get it to give me the R^2 solution. Note that the DSolve doesn't really give two solutions: it gives two expressions, which can match different initial conditions. Feb 5, 2020 at 22:23
• @mikado: I think this is because the function $f(g,R) = -2 \sqrt{R^2 - g^2} + 2R$ is not Lipschitz continuous in $g$, which means (in particular) that we shouldn't expect there to be a unique solution for a given set of initial conditions (i.e., the Picard-Lindelöf theorem doesn't apply.) Numerical algorithms such as NDSolve kind of implicitly assume that the solution is unique. Feb 6, 2020 at 14:45

I can no longer test V11.2, but in V13, DSolve has been improved. It still needs some help. The option IncludeSingularSolutions yields more solutions, but we need to rationalize the ode to get all of them. Of course, the solutions to the rationalized ode need to be checked on the original ode.

The constant of integration shows up in several places as Exp[2 C[1]]. One can apply constSimplify or just do the substitution C[1] -> Log[C[1]]/2 directly to get a nicer looking result. We check the solutions and find conditions for which they are valid.

dsol = DSolve[ode, g[R], R, IncludeSingularSolutions -> True] /.
C[1] -> Log[C[1]]/2 // Simplify
(*
{{g[R] -> -1 - 2 R + C[1] - 2 Sqrt[-(1 + R)^2 C[1]]},
{g[R] -> -1 - 2 R + C[1] + 2 Sqrt[-(1 + R)^2 C[1]]},
{g[R] -> 0},
{g[R] -> -1 - 2 R}}
*)

cond = ode /. DSolveDSolveToPureFunction@dsol // Map[Reduce[#, Reals] &]
(*
{R > -1 && C[1] <= 0,
(R < -1 && C[1] <= -1 - 2 R - R^2) ||
(R > -1 && -1 - 2 R - R^2 <= C[1] <= 0),
R >= 0,
R >= -1}
*)


Likewise, for the rationalization oderat of the original ode, we solve the differential equation and find conditions for which they are valid.

oderat = Eliminate[{ode /. Sqrt[R^2 - g[R]] -> dummy, dummy^2 == R^2 - g[R]},
dummy]

(*  (4 R - g'[R]) g'[R] == 4 g[R]  *)

dsolrat = DSolve[oderat, g[R], R, IncludeSingularSolutions -> True]

(*  {{g[R] -> R C[1] - C[1]^2/4}, {g[R] -> R^2}}  *)

condrat =
ode /. DSolveDSolveToPureFunction@dsolrat // Map[Reduce[#, Reals] &]

(*  {R >= C[1]/2, True}  *)


We can then assemble the final solutions, using ConditionExpression to restrict each solution (optional). All seem to be valid for some values of C[1] and interval for R.

allsols = {g[R] -> #} & /@ MapThread[ConditionalExpression,
{g[R] /. Join[dsol, dsolrat], Join[cond, condrat]}]
(*
{{g[R] -> ConditionalExpression[
-1 - 2*R + C[1] - 2*Sqrt[(-(1 + R)^2)*C[1]],
R > -1 && C[1] <= 0]},
{g[R] -> ConditionalExpression[
-1 - 2*R + C[1] + 2*Sqrt[(-(1 + R)^2)*C[1]],
(R > -1 && 2*R + R^2 + C[1] >= -1 && C[1] <= 0) ||
(R < -1 && 2*R + R^2 + C[1] <= -1)]},
{g[R] -> ConditionalExpression[ 0, R >= 0]},
{g[R] -> ConditionalExpression[ -1 - 2*R, R >= -1]},
{g[R] -> ConditionalExpression[ R*C[1] - C[1]^2/4,
R >= C[1]/2]},
{g[R] -> R^2}}
*)


Re-check solutions, assuming the conditions are satisfied for each one:

Map[
With[{conditions =
Cases[#, ConditionalExpression[_, assum_] :> assum, Infinity]},
Simplify[Reduce[Flatten@{ode /. #, conditions}], conditions]] &,
DSolveDSolveToPureFunction@allsols]

(*  {True, True, True, True, True, True}  *)

• Thanks for following up on this, even after so long :-) Jul 22, 2022 at 7:39

Usually a solution to an ODE needs as many constants as their degree. In this case, the solution

$$g(R) = R^2$$

doesn't requires such a constant. Analyzing the phase plane for this ODE we can observe that the initial conditions for a real solution should obey

$$R_0 \ge 0,\ \ g(R_0) \lt R_0^2$$

and for $$g(R)= R^2$$ we have a particular solution. Follows a plot showing the stream plot (phase plane) for this ODE, in red the particular solution and in green two solutions: for $$g(0) = 0$$ and for $$g(2+\epsilon^2) = 4$$

gr1 = StreamPlot[{1, -2 (Sqrt[x^2 - y] - x)}, {x, 0, 10}, {y, -1, 9}];
gr2 = ContourPlot[y == x^2, {x, 0, 10}, {y, -1, 9}, ContourStyle -> Red];
sol3 = NDSolve[{g'[r] == -2 Sqrt[-g[r] + r^2] + 2 r, g[0] == 0}, g, {r, 0, 10}][[1]];
gr3 = Plot[Evaluate[g[r] /. sol3], {r, 0, 10}, PlotStyle -> Green];
sol4 = NDSolve[{g'[r] == -2 Sqrt[-g[r] + r^2] + 2 r, g[2.01] == 4}, g, {r, 0, 10}][[1]];
gr4 = Plot[Evaluate[g[r] /. sol4], {r, 0, 10}, PlotStyle -> Green];
Show[gr1, gr4, gr3, gr2]


Not an answer, but might help getting help...

I evaluated the integral using the package NCAlgebra, which considers g as a noncommutative variable. You can see (it you click on the image...) that it transforms the DSolve into a Solve with two integration constants.

<< NC;
<< NCAlgebra;
DSolve[2 R - 2 Sqrt[R^2 - g[R]] - Derivative[1][g][R] == 0, g[R], R]
`

Maybe some gurus can infer what is happening from there (probably MMA transforms the ODE).

• Quite interesting indeed, but I'm not so sure how to proceed from there! thanks btw. Feb 6, 2020 at 15:24