2
$\begingroup$

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?

Improve this question
$\endgroup$
7
  • 2
    $\begingroup$ 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). $\endgroup$
    – Michael E2
    May 17, 2013 at 21:28
  • $\begingroup$ Oh, yes, this makes sense. Thank you! $\endgroup$
    – Kagaratsch
    May 17, 2013 at 21:44
  • $\begingroup$ 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? $\endgroup$
    – Kagaratsch
    May 17, 2013 at 21:53
  • $\begingroup$ BTW: The relevant information I want to acquire from the computation is the value of the solution at the right border y[1] $\endgroup$
    – Kagaratsch
    May 17, 2013 at 21:54
  • $\begingroup$ 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? $\endgroup$
    – Michael E2
    May 17, 2013 at 22:03

2 Answers 2

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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... $\endgroup$
    – Kagaratsch
    May 17, 2013 at 22:11
  • 1
    $\begingroup$ 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$. $\endgroup$
    – Andrew
    May 18, 2013 at 8:20
  • $\begingroup$ 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; $\endgroup$
    – Kagaratsch
    May 18, 2013 at 9:05
  • $\begingroup$ @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? $\endgroup$
    – Andrew
    May 18, 2013 at 10:37
1
$\begingroup$

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]]
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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