# 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? – anderstood Mar 6 '18 at 15:42
• @anderstood no, I run scipy in ipython – matheorem Mar 6 '18 at 15:44
• NIntegrate[Integrate[x y+Sin[x^y],{y,0,4}], {x,0,4}, PrecisionGoal->8, AccuracyGoal->8]//AbsoluteTiming – matrix89 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? – matheorem 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. – Henrik Schumacher 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... – Henrik Schumacher Mar 6 '18 at 15:56
• @HenrikSchumacher Kudos for being polite and constructive in your answer and comments, given OP's attitude. – Anton Antonov Mar 6 '18 at 15:57
• @AntonAntonov I am so sorry. I am absolutely not intend to be offensive. Sorry, sorry – matheorem Mar 6 '18 at 16:07