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