Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Yesterday I asked a question about the non converging integral. Woods told me that it is due to the function which has a singularity along a line which passes through the integration region. (Why does this numeric integral fail to converge?) Now the question is that how can I find a correct answer for this integral? Should I eliminate these singularities from the integral bound? How can I understand it is integrable singularity and gives convergent value or not? I read singularity handling in documentation center, but I couldn't solve the problem using Exclusions.(http://reference.wolfram.com/mathematica/tutorial/NIntegrateIntegrationStrategies.html#122144792)

function:

 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;

The region of the integral:

rpp = RegionPlot[
   Re[r2] < mphy - w1 - w2 < Re[r1], {w1, .11, .4}, {w2, .11, .4},  
   BoundaryStyle -> Blue, FrameLabel -> {"w1", "w2"}, 
   PlotRangePadding -> 0];

Where:

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;

Integration:

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

Where:

 w1min = w2min = Min@rpp[[1, 1, All, 1]];
 w1max = w2max = Max@rpp[[1, 1, All, 1]];

The green line is the singular curve which a part of it is in 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}]]]

enter image description here

With the following code I can find the singular line, but I don't know how to eliminate this line from integration bound:

sing1 = Solve[
    1/Msq2[w1, w2] == 0 && w1min < w1 < w1max && w2min < w2 < w2max &&
      Re[r2] < mphy - w1 - w2 < Re[r1], {w1, w2}][[1, 1]] /. 
  Rule -> Equal

ConditionalExpression[ w1 == 9.93345*10^-7 (474169.[VeryThinSpace]- 1.0067*10^6 w2), 0.209029 < w2 < 0.261985]

share|improve this question
add comment

1 Answer

Not an answer. I tried to use Exclusions but it doesn't seem to be reliable :

sing1 // Simplify
(* ConditionalExpression[w1 == 0.471013 - 1. w2, 0.209029 < w2 < 0.261985] *)

NIntegrate[Boole[Re[r2] < mphy - w1 - w2 < Re[r1]] 1/(64*Pi^3*mphy)*Msq2[w1, w2], 
  {w1, w1min, w1max}, {w2, w2min, w2max},  
  Exclusions -> {{1/Msq2[w1, w2] == 0, w1min < w1 < w1max && w2min < w2 < w2max && 
                  Re[r2] < mphy - w1 - w2 < Re[r1]}}]

(* 21.8701 *)

NIntegrate[Boole[Re[r2] < mphy - w1 - w2 < Re[r1]] 1/(64*Pi^3*mphy)*Msq2[w1, w2], 
  {w1, w1min, w1max}, {w2, w2min, w2max}, 
  Exclusions -> {{w1 == 0.4710132114830634` - 0.9999999999999999` w2, 
                  0.20902863354401882` < w2 < 0.2619845779390446`}}]

(* 4.32175*10^13 *)
share|improve this answer
1  
@gatessucks: You are true. It is not reliable and NIntegrate still has errors about some singular points although it seems they are excluded. It is strange. May be there is sth. wrong or we don't have integrable singularity? I don't know how to fix it. –  Soodeh Z. Feb 2 '13 at 9:42
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.