# NDSolve fails at the regular singular point of a second-order ODE

Here is the ODE I want to numerically integrate,

odey=-l (1 + l) R[y] + (k - y) (-2 Derivative[R][y] + (k - y) (R^\[Prime]\[Prime])[y])


If we rearrange it in the standard form,

$$R''(y)-\frac{2}{k-y}R'(y)-\frac{l(l+1)}{(k-y)^{2}}R(y)=0$$

we see that it has a regular singular point at $$y=k$$ where $$k<0$$. I also attempted to solve this using DSolve, and it gave me a solution that conforms with the form I expect, that is $$R(y)=c_{1}(y-k)^{a}+c_{2}(y-k)^{b}$$

DSolve[odey == 0, R[y], y]


Now, I want $$R(y)$$ to vanish for very large $$y$$ (boundary condition). So I make the transformation $$t=1/y$$ (where $$y=\infty$$ corresponds to $$t=0$$),

odet = odey /. R -> (R[1/#] &) /. y -> 1/t // Simplify


What I do next is to perform Frobenius method (series solution of $$R$$)

odetSeries = Series[odet /. R -> Function[{t}, t^p Sum[a[i] t^i, {i, 0, 7}]] // Simplify, {t, 0, 4}] // Simplify


where the indicial equation to solve is

Solve[-(l + l^2 + p - p^2) a == 0, p]


I chose $$p=1+l$$ to satisfy the boundary condition. After which, I obtained the first few coefficients of the expansion $$a_{0},a_{1},a_{2}$$

odetSeries /. p -> 1 + l /. a -> 0 // Simplify
Solve[(k^2 (2 + 3 l + l^2) a + 2 (3 + 2 l) a) == 0, a] //Flatten[#] & // Simplify
odetSeries /. p -> 1 + l /. a -> 0 /. % /. a -> 0 // Simplify
Solve[(-k^4 (24 + 50 l + 35 l^2 + 10 l^3 + l^4) a + 8 (15 + 16 l + 4 l^2) a) == 0, a] // Flatten[#] & // Simplify

a4[l_] := (k^4 (24 + 50 l + 35 l^2 + 10 l^3 + l^4) a0)/(8 (15 + 16 l + 4 l^2));
a2[l_] := -((k^2 (2 + 3 l + l^2) a0)/(6 + 4 l));

solInfPlus[l_] := t^(1 + l) (a0 + a2[l] t^2 + a4[l] t^4);
testsol = solInfPlus[l] /. t -> 1/y
Collect[testsol, y, Simplify] /. a0 -> 1
D[testsol, y] /. a0 -> 1 // Collect[#, y, Simplify] &

R[y_, l_] := (1/y)^(1 + l) - (k^2 (2 + 3 l + l^2) (1/y)^(3 + l))/(6 + 4 l) + (k^4 (24 + 50 l + 35 l^2 + 10 l^3 + l^4) (1/y)^(5 + l))/(8 (15 + 16 l + 4 l^2))
dR[y_, l_] := -(1 + l) (1/y)^(2 + l) + (k^2 (6 + 11 l + 6 l^2 + l^3) (1/y)^(4 + l))/(6 + 4 l) - (k^4 (120 + 274 l + 225 l^2 + 85 l^3 + 15 l^4 + l^5) (1/y)^(6 + l))/(8 (15 + 16 l + 4 l^2))


Now, I execute numerical calculations.

k=(Sqrt[\[Pi]] Gamma[-(1/3)])/(3 Gamma[1/6]);

rules = {AccuracyGoal -> Infinity, PrecisionGoal -> 20, WorkingPrecision -> 30, MaxSteps -> 10000};
rat = 10^-30;
yP = 10^3;
yM = -10^3;
y0 = 10^-2;

For[el = 0, el <= 8, el++, {
R0p = Rationalize[R[yP, el], rat];
dR0p = Rationalize[dR[yP, el], rat];
R0m = Rationalize[R[yM, el], rat];
dR0m = Rationalize[dR[yM, el], rat];

Rsolp = R /. First@NDSolve[{(odey /. l -> el) == 0, R[yP] == R0p, R'[yP] == dR0p}, {R}, {y, yP, y0}, rules];
Rsolm = R /. First@NDSolve[{(odey /. l -> el) == 0, R[yM] == R0m, R'[yM] == dR0m}, {R}, {y, yM, y0}, rules];

rp = Rationalize[ Rsolp[y0], rat];
rm = Rationalize[ Rsolm[y0], rat];
drp = Rationalize[ Rsolp'[y0], rat];
drm = Rationalize[ Rsolm'[y0], rat];

s = Rationalize[(-Q a Sqrt[4 Pi (2 el + 1)])/(a(y0-k))^2,rat];
c1f = (rm s)/(drp rm - drm rp);(*jump condition at y=y0 between Rsolp and Rsolm*)
baremode[el] = Sqrt[(2 el + 1)/(4 \[Pi])] (c1f Rsolp'[y0]);}]


But the solver stops evaluating exactly at $$y=k=-0.431188993916544683366628968310$$, this is the regular singular point of the ODE above. Specifically, the NDSolve fails in Rsolm since Rsolm covers the region from $$yM=-10^{3}$$ to $$y0=10^{-2}$$ in which the singularity $$y=k=-0.431188993916544683366628968310$$ belongs.

Is there a way to solve this problem? I am thinking of integrating near the singular point, bypass it and then continue integrating. But I do not have an idea on how to execute it properly.

• What equation do you want to solve - odey or  in the standard form? At first glance, these are different equations. Or did you make a mistake? Feb 2, 2019 at 13:53
• Hi @AlexTrounev, they are just the same. But you can just consider odey since it's a homogeneous ODE Feb 2, 2019 at 13:56
• Watch carefully compare degrees $(k-y)$. Feb 2, 2019 at 14:00
• Oh, sorry, I didn't see the brackets. Feb 2, 2019 at 14:09
• @AlexTrounev, I am thinking of a sort of numerical continuation methods. I guess it will help? Or, a good coordinate transformation so that the ode is free from singular points. I don't know exactly. I hope you can help me with this Feb 2, 2019 at 14:12

We write the equation in the form

odey = -l (1 + l) R[y] - 2*(k - y)*R'[y] + (k - y)^2*R''[y]


Make a replacement coefficient at the second derivative $$(k-y)^2\to(k-y)^2+h^2$$:

 odey = -l (1 + l) R[y] - 2*(k - y)*R'[y] + ((k - y)^2+h^2)*R''[y]


Choose h of the order of the minimum integration step.

Example

odey = -l (1 + l) R1[y] -
2*(k - y)*R1'[y] + ((k - y)^2 + h^2)*R1''[y] == 0;

k = (Sqrt[π] Gamma[-(1/3)])/(3 Gamma[1/6]); h = 10^-5;
R[y_, l_] := (1/y)^(1 + l) - (k^2 (2 + 3 l + l^2) (1/y)^(3 + l))/(6 +
4 l) + (k^4 (24 + 50 l + 35 l^2 + 10 l^3 +
l^4) (1/y)^(5 + l))/(8 (15 + 16 l + 4 l^2))
dR[y_, l_] := -(1 + l) (1/y)^(2 +
l) + (k^2 (6 + 11 l + 6 l^2 + l^3) (1/y)^(4 + l))/(6 +
4 l) - (k^4 (120 + 274 l + 225 l^2 + 85 l^3 + 15 l^4 +
l^5) (1/y)^(6 + l))/(8 (15 + 16 l + 4 l^2))
yP = 10^3;
yM = -10^3;
y0 = 10^-2;
a = 1; Q = 1; l = 4; el = l;
R0p = R[yP, el];
dR0p = dR[yP, el];
R0m = R[yM, el];
dR0m = dR[yM, el];

s1 = NDSolve[{odey, R1[yP] == R0p, R1'[yP] == dR0p}, R1, {y, y0, yP}]

s2 = NDSolve[{odey, R1[yM] == R0m, R1'[yM] == dR0m}, R1, {y, yM, y0}]

{Plot[R1[y] /. s1, {y, y0, yP}, PlotRange -> All,
PlotLabel -> Row[{"l=", l}]],
Plot[R1[y] /. s2, {y, yM, y0}, PlotRange -> All]} An example of using an exact solution

rl[ρ_, b0_,
q_] := ρ (1 - (Sqrt[π] Gamma[1/(q - 1)])/((1 - q) Gamma[
1/2 ((q + 1)/(q - 1))]) b0/ρ)
r1[ρ_] := rl[ρ, a1, -2];
f0 = 1 - a1^3/r1[a1 y0]^3;
R[y_, l_] := (1/y)^(1 +
l) - ((2 + 3 l + l^2) π (1/y)^(3 + l) Gamma[2/3]^2)/(2 (3 +
2 l) Gamma[1/6]^2) + ((24 + 50 l + 35 l^2 + 10 l^3 +
l^4) π^2 (1/y)^(5 + l) Gamma[2/3]^4)/(8 (15 + 16 l +
4 l^2) Gamma[1/6]^4);
dR[y_, l_] := -(1 + l) (1/y)^(2 +
l) + ((6 + 11 l + 6 l^2 + l^3) π (1/y)^(4 + l) Gamma[
2/3]^2)/(2 (3 + 2 l) Gamma[1/6]^2) - ((120 + 274 l +
225 l^2 + 85 l^3 + 15 l^4 + l^5) π^2 (1/y)^(6 + l) Gamma[
2/3]^4)/(8 (15 + 16 l + 4 l^2) Gamma[1/6]^4);
y0 = 10^-3;
yP = 10^3;
yM = -10^3;
a1 = 1;
b = a1 Sqrt[π] Gamma[2/3]/Gamma[1/6];
a = (a1 Sqrt[π] Gamma[2/3])/Gamma[1/6];
Q = 1;
R0p = R[yP, el];
dR0p = dR[yP, el];
R0m = R[yM, el];
dR0m = dR[yM, el];
sol = DSolve[-l (1 + l) R1[ρ] +
2 (ρ + b) R1'[ρ] + (ρ + a)^2*R1''[ρ] == 0,
R1, ρ];
Rs = First[ R1[ρ] /. sol] /. ρ -> y;
eq1 = Rs /. {y -> yP};
eq2 = D[Rs, y] /. {y -> yP};
sC1 = Solve[eq1 == R0p && eq2 == dR0p, {C, C}];
Rs1 = Rs /. sC1;
eq3 = Rs /. {y -> yM};
eq4 = D[Rs, y] /. {y -> yM};
sC2 = Solve[eq3 == R0m && eq2 == dR0m, {C, C}];
Rs2 = Rs /. sC2;
Rsolp = Rs1 /. l -> el;
Rsolm = Rs2 /. l -> el;
rp = Rsolp /. y -> y0;
rm = Rsolm /. y -> y0;
drp = D[Rsolp, y] /. y -> y0;
drm = D[Rsolm, y] /. y -> y0;
s = (-Q a1 Sqrt[4 Pi (2 el + 1)])/(rl[a1 y0, a1, -2])^2;
c1f = (rm s)/(drp rm - drm rp); baremode =
Sqrt[(2 el + 1)/(4 π)] (c1f D[Rsolp, y] /. y -> y0);

Table[{Plot[Rsolp /. el -> i, {y, y0, yP},
PlotLabel -> Row[{"l=", i}]],
Plot[Rsolm /. el -> i, {y, y0, yM}]}, {i, 1, 4}] The case $$l=0$$ is degenerate, therefore the exact solution is determined differently. To combine the two cases, divide the code into two blocks, and then combine them using If[]:

rl[\[Rho]_, b0_,
q_] := \[Rho] (1 - (Sqrt[\[Pi]] Gamma[1/(q - 1)])/((1 - q) Gamma[
1/2 ((q + 1)/(q - 1))]) b0/\[Rho])
r1[\[Rho]_] := rl[\[Rho], a1, -2];
f0 = 1 - a1^3/r1[a1 y0]^3;
R[y_, l_] := (1/y)^(1 +
l) - ((2 + 3 l + l^2) \[Pi] (1/y)^(3 + l) Gamma[2/3]^2)/(2 (3 +
2 l) Gamma[1/6]^2) + ((24 + 50 l + 35 l^2 + 10 l^3 +
l^4) \[Pi]^2 (1/y)^(5 + l) Gamma[2/3]^4)/(8 (15 + 16 l +
4 l^2) Gamma[1/6]^4);
dR[y_, l_] := -(1 + l) (1/y)^(2 +
l) + ((6 + 11 l + 6 l^2 + l^3) \[Pi] (1/y)^(4 + l) Gamma[
2/3]^2)/(2 (3 + 2 l) Gamma[1/6]^2) - ((120 + 274 l +
225 l^2 + 85 l^3 + 15 l^4 + l^5) \[Pi]^2 (1/y)^(6 + l) Gamma[
2/3]^4)/(8 (15 + 16 l + 4 l^2) Gamma[1/6]^4);
y0 = 10^-3;
yP = 10^3;
yM = -10^3;
a1 = 1;
b = a1 Sqrt[\[Pi]] Gamma[2/3]/Gamma[1/6];
a = (a1 Sqrt[\[Pi]] Gamma[2/3])/Gamma[1/6];
Q = 1;
Baremode[L_] := Block[{el = L}, R0p = R[yP, el];
dR0p = dR[yP, el];
R0m = R[yM, el];
dR0m = dR[yM, el];
sol = DSolve[-l (1 + l) R1[\[Rho]] +
2 (\[Rho] + b) R1'[\[Rho]] + (\[Rho] + a)^2*R1''[\[Rho]] == 0,
R1, \[Rho]];
Rs = First[ R1[\[Rho]] /. sol] /. \[Rho] -> y;
eq1 = Rs /. {y -> yP};
eq2 = D[Rs, y] /. {y -> yP};
sC1 = Solve[eq1 == R0p && eq2 == dR0p, {C, C}];
Rs1 = Rs /. sC1;
eq3 = Rs /. {y -> yM};
eq4 = D[Rs, y] /. {y -> yM};
sC2 = Solve[eq3 == R0m && eq2 == dR0m, {C, C}];
Rs2 = Rs /. sC2;
Rsolp = Rs1 /. l -> el;
Rsolm = Rs2 /. l -> el;
rp = Rsolp /. y -> y0;
rm = Rsolm /. y -> y0;
drp = D[Rsolp, y] /. y -> y0;
drm = D[Rsolm, y] /. y -> y0;
s = (-Q a1 Sqrt[4 Pi (2 el + 1)])/(rl[a1 y0, a1, -2])^2;
c1f = (rm s)/(drp rm - drm rp);
baremode = Sqrt[(2 el + 1)/(4 \[Pi])] (c1f D[Rsolp, y] /. y -> y0);
Re[baremode]];
Baremode0[L0_] := Block[{el = L0},
R0p = R[yP, el];
dR0p = dR[yP, el];
R0m = R[yM, el];
dR0m = dR[yM, el];
sol = DSolve[
2 (\[Rho] + b) R1'[\[Rho]] + (\[Rho] + a)^2*R1''[\[Rho]] == 0,
R1, \[Rho]];
Rs = First[ R1[\[Rho]] /. sol] /. \[Rho] -> y;
eq1 = Rs /. {y -> yP};
eq2 = D[Rs, y] /. {y -> yP};
sC1 = Solve[eq1 == R0p && eq2 == dR0p, {C, C}];
Rsolp = Rs /. sC1;
eq3 = Rs /. {y -> yM};
eq4 = D[Rs, y] /. {y -> yM};
sC2 = Solve[eq3 == R0m && eq2 == dR0m, {C, C}];
Rsolm = Rs /. sC2;
rp = Rsolp /. y -> y0;
rm = Rsolm /. y -> y0;
drp = D[Rsolp, y] /. y -> y0;
drm = D[Rsolm, y] /. y -> y0;
s = (-Q a1 Sqrt[4 Pi ])/(rl[a1 y0, a1, -2])^2;
c1f = (rm s)/(drp rm - drm rp);
baremode = Sqrt[1/(4 \[Pi])] (c1f D[Rsolp, y] /. y -> y0);
Re[baremode]];


Use the code with If[] for printing baremode (we dropped the zero imaginary part):

For[k = 0, k <= 4, k++,
Print[{Row[{"l=", k}],
Row[{"baremode=",
First[If[k == 0, N[Baremode0, 6], N[Baremode[k], 6]]]}]}]]
{l=0,baremode=-7.18650*10^6}

{l=1,baremode=-9.94531*10^10}

{l=2,baremode=-1.14027*10^20}

{l=3,baremode=-6.65090*10^24}

{l=4,baremode=-2.09589*10^33}

• What do you mean by 'of the order of the minimum integration step'? And what does the replacement coefficient do? Feb 2, 2019 at 14:50
• Put h=10^-5;, it is enough to get through y=k. Feb 2, 2019 at 14:58
• Ok. I'll try it. I will let you know what it will end up. Thanks @Alex Feb 2, 2019 at 15:02
• NDSolve still fails NDSolve::mxst: Maximum number of 10000 steps reached at the point y == -0.431183866743419284299032955175. Feb 2, 2019 at 15:07
• There is another problem. You try to calculate large numbers of order Exp[10^3] Feb 2, 2019 at 15:59

The solution of your ode is (l >= 0):

odey /. R -> Function[y, (c1 (k - y)^l + c2/(k - y)^(1 + l)) ]
(* 0 *)


The special solution Limit[R[y],y->Infinity]==0 requires c1->0:

 R~1/(k - y)^(1 + l)

• Right, so what we do with $y=k$? It blows up Feb 2, 2019 at 15:43
• Yes, and the asymptotic solution shows you how the singularity is "passed". Feb 2, 2019 at 15:56
• I don't quite understand what you really mean Feb 2, 2019 at 15:58
• The solution of odey shows a part vanishing for y->infinity (c1=0,c2!=0) and a part without singularity (c2==0,c1!=0) but increasing with y! That means the only solution without singularity, vanishing with y-> infinity is R==0 Feb 2, 2019 at 16:01
• I get that. But how do I actually evaluate $R$ at $y=k$? I just use $R\sim(k-y)^{l}$? Feb 2, 2019 at 16:05

I'm not sure what you're trying to do, but if you simply want to skip past the singularity in a numerical solver, you can deform your integration path by walking around the singularity in the complex plane. This will work for $$l\in\mathbb N_0$$ since the singularity will be a pole, so the continuation is unique. Try e.g. parametrizing your $$y$$ as

$$y(t)=\begin{cases}t+i\cos\left(\frac\pi{2\epsilon}(k-t)\right)^3 & \text{if } k-\epsilon\le t\le k+\epsilon, \\ t & \text{otherwise}. \end{cases}$$

Then express $$R'(y)$$ and $$R''(y)$$ in terms of another function $$g(t)=R(y(t))$$ and use NDSolve to find $$g$$. After this you can get $$R(y)$$ back (in our parametrization it's trivial for $$y\not\in(k-\epsilon,k+\epsilon)$$).

Example:

y[t_] = Piecewise[{{t + I*Cos[(Pi/(2*ϵ))*(k - t)]^3, k - ϵ <= t <= k + ϵ}, {t, True}}];
yp[t_] = y'[t] // Simplify;
ypp[t_] = y''[t] // Simplify;

k = (Sqrt[π] Gamma[-(1/3)])/(3 Gamma[1/6]);
l = 3;
ϵ = 0.1;

sol[t_] = NDSolveValue[{1/yp[t]^2 (g''[t] - ypp[t]/yp[t] g'[t]) -
2/(k - y[t])g'[t]/yp[t]-(l(l+1))/(k - y[t])^2 g[t] == 0,
g == 0, g' == -1}, g[t], {t, -1, 10}, WorkingPrecision -> 30];

(* Make sure that the solution is proportional to exact solution *)
Plot[sol[t]/(k - y[t])^(-1 - l), {t, -1, 1}, PlotRange -> All, AxesOrigin -> {0, 0}] • Hi @Ruslan, how can I actually retrieve $R(y)$ from the result of NDSolve for $g(t)$? Feb 4, 2019 at 10:36
• @user583893 by construction of y[t], $R(y)=g(y)=g(t)\;\forall y\not\in(k-\epsilon,k+\epsilon)$. Just make sure $\epsilon$ is small enough that you don't care about the excluded interval. Feb 4, 2019 at 11:22
• Can you update your answer including the actual code retrieving $R(y)$ @Ruslan? Thank you Feb 4, 2019 at 14:53
• @user583893 it's trivial: R[Y_]=Piecewise[{{Indeterminate, k-ϵ<=Y<=k+ϵ}, {sol[Y], True}}] Feb 4, 2019 at 16:51
• Thanks. By the way, what is the motivation for the form of the parametrization? Feb 4, 2019 at 16:55