# Very different results of NIntegrate using different methods: which one to believe

Consider the following function:

SpectrumEHth[th_, Eh_] = Exp[-(Eh/300)]*Cos[th]^20
DiracDeltaArgument[mS_, mh_, ES_, Eh_, ph_, cosa_] =
mh^2 - 2*ES*Eh + 2*Sqrt[ES^2 - mS^2]*ph*cosa;
Cosa[th_, tS_, cosfh_] = cosfh*Sin[th]*Sin[tS] + Cos[th]*Cos[tS];
Print["Solutions of the DiracDeltaArgument = 0 in terms of     \
\!$$\*SubscriptBox[\(E$$, $$h$$]\):"]
DiracDeltaEhSolution[ES_, cosa_, mS_, mh_] =
Eh /. Solve[
DiracDeltaArgument[mS, mh, ES, Eh, Sqrt[Eh^2 - mh^2], cosa] == 0,
Eh]
Print["\!$$\*SubscriptBox[\(E$$, $$h$$]\) must be real - restriction  \
on \!$$\*SubscriptBox[\(E$$, $$S$$]\):"]
DiracDeltaESSolution[cosa_, mS_, mh_] =
ES /. Solve[
Simplify[
DiracDeltaEhSolution[ES, cosa, mS, mh][[1]] -
DiracDeltaEhSolution[ES, cosa, mS, mh][[2]]] == 0, ES][[4]]
Abs[D[DiracDeltaArgument[mS, mh, ES, Eh, Sqrt[Eh^2 - mh^2], cosa],
Eh] /. {Eh -> DiracDeltaEhSolution[ES, cosa, mS, mh][[1]]}];
Abs[D[DiracDeltaArgument[mS, mh, ES, Eh, Sqrt[Eh^2 - mh^2], cosa],
Eh] /. {Eh -> DiracDeltaEhSolution[ES, cosa, mS, mh][[2]]}];
D4PDtSDESdthdfh1[mS_, tS_, th_, fh_, ES_] =
1/(2*Pi) Sin[tS]*
Sqrt[ES^2 -
mS^2]*(SpectrumEHth[th,
DiracDeltaEhSolution[ES, Cosa[th, tS, Cos[fh]], mS, 125][[1]]/
125] + SpectrumEHth[th,
DiracDeltaEhSolution[ES, Cosa[th, tS, Cos[fh]], mS, 125][[2]]/
125]);


I am interested in the integral

dPdtS1[mS_, tS_, method_] :=
NIntegrate[
D4PDtSDESdthdfh1[mS, tS, th, fh, ES], {th, 0, Pi}, {fh, 0,
2*Pi}, {ES, mS,
DiracDeltaESSolution[Cosa[th, tS, Cos[fh]], mS, 125]},
Method -> method]


The integral is defined for 0< mS < 62.5, 0 < ts < Pi. The problem is that using different methods of integration I get completely different results. For example, evaluating dPdtS1[61, Pi/7, method] once, I got

"MonteCarlo", 187 and 146 for the integral and error estimates;

"GlobalAdaptive", 15356 and 2780 for the integral and error estimates;

"QuasiMonteCarlo", 52.21 and 27.89 for the integral and error estimates;

"AdaptiveMonteCarlo", 4559 and 397 for the integral and error estimates.

Moreover, by relaunching the integral I get completely different results. What is a possible reason for such discrepancy, which method to believe and how to get correct result?

• Most likely the integral diverges. Could that be the reason? – AccidentalFourierTransform May 9 '19 at 20:08
• @AccidentalFourierTransform : it seems that the integrand is singular at $\cos(a) = \pm 1$. This is strange, as I expected it to be perfectly finite. – John Taylor May 9 '19 at 20:27
• @AccidentalFourierTransform : after isolating these points from the integral by adding UnitStep[cosa+0.99]*UnitStep[0.99-cosa] to the integrand I obtain similar answers with "GlobalAdaptive" and "MonteCarlo" methods. However, I am wondering whether this procedure is correct. – John Taylor May 9 '19 at 20:34

Somewhat changing your dPdtS1[mS_, tS_, method_]:

dPdtS1[mS_, tS_, method_] := NIntegrate[ D4PDtSDESdthdfh1[mS, tS, th, fh, ES],
{th, 0, Pi}, {fh, 0, 2*Pi}, {ES, mS, DiracDeltaESSolution[Cosa[th, tS, Cos[fh]], mS, 125]},
Method -> method, AccuracyGoal -> 3, AccuracyGoal -> 3, WorkingPrecision -> 12]


, I obtain results which are in accordance:

dPdtS1[61, Pi/7, "MonteCarlo"]


NIntegrate::maxp: The integral failed to converge after 50100 integrand evaluations. NIntegrate obtained 52.552699939412. and 22.945212575512. for the integral and error estimates. 52.5526999394

dPdtS1[61, Pi/7, "QuasiMonteCarlo"]


NIntegrate::maxp: The integral failed to converge after 50000 integrand evaluations. NIntegrate obtained 52.21045201512. and 27.891694754812. for the integral and error estimates.52.2104520150

dPdtS1[61, Pi/7, "AdaptiveMonteCarlo"]


33.9967243466

dPdtS1[61, Pi/7, "GlobalAdaptive"]


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.NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 27 recursive bisections in th near {th,fh,ES} = {0.14285714808198512200893945143233513185940146362578446566502126,8.9357332852024803908579020738908780656678045758304735988042582*10^-9,0.85551545010073033110760179589159365998596592914990664529781283}. NIntegrate obtained 7.061409092945544170324680114102101346280061026736880329252456462.*^6 and 149106.5113046128000673559420138321848074124348877185235094226262. for the integral and error estimates.7.06140909295*10^6

It seems "GlobalAdaptive" is not an appropriate method to this end. You may play with options.

• By relaunching the integral that uses "MonteCarlo" method with your set ups I get completely different results each time. – John Taylor May 9 '19 at 20:21
• @John Taylor: Yes, but each time the result is in accordance with other results, e.g. NIntegrate::maxp: The integral failed to converge after 50100 integrand evaluations. NIntegrate obtained 28.39392172612. and 7.8391810168512. for the integral and error estimates. for my second execution. – user64494 May 9 '19 at 20:44
• Since "MonteCarlo" samples the function randomly, you'll generally get different results from different runs. – John Doty May 9 '19 at 21:30