Re: numerics

*To*: mathgroup at smc.vnet.net*Subject*: [mg27285] Re: numerics*From*: "John Doty" <jpd at w-d.org>*Date*: Sat, 17 Feb 2001 03:31:03 -0500 (EST)*Organization*: Wampler-Doty Family*References*: <96iqph$d9j@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

In article <96iqph$d9j at smc.vnet.net>, Matt.Johnson at autolivasp.com wrote: > mathematica gurus-- > > Here's a frustrating problem, which seems to be simple. I can calculate > a matrix by hand and enter it into Solve or Reduce and get the correct > relationships: > > In[108]:= mat = Partition[{-0.2, 0.1, 0, 0.1, -0.3, 0.1, 0.1, 0.2, > -0.1}, 3]; Reduce[mat.{x1, x2, x3} == {0, 0, 0}, {x1, x2, x3}] Out[109]= > x1 == 0.5 x2 && x3 == 2.5 x2 > > However, if I try to manipulate the matrix in Mathematica then solve, it > doesn't work: > > In[110]:= mat1 = Partition[{0.8, 0.1, 0.1, 0.1, 0.7, 0.2, 0, 0.1, 0.9}, > 3]; imat = IdentityMatrix[Length[mat1]]; newmat = Transpose[mat1] - > imat; In[113]:= newmat == mat Out[113]= True In[114]:= newmat === mat > Out[114]= False In[115]:= Reduce[newmat.{x1, x2, x3} == {0, 0, 0}, {x1, > x2, x3}] Out[115]= x1 == 0. && x2 == 0. && x3 == 0. You can't reliably do this with approximate numbers: In[1]:= x = 0.2 Out[1]= 0.2 In[2]:= y = 1 - 0.8 Out[2]= 0.2 x and y are apparently equal and: In[3]:= x == y Out[3]= True appears to confirm this. However: In[4]:= x-y == 0 Out[4]= False In fact, because of rounding error: In[5]:= x-y -17 Out[5]= 5.55112 10 The equality test with "x == y" yields "True" because "==" allows for some rounding error. Your "imat" turned out to be exactly singular when reduced. However, "newmat" did not. In general, you should not count on singularity of approximate matrices: your result will depend on rounding error and the order in which the reduction is computed. This is not a particular difficulty with Mathematica: it is a problem that any approximate calculation can encounter. Mathematica, however, offers a way around it. If you replace your approximate numbers with exact numbers (e.g. 1/5 in place of 0.2), you will get the result you want. -- | John Doty "You can't confuse me, that's my job." | Home: jpd at w-d.org | Work: jpd at space.mit.edu