3
$\begingroup$

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

Output

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?

Improve this question
$\endgroup$
  • 1
    $\begingroup$ do not have time now, but quick comment: you have non-linear ode. A non-linear ode can admit a singular solution (which can't be obtained from the general solution by giving the constants of integration specific values). Another solution is y=-2. Another is y=(-4*((x^4-2*x^2+1)/(x^4+8*x^3+18*x^2+8*x+1))^(1/2)*x^4-20*((x^4-2*x^2+1)/(x^4+8*x^3+18*x^2+8*x+1))^(1/2)*x^3-4*x^4-24*((x^4-2*x^2+1)/(x^4+8*x^3+18*x^2+8*x+1))^(1/2)*x^2-16*x^3-20*((x^4-2*x^2+1)/(x^4+8*x^3+18*x^2+8*x+1))^(1/2)*x-32*x^2-4*((x^4-2*x^2+1)/(x^4+8*x^3+18*x^2+8*x+1))^(1/2)-16*x-4)/(3*x^4+8*x^3+14*x^2+8*x+3) as well as the.. $\endgroup$ – Nasser Sep 23 '14 at 8:54
  • $\begingroup$ ... the one you showed y=-4*x/(1+3*x). The point is, these are singular solutions. There is a singularity in your ode, which is why NDSolve gives a 1/0 error. So, there is nothing wrong with NDSolve. Notice that DSolve did not solve this. $\endgroup$ – Nasser Sep 23 '14 at 8:55
  • $\begingroup$ Thanks, that does clear things up a bit. With the singularity in the ODE, do you mean the fact that at $0$ we have that the ODE becomes $0*y''=0$? Can I conclude that for every value $y''(0)$ there is a different solution? If so, can Mathematica find these? $\endgroup$ – ScroogeMcDuck Sep 23 '14 at 9:56
  • $\begingroup$ @Nasser I'm also curious how you found your second solution $\endgroup$ – ScroogeMcDuck Sep 23 '14 at 16:12
3
$\begingroup$

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

Mathematica graphics

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]

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

Improve this answer
$\endgroup$
  • $\begingroup$ Thanks for the answer. Since it's been so long since I'd have to look up again what this was all about, but I promise to do so in a while and check your answer. $\endgroup$ – ScroogeMcDuck Jun 21 '16 at 11:52
  • $\begingroup$ @ScroogeMcDuck No problem. I was working on something else, searching the site, when I got some help that solved my problem (the "EquationSimplification" -> "Residual" trick). I thought I'd try it out on some the related questions I had found to see how robust it was. Of course, the trick failed here, but in an interesting way. I'm still not fully sure why it failed. You can't increase WorkingPrecision with this trick, and you can't analyze the numerics of internal code you don't have access to. $\endgroup$ – Michael E2 Jun 21 '16 at 12:56
2
$\begingroup$

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

Mathematica graphics

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

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

Mathematica graphics

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:

Mathematica graphics

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

Mathematica graphics

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

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

Mathematica graphics

For example, looking at the first one, change C1

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

Limit[First@sol2, C1 -> 1]

Mathematica graphics

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

Used Maple 18.01 and Mathematica 10.01 both on windows.

Improve this answer
$\endgroup$
  • $\begingroup$ Thanks for the answer, really helps. I'm a bit confused about the following comment: "But if y'(0)=-4, then it y''(0) should be zero. right?" Why is that? E.g. in the solution $y=-4x/(1+3x)$ we have y'(0)=-4 and y''(0)=24. $\endgroup$ – ScroogeMcDuck Sep 24 '14 at 9:04
  • $\begingroup$ Regardless of whether it's 24 or 0 or anything, I don't think my ODE contradicts this. I.e. I don't think it shows that $y''$ should be diverging. Rather, it's indeterminate because with $y'(0)=-4$ one gets $\infty-\infty$. If you replace rhs /. {x -> 0, y[x] -> 0, y'[0] -> -4} in your answer by rhs /. {x -> 0, y[x] -> 0, y'[x] -> -4}, Mathematica also returns Indeterminate, and imposing $y=-4x+ C x^2$ and Tayloring yields $y''(0)=2C$ as it should. Hence I don't think it's true that the initial conditions don't agree with the ODE, correct? $\endgroup$ – ScroogeMcDuck Sep 24 '14 at 9:23

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.