MathGroup Archive 2011

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

Search the Archive

Re: bug ?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg122363] Re: bug ?
  • From: Richard Fateman <fateman at cs.berkeley.edu>
  • Date: Wed, 26 Oct 2011 17:40:25 -0400 (EDT)
  • Delivered-to: l-mathgroup@mail-archive0.wolfram.com
  • References: <j862b6$5ls$1@smc.vnet.net>

On 10/25/2011 3:17 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?

As you know, you are dividing by zero. That's why you take a limit.
If you are not doing arithmetic NOT with the exact values you think you 
have, but some other numbers more or less nearby, you get different results.


Consider a slight variation on your input 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.0000000000000000])

which returns Indeterminate.

Although

-(((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.000000000000000])

(one fewer "0") in there, returns without comment, 0.09375.

Is this a bug? I am sure that some people will call it a feature.
It is a result of compromises made in numerics.


Your example on line 52 is more interesting.

Let q =
-(((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])

Then

Simplify[q]     returns, as you show,  1/16* (1+2*Sqrt[17])
FullSimplify[q] returns the same
N[q]            returns 0.
N[q,1]          returns Indeterminate  (with warning)
Expand[q]       returns 0

Map[N,Distribute[q]]  returns -1.49737* 10^14     [A HUGE number!]

Oh, for fun, take that q value and do this:

  qq=  (q/. 17 -> r)

That is, replace all the instances of 17 with r.

Limit[qq, r->17]  is  1/16*(1+2*Sqrt[17])

which is of course not the same as output 50, but is the result of Simplify.

Some of the results about can be explained by a feature/bug that is 
likely to be in Mathematica.  That is essentially a rule that says if 
one is computing a product like A * B, and one is able to simplify A
to zero, there is no need to look at B.  The product is zero. Probably 
other programs like Mathematica do the same thing, by the way.

In this case B =  1/0, but B wasn't examined.  So that is an 
explanation.   It is a result of a heuristic in the Simplify program.

Is it a bug or a feature?  Well, it wouldn't be the first time that 
something in Mathematica is simultaneously declared to be a feature and 
produces um, mathematically questionable results.  I guess that a fan of 
Mathematica could argue that it must be a feature because this behavior 
is well known to WRI and has not been changed, and WRI would change it 
if it were a bug, so it must be a feature. :)





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




  • Follow-Ups:
    • Re: bug ?
      • From: DrMajorBob <btreat1@austin.rr.com>
  • Prev by Date: Re: Finding distribution parameters for a mixture of product distributions
  • Next by Date: Re: runs slowly?
  • Previous by thread: Re: bug ?
  • Next by thread: Re: bug ?