# Symbolic integration and numerical integration yield very different results

This is a topic that has been discussed before, but after reading other posts I could not find a solution that applies for my specific case.

I have the following integral, which I need to evalutate at $$n=1,2,3,4...$$ and $$a=0.1$$:

$$I(n)=\int_{-\pi}^{\pi} {\rm{d}}x \frac{\Gamma(n,0,2ix)}{(1+iax)^n}.$$

The problem is that Integrate and NIntegrate seem to yield different results for sufficiently large $$n$$. Numerical integration with $$a=0.1$$ leads to

a = 1/10;
Table[NIntegrate[Gamma[n, 0, 2 I x]/(1 + I a x)^n, {x, -\[Pi], \[Pi]}], {n, 1,6}] // Chop
{5.82505, 9.23019, 42.4493, 117.67, -303.323, -4233.64}


However, integrating first and then substituting $$a=0.1$$ returns the following values:

T = Table[Integrate[Gamma[n, 0, 2 I x]/(1 + I a x])^n, {x, -\[Pi], \[Pi]}, Assumptions->a>0],
{n, 1, 6}];
a=1/10;
T//N//Chop
{5.82505, 9.22998, 41.7814, 0.951343, -26457.8, -3.06652*10^6}


which only agree with numerical integration for $$n=1$$ and then they become very different. Interestingly, using Integrate with $$a=0.1$$ yields another different set of values:

a = 1/10;
TT = Table[Integrate[Gamma[n, 0, 2 I x]/(1 + I a x)^n, {x, -\[Pi], \[Pi]}, Assumptions -> a > 0], {n, 1, 6}];
TT // N // Chop
{5.82505, 9.23058, 42.4145, 164.712 + 1.19844*10^-9 I, 0, 701973.}


All these evaluations occur without mathematica returning any errors. What is happening here? Which values are the correct ones?

Things that I know/ tried:

1. The integral is real, as the imaginary part of the integrand is an odd function of x.
2. Assumptions->a>0 is necessary, since for arbitrary $$a$$ the default solution mathematica returns is valid in the complex plane of $$a$$ minus the positive real axis.
3. From external arguments, the values that I am most inclined to trust are those from exact integration for arbitrary a>0 and then numerical evaluation at a=0.1 (the second set of numbers).
4. The integrand oscillates more and more the higher $$n$$ is, and I believe NIntegrate may have problems getting the right answer for highly oscillatory integrands. However, since it does not return any error I do not know if this is the issue.
5. The problem persists if I use GammaRegularized(n,0,2ix) to work with smaller numbers.
6. The problem persists if I increase numerical precision (I am not very familiar with how numerical precision should be managed, maybe I am doing it the wrong way: I increase $MaxExtraPrecision and WorkingPrecision in NIntegrate) I've read similar posts in which integration methods for NIntegrate such as ExtrapolatingOscillatory and DoubleExponentialOscillatory where suggested, but as far as I know those only work for integrals with infinite range. • Try Chop[N[N[T,100]]]. Oct 3, 2022 at 8:43 • @user293787 That transforms the symbolic integration values into those obtained by numerical integration. How can Chop transform something of order$10^6$into$4233$? Oct 3, 2022 at 8:51 • @user293787 I do not see why Chop has any effect. Screen shot !Mathematica graphics it is clear the results are different. Oct 3, 2022 at 8:55 • I understand now, this is what happens: the symbolic integration expression for lets say$n=6$is very, very long. When I evaluated it at$a=0.1\$ and used N, only 6 digits of precision were kept in each of the many terms and the total evaluation was wrong. With @user293787 command, 100 digits are kept in the evaluation, which drastically changes the result. Numerical integration was right all along. Oct 3, 2022 at 8:59
• The Chop[N[...]] is not important, that was just to produce a compact output. The main thing is N[T,100]. It seems that there are considerable cancellations among terms that lead to large errors with machine precision N. You could look at N[N[List@@@Expand[T],100]] to see the cancellations. Oct 3, 2022 at 9:01

The problem is not Integrate, but the numerical evaluation following Integrate.

Consider for example T[[4]]:

T[[4]]
(*
(1/(3 (100 + Pi^2)^3))80000 (48003000 Pi -
10520000000 E^20 Pi + 998810 Pi^3 - 315600000 E^20 Pi^3 +
5257 Pi^5 - 3156000 E^20 Pi^5 - 10520 E^20 Pi^7 +
5260000000 I E^20 ExpIntegralEi[-20 - 2 I Pi] +
157800000 I E^20 Pi^2 ExpIntegralEi[-20 - 2 I Pi] +
1578000 I E^20 Pi^4 ExpIntegralEi[-20 - 2 I Pi] +
5260 I E^20 Pi^6 ExpIntegralEi[-20 - 2 I Pi] -
5260 I E^20 (100 + Pi^2)^3 ExpIntegralEi[-20 + 2 I Pi])
*)


When one evaluates this numerically, the result depends strongly on how the expression is organized and on what precision is used:

N[T[[4]]]
(* 164.712 +1.19844*10^-9 I *)

N[Expand[T[[4]]]]
(* 256. +0. I *)

N[T[[4]],20]
(* 117.670138197013970506+0.*10^-19 I *)


The reason is that individual summands are many orders of magnitude larger than the final result:

Chop[N[N[List@@Expand[T[[4]]],100]]]
(*
3.03218*10^6
-3.22397*10^17
622685.
-9.5458*10^16
32346.2
-9.42133*10^15
-3.09949*10^14
1.61199*10^17-4.62852*10^6 I
4.7729*10^16-1.37045*10^6 I
4.71066*10^15-135258. I
1.54975*10^14-4449.81 I
1.61199*10^17+4.62852*10^6 I
4.7729*10^16+1.37045*10^6 I
4.71066*10^15+135258. I
1.54975*10^14+4449.81 I
*)


Therefore, higher than machine precision is needed to obtain the correct value, which in this case is the one near 117.67.