I have a problem with NIntegrate that I do not understand. I want to integrate a function f[x], which has a complicated analytic expression, but appears as a sum of narrow Lorentzian peaks. The function is defined so that $0 \leq f(x) \leq 1$ for all $x$:

Plot[f[x], {x, -0.9, 0}, PlotRange -> All, PlotLabel -> "f[x]"]


If I now integrate the function numerically in the plotted interval, I get an unrealistically high result:

NIntegrate[f[x], {x, -0.9, 0}, Method -> "MonteCarlo"]

NIntegrate::maxp: The integral failed to converge after 50100 integrand evaluations. NIntegrate obtained 3.111064646506622`*^24 and 9.60509125990105`*^22 for the integral and error estimates. >>


I tried to play around with different Nintegration methods, but the same order of magnitide shows up as a result. After some inspection, I found that the large contribution is located around $x= -0.81$. But here, $f(x)$ takes tiny values.

Does anyone know why I get this huge value? The same behavior of the function shows up when I use immediate definition:

g[x_] = f[x];
Plot[g[x], {x, -0.825, -0.8}, PlotPoints -> 150, PlotLabel -> "g[x]"]


  • 6
    $\begingroup$ I think you should share the definiton of your function too. $\endgroup$ – thedude Feb 15 '16 at 16:45
  • 3
    $\begingroup$ Why exactly do you insist on using Monte Carlo for a one-dimensional integral? There certainly are better methods of proceeding. $\endgroup$ – J. M.'s ennui Feb 15 '16 at 17:03
  • 2
    $\begingroup$ Hard to guess at why an unknown function evaluates in a way that is unexpected. $\endgroup$ – Daniel Lichtblau Feb 15 '16 at 17:58

The NIntegrate behavior you observe over your integrand is described and explained in the section "Examples of Pathological Behavior" in the chapter "NIntegrate Integration Rules" of the Advanced NIntegrate documentation .

Using larger values for the options MinRecursion, MaxRecursion, and PrecisionGoal would probably fix your issue.

NIntegrate[f[x], {x, 0, 1}, PrecisionGoal -> 12, MinRecursion -> 12, 
 MaxRecursion -> 20]

Longer answer with code

The estimates and plot shown below illustrate how the overestimation can happen.

enter image description here

The plot shows the integrand and the polynomials corresponding to the abscissas and weights of the Gauss-Kronrod integration rule.

Below is the code for the image above. I used the code from the cited document and an integrand function generated in a fashion similar to that of the answer of Dr. belisarius.

f[x_] := 318.31/(1 + 1.*^6*(-0.49158 + x)^2) + 
  318.31/(1 + 1.*^6*(-0.478 + x)^2)

Plot[f[x], {x, 0.3, 0.8}, PerformanceGoal -> "Quality", 
 PlotPoints -> 10000, PlotRange -> All]

{absc, weights, errweights} = 
  NIntegrate`GaussKronrodRuleData[12, MachinePrecision];

FitPlots[f_, {a_, b_}, abscArg_] := 
  Module[{absc = 
     Rescale[abscArg, {0, 1}, {a, 
       b}]},(*this finds the interpolating polynomial through the \
Gauss abscissas and the values of f over them*)
   polyGauss[x_] := 
      Transpose[{Take[absc, {2, -2, 2}], 
        f[#1] & /@ (Take[absc, {2, -2, 2}])}], x]];
   (*this finds the interpolating polynomial through the Gauss-
   Kronrod abscissas and the values of f over them*)
   polyGaussKronrod[x_] := 
     InterpolatingPolynomial[Transpose[{absc, f[#1] & /@ absc}], x]];
   (*plot of the Gauss interpolating points*)
   samplPointsGauss = 
    Graphics[{GrayLevel[0], PointSize[0.02], 
      Point /@ 
       Transpose[{Take[absc, {2, -2, 2}], 
         f[#1] & /@ Take[absc, {2, -2, 2}]}]}];
   (*plot of the Gauss-Kronrod interpolating points*)
   samplPointsGaussKronrod = 
    Graphics[{Red, PointSize[0.012], 
      Point /@ Transpose[{absc, f[#1] & /@ absc}]}];
   (*interpolating polynomials and f plots*)

   Block[{$DisplayFunction = Identity},
    funcPlots = Plot[{
        Tooltip[polyGauss[x], "Gauss"],
        Tooltip[polyGaussKronrod[x], "Gauss-Kronrod"],
        Tooltip[f[x], "Integrand"]}, {x, a, b},
       PlotRange -> All, PerformanceGoal -> "Quality", 
       PlotPoints -> 10000, 
       PlotLegends -> {"Gauss", "Gauss-Kronrod", "Integrand"}];];
   exact = Integrate[f[x], {x, a, b}];
   r1 = Integrate[polyGauss[x], {x, a, b}];
   r2 = Integrate[polyGaussKronrod[x], {x, a, b}];
   Print["estimated integral:" <> ToString@InputForm@r2, 
    "  exact integral:" <> ToString@InputForm@Re@exact];
   Print["estimated error:" <> ToString@InputForm@Abs[r1 - r2], 
    "  actual error:" <> ToString@InputForm@Abs[r2 - exact]];
   Show[{funcPlots, samplPointsGauss, samplPointsGaussKronrod}]

grRes = FitPlots[f[#1] &, {0.4, 0.55}, absc]

You may try for example:

(* Lorentzians*)

f[x_] := Tr@Evaluate@Table[PDF[CauchyDistribution[a, .0001], x], 
                           {a, RandomReal[{.1, .9}, 10]}]

Plot[Evaluate@f[x], {x, 0, 1}, Filling -> Axis, PlotRange -> All]

Mathematica graphics

NIntegrate[f[x], {x, 0, 1}, Method -> "LocalAdaptive"]
(* 9.99831 *) (*Quite good,should be 10 *)

Another possibility for integrating functions with wildly varying features is to avoid NIntegrate[] altogether, and consider reformulating as an initial value problem that can then be handled by NDSolve[]. Since DE solvers are often more careful than quadrature routines, they are less likely to miss sharp peaks or other features that usually give quadratures trouble.


BlockRandom[SeedRandom[0]; (* bel's/Michael's example *)
            f[x_] = Sum[PDF[CauchyDistribution[a, .0001], x],
                        {a, RandomReal[{.1, .9}, 10]}]];

AbsoluteTiming[NDSolveValue[{y'[x] == f[x], y[0] == 0}, y[1], {x, 1, 1},
                            Method -> "Extrapolation"]]
   {0.118919, 9.998274818835679}

(Thanks to Michael for pointing out a more compact way to evaluate the resulting InterpolatingFunction[] at the terminal value.)

  • $\begingroup$ +1 (earlier). By way of comparison with my timing, I get {0.015616, 9.998275411135387} with your code on my machine -- much faster in this case (V10.3.1, Macbook Pro). I also got a slightly different value for the integral. My experience is that NIntegrate is usually slower but more accurate. But a few extra digits probably doesn't matter in scientific applications (NDSolve gets almost 8, NIntegrate almost 12). Finally you can obtain the value more directly with NDSolveValue[..., y[1], {x, 1, 1},...]. $\endgroup$ – Michael E2 Feb 16 '16 at 13:06
  • $\begingroup$ @Michael, I hadn't adjusted AccuracyGoal/PrecisionGoal tho, so anyone who wants to squeeze out a few more digits would want to adjust those. However, as you say, 8-12 digits might be more than enough in most empirical applications. $\endgroup$ – J. M.'s ennui Feb 16 '16 at 13:33
  • $\begingroup$ Maybe I'm wrong, and it's a bit like comparing apples and oranges. I recall a couple of cases where I couldn't achieve the precision goal with NDSolve as easily as with NIntegrate (which will automatically extend the precision of the numbers when the going gets tough). Maybe they were edge cases. Seems to be different in this case. Perhaps I got the wrong general impression. $\endgroup$ – Michael E2 Feb 16 '16 at 14:12

If the function is like that of @Dr. b., one might divide the interval up programmatically, to help NIntegrate find the important regions and still use a "GlobalAdaptive" strategy.

Dr. belisarius's Lorentzians:

SeedRandom[0];               (* for reproducibility *)
f[x_] := Evaluate@Tr@
  Table[PDF[CauchyDistribution[a, .0001], x], {a, RandomReal[{.1, .9}, 10]}]        

Plot[Evaluate@f[x], {x, 0, 1}, Filling -> Axis, PlotRange -> All]

Mathematica graphics

We divide the integration regions at the local maxima, up to a tolerance (10^-12 below, but it should depend on the interval of integration, as well as perhaps f). The function new tries to find a new maximum at a starting point, usually the midpoint between two maxima. The function subdivide subdivides the interval if there is a local maximum between the end points, provided the new maximum is significantly different from the end points (governed by the tolerance 1*^-12 in this case).

new = FindArgMax[f[x], {x, #}, Method -> {"Newton", "StepControl" -> "TrustRegion"}] &;

subdivide[{a_, b_}] := With[{c = First@new[(a + b)/2]},
  If[a + 1*^-12 < c < b - 1*^-12,
   {subdivide[{a, c}], subdivide[{c, b}]},
   {a, b}]

The integral evaluates without difficulty.

   Sequence @@ Cases[subdivide[{0., 1.}], {a_Real, b_Real} :> a, Infinity],
   1.}] // AbsoluteTiming
(*  {0.30745, 9.99827}  *)

The division points, for the curious:

Cases[subdivide[{0., 1.}], {a_Real, b_Real} :> a, Infinity]
  {0., 0.180879, 0.290761, 0.553081, 0.606456, 0.61005,
   0.61642, 0.621974, 0.64625, 0.848162, 0.88095}

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.