3
$\begingroup$

I reach a dead end when trying to solve the following differential equation:

$$1+x-x^2=\frac{\partial}{\partial x}\Bigg\{\frac{2}{y(x)}\bigg[1+\bigg(\frac{1+0.001x}{2y(x)}\bigg)^2\bigg]^{-1}\Bigg\}$$

I use DSolve for this, with the boundary condition $y(0)=0$:

DSolve[{1 + x - x^2 == D[2/(y[x]*(((1 + 0.001*x)/(2*y[x]))^2 + 1)), {x}], y[0] == 0}, y[x], x]

For some reason, Mathematica gets stuck trying to calculate this, without obtaining any answer at all. The problem is clearly with the $(1+0.001x)$ part of the equation, because if I change this into simply $x$, meaning I solve the following equation:

$$1+x-x^2=\frac{\partial}{\partial x}\Bigg\{\frac{2}{y(x)}\bigg[1+\bigg(\frac{x}{2y(x)}\bigg)^2\bigg]^{-1}\Bigg\}$$

DSolve[{1 + x - x^2 == D[2/(y[x]*((x/(2*y[x]))^2 + 1)), {x}], y[0] == 0}, y[x], x]

It yields a solution within seconds.

Even if I get rid of the small floating number $0.001$, which intuitively seems to be the cause of the mishap, and simply replace it with $(1+x)$, it still gets stuck trying to return an answer.

Why is this happening?

$\endgroup$
7
$\begingroup$

This does not get stuck, it returns right away

ClearAll[x,y]
z=1/1000;
DSolve[{1+x-x^2==D[2/(y[x]*(((1+ z x)/(2*y[x]))^2+1)),{x}],y[0]==0},y[x],x]

Mathematica graphics

But your input makes it stuck

ClearAll[x,y]
z=0.001;
DSolve[{1+x-x^2==D[2/(y[x]*(((1+ z x)/(2*y[x]))^2+1)),{x}],y[0]==0},y[x],x]

As a general rules, use exact input with symbolic functions like DSolve, always safer.

To answer comment

You can see why it does not like y[0]=0

ClearAll[x, y, a]
z = 1/1000;
sol = y[x] /. 
  First@DSolve[{1 + x - x^2 == 
      D[2/(y[x]*(((1 + z x)/(2*y[x]))^2 + 1)), {x}], y[0] == a}, y[x],
     x]

sol2 = Limit[sol, a -> 0];

Mathematica graphics

And now at x=0 the above is

sol2/.x->0

Mathematica graphics

So you boundary condition y[0]=0 leads to a problem. May be this is why DSolve is having trouble with this. Changing the boundary conditions so it is little away from zero, now DSolve gives solution (not shown)

sol=y[x]/.First@DSolve[{1+x-x^2==D[2/(y[x]*(((1+z x)/(2*y[x]))^2+1)),
       {x}],y[0]==1/1000000},y[x],x]

Basically, your boundary condition at x=0 does not seem to be consistent with the ODE. But may be more analysis is needed, this is just a quick look.


update

Thanks to comment by bbgodfrey, changing sol2/.x->0 above by Limit[sol2,x->0] gives 0 instead of Indeterminate

Mathematica graphics

$\endgroup$
  • $\begingroup$ Thanks for the tip. But while you're right that it doesn't get stuck with exact input, it returns a blank output line... $\endgroup$ – VortexSheet Jul 29 '17 at 11:51
  • $\begingroup$ Ok, if I get rid of the boundary condition $y[0]==0$, I get an answer. Anybody knows why this is happening? $\endgroup$ – VortexSheet Jul 29 '17 at 11:58
  • 1
    $\begingroup$ Try Limit[sol2, x -> 0] instead of sol2/.x->0 to obtain 0, as desired. $\endgroup$ – bbgodfrey Jul 29 '17 at 15:37
  • $\begingroup$ @bbgodfrey thanks. Updated per your comment. $\endgroup$ – Nasser Jul 29 '17 at 19:18
5
$\begingroup$

The solution of the ODE without a boundary condition is

s = y[x] /. DSolve[1 + x - x^2 == D[2/(y[x]*(((1 + x z)/(2*y[x]))^2 + 1)), x], y[x], x]
(* {(-384 - √(147456 - 4 (-192 x - 96 x^2 + 64 x^3 - 12 C[1]) (-48 x - 24 x^2 + 
     6 x^3 - 96 x^2 z - 48 x^3 z + 32 x^4 z - 48 x^3 z^2 - 24 x^4 z^2 + 16 x^5 z^2 - 
     3 C[1] - 6 x z C[1] - 3 x^2 z^2 C[1])))/(2 (-192 x - 96 x^2 + 64 x^3 - 12 C[1])), 
    (-384 + √(147456 - 4 (-192 x - 96 x^2 + 64 x^3 - 12 C[1]) (-48 x - 24 x^2 + 
     16 x^3 - 96 x^2 z - 48 x^3 z + 32 x^4 z - 48 x^3 z^2 - 24 x^4 z^2 + 16 x^5 z^2 - 
     3 C[1] - 6 x z C[1] - 3 x^2 z^2 C[1])))/(2 (-192 x - 96 x^2 + 64 x^3 - 12 C[1]))} *)

Now consider the values of both solutions of y at x == 0.

Simplify[s /. x -> 0]
(* {(32 + Sqrt[1024 - C[1]^2])/(2 C[1]), -((-32 + Sqrt[1024 - C[1]^2])/(2 C[1]))} *)

The first clearly cannot equal zero for any value of the constant C[1]. However,

Series[Last@%, {C[1], 0, 3}] // Normal    
(* C[1]/128 + C[1]^3/524288 *)

So, C[1] == 0 causes the second solution for y[x] to vanish at the origin.

Plot[Last[s /. {C[1] -> 0, z -> 1/1000}], {x, -1, 1}, ImageSize -> Large, 
    AxesLabel -> {x, y}, LabelStyle -> Directive[12, Bold, Black]]

enter image description here

Evidently, DSolve did not recognize this. It is not uncommon for DSolve to have difficulties applying boundary conditions.

Addendum

Another way to see that the ODE is well behaved near y[0] == 0 is to expand the ODE at x == 0 and solve for higher-order derivatives. For instance,

Solve[Thread[CoefficientList[(Series[1 + x - x^2 - D[2/(y[x]*(((1 + z*x)/
    (2*y[x]))^2 + 1)), x], {x, 0, 2}] // Normal) /. y[0] -> 0, x] == 0], 
    {y'[0], y''[0], y'''[0]}]
(* {{Derivative[1][y][0] -> 1/8, 
     Derivative[2][y][0] -> (1 + 4*z)/8, 
     Derivative[3][y][0] -> (-13 + 48*z + 48*z^2)/64}} *)

It also is possible to integrate the ODE numerically with NDSolve, although some care must be taken, because the ODE is singular at x == 0, even though its solution is not.

$\endgroup$
  • $\begingroup$ That really helped, thanks! $\endgroup$ – VortexSheet Jul 30 '17 at 15:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.