6
$\begingroup$

Update June 25, 2023. The problem is still present in V 13.3

enter image description here


I am baffled why DSolve in V 13.1 unable to solve this standard first order ode. DSolve has improved a lot in recent version, so I was expecting it to be to solve this first order ode.

I am looking for suggestion what happend here, or a workaround.

I will show my hand solution. Mathematica verifies the solutions. Then show an attempt for workaround. DSolve seems to hang. Well, I waited 10 minutes and gave up. May be if I wait much more it will solve it. But from my experience with DSolve, if it takes more than few minutes, most likly it is stuck somewhere and no point waiting any more.

I do not understand why it needs that much time. Maple also solves this ode instantly. So I think something went wrong here with DSolve code path somewhere.

The ode is

ClearAll[y, x]
ode = y[x] == x + y'[x]^2*(1 - 2/3*y'[x])
DSolve[ode, y[x], x]

Waiting and waiting....but why?

Here is my hand solution, you can see there is nothing that should cause long time delay.

hand Solution

\begin{align*} y & =x+\left( y^{\prime}\right) ^{2}\left( 1-\frac{2}{3}y^{\prime}\right)\\ & =\phi\left( x,y^{\prime}\right) \end{align*} Let $y^{\prime}=p$ \begin{align} y & =x+\left( p^{2}-\frac{2}{3}p^{3}\right) \tag{1}\\ & =\phi\left( x,p\right) \nonumber \end{align} Differentiating w.r.t. $x$ gives \begin{align} y^{\prime} & =\frac{\partial\phi}{\partial x}+\frac{\partial\phi}{\partial p}\frac{dp}{dx}\nonumber\\ p & =1+\left( 2p-2p^{2}\right) \frac{dp}{dx}\nonumber\\ \frac{dp}{dx} & =\frac{p-1}{2\left( p-p^{2}\right) }\tag{2} \end{align} The singular solution is when $\frac{dp}{dx}=0$ which results in $p=1$. Substituting this in (1) gives \begin{align*} y & =x-\left( 1-\frac{2}{3}\right) \\ & =x+\frac{1}{3} \end{align*} The general solution is when $\frac{dp}{dx}\neq0$. Then (2) is now separable. Solving for $p$ gives \begin{align*} p & =-\sqrt{c_{1}-x}\\ p & =\sqrt{c_{1}-x} \end{align*} Substituting each one of the above solutions of $p$ in (1) gives \begin{align*} y_{1} & =x+\left( p^{2}-\frac{2}{3}p^{3}\right) \\ & =x+\left( \left( -\sqrt{c_{1}-x}\right) ^{2}-\frac{2}{3}\left( -\sqrt{c_{1}-x}\right) ^{3}\right) \\ & =x+\left( c_{1}-x+\frac{2}{3}\left( c_{1}-x\right) ^{\frac{3}{2}}\right) \\ & =c_{1}+\frac{2}{3}\left( c_{1}-x\right) ^{\frac{3}{2}}% \end{align*} And \begin{align*} y_{2} & =x+\left( p^{2}-\frac{2}{3}p^{3}\right) \\ & =x+\left( \left( \sqrt{c_{1}-x}\right) ^{2}-\frac{2}{3}\left( \sqrt{c_{1}-x}\right) ^{3}\right) \\ & =x+\left( c_{1}-x-\frac{2}{3}\left( c_{1}-x\right) ^{\frac{3}{2}}\right) \\ & =c_{1}-\frac{2}{3}\left( c_{1}-x\right) ^{\frac{3}{2}}% \end{align*} Therefore the solutions are \begin{align*} y & =x+\frac{1}{3}\\ y & =c_{1}+\frac{2}{3}\left( c_{1}-x\right) ^{\frac{3}{2}}\\ y & =c_{1}-\frac{2}{3}\left( c_{1}-x\right) ^{\frac{3}{2}} \end{align*}

Verification of hand solutions by Mathematica

ClearAll[y, x]
ode = y[x] == x + y'[x]^2*(1 - 2/3*y'[x])
ode /. y -> Function[{x}, x + 1/3]

Mathematica graphics

ode /. y -> Function[{x}, C[1] + 2/3 (C[1] - x)^(3/2)] // Simplify

Mathematica graphics

ode /. y -> Function[{x}, C[1] - 2/3 (C[1] - x)^(3/2)] // Simplify

Mathematica graphics

Attempts for workaround

I tried to see if solving for $y'$ first and solving each separate ode will help (even though this is the wrong approach to solving d'Alembert) and solve each generated ode (which is now of degree one), but this did not work. It also hanged

ClearAll[y, x]
ode = y[x] == x + y'[x]^2*(1 - 2/3*y'[x])
sol = y'[x] /. Solve[ode, y'[x]];
DSolve[y'[x] == First@sol, y[x], x]

Based on the above, my guess is that, this is what DSolve did internally? And that is why it hanged? Since the integrals generated are way too complicated to be solved as a quadrature ode.

Mathematica graphics

Is it possible this is what DSolve did?

And finally, this is Maple's result

ode:=y(x)=x+ diff(y(x),x)^2*(1-2/3*diff(y(x),x));
s := time[real]():
dsolve(ode);
time[real]()-s;

enter image description here

Any hints why DSolve can't solve this, or a workaround are welcome.

V 13.1 on windows 10.

$\endgroup$

2 Answers 2

2
$\begingroup$

Update notice: @lotus2019 pointed out that it was not working for them. The issue, as much as one can be sure of issues when fiddling with internals, is that the time constraint needs to be short enough. Tested on V14.0 and V14.1. I checked DSolve[ode, y[x], x] in V14.1, and it still results in an implicit equation with an unsolved integral.

Workaround:

ClearAll[y, x]
ode = y[x] == x + y'[x]^2*(1 - 2/3*y'[x]);
Internal`InheritedBlock[{Solve},
   (* ref: https://mathematica.stackexchange.com/a/63676 *)
   Unprotect[Solve];
   s_Solve /; ! TrueQ[$in] := Block[{$in = True},
     TimeConstrained[s, 0.1, $Failed] (* might need to adjust 0.1 *)
     ];
   Protect[Solve];
   DSolve[ode, y[x], x]
   ] //(* difficult simplification *)
  Collect[#, Sqrt[-x + C[1]], 
    ReplaceAll[c_ (x - C[1]) -> -c (-x + C[1])]@*Factor] & //
 AbsoluteTiming
(*
{84.2428 (* first time; 2.4242 sec. on subsequent runs *), 
{{y[x] -> C[1] - 2/3 (-x + C[1])^(3/2)},
 {y[x] -> C[1] + 2/3 (-x + C[1])^(3/2)}}}
*)

IncludeSingularSolutions -> True returns the same, after 50(!) additional seconds. So that seems a place for improvement.

Possibly the first strategies the parser tries waste a long time before returning unsolved. That would explain why a short time constraint speeds things up instead of causing a quick failure.


Addendum

Hackier but complete workaround as a proof of concept (V14.0.0 and the version might be more important here since we edit out the bug in the code):

ClearAll[y, x]
ode = y[x] == x + y'[x]^2*(1 - 2/3*y'[x]);
Internal`InheritedBlock[{DSolve`DSolveSingularSolutions, 
   TimeConstrained, Limit, Solve},
  DSolve`DSolveSingularSolutions;
  Unprotect[DSolve`DSolveSingularSolutions];
  DownValues@DSolve`DSolveSingularSolutions = 
   DownValues@DSolve`DSolveSingularSolutions /.
    HoldPattern[ (* BUG FIX *)
      If[! FreeQ[test_, False], s_ = {}]] :> (s = Pick[s, test]);
  Protect[DSolve`DSolveSingularSolutions];
  Unprotect[TimeConstrained];
  TimeConstrained[code_, tc_, fail___] /; ! TrueQ[$inTC] :=
   Block[{$inTC = True},
    If[FreeQ[Hold[code], DSolve`DSolveSingularSolutions],
     TimeConstrained[code, 0.1, fail], (* might need to adjust 0.1 *)
     (* Don't over-constrain DSolveSingularSolutions: *)
     TimeConstrained[code, tc, fail]
     ]
    ];
  Protect[TimeConstrained];
  Unprotect[Limit];
  lim_Limit /; ! TrueQ[$inLimit] := 
   Block[{$inLimit = True, $res1, $res2},
    TimeConstrained[lim, 1, Infinity]
    ];
  Protect[Limit];
  Unprotect[Solve];
  s_Solve /; ! TrueQ[$inSolve] := 
   Block[{$inSolve = True, $res1, $res2},
    TimeConstrained[s, 1, $Failed]
    ];
  Protect[Limit];
  
  DSolve[ode, y[x], x, IncludeSingularSolutions -> True]
  ] // AbsoluteTiming
(* unsimplified
{69.4 (* first time; 5.1 sec. on subsequent runs *), 
 {{y[x] -> 1/3 (3 C[1] + 2 x Sqrt[-x + C[1]] - 2 C[1] Sqrt[-x + C[1]])},
  {y[x] -> 1/3 (3 C[1] - 2 x Sqrt[-x + C[1]] + 2 C[1] Sqrt[-x + C[1]])},
  {y[x] -> 1/3 (1 + 3 x)}}}
*)

The time limits work here because everything important is easy, except the total time for DSolve`DSolveSingularSolutions. The small time limits seem to keep DSolve from going down several rabbit trails for too long. Beware, the time limits might cause DSolve to miss the solution on a different problem. Clearly the improved time on the second and subsequent runs means some work from the first run is re-used. One hopes that the short time limits did not cause failures to be cached. This partly what I mean by warning that this a proof of concept. This is not a general approach to be emulated by the casual user nor certainly by application developers.

The concept proved is that DSolve[] can solve the ODE. However, design choices make finding the general solution very slow, and an apparent bug means the singular solution is discarded.

$\endgroup$
7
  • $\begingroup$ are you using internal version of 14.1? Because for me on V 14. on windows 10, your code takes 84 seconds to complete. Not 2.4 seconds? Here is screen shot FYI i.sstatic.net/BHTK8Luz.png tried it 2 times., same timing. From clean kernel also. $\endgroup$
    – Nasser
    Commented Jul 20 at 18:38
  • $\begingroup$ @Nasser No, I'm using V14.0.0 (Mac ARM). It looks like some internal cache gets set. When I time it from a fresh kernel, it's 84 sec. for me, too. Second time, it's 2.4 sec. Probably something was cached when I tried your DSolve code without any changed. I can't explain why your 2nd run is still slow. I doubt it's a system thing. $\endgroup$
    – Michael E2
    Commented Jul 20 at 18:52
  • $\begingroup$ My second run was fast, 2.4 seconds. I only tried it twice, each from Fresh kernel. Yes, it looks like caching issues, where first time, some code gets loaded first time. Ok. $\endgroup$
    – Nasser
    Commented Jul 20 at 18:55
  • 1
    $\begingroup$ @Nasser BTW, it looks like there's a bug in the singular solver. It gets a list of candidates and tests them. It throws out all of them if any of them fail. In your case, it gets your solution and a spurious one, $y=x$. Consequently, it throws out the one works. $\endgroup$
    – Michael E2
    Commented Jul 20 at 18:57
  • 1
    $\begingroup$ @Nasser I coaxed (bullied?) the singular solution out of DSolve. A pretty severe hack! Not really of general use, I'm afraid. $\endgroup$
    – Michael E2
    Commented Jul 21 at 3:34
2
$\begingroup$

D'Alembert and Clairaut as trivial examples of first order ODE are not even mentioned by Arnold. With $y'=p$ make it explicit

  Solve[y == p^2 *(1 - 2/3) p, p]
   +-(-1)^(2/3) 3^(1/3) y^(1/3)

$$\text{DSolve}\left[y'(x)=\pm (-1)^{2/3} \sqrt[3]{3} \sqrt[3]{y(x)},y(x),x\right]$$ $$\left\{y(x)\to \pm \frac{2}{3} \sqrt{\frac{2}{3}} \left((-1)^{2/3} \sqrt[3]{3} x+c_1\right){}^{3/2}\right\}$$

$\endgroup$

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.