       Re: bug ?

• To: mathgroup at smc.vnet.net
• Subject: [mg122347] Re: bug ?
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Wed, 26 Oct 2011 17:37:30 -0400 (EDT)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com
• References: <201110251016.GAA05769@smc.vnet.net>

On 10/25/2011 05:16 AM, swiss wrote:
> Following rational function has a double pole matched by a double root at (4+\sqrt{17})/4. The first result is correct, the other ones are false, but come without warning. Can somebody explain this to me?
>
> In:= Limit[-(((s (4 s - 5)^2 - (5 + 3 Sqrt)/
>        8 (4 s + 1)) ((4 s - 5)^2 -
>       4 (5 + 3 Sqrt)/8 (4 s - 7)))/(16 s^2 - 32 s - 1)^2),
>   s ->  1/4 (4 + Sqrt)]
>
> Out= 1/544 (331 - 19 Sqrt)
>
> In:= 1/544 (331 - 19 Sqrt) // N
>
> Out= 0.46445
>
> In:= -(((s (4 s - 5)^2 - (5 + 3 Sqrt)/
>         8 (4 s + 1)) ((4 s - 5)^2 -
>        4 (5 + 3 Sqrt)/8 (4 s - 7)))/(16 s^2 - 32 s - 1)^2) /.
>    s ->  1/4 (4 + Sqrt) // Simplify
>
> Out= 1/16 (1 + 2 Sqrt)
>
> In:= 1/16 (1 + 2 Sqrt) // N
>
> Out= 0.577888
>
> In:= -(((s (4 s - 5)^2 - (5 + 3 Sqrt)/
>        8 (4 s + 1)) ((4 s - 5)^2 -
>       4 (5 + 3 Sqrt)/8 (4 s - 7)))/(16 s^2 - 32 s - 1)^2) /.
>   s ->  1/4 (4 + Sqrt[17.])
>
> Out= 0.09375

In:= expr = -(((s (4 s - 5)^2 - (5 + 3 Sqrt)/
8 (4 s + 1)) ((4 s - 5)^2 -
4 (5 + 3 Sqrt)/8 (4 s - 7)))/(16 s^2 - 32 s - 1)^2);

Just to verify the claims about the limit:

e2 = Together[expr, Extension -> Sqrt]
e3 = Simplify[e2 /. s -> 1/4 (4 + Sqrt)]

Out= ((17 + Sqrt - 8*s)*(31 - 7*Sqrt - 48*s +
8*Sqrt*s + 32*s^2))/
(16*(-4 + Sqrt + 4*s)^2)

Out= (1/544)*(331 - 19*Sqrt)

Fine so far.

The below numerical computation should give an idea of where that .09375
arises. It is simply a bady conditioned evaluation of a quotient of
approximate zeros. They would be exactly zero, where we doing exact
arithmetic.

In:= num = Numerator[expr] /. s -> (1/4.)*(4 + Sqrt)
den = Denominator[expr] /. s -> (1/4.)*(4 + Sqrt)
num/den

Out= 1.8932661725304283*^-29

Out= 2.0194839173657902*^-28

Out= 0.09375

The bad exact substitution is more of a problem. We see below that it is
certainly an indeterminate form.

In:= numex = Numerator[expr] /. s -> (1/4)*(4 + Sqrt);
denex = Denominator[expr] /. s -> (1/4)*(4 + Sqrt);
Simplify[numex]
Simplify[denex]
Simplify[numex/denex]

Out= 0

Out= 0

Out= (1/16)*(1 + 2*Sqrt)

The upshot is Mathematica is not recognizing a 0/0 situation and is,
somewhere, removing a zero rather than finding an indeterminate. This
can happen in internals of algebraic manipulation, wherein pieces can be
worked on in a nested fashion and some zeroesmight be more "hidden" than
others. In this case I do not know if it shows a bug or a feature (I
will have a closer look though).

Daniel Lichtblau
Wolfram Research



• References:
• bug ?
• From: swiss <gregoire.nicollier@hevs.ch>
• Prev by Date: Re: how do you prevent numerator expansion when using 'Together'?
• Next by Date: Re: Large control loops
• Previous by thread: Re: bug ?
• Next by thread: Re: bug ?