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

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.

```

• Prev by Date: Re: Setting up a probability problem
• Next by Date: Re: Function to transform lists