4
$\begingroup$

I've been trying to calculate an integral using NIntegrate as follows.

Q2 = 189^2;
f[x1_, x2_, Mph_] := -((1.2539284785732083*^-8*Mph^4*x1*x2)/
         ((x1 - 1.)^2*(x2 - 1.)^2)) + 
    (1.2539284785732083*^-8*Mph^4*x2)/((x1 - 1.)^2*(x2 - 1.)^2) + 
    (1.2539284785732083*^-8*Mph^4*x1)/((x1 - 1.)^2*(x2 - 1.)^2) - 
    (1.2539284785732083*^-8*Mph^4)/((x1 - 1.)^2*(x2 - 1.)^2) - 
    (0.0004479157918311357*Mph^2*x1^2*x2)/
      ((x1 - 1.)^2*(x2 - 1.)^2) + (0.0006718736877467036*Mph^2*x1^2)/
      ((x1 - 1.)^2*(x2 - 1.)^2) - 
    (0.0004479157918311357*Mph^2*x1*x2^2)/
      ((x1 - 1.)^2*(x2 - 1.)^2) + (0.0006718736877467036*Mph^2*x2^2)/
      ((x1 - 1.)^2*(x2 - 1.)^2) + 
    (0.0008958315836622714*Mph^2*x1*x2)/((x1 - 1.)^2*(x2 - 1.)^2) - 
    (0.0008958315836622714*Mph^2*x2)/((x1 - 1.)^2*(x2 - 1.)^2) - 
    (0.0008958315836622714*Mph^2*x1)/((x1 - 1.)^2*(x2 - 1.)^2) + 
    (0.0004479157918311357*Mph^2)/((x1 - 1.)^2*(x2 - 1.)^2) - 
    (8.*x1^3*x2)/((x1 - 1.)^2*(x2 - 1.)^2) + 
    (8.*x1^3)/((x1 - 1.)^2*(x2 - 1.)^2) + 
    (8.*x1^2*x2)/((x1 - 1.)^2*(x2 - 1.)^2) - 
    (8.*x1^2)/((x1 - 1.)^2*(x2 - 1.)^2) - 
    (8.*x1*x2^3)/((x1 - 1.)^2*(x2 - 1.)^2) + 
    (8.*x2^3)/((x1 - 1.)^2*(x2 - 1.)^2) + 
    (8.*x1*x2^2)/((x1 - 1.)^2*(x2 - 1.)^2) - 
    (8.*x2^2)/((x1 - 1.)^2*(x2 - 1.)^2); 
f1[m_] := 
 NIntegrate[
  f[x1, x2, m], {x1, 0, 1 - m^2/Q2}, {x2, 1 - x1 - m^2/Q2, 
   1 - 1/(1 - x1)*m^2/Q2}]
Plot[f1[10^x], {x, -5, 2}]

However, Mathematica keeps giving me errors like this:

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 18 recursive bisections in x2 near {x1,x2} = {0.429446,3.22008*10^-6}. NIntegrate obtained -5750.51 and 297.825 for the integral and error estimates.

Also, the plot looks fluctuating in the region of x below -2. enter image description here

I tried to increase WorkingPrecison but to no avail. Could anyone tell me how to fix this?

Improve this question
$\endgroup$

2 Answers 2

Reset to default
8
$\begingroup$

A completely analytical approach also works:

f[x1_, x2_, Mph_] := 
  1/((-1 + x1)^2 (-1 + x2)^2) (Mph^4 p (-1 + x1 + x2 - x1 x2) - 
     8 (-1 + x1) (-1 + x2) (x1^2 + x2^2) + 
     Mph^2 (q - s x1 + s (-1 + x1) x2 - q x1 x2 (x1 + x2) + 
        r (x1^2 + x2^2)));
y1 = Integrate[
    f[x1, x2, m], {x2, 1 - x1 - m^2/Q2, 1 - 1/(1 - x1)*m^2/Q2}, 
    Assumptions -> 
     Q2 > 0 && p > 0 && q > 0 && r > 0 && s > 0 && m > 0 && 
      1 > x1 > 0] // First;
y2 = Integrate[y1, x1];
g[m_] = (y2 /. {x1 -> 1 - m^2/Q2}) - (y2 /. {x1 -> 0});
Q2 = 189^2;
p = 1.2539284785732083*^-8;
q = 0.0004479157918311357;
r = 0.0006718736877467036;
s = 0.0008958315836622714;
Plot[g[10^x], {x, -5, 2}]
Improve this answer
$\endgroup$
3
  • $\begingroup$ After posting my answer, I saw that @ydd proposed essentially the same approach. $\endgroup$ Commented Oct 9 at 14:35
  • 2
    $\begingroup$ Your solution with exact symbolic parameters looks very nice (+1). $\endgroup$ Commented Oct 10 at 0:14
  • 1
    $\begingroup$ It works like a charm. Thank you very much! $\endgroup$ Commented Oct 10 at 9:23
8
$\begingroup$

We can simplify this problem by integrating on x2 as follows

Q2 = 189^2;
f[x1_, x2_, 
   Mph_] := -((1.2539284785732083*^-8*Mph^4*x1*
        x2)/((x1 - 1.)^2*(x2 - 1.)^2)) + (1.2539284785732083*^-8*
      Mph^4*x2)/((x1 - 1.)^2*(x2 - 1.)^2) + (1.2539284785732083*^-8*
      Mph^4*x1)/((x1 - 1.)^2*(x2 - 1.)^2) - (1.2539284785732083*^-8*
      Mph^4)/((x1 - 1.)^2*(x2 - 1.)^2) - (0.0004479157918311357*Mph^2*
      x1^2*x2)/((x1 - 1.)^2*(x2 - 1.)^2) + (0.0006718736877467036*
      Mph^2*x1^2)/((x1 - 1.)^2*(x2 - 1.)^2) - (0.0004479157918311357*
      Mph^2*x1*
      x2^2)/((x1 - 1.)^2*(x2 - 1.)^2) + (0.0006718736877467036*Mph^2*
      x2^2)/((x1 - 1.)^2*(x2 - 1.)^2) + (0.0008958315836622714*Mph^2*
      x1*x2)/((x1 - 1.)^2*(x2 - 1.)^2) - (0.0008958315836622714*Mph^2*
      x2)/((x1 - 1.)^2*(x2 - 1.)^2) - (0.0008958315836622714*Mph^2*
      x1)/((x1 - 1.)^2*(x2 - 1.)^2) + (0.0004479157918311357*
      Mph^2)/((x1 - 1.)^2*(x2 - 1.)^2) - (8.*x1^3*
      x2)/((x1 - 1.)^2*(x2 - 1.)^2) + (8.*
      x1^3)/((x1 - 1.)^2*(x2 - 1.)^2) + (8.*x1^2*
      x2)/((x1 - 1.)^2*(x2 - 1.)^2) - (8.*
      x1^2)/((x1 - 1.)^2*(x2 - 1.)^2) - (8.*x1*
      x2^3)/((x1 - 1.)^2*(x2 - 1.)^2) + (8.*
      x2^3)/((x1 - 1.)^2*(x2 - 1.)^2) + (8.*x1*
      x2^2)/((x1 - 1.)^2*(x2 - 1.)^2) - (8.*
      x2^2)/((x1 - 1.)^2*(x2 - 1.)^2);
F = Integrate[f[x1, x2, m], x2]

(*Out[]= 1/(-1. + x1)^2 ((
   m^2 (-0.000223958 + 0.000447916 x1 - 0.000223958 x1^2))/(-1. + 
    x2) - 0.000447916 (-17860.5 + 17860.5 x1 + 
      m^2 (-1.5 + 1. x1)) x2 - 4. (-1. + x1) x2^2 - 
   1.25393*10^-8 (-6.37995*10^8 + 6.37995*10^8 x1 - 
      6.37995*10^8 x1^2 + 6.37995*10^8 x1^3 + m^4 (-1. + 1. x1) + 
      m^2 (-35721. + 35721. x1^2)) Log[1. - 1. x2])*)

Then we define

F2 = F /. x2 -> 1 - 1/(1 - x1)*m^2/Q2 // Simplify;

F1 = F /. x2 -> 1 - x1 - m^2/Q2 // Simplify;
Int[m0_?NumericQ] := 
 NIntegrate[F2 - F1 /. m -> m0, {x1, 0, 1 - m0^2/Q2}];

Visualization

Plot[Int[10^x] // Quiet, {x, -5, 2}]

Figure1

Improve this answer
$\endgroup$
2
  • 2
    $\begingroup$ Nice answer. For what it's worth we can evaluate the definite integrals directly with some Assumptions (although it is pretty slow to do). $\endgroup$ Commented Oct 9 at 14:29
  • 1
    $\begingroup$ @ydd Yes, we can define integral directly. The problem is that the integral has complex part which is not small at $m=10^{−5}$. Therefore, we have same oscillations as for f1[m]. But it seems that your solution with rational parameters has no imaginary part. $\endgroup$ Commented Oct 10 at 0:11

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.