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}]

=================

I also check the questions in numerical-integration Tag but I didn’t find the answer. I know that the answer must be less than 1.

– Soodeh Z.

Feb 1 ’13 at 7:10

Your efforts to simplify the functions and isolate the problem would go a long way towards making this an answerable question.

– whuber

Feb 1 ’13 at 16:40

I tried to simplify the functions and edit my question. I hope it becomes better to understand. I really want to find the answer as soon as possible, so if there is any problem in my question please tell me. Thanks a lot.

– Soodeh Z.

Feb 1 ’13 at 17:00

1

You’re still asking us to work through and debug several pages of code for you, without providing any preliminary analysis yourself. That approach to problem-solving sometimes works, but it’s unreliable (and could take a long time) because it depends on some kind community member coming along who has the knowledge, time, and interest to do your work for you, as soon as possible. If you could do more to make this question accessible–shorten the code, isolate the problem, simplify the integrand, etc–then it’s much more likely someone will make the effort to provide a good answer.

– whuber

Feb 1 ’13 at 17:26

I eliminate the second part of my question. I didn’t simplify the integrand because by changing the parameters and functions this problem will not occure. I tried every method I knew and it is about 1 week I couldn’t find the answer. The question is that what is the problem with the integrand that none of numeric methods give convergent value.

– Soodeh Z.

Feb 1 ’13 at 17:37

=================

1 Answer

1

=================

Your Msq2 function has a singularity along a line which passes through the integration region:

With[{sing = Solve[1/Msq2[w1, w2] == 0, {w1, w2}] /. Rule -> Equal},

Show[rpp, ContourPlot[sing, {w1, 0, 1}, {w2, 0, 1}]]]

With Msq2b the singularity lies outside the region and so the integral converges:

With[{sing = Solve[1/Msq2b[w1, w2] == 0, {w1, w2}] /. Rule -> Equal},

Show[rpp, ContourPlot[sing, {w1, 0, 1}, {w2, 0, 1}]]]

Thanks a lot for the point which I didn’t see it at all. Now for integrating over this region should I remove the singularity from the bound or there is a better way to do it?

– Soodeh Z.

Feb 1 ’13 at 19:46

@soodeh, I don’t know how best to deal with it, sorry! Have you read this tutorial?

– Simon Woods

Feb 1 ’13 at 20:00

@Woods: As it is obvious there are many points in the region of integration which make Msq2 singular. So is it possible to omit them in the bound or not? How should I solve this problem? Thanks for your help.

– Soodeh Z.

Feb 1 ’13 at 20:03

@Woods: You helped me so much. Thanks a lot. I will read it.

– Soodeh Z.

Feb 1 ’13 at 20:05