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

Re: When does x/x not equal 1?

• To: mathgroup at smc.vnet.net
• Subject: [mg32426] Re: When does x/x not equal 1?
• From: bghiggins at ucdavis.edu (Brian Higgins)
• Date: Sat, 19 Jan 2002 19:03:24 -0500 (EST)
• References: <a2b4ib$lja$1@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

Kevin,
Would Chop do the trick?

In[202]:=N[Sqrt[1 - (1 + 1/10000000000000000000000000)]]

Out[202]=0. + 3.162277660168379*^-13*I
 
Now wrapping the expression with Chop

In[203]:=
Chop[N[Sqrt[1 - (1 + 1/10000000000000000000000000)]]]

Out[203]=0

Cheers

Brian



mjperson at mit.edu (Michael J Person) wrote in message news:<a2b4ib$lja$1 at smc.vnet.net>...
> Greetings all:
> 
> 	  I'm having a bit of a problem with Math 4.0.1.0 on Mac OS 9.1...
> 	  
> 	  Specifically we have a lot of large planetary calculations that only
> occasionally produce piles of complex numbers of the form:
> 
> 103.3 + 3.4432 *10^-39 I
> 
>       These tiny complex pieces then get carried all around in our 
> calculations and destroy the formatting of our outputs and things 
> even though the real part seems to be a quite reasonable answer to 
> our problems.
> 
>     Hunting it down piece by piece, we seem to have traced it to places
> where we're taking the square roots of things that include ratios of
> numbers that should be equaling 1 at the boundaries but aren't.
> 
> Sqrt[1-x/x] -> tiny complex number instead of 0!
> 
> 	    It seem to be that x/x that's causing the problem in our codes,
> but as we try to track it down, all sort of mystifying behavior starts
> to manefest.  Among the worst of it are the following sorts of things:
> 
> In[1]:=
> 	(x = (1.992+$MachineEpsilon))   //FullForm
> 	
> Out[1]//FullForm=
> 	1.9920000000000002`
> 	
> In[2]:=
> 	x/x //FullForm
> 	
> Out[2]//FullForm=
> 	0.9999999999999999`
> 	
> 	Now, I'd guess that wonky things start happening when you start
> fooling around with numbers near $MachineEpsilon, but shouldn't
> x/x = 1?  Worse, is when we tried to fix the problem in various ways,
> we came up with terribly confusing things like the following:
> 
> In[3]:=
> 	x/x //FullForm
> 	
> Out[3]//FullForm=
> 	0.9999999999999999`
> 	
> (*That's just like above.)
> 
> In[4]:=
> 	Divide[x,x] //FullForm
> 	
> Out[4]//FullForm=
> 	1.`
> 	
> 	Ah ha!  That's just the sort of thing we want in our calculations.
> So we try to see what the difference is....
> 
> In[5]:=
> 	Divide[a,b] //FullForm
> 
> Out[5]//FullForm=
> 	Times[a, Power[b, -1]]
> 	
> In[6]:=
> 	a/b //FullForm
> 	
> Out[6]//FullForm=
> 	Times[a, Power[b, -1]]
> 	
> 	They're handled exactly the same way! So, can anyone explain to me 
> what's going on here, and how I can keep very small complex numbers from
> arising in all of those cases when x/x != 1.  
> 
> 	(As I mentioned these problems are completely intermittent.  In the 
> examples given above, they do not occur for 1.991, or 1.994, just 1.992
> and 1.993.  Bah!)
> 
> Thanks for anything you can tell me.
>        
> -Michael Person
> mjperson at mit.edu
> Massachusetts Institute of Technology
> Department of Earth, Atmospheric, and Planetary Sciences