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