Re: [Maxima-discuss] simplify_sum leads to a wrong value

[Robert D. <rob...@gm...>]

On 2020-04-07, philippe <rou...@gm...> wrote:

> sum((k*sin((k*%pi)/3))/(4*k^2-1),k,1,inf)
>
> I know it converge to (sqrt(3)*%pi)/16=0.3400873807939158...
>
> simplify_sum find the limit in terms of erf(%i) but floating point
> approximation is far from the expected value : 2.564242426927912...

Hmm, I don't know where simplify_sum went off the tracks. I didn't look
carefully at what simplify_sum is doing, but if I understand correctly
one of the major algorithms it employs is the Zeilberger/Gosper
algorithm which works, from what I understand, on summands which are
rational functions. I speculate that it's the sin(k*%pi/3) which causes
trouble, but that's just a guess.

I find that I can make progress by replacing sin(k*%pi/3) with
+/- sqrt(3)/2 as follows and then applying simplify_sum later in the
process. First I separate the terms of the summation into three groups,
for k = 3*m + 1, 3*m + 2, and 3*m + 3.

  (%i3) sum((k*sin((k*%pi)/3))/(4*k^2-1),k,1,inf) $
  (%i4) summand: first(%);
  (%o4) (k*sin((%pi*k)/3))/(4*k^2-1)
  (%i5) summand1: summand, k = 3*m + 1;
  (%o5) ((3*m+1)*sin((%pi*(3*m+1))/3))/(4*(3*m+1)^2-1)
  (%i6) summand2: summand, k = 3*m + 2;
  (%o6) ((3*m+2)*sin((%pi*(3*m+2))/3))/(4*(3*m+2)^2-1)
  (%i7) summand3: summand, k = 3*m + 3;
  (%o7) ((3*m+3)*sin((%pi*(3*m+3))/3))/(4*(3*m+3)^2-1)

Now adding up the three summands is equivalent to the original.

  (%i8) sum (summand1 + summand2 + summand3, m, 0, inf);
  (%o8) 'sum(((3*m+3)*sin((%pi*(3*m+3))/3))/(4*(3*m+3)^2-1)
              +((3*m+2)*sin((%pi*(3*m+2))/3))/(4*(3*m+2)^2-1)
              +((3*m+1)*sin((%pi*(3*m+1))/3))/(4*(3*m+1)^2-1),m,0,
             inf)

Massage this a little to squeeze the trigonometric term out of it.

  (%i9) expand (%);
  (%o9) 'sum((3*m*sin(%pi*m+(2*%pi)/3))/(36*m^2+48*m+15)
              +(2*sin(%pi*m+(2*%pi)/3))/(36*m^2+48*m+15)
              +(3*m*sin(%pi*m+%pi/3))/(36*m^2+24*m+3)
              +sin(%pi*m+%pi/3)/(36*m^2+24*m+3)
              -(3*m*sin(%pi*m))/(36*m^2+72*m+35)
              -(3*sin(%pi*m))/(36*m^2+72*m+35),m,0,inf)
  (%i10) ''%;
  (%o10) 'sum((3^(3/2)*m*(-1)^m)/(2*(36*m^2+48*m+15))
               +(sqrt(3)*(-1)^m)/(36*m^2+48*m+15)
               +(3^(3/2)*m*(-1)^m)/(2*(36*m^2+24*m+3))
               +(sqrt(3)*(-1)^m)/(2*(36*m^2+24*m+3)),m,0,inf)

Now try simplify_sum again.

  (%i11) load (simplify_sum) $
  (%i12) simplify_sum (%o10);
  (%o12) (%pi/2-3*log(2))/(16*sqrt(3))
   -((-3*log(2))-%pi/2)/(16*sqrt(3))
   +(psi[0](11/12)+%gamma)/(32*sqrt(3))
   +(psi[0](7/12)+%gamma)/(32*sqrt(3))
   -(psi[0](5/12)+%gamma)/(32*sqrt(3))
   -(psi[0](1/12)+%gamma)/(32*sqrt(3))

Encourage Maxima to apply identities for rational arguments of psi[0].

  (%i13) maxpsifracnum: maxpsifracdenom: 12 $
  (%i14) ''%o12;
  (%o14) ((-(%pi*cot((11*%pi)/12))/2)-(sqrt(3)*log(sqrt(3)+2))/2
                                     +(sqrt(3)*log(2-sqrt(3)))/2
                                     -log(12)-log(4)/2-log(3)/2)
   /(32*sqrt(3))
   +((-(%pi*cot((7*%pi)/12))/2)+(sqrt(3)*log(sqrt(3)+2))/2
                               -(sqrt(3)*log(2-sqrt(3)))/2-log(12)
                               -log(4)/2-log(3)/2)
    /(32*sqrt(3))
   -((-(%pi*cot((5*%pi)/12))/2)+(sqrt(3)*log(sqrt(3)+2))/2
                               -(sqrt(3)*log(2-sqrt(3)))/2-log(12)
                               -log(4)/2-log(3)/2)
    /(32*sqrt(3))
   -((-(%pi*cot(%pi/12))/2)-(sqrt(3)*log(sqrt(3)+2))/2
                           +(sqrt(3)*log(2-sqrt(3)))/2-log(12)
                           -log(4)/2-log(3)/2)
    /(32*sqrt(3))+(%pi/2-3*log(2))/(16*sqrt(3))
   -((-3*log(2))-%pi/2)/(16*sqrt(3))

OK, let's collect common terms and call it a day.

  (%i15) ratsimp (%);
  (%o15) -(sqrt(3)*%pi*cot((11*%pi)/12)
   +sqrt(3)*%pi*cot((7*%pi)/12)-sqrt(3)*%pi*cot((5*%pi)/12)
   -sqrt(3)*%pi*cot(%pi/12)-4*sqrt(3)*%pi)
   /192
  (%i16) %, numer;
  (%o16) 0.3400873807939157

It's reassuring to see the expected numerical value.

It's still a bug that simplify_sum (or maybe some other simplification
step, I dunno) handled the original summation incorrectly. Maybe someone
will create a bug report about it.

best,

Robert Dodier