Re: Re: Integrate malfunction?

One thing is missing in my calculations: when I change a variable in an
integral (x->y/a, a = Sqrt[nitrogen/(k*t)]), I also need to change dx to
dy/a, but I forgot to do this.

These are the correct calculations:

In[1]:= amu = 1.661*10^-27;
nitrogen = 28*amu;
k = 1.381*10^-23;
t = 1000;

In[5]:= a = Sqrt[nitrogen/(k*t)];

In[6]:= expr = 4 Pi ((nitrogen/(2 Pi k t))^(3/2)) x^2 Exp[(-1 nitrogen x^2)/(2 k t)] /. x->y/a;

In[7]:= expr = expr/a  (* because dx = dy/a *)

2
0.797885 y
Out[7]= -----------
2
0.5 y
E

In[8]:= y0 = 11000*a

Out[8]= 20.1864

In[9]:= expr = SetPrecision[expr, 100];

In[10]:= y0 = SetPrecision[y0, 100];

In[11]:= Integrate[expr, {y, y0, Infinity}]

-88
Out[11]= 5.27525415 10

Zbigniew Leyk

zleyk at tca.net wrote:
>  What you get is so called roundoff error because all the calculations are
>  done in machine precision (this is the precision you get when you run a
>  program in C with type double variables). Because your problem is badly
>  scaled (you have small and large numbers in the same expression), therefore
>  your roundoff error is not zero. We can scale the problem introducing a new
>  variable, say y, such that, x = y/Sqrt[nitrogen/(k*t)]
>
>  In[1]:= amu = 1.661*10^-27;
>          nitrogen = 28*amu;
>          k = 1.381*10^-23;
>          t = 1000;
>  In[4]:= expr = 4 Pi ((nitrogen/(2 Pi k t))^(3/2)) x^2 Exp[(-1 nitrogen x^2)/(2 k t)] /. x->y/Sqrt[nitrogen/(k*t)]
>
>                      2
>          0.00146422 y
>  Out[4]= -------------
>                   2
>              0.5 y
>             E
>
>  We need to find the value of y corresponding to 11000:
>
>  In[6]:= y0 = 11000*Sqrt[nitrogen/(k t)]
>
>  Out[6]= 20.1864
>
>  And now we can do the calculation:
>
>  In[7]:= Integrate[expr, {y, y0, Infinity}]
>
>  Out[7]= 0.
>
>  This result is correct when you use machine precision. You cannot get
>  anything better.
>
>  But in Mathematica you can use arbitrary precision numbers. The question is
>  what precision should be used. Let's look at the value of Exp[-0.5 y^2] at
>  y = y0:
>
>  In[9]:= Exp[-0.5 y^2] /. y ->y0
>
>                    -89
>  Out[9]= 3.26725 10
>
>  So it seems that precision 100 should be enough.
>
>  In[10]:= expr = SetPrecision[expr, 100];
>
>  In[11]:= y0 = SetPrecision[y0, 100];
>
>  Then we can do the calculations:
>
>  In[12]:= Integrate[expr, {y, y0, Infinity}]
>
>                        -91
>  Out[12]= 9.68078067 10
>
>  The result is positive (although very small).
>
>  You can do the same computations with NIntegrate:
>
>  In[13]:= NIntegrate[expr, {y, y0, Infinity}]
>
>                     -91
>  Out[13]= 9.68078 10
>
>  The obtained result agrees with the previous result.
>
>  Zbigniew Leyk
>
>
>  Steve Story wrote:
>  >  I recently had to do an integral which should evaluate to be a very very
>  >  small quantity. The function is always positive, but the integral gives me a
>  >  negative number. I think this is a malfunction of some sort. Here's the
>  >  code:
>  >
>  >  amu = 1.661*10^-27
>  >  nitrogen = 28*amu
>  >  k = 1.381*10^-23
>  >  t = 1000
>  >
>  >  Integrate[
>  >    4 Pi ((nitrogen/(2 Pi k t))^(3/2)) x^2 Exp[(-1 nitrogen x^2)/(2 k t)], {x,
>  >  11000, Infinity}]
>  >
>  >  Out[13]= -5.878291776559734*10^-16
>  >
>  >  does anyone know why this is happening? i tried NIntegrate instead, but it
>  >  said:
>  >
>  >  General::"unfl": "Underflow occurred in computation."
>  >
>  >  thanks,
>  >  Steve Story
>  >
>  >
>  >  _______________
>  >
>  >  Steve Story
>  >  Polymer NEXAFS Research Group
>  >  411B Cox
>  >  North Carolina State University
>  >  1-919-515-8147
>  >  sbstory at unity.ncsu.edu
>  >  _______________
>  >

• Prev by Date: Re: Explicit Conjugate: a feature or a bug?
• Next by Date: RE: very simple plot format
• Previous by thread: Re: Integrate malfunction?
• Next by thread: Re: Integrate malfunction?