Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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}
share|improve this question
1  
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"] –  leibs Jan 27 at 20:11
1  
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 –  andre Jan 27 at 20:25
    
A sidenote, I believe polynomials can always be handled by Integrate very efficiently. –  Silvia Jan 27 at 20:46
add comment

5 Answers 5

AbsoluteTiming[
 NIntegrate[function[s, t], {t, 0, 1/2}, {s, 0, 1/2}, 
            PrecisionGoal -> 11, Method -> "LocalAdaptive"]]
(*
  {0.671875, -1.77636*10^-14}
*)
share|improve this answer
    
I wonder why "LocalAdaptive" strategy has a better performance here. –  M6299 Jan 28 at 11:02
    
@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). –  belisarius Jan 28 at 11:12
    
Thank you so much. –  M6299 Jan 28 at 14:44
    
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. –  Andrew Moylan Jan 28 at 23:53
add comment

For smooth functions like the one you are integrating, I often find it can be much faster to use Guassian 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]x[a,b].

share|improve this answer
    
great, thanks a lot –  sasa Jan 27 at 23:30
add comment

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>>]}
share|improve this answer
add comment

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.

share|improve this answer
add comment

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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