8
$\begingroup$

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}
Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ If you change the integration method you should have better luck. Try this: NIntegrate[function[s, t], {t, 0, 1/2}, {s, 0, 1/2}, Method -> "DoubleExponential"] $\endgroup$
    – leibs
    Commented Jan 27, 2014 at 20:11
  • 1
    $\begingroup$ You can reduce the absolute precision (default is Infinity). For example, with the option AccuracyGoal->5 it is 1000 times faster and there are no warning messages $\endgroup$
    – andre314
    Commented Jan 27, 2014 at 20:25
  • 1
    $\begingroup$ A sidenote, I believe polynomials can always be handled by Integrate very efficiently. $\endgroup$
    – Silvia
    Commented Jan 27, 2014 at 20:46

7 Answers 7

Reset to default
9
$\begingroup$

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.

Improve this answer
$\endgroup$
6
$\begingroup$

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

Improve this answer
$\endgroup$
1
  • $\begingroup$ great, thanks a lot $\endgroup$
    – sasa
    Commented Jan 27, 2014 at 23:30
5
$\begingroup$ 
AbsoluteTiming[
 NIntegrate[function[s, t], {t, 0, 1/2}, {s, 0, 1/2}, 
            PrecisionGoal -> 11, Method -> "LocalAdaptive"]]
(*
  {0.671875, -1.77636*10^-14}
*)
Improve this answer
$\endgroup$
3
  • $\begingroup$ I wonder why "LocalAdaptive" strategy has a better performance here. $\endgroup$
    – M6299
    Commented Jan 28, 2014 at 11:02
  • $\begingroup$ @M6299 From the docs (perhaps it helps): When "LocalAdaptive" is faster and performs better than "GlobalAdaptive", it is because the precision-goal-stopping criteria and partitioning strategy of "LocalAdaptive" are more suited for the integrand's nature. Another factor is the ability of "LocalAdaptive" to reuse the integral values of all points already sampled. "GlobalAdaptive" has the ability to reuse very few integral values (at most 3 per rule application, 0 for the default one-dimensional rule, the Gauss-Kronrod rule). $\endgroup$
    – Dr. belisarius
    Commented Jan 28, 2014 at 11:12
  • 5
    $\begingroup$ Unfortunately, this is actually an example of LocalAdaptive failing: It incorrectly believes it has converged to a precision of 11. But the correct answer is zero, so any finite answer (such as -1.77636*10^-14) is incorrect to any precision. $\endgroup$
    – Andrew Moylan
    Commented Jan 28, 2014 at 23:53
4
$\begingroup$

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>>]}
Improve this answer
$\endgroup$
4
$\begingroup$

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".

Improve this answer
$\endgroup$
3
$\begingroup$

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}
Improve this answer
$\endgroup$
1
  • $\begingroup$ (I changed the example a bit because a zero result here is a bit uninteresting.) $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Oct 18, 2018 at 18:15
1
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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