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]
Integrate
? Then there are no numerics issues with substituting values in for a and b. $\endgroup$b=0
$\endgroup$