NIntegrate::slwcon Problem

Question

I have a problem with numerical integration of this function. Integral value is zero, but NIntegrate[] needs a lot of time to calculate this. Is there any way to speed up this calculation?

Input:

function[s_, t_] :=100 (-2160 (1 - 2 s)^4 t^3 (-2 + 5 t) + 
 96 (1 - 2 s) t^3 (25 (-1 + 2 s) t^2 (-3 + 5 t) + 
    5/4 (-1 + 2 s)^3 (-1 + 5 t)) - 
 24 (-1 + 2 s)^3 t (5 (-1 + 2 s) t^2 (-3 + 5 t) + 
    25/4 (-1 + 2 s)^3 (-1 + 5 t)));

AbsoluteTiming[NIntegrate[function[s, t], {t, 0, 1/2}, {s, 0, 1/2}]]
AbsoluteTiming[Integrate[function[s, t], {t, 0, 1/2}, {s, 0, 1/2}]]

Output:

During evaluation of In[4]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>

During evaluation of In[4]:= NIntegrate::eincr: 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. NIntegrate obtained -1.49243*10^-14 and 1.0093478121591215`*^-12 for the integral and error estimates. >>

Out[5]= {43.421484, -1.49243*10^-14}

Out[6]= {3.963227, 0}

Answer 1

As the warning message you received suggests, this issue is caused by: "... the true value of the integral is 0".

When the true value is zero, the default PrecisionGoal -> 6 can never be satisfied. You need to set a finite AccuracyGoal in such cases:

NIntegrate[function[s, t], {t, 0, 1/2}, {s, 0, 1/2}, 
 AccuracyGoal -> 5]

gives -2.11892*10^-9.

A Cartesian product of Gaussian quadrature rules seems to be much faster than the default sparse multidimensional rule in this case too:

NIntegrate[function[s, t], {t, 0, 1/2}, {s, 0, 1/2}, 
 AccuracyGoal -> 10, Method -> "GaussKronrodRule"]

gives -2.98428*10^-13.

Answer 2

For smooth functions like the one you are integrating, I often find it can be much faster to use Gaussian quadrature methods. I'm sure someone with lots of experience using NIntegrate could find the right options to make Mathematica automatically do something like this.

Needs["NumericalDifferentialEquationAnalysis`"];
function = Compile[{{s, _Real, 0}, {t, _Real, 0}}, 100 (-2160 (1 - 2 s)^4 t^3(-2 + 5 t) + 
96 (1 - 2 s) t^3 (25 (-1 + 2 s) t^2 (-3 + 5 t) + 
5/4 (-1 + 2 s)^3 (-1 + 5 t)) - 
24 (-1 + 2 s)^3 t (5 (-1 + 2 s) t^2 (-3 + 5 t) + 
25/4 (-1 + 2 s)^3 (-1 + 5 t)))];

a = 0;
b = 1/2;
order = 6;
{x, w} = Transpose[GaussianQuadratureWeights[order, a, b]];
Map[function @@ # &, Flatten[Outer[List, x, x], 1]].Flatten[Outer[Times, w, w]] // Timing
(*{0.000154, -6.10727*10^-14}*)

You could easily modify the code for more general rectangular regions. Here I assumed that you were integrating the function function on the square domain $[a,b]\times[a,b]$.

Answer 3

AbsoluteTiming[
 NIntegrate[function[s, t], {t, 0, 1/2}, {s, 0, 1/2}, 
            PrecisionGoal -> 11, Method -> "LocalAdaptive"]]
(*
  {0.671875, -1.77636*10^-14}
*)

Answer 4

We can amplify and combine @Silva's and @Andrew Moylan's points, which were that the integrand was a polynomial and that the Gauss-Kronrod rule performed better than the multidimensional rule.

We can use the Cartesian product of any interpolatory rules such as "ClenshawCurtisRule" or "GaussKronrodRule", provided we use a sufficiently high order that the integration rule and the embedded error rule integrate the polynomial exactly. According to the docs, for a setting of "Points" -> n, the error rule for "GaussKronrodRule" will be exact for on a polynomial of degree up to 2n - 1 and for "ClenshawCurtisRule" will be exact on a polynomial of degree upto n (but they meant up to degree n or n - 1, according as n is odd or even resp.). The polynomial f[s, t] is of degree {6, 6}, so both the s and t integrals can use the same setting for "Points".

The principal issue, already pointed out by Andrew Moylan, is that zero cannot be approximated to any Precision, and so you have to set a finite AccuracyGoal -> acc, which indicates that a numerical result closer to zero than 10^-acc will be acceptable. More precisely, it means that the error between the integration rule and the embedded rule is to be less than 10^-acc + i0 * 10^-prec, where i0 is the result of the integration rule and prec is the PrecisionGoal.

With all that in mind and sufficient WorkingPrecision, any level of accuracy is obtainable (within practical limits). For instance, to get closer than 10^-100:

Block[{acc = 100, deg = 6, res},
 (res = NIntegrate[function[s, t], {t, 0, 1/2}, {s, 0, 1/2}, 
      Method -> {"ClenshawCurtisRule", "Points" -> deg + 1}, 
      MaxRecursion -> 0, AccuracyGoal -> acc, 
      WorkingPrecision -> acc]) // AbsoluteTiming // Print;
 Chop[res, 10^-acc]
 ]
(*
  {0.00845, -6.4369978618316750815140837884620002197320415488373585727...3*10^-114}
  0
*)

Block[{acc = 100, deg = 6, res},
 (res = NIntegrate[function[s, t], {t, 0, 1/2}, {s, 0, 1/2}, 
      Method -> {"GaussKronrodRule", 
        "Points" -> Ceiling[(deg + 1)/2]}, MaxRecursion -> 0, 
      AccuracyGoal -> acc, WorkingPrecision -> acc]) // 
   AbsoluteTiming // Print;
 Chop[res, 10^-acc]
 ]
(*
  {0.006585, 3.0164732651049401161481904591662713102211101606252755851...8*10^-113}
  0
*)

To use MachinePrecision, one has to account for rounding error and subtract a few digits off the AccuracyGoal.

NIntegrate[function[s, t], {t, 0, 1/2}, {s, 0, 1/2}, 
  Method -> {"ClenshawCurtisRule", "Points" -> 7}, MaxRecursion -> 0, 
  AccuracyGoal -> MachinePrecision - 2.7, 
  WorkingPrecision -> MachinePrecision] // AbsoluteTiming
(*  {0.003191, -1.55431*10^-14}  *)

NIntegrate[function[s, t], {t, 0, 1/2}, {s, 0, 1/2}, 
  Method -> {"GaussKronrodRule", "Points" -> 4}, MaxRecursion -> 0, 
  AccuracyGoal -> MachinePrecision - 2.7, 
  WorkingPrecision -> MachinePrecision] // AbsoluteTiming
(*  {0.003304, -2.45137*10^-13}  *)

One can also use "NewtonCotesRule" and "LobattoKronrodRule". Note also that MaxRecursion -> 0 means NIntegrate makes just a single application of the rule without recursive subdivision. @leibs was wondering whether classical Gauss quadrature could be done this way, which it can by using "GaussBerntsenEspelidRule".

Answer 5

This may be a little off-topic, but as OP's question is about integral of a certain function, here is a way without NIntegrate. That is, if all your concern is integrate a polynomial, and notice the indefinite version of Integrate is usually faster than the definite one, it should be safe and convenient to go with indefinite integral, than invoke the fundamental theorem of calculus. So

Clear[polynomialIntegrate]
polynomialIntegrate[poly_, varspecs : {Repeated[_, {3}]} ..] :=
    Fold[
         Function[{expr, varspec},
                  (Integrate[expr, #1] /. {{#1 -> #2}, {#1 -> #3}}).{-1, 1} & @@ varspec
                 ],
         poly, {varspecs}]

Use it like Integrate:

polynomialIntegrate[function[s, t], {s, 0, 1/2}, {t, 0, 1/2}] // AbsoluteTiming

{0.022002, 0}

polynomialIntegrate[function[s, t], {s, t, 1 - t^2}, {t, a, b}] // AbsoluteTiming

{0.023001, -(1000/7) Plus[<<16>>] + 1000/7 Plus[<<16>>]}

Answer 6

I am writing this answer to put on record a (relatively) quick and stable method for symbolically integrating polynomials, due to Velvel Kahan. The method is based on his procedure for computing divided differences of polynomials, and can be thought of as a modified Horner scheme:

polyInt[poly_, {x_, a_, b_}] /; PolynomialQ[poly, x] := 
    Module[{n = Exponent[poly, x], d, p, y},
           p = CoefficientList[poly, x]/Range[n + 1]; d = y = 0;
           Do[y = b y + p[[k]]; d = a d + y, {k, n + 1, 1, -1}];
           (b - a) d]

This can be used as a direct substitute for Integrate[poly, {x, a, b}] in the univariate case. For the multivariate case, one can use this with Fold[]; e.g.

AbsoluteTiming[Fold[polyInt, function[s, t], {{s, t, 1 - t}, {t, 0, 1/2}}]]
   {0.00633007, 375/28}

AbsoluteTiming[Integrate[function[s, t], {t, 0, 1/2}, {s, t, 1 - t}]]
   {5.72321, 375/28}

Answer 7

Another way is to manually split the integral up. First use FullSimplify to deal with a smaller integrand:

f[s_,t_]:=15000 (1 - 2 s)^2 t ((1 - 2 s)^4 - 5 (1 - 2 s)^4 t + 
32 (1 - 2 s)^2 t^2 - 80 (1 - 2 s)^2 t^3 + 48 t^4 - 80 t^5)

AbsoluteTiming[NIntegrate[f[s, t], {t, 1/8, 1/2}, {s, 0, 1/2}] + NIntegrate[f[s, t],
{t, 0, 1/8}, {s, 0, 1/2}]]

yields

{0.109375, -1.980108521593138*10^-9}

Not as good as the other methods, but sometimes it can pay off.