4
$\begingroup$

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). enter image description here

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: enter image description here

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.

$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – user9660 Apr 24 '15 at 4:24
  • $\begingroup$ In your first code-example the upper limit of integration is expressed in a Term involving b, also expression T is not given. Maybe something is missing in the code? $\endgroup$ – mgamer Apr 24 '15 at 10:41
  • $\begingroup$ Thank you mgamer, I think now everything is fine $\endgroup$ – Alex Apr 24 '15 at 13:57
  • $\begingroup$ If this is the actual you're using, why not get the symbolic result in terms of the parameters a and b with Integrate? Then there are no numerics issues with substituting values in for a and b. $\endgroup$ – Michael E2 Apr 24 '15 at 14:55
  • 1
    $\begingroup$ in the spirit of constructing minimal examples, I'll point out you get the same issue with b=0 $\endgroup$ – george2079 Apr 24 '15 at 17:23
4
$\begingroup$

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

Mathematica graphics


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]
$\endgroup$
  • $\begingroup$ Thank you so much Michael. Nevertheless, I still have two questions: First, I don't think I'm completely understanding what you are doing when inserting that fourth argument 0.5 Log[$MaxMachineNumber/(a^2 + 0.001^2)] in the integration limits part. Additionally and more as a curiosity, when trying LocalAdaptive I still get the coarse-graining problem. Do you know any reason that may prevent the LocalAdaptive method to do a successful integration? $\endgroup$ – Alex Apr 24 '15 at 16:39
  • $\begingroup$ @pepeperez I have tried to incorporate my thinking in the update. I hope to have answered your questions as well as explain my approach. $\endgroup$ – Michael E2 Apr 25 '15 at 2:54
  • $\begingroup$ This is actually a great answer full of details, thank you for the time you spent on it. $\endgroup$ – Alex Apr 25 '15 at 14:11
  • 1
    $\begingroup$ @pepeperez You're welcome. The only thing to add was touched on by george2079. A significant difference between a finite interval {0, 8000} and an infinite on {0, Infinity} in a non-oscillatory integral is this. The finite interval is treated symmetrically by the Gauss-Kronrod rule, with more or less evenly spaced sampling. The infinite interval {0, Infinity} is transformed to {0, 1} via the substitution u = 1/(1+x); Gauss-Kronrod is applied to {0, 1}, which means the sampling is denser near the finite end point x == 0. This works well with your example, which decays rapidly. $\endgroup$ – Michael E2 Apr 25 '15 at 19:58
6
$\begingroup$

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}]]

enter image description here

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

enter image description here

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]]]

enter image description here

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}]]]];

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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