# How to improve the calculation accuracy and speed of double numerical integration？

I want to solve the following double numerical integral:

$$A=\int _{0.52}^1\int _0^{\infty }\frac{\sqrt{s^2-0.52^2} \sqrt{1-s^2} \left(2 s^2-1\right) \left(\left(2 s^2-1\right) \text{T1}-2 s^2 \text{T2}\right)}{\left(2 s^2-1\right)^4+16 \left(s^2-0.52^2\right) s^4 \left(1-s^2\right)}d\xi ds$$

$$\text{T1}=\frac{\cos (-8.2 \xi ) \sin \left(\frac{\xi (0.936 t)}{s}\right)}{\xi }-\frac{\cos (11.8 \xi ) \sin \left(\frac{\xi (0.936 t)}{s}\right)}{\xi }$$

$$\text{T2}=\frac{\cos (-8.2 \xi ) \sin (0.877 \xi t)}{\xi }-\frac{\cos (11.8 \xi ) \sin (0.877 \xi t)}{\xi }$$

I just used the "Table" function to output the t vector：

Table[NIntegrate[A,{ξ ,0,infinity },{s,0.52,1}],{t,0,25,0.25}]


Adopting the "GlobalAdaptive" integration strategy takes a long time, and the error is large.And I get the following warning

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.939973715732701*^-6 and 0.0010362662139119024 for the integral and error estimates.


So I want to know how to improve the integration formula, or what integration strategy to use to improve the calculation speed and accuracy. Thanks!!!

• It is expected that you write out the Mathematica code for your integral such that it is easily copied and pasted by others. Can you, please, translate your latex formulas to Mathematica code (or, include the code that you have assumably already written!) so that other users can potentially answer this question without the burden of translating your latex formulas for you? Thanks! Aug 18, 2021 at 13:59

You can do the integration over xi with Integrate .

T1 = -((Sin[(0.936 t)/s \[Xi]] Cos[
11.8 \[Xi]])/\[Xi]) + (Sin[(0.936 t)/
s \[Xi]] Cos[-8.2 \[Xi]])/\[Xi];
T2 = -((Sin[0.877 t \[Xi]] Cos[11.8 \[Xi]])/\[Xi]) + (Sin[
0.877 t \[Xi]] Cos[-8.2 \[Xi]])/\[Xi];

A = ((Sqrt[s^2 - 0.52^2] Sqrt[
1 - s^2] (2 s^2 - 1) ((2 s^2 - 1) T1 -
2 s^2 T2))/((2 s^2 - 1)^4 + 16 s^4 (s^2 - 0.52^2) (1 - s^2)));
AR[ \[Xi]_, s_, t_] = A // Rationalize[#, 0] & // Together

int[s_, t_] =
Integrate[AR[\[Xi], s, t], {\[Xi], 0, \[Infinity]},
Assumptions -> {0 < t < 25, 52/100 < s < 1}]

(*   (25 \[Pi] (-1 + 2 s^2) Sqrt[-169 + 794 s^2 -
625 s^4] ((-1 + 2 s^2) Sign[1025 s - 117 t] + (1 - 2 s^2) Sign[
1475 s - 117 t] +
2 s^2 (-Sign[41/5 - (877 t)/1000] +
Sign[59/5 - (877 t)/1000])))/(4 (-625 + 5000 s^2 -
12296 s^4 + 7296 s^6))   *)

(tab = Table[{t, NIntegrate[int[s, t], {s, 52/100, 1}]}, {t, 0,
25, .25}]) // Timing

ListPlot[tab]


I testet it for t=625/100. (Cases with 3 equal values should be tested further.)

NIntegrate[int[s, 625/100], {s, 52/100, 1}]
(*   0.0455044   *)

NIntegrate[
AR[\[Xi], s, 625/100], {s, 52/100, 1}, {\[Xi], 0, \[Infinity]}]
(*   0.0454959   *)

• Thank you for your enthusiastic answer！! It’s really efficient to divide this integral into two parts. Aug 19, 2021 at 3:04
A=((Sqrt[s^2 - 0.52^2] Sqrt[1 - s^2] (2 s^2 - 1) ((2 s^2 - 1) T1 - 2 s^2 T2))/((2 s^2 - 1)^4 + 16 s^4 (s^2 - 0.52^2) (1 - s^2)))
T1=-((Sin[ (0.936 t)/s ξ] Cos[11.8 ξ ])/ξ) + (Sin[ (0.936 t)/s ξ] Cos[ -8.2 ξ])/ξ
T2=-((Sin[0.877 t ξ ] Cos[11.8 ξ ])/ξ) + (Sin[0.877 t ξ ] Cos[ -8.2 ξ])/ξ

• Please edit the original post to put this code in it (and thanks for adding the code). Aug 18, 2021 at 14:59