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