I have an integral over a region (rpp):
Msq2[w1_,
w2_] := (12.8228 + 10.518/(0.948338 - 2.0134 w1 - 2.0134 w2) -
6.69841/(1.72935 - 2.0134 w1 - 2.0134 w2) -
57.4434/(2.01348 - 2.0134 w1 - 2.0134 w2) -
13.4997/(3.45415 - 2.0134 w1 - 2.0134 w2) +
9.50782 (1/(-0.110612 + 2.0134 w1) +
1/(-0.110612 + 2.0134 w2)) -
82.5202 (1/(1.14046 + 2.0134 w1) + 1/(1.14046 + 2.0134 w2)))^2;
r1 = Sqrt[0.283 + (Sqrt[-0.018769 + w1^2] + Sqrt[-0.018769 + w2^2])^2];
r2 = Sqrt[0.283 + (Sqrt[-0.018769 + w1^2] - Sqrt[-0.018769 + w2^2])^2];
mphy = 1.007;
rpp = RegionPlot[
Re[r2] < mphy - w1 - w2 < Re[r1], {w1, .11, .4}, {w2, .11, .4},
BoundaryStyle -> Blue, FrameLabel -> {"w1", "w2"},
PlotRangePadding -> 0];
w1min = w2min = Min@rpp[[1, 1, All, 1]];
w1max = w2max = Max@rpp[[1, 1, All, 1]];
NIntegrate[
Boole[Re[r2] < mphy - w1 - w2 < Re[r1]] 1/(64*Pi^3*mphy)*
Msq2[w1, w2], {w1, w1min, w1max}, {w2, w2min, w2max},
AccuracyGoal -> 14] // Chop
With 2 errors mathematica give the following answer (Numerical integration converging too slowly):
Out[8]= 119.375
I tried to correct the problem by evaluating the following expression:
NIntegrate[
Piecewise[{{1/(64*Pi^3*mphy)*Msq2[w1, w2],
Re[r2] <= mphy - w1 - w2 <= Re[r1]}}], {w1, w1min, w1max}, {w2,
w2min, w2max}, MaxRecursion -> 50, AccuracyGoal -> 20,
Method -> {GlobalAdaptive, MaxErrorIncreases -> 10000}]
But the problem remains:
Out[12]= 3095.25
it seems it doesn't converge. I used some other NIntegrate integration strategies but none of them works. With which method does this integal converge and gives a precise and correct answer?
Is there any method which gives a correct answer? What is the real problem? I don't understand what's going on.
@Retay: If I replace Msq2b instead of Msq2, the integral converges. Msq2 and Msq2b are very close to each other ( compare their plots). So what is the problem with Msq2?
Msq2b[w1_,
w2_] := (10.8 + 11.85/(0.959 - 2.013*w1 - 2.013*w2) -
6.9768/(1.8059 - 2.01383*w1 - 2.01383*w2) -
60.9431/(2.04286 - 2.01383*w1 - 2.01383*w2) -
13.21819/(3.4901 - 2.013830*w1 - 2.01383*w2) +
10.461347*(1/(-0.072247 + 2.0138303*w1) +
1/(-0.072247 + 2.013830*w2)) -
82.309287*(1/(1.14002 + 2.013830*w1) +
1/(1.14002 + 2.01383*w2)))^2;
Plot3D[1/(64*Pi^3*mphy)*Msq2[w1, w2], {w1, w1min, w1max}, {w2, w2min,
w2max}, RegionFunction ->
Function[{w1, w2, z}, Re[r2] <= mphy - w1 - w2 <= Re[r1]] ,
BoxRatios -> {1, 1, 4}]
Plot3D[1/(64*Pi^3*mphy)*Msq2b[w1, w2], {w1, w1min, w1max}, {w2, w2min,
w2max}, RegionFunction ->
Function[{w1, w2, z}, Re[r2] <= mphy - w1 - w2 <= Re[r1]] ,
BoxRatios -> {1, 1, 4}]