3
$\begingroup$

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.

$\endgroup$
6
  • 1
    $\begingroup$ Try Chop[N[N[T,100]]]. $\endgroup$
    – user293787
    Oct 3, 2022 at 8:43
  • $\begingroup$ @user293787 That transforms the symbolic integration values into those obtained by numerical integration. How can Chop transform something of order $10^6$ into $4233$? $\endgroup$ Oct 3, 2022 at 8:51
  • $\begingroup$ @user293787 I do not see why Chop has any effect. Screen shot !Mathematica graphics it is clear the results are different. $\endgroup$
    – Nasser
    Oct 3, 2022 at 8:55
  • $\begingroup$ 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. $\endgroup$ Oct 3, 2022 at 8:59
  • 1
    $\begingroup$ 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. $\endgroup$
    – user293787
    Oct 3, 2022 at 9:01

1 Answer 1

7
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.