# Mathematica's NIntegrate is so slower compare to scipy

In this post, Antonov said Mathematica is better than scipy in integration.

However, I've tried several examples, Mathematica is always several times slower than scipy.

For example,

In[1]:= NIntegrate[x y + Sin[x^y], {x, 0, 4}, {y, 0, 4}, PrecisionGoal -> 8,
AccuracyGoal -> 8] // AbsoluteTiming

Out[1]= {2.98762, 68.8117}

In[2]:= NIntegrate[x y + Sin[x^y], {x, 0, 4}, {y, 0, 4}] // AbsoluteTiming

Out[2]= {2.10989, 68.8117}


According to this post, scipy's integration corresponds to PrecisionGoal -> 8,AccuracyGoal -> 8

Now look at scipy.

In [1]: import scipy.integrate

In [2]: import numpy as np

In [3]: %time scipy.integrate.dblquad(lambda x,y:x*y+np.sin(x**y), 0, 4,lambda x
...: :0,lambda x:4)
CPU times: user 593 ms, sys: 45.8 ms, total: 639 ms
Wall time: 638 ms
Out[3]: (68.8116996894142, 4.805091491642489e-07)


It is almost 5 times faster than 2.98762s.

What magic made scipy this fast? Any idea to make NIntegration faster?

What I actually want to is to plot integration result like this

Plot[NIntegrate[-2 Im[((-0.0006250000000000001 + ((3.5 +
0.02 I) + \[Omega] - 1.9 Cos[kx] -
2.1 Cos[ky]) ((3.5 + 0.02 I) + \[Omega] -
2.1 Cos[kx] - 1.9 Cos[ky])) ((-3.5 +
0.02 I) + \[Omega] + 2.1 Cos[kx] +
1.9 Cos[
ky]))/(-0.0006250000000000001 (-0.0006250000000000001 + \
((-3.5 + 0.02 I) + \[Omega] + 2.1 Cos[kx] +
1.9 Cos[ky]) ((-3.5 + 0.02 I) + \[Omega] +
1.9 Cos[kx] + 2.1 Cos[ky])) + ((3.5 +
0.02 I) + \[Omega] - 1.9 Cos[kx] -
2.1 Cos[ky]) ((3.5 + 0.02 I) + \[Omega] - 2.1 Cos[kx] -
1.9 Cos[
ky]) (-0.0006250000000000001 + ((-3.5 +
0.02 I) + \[Omega] + 2.1 Cos[kx] +
1.9 Cos[ky]) ((-3.5 + 0.02 I) + \[Omega] +
1.9 Cos[kx] +
2.1 Cos[
ky])))], {kx, -\[Pi], \[Pi]}, {ky, -\[Pi], \[Pi]}], {\
\[Omega].-0.5, 0.5}]


each integration took more than 3 seconds, the plotting is slow. Do I really need to call scipy inside Mathematica. This way is awkward, and there is no standard conversion function to conveniently convert long expression to python form.

• Did you call scipy from MMA? If so, could you post the full code? Mar 6 '18 at 15:42
• @anderstood no, I run scipy in ipython Mar 6 '18 at 15:44
• NIntegrate[Integrate[x y+Sin[x^y],{y,0,4}], {x,0,4}, PrecisionGoal->8, AccuracyGoal->8]//AbsoluteTiming Sep 20 '18 at 6:49

There are various ways to fine tune NIntegrate. The difference of the following result to scipy's is of order 10^-7 and it needs less than a thenth of the computation time of NIntegrate without specified Method.

a = 68.811699689414;
{t, b} = NIntegrate[x y + Sin[x^y], {x, 0, 4}, {y, 0, 4},
PrecisionGoal -> 8, AccuracyGoal -> 8,
Method -> {"GaussKronrodRule", "Points" -> 9}}
] // AbsoluteTiming
a - b


{0.176262, 68.81169993374982}

-2.44336*10^-7

• Thank you so much! Seems great. Your fine tuning also made my last complex integration 3 times faster. But this kind of find tuning requires so much knowledge... And I am curious, why Automatic option is so stupid, it has symbolicprocessing! Doesn't help to automatically decide a better method? Mar 6 '18 at 15:50
• Well, NIntegrate is build as a swiss army knive that is able to decide which tool to use. That requires some overhead. In this case, I know that the function is very oscillatory in a localized region, so "LocalAdaptive" might be a good idea. On the other hand, the function is very smooth (it is even analytic), so high order Gauss quadratures should help. The latter are less useful when there is a discontinuity or a pole in the integrant. Mar 6 '18 at 15:54
• Without any a priori knowledge, Mathematica has to analyze the integrant and perform some checks to chose an integration strategy. This choice may be however suboptimal in certain cases... Mar 6 '18 at 15:56
• @HenrikSchumacher Kudos for being polite and constructive in your answer and comments, given OP's attitude. Mar 6 '18 at 15:57
• @AntonAntonov I am so sorry. I am absolutely not intend to be offensive. Sorry, sorry Mar 6 '18 at 16:07