MathGroup Archive 2002

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

Search the Archive

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?