7.91 bug in integration in Maple 7 and correction (15.4.02)

7.91.1 Marko Horbatsch
7.91.2 Helmut Kahovec (25.4.02)

7.91.1 Marko Horbatsch

In case this bug coincides with a previous report, my apologies, but it seems that the developers are not checking their software updates with an appropriate suite of integrals, they should add the following:

int(ln(x)/(1-x),x=0..1);
 

should result in -Pi^2/6

but gives infinity (maple6 is OK on this one).

If one uses the anti-derivative, and takes the limit for the lower limit, the correct result is returned; it is however the right-hand integration limit which seems to throw maple7, because changing that to something else than 1, makes the integral work.

Of course, this hits a class of integrals, i.e., including powers of ln(x) produces similarly erroneous results.

for

f:=ln(x)^3/(1-x)^3; 
 
int(f,x=0..1);
 

isn’t wrong, but doesn’t evaluate, even though taking the limit on the anti-derivative works,

limit(int(f,x=a..1),a=0,right);
 

works out. This complaint applies to both maple7 and maple6, i.e., apparently maple does not try on these improper integrals taking the antiderivative and then performing limits to obtain the definite integral.

It is corrected with Maple 8 (U. Klein)

7.91.2 Helmut Kahovec (25.4.02)

There seems to be a bug in ‘int/cook/IIntd1c‘ of Maple7, which defines two different versions of ‘int/cook/ngritty‘. The second one evaluates the derivative of the Beta(z,w) function at w=0 although the value of w may be zero. This evaluation should be replaced by the limit lim w->0. Then we get the following results:

> `int/cook/IIntd1c`:=proc( 
    t,ta,tb,param,model,triglabel,ans,fltype,fail 
  ) 
  local cof,a,b,c,d,dl,a0,a1,f,m,w,np1,p,u,u2,r,s,s2,sp,ucomplex, 
        v,z,btype,ptype,a0type,a1type,wtype,dtype,dltype,a0fl, 
        a1fl,bfl,wfp1l,dfl,dlfl,pfm1l,B,S1,S2,D,DL,A0,A1,C1,M, 
        N,P,R,Ucplex,U2,MODEL,failtemp,fltypetemp,res,res1,res2, 
        i,sig,dc1,dc2,dc3,dc4,`int/cook/ngritty`; 
    failtemp:=0; 
    cof:=param[1]; 
    b:=param[2]; 
    u:=-op(1,param[3]); 
    u2:=-op(2,param[3]); 
    s:=op(1,param[4]); 
    s2:=op(2,param[4]); 
    dl:=param[5]; 
    m:=param[6]; 
    w:=param[7]; 
    p:=param[8]; 
    a0:=param[9]; 
    a1:=param[10]; 
    c:= param[11]; 
    r:=param[12]; 
    d:=param[13]; 
    sp:=param[14]; 
    f:=param[15]; 
    MODEL:=exp(-Ucplex*t^S1)*t^N*ln(B*t^DL)^M*cos(C1*t^R)/ 
      ((A0+A1*f^D)^P); 
    d:=`int/cook/secure`(d,'dfl','dtype'); 
    p:=`int/cook/secure`(p-1,'pfm1l','ptype')+1; 
    w:=`int/cook/secure`(w+1,'wfp1l','wtype')-1; 
    b:=`int/cook/secure`(b,'bfl','btype'); 
    dl:=`int/cook/secure`(dl,'dlfl','dltype'); 
    a0:=`int/cook/secure`(a0,'a0fl','a0type'); 
    a1:=`int/cook/secure`(a1,'a1fl','a1type'); 
    if d<>0 then 
      np1:=(w+1)/d 
    else 
      fail:=2; 
      `int/cook/nogo1`(MODEL,model,t,ta,tb,s); 
      return 
    end if; 
    if not type(pfm1l,constant) then 
      fail:=2; 
      `int/cook/nogo2`(MODEL,model,t,ta,tb,p); 
      return 
    end if; 
    if not type(wfp1l,constant) then 
      fail:=2; 
      `int/cook/nogo2`(MODEL,model,t,ta,tb,w+1); 
      return 
    end if; 
    if type([a0,a1],[constant,constant]) and 
       type([a1fl,a0fl],[numeric,numeric]) and 
       abs(a0fl)<abs(a1fl) 
    then 
      `int/cook/nogo1`(MODEL,model,t,ta,tb,abs(a1)-abs(a0)); 
      return 
    end if; 
    if not type(dfl,constant) then 
      fail:=2; 
      `int/cook/nogo2`(MODEL,model,t,ta,tb,d); 
      return 
    end if; 
    if u=0 and u2=0 and m=0 and c=0 and sp=0 and 
       f=t and pfm1l<1 and 0<dfl and abs(a0/a1)<>1 
    then 
      sig:=signum(a0); 
      if not type(p,integer) and (type(sig,'function') and 
         op(0,sig)=('signum') or sig=-1) 
      then 
        fail:=2; 
        if 2<printlevel then 
          `int/cook/nogo1`(MODEL,model,t,ta,tb); 
          printf("--> Does not fit into this sub-class\n") 
        end if; 
        return 
      end if; 
      s:=d; 
      `int/cook/ngritty`:=proc(s,dc1,np1,p,dc2,a0,a1,dc3,dc4) 
      local lam,b,Kt,F; 
        lam:=np1; 
        b:=p; 
        Kt:=-a1/a0; 
        F:=Beta(lam,1)*hypergeom([b,lam],[lam+1],Kt); 
        simplify(a0^(-p))*simplify(F)/s 
      end proc 
    elif u=0 and u2=0 and a1=-a0 and c=0 and sp=0 and f=t and 
         pfm1l<=0 and 0<dfl and type(m,integer) and 0<=m and 
         type(d,freeof(t)) 
    then 
      if ta=0 and wfp1l<=0 then 
        z:=limit(t*model,t=0,'right'); 
        if hastype(z,'infinity') then 
          v:=limit((t-1)*model,t=1,'left'); 
          if signum(z)=signum(v) or not hastype(v,'infinity') then 
            ans:=z; 
            fail:=4; 
            return 
          elif has(v,'infinity') then 
            fail:=3; 
            if 2<printlevel then 
              `int/cook/nogo1`(MODEL,model,t,ta,tb); 
              printf("integral is most likely divergent\n"); 
              printf("Cauchy principal value ?\n") 
            end if; 
            return 
          end if 
        end if 
      end if; 
      s:=d; 
      z:=np1; 
      v:=-p+1; 
      `int/cook/ngritty`:=proc(s,m,np1,p,z0,a0,a1,dc2,v) 
      local z,w,F; 
        F:=Beta(z,w); 
        if 0<m then F:=diff(F,`$`(z,m)) end if; 
        try 
          if v=0 then                                             #<== 
            simplify(                                             #<== 
              a0^(-p))*simplify(limit(eval(F,z=z0),w=v))/(s^(m+1) #<== 
            )                                                     #<== 
          else                                                    #<== 
            simplify(                                             #<== 
              a0^(-p))*simplify(eval(F,{z=z0,w=v}))/(s^(m+1)      #<== 
            )                                                     #<== 
          end if                                                  #<== 
        catch "numeric exception: division by zero": 
          a0^(-p)/(s^(m+1))*infinity 
        end try 
      end proc: 
    else 
      fail:=2; 
      if 2<printlevel then 
        `int/cook/nogo1`(MODEL,model,t,ta,tb); 
        printf("--> Does not fit into sub-classes\n") 
      end if; 
      return 
    end if; 
    if (bfl<>`int/cook/insecure`(1,btype) or 3<btype) and 
        bfl<>0 and dlfl<>0 
    then 
      res:=0; 
      for i from 0 to m do 
        res:=res+ 
          binomial(m,i)*ln(b)^(m-i)*dl^i* 
          `int/cook/ngritty`(s,i,np1,p,z,a0,a1,ucomplex,v) 
      end do 
    elif bfl=`int/cook/insecure`(1,btype) and dlfl<>0 then 
      res:=dl^m*`int/cook/ngritty`(s,m,np1,p,z,a0,a1,ucomplex,v) 
    else 
      res:=`int/cook/ngritty`(s,m,np1,p,z,a0,a1,ucomplex,v) 
    end if; 
    res:=cof*res; 
    if has(res,abs) then res:=`int/cook/abs`(res) end if; 
    if has(res,'infinity') then failtemp:=4 end if; 
    if fltypetemp=1 then res:=evalf(res) end if; 
    ans:=eval(res); 
    fltype:=fltypetemp; 
    fail:=failtemp 
  end proc:
 
> int(ln(x)/(1-x),x=0..1); 
 
                                      2 
                              - 1/6 Pi 
 
> int(ln(x)^2/(1-x)^2,x=0..1); 
 
                                     2 
                               1/3 Pi 
 
> j:=ln(x)^3/(1-x)^3; 
 
                                       3 
                                  ln(x) 
                            j := -------- 
                                        3 
                                 (1 - x) 
 
> J:=applyop(expand,1,int(j,x=0..1)); 
 
                                   2   1/2       ln(1 - exp(1/y)) 
  J :=   lim    -3 Zeta(3) + 1/2 Pi  + --- - 3/2 ---------------- 
       y -> 0+                          3                2 
                                       y                y 
 
           3 polylog(2, exp(1/y))                            3/2 
         - ---------------------- + 3 polylog(3, exp(1/y)) + --- 
                     y                                        2 
                                                             y 
 
           3 ln(1 - exp(1/y)) 
         - ------------------ - 3 polylog(2, exp(1/y)) 
                   y
 

Maple cannot compute that limit symbolically. However, it will return the numerical value of J:

> simplify(fnormal(evalf(J)),zero); 
 
                             -8.540972910
 

Let us check this result:

> int(j,x=a..1); 
 
                              1 
                             /        3 
                            |    ln(x) 
                            |   -------- dx 
                            |          3 
                           /    (1 - x) 
                             a 
 
> limit(%,a=0,right); 
 
                                2 
                        - 1/2 Pi  - 3 Zeta(3) 
 
> evalf(%); 
 
                             -8.540972911