Re: Re: Solve Feature?

*To*: mathgroup at smc.vnet.net*Subject*: [mg52876] Re: [mg52846] Re: Solve Feature?*From*: DrBob <drbob at bigfoot.com>*Date*: Wed, 15 Dec 2004 04:27:20 -0500 (EST)*References*: <200412100123.UAA18967@smc.vnet.net> <opsir6lsn2iz9bcq@monster.ma.dl.cox.net> <002401c4dfb5$809b6a00$1802a8c0@Pentium4> <opsiu9j4nmiz9bcq@monster.ma.dl.cox.net> <opsiu90xstiz9bcq@monster.ma.dl.cox.net> <41BE2889.8030108@wolfram.com> <200412141100.GAA24711@smc.vnet.net> <28F087EF-4E26-11D9-A0F0-000A95B4967A@mimuw.edu.pl>*Reply-to*: drbob at bigfoot.com*Sender*: owner-wri-mathgroup at wolfram.com

Is there ever any reason not to use Rationalize before passing equations to Solve or NSolve? Bobby On Wed, 15 Dec 2004 08:16:18 +0900, Andrzej Kozlowski <akoz at mimuw.edu.pl> wrote: > *This message was transferred with a trial version of CommuniGate(tm) Pro* > FindRoot will not give you a "closed formula", which certainly has its > uses, at least in algebraic cases. However, the point seems to me a > fair one, numerical algebra based on machine arithmetic seems > inherently unreliable and I am also not sure it is worth the trouble. > Exact methods are much too slow except for relatively "small problems". > But significance arithmetic works well here and seems reasonably fast. > This, in fact, looks like the best justification I know of for having > it at all. In fact the following works perfectly: > > > {c1, c5} = {(-50 + x)^2 + (-50 + y)^2 == > 156.25`16, > (-50.00000000000002`16 + x)^2 + > (-37.49999999999999`16 + y)^2 == > 156.25`16}; > > > NSolve[%, {x, y}] > > > { > {x -> > 39.1746824526924698688078348607`14.855117414427017, > y -> > 43.749999999999977679491924308049646499`29.6989700043\ > 3602}, > {x -> 60.8253175473075501311921649336`15.04619631453911\ > 1, y -> > 43.750000000000012320508075692150353501`29.6989700043\ > 3602}} > > > {c1, c5} /. % > > > {{True, True}, {True, True}} > > > What's more, going down to much lower precision still gives perfect > answers: > > In[14]:= > {c1, c5} = {(-50 + x)^2 + (-50 + y)^2 == > 156.25`8, > (-50.00000000000002`8 + x)^2 + > (-37.49999999999999`8 + y)^2 == > 156.25`8}; > > In[15]:= > NSolve[%, {x, y}] > > Out[15]= > {{x -> 39.1746824526945169247`6.698970004336016, > y -> 43.75`6.999999999999997}, > {x -> 60.8253175473054829832`6.698970004336016, > y -> 43.75`6.999999999999997}} > > In[13]:= > {c1, c5} /. % > > Out[13]= > {{True, True}, {True, True}} > > > This is making me almost enthusiastic about significance arithmetic. > Wouldn't it be a good idea perhaps to dispense with machine arithmetic > for this type of problems entirely? > > Andrzej Kozlowski > > > On 14 Dec 2004, at 20:00, DrBob wrote: > >> *This message was transferred with a trial version of CommuniGate(tm) >> Pro* >> OK, but it's a very simple polynomial system; if Solve and NSolve give >> wildly wrong answers without notifying us, it's difficult to justify >> ever using them at all (with machine-precision coefficients). >> >> I think I'm forced to use FindRoot virtually every time, already; I >> just hadn't thought about it. >> >> Bobby >> >> On Mon, 13 Dec 2004 17:40:57 -0600, Daniel Lichtblau >> <danl at wolfram.com> wrote: >> >>> DrBob wrote: >>>> By the way, it gets even worse if I use NSolve: >>>> >>>> solution = NSolve[{c1, c5}] >>>> {c1, c5} /. solution >>>> Apply[Subtract, {c1, c5}, {1}] /. solution >>>> >>>> {{x -> -4.27759132*^8, >>>> y -> 43.74999927054117}, >>>> {x -> 4.2775924*^8, >>>> y -> 43.75000072945884}} >>>> {{False, False}, {False, False}} >>>> {{1.82977917785309*^17, >>>> 1.82977917785309*^17}, >>>> {1.8297792462945597*^17, >>>> 1.8297792462945597*^17}} >>>> >>>> These are simple polynomial equations, so this is really puzzling. >>>> >>>> Bobby >>>> >>>> On Sat, 11 Dec 2004 17:15:30 -0600, DrBob <drbob at bigfoot.com> wrote: >>>> >>>>> HELP!! Here are two very simple quadratic equations (circles): >>>>> >>>>> {c1, c5} = {(-50 + x)^2 + (-50 + y)^2 == 156.25, >>>>> (-50.00000000000002 + x)^2 + >>>>> (-37.49999999999999 + y)^2 == 156.25}; >>>>> >>>>> If we rationalize before solving, we get accurate solutions: >>>>> >>>>> solution = Solve@Rationalize@{c1, c5} >>>>> solution // N >>>>> {c1, c5} /. solution >>>>> Subtract @@@ {c1, c5} /. solution >>>>> >>>>> {{x -> (25/4)*(8 - Sqrt[3]), y -> 175/4}, >>>>> {x -> (25/4)*(8 + Sqrt[3]), y -> 175/4}} >>>>> {{x -> 39.17468245269452, y -> 43.75}, >>>>> {x -> 60.82531754730548, y -> 43.75}} >>>>> {{True, True}, {True, True}} >>>>> {{-4.263256414560601*^-14, 4.973799150320701*^-13}, >>>>> {-4.263256414560601*^-14, -4.263256414560601*^-13}} >>>>> >>>>> You can check this visually with ImplicitPlot: >>>>> >>>>> ImplicitPlot[{c1, c5}, {x, -40, 70}] >>>>> >>>>> But if we solve without Rationalize, we get wildly inaccurate >>>>> results: >>>>> >>>>> solution = Solve@{c1, c5} >>>>> {c1, c5} /. solution >>>>> Subtract @@@ {c1, c5} /. solution >>>>> >>>>> {{x -> 16., y -> 43.74999999999993}, >>>>> {x -> 84., y -> 43.75000000000004}} >>>>> {{False, False}, {False, False}} >>>>> {{1038.812500000001, 1038.8125000000005}, >>>>> {1038.8124999999995, 1038.8124999999993}} >>>>> >>>>> I'm using version 5.1. >>>>> >>>>> Bobby >>> >>> >>> The behavior of Solve with these approximate inputs is subject to the >>> vagaries of how it creates an infinite precision problem. As this one >>> is >>> poorly conditioned it is not too surprising that it does not give a >>> good >>> result. I'll point out that with VerifySolutions->True Solve will use >>> significance arithmetic throughout, and this is sufficient (in this >>> case) to get a viable result. >>> >>> The NSolve result is more troublesome. I tracked it to some difficulty >>> in assessing either of two important facts. >>> >>> (1) The precision used in an eigen decomposition is not adequate for >>> the >>> problem at hand. >>> >>> (2) The result gives a large residual for the inputs when normalized >>> in >>> what I believe to be a reasonable way (obviously we can make the >>> residual as large or small as we please by suitably rescaling the >>> input). >>> >>> It is not too hard to fix the code for either of these but there are >>> drawbacks. >>> >>> The second issue came about from a poor choice on my part of what to >>> check: I used a numerical Groebner basis instead of the original >>> inputs. >>> So this can (and will) be addressed. In those cases where it catches >>> trouble (as in this example) the fix is to redo certain computations >>> at >>> higher precision. As the code for this is already present the only >>> issue >>> is in the detection itself. The down side is that I am not certain it >>> will be adequate to assess poor solutions in general, though it works >>> reasonably well for this example. >>> >>> For the first issue I am right now at a loss. It seems that the >>> numeric >>> EigenXXX code gives a very poor result for a particular 2x2 matrix >>> that >>> appears in the course of the NSolve code. It is not too hard to figure >>> out for a given matrix that the numeric eigen decomposition might be >>> troublesome. Unfortunately that involves finding numeric singular >>> values, and I believe that tends to be typically more expensive than >>> the >>> EigenXXX computation itself. So the problem becomes one of efficiency: >>> detection of the problem cases could make all NSolve calls slower just >>> to find the rare bad apples. >>> >>> I'll give this some more thought. >>> >>> Another possible issue is in how the numeric Groebner basis gets >>> computed. This is something I also need to look at more closely. It >>> appears that for one variable order the result is quite tame, but not >>> for the other (NSolve[polys, {y,x}] works fine). I'm not yet certain >>> but >>> I think the nuemric instability at this step is essentially >>> unavoidable. >>> I believe this to be the case because I can find "nearby" polynomial >>> systems for which a problematic Groebner basis will emerge (and in >>> fact >>> that happens in this example). It is not actually incorrect, but >>> causes >>> numeric problems at low precision for the eigen decomposition on a >>> matrix formed from that Groebner basis. >>> >>> Whatever the outcome of investigating that numeric Groebner basis >>> part, >>> EigenXXX computation needs to sort out how to handle the resulting >>> matrix. So this part at least is a bug to be fixed. >>> >>> >>> Daniel Lichtblau >>> Wolfram Research >>> >>> >>> >>> >> >> >> >> -- >> DrBob at bigfoot.com >> www.eclecticdreams.net >> > > > > -- DrBob at bigfoot.com www.eclecticdreams.net

**References**:**Solve bug?***From:*paul@selfreferral.com

**Re: Solve Feature?***From:*DrBob <drbob@bigfoot.com>

**Re: Re: Solve Feature?**

**Re: Re: multiple outputs from a function**

**Re: Re: Solve Feature?**

**Re: Re: Solve Feature?**