2
$\begingroup$

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,

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

$\endgroup$
10
  • $\begingroup$ 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? $\endgroup$ Feb 2, 2019 at 13:53
  • $\begingroup$ Hi @AlexTrounev, they are just the same. But you can just consider odey since it's a homogeneous ODE $\endgroup$
    – user583893
    Feb 2, 2019 at 13:56
  • $\begingroup$ Watch carefully compare degrees $(k-y)$. $\endgroup$ Feb 2, 2019 at 14:00
  • $\begingroup$ Oh, sorry, I didn't see the brackets. $\endgroup$ Feb 2, 2019 at 14:09
  • $\begingroup$ @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 $\endgroup$
    – user583893
    Feb 2, 2019 at 14:12

3 Answers 3

2
$\begingroup$

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

fig1

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[1], C[2]}];
Rs1 = Rs /. sC1;
eq3 = Rs /. {y -> yM};
eq4 = D[Rs, y] /. {y -> yM};
sC2 = Solve[eq3 == R0m && eq2 == dR0m, {C[1], C[2]}];
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}]

fig2

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[1], C[2]}];
   Rs1 = Rs /. sC1;
   eq3 = Rs /. {y -> yM};
   eq4 = D[Rs, y] /. {y -> yM};
   sC2 = Solve[eq3 == R0m && eq2 == dR0m, {C[1], C[2]}];
   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[1], C[2]}];
   Rsolp = Rs /. sC1;
   eq3 = Rs /. {y -> yM};
   eq4 = D[Rs, y] /. {y -> yM};
   sC2 = Solve[eq3 == R0m && eq2 == dR0m, {C[1], C[2]}];
   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[0], 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}
$\endgroup$
14
  • $\begingroup$ What do you mean by 'of the order of the minimum integration step'? And what does the replacement coefficient do? $\endgroup$
    – user583893
    Feb 2, 2019 at 14:50
  • $\begingroup$ Put h=10^-5;, it is enough to get through y=k. $\endgroup$ Feb 2, 2019 at 14:58
  • $\begingroup$ Ok. I'll try it. I will let you know what it will end up. Thanks @Alex $\endgroup$
    – user583893
    Feb 2, 2019 at 15:02
  • $\begingroup$ NDSolve still fails NDSolve::mxst: Maximum number of 10000 steps reached at the point y == -0.431183866743419284299032955175. $\endgroup$
    – user583893
    Feb 2, 2019 at 15:07
  • $\begingroup$ There is another problem. You try to calculate large numbers of order Exp[10^3] $\endgroup$ Feb 2, 2019 at 15:59
0
$\begingroup$

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) 
$\endgroup$
6
  • $\begingroup$ Right, so what we do with $y=k$? It blows up $\endgroup$
    – user583893
    Feb 2, 2019 at 15:43
  • $\begingroup$ Yes, and the asymptotic solution shows you how the singularity is "passed". $\endgroup$ Feb 2, 2019 at 15:56
  • $\begingroup$ I don't quite understand what you really mean $\endgroup$
    – user583893
    Feb 2, 2019 at 15:58
  • $\begingroup$ 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 $\endgroup$ Feb 2, 2019 at 16:01
  • $\begingroup$ I get that. But how do I actually evaluate $R$ at $y=k$? I just use $R\sim(k-y)^{l}$? $\endgroup$
    – user583893
    Feb 2, 2019 at 16:05
0
$\begingroup$

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[10] == 0, g'[10] == -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}]

output from Plot

$\endgroup$
17
  • $\begingroup$ Hi @Ruslan, how can I actually retrieve $R(y)$ from the result of NDSolve for $g(t)$? $\endgroup$
    – user583893
    Feb 4, 2019 at 10:36
  • $\begingroup$ @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. $\endgroup$
    – Ruslan
    Feb 4, 2019 at 11:22
  • $\begingroup$ Can you update your answer including the actual code retrieving $R(y)$ @Ruslan? Thank you $\endgroup$
    – user583893
    Feb 4, 2019 at 14:53
  • $\begingroup$ @user583893 it's trivial: R[Y_]=Piecewise[{{Indeterminate, k-ϵ<=Y<=k+ϵ}, {sol[Y], True}}] $\endgroup$
    – Ruslan
    Feb 4, 2019 at 16:51
  • $\begingroup$ Thanks. By the way, what is the motivation for the form of the parametrization? $\endgroup$
    – user583893
    Feb 4, 2019 at 16:55

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.