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$
  • $\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$ – Alex Trounev Feb 2 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 at 13:56
  • $\begingroup$ Watch carefully compare degrees $(k-y)$. $\endgroup$ – Alex Trounev Feb 2 at 14:00
  • $\begingroup$ Oh, sorry, I didn't see the brackets. $\endgroup$ – Alex Trounev Feb 2 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 at 14:12
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$
  • $\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 at 14:50
  • $\begingroup$ Put h=10^-5;, it is enough to get through y=k. $\endgroup$ – Alex Trounev Feb 2 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 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 at 15:07
  • $\begingroup$ There is another problem. You try to calculate large numbers of order Exp[10^3] $\endgroup$ – Alex Trounev Feb 2 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$
  • $\begingroup$ Right, so what we do with $y=k$? It blows up $\endgroup$ – user583893 Feb 2 at 15:43
  • $\begingroup$ Yes, and the asymptotic solution shows you how the singularity is "passed". $\endgroup$ – Ulrich Neumann Feb 2 at 15:56
  • $\begingroup$ I don't quite understand what you really mean $\endgroup$ – user583893 Feb 2 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$ – Ulrich Neumann Feb 2 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 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$
  • $\begingroup$ Hi @Ruslan, how can I actually retrieve $R(y)$ from the result of NDSolve for $g(t)$? $\endgroup$ – user583893 Feb 4 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 at 11:22
  • $\begingroup$ Can you update your answer including the actual code retrieving $R(y)$ @Ruslan? Thank you $\endgroup$ – user583893 Feb 4 at 14:53
  • $\begingroup$ @user583893 it's trivial: R[Y_]=Piecewise[{{Indeterminate, k-ϵ<=Y<=k+ϵ}, {sol[Y], True}}] $\endgroup$ – Ruslan Feb 4 at 16:51
  • $\begingroup$ Thanks. By the way, what is the motivation for the form of the parametrization? $\endgroup$ – user583893 Feb 4 at 16:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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