2
$\begingroup$

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?

Improve this question
$\endgroup$
4
  • 4
    $\begingroup$ 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. $\endgroup$
    – J. M.'s eventual burnout ♦
    Feb 7, 2021 at 21:53
  • 1
    $\begingroup$ You ought to post an example that goes astray, in copy-pastable format. $\endgroup$
    – Daniel Lichtblau
    Feb 7, 2021 at 22:28
  • 1
    $\begingroup$ My understanding in this situation is to indent the contour around the pole over a very small semi-circle. $\endgroup$
    – Dominic
    Feb 8, 2021 at 11:12
  • $\begingroup$ Since Dominic and Carl mention it: use their proposal if you are willing and able to evaluate at complex arguments, and mine otherwise. $\endgroup$
    – J. M.'s eventual burnout ♦
    Feb 9, 2021 at 12:42

2 Answers 2

Reset to default
2
$\begingroup$

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

Improve this answer
$\endgroup$
1
$\begingroup$

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}  *)
Improve this answer
$\endgroup$
3
  • $\begingroup$ I can't help but feel that with a route like this, one should just be doing Clenshaw-Curtis directly... (+1 of course) $\endgroup$
    – J. M.'s eventual burnout ♦
    Feb 9, 2021 at 18:14
  • $\begingroup$ @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. $\endgroup$
    – Michael E2
    Feb 9, 2021 at 18:31
  • 1
    $\begingroup$ In case someone wants to bypass NIntegrate: {abscissae, weights} = Most@NIntegrate`GaussRuleData[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 NIntegrate`GaussKronrodRuleData to help with estimating the error. $\endgroup$
    – Michael E2
    Feb 9, 2021 at 18:38

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.