0
$\begingroup$

I'm having trouble to do the following numerical integrations:

iEAs[Q_?NumericQ, m1_?NumericQ, m2_?NumericQ, mh_?NumericQ, V_?NumericQ, Vc_?NumericQ, A_?NumericQ, Ac_?NumericQ] := m1 Q NIntegrate[ x (A Ac (m2 + m1 (-1 + x)) + V Vc (m1 + m2 -          m1 x))/(8 \[Pi]^2 (-mh^2 (-1 + x) + (m2^2 + m1^2 (-1 + x)) x)), {x, 0, 1}, MaxRecursion -> 100, WorkingPrecision -> 13];

iEAf[Q_?NumericQ, m1_?NumericQ, m2_?NumericQ, mh_?NumericQ, V_?NumericQ, Vc_?NumericQ, A_?NumericQ, Ac_?NumericQ] := m1 Q NIntegrate[-(  (A Ac (12 m2 (1 + x) + m1 (1 - 9 x - 12 x^2)) + V Vc (12 m2 (1 + x) + m1 (-1 + 9 x + 12 x^2))))/(96 \[Pi]^2 (-m2^2 (-1 +x) + (mh^2 + m1^2 (-1 + x)) x)), {x, 0, 1}, MaxRecursion -> 100, WorkingPrecision -> 13];

iVAv[G_?NumericQ, m1_?NumericQ, m2_?NumericQ, M_?NumericQ, V_?NumericQ, Vc_?NumericQ, A_?NumericQ, Ac_?NumericQ] := G m1 NIntegrate[- (A Ac (-3 m2 (1 + x) + m1 (-3 + x + 2 x^2)) +      V Vc (3 m2 (1 + x) + m1 (-3 + x + 2 x^2)))/(8 \[Pi]^2 (M^2 (-1 + x) - (m2^2 + m1^2 (-1 + x)) x)), {x, 0, 1}, MaxRecursion -> 100, WorkingPrecision -> 13];

iVAf[Q_?NumericQ, m1_?NumericQ, m2_?NumericQ, M_?NumericQ, V_?NumericQ, Vc_?NumericQ, A_?NumericQ, Ac_?NumericQ] := m1 Q NIntegrate[ -( (V Vc (m1 - 2 m2 + m1 x) + A Ac (m1 + 2 m2 + 
m1 x))/(4 \[Pi]^2 (M^2 + (m1^2 + m2^2) (-1 + x)))), {x, 0, 1}, MaxRecursion -> 100, WorkingPrecision -> 13];

Then when I try to calculate my other function, that is just a sum of several of those integrals taken with diferent arguments, I get the following errror:

During evaluation of In[1]:= NIntegrate::zeroregion: Integration region {{1.000000000000,0.999999999999999999999999999972752437499605743914855469900697130}} cannot be further subdivided at the specified working precision. NIntegrate assumes zero integral there and on any further indivisible regions. >>

During evaluation of In[1]:= NIntegrate::inumri: The integrand -((3.94734084367020280294200477311377373096717807463275517933069283*10^-10 (1.2679004580000000000 (1+x)+0.10565837150000000000 (1-9 x-12 x^2)))/(-0.0111636914680320122499999999999999999999999999999999999999121163 (-1+x)+0.011163691468032012250 (-1+x) x)) has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{0.999999999999999999999999999972752437499605743914855469900697130,0.999999999999999999999990130457402984116116415865423076016809530}}.

>

During evaluation of In[1]:= NIntegrate::zeroregion: Integration region {{1.000000000000,0.999999999999999999999999999972752437499605743914855469900697130}} cannot be further subdivided at the specified working precision. NIntegrate assumes zero integral there and on any further indivisible regions. >>

During evaluation of In[1]:= NIntegrate::inumri: The integrand -((3.94734084367020280294200477311377373096717807463275517933069283*10^-10 (1.2679004580000000000 (1+x)+0.10565837150000000000 (1-9 x-12 x^2)))/(-0.0111636914680320122499999999999999999999999999999999999999121163 (-1+x)+0.011163691468032012250 (-1+x) x)) has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{0.999999999999999999999999999972752437499605743914855469900697130,0.999999999999999999999990130457402984116116415865423076016809530}}.

>

During evaluation of In[1]:= 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. >>

During evaluation of In[1]:= NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 400 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained -1.102332818474*10^-28 and 1.773650837526`13.*^-29 for the integral and error estimates. >>

During evaluation of In[1]:= 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. >>

During evaluation of In[1]:= NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 400 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 1.6842750826060307585060318968899777955348273049856086641756531763.*^-25 and 6.2565249622956313692877327314984953630958701988698617625287317863.*^-26 for the integral and error estimates. >>

During evaluation of In[1]:= 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. >>

During evaluation of In[1]:= General::stop: Further output of NIntegrate::slwcon will be suppressed during this calculation. >>

During evaluation of In[1]:= NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 400 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained -1.343930509261*10^-26 and 7.504336838898`13.*^-27 for the integral and error estimates. >>

During evaluation of In[1]:= General::stop: Further output of NIntegrate::eincr will be suppressed during this calculation. >>

During evaluation of In[1]:= NIntegrate::zeroregion: Integration region {{1.000000000000,0.999999999999999999999999999972752437499605743914855469900697130}} cannot be further subdivided at the specified working precision. NIntegrate assumes zero integral there and on any further indivisible regions. >>

During evaluation of In[1]:= General::stop: Further output of NIntegrate::zeroregion will be suppressed during this calculation. >>

During evaluation of In[1]:= NIntegrate::inumri: The integrand -((3.94734084367020280294200477311377373096717807463275517933069283*10^-10 (1.2679004580000000000 (1+x)+0.10565837150000000000 (1-9 x-12 x^2)))/(-0.0111636914680320122499999999999999999999999999999999999999121163 (-1+x)+0.011163691468032012250 (-1+x) x)) has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{0.999999999999999999999999999972752437499605743914855469900697130,0.999999999999999999999990130457402984116116415865423076016809530}}.

>

During evaluation of In[1]:= General::stop: Further output of NIntegrate::inumri will be suppressed during this calculation. >>

Out[56]= -0.1054127670153 - 0.10565837150000000000 NIntegrate[-(0.0006115582135884999042 \ 0.0006115582135884999042 (12 0.10565837150000000000 (1 + x) + 0.10565837150000000000 (1 - 9 x - 12 x^2)) + 0 0 (12 0.10565837150000000000 (1 + x) + 0.10565837150000000000 (-1 + 9 x + 12 x^2)))/(96 [Pi]^2 (-0.10565837150000000000^2 (-1 + x) + (0^2 + 0.10565837150000000000^2 (-1 + x)) x)), {x, 0, 1}, MaxRecursion -> 100, WorkingPrecision -> 13]

I have tried many methods of integration. For most of them I get the kind of error above, for local adaptive methods I don't get any error, but I left Mathematica running for a whole afternoon and it didn't finish the calculation (it just kept running).

I need the WorkingPrecision->13 because I'm going to compare this result with the electron anomalous magnetic moment (which has got 12 significant digits for it's experimental value). My final goal is to use thse integrals for several values of it's parameters, and then make RegionPlots comparing my calculations with the experimental value.

EDIT: If I get divergences it is fine for me to discard the parameter values (m2, mh and M) where they happen.

So, if anyone can find a way to solve my problem, I would be very grateful.

Improve this question
$\endgroup$
5
  • $\begingroup$ First error, Pi^2,not [Pi]^2. $\endgroup$
    – Apple
    Jul 2, 2014 at 14:27
  • $\begingroup$ @Chenminqi Or, perhaps \[Pi]^2, as it now appears in the code block. $\endgroup$
    – Mark McClure
    Jul 2, 2014 at 14:27
  • $\begingroup$ the value of Q, m1, m2, mh, V, Vc, A, Ac ? $\endgroup$
    – Apple
    Jul 2, 2014 at 14:38
  • $\begingroup$ Can you edit your question to show whether or not the first function, iEAs is behaving properly? $\endgroup$
    – bobthechemist
    Jul 2, 2014 at 15:10
  • $\begingroup$ Q is 1 or -1. m1 is either 10^-1 or 10^-4. m2 and mh I will vary from 0 to 2000. V, Vc, A, Ac and G are of order 10^-3. Also, when I ploted the arguments of those integrals, everything went well. But I know that for some values of m1, m2 and mh I may get divergences. $\endgroup$
    – George
    Jul 2, 2014 at 16:14

2 Answers 2

Reset to default
9
$\begingroup$

Just consider the first integral.

expr = x (A Ac (m2 + m1 (-1 + x)) + V Vc (m1 + m2 - 
     m1 x))/(8 π^2 (-mh^2 (-1 + x) + (m2^2 + m1^2 (-1 + x)) x));
denominator = Collect[Denominator[expr], x]

8 mh^2 π^2 + 8 (-m1^2 + m2^2 - mh^2) π^2 x + 8 m1^2 π^2 x^2

It has two singular points.

sol = Solve[denominator == 0, x] // Simplify

enter image description here

If the singular points are between 0 and 1:

Reduce[{0 < # < 1, mh > 0, m2 > 0, m1 > 0}] & /@ (sol[[All, 1, 2]])

{m2 > 0 && mh > 0 && m1 >= m2 + mh, m2 > 0 && mh > 0 && m1 >= m2 + mh}

So if m1 >= m2 + mh, the singular points are between 0 and 1. For example:

expr2=expr /. {A -> 1, V -> 10, Vc -> 10, Ac -> 10, m1 -> 10^-1, m2 -> 4.5*10^-3, 
    mh -> 4.5*10^-3};
Plot[expr2, {x, 0, 1}]

enter image description here

then Integrate[expr2, {x, 0, 1}] will diverge, just like Integrate[1/x, {x, -1, 2}]. So we need PrincipalValue -> True.

Plot[1/x, {x, -1, 2}]

enter image description here

Integrate[1/x, {x, -1, 2}, PrincipalValue -> True]

Log[2]

N@Integrate[expr2, {x, 0, 1}, PrincipalValue -> True]

-34.7311

If m1 < m2 + mh, there are no singular points between 0 and 1. For example:

expr3 = expr /. {A -> 1, V -> 10, Vc -> 10, Ac -> 10, m1 -> 10^-1, m2 -> 1/10,
     mh -> 1/15};
Plot[expr3, {x, 0, 1}]

enter image description here

Then just integrate it.

NIntegrate[expr3, {x, 0, 1}]
N@Integrate[expr3, {x, 0, 1}]

15.5586

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks Chenminqi, your solution worked. Also, thanks to everyone else that tried to solve it. I was stuck on this problem for quite a while. $\endgroup$
    – George
    Jul 2, 2014 at 17:13
  • $\begingroup$ Nice answer +1. I fixed quite a few typos for you. I hope I made no mistakes. $\endgroup$
    – Michael E2
    Jul 2, 2014 at 17:20
2
$\begingroup$

All of your integrals can be done with Integrate. These initial calculations are slow but by using Set rather than SetDelayed their subsequent use will be much quicker. You will also get better precision.

Since you are comparing with experimental results presumably you satisfy the conditions suppressed by GenerateConditions -> False. If you want to see the ConditionalExpression just eliminate GenerateConditions -> False. 

iEAs[Q_, m1_, m2_, mh_, V_, Vc_, A_, Ac_] =
  m1 Q Integrate[
    x (A Ac (m2 + m1 (-1 + x)) + V Vc (m1 + m2 - m1 x))/
      (8 \[Pi]^2 (-mh^2 (-1 + x) + (m2^2 + m1^2 (-1 + x)) x)),
    {x, 0, 1}, GenerateConditions -> False];

iEAf[Q_, m1_, m2_, mh_, V_, Vc_, A_, Ac_] =
  m1 Q Integrate[
    -((A Ac (12 m2 (1 + x) + m1 (1 - 9 x - 12 x^2)) +
         V Vc (12 m2 (1 + x) + m1 (-1 + 9 x + 12 x^2))))/
     (96 \[Pi]^2 (-m2^2 (-1 + x) + (mh^2 + m1^2 (-1 + x)) x)),
    {x, 0, 1}, GenerateConditions -> False];

iVAv[G_, m1_, m2_, M_, V_, Vc_, A_, Ac_] =
  G m1 Integrate[
    -(A Ac (-3 m2 (1 + x) + m1 (-3 + x + 2 x^2)) +
        V Vc (3 m2 (1 + x) + m1 (-3 + x + 2 x^2)))/
     (8 \[Pi]^2 (M^2 (-1 + x) - (m2^2 + m1^2 (-1 + x)) x)),
    {x, 0, 1}, GenerateConditions -> False];

iVAf[Q_, m1_, m2_, M_, V_, Vc_, A_, Ac_] =
  m1 Q Integrate[
    -((V Vc (m1 - 2 m2 + m1 x) + A Ac (m1 + 2 m2 + m1 x))/
       (4 \[Pi]^2 (M^2 + (m1^2 + m2^2) (-1 + x)))),
    {x, 0, 1}, GenerateConditions -> False];
Improve this answer
$\endgroup$
3
  • $\begingroup$ I've tried your idea, using Integrate instead of NIntegrate (with and without the GenerateConditions->False).In both cases I get convergence problems and a complex result. I was using NIntegrate to make sure that I would get a real result. $\endgroup$
    – George
    Jul 2, 2014 at 16:10
  • $\begingroup$ @George - what are the values of the parameters? $\endgroup$
    – Bob Hanlon
    Jul 2, 2014 at 16:21
  • $\begingroup$ Q is 1 or -1. m1 is either 10^-1 or 10^-4. M, m2 and mh I will vary from 0 to 2000. V, Vc, A, Ac and G are of order 10^-3. $\endgroup$
    – George
    Jul 2, 2014 at 16:22

Your Answer

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

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