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>
- Solve bug?