Strange Behavior of NDSolve

I am trying to evaluate the following ODE numerically:

NDSolve[(2 + Sqrt[2] +s^2 -(-2 + Sqrt[2]) s^4 + s^6) /(1 + s^2)^2 yy[s] -(1 - s^2)^2 yy''[s] == 0 &&
yy[-1] == 0 &&   yy'[-1] == 1, yy[s], {s, -1, 1}]


and I surprisingly get the following error:

Power::infy: "Infinite expression 1/0.^2 encountered."
Infinity::indet: "Indeterminate expression 0.\ ComplexInfinity encountered."
NDSolve::ndnum: "Encountered non-numerical value for a derivative at s == -1.."


One can verify directly that no infinities should appear at s=-1:

D[(2 + Sqrt[2] + s^2 - (-2 + Sqrt[2]) s^4 + s^6)/(1 + s^2)^2, s] /.  s -> -1

gives: 3 + 1/4 (-8 + 4 (-2 + Sqrt[2]))


But apparently, the ODE solver automatically divides by the factor in front of the yy''[s] term before computation. in this case we indeed find a singular

D[-(2 + Sqrt[2] + s^2 - (-2 + Sqrt[2]) s^4 + s^6)/((1 + s^2)^2 (1 - s^2)^2), s] /. s -> -1


Clearly, if one does not divide by the factor of the kinetic term, the ODE is well defined on the whole interval including the boundaries. Why does the solver divide? What can I do to prevent this? Or maybe there is some workaround?

• It cannot solve for yy''[s] at s == -1, since the coefficient vanishes. (At a naive level, it cannot then approximate the next value of yy'[s] for s a little greater than -1, etc.) What is being detected is that s == -1 is a singular point, and there are infinitely many solutions to this IVP (one for each value of yy''[-1] you care to assign). Commented May 17, 2013 at 21:28
• Oh, yes, this makes sense. Thank you! Commented May 17, 2013 at 21:44
• So, now that the problem with the ODE has been named, what would be an appropriate workaround? I can think of two: 1) Andrew's answer below, 2) Introduce a regulator eps=10^-14 as a summand to the factor in front of yy''[s]. Are these any good? Commented May 17, 2013 at 21:53
• BTW: The relevant information I want to acquire from the computation is the value of the solution at the right border y[1] Commented May 17, 2013 at 21:54
• Do you know the values of y[s], y'[s] at s == 0 or some other point? Or is knowing that y approaches infinity at s == 1 enough? Commented May 17, 2013 at 22:03

As an illustration of my initial comment, let's look at various starting points, following @Andrew's method:

sols = Table[
NDSolveValue[(2 + Sqrt[2] + s^2 - (-2 + Sqrt[2]) s^4 + s^6)/(1 + s^2)^2 yy[s] -
(1 - s^2)^2 yy''[s] == 0 && yy[start] == yy'[start] (start + 1) && yy'[start] == 1,
yy[s], {s, start, 1}], {start, -1 + 10^-9, -1 + 2 10^-8, 10^-9}];
Plot[sols, {s, -1, 1}]


Expanded comment

I may have been hasty in asserting there are infinitely many solutions to the IVP in the question, although it looks plausible. The point s == -1 is a a regular singular point whose indicial equation has complex roots (if I didn't make a mistake). Unfortunately I do not know much about such situations. Certainly it should be clear that one cannot pick a particular starting point as a workaround.

• Similarly, sols = Table[ NDSolveValue[(2 + Sqrt[2] + s^2 - (-2 + Sqrt[2]) s^4 + s^6)/(1 + s^2)^2 yy[s] - ((1 - s^2)^2 + shift) yy''[s] == 0 && yy[-1] == 0 && yy'[-1] == 1, yy[s], {s, -1, 1}], {shift, 10^-9, 2 10^-8, 10^-9}]; Plot[sols, {s, 1 - 0.000001, 1}] also displays a shift... Commented May 17, 2013 at 22:11
• The problem here can be that the coefficient at $y''$ vanishes at $s=1$ too. And looking at the plot it seems not too implausible that $\lim_{s\to1-0}y(s)=+\infty$. Also Mma gives a warning that the solution may be stiff near $s=1$. Commented May 18, 2013 at 8:20
• Yes, considering that the potential is symmetric it is clear that one should regularize both ends. e.g. leftend = -1 + 10^-8; and rightend = 1 - 10^-8; Commented May 18, 2013 at 9:05
• @Kagaratsch But as far as I understood from the comments to the question the value of $y(1)$ is required. If it is equal to $+\infty$, what else information of solution is needed? Commented May 18, 2013 at 10:37

As a workaround it is possible to step up a bit from the left end:

leftend = -1 + 10^-8;
NDSolve[(2 + Sqrt[2] + s^2 - (-2 + Sqrt[2]) s^4 + s^6)/(1 + s^2)^2 yy[s]
- (1 - s^2)^2 yy''[s] == 0 && yy[leftend] == 0 && yy'[leftend] == 1, yy[s],
{s, leftend, 1}][[1]]
`