Re: solving third order equation: bad result
- To: mathgroup at smc.vnet.net
- Subject: [mg22163] Re: [mg22152] solving third order equation: bad result
- From: Andrzej Kozlowski <andrzej at platon.c.u-tokyo.ac.jp>
- Date: Thu, 17 Feb 2000 01:23:51 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
on 00.2.16 4:35 PM, Wolter Kaper at kaper at chem.uva.nl wrote:
> Dear Mathgroup members,
>
> When I solve the following equation for V, in terms of symbolic P and then
> substitute a value for P, I get a result that is wrong. (The right result is
> returned when I solve the equation with P substituted beforehand.)
>
> P == RT /(V - b) - a/V^2 /. {a -> 0.142666, b -> 3.913*10^(-5), RT -> 2494})
>
> The following notebook demonstrates the problem and compares results after
> the equation has been rewritten in standard (third order polynomial ==0)
> form. In that case the right result is returned.
> My questions are: does anybody know why this behaviour occurs, has anybody
> seen it before, should I report a bug to Wolfram research?
>
> Thanks beforehand,
> Wolter Kaper
> Univ. Of Amsterdam, dept. of Chemistry
>
> Notebook follows.
I can't really explain the exact cause of the problem, but here are a
couple of observations. First, if you convert all the numbers to arbitrary
precision numbers using Rationalize and then use Solve you will get the
right answer.
This is what happens:
In[65]:=
Clear[eq1]
In[66]:=
eq1 = P == T/(V - b) - a/V^2 /.
{a -> Rationalize[0.142666, 0],
b -> Rationalize[3.913/10^5, 0], T -> 2494}
Out[66]=
2494 71333
P == ---------------- - ---------
3913 2
-(---------) + V 500000 V
100000000
In[67]:=
OplA = N[Solve[eq1, V]];
In[68]:=
OplA /. P -> 4*10^8
Out[68]=
{{V -> 0.0000444101}, {V ->
-7
4.77432 10 + 0.0000177209 I},
-7
{V -> 4.77432 10 - 0.0000177209 I}}
So this suggests that this is indeed a problem with numerical accuracy.
But now, here is something interesting. Suppose we replace the equation eq1
by a trivially equivalent system of two simultaneous equations:
In[69]:=
Clear[eq1]
In[70]:=
eq1 = P == T/(V - b) - a/V^2 /. {a -> 0.142666,
b -> 3.913/10^5, T -> 2494}
Out[70]=
2494 0.142666
P == --------------- - --------
-0.00003913 + V 2
V
In[71]:=
OplA = Solve[{eq1, T == P}, {V, T}];
In[72]:=
OplA /. P -> 4*10^8
Out[72]=
8 -7
{{T -> 4. 10 , V -> 4.77432 10 - 0.0000177209 I},
8 -7
{T -> 4. 10 , V -> 4.77432 10 + 0.0000177209 I},
8
{T -> 4. 10 , V -> 0.0000444101}}
Again the right answer. Now I know that when dealing with a system of
algebraic equations Mathematica uses Groebner basis. I believe that Groebner
basis works badly over inexact numbers so perhaps Mathematica treats in this
case converts the finite precision numbers to infinite precision ones and
the two above approaches are really the equivalent. That is only my
hypothesis. Daniel Lichtblau of course knows the true answer and I am sure
he will let us know.