Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Integrate malfunction?

• To: mathgroup at smc.vnet.net
• Subject: [mg33886] Re: [mg33860] Integrate malfunction?
• From: zleyk at tca.net
• Date: Sun, 21 Apr 2002 06:12:43 -0400 (EDT)
• References: <200204200649.CAA22309@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

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

• Follow-Ups:
  • Re: Re: Integrate malfunction?
    • From: zleyk@tca.net

• References:
  • Integrate malfunction?
    • From: "Steve Story" <sbstory@unity.ncsu.edu>