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
>
> 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

```

• Prev by Date: Re: Linux: how compatible with add-ons ?
• Next by Date: Re: What is happening here?