# Mathematica gets stuck trying to solve a simple differential equation

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?

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]


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.

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];


And now at x=0 the above is

sol2/.x->0


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

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

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]]


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

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.

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