[MUG] Integration bug dependent on N

[MUG] Re: Integration bug dependent on N
Author: Helmut Kahovec    Posted: Fri, 12 Apr 2002 01:26:21 +0200
>> From: Helmut Kahovec /> Glenn Sowell wrote: >| 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. 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. > restart; > 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 Kind regards Helmut

Previous by date: [MUG] Re: Parallel use?,  Jerome BENOIT
Next by date: [MUG] Stop reducing & data points.,  Jason C Leach
Previous thread: [MUG] is there an easy way to compute a Fourier series spectrum of raw data taken from text file,  Deveney, Edward
Next thread: [MUG] Stop reducing & data points.,  Jason C Leach