# How to speed up numerical integration?

I need to calculate the integral numerically

NIntegrate[ StruveH[0, Sqrt[r1^2 + r2^2 - 2*r1*r2*Cos[\[Theta]]]] - BesselY[0, Sqrt[r1^2 + r2^2 - 2*r1*r2*Cos[\[Theta]]]], {\[Theta], 0, 2*Pi}, {r1, 0, 7}, {r2, 0, 7}] // Timing

but Mathematica takes a very long time to calculate it. Could you tell me, please, is there any way to speed up the calculation of the integral?

• You can try Method -> "LocalAdaptive" for a bit of speedup. Commented Sep 24, 2022 at 19:34
• Are you sure that your integrand is correct? Intuition tells me you're missing a Jacobian, possibly something like $8\pi^2r_1^2r_2^2\sin\theta$. If that is so, then analytic integration may be of use for speeding up the calculation. Commented Sep 24, 2022 at 19:45
• It's correct integrand Commented Sep 24, 2022 at 20:06
• If you have many such integrals, you should consider rewriting them as 1-dimensional integrals: There is a function $h(x)$ such that for all functions $f(x)$ one has $\iiint f(\sqrt{r_1^2 + r_2^2 - 2r_1r_2 \cos \theta}) dr_1 dr_2 d\theta = \int_0^{14} f(x) h(x) dx$. To do this you have to construct or compute $h$ first, but that would make sense if you need it for many $f$. Commented Sep 25, 2022 at 11:43

We can speed up with using precision and accuracy options. First, we compute without options for comparison.

NIntegrate[
StruveH[0, Sqrt[r1^2 + r2^2 - 2*r1*r2*Cos[\[Theta]]]] -
BesselY[0, Sqrt[r1^2 + r2^2 - 2*r1*r2*Cos[\[Theta]]]], {\[Theta],
0, 2*Pi}, {r1, 0, 7}, {r2, 0, 7}] // Timing

During evaluation of In[1]:= 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[1]:= 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 52.87383844603967 and 0.0006655311049737613 for the integral and error estimates.

Out[1]= {192.047, 52.8738}


There are 2 messages about working precision and error estimate. Last is very important, it tells us that error is about 0.0006655311049737613. Therefore, we don't need automatic accuracy and precision and we can put these parameters as follows

NIntegrate[
StruveH[0, Sqrt[r1^2 + r2^2 - 2*r1*r2*Cos[\[Theta]]]] -
BesselY[0, Sqrt[r1^2 + r2^2 - 2*r1*r2*Cos[\[Theta]]]], {\[Theta],
0, 2*Pi}, {r1, 0, 7}, {r2, 0, 7}, AccuracyGoal -> 5,
PrecisionGoal -> 4] // Timing

During evaluation of In[2]:= 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.

Out[2]= {41.375, 52.8736}


We compute integral with same error but much faster. But we have message about working precision. Therefore, we can decrease precision and computation time as well

NIntegrate[
StruveH[0, Sqrt[r1^2 + r2^2 - 2*r1*r2*Cos[\[Theta]]]] -
BesselY[0, Sqrt[r1^2 + r2^2 - 2*r1*r2*Cos[\[Theta]]]], {\[Theta],
0, 2*Pi}, {r1, 0, 7}, {r2, 0, 7}, AccuracyGoal -> 2,
PrecisionGoal -> 2] // Timing

Out[3]= {2.09375, 53.0089}


Finally, we decrease computation time in 100 times, but we can improve result as follows

NIntegrate[
StruveH[0, Sqrt[r1^2 + r2^2 - 2*r1*r2*Cos[\[Theta]]]] -
BesselY[0, Sqrt[r1^2 + r2^2 - 2*r1*r2*Cos[\[Theta]]]], {\[Theta],
0, 2*Pi}, {r1, 0, 7}, {r2, 0, 7}, AccuracyGoal -> 3,
PrecisionGoal -> 3] // Timing

Out[5]= {7.67188, 52.8642}



This is maybe the best result for this kind of integrals.

• Thanks, this is a good option. Is there any other way to reduce the calculation time? Commented Sep 25, 2022 at 16:10
• @MamMam For some problem we use compilation - see this post mathematica.stackexchange.com/questions/233509/… Commented Sep 26, 2022 at 4:12
• Yes, I already used compilation. Commented Sep 26, 2022 at 9:03