MathGroup Archive 2011

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

Search the Archive

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[50]:= Limit[-(((s (4 s - 5)^2 - (5 + 3 Sqrt[17])/
>        8 (4 s + 1)) ((4 s - 5)^2 -
>       4 (5 + 3 Sqrt[17])/8 (4 s - 7)))/(16 s^2 - 32 s - 1)^2),
>   s ->  1/4 (4 + Sqrt[17])]
>
> Out[50]= 1/544 (331 - 19 Sqrt[17])
>
> In[51]:= 1/544 (331 - 19 Sqrt[17]) // N
>
> Out[51]= 0.46445
>
> In[52]:= -(((s (4 s - 5)^2 - (5 + 3 Sqrt[17])/
>         8 (4 s + 1)) ((4 s - 5)^2 -
>        4 (5 + 3 Sqrt[17])/8 (4 s - 7)))/(16 s^2 - 32 s - 1)^2) /.
>    s ->  1/4 (4 + Sqrt[17]) // Simplify
>
> Out[52]= 1/16 (1 + 2 Sqrt[17])
>
> In[53]:= 1/16 (1 + 2 Sqrt[17]) // N
>
> Out[53]= 0.577888
>
> In[54]:= -(((s (4 s - 5)^2 - (5 + 3 Sqrt[17])/
>        8 (4 s + 1)) ((4 s - 5)^2 -
>       4 (5 + 3 Sqrt[17])/8 (4 s - 7)))/(16 s^2 - 32 s - 1)^2) /.
>   s ->  1/4 (4 + Sqrt[17.])
>
> Out[54]= 0.09375


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

Just to verify the claims about the limit:

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

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

Out[71]= (1/544)*(331 - 19*Sqrt[17])

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[66]:= num = Numerator[expr] /. s -> (1/4.)*(4 + Sqrt[17])
den = Denominator[expr] /. s -> (1/4.)*(4 + Sqrt[17])
num/den

Out[66]= 1.8932661725304283*^-29

Out[67]= 2.0194839173657902*^-28

Out[68]= 0.09375

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

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

Out[80]= 0

Out[81]= 0

Out[82]= (1/16)*(1 + 2*Sqrt[17])

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 ?