There is a bug in the definite integration algorithm which only appears when the limits of
integration are floating point numbers. Expressions of the form "limit/X" appear uninvited
in the result. I had thought that this problem had been fixed in the transition from V4
to V5, but it is still around. In the first example below, the limits of integration
are rational numbers, so the problem does not appear. In the second, the same
limits have been written as floating point numbers, and the result is suddenly
useless.
|\^/| Maple V Release 5 (WMI Campus Wide License) ._|\| |/|_. Copyright (c) 1981-1997 by Waterloo Maple Inc. All rights \ MAPLE / reserved. Maple and Maple V are registered trademarks of <____ ____> Waterloo Maple Inc. | Type ? for help. > int( x/sin(1-x^2), x=1/10..9/10 ); 19 19 19 - 1/2 ln(1 - cot(---) sin(---)) + 1/2 ln(sin(---)) 100 100 100 99 99 99 + 1/2 ln(1 - cot(---) sin(---)) - 1/2 ln(sin(---)) 100 100 100 > int( x/sin(1-x^2), x=0.1..0.9 ); 2 .8671793960 + 1.570796327 I signum(csc(-.19 - 1.8 limit/X + 1. limit/X ) 2 - 1. cot(-.19 - 1.8 limit/X + 1. limit/X )) - 1.570796327 I signum( 2 2 csc(-.99 + .2 limit/X + limit/X ) - 1. cot(-.99 + .2 limit/X + limit/X ))
It is corrected with Maple 6. (U. Klein)
Yech. It does clear up, however, if you add abs() around the argument of sin():
> int( x/sin(abs(1-x^2)), x=0.1..0.9 ); .8671793961
which agrees with the rational version.
This bug occurs in Release 4 and Release 5 because in your example
- int calculates the antiderivative of the integrand and then calls int/DefInt.
- int/Defint calls limit to calculate the value of the antiderivative at the limits of integration.
- limit substitutes L-`limit/X` for the indeterminate in the antiderivative. L is either limit
of integration.
- limit does a series expansion of the antiderivative w.r.t. `limit/X` and calls
series/leadterm since L is a float and there isn’t any float in the antiderivative.
- series/leadterm does not work correctly if the series expansion of the antiderivative doesn’t have an explicit leading term. It always takes the first term of the series expansion.
Thus the real bug is in series/leadterm. Look at the following:
> restart; > j:=x/sin(1-x^2); x j := - ------------ 2 sin(-1 + x ) > Int(j,x): %=value(%); / | x 2 2 | - ------------ dx = - 1/2 ln(csc(-1 + x ) - cot(-1 + x )) | 2 / sin(-1 + x ) > J:=rhs(%); 2 2 J := - 1/2 ln(csc(-1 + x ) - cot(-1 + x )) > limit(J,x=0.1,right); .3082524570 - 1.570796327 I csgn(I ( 2 -csc(-.99 + .2 limit/X + limit/X ) 2 + 1. cot(-.99 + .2 limit/X + limit/X ))) > limit(J,x=0.9,left); 1.175431853 + 1.570796327 I csgn(I ( 2 csc(-.19 - 1.8 limit/X + 1. limit/X ) 2 - 1. cot(-.19 - 1.8 limit/X + 1. limit/X )))
Compare with the following:
> limit(J,x=1/10,right); 99 99 99 - 1/2 ln(1 - cot(---) sin(---)) + 1/2 ln(sin(---)) - 1/2 I Pi 100 100 100 > limit(J,x=9/10,left); 19 19 19 - 1/2 ln(1 - cot(---) sin(---)) + 1/2 ln(sin(---)) - 1/2 I Pi 100 100 100
Actually, the order of the series expansion is 6. We use 4 to save space:
> series(J,x,4); / 1 |- 1/2 ln(------ - cot(1)) \ sin(1) 2 2 \ + 1/2 I csgn(I (csc(-1 + x ) - cot(-1 + x ))) Pi| - / cos(1) 2 ------- - 1 - cot(1) 2 sin(1) 2 4 1/2 --------------------- x + O(x ) 1 ------ - cot(1) sin(1) > series(leadterm(J),x,4); / 1 |- 1/2 ln(------ - cot(1)) \ sin(1) 2 2 \ - 1/2 I csgn(I (-csc(-1 + x ) + cot(-1 + x ))) Pi| /
Note that if L is a float, then series/leadterm seems to be called if there also is a float in the integrand. But if you manually introduce a float in the antiderivative then series/leadterm isn’t called by limit and the antiderivative is calculated correctly at the limits of integration:
> subs(-1+x^2=-1.+x^2,J); 2 2 - 1/2 ln(csc(-1. + x ) - cot(-1. + x )) > limit(%,x=0.1,right); .3082524570 - 1.570796327 I > limit(%%,x=0.9,left); 1.175431854 - 1.570796327 I
Now we may easily construct another example of this bug:
> -ln(sin(-1+x^2)); 2 -ln(sin(-1 + x )) > series(%,x,6); 2 cos(1) 2 (-ln(sin(1)) + I csgn(I sin(-1 + x )) Pi) + ------ x + sin(1) / 2\ | cos(1) | 4 6 |1/2 + 1/2 -------| x + O(x ) | 2| \ sin(1) / > series(leadterm(%%),x,6); 2 (-ln(sin(1)) + I csgn(I sin(-1 + x )) Pi) > diff(%%%,x); 2 cos(-1 + x ) x -2 -------------- 2 sin(-1 + x ) > int(%,x); 2 -ln(sin(-1 + x )) > int(%%,x=0.1..0.9); 1.487659539 2 + 3.141592654 I csgn(I sin(-.19 - 1.8 limit/X + 1. limit/X )) 2 - 1. I csgn(I sin(-.99 + .2 limit/X + limit/X )) Pi > int(%%%,x=1/10..9/10); 19 99 -ln(sin(---)) + ln(sin(---)) 100 100