Best way to compute the following numerical integral?

I want to perform

NIntegrate[-(2 u)^1000 E^(-2000 u) 1/(1 - 2 u), {u, 0, Infinity}]

When I attempt this on my machine, Mathematica warns me that

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections

I am a very pedestrian user of numerical integration, but I am aware that Mathematica has different strategies for this purpose. Can anybody recommend me a proper approach to tackle this numerical integral?

• NIntegrate cannot integrate stuff that is not purely numerical: you seem to want this integral for some unspecified N. Another problem is that N is a reserved Mathematica function name. Trying Integrate[-(2 u)^n E^(-2 n u) 1/(1 - 2 u), {u, 0, Infinity}] instead, Mathematica tells me after a few seconds that the integral does not converge. – Marius Ladegård Meyer Oct 31 '16 at 10:59
• @MariusLadegårdMeyer corrected – PhoenixPerson Oct 31 '16 at 11:23
• Since you have a singularity at u == 1/2, you need to specify that you want a principal value. The symbolic answer is more useful because you have such crazy huge exponents: Integrate[-(2 u)^n E^(-2 n u) 1/(1 - 2 u), {u, 0, Infinity}, PrincipalValue -> True] gives something... – Marius Ladegård Meyer Oct 31 '16 at 12:29
• @MariusLadegårdMeyer I know that I can perform it analytically. I have that. I want a numerical check. – PhoenixPerson Oct 31 '16 at 12:33

This is how the integrand looks like:

GraphicsGrid[
Partition[
Table[Plot[-((2*u)^n/(E^(2*n*u)*(1 - 2*u))),
{u, 0.4, 0.6}, Frame -> True, PlotRange -> {All, {-1, 5}}, MaxRecursion -> 12,
PlotLabel -> StringJoin["n = ", ToString[n]]],
{n, 2, 12, 2}], 3], ImageSize -> 600, Spacings -> 0] It has an asymptote at $u=0.5$. Let's split the integral into two parts excluding the singularity (for illustrative purposes - see also the comment at the end of this post):

ε = 0.00001;

int[n_] :=
NIntegrate[-(2 u)^n E^(-2 n u) 1/(1 - 2 u), {u, 0, 0.5 - ε}] +
NIntegrate[-(2 u)^n E^(-2 n u) 1/(1 - 2 u), {u, 0.5 + ε, Infinity}]

and

ListPlot[Table[{n, int[n]}, {n, 1, 20}], Frame -> True,
PlotRange -> All, PlotStyle -> {Red, PointSize[Large]},
FrameLabel -> {"n", "int[n]"}] int is $1.92\cdot 10^{-10}$, and int - $1.55\cdot 10^{-45}$. The value n = 1000 gives a numerical zero.

Those values are the same if obtained with Integrate performed with PrincipalValue -> True - see Marius's comment; in general

Normal @ Integrate[-(2 u)^n E^(-2 n u) 1/(1 - 2 u), {u, 0, Infinity},
PrincipalValue -> True] with a condition Re[n] >= 0.

Explicitly listing the singularity at 1/2 gives an error-free result:

NIntegrate[-(2 u)^1000 E^(-2000 u) 1/(1 - 2 u), {u, 0, 1/2, Infinity},
MaxRecursion -> 20, WorkingPrecision -> 16,
Method -> "PrincipalValue"]
(*  6.705858148612409*10^-437  *)

It is stable under increasing WorkingPrecision, PrecisionGoal and MinRecursion. Update: It agrees with the following symbolic result:

int = Integrate[-(2 u)^1000 E^(-2000 u) 1/(1 - 2 u), u];
Block[{\$MaxExtraPrecision = 500},
N[Limit[int, u -> Infinity] - Limit[int, u -> 0], 16]
]
(*  6.705858148607688*10^-437  *)

The agreement between the numeric result above and the exact result improves as the WorkingPrecision increases. So this now seems to be the correct result. (Compare with the result below.)

Original: The symbolic result of Marius Ladegård Meyer is somewhat different from the answer above:

val = Integrate[-(2 u)^n E^(-2 n u) 1/(1 - 2 u), {u, 0, Infinity},
PrincipalValue -> True];
N[Limit[val, n -> 1000], 8]
(* 1.26877823*10^-438 + 0.*10^-447 I  *)

Apparently this symbolic result is erroneous.