NDSolve giving the wrong solution?

Question

I'm considering the non-linear second order ODE DE $=0$, with DE given by

DE = y[x] (2 + y[x]) y''[x] + 2 x y[x] (2 + y[x]) y'[x]/(-1 + x^2) - (1 + 2 y[x]) y'[x]^2  
+ 4 (1 + y[x]) (2 + y[x])^2/(-1 + x^2)^2.

on the domain $x\geq 0$. The initial conditions are $y[0]=0,$ $y'[0]=-4$. Now I know that a solution is given by $ysol=\frac{-4x}{1+3x}$, which Mathematica can indeed verify. Now I'm trying to acquire this solution with NDSolve. When I use the following input

NDSolve[{DE == 0, y[0] == 0, y'[0] == -4}, y, {x, 0, 10}]

Mathematica gives the error

Power::infy: "Infinite expression 1/0. encountered.

Now I am assuming out that the problem is the fact that the coefficient in the ODE in front of the second derivative vanishes at $x=0$. Therefore I tried using a slightly different initial condition. The input I used is:

eps =10^-10
Sol = NDSolve[{DE == 0, y[0] == -eps, y'[0] == -4}, y, {x, 0, 10}];
Plot[{y[x] /. Sol, ysol[x]}, {x, 0, 1}]

Here $ysol$ is the known solution I gave above. The output is

hence Mathematica does not give the solution I already know.

Questions:

• Is the solution Mathematica is giving me correct because of the slightly different initial conditions?

• How can I make Mathematica give the solution I already know? I.e. what should I change in my input?

• It seems Mathematica cannot solve the ODE with the initial conditions that I know. Is this an error caused by the numerical solution method, or does this mean that the ODE actually does not have a unique solution?

Answer 1

The DE has problems at y[x] == 0 as noted, but also at x == 1, where there is a pole of order 2. This suggests that numerical issues near the singularities could cause substantial error to accumulate. One can use Piecewise to substitute the limiting value at a discontinuity, but this won't cure any numerical instability in the neighborhood of the discontinuity.

One can remove the 1/0 errors with "EquationSimplification" -> "Residual", but the numerical trouble persists. The code below results in a plot similar to, if not worse than, the OP's workaround.

sol = NDSolve[{DE == 0, y[0] == 0, y'[0] == -4, y''[0] == 24}, 
  y, {x, 0, 10}, Method -> {"EquationSimplification" -> "Residual"}]

Plot[{(-4 x)/(1 + 3 x), y[x] /. First[sol]}, {x, 0, 10}, 
 PlotStyle -> {AbsoluteThickness[5], Automatic}]

One can see that the second derivative is not calculated very well near x == 0 by this method:

foo = D[(-4 x)/(1 + 3 x), x, x];
Plot[{y''[x] /. First[sol], foo}, {x, 0, 0.001}, PlotRange -> All]

Using Piecewise helps near x == 0, but it fails at x == 1:

sys = {y''[x] == 
   Piecewise[{{y''[x] /. First@Solve[DE == 0, y''[x]] // Expand // 
       Simplify, y[x] != 0}}, 24], y[0] == 0, y'[0] == -4};
sol = NDSolve[sys, y, {x, 0, 10}];

NDSolve::ndsz: At x == 1.00008646582443`, step size is effectively zero; singularity or stiff system suspected. >>

Plot[{(-4 x)/(1 + 3 x), y[x] /. First[sol]}, {x, 0, 1.}, 
 PlotStyle -> {AbsoluteThickness[5], Automatic}]

To get a solution that crosses the discontinuity is difficult. The numerical evidence* suggests the discontinuity is unstable and virtually impossible without some preconditioning. The best thing I have to offer is to start the integration at x == 1. While there are some numerical issues near x == 0, integrating toward it seems stable.
(*I tried WhenEvent[x == 1, "CrossDiscontinuity"], increasing WorkingPrecision, and such stuff. Whenever it integrated past x == 1, there was a large error.)

To start the integration at x == 1, I tried to analyze the DE and come up with a solution (pretending not to know the OP's desired outcome). It turns out there are two choices for y[1], so that didn't work. And for the desired y[1] == -1, the choice for y'[1] appears to be free, so again I had to use the OP's formula to determine the initial value. I determined these constraints by examining the series expansion of the DE. The default value for y''[1] is also determined this way.

SeriesCoefficient[DE, {x, 1, -2}] == 0 // Solve
SeriesCoefficient[DE, {x, 1, -1}] == 0 /. y[1] -> -1 // Solve
SeriesCoefficient[DE, {x, 1, 0}] == 0 /. y[1] -> -1 // Solve
(*
  {{y[1] -> -2}, {y[1] -> -2}, {y[1] -> -1}}
  {{}}
  {{y''[1] -> y'[1] (-1 + 2 y'[1])}}
*)

Note that "StiffnessSwitching" is needed to step through the muck around x == 1:

ypp[x0_?NumericQ, y0_?NumericQ, yp0_?NumericQ] = Piecewise[
   {{y''[x] /. First@Solve[DE == 0, y''[x]] /. {x -> x0, y[x] -> y0, 
       y'[x] -> yp0}, x0 != 1}},
   yp0 (-1 + 2 yp0)];
sol = NDSolve[{y''[x] == ypp[x, y[x], y'[x]], y[1] == -1, 
   y'[1] == -1/4}, y, {x, 0, 10}, Method -> "StiffnessSwitching"]

Plot[{(-4 x)/(1 + 3 x), y[x] /. First[sol]}, {x, 0, 10}, 
 PlotStyle -> {AbsoluteThickness[5], Automatic}]

Answer 2

I am not an expert on this (solving non-linear odes'). But to put all what I saw here in stead of in comment is easier. Basically, you have a non-linear ode with singularity when y(0)=0 as well due to the initial conditions. You can see this part like this:

de = y[x] (2 + y[x]) y''[x] + 2 x y[x] (2 + y[x]) y'[x]/(-1 + x^2) - 
  (1 + 2 y[x]) y'[x]^2 + 4 (1 + y[x]) (2 + y[x])^2/(-1 + x^2)^2;
sol=Solve[de == 0, y''[x]]

What happens in the above when y(0)=0?

rhs = y''[x] /. sol
rhs /. {x -> 0, y[x] -> 0, y'[0] -> -4}

So, your initial conditions leads to singularity, without even solving the ode. i.e. the above is saying y''(0) is infinity. But if y'(0)=-4, then it y''(0) should be zero. right? So, your initial conditions and the ode do not agree, at least from initial look as in the above.

But since this is a non-linear ode, they have what is called a singular solutions, if you google singular solution of nonlinear ode you see articles on this. These do not come from setting initial conditions to the general solution. Mathematica DSolve can't solve this in closed form, and NDsolve can't due to the above problem with initial conditions (after all, it is a numerical solver).

Maple can give a singular solution but parametrized. By changing the parameters, it gave me back the solutions I gave. Here is the code if you are interested:

Taking this solution to Mathematica and looking at it:

sol = -ArcTanh[x] + (1/2)*ArcTanh[(1/4)*(4*C1*(2 + y[x]) + 8)/(4 - 2*C1*y[x]^2 -
  4*C1*y[x])^(1/2)] - (1/2)*ArcTanh[C1 + 1] == 0

Now solving for y[x] and changing C1 gives the different solutions I showed

sol2 = y[x] /. Solve[sol, y[x]];  (*2 solutions*)

For example, looking at the first one, change C1

Limit[First@sol2, C1 -> Infinity]
(*  -2  *)

Limit[First@sol2, C1 -> 1]

Notice that at C1=0 the solution is not defined.

Used Maple 18.01 and Mathematica 10.01 both on windows.