2
$\begingroup$

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.

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

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 -> {"LocalAdaptive", 
     Method -> {"GaussKronrodRule", "Points" -> 9}}
   ] // AbsoluteTiming
a - b

{0.176262, 68.81169993374982}

-2.44336*10^-7

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – matheorem Mar 6 '18 at 15:50
  • 3
    $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Mar 6 '18 at 15:54
  • $\begingroup$ 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... $\endgroup$ – Henrik Schumacher Mar 6 '18 at 15:56
  • 2
    $\begingroup$ @HenrikSchumacher Kudos for being polite and constructive in your answer and comments, given OP's attitude. $\endgroup$ – Anton Antonov Mar 6 '18 at 15:57
  • 1
    $\begingroup$ @AntonAntonov I am so sorry. I am absolutely not intend to be offensive. Sorry, sorry $\endgroup$ – matheorem Mar 6 '18 at 16:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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