List Archives >
Maple User Group List Archive >
Archive by date >
This Month By Date >
This Month By Topic
[MUG] Re: 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 "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
|
[View Complete Thread]
Previous by date: [MUG] results differ in plot3d(), Zlatko Savic
Next by date: [MUG] Re: simplifying a sum of integrals, Maple User Group
Previous thread: [MUG] Problems with solving systems of equations, Roush, Craig Ryan UMKC-Student
Next thread: [MUG] simplifying a sum of integrals, Sam H Davis
|
(NAG)
