# How to obtain these ODE solutions using DSolve?

V 12.1.1 on windows 10

Any one has any suggestions how to persuade Mathematica to obtain these solutions below?

I have small collection of such ODE's from textbooks. But do not know now how to obtain these special solutions.

I'll just show 2 for now. I think one method should works for all.

The problem is that one can not just follow the standard method on these, which is to obtain the general solution and then solve for the constant of integration using the initial conditions, since that leads to singularity.

## First

ode = y'[x] == (x^2 + y[x]^2)/(2 x^2)
ic = y[-1] == -1;
DSolve[ode, y[x], x] (*no problem finding general solution*)


DSolve[{ode, ic}, y[x], x]
(* {} *)


But a particular solution exist, which is y[x]==x:

 sol = y -> Function[{x}, x];
ic /. sol
(* True *)
ode /. sol
(* True *)


## Second

ode = (x + y[x]) + (x - y[x])*y'[x] == 0;
ic = y[0] == 0;
DSolve[ode, y[x], x]  (*no problem finding general solution*)


DSolve[{ode, ic}, y[x], x]
(* {} *)


But a particular solution exist, which is y[x]==(1+Sqrt[2])x:

sol = y -> Function[{x}, (1 + Sqrt[2]) x];
ic /. sol
(* True *)
ode /. sol // Simplify
(* True*)


ps. I tried the method to find singular solutions given in nonlinear-first-order-differential-equation but it did not find these.

• Thanks for the accept. Just to let you know, I woke up this morning and realized I hadn't updated projSolve to solve over the projective complex line. It solves over the projective real line. I was going to update the answer later when I get a chance. I don't think it matters in practice, but it might be less robust. Conceptually it bothers me. :) – Michael E2 Jul 4 '20 at 16:05
• Related (similar issue): mathematica.stackexchange.com/questions/57910/… – Michael E2 Jul 4 '20 at 21:18

Let $$F(x,y,y')=0$$ be the differential equation, and suppose $$y_C = y(x; C)$$ is a solution for any complex number $$C$$. Then $$F(x,y_C,y'(C))=0$$ for all $$C$$. Then $$0=\lim_{C\rightarrow\infty}F(x,y_C,y'(C)) \buildrel ? \over = F(x,\lim_{C\rightarrow\infty} y_C, \lim_{C\rightarrow\infty} y'_C)\,.$$ So if the limits of $$y_C$$ and $$y'_C$$ exist, the limit of $$y'_C$$ is the derivative of the limit of $$y_C$$, and the limit may be brought inside $$F$$, then the limit of $$y_C$$ will be a solution. The hypotheses are often true if the limit of $$y_C$$ exists.

Here is a way to hack into DSolve and try to solve for the solution at $$C=\infty$$. It does so projectively: Treat the parameter $$C = [v:w]$$ as an element of the (complex) projective line, so that $$\infty = [1:0]$$. We convert the infinite solution to a limit by defining an upvalue for a special head limitRule that invokes Limit when the rule is used in ReplaceAll. It took a little experimentation to determine the form of ReplaceAll that DSolve calls. This is specifically restricted to single-parameter solutions, that is, first-order equations. It could be extended to multiple parameters.

ClearAll[projSolve, limitRule, withProjectiveParameters];

projSolve[eq_, {v_}, rest___] :=
Module[{w, sol},
sol = Solve[eq, {v}, rest]; (* could skip to proj. solver *)
If[sol === {},(* solve over projective line if regular Solve[] fails *)
sol = Solve[
Flatten@{eq /.
v -> v/w, (w == 1 && -1 <= v <= 1) || (v == 1 &&
0 <= w < 1) || (v == -1 && 0 < w < 1)}, {v, w}, rest];
If[sol =!= {}, (* the only solution should be ComplexInfinity *)
sol = limitRule @@ (List /@ Thread[v -> (v/w /. sol)])]
];
sol
];

limitRule /:
ReplaceAll[HoldPattern[{v_ -> body_}],
limitRule[rules : {_ -> _} ..]] := v -> Limit[body, #] & /@ rules;

SetAttributes[withProjectiveParameters, HoldFirst];
withProjectiveParameters[ds_DSolve] :=
InternalInheritedBlock[{Solve}, Unprotect[Solve];
call : Solve[eq_, v_, opts___] /; ! TrueQ[$$in] := Block[{$$in = True},
Hold[call] /.
Hold[Solve[e_, {c_}, o___]] /; ! FreeQ[e, _C] :>
projSolve[e, {c}, o] // ReleaseHold
];
Protect[Solve];
ds
]


Example 1:

ode = y'[x] == (x^2 + y[x]^2)/(2 x^2);
ic = y[-1] == -1;
withProjectiveParameters[DSolve[{ode, ic}, y, x]]

(*  {y -> Function[{x}, x]}  *)


Example 2:

ode = (x + y[x]) + (x - y[x])*y'[x] == 0;
ic = y[0] == 0;
withProjectiveParameters[DSolve[{ode, ic}, y, x]]


Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

DSolve::bvnul: For some branches of the general solution, the given boundary conditions lead to an empty solution.

(*  {{y -> Function[{x}, x - Sqrt[2] Sqrt[x^2]]}}  *)

{ode, ic} /. % // Simplify

(*  {{True, True}}  *)

Clear["Global*"]


Example 1

ode1 = y'[x] == (x^2 + y[x]^2)/(2 x^2);
ic1 = y[-1] == -1;


The general solution is

solg1 = DSolve[ode1, y, x][[1]]

(* {y -> Function[{x}, (x (-2 + 2 C[1] + Log[x]))/(2 C[1] + Log[x])]} *)


Verifying the general solution,

ode1 /. solg1 // Simplify

(* True *)


The particular solution is the limiting case as C[1] -> Infinity

test1 = Limit[{y[x], y[-1]} /. solg1, C[1] -> Infinity]

(* {x, -1} *)

solp1 = y -> Function[{x}, Evaluate@#[[1]]] &@test1

(* y -> Function[{x}, x] *)


Verifying the particular solution,

{ode1, ic1} /. solp1

(* {True, True} *)


EDIT: Alternatively, temporarily generalize the initial condition

ic1r = y[-1] == k;

solp1r = y -> Function[{x},
Evaluate[y[x] /. DSolve[{ode1, ic1r}, y, x][[1]] /. k -> -1]]

(* y -> Function[{x}, x] *)


This is identical to solp1

solp1r === solp1

(* True *)


Example 2

ode2 = (x + y[x]) + (x - y[x])*y'[x] == 0;
ic2 = y[0] == 0;


The general solutions are

solg2 = DSolve[ode2, y, x]

(* {{y -> Function[{x}, x - Sqrt[E^(2 C[1]) + 2 x^2]]}, {y ->
Function[{x}, x + Sqrt[E^(2 C[1]) + 2 x^2]]}} *)


Verifying the general solutions,

ode2 /. solg2 // Simplify

(* {True, True} *)


The particular solutions are the limiting cases as C[1] -> -Infinity

test2 = Limit[{y[x], y[0]} /. solg2, C[1] -> -Infinity]

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

solp2 = {y -> Function[{x}, Evaluate@#[[1]]]} & /@ test2

(* {{y -> Function[{x}, x - Sqrt[2] Sqrt[x^2]]}, {y ->
Function[{x}, x + Sqrt[2] Sqrt[x^2]]}} *)


Verifying the particular solutions,

{ode2, ic2} /. solp2 // Simplify

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


EDIT: Alternatively, temporarily generalize the initial condition

ic2r = y[0] == k;

solp1r = y  -> Function[{x}, Evaluate[y[x] /.
DSolve[{ode2, ic2r}, y, x][[1]] /. k -> 0]]

(* Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

y -> Function[{x}, x - Sqrt[2] Sqrt[x^2]] *)


Note that this misses one of the particular solutions.

• Thanks Bob. But how, in a program, would one know to take such a specific limit of C[1] to infinity? I mean, I am doing this in a program, and need an algorithm. How did you know that limit of C[1] to infinity will give the correct solution? If this was by trying different things, that is OK, but in a program, I would need a mathemtical reason or algorithm in order to check for such limit. – Nasser Jul 3 '20 at 18:34
• See edit. Try temporarily generalizing the initial condition. – Bob Hanlon Jul 3 '20 at 19:44
• @Nasser I tried to address your question in the opening paragraph of my answer. In general one would know to try the limit at infinity if the condition for the parameter C[1] has no solution. Then there might or might not be a solution at infinity. If the IVP has a solution and there is no solution for finite C[1], then there has to be one for infinite C[1]. Sometimes the limits at $\pm\infty$ are different, but my answer ignores that case and assumes the ode is sufficiently nice over the complexes. – Michael E2 Jul 3 '20 at 20:50