I have an integral that does not converge satisfactorily for some values and appears random. I have included the code (and plot). However, since the integral depends of imported data, it cannot be tested. Nevertheless, any help through observation alone will be highly appreciated

II[t_, s_] := (288*(-5 + s^2 + 2*t + 
     t^2)^2)/(((1 - s + t)*(1 + s + t))^6)*(0.25*
  N[Pi]*((-5 + s^2 + 2*t + t^2)^2)*
  UnitStep[t - Sqrt[3] + 1] + (-(t - s + 1)*(t + s + 1) + 
   0.5*(-5 + s^2 + t*(2 + t))*
    Log[Abs[(-2 + 2*t + t^2)/(3 - s^2)]])^2);
Omega[k_] := (1./12)*
NIntegrate[((t*(2 + t)*(s^2 - 1))/((1 - s + t)*(1 + s + t)))^2*
   PSpecFunc2[0.5*(t - s + 1)*k]*PSpecFunc2[0.5*(t + s + 1)*k] // 
  Rationalize[#, 0] & // Evaluate, {s, -1, 1}, {t, 0, Infinity}, 
WorkingPrecision -> 40, PrecisionGoal -> 8, MaxRecursion -> 20, 
Method -> {"MonteCarlo", "SymbolicProcessing" -> 0}];
k = Table[i, {i, 10^Range[3, 16, 13./500]}];
GWdata = Table[
Omega[k]*0.68^2, {k, 
10^Range[3, 16, 13./500]}];

In this code the "PSpecFunc2" is an interpolating function created from data that, unfortunately, cannote be imported here. Regardless, when I ran the code, I obtain the following plot enter image description here It is mostly alright except for the regions where the data appears random. Mathematica throws NIntegrate::maxp error stating that the integral failed to converge after 50100 integrand evaluations. I know it is impossible to diagnose this without being able to run it, however any hints will be extremely helpful.

  • $\begingroup$ Your integrand contains singularities. For instance for 1-s+t==0. I guess you need to treat them in the PrincipalValue sense by providing Exclusions. Something along these lines NIntegrate[....,Method -> PrincipalValue, Exclusions -> {1-s+t==0}]. You are using MonteCarlo for speed, however, unless you regularize the integrand somehow (by adding small imaginary part to the denominator), it has no chances to converge. Rather, I would go with above option and reduce PrecisionGoal and AccuracyGoal. $\endgroup$
    – yarchik
    Commented Oct 27, 2019 at 6:27
  • $\begingroup$ Thanks. I will try to take care of the singularity issue. However, should I not use Monte Carlo in this case? $\endgroup$
    – Judas503
    Commented Oct 27, 2019 at 6:36
  • $\begingroup$ Integration in the principal value sense assumes a very specific way how the singularity is approached from both sides. Monte-Carlo is a stochastic method that "knows" nothing about these intricacies. However, if you regularize the function, MC can be very efficient. $\endgroup$
    – yarchik
    Commented Oct 27, 2019 at 6:40
  • $\begingroup$ I tried adding the principal value method. Mathematica states that it can only be used for 1d integrals. $\endgroup$
    – Judas503
    Commented Oct 27, 2019 at 6:44
  • $\begingroup$ I was awaiting this question. Write 2d integral as two 1d integrals. Have a look also here reference.wolfram.com/language/tutorial/… $\endgroup$
    – yarchik
    Commented Oct 27, 2019 at 6:53


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.