There seems to be a bug in Maple 6. When I do the following computation I get different
results for the numerical Intergration. In the first case I took as limits q=0..2*Pi
and in the
second q=0..2*evalf(Pi)
. In the first case the result is wrong. In Maple V both
computations give the same result.
Is this a bug ?
> restart: > with(linalg): > h:=matrix(9,1,[a1*cos(q),a1*sin(q),q,r*cos(q)+a2*sqrt(1-lambda^2*(sin(q))^2), (r-lambda*a2)*sin(q),-arcsin(lambda*sin(q)), r*cos(q)+l*sqrt(1-lambda^2*(sin(q))^2),0,0]); > hst:=map(diff,h,q); > M:=diag(m1,m1,I1,m2,m2,I2,m3,m3,I3); > Theta:=evalm(transpose(hst)&* M &* hst); > l:=0.25; > r:=0.03; > lambda:=r/l; > a1:=21.877e-3; > m1:=50.32; > I1:=1.0267; > a2:=0.034345; > m2:=16.645; > I2:=0.23251; > m3:=44.55; > Theta_m:=evalf(int(1/2/Pi*Theta[1,1],q=0..2*Pi)); Theta_m := 1.058110253 > Theta_m:=evalf(int(1/2/evalf(Pi)*Theta[1,1],q=0..2*evalf(Pi))); Theta_m := 1.085648004
es, this is a bug. But it’s not numerical integration, it’s symbolic integration. If you want
numerical integration, use evalf(Int(...))
not evalf(int(...))
.
Here is this bug in a simpler form:
> h:= (1-cos(t)/sqrt(1-a*sin(t)^2))^2; > int(expand(h),t) - int(h,t); t
The integral of the expanded form of \(h\) would be correct (at least when \(a\) is small and positive).
Somehow it’s as if int(h,t)
expands \(h\) but forgets about the 1. In your example, if you use
expand(Theta[1,1])
instead of Theta[1,1]
you should get the correct answer 1.085648004
both times.
All I found out is that if the (trivial) changes described below are made in
`int/indefinite`()
and `int/definite`()
then the example input will give the correct
results.
Unfortunately, that statement is not correct, it is a bug, really. I am sorry and congratulate Robert Israel!
Well, I did not at all expect the symbolic integration to be wrong. Note that releases 4 and 5 of Maple V cannot evaluate that integral symbolically and thus numerically integrate definite integrals like
> lambda:=1/2: > j:=(1+lambda^2*cos(q)/(sqrt(1-lambda^2*sin(q)^2)))^2; > evalf(int(j,q=0..2*Pi));
maple6, however, can indeed integrate
> int(j,q);
and
> int(j,q=0..2*Pi);
symbolically but, as Robert has shown, will give the wrong result. `int/indefinite`()
and
`int/definite`()
should expand their first argument in this case. I cannot pinpoint the
bug more closely than that. If we make the following trivial changes to `int/indefinite`()
and `int/definite`()
`int/indefinite` := proc(FF, zz) local ff, opts, answer; global StandardFunctions; 1 ff := expand(FF); 2 if nargs = 2 then ... 19 if answer = FAIL then 20 answer := ('int')(ff,zz) end if; 21 answer end proc `int/definite` := proc(FF, zz, aa, bb) local ff, answer, opts; 1 ff := expand(FF); 2 if nargs = 5 then ... 33 if answer = FAIL then 34 answer := ('int')( ff,zz = aa .. bb,op(subs(('formula') = NULL,opts)) ) end if; 35 RETURN(answer) end proc
then the following integrals will give the correct results:
> int(j,q); q + arcsin(1/2 sin(q)) - 1/8 sqrt(3) arctan(1/3 (2 tan(1/2 q) - 1) sqrt(3)) - 1/8 sqrt(3) arctan(1/3 (2 tan(1/2 q) + 1) sqrt(3)) + 1/2 arctan(tan(1/2 q)) > int(j,q=0..2*Pi); 5/2 Pi - 1/4 sqrt(3) Pi
Note, too, that maple6 integrates
> int(j,q=0..2*evalf(Pi));
numerically right from the beginning, cf. lines 1-4
of the procedure int()
.
Please note that V 5.1 will integrate it symbolically with the answer given in terms of elliptic integrals.