# NDSolve divide-by-zero trouble

So I have a set of differential equations (which I know are analytically integrable) but I want to integrate them numerically to get the form of the solution (the solution is not explicit and hence I want to numerically integrate it).

g=7/5

v = Function[r, 1/(n[r] r^2)]

eulerEq = v[r]*v'[r] + n'[r]*n[r]^(g - 2) + 2/r^2 == 0

initcond = n == 1 (*These are my equation and initial condition *)


Now to numerically solve it.

S = NDSolve[{eulerEq, initcond}, n[r], {r, 0.5, 10}]


When I hit shift+Enter, it shows me the following error:

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

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

NDSolve::ndnum: Encountered non-numerical value for a derivative at r == 1.. >>


I tried using the debugger and I realized that it is computing $n'$ in the process and it is getting a 'division-by-zero'. I can integrate the equation easily and get the constant of integration using the initial condition (it is getting $n$ explicitly in terms of $r$ that is troublesome). Is there a way to do this still using NDSolve or is there an alternative? I want to do it numerically only because this equation is just a special case which is analytically integrable and I need to get the other solutions also.

Determining n'[r] from

Solve[eulerEq, Derivative[n][r]][[1, 1, 2]]
(* -((2 n[r] (-1 + r^3 n[r]^2))/(r (-1 + r^4 n[r]^(12/5)))) *)


we see that it is indeterminate at r = 1, and it is for this reason that NDSolve fails there. Nonetheless, the limiting value can be obtained as follows,

Limit[% /. n[r] -> 1 + a (r - 1), r -> 1]
(* -((5 (3 + 2 a))/(10 + 6 a)) *)
sl = Solve[% == a]
(* {{a -> 1/6 (-10 - Sqrt)}, {a -> 1/6 (-10 + Sqrt)}} *)


and used to define boundary conditions very near r = 1.

dr = .0001;
s1 = NDSolveValue[{eulerEq, n[1 + dr] == 1 + dr a /. First[sl]}, n, {r, 1 + dr, 10}];
s2 = NDSolveValue[{eulerEq, n[1 + dr] == 1 + dr a /. Last[sl]}, n, {r, 1 + dr, 10}];
s3 = NDSolveValue[{eulerEq, n[1 - dr] == 1 - dr a /. First[sl]}, n, {r, .5, 1 - dr}];
s4 = NDSolveValue[{eulerEq, n[1 - dr] == 1 - dr a /. Last[sl]}, n, {r, .5, 1 - dr}];

p1 = Plot[{s1[r], s2[r]}, {r, 1 + dr, 10}, PlotRange -> {{0, 10}, {0, 5}}];
p3 = Plot[{s3[r], s4[r]}, {r, .5, 1 - dr}, PlotRange -> {{0, 10}, {0, 5}}];
Show[p1, p3, AxesLabel -> {r, n}] • good solution you managed to get it. Maple does actually solve this analytically for same initial condition, and I get the same plot as you show. The solution n(r) it gives is in terms of a complicated non-linear function in n(r) – Nasser Dec 18 '15 at 3:51
• @Nasser so Maple uses a different algorithm to compute the numerical integral? Is it possible to implement it here? Because it would be tedious to compute the limit every time such a situation pops up. – Sunil S. Dec 18 '15 at 4:17
• I think you can explain a bit more about why you choose n[r] -> 1 + a (r - 1) but not n[r] -> 1 + a (r - 1)^2, n[r] -> 1 + a (r - 1)^3, etc. :) – xzczd Dec 18 '15 at 4:23
• @SunilS. I do not know how Maple solved it, but here is the solution it gives. There is no numerical anything. it solved the analytical ode. Not numerical. The plots seems similar to what in the above answer: !Mathematica graphics here is plain text of the solution of the ODE itself. 1+5*n(r)^(12/5)*r^4+(-2*r^4-4*r^3)*n(r)^2 = 0 so it is non-linear in n(r) – Nasser Dec 18 '15 at 4:32
• @xzczd Linear was my first guess, and it worked. Once assumed, it is easy to show that the guess is self-consistent, as is apparent from the plot.. – bbgodfrey Dec 18 '15 at 5:06

bbgodfrey has already given a good answer, here's just a not that rigorous but cheap and accurate enough way to resolve the problem:

S = NDSolveValue[{eulerEq, n@1 == 1 + #}, n[r], {r, 0.5, 10},
Method -> "StiffnessSwitching"] & /@ ({1, -1} $MachineEpsilon) Plot[S, {r, 0.5, 10}, PlotRange -> All] • Why does {1.7, -1.7} work as the perturbation, when {1.8, -1.8} does not? Perhaps, the boundary of the basins of attraction. In any case, +1. – bbgodfrey Dec 18 '15 at 5:02 • @bbgodfrey Interesting observation. I guess it might be related to how "StiffnessSwitching" works (try Manipulate[With[{S = NDSolveValue[{eulerEq, n@1 == 1 + #}, n[r], {r, 0.5, 10} , Method -> "StiffnessSwitching"] & /@ ({1.7, -a}$MachineEpsilon)}, Plot[S, {r, 0.5, 10}, PlotRange -> {0, 5}, Epilog -> {PointSize@Medium, Point[{1, 1}]}]], {a, 1, 10}]), but currently I can't find a better Method to handle the singularity. – xzczd Dec 18 '15 at 5:58

Using the limiting values from bbgodfrey's answer, we can take a more direct approach with NDSolve using "EquationSimplification" -> "Residual":

sol = Join @@@
(NDSolve[{eulerEq, initcond, n' == #},
n[r], {r, 0.5, 10},
Method -> {"EquationSimplification" -> "Residual"}] & /@
{1/6 (-10 + Sqrt), 1/6 (-10 - Sqrt)});

Plot[n[r] /. sol // Evaluate, {r, 0.5, 10}, PlotRange -> All,
AxesOrigin -> {0, 0}]
` 