Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2005
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2005

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Num. integration problem in Mathematica

  • To: mathgroup at smc.vnet.net
  • Subject: [mg53826] Re: Num. integration problem in Mathematica
  • From: Urijah Kaplan <uak at sas.upenn.edu>
  • Date: Fri, 28 Jan 2005 02:44:37 -0500 (EST)
  • Organization: University of Pennsylvania
  • References: <ctahqm$apt$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Using Mtm 5.1 and making some small changes (don't use floating point 
values!!!) I received a (hopefully) correct answer.

Ma = 10^10;
m1 = Rationalize[0.06*10^(-9)];
mPl = Rationalize[1.221*10^19];

v = 174;
lambda = Rationalize[0.4];
yt = Rationalize[0.6];
gy = ((1/81)*4*Pi)^(1/2);
g2 = ((1/38)*4*Pi)^(1/2);
g3 = ((1/26)*4*Pi)^(1/2);

mH = ((3/16)*g2^2 + (1/16)*gy^2 + (1/4)*yt^2 + (1/2)*lambda)^(1/2)*Ma;
mL = ((3/32)*g2^2 + (1/32)*gy^2)^(1/2)*Ma;
mQ = ((1/6)*g3^2 + (3/32)*g2^2 + (1/288)*gy^2 + (1/16)*yt^2)^(1/2)*Ma;
mU = ((1/6)*g3^2 + (1/18)*gy^2 + (1/8)*yt^2)^(1/2)*Ma;

aH = mH^2/Ma^2;
aL = mL^2/Ma^2;
aQ = mQ^2/Ma^2;

x = s/Ma^2;

smin = Max[(mQ + mU)^2, (mL + Ma)^2];
s = smin/y;

NIntegrate[
     s^(1/2)*BesselK[1, s^(1/2)/Ma]*(3/(4*Pi))*(Ma*m1/v^2)*yt^2*((x - 1 -
       aL)*(x - 2*aQ)/(x*(
     x - aH)^2))*(((1 + aL - x)^2 - 4*aL)*(1 - 4*aQ/x))^(1/
             2)*(1/y^2), {y, 0, 1}, MinRecursion -> 10, MaxRecursion ->
         50, WorkingPrecision -> 35]



23573.56641066495124896399363987458797535501`25.202227120473328


(`25.202227120473328 means that the answer is given to about 25 digits of 
precision)

--Urijah Kaplan


Antonio Cardoso wrote:
> Hello, I'm trying to solve this numerical integration in Mathematica:
> 
> Ma = 10^10;
> m1 = 0.06*10^(-9);
> mPl = 1.221*10^19;
> 
> v = 174;
> lambda = 0.4;
> yt = 0.6;
> gy = ((1/81)*4*Pi)^0.5;
> g2 = ((1/38)*4*Pi)^0.5;
> g3 = ((1/26)*4*Pi)^0.5;
> 
> mH = ((3/16)*g2^2 + (1/16)*gy^2 + (1/4)*yt^2 + (1/2)*lambda)^0.5*Ma;
> mL = ((3/32)*g2^2 + (1/32)*gy^2)^0.5*Ma;
> mQ = ((1/6)*g3^2 + (3/32)*g2^2 + (1/288)*gy^2 + (1/16)*yt^2)^0.5*Ma;
> mU = ((1/6)*g3^2 + (1/18)*gy^2 + (1/8)*yt^2)^0.5*Ma;
> 
> aH = mH^2/Ma^2;
> aL = mL^2/Ma^2;
> aQ = mQ^2/Ma^2;
> 
> x = s/Ma^2;
> 
> smin = Max[(mQ + mU)^2, (mL + Ma)^2];
> s = smin/y;
> 
> NIntegrate[s^0.5*BesselK[1,s^0.5/Ma]*(3/(4*Pi))*(Ma*m1/v^2)*yt^2*((x-1-aL)*(x-2*aQ)/(x*(x-aH)^2))*(((1+aL-x)^2-4*aL)*(1-4*aQ/x))^0.5*(1/y^2),{y,0,1},MinRecursion->10,MaxRecursion->35,WorkingPrecision->35]
> 
> but it can't give me the right answer. I't gives the error:
> 
> NIntegrate::slwcon:
> Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration being 0, oscillatory integrand, or insufficient WorkingPrecision. If your integrand is oscillatory try using the option Method->Oscillatory in NIntegrate.
> 
> I have already tried to change the options, like WorkingPrecision, PrecisionGoal, SingularityDepth, etc., but it didn't work.
> 
> Can someone help me to solve this problem?
> 
> Thank you,
> 
> Antonio Cardoso


  • Prev by Date: Re: Algebra of Einstein velocity addition
  • Next by Date: Re: Solving a weakly singular integral equation
  • Previous by thread: Re: Num. integration problem in Mathematica
  • Next by thread: Algebra of Einstein velocity addition