Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
2  
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). –  Michael E2 May 17 '13 at 21:28
    
Oh, yes, this makes sense. Thank you! –  Kagaratsch May 17 '13 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? –  Kagaratsch May 17 '13 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] –  Kagaratsch May 17 '13 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? –  Michael E2 May 17 '13 at 22:03
show 2 more comments

2 Answers

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

Output

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.

share|improve this answer
    
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... –  Kagaratsch May 17 '13 at 22:11
1  
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$. –  Andrew May 18 '13 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; –  Kagaratsch May 18 '13 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? –  Andrew May 18 '13 at 10:37
add comment

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]]
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.