Re: solving third order equation: bad result
- To: mathgroup at smc.vnet.net
- Subject: [mg22190] Re: [mg22152] solving third order equation: bad result
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Thu, 17 Feb 2000 01:24:27 -0500 (EST)
- References: <B4D0A1C6.4E1B%andrzej@platon.c.u-tokyo.ac.jp>
- Sender: owner-wri-mathgroup at wolfram.com
Andrzej Kozlowski wrote:
>
> 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.
While such equations can suffer from numeric stability considerations on
back-substitution for parameters, thus was not the issue in this case.
The behavior of the one-equation-one-variable example is indeed a bug,
one that lies in the internals of root-extraction code. This bug is a
bit subtle and was basically elluded in going to the essentially
equivalent two-equation-two-variable problem. In any case I expect to
have it fixed for the next kernel release.
As surmised, Solve will rationalize its inputs.
While I'm on the podium I'll note that I'm a bit crunched for time
nowadays so I do not look at notebooks to test candidate kernel bugs. It
simply takes too long to find the appropriate cells and convert to
InputForm that I can use within a kernel session run inside a debugger.
Had Andrzej not done this conversion it is doubtful whether the bug
would have been looked at for days to weeks at a minimum (it would get
reported to WRI, go through our SQA department, someone there would
perform said conversion, they'd file a report, it would get forwarded to
me....). Moral: Provide all InputForm needed to replicate the
misbehavior, or else give a convincing reason why this cannot be done
(or else be prepared to exercise alot of patience).
Daniel Lichtblau
Wolfram Research