#### 7.85 bug in integration dependent on N, Maple 6 to Maple 8 (5.4.02)

##### 7.85.1 Glenn Sowell

I have seen the following integration bug in Maple 6 & 7, but not in 5.1.

The bug seems to be N dependent, where N is the power of the denominator of the integrand.

 > interface(version);
Maple Worksheet Interface, Maple 7.00,
APPLE_PPC_MAC, Tue, May 29, 2001 Build ID 96223
> assume(a>0);
> f := 1/(x^2+a^2)^N;
> # N=6 This does not work.
> N:=6:
> Df := diff(f,x,x);
> fDf := simplify( f * Df );
> Int1 := int( fDf, x=0..infinity );

Int1 := 0



A surprising result!

Try expanding fDf before integrating:

 > Int2 := int( expand(fDf), x=0..infinity );

468027 Pi
Int2 := - ------------
25
1048576 a~



A diﬀerent result?!

he same dichotomy occurs for N>=6. For N<=5, both results are the same.

It gets stranger...

Do the indeﬁnite integral ﬁrst

 > Int3 := int( fDf, x);



Now apply the limit by hand:

 > Int2 := limit( Int3, x=infinity ) - limit( Int3, x=0 );

468027 Pi
Int2 := - ------------
25
1048576 a~



Another surprising result!

Evidently, the deﬁnite integral algorithm is doing something wrong!

Why does it occur for N>=6 and not for N<=5?

It is corrected with Maple 9 (U. Klein)

##### 7.85.2 Helmut Kahovec (12.4.02)

Actually, the procedure int/definite/contour/residue, which computes residues during deﬁnite integration, is ﬂawed. Additionally, the procedure 'residue' cannot compute a residue if the order of the pole is suﬃciently high. Following is a Maple7 session on a PC running MS Windows NT 4.0.

> assume(a>0);
> num,den,alpha,m:=(13*_X^2-a^2)*ln(-_X),(_X^2+a^2)^14,I*a,14:

> residue(eval(num/den,_X=x),x=alpha);

2    2
(13 x  - a ) ln(-x)
residue(-------------------, x = I a)
2    2 14
(x  + a )



’residue’ returns unevaluated. We now modify ’residue’ as indicated:

> unprotect(residue);

> residue:=proc(f,a)
local g,i,t,x,t1;
options remember;
if nargs<>2 or not type(a,equation) or not type(op(1,a),name) then
error "invalid arguments"
end if;
x:=op(1,a);
if has(op(2,a),x) then
error "invalid point"
elif type(op(2,a),infinity) then
g:=-1/x^2*subs(x=1/x,f)
else
g:=normal(subs(x=x+op(2,a),f),expanded)
end if;
#    for i from 0 to 5 do   <=== replace the number 5 by 6 ===
for i from 0 to 6 do
t:=traperror(series(g,x,i^2+2));
if t=0 then return 0 end if;
if t=lasterror or not type(t,series) then next end if;
t1:=order(t);
if t1=infinity or -1<t1 then return coeff(t,x,-1) end if
end do;
'residue(args)'
end proc:

> protect(residue);

> residue(eval(num/den,_X=x),x=alpha): R11:=simplify(expand(%));

R11 :=

669278610 I ln(a) - 2269536821 I + 334639305 Pi
1/17993564160 -----------------------------------------------
25
a



This seems to be the correct value. During deﬁnite integration the procedure int/definite/contour/residue gets called with the appropriate arguments. However, since int/definite/contour/residue is ﬂawed, it will return an incorrect residue if the order of the pole is high enough.

> int/definite/contour/residue(num,den,alpha,m);

-79
-------- I
14057472
----------
25
a



We ﬁx the incorrect lines in int/definite/contour/residue as indicated:

> int/definite/contour/residue:=proc(n,d,a,m)
local i,r,v1;
if nargs=4 and m=1 then
i:=n/diff(d,_X);
i:=traperror(eval(subs(_X=a,i)));
if i<>lasterror then return i end if
end if;
for i from 0 by 3 to 12 do
r:=traperror(int/definite/contour/seriesc(n/d,_X=a,i));
if r=FAIL then
return FAIL
elif r=lasterror or
not type(r,series) or
op(nops(r),r)<0 and op(nops(r)-1,r)=O(1)
then
next
else
break
end if
end do;
if i<=12 then
if hastype(n/d,'nonreal') or hastype(a,'nonreal') then
return coeff(r,evalc(_X-a),-1)
else
return coeff(r,_X-a,-1)
end if
end if;
if type(d,polynom(anything,_X)) and nargs=4 and 1<m then
v1:=quo(d,(_X-a)^m,_X,'r');
if r<>0 then return FAIL end if;
#
#     One has to differentiate n/v1 instead of differentiating n
#     alone and dividing by v1 afterwards!
#
#                                       |||
#                                       vvv
i:=traperror(eval(subs(_X=a,diff(n/v1,$(_X,m-1))/(m-1)!))); if i<>lasterror then return i else return FAIL end if elif nargs=4 and 1<m then # # One has to differentiate (_X-a)^m*(n/d) instead of (_X-a)^m # alone! # |||| # vvvv i:=limit(diff((_X-a)^m*n/d,$(_X,m-1))/(m-1)!,_X=a);
if type(i,function) and op(0,i)=limit or
type(i,undefined) or
type(i,range)
then
FAIL
else
i
end if
else
FAIL
end if
end proc:



Now, int/definite/contour/residue returns the same result as ’residue’:

> int/definite/contour/residue(num,den,alpha,m):
> R12:=simplify(expand(%));

R12 :=

669278610 I ln(a) - 2269536821 I + 334639305 Pi
1/17993564160 -----------------------------------------------
25
a



The other pole of your integrand gives the following residue:

> num,den,alpha,m:=(13*_X^2-a^2)*ln(-_X),(_X^2+a^2)^14,-I*a,14:
> residue(eval(num/den,_X=x),x=alpha): R21:=simplify(expand(%));

R21 := - 1/17993564160

669278610 I ln(a) - 334639305 Pi - 2269536821 I
-----------------------------------------------
25
a

> int/definite/contour/residue(num,den,alpha,m):
> R22:=simplify(expand(%));

R22 := - 1/17993564160

669278610 I ln(a) - 334639305 Pi - 2269536821 I
-----------------------------------------------
25
a



Finally, Maple will correctly compute your deﬁnite integral, too:

> f:=1/(x^2+a^2)^N:
> N:=6:
> Df:=diff(f,x,x):
> fDf:=simplify(f*Df);

2    2
13 x  - a
fDf := 12 -----------
2    2 14
(x  + a )

> int(fDf,x=0..infinity): normal(%,expanded);

468027  Pi
- ------- ---
1048576  25
a



If we expand the integrand then Maple ﬁrst computes the antiderivatives and takes the limits afterwards.

> int(expand(fDf),x=0..infinity): normal(%);

468027  Pi
- ------- ---
1048576  25
a