MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Re: Solve Feature?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg52877] Re: [mg52846] Re: Solve Feature?
  • From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
  • Date: Wed, 15 Dec 2004 04:27:26 -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> <opsi0yvxd6iz9bcq@monster.ma.dl.cox.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Certianly, if you raionalize and try to use Solve any really complex 
system will take longer than you can wait. Even a single polynomial of 
high enough degree.
NSolve it is another matter; Daniel can explain this better. It seems 
to me that rationalizing first and then using NSolve may sometiems be a 
good idea, but of ocurse you will still be using MachiePrecision or 
whatever precision you specify in NSolve.

Andrzej


On 15 Dec 2004, at 10:10, DrBob wrote:

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


  • Prev by Date: Re: Re: multiple outputs from a function
  • Next by Date: Re: Remove the higher order terms in series expansion
  • Previous by thread: Re: Re: Solve Feature?
  • Next by thread: Re: Solve Feature?