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.
