# Most efficient strategy for integrating over removable poles?

I am finding many situations where I have to numerically integrate some function $$f(x)$$ of the form:

$$f(x)=f_{s}(x)-ax^{-n},$$

where $$f_s$$ is a special function with a finite-order pole that is canceled by the $$-ax^{-n}$$ term. I.e., $$f_s(x)$$ has Laurent expansion

$$f_s(x)=ax^{-n}+a_0+a_1x+a_2x^2+\cdots.$$

As an example, consider $$f_s(x)=E_{-n}(x),$$ or fs[x_]:=ExpIntegralE[-n,x], and $$f(x)=E_{-n}(x)-(n!)x^{-n},$$ which has a removable singularity at $$x=0.$$

In most cases I have found that integrating directly over $$f(x)$$ works, but sometimes I have to crank up WorkingPrecision fairly high in order for the cancellation near the pole to work out. And I think this method of computation must be inefficient, because Mathematica must compute many digits of precision that only end up canceling.

So is there a better/more computationally efficient strategy for solving these integrals in general? In particular, I wonder if there is a more efficient way of computing $$f(x)$$ near the removable poles in general, even before integration comes into play. Since Mathematica has built-in algorithms to calculate special functions to arbitrary precision (presumably as efficiently as possible), I wonder if it has smart ways of computing only the "non-pole" part of special functions near their poles?

• Define a new function to integrate, such that you evaluate an e.g. Padé approximant in the neighborhood of the singularity, and the function itself otherwise. Feb 7 at 21:53
• You ought to post an example that goes astray, in copy-pastable format. Feb 7 at 22:28
• My understanding in this situation is to indent the contour around the pole over a very small semi-circle. Feb 8 at 11:12
• Since Dominic and Carl mention it: use their proposal if you are willing and able to evaluate at complex arguments, and mine otherwise. Feb 9 at 12:42

One idea is to deform the integration contour around the singularity. For your example:

f[x_] := ExpIntegralE[-5, x] - 5! x^-6
NIntegrate[f[x], {x, -1, I, 1}]


-0.3767 - 2.02699*10^-13 I

We can check by integrating the series approximation:

g[x_] = Normal @ Series[f[x], {x, 0, 12}]
NIntegrate[g[x], {x, -1, 1}]


-(1/6) + x/7 - x^2/16 + x^3/54 - x^4/240 + x^5/1320 - x^6/8640 + x^7/65520 -
x^8/564480 + x^9/5443200 - x^10/58060800 + x^11/678585600 - x^12/8622028800

-0.3767

Compare this to integrating along the real axis:

NIntegrate[f[x], {x, -1, 1}]


NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {-0.0000310833}. NIntegrate obtained 2.926094453677552*^9 and 2.926342889896642*^9 for the integral and error estimates.

2.92609*10^9

I think a Chebyshev method could be adapted to your workflow. I don't know what your workflow is, so I don't have any advice about that. Here's a comparison with NIntegrate and @Carl Woll's example.

(ClearAll[f];
f0[x_] := ExpIntegralE[-5, x] - 5! x^-6;
f[0] = SeriesCoefficient[f0[x], {x, 0, 0}];
f[0.] = N@f[0];
f[x_] = f0[x];
pp = 16 (* order*);
wp = MachinePrecision;
xx = N[Sin[Pi/2 Range[pp, -pp, -2]/pp], wp];
yy = f /@ xx /. y_ /; Im[y] == 0 :> Re[y];
cc = Sqrt[2/pp] FourierDCT[yy, 1];
cc[[{1, -1}]] /= 2;
cc[[;; ;; 2]] .
Table[(1 + Cos[n \[Pi]])/(
1 - n^2), {n, 0, pp, 2}]) // RepeatedTiming
(*  {0.000521951, -0.3767 + 3.44643*10^-16 I}  *)

NIntegrate[f0[x], {x, -1, I, 1}] // RepeatedTiming
(*  {0.00537676, -0.3767 - 2.02699*10^-13 I}  *)


High-precision NIntegrate does not complain and is just a bit slower than machine-precision NIntegrate:

NIntegrate[f0[x], {x, -1, 1}, PrecisionGoal -> 6,
WorkingPrecision -> 32] // RepeatedTiming
(*  {0.00788756, -0.37670013227515802195163230745995}  *)


The Chebyshev code with working precision wp = 32:

(*  {0.00169116, -0.37670013227515802195163}  *)

• I can't help but feel that with a route like this, one should just be doing Clenshaw-Curtis directly... (+1 of course) Feb 9 at 18:14
• @J.M. You're right, of course. In fact, one can bypass NIntegrate and manually do Gauss or Clenshaw-Curtis directly over the complex path, if speed is important. I was hoping the OP would update with more detail, because I thought having the Chebyshev approximation might turn out to be useful. Since the integrand is in fact analytic, either method over either path should work well if the function can be computed accurately. Probably NIntegrate over a complex path is the most convenient way, but I got tired of waiting and posted what I had before my computer crashes again. Feb 9 at 18:31
• In case someone wants to bypass NIntegrate: {abscissae, weights} = Most@NIntegrateGaussRuleData[5, MachinePrecision]; ((I - (-1)) (f0 /@ Rescale[abscissae, {0, 1}, {-1, I}]) . weights + (1 - I) (f0 /@ Rescale[abscissae, {0, 1}, {I, 1}]) . weights) // RepeatedTiming --> {0.000248553, -0.3767 + 1.32808*10^-12 I}. A drawback is that the programmer is responsible for the managing the error of approximation. One could use NIntegrateGaussKronrodRuleData to help with estimating the error. Feb 9 at 18:38