# Match NIntegrate vs Integrate with HighPrecision

This a pretty basic question but I am not sure if the answer is simple.

I have a function n-dimensional and lets assume that I want to integrate with 50 digit accuracy . Integrate is easy but I am failing to match NIntegrate with Integrate. Actually, I am failing to get a converged result with NIntegrate.

AbsoluteTiming[
N[Integrate[
z^2 Sin[π x] + I z^4 Cos[π y] + I y^6 Cos[2 π y] +
w Sin [π y^2], {x, -12/10, 1}, {y, 1, 2}, {z, -21/10,
1}, {w, -3, 3}], 30]]
AbsoluteTiming[
NIntegrate[
z^2 Sin[π x] + I z^4 Cos[π y] + I y^6 Cos[2 π y] +
w Sin [π y^2], {x, -12/10, 1}, {y, 1, 2}, {z, -21/10,
1}, {w, -3, 3}, WorkingPrecision -> 30, PrecisionGoal -> 30]]


EDIT 1

Output of the NIntegrate

The global error of the strategy GlobalAdaptive has increased more \ than 2000 times. The global error is expected to decrease \ monotonically after a number of integrand evaluations. Suspect one of \ the following: the working precision is insufficient for the \ specified precision goal; the integrand is highly oscillatory or it \ is not a (piecewise) smooth function; or the true value of the \ integral is 0. Increasing the value of the GlobalAdaptive option \ MaxErrorIncreases might lead to a convergent numerical integration. \

And result was matching only up to 6 digits.

• What was your output? – Michael E2 Aug 6 '16 at 21:37

First, let's compute a somewhat higher precision approximation to the exact integral for the sake of comparison with the numerical approximation below:

exact[50] =
N[Integrate[
z^2 Sin[π x] + I z^4 Cos[π y] + I y^6 Cos[2 π y] +
w Sin[π y^2], {x, -12/10, 1}, {y, 1, 2}, {z, -21/10,
1}, {w, -3, 3}], 50];


High dimensional integrals are reputed to be hard to numerically compute. This one has a pretty easy integrand (at least over the domain). You can get it rather accurately with one go at it with a sufficiently high-order Gauss-Kronrod rule.

AbsoluteTiming[
approx[11] =
NIntegrate[
z^2 Sin[π x] + I z^4 Cos[π y] + I y^6 Cos[2 π y] +
w Sin[π y^2], {x, -12/10, 1}, {y, 1, 2}, {z, -21/10,
1}, {w, -3, 3},
Method -> {"GaussKronrodRule", "Points" -> 11}, WorkingPrecision -> 50,
PrecisionGoal -> 30, MaxRecursion -> 0]]
(* ignore the warning
{18.4073,
1.2475688841959872224585179766063977904579143591081 +
171.21635046746085306563299746868711863872735418560 I}
*)


It has greater than 30 digits of accuracy:

exact[50] - approx[11]
(*  1.9088394186586*10^-34 + 1.0196130029198730*10^-31 I  *)


If one thinks about the order of the GK rule and the integrand, one can figure out a faster way:

AbsoluteTiming[
gkr = Quiet@
NIntegrate[
z^2 Sin[π x] + I z^4 Cos[π y] + I y^6 Cos[2 π y] +
w Sin[π y^2], {x, -12/10, 1}, {y, 1, 2}, {z, -21/10,
1}, {w, -3, 3},
Method -> {"CartesianRule", Method -> {
{"GaussKronrodRule", "Points" -> 11},
{"GaussKronrodRule", "Points" -> 11},
{"GaussKronrodRule", "Points" -> 3},
{"GaussKronrodRule", "Points" -> 3}}},
WorkingPrecision -> 50, PrecisionGoal -> 30, MaxRecursion -> 0]]
exact[50] - gkr
(*
{1.69512,
1.2475688841959872224585179766063977904579143591081 +
171.21635046746085306563299746868711863872735418560 I}
1.9088394186586*10^-34 + 1.0196130029198730*10^-31 I
*)


Update: Some of the theory behind the answer

First, when I described the integral as "easy," I meant the integrand was (1) analytic over a neighborhood of the interval of integration ("analytic" = "represented its power series", not the numerical analysis meaning of "symbolic") and (2) does not oscillate much over the interval. The remarks below concern such functions and integrals.

1. The default "MultidimensionalRule" does rather sparse sampling with poor accuracy - a compromise between speed and accuracy. For a high dimensional integral, most users are satisfied with a few digits of precision. In such a case it is a good default method. On the other hand, the "CartesianRule" uses a tensor-product grid. (The "CartesianRule" is also applied when a one-dimensional method is indicated, such as Method -> "GaussKronrodRule".) Such a grid is relatively dense. High precision should be possible, but sampling over it, especially a 4D one, is expensive. If recursion can be minimized or avoided and the grid is not too large, it can be effective. That in turn will depend on the convergence rate. For an analytic function, the "GaussKronrodRule" (as well as for "GaussBerntsenEspelidRule", "LobattoKronrodRule" and "ClenshawCurtisRule") ultimately converges exponentially as the number of "Points" increases.

The plot below shows the logarithm of the absolute error of the Gauss-Kronrod rule for an increasing n in the setting of the option "Points" -> n, in the graphic below for the 1D integral of 50 Sin[1 + 3 x/2]^2 Cos[x]^2 over {x, 0, 10}, which has 6 maxima and 6 minima over the interval. Keep in mind the number of sampling nodes is 2n + 1 and has an order of 3n - 1 or 3n - 2 depending on whether n is odd or even. Normally, then, to get past the pre-convergent domain, you would expect the order to have to be at least the number of inflection points plus one, and in practice it usually needs to be a bit greater.

There are three kinds of behavior. In the first,"pre-convergent" domain (n < 7), the error bounce around their maximum values. In the second, "exponential convergence" domain (7 < n < 17), we see the steady decline of the error. In the third, "precision limit" domain (n > 17), the decline stops, due to the limitations of the working precision. For machine precision the range from the maximum to the minimum is roughly 16 (digits), with exceptions typically due to random rounding error.

The take-away point is that the real advantage of the rule comes into play when the "Points" land in the exponential convergence domain. Frequently, in my experience, for an analytic function, the default "Points" -> 5 is not quite high enough to reach exponential convergence for a given integrand. Doubling the setting to "Points" -> 9 or 11 improves the performance. (On the other hand, when the function is not analytic and there is not exponential convergence, it can be more efficient to subdivide.) The example above doesn't enter exponential convergence until the number of points is 6 or 7, but it is contrived to have an "obviously" pre-convergent phase. More important is that it doesn't reach single precision (~8 digits) until n is at least 12. In this region, we see that the number of digits of precision increases linearly with n. If n is not already too big, then it would be more efficient to increase n than to use recursive subdivision. Except when integrating polynomials of degree less than 15, the rule is rarely near the end of exponential convergence with the default "Points" -> 5. A rule of thumb for analytic functions is that doubling the setting to "Points" -> 9 or 11 usually improves performance.

2. Another complication is that the cost of increasing n is multiplied by the numbers of sample nodes in the other dimensions. This prospect ought to be discouraging, but in this case, we were lucky. As a function of each of z and w, the integrand is a polynomial of degree 4 and 1 respectively. So "Points" -> 2 is sufficient (I used 3 thoughtlessly above, out of habit, because it's odd). That makes increasing the points in the difficult dimensions less costly. In fact, the integral over w cancels out the Sin[π y^2] term, which reduces the need for more points.

3. As @bbgodfrey noticed, MaxRecursion -> 0 is critical to speed. The error estimator tends to overestimate the error. MaxRecursion -> 0 prevents subdivision even if the estimated error calls for it. I realized later that the Gauss rule, "GaussBerntsenEspelidRule", would perform better than Gauss-Kronrod. It uses n sample points for "Points" -> n and as an order of 2n - 1. The ratio of order to points is roughly 2:1 while for GK, it is 3:2. The method

Method -> {"CartesianRule", Method -> {
{"GaussBerntsenEspelidRule", "Points" -> 8},
{"GaussBerntsenEspelidRule", "Points" -> 9},
{"GaussBerntsenEspelidRule", "Points" -> 2},
{"GaussBerntsenEspelidRule", "Points" -> 2}}}


produces an answer in 0.543138 sec., with a precision of more than 31 digits. One can estimate these numbers for "Points" from the convergence of the one-dimensional integrals:

Table[
With[{iter = i},
Block[{x, y, z, w},
x = RandomReal[{-12/10, 1}];
y = RandomReal[{1, 2}];
z = RandomReal[{-21/10, 1}];
w = 0;(*w=RandomReal[{-3,3}];*)   (* the integral of w cancels out! *)
-1 + First@NestWhile[
{1 + First@#,
Quiet@NIntegrate[
z^2 Sin[π x] + I z^4 Cos[π y] + I y^6 Cos[2 π y] + w Sin[π y^2],
iter,
Method -> {"GaussBerntsenEspelidRule", "Points" -> 1 + First[#]},
WorkingPrecision -> 50, MaxRecursion -> 0]} &,
{1, 1},
RealExponent[{1, -1}.{##}[[All, 2]]/Last[#2]] > -30 &,
2,
20
]]
],
{i, {{x, -12/10, 1}, {y, 1, 2}, {z, -21/10, 1}, {w, -3, 3}}}] // AbsoluteTiming
(*
{0.065465, {8, 9, 2, 2}}
*)


Note that it only probably gives good settings.

• @bbgodfrey The default (Automatic) "MultidimensionalRule" does rather sparse sampling with poor accuracy - a compromise between speed and accuracy. (For a high dimensional integral, most are satisfied with a few digits of precision.) OTOH, the "CartesianRule", which is also invoked in your Method -> "GaussKronrodRule" (+1, btw), uses a tensor-product grid. This is expensive, but it can be effective, if recursion can be minimized or avoided. By increasing "Points", it can be avoided, except that the error estimator overestimates the error. MaxRecursion -> 1 overrides it.... – Michael E2 Aug 7 '16 at 1:41
• ...The default for "GaussKronrodRule" has "Points" -> 5, which has 11 sample points (per dimension) and gives a rule of order 17 (zero error up to degree 17), so it probably needs a few recursive subdivisions to reach the desired precision (at least in the x and y dimensions, given the overestimation of the error). – Michael E2 Aug 7 '16 at 1:47
• Oops, I meant "MaxRecursion -> 0 overrides it", of course. – Michael E2 Aug 7 '16 at 1:52
• @bbgodfrey Another thing - it might even be a rule of thumb - I've noticed: The Gauss-Kronrod (& Clenshaw-Curtis) rule converges exponentially as the number of "Points" increases for analytic functions ("analytic" = "equals its power series"). The default setting often is not high enough to reach exponential convergence for a given integrand. Frequently for analytic functions, doubling the setting to "Points" -> 9 or 11 improves the performance. (By habit, I pick odd numbers. In fact you can halve the time by setting "Points" -> 2 for z & w above & get the same answer.) – Michael E2 Aug 7 '16 at 11:23
• @bbgodfrey Done. And then some. :) Perhaps overkill. But in fact, I sometimes find it easier to locate my code by searching SE than searching my computer. – Michael E2 Aug 8 '16 at 16:00

The following converges:

AbsoluteTiming[
NIntegrate[z^2 Sin[π x] + I z^4 Cos[π y] + I y^6 Cos[2 π y] +
w Sin[π y^2] + w Sin[π y^2], {x, -12/10, 1}, {y, 1, 2}, {z, -21/10, 1},
{w, -3, 3}, WorkingPrecision -> 60, PrecisionGoal -> 30,
Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 10000,
Method -> "GaussKronrodRule"}]]

(* {159.524,
1.24756888419598722245851797660639798134185622496647479523567 +
171.216350467460853065632997468687220600027646172902541171176 I} *)


which agrees well with the Integrate result

AbsoluteTiming[
N[Integrate[z^2 Sin[π x] + I z^4 Cos[π y] + I y^6 Cos[2 π y] +
w Sin[π y^2], {x, -12/10, 1}, {y, 1, 2}, {z, -21/10, 1}, {w, -3, 3}], 70]]

(* {1.64704,
1.24756888419598722245851797660639798134185622496647479523567 +
171.216350467460853065632997468687220600027646172902541171176 I} *)

Abs[Last@% - Last@%%]
(* 0.*10^-58 *)


As suggested in the error message from NIntegrate as used in the question, increasing "MaxErrorIncreases" solves the problem. See NIntegrate Integration Strategies for further discussion.