MathGroup Archive 2002

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

Search the Archive

Re: Integrate malfunction?

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)]

        0.00146422 y
Out[4]= -------------
            0.5 y

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

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}]

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}]

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
>  _______________

  • Prev by Date: Re: Bug in Sum causes index variable to remain set?
  • Next by Date: Why No Solution Using Solve?
  • Previous by thread: Integrate malfunction?
  • Next by thread: Re: Re: Integrate malfunction?