Here is the ODE I want to numerically integrate,

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

If we rearrange it in the standard form,


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 the relevant solution (and its derivative) is

R[y_, l_] := 1/(y - k)^((1 + Sqrt[1 + 4 l (l + 1)])/2)
dR[y_, l_] := -(1/2) (1 + Sqrt[(1 + 2 l)^2]) (-k + y)^(-(3/2) - 1/2 Sqrt[(1 + 2 l)^2])

K[q_] := (Sqrt[\[Pi]] Gamma[1/(q - 1)])/((1 - q) Gamma[1/2 ((q + 1)/(q - 1))])
k = K[-2];

where I just set the constant to 1. Now, I get my initial condition from R itself and then numerically solve the ODE across the region $yM$ to $y0$,

rules = {AccuracyGoal -> Infinity, PrecisionGoal -> 20, WorkingPrecision -> 200, MaxSteps -> 10000};
yM = - 10^3;
y0 = 1; 
R0m = Rationalize[R[yM, el], rat];
dR0m = Rationalize[dR[yM, el], rat];
BCm = {R[yM] == R0m, R'[yM] == dR0m};
EQm = {(odey /. l -> el) == 0};
Rsolm = R /. First@NDSolve[Union[EQm, BCm], {R}, {y, yM, y0}, rules, Method -> "StiffnessSwitching"];

With the range of integration of the ODE, the NDSolve surely fails since it passes through the singular point at $y=k$. I want to solve the ODE for all $l=(0,30)$.

My question is that, is there a convenient numerical continuation method to bypass the singularity so that NDSolve does not fail?


1 Answer 1


We can take the approach in Extending NDSolve beyond a singularity and bump it up to order 2. The approach is to solve the ODE in projective space, with $R(z) = [p(x):q(z)]$ in projective coordinates. The real affine solution is then given by $R(z) = p(z)/q(z)$. Since we're adding another variable, we have to add an equation here and there. They are all based on some normalization of the projective coordinates. In this case, we chose $$p(z)^2 + q(z)^2 = \text{constant}\,.$$ By differentiating this relation, we can obtain the extra equations that we need.

projODE = {
    EQm[[1, 1]] /. R -> (p[#]/q[#] &) // Together // Numerator,
    D[p[y]^2 + q[y]^2, {y, 2}]} == 0 // Thread;
projICS = BCm /. {
    R[y0_] == r0_ :> {p[y0], q[y0]} == {Numerator[r0], Denominator[r0]},
    R'[y0_] == rp0_ :> 
     ({D[p[y]/q[y], y] == rp0, D[p[y]^2 + q[y]^2, y] == 0} /. y -> y0)};

projSOL = NDSolve[{projODE, projICS}, {p, q}, {y, yM, y0},
   rules, Method -> "StiffnessSwitching"];

The solution p[y]/q[y] tracks Rsolm[y] up to the singularity, and then continues on past it.

Plot[{Rsolm[y], p[y]/q[y] /. First@projSOL}, {y, yM, y0}, 
 PlotStyle -> {AbsoluteThickness[4], AbsoluteThickness[2]},
 PlotLegends -> {R, p/q}]

enter image description here

Plot[{Rsolm[y], p[y]/q[y] /. First@projSOL}, {y, -1, y0}, 
 PlotStyle -> {AbsoluteThickness[4], AbsoluteThickness[2]}, 
 Exclusions -> (q[y] == 0 /. First@projSOL), PlotLegends -> {R, p/q}]

enter image description here


Your Answer

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

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