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 different result?!
he same dichotomy occurs for N>=6. For N<=5, both results are the same.
It gets stranger...
Do the indefinite integral first
> 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 definite 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)
Actually, the procedure `int/definite/contour/residue`, which computes residues
during definite integration, is flawed. Additionally, the procedure 'residue' cannot compute
a residue if the order of the pole is sufficiently 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 definite integration the procedure
`int/definite/contour/residue` gets called with the appropriate arguments. However,
since `int/definite/contour/residue` is flawed, 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 fix 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 definite 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 first computes the antiderivatives and takes the limits afterwards.
> int(expand(fDf),x=0..infinity): normal(%); 468027 Pi - ------- --- 1048576 25 a