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?