1
$\begingroup$

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]]
DiracDeltaDerivative1[ES_, cosa_, mS_, mh_] = 
  Abs[D[DiracDeltaArgument[mS, mh, ES, Eh, Sqrt[Eh^2 - mh^2], cosa], 
     Eh] /. {Eh -> DiracDeltaEhSolution[ES, cosa, mS, mh][[1]]}];
DiracDeltaDerivative2[ES_, cosa_, mS_, mh_] = 
  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]/DiracDeltaDerivative1[ES, Cosa[th, tS, Cos[fh]], mS, 
       125] + SpectrumEHth[th, 
       DiracDeltaEhSolution[ES, Cosa[th, tS, Cos[fh]], mS, 125][[2]]/
        125]/DiracDeltaDerivative2[ES, Cosa[th, tS, Cos[fh]], mS, 
       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?

$\endgroup$
  • 2
    $\begingroup$ Most likely the integral diverges. Could that be the reason? $\endgroup$ – AccidentalFourierTransform May 9 '19 at 20:08
  • $\begingroup$ @AccidentalFourierTransform : it seems that the integrand is singular at $\cos(a) = \pm 1$. This is strange, as I expected it to be perfectly finite. $\endgroup$ – John Taylor May 9 '19 at 20:27
  • $\begingroup$ @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. $\endgroup$ – John Taylor May 9 '19 at 20:34
-1
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ By relaunching the integral that uses "MonteCarlo" method with your set ups I get completely different results each time. $\endgroup$ – John Taylor May 9 '19 at 20:21
  • $\begingroup$ @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. $\endgroup$ – user64494 May 9 '19 at 20:44
  • $\begingroup$ Since "MonteCarlo" samples the function randomly, you'll generally get different results from different runs. $\endgroup$ – John Doty May 9 '19 at 21:30

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.