Re: Solve / NSolve
- To: mathgroup at smc.vnet.net
- Subject: [mg95386] Re: Solve / NSolve
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Sat, 17 Jan 2009 05:29:27 -0500 (EST)
- References: <gkgouo$3sg$1@smc.vnet.net> <gkppsb$dot$1@smc.vnet.net>
On Jan 16, 5:08 am, SigmundV <sigmu... at gmail.com> wrote:
> Thanks to Andrzej and Jean-Marc for the discussion on the differences
> between Solve and NSolve. Have any of you tried to run the code i
> provided (including the missing function definition which I provided
> later)? I'm aware that Solve and NSolve do not use the same algorithm
> to solve an equation, but it does not seem "logical" to me that they
> don't give the same solution to the same system of algebraic
> equations. What could be the reason for this? I'm also concerned witht
> he speed. Solve is much faster than NSolve for my problem. How come?
> In addition to the code I provided earlier, I just tried with another
> example of a system of two algebraic equations:
>
> {1.1 x^2 + y^2 == 1.0, x + 1.1 y == 0.0}
> Solve[%, {x, y}] // Timing
> NSolve[%%, {x, y}] // Timing
>
> When evaluating this cell you'll see that Solve is substantially
> faster than NSolve, but they yield the same solution. However, in my
> first example they do not yield the same solution in all cases.
>
> Hopefully some of you can shed some light on all this. It could also
> be me comparing apples and oranges. Let me hear your thoughts on all
> this.
>
> Kind regards,
> Sigmund
This is a bit of a tangle. One thing you probably do not realize is
that people using version 7 get different behavior from you.
Specifically, the Solve variant will give a slew of error messages due
to badness in your Drop code. I isolated a simple example that shows
this difference, and I'll cover it below because it also serves, to an
extent, to show why NSolve is safer for a problem like this even in
absence of such error messages.
polys = {100*(-0.11 + x)^2 + 25*(-0.21 + y)^2 - 1,
(-1.8774107494749928*(-0.11 + x))/10^13 -
10.971072010497574*(-0.21 + y)};
Mathematica 6:
In[12]:= Solve[polys==0, {x, y}]
Out[12]= {{x -> 0.0785768, y -> 0.21}, {x -> 0.139692, y -> 0.21}}
Mathematica 7:
In[13]:= Solve[polys==0, {x, y}]
Out[13]= {{x -> 0.109135, y -> 0.21}, {x -> 0.109135, y -> 0.21}}
Observe that the latter claims a double root whereas the x coordinates
in the former are well separated.
So let me say a bit about what Solve does. It starts by rationalizing
all coefficients and computing an exact Groebner basis. In both
versions this is similarly handled, although the rationalization is
different in a way that is irrelevant for the present discussion. It
then feeds a univariate polynomial to a numerical root finder. In
version 6 and prior this was Jenkins-Traub. In version 7 I believe
this is no longer the case, at least by default.
Here is the univariate polynomial produced by one of the versions. I
will show the resulting roots in both versions.
upoly = -1.049020730058627*^75*y + 2.4976684049014928*^75*y^2 +
1.1014717665615586*^74;
Mathematica 6:
In[17]:= InputForm[Roots[upoly==0, y]]
Out[17]//InputForm= y == 0.20999999999999952 || y ==
0.21000000000000058
Mathematica 7:
In[16]:= InputForm[Roots[upoly==0, y]]
Out[16]//InputForm= y == 0.21000000000000002 || y ==
0.21000000000000002
In the first case the roots are separated slightly, and in the second
they are not. This extremely small difference is what messes up later
computations of x; the small root difference gets magnified by some
large coefficients in the polynomial(s) that relate x to y, so the x
values get separated in version 6. Clearly when the y values are
identical this will not happen.
I do not view this as a bug in the new root finding code. It may well
be a bug fix, in that the code now only tries to return digits that
are justifiable by the input precision. The bad handling, if any, is
in Solve, in rationalizing and then finding a Groebner basis. Reason
being, this can make numerator coefficients huge, thus introducing
numeric instability. Huge coefficients pose no problem in principle,
when doing exact computation (the domain of Solve). For approximate
computation this introduces all sorts of trouble with cancellation
error downstream (e.g. in numeric handling of Roots).
As for speed, I can say a bit. This particular equation set is fast to
handle, and Solve manages to beat NSolve in timing. That is by no
means always the case, and is quite rare for most polynomial systems
that come my way. Moreover even in version 6, the Solve solution is
garbage (check the residuals), hence your bad pictures. It is just not
garbage with spurious multiplicity, hence does not run afoul of the
Drop code. In contrast, the NSolve solution gives small residuals.
Probably worth spending a bit more time in that code, if it is better
suited to handle the problem.
Daniel Lichtblau
Wolfram Research