Coarse-graining in numerical integrations

Question

I have been working recently in a coarse-graining problem I found when using NIntegrate:

I am trying to evaluate the function

$$f(a)=\int_0^{\infty}x\,e^{-(a^2+b^2)x^2}\text{d}x$$

numerically as a function of $a$, where $b=0.001$. To do so, I am writing the following code:

ListPlot[
   ParallelTable[
     {a,NIntegrate[x Exp[-(a^2 + 0.001^2) x^2], {x, 0, 8000}]},{a, 0.001, 1, 1/100}],
        Joined->True]]

Where the upper limit is just a cuttoff. Here is the result of integration (in red), compared with the analytical form of the integral, $f(a)=\frac{1}{2(a^2+b^2)}$ (in blue).

As you can see, there is a coarse-graining problem starting in $a\sim0.86$. I also tried to split the integral in parts with the following function

SplitInt[f_, x_, xmin_, xmax_, div_] := 
 Sum[NIntegrate[f, {x, xmin + (xmax - xmin)/div i, xmin + (xmax - xmin)/div (i + 1)},
AccuracyGoal -> 9, MaxRecursion -> 22, Method -> "LocalAdaptive"], {i, 0, div - 1}]

And now

ListPlot[ParallelTable[{a, 
  SplitInt[x Exp[-(a^2 + 0.001^2) x^2], x, 0, 8000, 
   16]}, {a, 0.001, 0.9, 1/100}]]

moves the problem to another region, as seen below:

As before, the blue curve is the analytic integration and the red one is the numerical one.

I have tried also different integration strategies, namely GlobalAdaptive, LocalAdaptive and GaussKronrodRule, but none of them seem to fix the coarse-graining.

Is there any way I can attempt to fix this issue?

Thank you in advance.

Answer 1

Treat the maximum machine number as a singularity:

ListPlot[Table[{a, 
   NIntegrate[x Exp[-(a^2 + 0.001^2) x^2],
 {x, 0, 0.5 Log[$MaxMachineNumber/(a^2 + 0.001^2)], 8000}]}, {a, 0.001, 1, 1/100}], 
 Joined -> True

---

Update

[Forgive me, I actually have a job, and, while I could solve the problem quickly over breakfast, I could not compose a complete answer.]

The problem has to do with the sampling over the domain {x, 0, 8000}. The exponential function decays rather quickly for larger values of a, to the point that it is nearly zero at all initial sample points. This makes NIntegrate think the integral is zero, whether the GlobalAdaptive or LocalAdaptive strategy is chosen. For instance, the maximum value is around x == 1/(Sqrt[2] a), which is extremely close to 0 compared to 8000, when a is close to 1. The initial sampling for the Gauss-Kronrod rule is

8000 First@NIntegrate`GaussKronrodRuleData[5, MachinePrecision]
(*
  {63.6586, 375.281, 983.333, 1846.12, 2881.48, 4000.,
   5118.52, 6153.88, 7016.67, 7624.72, 7936.34}
*)

The trick is to get NIntegrate to sample nearer x == 0 than x == 63.6586.

The code

NIntegrate[f[x], {x, x0, x1, x2, ..., xn}]

effectively partitions the interval {x0, xn} into subintervals {x0, x1}, {x1, x2}, etc., each of which NIntegrate will sample independently according to the Gauss-Kronrod rule. This is one way to specify singularities, as it will ensure that NIntegrate finds the points x1, x2, etc. in calculating the integral; see tutorial/NIntegrateIntegrationStrategies for more. In this case we can use to ensure sampling the region where the integrand has strong support.

I somewhat hastily chose 1/$MaxMachineNumber for the cut-off value of the function Exp[-(a^2 + 0.001^2) x^2] (i.e. x == 0.5 Log[$MaxMachineNumber/(a^2 + 0.001^2)]), since it is about 10^-308. When I tested the code, this turned out to be a good choice. But it's not the only choice. Any number between 1/10 and 10 times this number produces an accurate plot.

Given that the maximum of the integrand is around x == 1/(Sqrt[2] a) or 0.71 / a, a simpler choice, would be something like 1 / a, which would range in value from 1000 down to 1 over the range for a in the Table. So the NIntegrate call would be the slightly simpler

NIntegrate[x Exp[-(a^2 + 0.001^2) x^2], {x, 0, 1/a, 8000}]

It can be hard to find a single approach that works for functions in general. I suppose one guideline would be to make sure the maxima (of the absolute value) get sampled. If these extrema occur at a1, a2, etc., then the following would ensure they are not missed:

NIntegrate[x Exp[-(a^2 + 0.001^2) x^2], {x, 0, a1, a2, ..., 8000}]

Alternatively, I often check a test calculation by bumping up the working precision to 20 or 30. This not only uses arbitrary-precision numbers but also bumps up PrecisionGoal. In fact, it fixes the plot in this case. The following produces the plot above:

ListPlot[Table[{a, 
   NIntegrate[x Exp[-(a^2 + (0.001`20)^2) x^2], {x, 0, 8000}, 
    WorkingPrecision -> 20]}, {a, 0.001`20, 1, 1/100}], 
 Joined -> True]

Answer 2

To illustrate whats going on, your function is 0 at x=0, rises to a max and becomes essentially zero very quickly.

 With[{a = .9}, Plot[x Exp[-(a^2+.001^2) x^2], {x, 0, 3}]]

Now look at the values NIntegrate computes:

 res = Reap[With[{a = .9},
         NIntegrate[y = x Exp[-(a^2 + 0.001^2) x^2], {x, 0, 8000},
         EvaluationMonitor :> Sow[{x, y}]]]];
         Sort@Last@Last@res // MatrixForm

The lowest x it even tries is too large, so it sees essentially all zeros and concludes the result is zero.

You can 'brute force' fix this with MinRecursion->9 for example.

Now look at what happens with Michaels's extra values in the x list:

 res = Reap[With[{a = .9},
      NIntegrate[y = x Exp[-(a^2 + 0.001^2) x^2], {x, 0, 300, 8000},
      EvaluationMonitor :> Sow[{x, y}]]]];
             ListLogPlot[(Last@Last@res)[[All, 1]]]

now plot the evaluation x points in order of evaluation:

 ListLogPlot[(Last@Last@res)[[All, 1]]]

You see it forces eval at that point and evenly spaced points in each interval. The nonzero part of the curve is picked up and the adaptive scheme takes over well.

Finally, note if you actually use Infinity as the upper limit, the initial sampling is spaced exponentially and you get the correct result with no fuss..

 res = Reap[
    With[{a = .9}, 
      NIntegrate[y = x Exp[-(a^2 + 0.001^2) x^2], {x, 0, Infinity}, 
         EvaluationMonitor :> Sow[{x, y}]]]];