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

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

• I don't think this is a power series because your $c_n$ depend on $a$ as well. Is this intended? Mar 27 at 11:50
• If your input contains machine numbers such as 0.2 then NIntegrate will not be able to use higher precision. So first thing is to give either exact input or else input at a precision at least as high as the WorkingPrecision setting. Mar 27 at 14:49
• @Roman Yes, the $a$ dependence of the $c_n$ coefficients is intended. I am comparing the behaviour of the fully expanded (in $a$) power series with this 'rearrangment'. If you are familiar with particle physics, results are usually given as power series in some variable $a(\mu)$ which satisfies a first-order differential equation with an initial condition $a(\mu_0)\equiv a$; schematically, it is of the form $a(\mu)=a/(1+a \log(\mu/\mu_0))$. That leads to power series such as mine. Mar 28 at 8:16
• @DanielLichtblau understood, thanks! Mar 28 at 9:22

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."

• Thank you so much, the approach with Gamma[n,0,2 I f] vastly improves the computing time of my code. Anallytically it is possible to proof that the integrand is an odd function of $\varphi$, so indeed the imaginary part must be discarded. Mar 28 at 18:34
• @NéstorGonzálezGracia You're welcome. I'm glad it helps. Mar 28 at 18:35
• Cool stuff Michael. Just a quick question: since we are dealing with special functions, shouldn't methods of NIntegrate help with the performance? Like GlobalAdaptive or the local version. I did not check for myself, so this is based on past experience. Maybe too naive for this example to be honest.
– bmf
Mar 28 at 18:45
• @bmf As a gets bigger, the singularity at I/a gets closer to the path of integration, so that might interfere with convergence and "LocalAdaptive" might help. For some reason, I've never figured out how to predict when local adaptive will save time and when it will waste it. It's not simply "special functions," I don't think. (The default is "GlobalAdaptive" with "GaussKronrodRule" and 5 Gauss points.) Mar 28 at 19:08
• Thanks for your comment. I, too, never figured out how to predict whether it's useful or not. Usually I apply the good old "trial and error" method, so I thought to ask. Thanks for your thoughts :-)
– bmf
Mar 28 at 19:10

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:

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.

0.2 // Rationalize[#, 0] &


gives back

1/5

1. 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.

• @NéstorGonzálezGracia Just as a reference, I executed again the table with Rubi and I am providing a screenshot. Is it possible you have an old Mathematica version? Did you get the latest Rubi version?
– bmf
Mar 27 at 10:16
• Thanks for your answer! I know mathematica can solve the integral analytically for $n=1,2,3...$, but it simply takes too long for large $n$. Regarding Rubi, I didn't know about it, so I installed it and gave it a try. Copying what you had in your screenshot I got 66 seconds of computation time (!) and also the results are given explicitely in terms of the limits of the antiderivative when $f\mapsto\pm\pi$ (I don't know how to add a screenshot to a comment). I am using Mathematica 12.0 and the latest version of Rubi. Mar 27 at 10:25
• @NéstorGonzálezGracia no need to add a screenshot, I understand what you are seeing. That's odd. We have the same Mathematica and Rubi versions. It's confusing. Could you please update your post with the code you tried? I have to sleep now, but I will try to have a look tomorrow, and maybe others can trouble-shoot what's going on before me :-)
– bmf
Mar 27 at 10:27
• For some reason adding Quit at the beginning of the code solved both problems. I am going to update my post now. Mar 27 at 10:53
• @NéstorGonzálezGracia oh yes. sure. You probably had some interfering definitions before. Glad it's working. Take into consideration that the Rubi results are without restrictions on a` hence some inspection and caution are required :-)
– bmf
Mar 27 at 11:02