Why NDSolve has problem with Gamma[5/2,0,x]?

I was trying to solve three arrays of Poisson's equation with three different incomplete gamma functions and 10 different initial conditions for each. I did the following (taking my first array as an example):

ϵ = 10^-6
G1[x_] := Exp[x] Gamma[3/2, 0, x]
T1 = {ψ[i]''[r] + 2/r ψ[i]'[r] == -9*G1[ψ[i][r]]/G1[i]}
T2 = {ψ[i][ϵ] == i}
T3 = {ψ[i]'[ϵ] == 0}
eqn1 = Join[T1, T2, T3]
psi1 = Table[NDSolve[eqn1, ψ[i], {r, ϵ, 20}], {i, 1, 10}]


Then I just repeated the above for Gamma[5/2, 0, x] and Gamma[7/2, 0, x]. The first and third array worked fine, but I couldn't get the equations solved for Gamma[5/2, 0, x]. There was no warnings, the programme just kept running. It seems like I have to narrow down my domain to around {r, ϵ, 2}. Is there something special about this gamma function that made the computation complicated, or is there any flaw with my codes that created a pitfall at this point?

• It should be :G1[x_] := Exp[x] Gamma[3/2, 0, x] – Mariusz Iwaniuk Feb 27 '17 at 21:22
• There is a discontinuity in the imaginary part of the principal value of Gamma[5/2, 0, z] along the negative real axis. See dlmf.nist.gov/8.3.F15.mag or Plot3D[Im@Gamma[5/2, 0, x + I y], {x, -1, 1}, {y, -1, 1}]. When the solution becomes negative, the approximate floating point calculations might jump back and forth because of the sign change in the discontinuity, which for 5/2 is the opposite of 3/2, 7/2. – Michael E2 Mar 2 '17 at 12:13

For $a = (2n+1)/2$, the incomplete gamma function $\gamma(a,z)$ has two branches,

Gamma[a, 0, x + I y]
Exp[-2 Pi I * a] Gamma[a, 0, x + I y] (* yes, this is just  -Gamma[a, 0 x + I y],
but I kept the general form to play with  a  *)


which are branched along the negative real axis (see DLMF (8.2.8)):

Real and imaginary parts of $\gamma(a,z)$ for $a=\textstyle{3\over2},{5\over2},{7\over2}$.

GraphicsGrid[
Table[
Plot3D[
{part[Gamma[a, 0, x + I y]],
part[Exp[-2 Pi I *a] Gamma[a, 0, x + I y]]},
{x, -1, 1}, {y, -1, 1},
PlotLabel ->
With[{p = part, a = a},
HoldForm[p@γ[Style[a, Small], x + I y]]],
ViewPoint -> {-2.5, 2, 1.2}, AxesLabel -> {x, y, part@γ}],
{part, {Re, Im}}, {a, {3/2, 5/2, 7/2}}],
ImageSize -> 800
]


The biggest problem is that each IVP drives the solution ψ[r] down the positive real axis through though the branch point in γ[a, ψ[r]] at ψ[r] == 0. The question is on what branch to exit the branch point? Either seems equally valid. A way to choose the branch is presented further down, after we introduce the notion of a winding number.

Assuming one wants to avoid the discontinuity in the principal value (yellow branch) when crossing the negative real axis, one can maintain a winding number w using WhenEvent and use the general value of gamma:

γ[a_, z_?NumericQ, w_] := Exp[-2 Pi*I*w*a] Gamma[a, 0, z];

WhenEvent[Im[ψ[i][r]] > 0 && Re[ψ[i][r]] < 0, (* cross - real axis counterclockwise *)
w[r] -> w[r] + 1,
"LocationMethod" -> "StepBegin"],
WhenEvent[Im[ψ[i][r]] < 0 && Re[ψ[i][r]] < 0, (* cross - real axis clockwise *)
w[r] -> w[r] - 1,
"LocationMethod" -> "StepBegin"]


This works, but crossing the branch introduces error, because the winding number is changed only at a step taken by NDSolve. In computing a step when the branch is crossed, it is very likely that $\gamma(a,z)$ is evaluated at points on either side of the branch.

One can move the branch cut to the positive real axis:

γ[a_, z_?NumericQ, w_] := Exp[-2 Pi*I*(w - UnitStep@Sign@Im[z] + 1)*a] Gamma[a, 0, z]


This seems to produce a better result, perhaps because the difference between the branches seems less in magnitude there. because it does not cross the positive real axis after the solution becomes complex.

To address the "big problem" of choosing a branch when the solution passes through the origin ne could select the branch by setting the winding number when the origin is crossed. For example, this choose the non-principal value of gamma:

WhenEvent[Re[ψ[i][r]] < 0 && Abs@Im[ψ[i][r]] < 1*^-8, w[r] -> 1]


Finally, if you want to move the initial condition to r == 0, that is possible, too.
\eqalign{ v_i'(r) &= -9\,r^2\, G1(\psi_i(r))/G1(i)\strut \cr \psi_i'(r) &= \cases{v_i(r)/r^2 & r \ne 0 \cr 0 & otherwise} \cr }

Code:

T1[i_] := {
v'[r] == -9*r^2 G1[ψ[i][r], w[r]]/G1[i, 0],
ψ[i]'[r] == Piecewise[{{v[r]/r^2, r != 0}}]};


Here's the complete approach:

Clear[T1, T2, T3, evts, eqn1, γ, psi2];
ϵ = 0;
a = 5/2;
γ[a_, z_?NumericQ, w_] :=          (* move branch cut to +real axis *)
Exp[-2 Pi*I*(w - UnitStep@Sign@Im[z] + 1)*a] Gamma[a, 0, z];
G1[x_, w_] := Exp[x] γ[a, x, w];
T1[i_] := {v'[r] == -9*r^2 G1[ψ[i][r], w[r]]/G1[i, 0],
ψ[i]'[r] == Piecewise[{{v[r]/r^2, r != 0}}]};
T2[i_] := {ψ[i][ϵ] == i};
T3[i_] := {v[ϵ] == 0, w[ϵ] == 0};
evts[i_] :=                        (* switch branch index on crossing +real axis *)
{WhenEvent[Im[ψ[i][r]] < 0 && Re[ψ[i][r]] > 0,
w[r] -> w[r] + 1,
"LocationMethod" -> "StepBegin"],
WhenEvent[Im[ψ[i][r]] > 0 && Re[ψ[i][r]] > 0,
w[r] -> w[r] - 1,
"LocationMethod" -> "StepBegin"]};
branchpt[i_, w0_] :=               (* choose branch on passing through origin *)
WhenEvent[Re[ψ[i][r]] < 0 && Abs@Im[ψ[i][r]] < 1*^-8, w[r] -> w0];
eqn1[i_] := Join @@ Through[{T1, T2, T3, evts}[i]];

mem : psi2[bp_] := mem = Flatten@Table[
NDSolve[{eqn1[i], branchpt[i, bp]}, ψ[i], {r, ϵ, 20},
DiscreteVariables -> {w}, Method -> "Extrapolation"],
{i, 1, 10}];
psi2 /@ {0, 1};


Imaginary parts of the two sets of solutions track opposite branches:

Real and imaginary parts of ψ[i][r] for each choice of branch (psi2[0], psi2[1]); the real parts (blue, green) coincide, but the imaginary parts (gold, orange) have opposite signs.

Table[
Plot[Flatten@{ReIm[ψ[i][r] /. psi2[0]], ReIm[ψ[i][r] /. psi2[1]]} // Evaluate,
{r, ϵ, 20}, PlotLabel -> ψ[i], PlotRange -> {-0.5, 0.5}, Frame -> True],
{i, 10}] // Multicolumn[#, 3] &


Can't reproduce your problem in v9.0.1, the calculation finishes in about 30 seconds, but NDSolve uses up its 10000 steps before r = 20 for i < 6. Further check shows this problem is actually similar to this one:

Orbit followed by a particle around Schwarzschild Black Hole

i.e. the culprit is the very small imaginary part of G1[x] caused by numerical error, so it can be fixed by adding a Re to chop the imaginary part:

T1 = {ψ[i]''[r] + 2/r ψ[i]'[r] == -9 Re@G1[ψ[i][r]]/G1[i]};
T2 = {ψ[i][ϵ] == i};
T3 = {ψ[i]'[ϵ] == 0};
eqn1 = {T1, T2, T3};

psi1 =
Table[NDSolve[eqn1, ψ[i], {r, ϵ, 20}], {i, 1, 10}]; // AbsoluteTiming
(* {0.152637, Null} *)

Plot[#[[1, 1, -1]][x] & /@ psi1 // Evaluate, {x, ϵ, 20},
PlotRange -> All]