12
$\begingroup$

[Note that in the cited duplicate, DSolve not finding solution I expected, the general solution returned by DSolve is missing a solution for quite different reasons than this one. In the duplicate question, the missing solution may be obtained as the limiting function (or closure) of the solution set. In this question, the missing solution is an isolated point and not the limiting value of the general solution returned by DSolve. The reason DSolve misses a solution here is different than the reasons given in the answers to the duplicate. Further the obstruction to obtaining the missing solution is deeper and more difficult to bridge. -Michael E2]


I am trying to find the general solution to the differential equation

$y'=y/x-2/x^2+(2\sqrt{1-xy})/x^2$.

DSolve[y'[x] == y[x]/x - 2/x^2 + (2*Sqrt[1 - x y[x]])/x^2, y[x], x, 
GeneratedParameters -> (Subscript[c, #] &)]

I used DSolve, it gives the linear functions:

$y=-\frac{4 \left(c_1+x\right)}{c_1^2}$

I know, that $y=1/x$ is also a solution, which is not found by DSolve. I understand that there are limits of DSolve, but there was no indication in the output that these limits are reached. In this case, I can work out the solution, but in other cases, how can I be sure that DSolve returned all possible solutions? Is there an option that I can give to DSolve to explicitly request a note that maybe not all solutions are listed?

$\endgroup$
4
  • $\begingroup$ Please add your actual Mathematica code to the Question. $\endgroup$
    – bbgodfrey
    Commented Jun 17, 2015 at 15:31
  • 1
    $\begingroup$ From the documentation: DSolve can find general solutions for linear and weakly nonlinear partial differential equations. Truly nonlinear partial differential equations usually admit no general solutions. Since your equation is not linear there should be no expectation of finding all possible solutions. $\endgroup$
    – Bob Hanlon
    Commented Jun 17, 2015 at 16:14
  • 3
    $\begingroup$ The equation is weakly non-linear, though (i.e., the equation is linear in its highest-order derivatives), and Mathematica does find a general solution. It just doesn't find the singular solution. (Note that $y = 1/x$ is the envelope of the family of solutions.) $\endgroup$ Commented Jun 17, 2015 at 18:07
  • $\begingroup$ Guys i cant find what is shortcut of OP's at differential equations ? $\endgroup$ Commented Nov 30, 2018 at 23:32

2 Answers 2

31
$\begingroup$

$ \def\figA#1{\def\nextfig##1{\figB##1}\def#1{1}#1} \def\figB#1{\def\nextfig##1{\figC##1}\def#1{2}#1} \def\figC#1{\def\nextfig##1{\figD##1}\def#1{3}#1} \def\figD#1{\def\nextfig##1{\figE##1}\def#1{4}#1} \def\figE#1{\def\nextfig##1{\figF##1}\def#1{5}#1} \def\figF#1{\def\nextfig##1{\figG##1}\def#1{6}#1} \def\figG#1{\def\nextfig##1{\figH##1}\def#1{7}#1} \def\figH#1{\def\nextfig##1{\figI##1}\def#1{8}#1} \def\figI#1{\def\nextfig##1{\figJ##1}\def#1{9}#1} \def\figJ#1{\def\nextfig##1{\figK##1}\def#1{10}#1} \def\figK#1{\def\nextfig##1{\figL##1}\def#1{11}#1} \def\nextfig#1{\figA#1} $

Introduction

In this answer we will develop a geometric picture of the OP's differential equation. It will illustrate the behavior of DSolve on such equations. We will show that the singular solution arises as the discriminant of the ODE viewed as an equation in three variables, and that the discriminant curve necessarily is the same as the envelope of the general solution shown in the answer by Michael Seifert. The approach will reveal two more solutions, one of which might be considered spurious according to one's point of view. We will see that the general solution is valid only over a restricted domain. Finally, our analysis will suggest reasons why DSolve does not, and probably does not attempt to, find the singular solution.

We will give a general development of the ideas and return from time to time to the OP's ODE as a particular example.

Set-up

A solution $y=y(x)$ of a first-order ODE $F(x,y(x),y'(x))=0$ parametrizes a curve on the surface $$F(x,y,p)=0\,,$$ where we think of the coordinate $p$ as $p = dy/dx$. The projection $(x,y,p)\mapsto(x,y)$ maps the curve on the surface back to the solution curve in the $xy$ plane. In what follows, we will assume $F$ is sufficiently smooth and, for the code, that $F$ is a polynomial (or rational). That is not required by mathematics, but if the ODE is a polynomial at least in $p$, then there is a simple function in Mathematica that can yield the answer.

In the OP's case, we need to rationalize the ODE to eliminate the radical. There are many ways to do this (by hand is probably faster in this case than trying to figure out how to do it in Mathematica), such as GroebnerBasis.

Clear[x, y, p, F];
ode = y'[x] == y[x]/x - 2/x^2 + (2*Sqrt[1 - x y[x]])/x^2;
odeToFn = {y[x] -> y, y'[x] -> p, Equal -> Subtract};
fn = ode /. odeToFn;
F[x_, y_, p_] = Module[{u},
  With[{rads = DeleteDuplicates@
      Cases[fn, Power[a_, b_Rational] :> u[a, b], Infinity]},
   First@GroebnerBasis[
     Flatten@{fn /. Power[a_, b_Rational] :> u[a, b],
       Map[
        Numerator@ Together[  (* gets rid of denominators in neg. powers *)
          #^Denominator[Last@#] - First[#]^Numerator[Last@#]] &,
        rads]},
     rads]
   ]]
F[x_, y_, p_] = Numerator@ Together@ F[x, y, p]
(*  p^2 + (4 p)/x^2 - (2 p y)/x + y^2/x^2  *)
(*  4 p + p^2 x^2 - 2 p x y + y^2  *)

Clearing the denominators with Numerator@Together@F[x, y, p] can help with numerics.

Rationalizing usually introduces extra solutions arising from conjugate equations. We will return to this issue after we have developed a clearer picture of the ODE.

Terminology

For convenience we will call the solution returned by DSolve the general solution. A particular solution of the general solution we will call a regular solution. The differential equation ode we will call the OP's ODE. The rationalized equation $F(x, y(x), y'(x))=0$ we will call the differential equation $F = 0$; the surface $F(x,y,p)=0$ we will call the surface $F=0$. We will use these terms generally; the function F[x, y, p] defined above is a particular example.

Simplest answer

When there is a singular solution, it is simply the discriminant of the equation:

Discriminant[F[x, y, p], p] == 0 // Simplify
(*  x y == 1  *)

Now, every surface that folds over itself will have such a discriminant, but it will not turn always out to be a solution to the corresponding differential equation $F(x,y,y')=0$. We will discuss why and when the discriminant yields a singular solution below.

What's going on?

The singular points of the projection from the surface $F=0$ to the $xy$ plane form the so-called criminant. It is the set of points on the surface for which $\partial F/\partial p = 0$. The discriminant curve is the projection of the criminant onto the $xy$ plane. (See, for instance, Arnold, Geometric Methods in the Theory of Differential Equations.) The discriminant may be obtained by eliminating $p$ from the system $$(1)\ \ F(x,y,p)=0,\quad\text{and}\quad(2)\ \ F_p(x,y,p)=0\,.$$

Eliminate[{F[x, y, p] == 0, D[F[x, y, p], p] == 0}, p]
(*  x y == 1  *)

Geometrically, this corresponds to where the tangent plane is vertical (parallel to the $p$ axis):

Mathematica graphics
Fig. $\nextfig\criminant$. The criminant (left) and its projection onto the discriminant curve (right).

GraphicsRow[{#, Show[#, ViewPoint -> {0, 0, Infinity}]}] &[
 ContourPlot3D[F[x, y, p] == 0,
  {x, 0.02, 2}, {y, 0.02, 2}, {p, -4, 0},
  ContourStyle -> Opacity[0.5],
  MeshFunctions -> {Function[{x, y, p}, x^2 p + 1]},
  Mesh -> {{0}}, MeshStyle -> {Red, Thick]},
  AxesLabel -> Automatic]
 ]
(* Because x y <= 1 always, we use a different form x^2 p == -1 in
   the MeshFunctions to draw the criminant. *)

If the discriminant defines $y$ as a smooth function of $x$, then we can deduce a relation for its derivative: $$ 0={d \over dx}\,F{(x,y,y')}=F_x+F_y\,y'+F_p\,y''=F_x+y'F_y \implies y'=-F_x/F_y\,. \tag{1} $$ Before discussing when the discriminant is a solution to the ODE, let us first look at its connection to the envelope of the general solution.

Why is the discriminant the same as the envelope of the general solution?

Suppose we have a smooth parametrization of the general solution in the form $y = y(x,c)$. Since $$ F(x,\, y(x,c),\, y'(x,c))=0 $$ (where $y'(x,c)$ is the derivative with respect to $x$), then we have at the (dis)criminant where $F_p=0$ $$ {\partial \over \partial c}\,F(x,\, y(x,c),\, y'(x,c)) = F_y\,{\partial\over\partial c}\,y(x,c) = 0 $$ Thus, provided $F_y \ne 0$, which is generally true, the discriminant satisfies the conditions $y = y(x,c),\ 0 = {\partial y / \partial c}$, which are the conditions for the envelope.

[Otherwise $F_y=F_p=0$ and $dF=0$ would imply $F_x=0$ and the surface would have a singular point.]

Mathematica graphics
Fig. $\nextfig\family$. The family of solutions for different values of the constant, plotted on the surface $F=0$ (left) and showing the envelope (right).

{sol} = DSolve[ode, y, x];

GraphicsRow[{#, Show[#, ViewPoint -> {0, 0, Infinity}]}] &[
  ParametricPlot3D[{x, y[x], y'[x]} /. sol /. C[1] -> 1/c // Simplify // Evaluate,
    {c, -1, -0.01}, {x, Max[0.02, -((1 + 2 c)/(2 c^2))], Min[2, (-0.005 - c)/c^2]},
    PlotPoints -> 35, PlotStyle -> Opacity[0.5],
    MeshFunctions -> {Function[{x, y, p, c, x2}, c]},
    MeshStyle -> {Thick}, BoxRatios -> {1, 1, 1},
    AxesLabel -> {x, y, p}] /. {foo = 0; mesh_Line :> {ColorData[97, ++foo], mesh}}
  ]

What is more, if any curve on the surface $F=0$ passes through the criminant, its projection onto the $xy$ plane intersects the discriminant in one of two ways. Either it is tangent or it has a singularity at the point of intersection. For at a point on the criminant we have exactly as in equation (1) $$0={dF}=F_x\;dx+F_y\;dy+F_p\;dp \implies dy/dx=-F_x/F_y\,,$$ and the solution will be tangent to the discriminant, provided it is nonsingular.

What's really going on?

We can construct a direction field on the surface $F=0$ by a system of planes, called a contact structure, as follows. At a point $(x,y,p)$ on the surface, the contact plane is given by $dx = p\; dy$, i.e., the plane through $(x,y,p)$ with normal $(1,-p,0)$. It is vertical and it intersects the $xy$ plane in a line with a slope $dy/dx$ that is given by the height $p$. A curve on the surface $F=0$ that is tangent to the contact plane at each point it passes through is called an integral curve of the contact structure. Its projection onto the $xy$ plane is an integral curve of the ODE $F=0$.

Mathematica graphics
Mathematica graphics
Fig. $\nextfig\contact$. Contact structure on the surface: $F=0$ (top), $p^2=y$ (left), and $p^2=x$ (right). The integral curves are shown on each surface and their projections, which are integral curves of the ODE are shown below.

At a point on the criminant there are two possibilities. The contact plane intersects the tangent plane in a vertical line, or the contact plane is identical to the tangent plane. The second case implies that the slope of the discriminant is given by $p = dy/dx$ and satisfies the ODE. If each contact plane of the criminant is tangent to the criminant, the discriminant curve will be an integral curve of the ODE. It should be clear that the contact planes coinciding with the tangent planes is an exceptional situation and to have all of them line up along the criminant is really unlikely. In that case, the contact planes do not define a direction on the surface, and technically, the criminant is not an integral curve, but a singular locus of the direction field.

A generic, small, smooth perturbation of the surface $F=0$ will transform the contact system to the first case. For instance, if the surface is translated up, the slope $p$ increases and the contact planes rotate but the tangent planes stay the same (Fig. $\nextfig\deformation$). If the surface is rotated, then the tangent planes rotate but not the contact planes.

Mathematica graphics
Fig. $\deformation$. Deformation of the contact structure (via translation). A slight raising or lowering changes the direction field, which is given by the height $p$ at each point. In turn the integral curves are altered, most dramatically near the discriminant curve, which itself is destroyed as a solution.

After such a deformation, the discriminant curve will no longer be solution of the differential equation. The contact plane of the transformed surface will define a vertical direction at the criminant. The integral curves of the contact structure will pass through the criminant, and their projections onto the $xy$ plane will have a cusp on the discriminant. Indeed this is the norm. Generically at the criminant an ODE can locally be put in the normal form (Fig. $\contact$, right) $$P^2=X \quad (P = dY/dX)$$ by suitable change of coordinates. Only exceptionally can it be put in the form (Fig. $\contact$, left) $$P^2=Y$$ for which the $X$ axis is the criminant and the exceptional/singular solution. Indeed the OP's ODE may be put in this form by the change of coordinates $$y= {1 \over x}-x\,Y,\quad x = 2/X\,.$$

(* change of variables *)
alg2Fn = HoldPattern[y_[x_] -> fx_] :> (y -> Function[{x}, fx]);
x2X = First@Solve[x == 2/X, x];                            (* change x to X *)
X2x = First@Solve[x == 2/X, X];                           (* change X to x *)
y2Y = First@Solve[y[x] == 1/x - x Y[2/x], y[x]] /. alg2Fn; (* change y to Y *)
Y2y = First@Solve[y[x] == 1/x - x Y[2/x] /. x2X, Y[X]] /. alg2Fn; (* change Y to y *)

(* mapping ode to P^2 == Y *)
odeY =  F[x, y[x], y'[x]] == 0 /. y2Y /. x2X // Simplify
(*  Y[X] == Y'[X]^2  *)

Note that DSolve does not find the singular solution Y -> Function[{X}, 0] to this equation either, even though it is "obvious." However, the general and singular solution transform under the change of coordinates back to the general and singular solution of of the ODE $F=0$, after adjusting the constants of integration:

Ysol = DSolve[odeY, Y, X]
y[x] /. y2Y /. x2X                      (* change y[x] to terms of Y, X *)
% /. Ysol /. X2x /. C[1] -> 2 C[1] // Expand
y[x] /. DSolve[F[x, y[x], y'[x]] == 0, y[x], x] /. C[1] -> 2/C[1] // Expand
y[x] /. DSolve[ode, y[x], x] /. C[1] -> 2/C[1] // Expand
(*
  {{Y -> Function[{X}, 1/4 (X^2 - 2 X C[1] + C[1]^2)]},
   {Y -> Function[{X}, 1/4 (X^2 + 2 X C[1] + C[1]^2)]}}

  X/2 - (2 Y[X])/X                          <-- y[x] in terms of Y, X

  {2 C[1] - x C[1]^2, -2 C[1] - x C[1]^2}   <-- Ysol in terms of y, x
  {2 C[1] - x C[1]^2, -2 C[1] - x C[1]^2}   <-- DSolve @ F == 0
  {-2 C[1] - x C[1]^2}                      <-- DSolve @ ode (OP's)
*)

There are a couple of things worth observing about these results. First, the DSolve solution of the ODE $F = 0$ returns two formulas, while the solution to ode has only one. This is due to the rationalization of ode that produced F. The two formulas are equivalent and may be reconciled by replacing C[1] -> - C[1]. Second, the missing singular solution of odeY is mapped to the missing singular solution of the OP's ode:

Y[X] == 0 /. Y2y /. X2x // Simplify
(*  y[x] == 1/x  *)

On solving for the singular solutions

The singular solution begins to seem special. Let us consider how one might solve for it. The simplest way to find the singular solution is to check whether the discriminant curve satisfies the differential equation. Another approach is to integrate the contact structure. At the end we will review the desirability of automatically solving for the singular solution.

Checking the discriminant

Instead of Eliminate, we can use three-argument form of Solve in which the third argument is a List of variables to be eliminated. This will return the discriminant as a list of formulas for y. We then can pick the ones that satisfy the ODE.

ClearAll[singularSolve];
singularSolve[de_Equal, form : y_[x_] | y_, x_] :=
 Module[{p, q, u, fn, rads, a, b, ode, discriminant},
  fn = de /. {y[x] -> y, y'[x] -> p, Equal -> Subtract};
  rads = DeleteDuplicates@ Cases[fn, Power[a_, b_Rational], Infinity];
  ode = First@ GroebnerBasis[
     Flatten@ {
       fn /. Power[a_, b_Rational] :> u[a, b],
       rads /. Power[a_, b_Rational] :> u[a, b]^Denominator[b] - a^Numerator[b]},
     rads /. Power[a_, b_Rational] :> u[a, b]];
  discriminant = Solve[{ode == 0, D[ode, p] == 0}, {y}, {p}];
  Pick[
   Which[
      MatchQ[form, y[x]], {y[x] -> (y /. #)},
      MatchQ[form, y], {y -> Function @@ {{x}, y /. #}},
      True, # /. p -> y'[x]
      ] & /@ discriminant,
   TrueQ@Simplify[ode == 0 /. Flatten@{#, p -> D[y /. #, x]}] & /@ discriminant
   ]
  ]

Examples:

singularSolve[ode, y, x]                 (* OP's ode *)
(*  {{y -> Function[{x}, 1/x]}}  *)

singularSolve[y'[x]^2 == y[x], y[x], x]  (* generic singular solution *)
(*  {{y[x] -> 0}}  *)

singularSolve[y'[x]^2 == x, y, x]        (* generic normal form (no singular curve) *)
(*  {}  *)

Integrating the contact structure

By solving the contact structure we would find four solutions to the ODE $F=0$: $$ \begin{align} y &= 1/x & &\text{singular solution} \\ y &= -4 (c + x)/c^2 & &\text{general solution} \\ y &= 0 & &\text{limit of general solution} \\ x &= 0 & &\text{projective closure of general solution} \\ \end{align} $$ We will explain how these arise below. The fourth one may or may not be admitted depending on one's viewpoint.

Let's illustrate how the singular curve arises in the simple case of $(dY/dX)^2=Y$. Practically by inspection one can see that the singular solution is $Y=0$, which is not represented by the general solution $Y={1\over4}(X-c)^2$ returned by DSolve nor a limit of it. If we solve the contact system, we will see how the singular solution falls out. The corresponding equations are $$ \begin{align} P^2 &= Y & &\text{contact manifold} \\ 2P \; dP &= dY & &\text{tangent field} \\ dY &= P \; dX & &\text{contact field} \\ \end{align} $$ whence $P(2\;dP-dX)=0$. One factor gives $P=0$, from which follows the singular curve $Y=0$. Integrating the other factor yields $2P=X-c$ or after squaring, $Y={1\over4}(X-c)^2$.

In general, the contact system to be solved is $$ \begin{align} F(x,y,p) &= 0 & &\text{contact manifold} \\ dF(x,y,p) &= 0 & &\text{tangent field} \\ dy &= p \; dx & &\text{contact field} \\ \end{align} $$ The general approach is to reduce the direction field equations, that is, the tangent and contact fields, to disjoint systems of algebraic and differential equations (via Reduce). The contact manifold $F(x,y,p)=0$ is used to eliminate variables where necessary. The differential equations are integrated. Then we can get $y$ in terms of $x$ and a constant of integration, eliminating $p$ as needed. It turns out the Reduce may produce a component of the system with more than one differential equation. For instance, in the OP's case I obtained four disjoint systems, the first two of which have two differential equations:

(p == 0 && dy == 0 && -2 + x y != 0 && dp == 0 && dx != 0) ||            (*=> y == 0 *)
 (p != 0 && dx == dy/p && 2 + p x^2 - x y != 0 && dp == 0 && dy != 0) || (*=> y == (4(c-x))/c^2 *)
 (x != 0 && p == (-2 + x y)/x^2 && p != 0 && dx == dy/p && dy != 0) ||   (*=> y == 1/x *)
 (y != 0 && x == 2/y && p == 0 && dy == 0 && dx != 0)                 (* => not a solution *)

Their integrations produce independent constants that need to be reconciled.

It is worth emphasizing that solving this system yields both the general and the singular solutions. In the OP's case, it turned up two solutions not yet mentioned. One is

y -> Function[{x}, 0]

This is the limit of the general solution returned by DSolve as C[1] -> Infinity. Another is x == 0: It can be considered an integral curve of the projective closure of the contact structure at p == ∞. It makes no sense as a solution for y, but if we think of the differential equation in terms of (1) x as a function of y and (2) p == 1/x'[y], then it is a solution of the differential equation F[x[y], y, 1/x'[y]] == 0 (after clearing denominators). In fact, F is symmetric in this way and this differential equation is essentially the same as the original with x and y switched.

DSolve[Numerator@Together@F[x[y], y, 1/x'[y]] == 0, x, y]
(*
  {{x -> Function[{y}, (4 (-y + C[1]))/C[1]^2]},
   {x -> Function[{y}, -((4 (y + C[1]))/C[1]^2)]}}
*)

One can also view the $y$ axis $x == 0$ as a limit of the solution as C[1] -> Infinity, or in terms of the original general solution as C[1] -> 0, which is indeed how it appears in Michael Seifert's animation.

Note: The results obtained depended on how the Reduce command was set up. I have had trouble developing an approach that did not need tweaking on a case by case basis. In the big picture it does not seem that important to do so.

The OP's ODE revisited

We promised to say something about the OP's original ode and the rationalized F[x, y[x], y'[x]] == 0. Two observations raise questions about the DSolve general solution and its relationship to the contact structure. One is that the direction field does not exist along the criminant. The other is that the surface $F=0$ folds over itself along the criminant, and the surface that lies above criminant corresponds to the OP's ode while the surface below corresponds to the conjugate ODE; and yet the traces of the general solution on $F=0$ cross the criminant.

The general solution is valid only over a restricted domain

In fact the general solution of the OP's ode is valid only up to the point it contacts the discriminant.

{sol} = DSolve[ode, y, x];

Reduce[ode && $y == y[x] /. sol, {y, x}, {C[1]}, Reals]   (* compute domain *)
(*  ($y < 0 && (1/$y <= x < 0 || x > 0)) || ($y > 0 && (x < 0 || 0 < x <= 1/$y))  *)

If we check the residual error of ode for the solution that passes through $(1,1)$, we see it is zero up until x == 1 at which point y == 1 and the solution contacts the discriminant.

Mathematica graphics
Fig. $\nextfig\residual$. After the solution touches the discriminant at x == 1, it no longer satisfies the ODE ode (left), but it does satisfy the ODE F[x, y[x], y'[x]] == 0 (right).

Block[{C},
 C[1] = C[1] /. First@Solve[y[1] == 1 /. sol, C[1]];
 GraphicsRow[{
   Plot[ode /. Equal -> Subtract /. sol, {x, 0.5, 2},
    PlotRange -> All, Frame -> True],
   Plot[F[x, y[x], y'[x]] /. sol, {x, 0.5, 2}, PlotRange -> 0.8*^-15,
    Frame -> True, WorkingPrecision -> 16]
   }, ImageSize -> 533]
 ]

Using Reduce instead of Solve within DSolve

One can get Mathematica to produce this restriction by setting the option Method -> Reduce to Solve as I did here.

With[{opts = Options[Solve]},
 Internal`WithLocalSettings[
  SetOptions[Solve, Method -> Reduce],
  {redsol} = DSolve[ode, y, x],
  SetOptions[Solve, opts]
  ]]

Solve::useq: The answer found by Solve contains equational condition(s)....A likely reason for this is that the solution set depends on branch cuts of Wolfram Language functions. >>

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

We simplify the condition further.

redsol /. Cases[
  redsol,
  ConditionalExpression[val_, cond_] :>
   (HoldPattern[ConditionalExpression[val, cond]] ->
     With[{c = Reduce[cond, x, Reals]}, ConditionalExpression[val, c]]),
  Infinity]
(*
  {y -> Function[{x},
    ConditionalExpression[-((4 (x + C[1]))/C[1]^2),
    (2 x + C[1] >= 0 && C[1] > 0) || (2 x + C[1] <= 0 && C[1] < 0)]]}
*)

The condition 2 x + C[1] == 0 is equivalent to the solution touching the discriminant y = 1/x:

Reduce[y[x] == 1/x /. redsol]
(*  C[1] != 0 && x == -(C[1]/2)  *)

Fig. $\nextfig\domain$ shows a plot of the valid solutions to ode. (One may observe the uniqueness of a solution with a given initial condition.)

Mathematica graphics
Fig. $\nextfig\domain$. The general solution from DSolve is valid only up until it contacts the singular solution.

Show[
 Plot[1/x, {x, 0, 2}, PlotStyle -> Gray, PlotRange -> {0, 2}, AspectRatio -> Automatic],
 Table[Plot[
   Evaluate@
    Block[{C},
     C[1] = C[1] /. First@Solve[y[x0] == 1/x0 /. redsol, C[1]];
     y[x] /. redsol],
   {x, 0, 2}, PlotStyle -> ColorData[97, Round[(x0 - 0.4)/0.2]],
   MeshFunctions -> {Function[{x, y}, x y]}, Mesh -> {{0.99999}},
   MeshStyle -> {PointSize[Medium], ColorData[97, Round[(x0 - 0.4)/0.2]]}],
  {x0, 0.6, 1.8, 0.2}]
 ]

Connection to the perturbed ODE

One can see in the animation above (Fig. $\deformation$) how the general solution is "cut" at the discriminant and reflected back along the discriminant as the translated ODE $F(x,y,y'-a)=0$ moves through $a=0$. For a small perturbation, two branches of a regular solution meet in a cusp at the discriminant, one of them closely tracking a regular solution (line) of the orginal ODE ($a=0$) and the other closely tracking the discriminant.

Mathematica graphics

Fig. $\nextfig\perturbed$. The perturbed ODE $F(x,y,y'-a)=0$ for $a = \pm0.01$. Each solution has a cusp on the discriminant. One branch follows closely a line of the original general solution (shown lightly shaded); the other follows closely the discriminant.

It is interesting that one can see the singular solution as a limit, not of the general solution to the OP's ODE, but of the ODEs in a neighborhood of $F(x,y,y')=0$. When the criminant aligns with the contact field, the singular solution appears as an independent solution and the solutions pass by the discriminant (or through the criminant).

Conclusion

The upshot is that the rationalized ODE admits as solutions the complete lines in the general solution to the OP's ode, which only admits rays starting at the singular solution.

Aside: Connection to numerical integration

The picture of the ODE as a contact system also illustrates the numerical integration of the ODE.

  • When the integration of a regular solution by NDSolve reaches the discriminant/singular solution, the system becomes stiff and integration can continue only along the singular solution. This is true for both ode and $F=0$; in the case of F[x, y[x], y'[x]] == 0, NDSolve first solves for y'[x] and effectively creates two ODEs, equivalent to ode and its conjugate.
  • The discriminant appears as boundary in ode in the form Sqrt[1 - x y[x]], which causes complex numbers to enter the calculations as round-off error and approximate steps lead to the boundary being violated.
  • Whether NDSolve tracks the singular solution on the boundary or leaves it to follow a regular solution depends on which branch of the ODE $F=0$ is being integrated and the direction of integration.
  • Rectifying (straightening) the singular solution (by transforming the ODE to the form Y'[X]^2 == Y[X]) aligns the steps taken by NDSolve with the direction of the singular solution Y[X] == 0. In this case, NDSolve naturally follows the singular solution and does not leave it, provided the solution lies exactly on Y == 0; otherwise, the effects of error will be similar to those in integrating ode.
  • If we integrate the contact structure to construct parametrizations of the integral curves, NDSolve seems to always integrate across the singularity at the criminant.

Mathematica graphics
Fig. $\nextfig\ndsolve$. Integrating the ODE $F=0$ with an initial condition on the singular solution. NSolve computes two solutions (dark green, dark red). The integral curve (light yellow) is traced until it meets the singular solution (light blue). Then it traces the discriminant, with increasing error the further from the initial condition. The values along the singular solution have a small imaginary component; only the real part is plotted.

solF = NDSolve[
   {F[x, y[x], y'[x]] == 0, y[1] == 1}, y, {x, 0.5, 2},
   "ExtrapolationHandler" -> {Indeterminate &, "WarningMessage" -> False},
   Method -> "StiffnessSwitching" (* aids integrations along singular solution *)
   ];
Plot @@ {
  Flatten[{1/x, 2 - x, Re[y[x] /. solF]}], Flatten[{x, 0.5, 2}],
  PlotStyle -> {Directive[Thickness[0.015], Opacity[0.5]],
    Directive[Thickness[0.015], Opacity[0.5]],
    Directive[Darker[Red, 0.5], Opacity[0.8]],
    Directive[Darker[Green, 0.66], Opacity[0.8]]},
  GridLines -> {{1}, {1}}
  }

The numerical integration of a related problem is discussed here: NDSolve solves this ordinary differential equation only "half-way"

DSolve and the complete solution

Why is it that DSolve misses a solution of this differential equation?

A little reflection will suggest why DSolve does not solve for the singular solution.

Why singular solutions are missing in the first place

  • Singular solutions are not a member or a limit of the general solution.
  • A singular solution is a rare coincidence.
  • The singular solution in the example fell out from algebraic equations obtained from differentiation and not from integrating differential equations.
  • Using Reduce does not guarantee all solutions to the contact system will be found. (This is an observed fact. I do not understand why it happened or whether another approach might always be successful.)

Why they are not automatically sought by DSolve

  • Singular solutions are rare, so one has to know when to look for them, or it will be mainly a waste of time.
  • Reducing systems can take a lot of time. Is it really worth it if most of the time the result will be the null set?

How to be sure that DSolve returned all possible solutions?

This is the OP's primary question. Leaving the difficulties of symbolic integration aside, there seem to be three steps to getting the complete answer (or as complete as possible in Mathematica.)

  • Use DSolve to get the general solution.
  • Check the limits of the general solution as the constants of integration approach the boundaries of their natural domains. (This is the issue with the ODE in DSolve not finding solution I expected.)
  • Check the discriminant or the envelope of the general solution.
$\endgroup$
7
  • 4
    $\begingroup$ Are you sure you didn't mean to publish this as an expository journal article instead of as a mere SE answer? :) $\endgroup$ Commented Jul 29, 2015 at 2:10
  • $\begingroup$ @J. M. I considered it. But it's fairly simple, in a way, even though "contact structure" is typically grad-level. It is actually fairly rich and simple, which makes it a nice example. -- I dread the corrections to such a long post, though. Although maybe no one will read it. $\endgroup$
    – Michael E2
    Commented Jul 29, 2015 at 2:19
  • $\begingroup$ Well, I read it and liked it; hence the joke. :) Mentioning Arnol'd was a nice touch. $\endgroup$ Commented Jul 29, 2015 at 2:44
  • $\begingroup$ @J. M. Thanks. I saw Arnol'd lecture a few times and got to meet him once. I really like(d) him mathematically. $\endgroup$
    – Michael E2
    Commented Jul 29, 2015 at 3:07
  • 3
    $\begingroup$ @MichaelE2, this is an amazing post, and deserves much more recognition than the 7 upvotes it currently has. $\endgroup$
    – SPPearce
    Commented Jul 25, 2016 at 8:52
7
$\begingroup$

The following isn't quite an answer to the question—rather than having DSolve remind you that you might not have all the solutions, the procedure below tries to find other ones (and succeeds in the case you've outlined.)


As I mentioned in my comment above, the "missing" solution is the envelope of the general family of curves that you found:

enter image description here

In general, if a differential equation generates a one-parameter family of curves with a continuous envelope, then the envelope is a solution. This is because at every point along the envelope, $y'$, $y$, and $x$ agree with some solution of the original differential equation.

To find the envelope of a family of curves given by $F(x,y,c) = 0$, where $c$ is our parameter, we require a non-zero solution to the system of equations $F = 0$ and $\partial F/\partial c = 0$. If a non-trivial solution to these equations exists, it will be a solution of the original ODE that is not in the family of curves. As an example, let's apply this to the family of functions you found:

soln = DSolve[y'[x] == y[x]/x - 2/x^2 + (2*Sqrt[1 - x y[x]])/x^2, y[x], x, GeneratedParameters -> (Subscript[c, #] &)]
enveqns = {ys[x] == (y[x] /. First[soln]), D[ys[x] == (y[x] /. First[soln]), Subscript[c, 1]]};
Eliminate[enveqns, Subscript[c, 1]]

(* x ys[x] == 1 *)

Thus, $y_s(x) = 1/x$ is the envelope of the curves, and as such is a singular solution to the differential equations.

As another example, let's do a family of non-linear equations that doesn't have an envelope. The solutions to $y' = y^2/x^2$ are $y=x/(1 - c_1 x)$

soln = DSolve[y'[x] == y[x]^2/x^2, y[x], x, GeneratedParameters -> (Subscript[c, #] &)]
(* rest of code as above *)

(* x^2 == 0 && ys[x] == x *)

In this case, the envelope consists of a single point $x = y = 0$, and so the envelope is not a valid solution.

I don't know whether all singular solutions (in the sense of not belonging to a general family of solutions) can always be found in this way. Still, it could be a useful technique to look for singular solutions, and I had fun coding it up.

$\endgroup$
1
  • $\begingroup$ A similar situation happens with the simpler Clairaut equation; DSolve[] will happily give you a general solution, but it will take some work to obtain the singular solution. $\endgroup$ Commented Jun 18, 2015 at 0:18

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.