Best way to handle numerical integration and power series with large numbers

Question

Due to the fact that almost everything I do in my research is analytic, I am quite unfamiliar with numeric calculus, so I was wondering if anyone could give me some advice on the most efficient way to face the following problem.

I have a power series of the form $S\equiv\sum_{n=1}^\infty c_n a^n$, where $a\in(0,1)$ and the coefficients are:

$$c_n\equiv \frac{1}{2^{n+1}\pi} \int_{-\pi}^\pi d\varphi\,{}\frac{\Gamma(n)-\Gamma(n,2i\varphi)}{(1+ia\varphi)^n},$$ where of course $i$ is the imaginary unit. The problem here is that the integral defining $c_n$ , despite convergent, does not seem to have a solution in terms of the known functions. This means that any analysis on the series $S$ must be performed numerically.

My main concern is to get an idea of the convergence radius of the series, which people claim (with no proof) that is somewhere between $1/\pi$ and $1$. To do this I need to extract the numerical values of $c_n$ for large $n$, as big as it is possible, and perform things such as the root and ratio tests. The only trick I could think of to improve the numeric computations is to write $$c_n= \frac{\Gamma(n)}{2^{n+1}\pi} \int_{-\pi}^\pi d\varphi\,{}\frac{1-\tilde{\Gamma}(n,2i\varphi)}{(1+ia\varphi)^n},$$ where $\tilde{\Gamma}$ is the regularized version of the (incomplete) gamma function: $\tilde{\Gamma}(a,z)=\Gamma(a,z)/\Gamma(a)$. This, I believe, avoids cancelations between large numbers.

The main issues I find are (I am using the NIntegrate mathematica function):

• I can only get to $n\approx 200$. For higher values I get the folllowing error:

NIntegrate::errprec: Catastrophic loss of precision in the global error estimate due to insufficient WorkingPrecision or divergent integral.

• Increasing WorkingPrecision to 20 returns two more errors:

NIntegrate::precw: The precision of the argument function ((1-GammaRegularized[241,2 I \[Phi]])/(1+(0. +0.2 I) \[Phi])^241) is less than WorkingPrecision (20.`).

NIntegrate::izero: Integral and error estimates are 0 on all integration subregions. Try increasing the value of the MinRecursion option. If value of integral may be 0, specify a finite value for the AccuracyGoal option.

• Finally, I have observed that several values of $n$ (such as $n=32$ to $35$ for $a=0.2$) return the WorkingPrecision error despite the higher values don't.

In sum, ony being able to extract the first 200 $c_n$ coefficients is not enough to get a good approximation for the convergence radius, so I need to find a way of efficiently computing more of them.

Edit1: Setting WorkingPrecision->10 and AccuracyGoal->10 and in NIntegrate seems to get rid of all the errors and lets me go to higher $n$. This is because the $c_n$ get very small very quickly, and, if I understand correctly, AccuracyGoal->10 sets them to 0 as soon as the first 10 digits are 0.

Edit2: Following the suggestion by @bmf of using Rubi, the integrals can be now solved anallytically under acceptable computation time. I have yet to test how far in $n$ is it possible to go, but this path seems promising. The questions that I now have about Rubiare

• Is it possible to add assumptions to the Intfunction?. I am only interested on $a>0$.
• Why the code provided in the comments, i.e.,

    Get["Rubi`"]
    d[n_] := 1/(2^(
         n + 1) \[Pi]) Int[(Gamma[n] - Gamma[n, 2 I f])/(1 + I a f)^
         n, {f, -Pi, Pi}]
    rubitable = Table[d[n], {n, 1, 5}]; // AbsoluteTiming
    rubitable // TableForm

works fine for generic $a$ but gets many, many errors for numeric $a$?

    Get["Rubi`"]
    a=0.2;
    d[n_] := 1/(2^(
         n + 1) \[Pi]) Int[(Gamma[n] - Gamma[n, 2 I f])/(1 + I a f)^
         n, {f, -Pi, Pi}]
    rubitable = Table[d[n], {n, 1, 5}]; // AbsoluteTiming
    rubitable // TableForm

Answer 1

Some of the numerical problem probably comes from subtracting the two gamma terms. One can use the form Gamma[n, 0, 2 I f] to avoid that. There is also GammaRegularized[n, 0, 2 I f] which gave a slight advantage in timing that I can't explain. Finally, it appears the imaginary parts of the integral cancel, which might be another source of numerical trouble. Certainly one can speed things up by taking the real part of the integrand.

(* using Gamma[n, 0, 2 I f] and high precision to tame the numerics --
   not too slow, in fact *)
cc[n_, a_] := 
 1/(2^(n + 1) Pi) NIntegrate[
   Gamma[n, 0, 2 I f]/(1 + I a f)^n, {f, -Pi, Pi}, 
   WorkingPrecision -> 50, PrecisionGoal -> 8]

(* using GammaRegularized[] and the real part to tame the numerics -- 
   a bit faster *)
ccc[n_, a_] := 
 1/(2^(n + 1) Pi) Gamma[n] NIntegrate[
   GammaRegularized[n, 0, 2 I f]/(1 + I a f)^n // Re, {f, -Pi, Pi}]

cc[200, 1/5] // AbsoluteTiming
Style[ccc[200, 1/5], PrintPrecision -> 17] // AbsoluteTiming

(*
{0.039149, 1.4619894219240324987154273848708784861578925349920*10^79}
{0.016952, 1.461989421924107`*^79)}
*)

If you remove Re, you get an imaginary component that is relatively small (rel. err. ~ 10^-14). In this case (n = 200, a = 1/5) that means it's around 10^65 I, which is sometimes disturbing to those who are not used to thinking of such numbers as "small."

Answer 2

This is just a comment, but it's too extended to be left under the OP. Hopefully interesting and helpful.

My main concern is that it does not seem that you need to do numerics. At least not immediately. Is it possible that you didn't code something properly? Or maybe you have an old version? Please see below.

I am currently using version 12.

With regards to the first formulation of the integral, i.e

$$c_n\equiv \frac{1}{2^{n+1}\pi} \int_{-\pi}^\pi d\varphi\,{}\frac{\Gamma(n)-\Gamma(n,2i\varphi)}{(1+ia\varphi)^n},$$

I defined

c[n_] := 1/(2^(n + 1) Pi)
   Assuming[0 < a < 1, 
   Integrate[(Gamma[n] - Gamma[n, 2 I f])/(1 + I a f)^
    n, {f, -Pi, Pi}]]

and then called

c[1]

to get

So, Mathematica seems to be willing to do it analytically.

For completeness and clarity, I tried a couple of values

mmatable = Table[c[n], {n, 1, 5}];

to see if it works. The output is

mmatable // TableForm

Then just for fun, I wanted to see if Mathematica would do the integral without any assumption on a. It was taking too long for my taste to be honest,so I turned to Rubi

Get["Rubi`"]

d[n_] := 1/(2^(n + 1) Pi)
   Int[(Gamma[n] - Gamma[n, 2 I f])/(1 + I a f)^n, {f, -Pi, Pi}]

rubitable = Table[d[n], {n, 1, 5}];

rubitable // TableForm

---

Edit: Large values of $n$

For $n=54$, performed the integration in a couple of seconds. Keep in mind that we don't use any restrictions on a here.

d[54]

Mathematica needed some more time. Significantly more and it's very late here. So, I tried $n=11$ just to see if she does it.

c[11]

---

Edit:

• 1. Addressing the numerical issue:

When Rubi complained to you the error was self-explanatory. What you need to avoid such situations, is the command Rationalize. From the documentation, if you check the Details you can see that

Rationalize[x,0] converts any inexact number x to rational form.

For instance, in your example

0.2 // Rationalize[#, 0] &

gives back

1/5

2. Tension between Mathematica and Rubi

I mentioned that some trouble-shooting and checks were necessary, so I am providing two.

Before I proceed have a look at what Chop does.

Results from Mathematica:

Table[c[nn] /. a -> 0.2, {nn, 1, 7}] // Chop

Table[c[nn] /. a -> 0.8, {nn, 1, 7}] // Chop

Results from Rubi:

Table[d[nn] /. a -> 0.2, {nn, 1, 7}] // Chop

Table[d[nn] /. a -> 0.8, {nn, 1, 7}] // Chop

Punchline: There seems to be some tension for the $n=1$ term between the two. I am not sure which one is correct, though intuitively -following the pattern of the numbers- I am inclined to say that Rubi got wrong the $n=1$ for some reason. Even if this is the case, for $n>1$ they both give the same results -as it seems so far- and Rubi is just much more efficient.