# How to resolve the singularity in NDSolve

I ran into a problem in my Mathematica code. There is a singularity or stiff system in NDSolve. The explicit code is as follows:

zslist = Table[9995/10000 - i/480 9995/10000, {i, 0, 479}];
ϵ := 1/10^6;
zc = 10^-2;
d = 4;
Do[z0[i_, d_] :=
zslist[[i]] - ((1 - zslist[[i]]^d) ϵ^2)/(2 zslist[[i]]);
z1[i_, d_] := -(((1 - zslist[[i]]^d) ϵ)/zslist[[i]]);
s[i, d] =
NDSolve[{z''[ρ] == -((2 z[ρ]^3
z'[ρ]^2)/(1 - z[ρ]^4)) + (-(2/z[ρ]) -
z'[ρ]/(ρ (1 - z[ρ]^4))) (1 -
z[ρ]^4 + z'[ρ]^2),
z[ϵ] == z0[i, d], z'[ϵ] == z1[i, d]},
z, {ρ, ϵ, 10}];
f[ρ_, i_, d_] := s[i, d][[1, 1]][[2]][ρ];
r[i_, d_] :=
FindRoot[f[ρ, i, d] == zc, {ρ,
s[i, d][[1, 1]][[2]][[1, 1, 2]]}][[1, 2]];,
{i, 1, Length[zslist]}]


I don't know how to resolve this puzzle. Could you tell me some methods to make it work?

• While you can use HTML tags like <code> formatting, I changed it to the standard Markdown method to fix a couple little things. (The two methods seem to behave slightly differently, and I knew how to do what I wanted in Markdown. You format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers.) – Michael E2 Jul 26 '15 at 11:06
• And welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Michael E2 Jul 26 '15 at 11:06
• There is something about your code that is unclear to me. You sometimes use a pattern i_ as in z0[i_, d_] := ... and sometimes you use the variable i as in s[i, d] = .... These mean very different things. For instance, the definition z0[i_, d_] := ... is the exactly the same at each iteration of the Do loop; instead, it should be placed outside the loop. If you need further explanation, please say so. Someone can help -- there may be a Q&A on site that explains the difference. – Michael E2 Jul 26 '15 at 11:13
• Thank you very much! Let me modify the errors. However I don't think that this is the origin of the problem in my code. – amon xu Jul 26 '15 at 11:16
• I agree. I'd be surprised if that had something to do with the NDSolve::ndsz error. – Michael E2 Jul 26 '15 at 11:24

I think if you examine the solutions, you will see that z'[t] -> - Infinity near the point where the integration ends.

Manipulate[
With[{ρminmax = Flatten[z["Domain"] /. s[i, d]]}, (* start/stop values *)
Plot[
{z[ρ], z'[ρ]} /. s[i, d] // Flatten // Evaluate,  (* fn. & deriv. *)
{ρ, ρminmax[[1]], ρminmax[[2]]},
PlotLabel ->
Row[{z'[Subscript[ρ, final]], " = ", z'[ρminmax[[2]]] /. First@s[i, d]}]
]],
{i, 1, Length[zslist], 1}
]


This behavior persists even if the option Method -> "StiffnessSwitching" is used.

My first thought at this point is that system is bound to run into a singularity. Is there some reason to think the system can be integrated past this singularity, say if z were like ρ^(1/3)?

• I think that the singularity always exist, for example we can choose $zslist = Table[9995/10000 - i/480 9995/10000, {i, 0, 99}];$ – amon xu Jul 26 '15 at 12:08
• @amonxu So the result makes sense to you? Note that the NDSolve::ndsz message is not really an "error"; it is a warning that the integration was successful only up to a point and not all the way to the end of the interval {ρ, ϵ, 10}. – Michael E2 Jul 26 '15 at 12:12
• Maybe, because I want to replicate figure.7 in the article 1103.2683 which is plot in the interval '{ρ,0,2.6} ' . The code is complicated, if we get a smooth solution, we need to use it integrating numerically. In the step, there are also warning messages, which leads to a wrong figure which contains many divergence points near the singularity. So I am still not sure that this message isn't harmless to our final goal. – amon xu Jul 26 '15 at 12:21
• @amonxu I'm suggesting only that Mathematica has calculated accurately the system as it is set up in the code; but perhaps the way it is set up needs some adjustment to bring the results closer to the paper? – Michael E2 Jul 26 '15 at 13:38
• I am grateful for your help. I hope you can help me resolve the problems in the code. I wrote emails to one of the authors, he told me the code is almost the same as theirs. He tell me that my code need tuning the workprecison and there are some noises in the data. – amon xu Jul 26 '15 at 14:00