Re: Re: Integrate malfunction?
- To: mathgroup at smc.vnet.net
- Subject: [mg33909] Re: [mg33886] Re: [mg33860] Integrate malfunction?
- From: zleyk at tca.net
- Date: Mon, 22 Apr 2002 00:57:49 -0400 (EDT)
- References: <200204200649.CAA22309@smc.vnet.net> <200204211012.GAA28553@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
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
> > _______________
> >
- References:
- Integrate malfunction?
- From: "Steve Story" <sbstory@unity.ncsu.edu>
- Re: Integrate malfunction?
- From: zleyk@tca.net
- Integrate malfunction?