Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Re: bug ?

• To: mathgroup at smc.vnet.net
• Subject: [mg122484] Re: Re: bug ?
• From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
• Date: Sat, 29 Oct 2011 07:15:00 -0400 (EDT)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com
• References: <j862b6$5ls$1@smc.vnet.net> <201110280933.FAA20663@smc.vnet.net>

On 28 Oct 2011, at 11:33, Peter Pein wrote:

> Am 25.10.2011 12:17, schrieb swiss:
>> 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
>>
>
> Well, if Cancel[], Factor[] and/or Apart[] can't do that (and I would
> consider this a bug), let's do it by hand:
>
> In[1]:= {n, d} = Through[{Numerator, Denominator}[
>    Together[-(((s*(4*s - 5)^2 - (1/8)*(5 + 3*Sqrt[17])*(4*s + 1))*
>        ((4*s - 5)^2 - (4/8)*(5 + 3*Sqrt[17])*(4*s - 7)))/
>       (16*s^2 - 32*s - 1)^2)]]]
>
>    s0 = (1/4)*(4 + Sqrt[17]);
>
> Out[1]= {(-(85 + 21*Sqrt[17] - 100*s - 12*Sqrt[17]*s + 32*s^2))*
>   (-5 - 3*Sqrt[17] + 180*s - 12*Sqrt[17]*s - 320*s^2 + 128*s^3),
>  16*(-1 - 32*s + 16*s^2)^2}
>
> In[3]:= gcd = PolynomialGCD[n, d, Extension -> Automatic]
> Out[3]= (4 + Sqrt[17] - 4*s)^2
>
> In[4]:= pqr = (PolynomialQuotientRemainder[#1, gcd, s] & ) /@ {n, d}
> Out[4]= {{408 - 88*Sqrt[17] + (-928 + 144*Sqrt[17])*s +
>         (928 - 32*Sqrt[17])*s^2 - 256*s^3, 0},
>        {528 - 128*Sqrt[17]    +(-512 + 128*Sqrt[17])*s + 256*s^2, 0}}
>
> In[5]:= y = Divide @@ pqr[[All,1]];
> In[6]:= Limit[y, s -> s0]
>        N[%]
>
> Out[6]= (1/544)*(331 - 19*Sqrt[17])
> Out[7]= 0.4644503549876185
>
> In[8]:= y /. s -> s0
>        N[%]
>        Simplify[%%]
>
> Out[8]= (408 - 88*Sqrt[17] +(1/16)*(928 - 32*Sqrt[17])*(4 + Sqrt[17])^2
> -  4*(4 + Sqrt[17])^3 + (1/4)*(4 + Sqrt[17])*(-928 + 144*Sqrt[17]))/
>  (528 - 128*Sqrt[17] + 16*(4 + Sqrt[17])^2 + (1/4)*(4 + Sqrt[17])*
>  (-512 + 128*Sqrt[17]))
> Out[9]= 0.4644503549876185
> Out[10]= (1/544)*(331 - 19*Sqrt[17])
>
> In[11]:= y /. s -> N[s0]
> Out[11]= 0.46445035498761833
>
> In[12]:= RootApproximant[%, 2]
> Out[12]= (1/544)*(331 - 19*Sqrt[17])
>
>
> IMHO these steps should belong to Cancel[].
>
> Peter
>

I have not read this thread carefully so maybe I am have missed something but the simplest thing to do seems to me to be:

ff = -(((s*(4*s - 5)^2 - (1/8)*(3*Sqrt[17] + 5)*(4*s + 1))*((4*s - 5)^2 - (4/8)*(3*Sqrt[17] + 5)*(4*s - 7)))/(16*s^2 - 32*s - 1)^2);
gg = Factor[ff, Extension -> (s /. {ToRules[Roots[ff == 0, s]]})];

Simplify[gg /. s -> (1/4)*(Sqrt[17] + 4)]
(1/544)*(331- 19*Sqrt[17] )

Of course I agree that Cancel should do this automatically.

Andrzej Kozlowski

• References:
  • Re: bug ?
    • From: Peter Pein <petsie@dordos.net>